From f204496e1aad57a97de316a0407cb9935e65ebd1 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sat, 15 Mar 2025 14:12:05 +0100 Subject: [PATCH] detect propagation past a cutoff MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit There are two cases in which a wave in GRAY can reach and pass through a cutoff: 1. a discontinuity or too large step size causes the integrator to jump into region with N⊥² < 0, according to the cold plasma theory. 2. a wave gets close to a cutoff region is still propagating according to the cold plasma theory, but the warm disperion relation already gives re(N⊥²) < 0. For case 1. we raise an error and immediately stop the raytracing as it must be completely invalid. For case 2. we force N⊥² = 0 to avoid interpreting the wave reflection as a strong energy absorption, but allow the raytracing to continue. --- src/dispersion.f90 | 18 +++++++++++--- src/gray_core.f90 | 38 ++++++++++++++++++++---------- src/gray_errors.f90 | 57 ++++++++++++++++++++++++++++++--------------- 3 files changed, 79 insertions(+), 34 deletions(-) diff --git a/src/dispersion.f90 b/src/dispersion.f90 index a0d8744..b5777ae 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -26,7 +26,7 @@ module dispersion real(wp_), parameter :: extv(*) = exp(-ttv**2) private - public colddisp, warmdisp + public colddisp, warmdisp, dielectric_tensor public zetac, harmnumber contains @@ -438,11 +438,23 @@ subroutine warmdisp(X, Y, mu, Npl, Npr_cold, sox, & error = raise_error(error, warmdisp_result, 0) end if - if (Npr2%re < 0 .and. Npr2%im < 0) then - Npr2 = sqrt(0.5_wp_) ! Lorenzo, what should we have here?? + ! Catch unphysical growing waves + if (Npr2%re < 0 .and. Npr2%im > 0) then + ! Set N⊥²=0 to avoid altering the wave residual power + ! Note: This is impossible in a Maxwellian plasma, so + ! it must be a numerical problem. + Npr2 = 0 error = raise_error(error, warmdisp_result, 1) end if + ! Catch propagation past a cut-off + if (Npr2%re < 0 .and. Npr2%im < 0) then + ! Set N⊥²=0 to avoid spurious absorption + ! Note: im(N⊥)<0 here just means the wave is being reflected. + Npr2 = 0 + error = raise_error(error, warmdisp_result, 2) + end if + ! Final result Npr = sqrt(Npr2) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 140992c..a1d0ed2 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -295,10 +295,6 @@ contains if (igrad_b == 1) call gradi_upd(params%raytracing, yw, ak0, xc, & du1, gri, ggri, curr_error) - if (has_new_errors(prev_error, curr_error)) then - call print_err_raytracing(curr_error, i, jk) - end if - iopmin = 10 ! =========== rays loop BEGIN =========== @@ -312,7 +308,7 @@ contains xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & xgcn,ypw(:,jk), psinv, dens, btot, bv, xg, yg, & derxg, anpl, anpr, ddr, ddi, dersdst, derdnm, & - igrad_b) + igrad_b, curr_error) ! check entrance/exit plasma/wall zzm = xv(3)*0.01_wp_ @@ -413,7 +409,7 @@ contains xv, anv, gri(:,jk), ggri(:,:,jk), sox, bres, & xgcn, ypw(:,jk), psinv, dens, btot, bv, xg, & yg, derxg, anpl, anpr, ddr, ddi, dersdst, & - derdnm, igrad_b) ! * update derivatives after reflection + derdnm, igrad_b, curr_error) ! * update derivatives after reflection if(jk == 1 .and. ip == 1) then ! * 1st pass, polarization angles at reflection for central ray psipv(index_rt) = psipol @@ -525,6 +521,10 @@ contains end do rays_loop ! ============ rays loop END ============ + if (has_new_errors(prev_error, curr_error)) then + call print_err_raytracing(curr_error, i, jk) + end if + if(i==params%raytracing%nstep) then ! step limit reached? do jk=1,params%raytracing%nray if(iop(jk)<3) call turnoffray(jk,ip,npass,ib,iroff) ! * ray hasn't exited+reentered the plasma by last step => stop ray @@ -945,12 +945,13 @@ contains subroutine ywppla_upd(params, equil, plasma, & xv, anv, dgr, ddgr, sox, bres, xgcn, dery, & psinv, dens, btot, bv, xg, yg, derxg, anpl, & - anpr, ddr, ddi, dersdst, derdnm, igrad) + anpr, ddr, ddi, dersdst, derdnm, igrad, error) ! Compute right-hand side terms of the ray equations (dery) ! used after full R-K step and grad(S_I) update use gray_params, only : gray_parameters use gray_equil, only : abstract_equil use gray_plasma, only : abstract_plasma + use gray_errors, only : gray_error ! subroutine rguments type(gray_parameters), intent(in) :: params @@ -968,6 +969,7 @@ contains real(wp_), intent(out) :: bv(3) real(wp_), intent(out) :: derxg(3) integer, intent(in) :: igrad + integer(kind=gray_error), intent(inout) :: error ! local variables real(wp_):: deryg(3) @@ -975,8 +977,9 @@ contains call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg, psinv) - call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, dgr, & - ddgr, igrad, dery, anpl, anpr, ddr, ddi, dersdst, derdnm) + call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, & + dgr, ddgr, igrad, dery, anpl, anpr, ddr, ddi, dersdst, & + derdnm, error) end subroutine ywppla_upd @@ -1369,7 +1372,7 @@ contains subroutine disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, & dgr, ddgr, igrad, dery, anpl_, anpr, ddr, ddi, & - dersdst, derdnm_) + dersdst, derdnm_, error) ! Computes the dispersion relation, derivatives and other ! related quantities ! @@ -1381,6 +1384,7 @@ contains use const_and_precisions, only : zero, half use gray_params, only : gray_parameters, STEP_ARCLEN, & STEP_TIME, STEP_PHASE + use gray_errors, only : gray_error, cutoff, raise_error ! subroutine arguments @@ -1416,6 +1420,8 @@ contains real(wp_), optional, intent(out) :: dersdst ! |∂Λ/∂N̅| real(wp_), optional, intent(out) :: derdnm_ + ! raytracing error + integer(kind=gray_error), optional, intent(inout) :: error ! Note: assign values to missing optional arguments results in a segfault. ! Since some are always needed anyway, we store them here and copy @@ -1620,11 +1626,19 @@ contains dery(1:3) = derdnv(:)/denom ! +∂Λ/∂N̅ dery(4:6) = -derdxv(:)/denom ! -∂Λ/∂x̅ + ! Ns², solution of the cold dispersion relation + an2s = 1 - xg - half*xg*yg2*(1 + anpl2 + sox*del)/duh + + if (present(error) .and. an2s-anpl2 < 0) then + ! The wave is past a cutoff (N²⊥ < 0) + ! Note: this is likely due to an integration error + error = raise_error(error, cutoff) + end if + if (present(ddr) .or. present(ddi)) then ! Dispersion relation (geometric optics) ! ddr ← Λ = N² - N²s(X,Y,N∥) = 0 - an2s = 1 - xg - half*xg*yg2*(1 + anpl2 + sox*del)/duh - ddr = an2 - an2s + ddr = an2 - an2s end if if (present(ddr) .and. igrad > 0) then diff --git a/src/gray_errors.f90 b/src/gray_errors.f90 index a5f7a79..8c1f2a8 100644 --- a/src/gray_errors.f90 +++ b/src/gray_errors.f90 @@ -5,7 +5,7 @@ ! function and combined simply with the intrinsic `ior` function. module gray_errors - use logger, only : log_warning + use logger, only : log_warning, log_error use, intrinsic :: iso_fortran_env, only : int32 implicit none @@ -33,11 +33,16 @@ module gray_errors type(error_spec), parameter :: unstable_beam = & error_spec(offset=0, subcases=1, & - mod='gray_core', proc='gray_main', & + mod='gray_core', proc='gradi_upd', & msg=[str :: 'beamtracing may be unstable', pad1]) + type(error_spec), parameter :: cutoff = & + error_spec(offset=after(unstable_beam), subcases=1, & + mod='gray_core', proc='disp_deriv', & + msg=[str :: 'N⊥²<0, passed through a cutoff', pad1]) + type(error_spec), parameter :: dielectric_tensor = & - error_spec(offset=after(unstable_beam), subcases=2, & + error_spec(offset=after(cutoff), subcases=2, & mod='gray_core', proc='gray_main', & msg=[str :: 'ε tensor, overflow in `fsup`', & 'ε tensor, integration error in `hermitian_2`', & @@ -50,12 +55,13 @@ module gray_errors 'failed to converge, returned last value', & pad2]) - type(error_spec), parameter :: warmdisp_result = & - error_spec(offset=after(warmdisp_convergence), subcases=2, & - mod='dispersion', proc='warmdisp', & - msg=[str :: 'final N⊥² is NaN or ±Infinity', & - 'final N⊥² in 3rd quadrant', & - pad2]) + type(error_spec), parameter :: warmdisp_result = & + error_spec(offset=after(warmdisp_convergence), subcases=3, & + mod='dispersion', proc='warmdisp', & + msg=[str :: 'final N⊥² is NaN or ±Infinity', & + 'final N⊥² had re(N⊥²)<0, wave past cutoff', & + 'final N⊥² had im(N⊥²)>0, unphysical growing wave', & + pad3]) type(error_spec), parameter :: negative_absorption = & error_spec(offset=after(warmdisp_result), subcases=1, & @@ -76,7 +82,7 @@ module gray_errors pad3]) ! Errors occuring during raytracing - type(error_spec), parameter :: raytracing_errors(*) = [unstable_beam] + type(error_spec), parameter :: raytracing_errors(*) = [unstable_beam, cutoff] ! Errors occuring during absorption and current drive computations type(error_spec), parameter :: ecrh_cd_errors(*) = & @@ -96,7 +102,8 @@ contains integer(kind=gray_error), intent(in) :: error logical :: is_critical - is_critical = has_error(error, negative_absorption) + is_critical = has_error(error, negative_absorption) .or. & + has_error(error, cutoff) end function is_critical @@ -154,10 +161,10 @@ contains character(256) :: line ! formatted log line type(error_spec) :: spec integer :: i, j + integer(kind=gray_error) :: tmp ! format specifier of the log line - character(*), parameter :: fmt = & - '(a,": ","error=",g0," step=",g0," ray=",g0)' + character(*), parameter :: fmt = '(a,": ","step=",g0," ray=",g0)' ! iterate on the known errors do i = 1, size(raytracing_errors) @@ -168,8 +175,14 @@ contains ! iterate on the subcases do j = 0, spec%subcases - 1 if (ibits(slice, j, 1) == 0) cycle - write(line, fmt) trim(spec%msg(j)), slice * 2**j, step, ray - call log_warning(line, mod=spec%mod, proc=spec%proc) + write(line, fmt) trim(spec%msg(j)), step, ray + + if (is_critical(raise_error(tmp, spec, j))) then + call log_error(line, mod=spec%mod, proc=spec%proc) + else + call log_warning(line, mod=spec%mod, proc=spec%proc) + end if + end do end do end subroutine print_err_raytracing @@ -193,10 +206,11 @@ contains character(256) :: line ! formatted log line type(error_spec) :: spec integer :: i, j + integer(kind=gray_error) :: tmp ! format specifier of the log line character(*), parameter :: fmt = & - '(a,": ", "error=",g0, " N⊥=",f0.0,sp,f0.0,"i", " α=",g0, " step=",g0," ray=",g0)' + '(a,": ", "N⊥=",g0.3,sp,g0.3,"i", " α=",g0.3, " step=",g0," ray=",g0)' ! iterate on the known errors do i = 1, size(ecrh_cd_errors) @@ -207,9 +221,14 @@ contains ! iterate on the subcases do j = 0, spec%subcases - 1 if (ibits(slice, j, 1) == 0) cycle - write(line, fmt) trim(spec%msg(j)), slice * 2**j, & - Npr%re, Npr%im, alpha, step, ray - call log_warning(line, mod=spec%mod, proc=spec%proc) + write(line, fmt) trim(spec%msg(j)), Npr%re, Npr%im, alpha, step, ray + + if (is_critical(raise_error(tmp, spec, j))) then + call log_error(line, mod=spec%mod, proc=spec%proc) + else + call log_warning(line, mod=spec%mod, proc=spec%proc) + end if + end do end do end subroutine print_err_ecrh_cd