detect propagation past a cutoff
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.
This commit is contained in:
Michele Guerini Rocco
parent
e3eca72dff
commit
f204496e1a
GPG Key ID:
BFBAF4C975F76450
@ -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)
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
|
||||
Loading…
Reference in New Issue
Block a user