diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 5a2add2..f71dcb7 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -2,30 +2,69 @@ module gray_core use const_and_precisions, only : wp_ + use gray_params, only : gray_parameters + use gray_equil, only : abstract_equil + use gray_plasma, only : abstract_plasma + use gray_errors, only : gray_error, raise_error, has_error, has_new_errors implicit none + type beam_state + ! Note: this should be converted into a PDT once gfotran can handle it + + integer :: nray ! total number of rays + integer :: nrayr, nrayth ! radial/angular number of rays + + ! Main beam variables + real(wp_), allocatable :: y(:, :) ! position and refractive index vector (6, nray) + real(wp_), allocatable :: yp(:, :) ! dy/dσ, with σ integration var (6, nray) + integer :: step ! current number of integration steps + real(wp_), allocatable :: dist(:) ! distance travelled in integration var (nray) + logical :: beamtracing ! whether beamtracing is enabled + + ! Complex eikonal + real(wp_), allocatable :: pos(:, :, :) ! ray positions (3, nrayth, nrayr) + real(wp_), allocatable :: grad_u(:, :, :) ! variations Δu (3, nrayth, nrayr) + real(wp_), allocatable :: grad_S(:, :) ! gradient of S_I (3, nray) + real(wp_), allocatable :: hess_S(:, :, :) ! Hessian of S_I (3, 3, nray) + + ! Local quantities of each ray + real(wp_), allocatable :: Xg(:), Yg(:) ! normalised density, magnetic field + real(wp_), allocatable :: N_pl(:), N_pr(:) ! N∥, Nَ⊥ solution of dispersion rel + real(wp_), allocatable :: dsdvar(:) ! ds/dσ, with σ integration var + real(wp_), allocatable :: dlambdadN(:) ! |∂Λ/∂N|, with Λ dispersion rel + + ! Polarisation of each ray + integer, allocatable :: mode(:) ! -1 ⇒ O mode, +1 ⇒ X mode + + ! ECRH & CD quantities of each ray + real(wp_), allocatable :: psi_n(:) ! normalised poloidal flux + real(wp_), allocatable :: tau(:) ! optical depth + real(wp_), allocatable :: I_cd(:) ! total driven current + + ! Beam constants + real(wp_) :: B_res, X_norm, k0, rwmax + end type + + interface + function rhs(y) + import :: wp_ + real(wp_), intent(in) :: y(:, :) + real(wp_) :: rhs(size(y, dim=1), size(y, dim=2)) + end function rhs + end interface + contains subroutine gray_main(params, equil, plasma, limiter, results, error) - use const_and_precisions, only : zero, one, comp_tiny use polarization, only : ellipse_to_field use types, only : contour, table, wrap - use gray_params, only : gray_parameters, gray_results, EQ_VACUUM, ABSORP_OFF - use gray_equil, only : abstract_equil - use gray_plasma, only : abstract_plasma + use gray_params, only : gray_results, EQ_VACUUM use gray_project, only : ray_projector use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, & store_beam_shape, find_flux_surfaces, & flux_surfaces, kinetic_profiles, flux_averages, & ec_resonance, inputs_maps - use gray_errors, only : gray_error, is_critical, has_new_errors, & - print_err_raytracing, print_err_ecrh_cd - use beams, only : xgygcoeff, launchangles2n - use beamdata, only : pweight - use utils, only : vmaxmin - use multipass, only : initbeam, initmultipass, turnoffray, & - plasma_in, plasma_out, wall_out use logger, only : log_info, log_debug ! subroutine arguments @@ -37,86 +76,14 @@ contains integer(kind=gray_error), intent(out) :: error ! local variables - real(wp_), parameter :: taucr = 12._wp_, etaucr = exp(-taucr) - character, dimension(2), parameter :: mode=['O','X'] - - integer :: sox - real(wp_) :: ak0,bres,xgcn,xg,yg,rrm,zzm,alpha,didp,anpl,anpr,anprim,anprre - real(wp_) :: chipol,psipol,btot,psinv,dens,tekev,dersdst,derdnm - real(wp_) :: tau,pow,dids,ddr,ddi,taumn,taumx - real(wp_) :: rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip - real(wp_) :: rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx - real(wp_) :: pabs_beam,icd_beam,cpl_beam1,cpl_beam2,cpl_cbeam1,cpl_cbeam2 - - integer :: iox,nharm,nhf,iokhawa,index_rt, parent_index_rt - integer :: ip,ib,iopmin,ipar,child_index_rt - integer :: igrad_b,istop_pass,nbeam_pass - logical :: ins_pl,ins_wl,ent_pl,ext_pl,ent_wl,ext_wl,iboff - - ! Error handing - integer(kind=gray_error) :: prev_error ! from previous step - integer(kind=gray_error) :: curr_error ! current step - - ! i: integration step, jk: global ray index - integer :: i, jk - - ! buffer for formatting log messages character(256) :: msg - type(ray_projector) :: projector - - real(wp_), dimension(2) :: pabs_pass,icd_pass,cpl - real(wp_), dimension(3) :: xv,anv0,anv,bv,derxg - - associate (nray => params%raytracing%nray, & - nrayr => params%raytracing%nrayr, & - nrayth => params%raytracing%nrayth, & - nstep => params%raytracing%nstep, & - npass => params%raytracing%ipass, & - nbeam_tot => 2**(params%raytracing%ipass+1)-2, & - nbeam_max => 2**(params%raytracing%ipass)) - block - - ! ray variables - real(wp_), dimension(6, nray) :: yw, ypw - real(wp_), dimension(3, nray) :: gri - real(wp_), dimension(3, 3, nray) :: ggri - real(wp_), dimension(3, nrayth, nrayr) :: xc, du1 - real(wp_), dimension(nray) :: ccci0, p0jk - real(wp_), dimension(nray) :: tau0, alphaabs0 - real(wp_), dimension(nray) :: dids0 - integer, dimension(nray) :: iiv - real(wp_), dimension(nray, nstep) :: psjki, ppabs, ccci - complex(wp_), dimension(nray) :: ext, eyt - - ! multipass variables - logical, dimension(nray) :: iwait - integer, dimension(nray) :: iow, iop - logical, dimension(nray, nbeam_tot) :: iroff - real(wp_), dimension(6, nray, nbeam_max-2) :: yynext, yypnext - real(wp_), dimension(6, nray) :: yw0, ypw0 - real(wp_), dimension(nray, nbeam_tot) :: stnext - real(wp_), dimension(nray) :: stv, p0ray - real(wp_), dimension(nray, nbeam_tot) :: taus, cpls - real(wp_), dimension(nray) :: tau1, etau1, cpl1, lgcpl1 - real(wp_), dimension(0:nbeam_tot) :: psipv, chipv - - ! beam variables - real(wp_), dimension(params%output%nrho) :: jphi_beam, pins_beam - real(wp_), dimension(params%output%nrho) :: currins_beam, dpdv_beam, jcd_beam - real(wp_), dimension(2) :: mode_ellipse ! ψ, χ of the current polarisation mode - - ! ======== set environment BEGIN ======== - ! Compute X=ω/ω_ce and Y=(ω/ω_pe)² (with B=1) - call xgygcoeff(params%antenna%fghz, ak0, bres, xgcn) + type(beam_state) :: beam + type(beam_state), allocatable :: trace(:) ! Initialise the output tables call init_tables(params, results%tables) - ! Compute the initial cartesian wavevector (anv0) - ! from launch angles α,β and the position - call launchangles2n(params%antenna, anv0) - if (params%equilibrium%iequil /= EQ_VACUUM) then ! Initialise the output profiles call projector%init(params%output, equil) @@ -131,14 +98,13 @@ contains results%icd = 0 results%dpdv = 0 results%jcd = 0 - ! ========= set environment END ========= ! Pre-determinted tables results%tables%kinetic_profiles = kinetic_profiles(params, equil, plasma) results%tables%flux_averages = flux_averages(params, equil) results%tables%flux_surfaces = flux_surfaces(params, equil) - results%tables%ec_resonance = ec_resonance(params, equil, bres) - results%tables%inputs_maps = inputs_maps(params, equil, plasma, bres, xgcn) + results%tables%ec_resonance = ec_resonance(params, equil, beam%B_res) + results%tables%inputs_maps = inputs_maps(params, equil, plasma, beam%B_res, beam%X_norm) ! print ψ rational surfaces call find_flux_surfaces(qvals=[1.0_wp_, 1.5_wp_, 2.0_wp_], & @@ -152,553 +118,490 @@ contains 'α', params%antenna%alpha, 'β', params%antenna%beta call log_info(msg, mod='gray_core', proc='gray_main') - ! =========== main loop BEGIN =========== - call initmultipass(params, iroff, yynext, yypnext, yw0, ypw0, & - stnext, p0ray, taus, tau1, etau1, cpls, & - cpl1, lgcpl1, psipv, chipv) + call compute_initial_conds(params%raytracing, params%antenna, beam) + beam%mode = -1 + allocate(trace(0:params%raytracing%nstep)) + call propagate(params, beam, equil, plasma, trace, error) - sox=1 ! mode inverted for each beam - iox=2 ! start with O: sox=-1, iox=1 - - call pweight(params%raytracing, params%antenna%power, p0jk) - - ! Set original polarisation - psipv(0) = params%antenna%psi - chipv(0) = params%antenna%chi - call ellipse_to_field(psipv(0), chipv(0), ext, eyt) - - ! Error value for the whole simulation - error = 0 - - nbeam_pass=1 ! max n of beam per pass - index_rt=0 ! global beam index: 1,O 2,X 1st pass - ! | | | | - call log_debug('pass loop start', mod='gray_core', proc='gray_main') ! 3,O 4,X 5,O 6,X 2nd pass - main_loop: do ip=1,params%raytracing%ipass - write (msg, '("pass: ",g0)') ip - call log_info(msg, mod='gray_core', proc='gray_main') - - pabs_pass = 0 - icd_pass = 0 - istop_pass = 0 ! stop flag for current pass - nbeam_pass = 2*nbeam_pass ! max n of beams in current pass - - if(ip > 1) then - du1 = 0 - gri = 0 - ggri = 0 - if(ip == params%raytracing%ipass) cpl = [zero, zero] ! no successive passes - end if - - ! =========== beam loop BEGIN =========== - call log_debug('beam loop start', mod='gray_core', proc='gray_main') - - beam_loop: do ib=1,nbeam_pass - - sox = -sox ! invert mode - iox = 3-iox ! O-mode at ip=1,ib=1 - index_rt = index_rt +1 - child_index_rt = 2*index_rt + 1 ! * index_rt of O-mode child beam - parent_index_rt = ceiling(index_rt / 2.0_wp_) - 1 ! * index_rt of parent beam - - call initbeam(index_rt,iroff,iboff,iwait,stv,jphi_beam, & - pins_beam,currins_beam,dpdv_beam,jcd_beam) - - write(msg, '(" beam: ",g0," (",a1," mode)")') index_rt, mode(iox) - call log_info(msg, mod='gray_core', proc='gray_main') - - if(iboff) then ! no propagation for current beam - istop_pass = istop_pass +1 ! * +1 non propagating beam - call log_info(" beam is off", mod='gray_core', proc='gray_main') - cycle - end if - - psjki = 0 - ppabs = 0 - ccci = 0 - tau0 = 0 - alphaabs0 = 0 - dids0 = 0 - ccci0 = 0 - iiv = 1 - - if(ip == 1) then ! 1st pass - igrad_b = params%raytracing%igrad ! * input value, igrad_b=0 from 2nd pass - - tau1 = 0 ! * tau from previous passes - etau1 = 1 - cpl1 = 1 ! * coupling from previous passes - lgcpl1 = 0 - p0ray = p0jk ! * initial beam power - - call compute_initial_conds(params%raytracing, params%antenna, & ! * initial conditions - anv0, ak0, yw, ypw, stv, xc, du1, gri, ggri) - - do jk=1,params%raytracing%nray - ! Save the step "zero" data - call store_ray_data(params, equil, results%tables, & - i=0, jk=jk, s=stv(jk), P0=p0jk(jk), pos=yw(1:3,jk), & - psi_n=-one, B=zero, b_n=[zero,zero,zero], k0=ak0, & - N_pl=zero, N_pr=zero, N=yw(4:6,jk), N_pr_im=zero, & - n_e=zero, T_e=zero, alpha=zero, tau=zero, dIds=zero, & - nhm=0, nhf=0, iokhawa=0, index_rt=index_rt, & - lambda_r=zero, lambda_i=zero, X=zero, Y=zero, & - grad_X=[zero,zero,zero]) - - zzm = yw(3,jk)*0.01_wp_ - rrm = sqrt(yw(1,jk)*yw(1,jk)+yw(2,jk)*yw(2,jk))*0.01_wp_ - - if (limiter%contains(rrm, zzm)) then ! * start propagation in/outside vessel? - iow(jk) = 1 ! + inside - else - iow(jk) = 0 ! + outside - end if - end do - - else ! 2nd+ passes - ipar = (index_rt+1)/2-1 ! * parent beam index - yw = yynext(:,:,ipar) ! * starting coordinates from - ypw = yypnext(:,:,ipar) ! parent beam last step - stv = stnext(:,ipar) ! * starting step from parent beam last step - iow = 1 ! * start propagation inside vessel - - tau1 = taus(:,index_rt) ! * tau from previous passes - etau1 = exp(-tau1) - cpl1 = cpls(:,index_rt) ! * coupling from previous passes - lgcpl1 = -log(cpl1) - p0ray = p0jk * etau1 * cpl1 ! * initial beam power - end if - iop = 0 ! start propagation outside plasma - - if(params%raytracing%nray>1 .and. all(.not.iwait)) & ! iproj=0 ==> nfilp=8 - call store_beam_shape(params%raytracing, results%tables, stv, yw) - - ! ======= propagation loop BEGIN ======= - call log_debug(' propagation loop start', mod='gray_core', proc='gray_main') - - ! error value for the this step - curr_error = 0 - - propagation_loop: do i=1,params%raytracing%nstep - ! rotate current and previous step errors - prev_error = curr_error - curr_error = 0 - - ! advance one step with "frozen" grad(S_I) - do jk=1,params%raytracing%nray - if(iwait(jk)) cycle ! jk ray is waiting for next pass - stv(jk) = stv(jk) + params%raytracing%dst ! current ray step - call rkstep(params, equil, plasma, sox, bres, xgcn, yw(:,jk), & - ypw(:,jk), gri(:,jk), ggri(:,:,jk), igrad_b) - end do - ! update position and grad - 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 =========== - rays_loop: do jk=1,params%raytracing%nray - if(iwait(jk)) cycle ! jk ray is waiting for next pass - - ! compute derivatives with updated gradient and local plasma values - xv = yw(1:3,jk) - anv = yw(4:6,jk) - call ywppla_upd(params, equil, plasma, & - 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) - - ! check entrance/exit plasma/wall - zzm = xv(3)*0.01_wp_ - rrm = sqrt(xv(1)*xv(1)+xv(2)*xv(2))*0.01_wp_ - - ins_pl = (psinv>=zero .and. psinv continue current pass - - if(cpl(iox) < etaucr) then ! + IF low coupled power for current mode => de-activate derived rays - call turnoffray(jk,ip+1,npass,2*ib+2-iox,iroff) - iwait(jk) = .true. ! . stop advancement and H&CD computation for current ray - if(cpl(iox).le.comp_tiny) cpl(iox)=etaucr - else ! + ELSE assign coupled power to current ray - p0ray(jk) = p0ray(jk)*cpl(iox) - end if - cpls(jk,index_rt) = cpl(iox) - - if(jk == 1) then - write (msg,'(" 1st pass - central ray (",a1,"-mode) c=",g0.4)') & - mode(iox), cpl(iox) - call log_info(msg, mod='gray_core', proc='gray_main') - write (msg,'(" polarisation: ψ=",g0.5,"°, χ=",g0.5,"°")') & - psipol, chipol - call log_debug(msg, mod='gray_core', proc='gray_main') - psipv(index_rt) = psipol ! + polarization angles at plasma boundary for central ray - chipv(index_rt) = chipol - end if - - else if(iop(jk) > 2) then ! * 2nd entrance on 1st pass / entrance on 2nd+ pass => end of current pass for ray jk - igrad_b = 0 ! + switch to ray-tracing - iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray - - if(ip < params%raytracing%ipass) then ! + not last pass - yynext(:,jk,index_rt) = yw0(:,jk) ! . copy starting coordinates - yypnext(:,jk,index_rt) = ypw0(:,jk) ! for next pass from last step - stnext(jk,index_rt) = stv(jk) - params%raytracing%dst ! . starting step for next pass = last step - - if(cpl(1) < etaucr) then ! . low coupled power for O-mode => de-activate derived rays - call turnoffray(jk,ip+1,npass,2*ib-1,iroff) - if(cpl(1).le.comp_tiny) cpl(1)=etaucr - end if - if(cpl(2) < etaucr) then ! . low coupled power for X-mode => de-activate derived rays - call turnoffray(jk,ip+1,npass,2*ib,iroff) - if(cpl(2).le.comp_tiny) cpl(2)=etaucr - end if - - taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass - cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass - cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass - - if(jk == 1) then ! . polarization angles at plasma boundary for central ray - psipv(child_index_rt:child_index_rt+1) = psipol - chipv(child_index_rt:child_index_rt+1) = chipol - end if - else ! * 1st entrance on 2nd+ pass (ray hasn't entered in plasma since end of previous pass) => continue current pass - cpl = [zero, zero] - end if - end if - - else if(ext_pl) then ! ray exits plasma - write (msg, '(" ray ", g0, " left plasma")') jk - call log_debug(msg, mod='gray_core', proc='gray_main') - call plasma_out(jk, equil, xv, anv, bres, sox, iop, ext, eyt) - end if - - if(ent_wl) then ! ray enters vessel - iow(jk)=iow(jk)+1 ! * out->in - - else if(ext_wl) then ! ray exit vessel - call wall_out(jk, equil, limiter, ins_pl, xv, anv, & - params%raytracing%dst, bres, sox, psipol, chipol, & - iow, iop, ext, eyt) - yw(:,jk) = [xv, anv] ! * updated coordinates (reflected) - igrad_b = 0 ! * switch to ray-tracing - - call ywppla_upd(params, equil, plasma, & - 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 - - if(jk == 1 .and. ip == 1) then ! * 1st pass, polarization angles at reflection for central ray - psipv(index_rt) = psipol - chipv(index_rt) = chipol - end if - - if(ins_pl) then ! * plasma-wall overlapping => wall+plasma crossing => end of current pass - iwait(jk) = .true. ! + stop advancement and H&CD computation for current ray - - if(ip < params%raytracing%ipass) then ! + not last pass - yynext(:,jk,index_rt) = [xv, anv] ! . starting coordinates - yypnext(:,jk,index_rt) = ypw(:,jk) ! for next pass = reflection point - stnext(jk,index_rt) = stv(jk) ! . starting step for next pass = step after reflection - - call plasma_in(jk, equil, xv, anv, bres, sox, cpl, & ! . ray re-enters plasma after reflection - psipol, chipol, iop, ext, eyt, perfect=.false.) - - if(cpl(1) < etaucr) then ! . low coupled power for O-mode? => de-activate derived rays - call turnoffray(jk,ip+1,npass,2*ib-1,iroff) - if(cpl(1).le.comp_tiny) cpl(1)=etaucr - end if - if(cpl(2) < etaucr) then ! . low coupled power for X-mode? => de-activate derived rays - call turnoffray(jk,ip+1,npass,2*ib,iroff) - if(cpl(2).le.comp_tiny) cpl(2)=etaucr - end if - - taus(jk,child_index_rt:child_index_rt+1) = tau1(jk) + tau0(jk) ! . starting tau for next O-mode pass - cpls(jk,child_index_rt) = cpl1(jk) * cpl(1) ! . cumulative coupling for next O-mode pass - cpls(jk,child_index_rt+1) = cpl1(jk) * cpl(2) ! . cumulative coupling for next X-mode pass - - if(jk == 1) then ! + polarization angles at plasma boundary for central ray - psipv(child_index_rt:child_index_rt+1) = psipol - chipv(child_index_rt:child_index_rt+1) = chipol - end if - end if - end if - end if - - iopmin = min(iopmin,iop(jk)) - if(ip < params%raytracing%ipass) then ! not last pass - yw0(:,jk) = yw(:,jk) ! * store current coordinates in case - ypw0(:,jk) = ypw(:,jk) ! current pass ends on next step - end if - - ! Compute ECRH&CD if (inside plasma & power available>0 & ray still active) - ! Note: power check is τ + τ₀ + log(coupling O/X) < τ_critic - - if (params%ecrh_cd%iwarm /= ABSORP_OFF .and. ins_pl .and. & - (tau1(jk)+tau0(jk)+lgcpl1(jk))<=taucr .and. .not.iwait(jk)) then ! H&CD computation check - tekev = plasma%temp(psinv) - block - complex(wp_) :: Npr_warm - call alpha_effj(params%ecrh_cd, equil, plasma, & - psinv, xg, yg, dens, tekev, ak0, bres, & - derdnm, anpl, anpr, sox, Npr_warm, alpha, & - didp, nharm, nhf, iokhawa, curr_error) - anprre = Npr_warm%re - anprim = Npr_warm%im - - if (has_new_errors(prev_error, curr_error)) then - call print_err_ecrh_cd(curr_error, i, jk, Npr_warm, alpha) - end if - end block - else - tekev=0 - alpha=0 - didp=0 - anprim=0 - anprre=anpr - nharm=0 - nhf=0 - iokhawa=0 - end if - if(nharm>0) iiv(jk)=i - - psjki(jk,i) = psinv - - ! Computation of the ray τ, dP/ds, P(s), dI/ds, I(s) - - ! optical depth: τ = ∫α(s)ds using the trapezoid rule - tau = tau0(jk) + 0.5_wp_*(alphaabs0(jk) + alpha) * dersdst * params%raytracing%dst - - pow = p0ray(jk) * exp(-tau) ! residual power: P = P₀exp(-τ) - ppabs(jk,i) = p0ray(jk) - pow ! absorbed power: P_abs = P₀ - P - dids = didp * pow * alpha ! current driven: dI/ds = dI/dP⋅dP/ds = dI/dP⋅P⋅α - - ! current: I = ∫dI/ds⋅ds using the trapezoid rule - ccci(jk,i) = ccci0(jk) + 0.5_wp_*(dids0(jk) + dids) * dersdst * params%raytracing%dst - - tau0(jk) = tau - alphaabs0(jk) = alpha - dids0(jk) = dids - ccci0(jk) = ccci(jk,i) - - if(iwait(jk)) then ! copy values from last pass for inactive ray - ppabs(jk,i:params%raytracing%nstep) = ppabs(jk,i-1) - ccci(jk,i:params%raytracing%nstep) = ccci(jk,i-1) - psjki(jk,i:params%raytracing%nstep) = psjki(jk,i-1) - else - call store_ray_data(params, equil, results%tables, & - i, jk, stv(jk), p0ray(jk), xv, psinv, btot, bv, ak0, & - anpl, anpr, anv, anprim, dens, tekev, alpha, tau, dids, & - nharm, nhf, iokhawa, index_rt, ddr, ddi, xg, yg, derxg) ! p0ray/etau1 [dids normalization] = fraction of p0 coupled to this ray (not including absorption from previous passes) - end if - - ! Accumulate errors from the latest step - error = ior(error, curr_error) - - end do rays_loop - ! ============ rays loop END ============ - - 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 - end do - end if - - ! print ray positions for j=nrayr in local reference system - if(mod(i,params%output%istpr) == 0) then - if(params%raytracing%nray > 1 .and. all(.not.iwait)) & - call store_beam_shape(params%raytracing, results%tables, stv, yw) - end if - - ! test whether further trajectory integration is unnecessary - call vmaxmin(tau1+tau0+lgcpl1, taumn, taumx) ! test on tau + coupling - - if (is_critical(curr_error)) then ! stop propagation for current beam - istop_pass = istop_pass +1 ! * +1 non propagating beam - if(ip < params%raytracing%ipass) call turnoffray(0,ip,npass,ib,iroff) ! * de-activate derived beams - exit - else if(all(iwait)) then ! all rays in current beam are waiting for next pass => do not increase istop_pass - exit - end if - end do propagation_loop - call log_debug(' propagation loop end', mod='gray_core', proc='gray_main') - ! ======== propagation loop END ======== - - ! print all ray positions in local reference system - if(params%raytracing%nray > 1 .and. all(.not.iwait)) & - call store_beam_shape(params%raytracing, results%tables, & - stv, yw, full=.true.) - - ! =========== post-proc BEGIN =========== - ! compute total absorbed power and driven current for current beam - if(i>params%raytracing%nstep) i=params%raytracing%nstep - pabs_beam = sum(ppabs(:,i)) - icd_beam = sum(ccci(:,i)) - call vmaxmin(tau0, taumn, taumx) ! taumn,taumx for print - - if (params%equilibrium%iequil /= EQ_VACUUM) then - ! compute power and current density profiles for all rays - call projector%project(psjki, ppabs, ccci, iiv, dpdv_beam, & - jphi_beam, jcd_beam, pins_beam, currins_beam) - end if - - pabs_pass(iox) = pabs_pass(iox) + pabs_beam ! 0D results for current pass, sum on O/X mode beams - icd_pass(iox) = icd_pass(iox) + icd_beam - - if(ip < params%raytracing%ipass .and. iopmin > 2) then ! not last pass AND at least one ray re-entered plasma - cpl_beam1 = sum(& - p0ray * exp(-tau0) * cpls(:,child_index_rt)/cpl1, MASK=iop > 2) / & - sum(p0ray * exp(-tau0), MASK=iop > 2) ! * average O-mode coupling for next beam (on active rays) - cpl_beam2 = 1 - cpl_beam1 ! * average X-mode coupling for next beam - - if(iop(1) > 2) then ! * central ray O/X-mode coupling for next beam - cpl_cbeam1 = cpls(1,child_index_rt)/cpl1(1) - cpl_cbeam2 = 1 - cpl_cbeam1 - end if - else ! last pass OR no ray re-entered plasma - cpl_beam1 = 0 - cpl_beam2 = 0 - end if - - ! print final results for pass on screen - call log_info(' partial results:', mod='gray_core', proc='gray_main') - write(msg, '(" final step:", x,a,g0, x,a,g0.4)') 'n=', i, '(s, ct, Sr)=', stv(1) - call log_info(msg, mod='gray_core', proc='gray_main') - - write(msg, '(3x,a,2(x,a,"=",g0.4))') 'optical depth:', 'τ_min', taumn, 'τ_max', taumx - call log_info(msg, mod='gray_core', proc='gray_main') - - write(msg, '(3x,a,g0.3," MW")') 'absoption: P=', pabs_beam - call log_info(msg, mod='gray_core', proc='gray_main') - - write(msg, '(3x,a,g0.3," kA")') 'current drive: I=', icd_beam*1.0e3_wp_ - call log_info(msg, mod='gray_core', proc='gray_main') - - if(ip < params%raytracing%ipass) then - write (msg,'(3x,a,(g0.4,", ",g0.4))') & ! average coupling for next O/X beams (=0 if no ray re-entered plasma) - 'next couplings [O,X mode]: c=', cpl_beam1, cpl_beam2 - call log_info(msg, mod='gray_core', proc='gray_main') - if(iop(1) > 2) then - write(msg, '(3x,a,(g0.4,", ",g0.4))') & - 'coupling [ctr ray, O/X]:', cpl_cbeam1, cpl_cbeam2 ! central ray coupling for next O/X beams - call log_info(msg, mod='gray_core', proc='gray_main') - end if - end if - - if (params%equilibrium%iequil /= EQ_VACUUM) then - call store_ec_profiles( & - results%tables%ec_profiles, projector%rho_p, projector%rho_t, & - jphi_beam, jcd_beam, dpdv_beam, currins_beam, pins_beam, index_rt) - - call projector%statistics(equil, dpdv_beam, jphi_beam, & - rhotpav, drhotpav, rhotjava, drhotjava, dpdvp, jphip, & - rhotp, drhotp, rhotj, drhotj, dpdvmx, jphimx, ratjamx, ratjbmx) ! *compute profiles width for current beam - - if (results%tables%summary%active) then - call results%tables%summary%append([ & - wrap(icd_beam), wrap(pabs_beam), wrap(jphip), wrap(dpdvp), & - wrap(rhotj), wrap(rhotjava), wrap(rhotp), wrap(rhotpav), & - wrap(drhotjava), wrap(drhotpav), wrap(ratjamx), wrap(ratjbmx), & - wrap(stv(1)), wrap(mode_ellipse(1)), wrap(mode_ellipse(2)), & - wrap(index_rt), wrap(jphimx), wrap(dpdvmx), wrap(drhotj), & - wrap(drhotp), wrap(sum(p0ray)), wrap(cpl_beam1), wrap(cpl_beam2)]) ! *print 0D results for current beam - end if - - end if - - ! ============ post-proc END ============ - - end do beam_loop - call log_debug('beam loop end', mod='gray_core', proc='gray_main') - ! ============ beam loop END ============ - - ! ======= cumulative prints BEGIN ======= - results%pabs = results%pabs + sum(pabs_pass) ! *final results (O+X) [gray_main output] - results%icd = results%icd + sum(icd_pass) - - ! print final results for pass on screen - call log_info(' comulative results:', mod='gray_core', proc='gray_main') - - write(msg, '(" absoption [O,X mode] P=",g0.4,", ",g0.4," MW")') & - pabs_pass(1), pabs_pass(2) - call log_info(msg, mod='gray_core', proc='gray_main') - - write(msg, '(" current drive [O,X mode] I=",g0.4,", ",g0.4," kA")') & - icd_pass(1)*1.0e3_wp_, icd_pass(2)*1.0e3_wp_ - call log_info(msg, mod='gray_core', proc='gray_main') - ! ======== cumulative prints END ======== - - if(istop_pass == nbeam_pass) exit ! no active beams - end do main_loop - call log_debug('pass loop end', mod='gray_core', proc='gray_main') - ! ============ main loop END ============ - - end block - end associate end subroutine gray_main - subroutine compute_initial_conds(rtx, beam, N_c, k0, y, yp, dist, & - pos, grad_u, grad, hess) - ! Computes the initial conditions for tracing a beam + subroutine propagate(params, beam, equil, plasma, trace, error) + ! Propagates a beam producing a trace + + ! subroutine arguments + type(gray_parameters), intent(in) :: params ! gray parameters + type(beam_state), intent(inout) :: beam ! initial beam state + class(abstract_equil), intent(in) :: equil ! equilibrium object + class(abstract_plasma), intent(in) :: plasma ! plasma object + type(beam_state), intent(out) :: trace(params%raytracing%nstep) ! beam trajectory + integer(kind=gray_error), intent(out) :: error + + ! local variables + integer :: step, var + real(wp_) :: h + real(wp_) :: alpha_0(beam%nray), alpha_1(beam%nray), dIdP(beam%nray) + real(wp_) :: dIds_0(beam%nray), dIds_1(beam%nray) + + h = params%raytracing%dst ! step size + var = params%raytracing%idst ! integration variable + + error = 0 + + alpha_0 = 0 + dIds_0 = 0 + + do step = 1, params%raytracing%nstep + ! Integrate by one step + call integrate(f, h, beam%y) + + ! Update travelled distance + beam%step = step + beam%dist = beam%dist + h + + ! Update eikonal gradient + if (beam%beamtracing) call update_gradient(beam, error) + + ! Update local quantities + call update_locals(equil, plasma, var, beam) + + ! Compute absorption and current drive + call compute_absorption(params%ecrh_cd, equil, plasma, beam, & + alpha_1, dIdP, error) + + ! optical depth: τ = ∫α(s) ds/dσ dσ + beam%tau = beam%tau + (alpha_1 + alpha_0)/2 * beam%dsdvar * h + + ! driven current: dI/ds = dI/dP⋅dP/ds = dI/dP α exp(-τ) + dIds_1 = dIdP * alpha_1 * exp(-beam%tau) + + ! total driven current, I_cd = ∫dI/ds ds + beam%I_cd = beam%I_cd + (dIds_1 + dIds_0)/2 * beam%dsdvar * h + + if (dIds_1(1) > 0) & + print *, "s:", beam%dist, "ψ:", beam%psi_n(1), "τ", beam%tau(1) + + ! Save current state + dIds_0 = dIds_1 + alpha_0 = alpha_1 + trace(step) = beam + end do + + contains + + function f(y) result(dy) + ! function argument + real(wp_), intent(in) :: y(:, :) + real(wp_) :: dy(size(y, dim=1), size(y, dim=2)) + + ! local variables + integer :: ray + real(wp_) :: n_e, B_tot, X, Yg + real(wp_) :: b(3), grad_b(3,3), grad_X(3), grad_Y(3) + + do concurrent (ray = 1 : beam%nray) + call plas_deriv(equil, plasma, beam%y(1:3,ray), beam%B_res, beam%X_norm, & + n_e, B_tot, b, grad_b, X, Yg, grad_X, grad_Y) + call disp_deriv(var, beam%y(4:6,ray), beam%mode(ray), X, Yg, grad_X, grad_Y, & + b, grad_b, beam%grad_S(:,ray), beam%hess_S(:,:,ray), & + beam%beamtracing, dy(:,ray)) + end do + + end function f + + end subroutine propagate + + + subroutine integrate(f, h, y) + ! Runge-Kutta integrator + + ! subroutine arguments + procedure(rhs) :: f + real(wp_), intent(in) :: h + real(wp_), intent(inout) :: y(:,:) + + ! local variables + real(wp_), dimension(size(y, dim=1), size(y, dim=2)) :: k1, k2, k3, k4 + + k1 = f(y) + k2 = f(y + k1*h/2) + k3 = f(y + k2*h/2) + k4 = f(y + k3*h) + y = y + h * (k1/6 + k2/3 + k3/3 + k4/6) + end subroutine integrate + + + subroutine update_locals(equil, plasma, var, beam) + ! Computes the values of X, Y, N∥, N⊥, ds/dσ, |∂Λ/∂N| of each ray. + ! + ! Note: these require solving the dispersion relation again. + use gray_params, only : step_enum + + ! subroutine argument + class(abstract_equil), intent(in) :: equil ! equilibrium object + class(abstract_plasma), intent(in) :: plasma ! plasma object + integer(kind(step_enum)), intent(in) :: var ! integration variable + type(beam_state(*,*,*)), intent(inout) :: beam ! current beam state + + ! local variables + integer :: ray + real(wp_) :: n_e, B_tot, b(3) + real(wp_) :: grad_b(3,3), grad_X(3), grad_Y(3) + + do concurrent (ray = 1 : beam%nray) + call plas_deriv(equil, plasma, beam%y(1:3,ray), beam%B_res, beam%X_norm, & + n_e, B_tot, b, grad_b, beam%Xg(ray), beam%Yg(ray), & + grad_X, grad_Y, beam%psi_n(ray)) + call disp_deriv(var, beam%y(4:6,ray), beam%mode(ray), beam%Xg(ray), & + beam%Yg(ray), grad_X, grad_Y, b, grad_b, & + beam%grad_S(:,ray), beam%hess_S(:,:,ray), & + beam%beamtracing, beam%yp(:,ray), & + anpl_=beam%N_pl(ray), anpr=beam%N_pr(ray), & + dersdst=beam%dsdvar(ray), derdnm_=beam%dlambdadN(ray)) + end do + end subroutine update_locals + + + subroutine update_gradient(beam, error) + ! Computes ∇u, ∇S_I, HS_I after a new integration step + use gray_errors, only : unstable_beam, print_err_raytracing + + ! subroutine parameters + type(beam_state(*,*,*)), intent(inout) :: beam + integer(kind=gray_error), intent(inout) :: error + + ! local variables + real(wp_), dimension(3, beam%nrayth, beam%nrayr) :: xco, grad_uo + integer :: jk, j, jm, jp, k, km, kp + real(wp_) :: ux, uxx, uxy, uxz, uy, uyy, uyz, uz, uzz + real(wp_) :: dr, dfuu, dffiu + real(wp_), dimension(3) :: dxv1, dxv2, dxv3, dgu + real(wp_), dimension(3,3) :: dgg, dff + integer(kind=gray_error) :: new_error + + new_error = 0 + + ! update position and beam%grad_u vectors + xco = beam%pos + grad_uo = beam%grad_u + + jk = 1 + do j=1,beam%nrayr + do k=1,beam%nrayth + if(j>1) jk=jk+1 + beam%pos(1:3,k,j)=beam%y(1:3,jk) + end do + end do + + ! compute grad u1 for central ray + j = 1 + jp = 2 + do k=1,beam%nrayth + if(k == 1) then + km = beam%nrayth + else + km = k-1 + end if + if(k == beam%nrayth) then + kp = 1 + else + kp = k+1 + end if + dxv1 = beam%pos(:,k ,jp) - beam%pos(:,k ,j) + dxv2 = beam%pos(:,kp,jp) - beam%pos(:,km,jp) + dxv3 = beam%pos(:,k ,j) - xco(:,k ,j) + call solg0(dxv1,dxv2,dxv3,dgu) + beam%grad_u(:,k,j) = dgu + end do + beam%grad_S(:,1) = 0 + + ! compute grad u1 and grad(S_I) for all the other rays + if (beam%nrayr > 1) then + dr = beam%rwmax / (beam%nrayr - 1) + else + dr = 1 + end if + dfuu=2*dr**2/beam%k0 + jm=1 + j=2 + k=0 + dffiu = dfuu + do jk=2,beam%nray + k=k+1 + if(k > beam%nrayth) then + jm = j + j = j+1 + k = 1 + dffiu = dfuu*jm + end if + kp = k+1 + km = k-1 + if (k == 1) then + km=beam%nrayth + else if (k == beam%nrayth) then + kp=1 + end if + dxv1 = beam%pos(:,k ,j) - beam%pos(:,k ,jm) + dxv2 = beam%pos(:,kp,j) - beam%pos(:,km,j) + dxv3 = beam%pos(:,k ,j) - xco(:,k ,j) + + call solg0(dxv1, dxv2, dxv3, dgu) + + ! Check for numerical instability in the computation of ∇S_I. + ! This can happen when the distances between the rays and the + ! integration step size differ wildly; in such cases the linear + ! system Δx⋅∇u = Δu becomes bad conditioned. + ! + ! Source: https://discourse.julialang.org/t/what-is-the-best-way-to-solve-a-ill-conditioned-linear-problem/95894/10 + block + use const_and_precisions, only : comp_eps + use utils, only : singvals + real(wp_) :: A(3,3), x(3), b(3), sigma(3) + real(wp_) :: R, normA, back_error, forw_error + A = reshape([dxv1, dxv2, dxv3], shape=[3,3], order=[2,1]) + sigma = singvals(A) + b = [1, 0, 0] + x = dgu + R = maxval(sigma) / minval(sigma) ! condition number + normA = maxval(sigma) ! L₂ norm of A + back_error = norm2(matmul(A, x) - b) / (normA * norm2(b)) + forw_error = back_error * R * norm2(x) + ! TODO: figure out how to tie this threshold + ! to the derivatives of the dispersion relation + if (forw_error / comp_eps > 5000) then + new_error = raise_error(new_error, unstable_beam) + end if + end block + + if (has_new_errors(error, new_error)) then + call print_err_raytracing(new_error, beam%step, jk) + end if + + beam%grad_u(:,k,j) = dgu + beam%grad_S(:,jk) = dgu(:)*dffiu + end do + + ! compute derivatives of grad u and grad(S_I) for rays jk>1 + beam%hess_S(:,:,1) = 0 + jm=1 + j=2 + k=0 + dffiu = dfuu + do jk=2,beam%nray + k=k+1 + if(k > beam%nrayth) then + jm=j + j=j+1 + k=1 + dffiu = dfuu*jm + end if + kp=k+1 + km=k-1 + if (k == 1) then + km=beam%nrayth + else if (k == beam%nrayth) then + kp=1 + end if + dxv1 = beam%pos(:,k ,j) - beam%pos(:,k ,jm) + dxv2 = beam%pos(:,kp,j) - beam%pos(:,km,j) + dxv3 = beam%pos(:,k ,j) - xco(:,k ,j) + dff(:,1) = beam%grad_u(:,k ,j) - beam%grad_u(:,k ,jm) + dff(:,2) = beam%grad_u(:,kp,j) - beam%grad_u(:,km,j) + dff(:,3) = beam%grad_u(:,k ,j) - grad_uo(:,k ,j) + call solg3(dxv1, dxv2, dxv3, dff, dgg) + + ! derivatives of u + ux = beam%grad_u(1,k,j) + uy = beam%grad_u(2,k,j) + uz = beam%grad_u(3,k,j) + uxx = dgg(1,1) + uyy = dgg(2,2) + uzz = dgg(3,3) + uxy = (dgg(1,2) + dgg(2,1))/2 + uxz = (dgg(1,3) + dgg(3,1))/2 + uyz = (dgg(2,3) + dgg(3,2))/2 + + ! Hessian of S_I + beam%hess_S(1,1,jk) = dfuu*ux*ux + dffiu*uxx + beam%hess_S(2,1,jk) = dfuu*ux*uy + dffiu*uxy + beam%hess_S(3,1,jk) = dfuu*uy*uz + dffiu*uyz + beam%hess_S(1,2,jk) = dfuu*uy*uz + dffiu*uyz + beam%hess_S(2,2,jk) = dfuu*uy*uy + dffiu*uyy + beam%hess_S(3,2,jk) = dfuu*uy*uz + dffiu*uyz + beam%hess_S(1,3,jk) = dfuu*ux*uz + dffiu*uxz + beam%hess_S(2,3,jk) = dfuu*uy*uz + dffiu*uyz + beam%hess_S(3,3,jk) = dfuu*uz*uz + dffiu*uzz + end do + + error = ior(error, new_error) + end subroutine update_gradient + + + subroutine compute_absorption(params, equil, plasma, beam, alpha, dIdP, error) + ! Computes the absorption coefficient and effective current + + use const_and_precisions, only : pi, mc2=>mc2_ + use gray_params, only : ecrh_cd_parameters, CD_OFF, CD_COHEN, & + CD_NO_TRAP, CD_NEOCLASSIC, ABSORP_OFF + use dispersion, only : harmnumber, colddisp, warmdisp + use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl + use gray_errors, only : negative_absorption, print_err_ecrh_cd + + ! subroutine arguments + type(ecrh_cd_parameters), intent(in) :: params ! ECRH & CD parameters + class(abstract_equil), intent(in) :: equil ! equilibrium object + class(abstract_plasma), intent(in) :: plasma ! plasma object + type(beam_state(*,*,*)), intent(in) :: beam ! beam state + real(wp_), intent(out) :: alpha(beam%nray) ! absorption coefficient + real(wp_), intent(out) :: dIdP(beam%nray) ! current driver efficiency + integer(kind(gray_error)), intent(inout) :: error + + ! local variables + real(wp_) :: psi_n, n_e, T_e, X, Y, mu, N_pl + real(wp_) :: N_pr_cold + complex(wp_) :: N_pr + integer(kind(gray_error)) :: new_error + + integer :: ray, nhmin, nhmax + integer :: nlarmor, ithn, ierrcd, iokhawa + real(wp_) :: cst2, Z_eff, B, B_avg, B_max, B_min, f_c, R_J + real(wp_), dimension(:), allocatable :: eccdpar(:) + real(wp_) :: effjcd, effjcdav + complex(wp_) :: e(3) + + new_error = 0 + + alpha = 0 + dIdp = 0 + + do ray = 1, beam%nray + ! Skip when power is too low + if (beam%tau(ray) > 12) cycle + + ! Skip when the temperature is zero (no plasma) + psi_n = beam%psi_n(ray) + T_e = plasma%temp(psi_n) + if (T_e <= 0) cycle + mu = mc2 / T_e + + ! Skip when the lowest resonant harmonic number is zero + X = beam%Xg(ray) + Y = beam%Yg(ray) + N_pl = beam%N_pl(ray) + call harmnumber(Y, mu, N_pl**2, params%iwarm == 1, nhmin, nhmax) + if (nhmin <= 0) cycle + + ! Solve the full dispersion only when needed + N_pr_cold = beam%N_pr(ray) + if (params%iwarm /= ABSORP_OFF .or. params%ieccd /= CD_OFF) then + nlarmor = max(params%ilarm, nhmax) + call warmdisp(X, Y, mu, N_pl, N_pr_cold, beam%mode(ray), & + new_error, N_pr, e, & + model=params%iwarm, nlarmor=nlarmor, & + max_iters=abs(params%imx), & + fallback=params%imx < 0) + end if + + ! Compute α from the solution of the dispersion relation + ! The absoption coefficient is defined as + ! + ! α = 2 Im(k̅)⋅s̅ + ! + ! where s̅ = v̅_g/|v_g|, the direction of the energy flow. + ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion + ! relation, we have that + ! + ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| + ! = [2N̅ - ∂(N²s)/∂N∥ b̅ + ! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅| + ! + ! Assuming Im(k∥)=0: + ! + ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| + ! + block + real(wp_) :: k_im + k_im = beam%k0 * N_pr%im ! imaginary part of k⊥ + alpha(ray) = 4 * k_im*N_pr_cold / beam%dlambdadN(ray) + end block + + if (alpha(ray) < 0) new_error = raise_error(new_error, negative_absorption) + + if (has_new_errors(error, new_error)) then + call print_err_ecrh_cd(new_error, beam%step, ray, N_pr, alpha(ray)) + end if + + if (has_error(error, negative_absorption)) cycle + + ! Current drive computation + if (params%ieccd == CD_OFF) cycle + + Z_eff = plasma%zeff(psi_n) + ithn = 1 + if (nlarmor > nhmin) ithn = 2 + call equil%flux_average(rho_p=sqrt(psi_n), R_J=R_J, B_avg=B_avg, & + B_min=B_min, B_max=B_max, f_c=f_c) + B = Y * beam%B_res + n_e = X * beam%X_norm + + select case(params%ieccd) + case (CD_COHEN) + call setcdcoeff(Z_eff, B/B_min, B/B_max, cst2, eccdpar) + call eccdeff(Y, N_pl, N_pr%re, n_e, mu, e, nhmin, nhmax, & + ithn, cst2, fjch, eccdpar, effjcd, iokhawa, ierrcd) + + case (CD_NO_TRAP) + call setcdcoeff(Z_eff, cst2, eccdpar) + call eccdeff(Y, N_pl, N_pr%re, n_e, mu, e, nhmin, nhmax, & + ithn, cst2, fjch0, eccdpar, effjcd, iokhawa, ierrcd) + + case (CD_NEOCLASSIC) + call setcdcoeff(Z_eff, B/B_max, f_c, mu, sqrt(psi_n), & + equil%spline_cd_eff, cst2, eccdpar) + call eccdeff(Y, N_pl, N_pr%re, n_e, mu, e, nhmin, nhmax, & + ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd) + end select + error = error + ierrcd + + ! current drive efficiency R* = / [A⋅m/W] + effjcdav = (B_avg / B_min) * effjcd + dIdp(ray) = equil%sgn_bphi * effjcdav / (2*pi*R_J) + + end do + + error = ior(error, new_error) + + end subroutine compute_absorption + + + subroutine compute_initial_conds(rtx, ant, beam) + ! Initialises the beam and computes: + ! - beam%y, (x̅, N̅) for each ray + ! - beam%yp, (dx̅/dσ, dN̅/dσ) for each ray + ! - beam%dist, distance travelled in σ for each ray + ! - beam%pos, ray positions (used in update_gradient) + ! - beam%grad_u, variations Δu (used in update_gradient) + ! - beam%grad_S, ∇S_I, gradient of imaginary eikonal + ! - beam%hess_S, HS_I, Hessian of imaginary eikonal use const_and_precisions, only : zero, one, pi, im, degree use gray_params, only : raytracing_parameters, antenna_parameters use gray_params, only : STEP_ARCLEN, STEP_TIME, STEP_PHASE use utils, only : diag, rotate - use beams, only : gaussian_eikonal + use beams, only : gaussian_eikonal, xgygcoeff, launchangles2n use beamdata, only : fold_indices, tokamak2beam ! subroutine arguments - - ! Inputs - type(raytracing_parameters), intent(in) :: rtx ! simulation parameters - type(antenna_parameters), intent(in) :: beam ! beam parameters - real(wp_), intent(in) :: N_c(3) ! vacuum refractive index - real(wp_), intent(in) :: k0 ! vacuum wavenumber - - ! Outputs - real(wp_), intent(out) :: y(6, rtx%nray) ! (x̅, N̅) for each ray - real(wp_), intent(out) :: yp(6, rtx%nray) ! (dx̅/dσ, dN̅/dσ) for each ray - real(wp_), intent(out) :: dist(rtx%nrayth * rtx%nrayr) ! distance travelled in σ for each ray - real(wp_), intent(out) :: pos(3, rtx%nrayth, rtx%nrayr) ! ray positions (used in gradi_upd) - real(wp_), intent(out) :: grad_u(3, rtx%nrayth, rtx%nrayr) ! variations Δu (used in gradi_upd) - real(wp_), intent(out) :: grad(3, rtx%nray) ! ∇S_I, gradient of imaginary eikonal - real(wp_), intent(out) :: hess(3, 3, rtx%nray) ! HS_I, Hessian of imaginary eikonal + type(raytracing_parameters), intent(in) :: rtx ! simulation parameters + type(antenna_parameters), intent(in) :: ant ! beam launcher parameters + type(beam_state), intent(out) :: beam ! initial beam state ! local variables real(wp_) :: x0_c(3), x0(3) ! beam centre, current ray position real(wp_) :: xi(2) ! position in the principal axes basis (ξ, η) integer :: j, k, jk ! ray indices, j: radial, k: angular real(wp_) :: dr, da ! ray grid steps, dρ: radial, dα: angular - real(wp_) :: N(3) ! wavevector + real(wp_) :: N_c(3), N(3) ! wavevector complex(wp_) :: Q(2,2) ! complex cuvature tensor real(wp_) :: M(3,3) ! local beam → tokamak change of basis matrix complex(wp_) :: S ! complex eikonal @@ -706,11 +609,42 @@ contains complex(wp_) :: hess_S(3,3) ! its Hessian matrix real(wp_) :: phi_k, phi_w ! ellipses rotation angles + ! Set beam constants + beam%nrayr = rtx%nrayr + beam%nrayth = rtx%nrayth + beam%nray = rtx%nray + beam%rwmax = rtx%rwmax + call xgygcoeff(ant%fghz, beam%k0, beam%B_res, beam%X_norm) + beam%beamtracing = rtx%igrad + + ! Allocate arrays + allocate(beam%y(6, beam%nray), beam%yp(6, beam%nray)) + allocate(beam%Xg(beam%nray), beam%Yg(beam%nray)) + allocate(beam%N_pl(beam%nray), beam%N_pr(beam%nray)) + allocate(beam%dsdvar(beam%nray), beam%dlambdadN(beam%nray)) + allocate(beam%pos(3, beam%nrayth, beam%nrayr)) + allocate(beam%grad_u(3, beam%nrayth, beam%nrayr)) + allocate(beam%grad_S(3, beam%nray)) + allocate(beam%hess_S(3, 3, beam%nray)) + allocate(beam%psi_n(beam%nray)) + allocate(beam%tau(beam%nray)) + allocate(beam%I_cd(beam%nray)) + allocate(beam%mode(beam%nray)) + + beam%step = 0 + beam%dist = [(0, j = 1,beam%nray)] + beam%tau = [(0, j = 1,beam%nray)] + beam%I_cd = [(0, j = 1,beam%nray)] + + ! Compute the wavevector of central ray given + ! the beam launcher position and orientation + call launchangles2n(ant, N_c) + ! Compute the complex curvature tensor Q from the beam parameters block real(wp_) :: K(2,2), W(2,2) - K = diag(beam%ri) ! curvature part - W = diag(2/k0 * 1/beam%w**2) ! beam width part + K = diag(ant%ri) ! curvature part + W = diag(2/beam%k0 * 1/ant%w**2) ! beam width part ! Switch from eigenbasis to local beam basis ! @@ -723,8 +657,8 @@ contains ! 2. The φs are nothing more that the angles of the rotations that ! diagonalises the respective matrix. ! - phi_k = beam%phi(2) * degree ! curvature rotation angle - phi_w = beam%phi(1) * degree ! beam width rotation angle + phi_k = ant%phi(2) * degree ! curvature rotation angle + phi_w = ant%phi(1) * degree ! beam width rotation angle K = matmul(rotate(phi_k), matmul(K, rotate(-phi_k))) W = matmul(rotate(phi_w), matmul(W, rotate(-phi_w))) @@ -773,21 +707,21 @@ contains ! For the central ray (j=1) r=0, so ∇S_I=HS_I=0 jk = 1 - grad(:,1) = 0 - hess(:,:,1) = 0 - x0_c = beam%pos - y(:,1) = [x0_c, N_c] - yp(:,1) = [N_c, 0*N_c] + beam%grad_S(:,1) = 0 + beam%hess_S(:,:,1) = 0 + x0_c = ant%pos + beam%y(:,1) = [x0_c, N_c] + beam%yp(:,1) = [N_c, 0*N_c] ! We consider the central ray as k degenerate rays ! with positions x̅=x₀_c, but ∇u computed at the points ! ξ̅_1k = Δρ [w_ξcos(α_k), w_ηcos(α_k)]. This choice ! simplifies the computation of ∇u in `gradi_upd`. do k = 1, rtx%nrayth - xi = beam%w * dr * [cos((k-1)*da), sin((k-1)*da)] + xi = ant%w * dr * [cos((k-1)*da), sin((k-1)*da)] x0 = [matmul(rotate(phi_w), xi), zero] - pos(:,k,1) = x0_c - grad_u(:,k,1) = [-matmul(Q%im, x0(1:2) * k0/(2*dr**2)), zero] + beam%pos(:,k,1) = x0_c + beam%grad_u(:,k,1) = [-matmul(Q%im, x0(1:2) * beam%k0/(2*dr**2)), zero] end do ! Compute x̅, ∇S_I, ∇u for all the other rays @@ -795,7 +729,7 @@ contains jk = fold_indices(rtx, j, k) ! Compute the ray position - xi = beam%w * (j-1)*dr * [cos((k-1)*da), sin((k-1)*da)] + xi = ant%w * (j-1)*dr * [cos((k-1)*da), sin((k-1)*da)] x0 = [matmul(rotate(phi_w), xi), zero] ! Compute the Eikonal and its derivatives @@ -803,12 +737,12 @@ contains S=S, grad_S=grad_S, hess_S=hess_S) ! Compute ∇u using the formula above - grad_u(:,k,j) = matmul(M, k0/(2 * dr**2 * (1-j)) * grad_S%im) + beam%grad_u(:,k,j) = matmul(M, beam%k0/(2 * dr**2 * (1-j)) * grad_S%im) ! Switch from local beam to tokamak basis - x0 = x0_c + matmul(M, x0) ! ray position - grad(:,jk) = matmul(M, grad_S%im) ! gradient of S_I - hess(:,:,jk) = matmul(M, matmul(hess_S%im, transpose(M))) ! Hessian of S_I + x0 = x0_c + matmul(M, x0) ! ray position + beam%grad_S(:,jk) = matmul(M, grad_S%im) ! gradient of S_I + beam%hess_S(:,:,jk) = matmul(M, matmul(hess_S%im, transpose(M))) ! Hessian of S_I ! Compute the refractive index N such that: ! @@ -840,8 +774,8 @@ contains end block ! Compute the actual initial conditions - pos(:,k,j) = x0 - y(:,jk) = [x0, matmul(M, N)] + beam%pos(:,k,j) = x0 + beam%y(:,jk) = [x0, matmul(M, N)] ! Compute the r.h.s. of the raytracing eqs. ! dx̅/dσ = + ∂Λ/∂N̅ / denom @@ -852,317 +786,24 @@ contains ! Compute travelled distance in σ and r.h.s. normalisation select case(rtx%idst) case (STEP_ARCLEN) ! σ=s, arclength - dist(jk) = 0 + beam%dist(jk) = 0 denom = norm2(N) case (STEP_TIME) ! σ=ct, time - dist(jk) = 0 + beam%dist(jk) = 0 denom = 1 case (STEP_PHASE) ! σ=S_R, phase - dist(jk) = S%re + beam%dist(jk) = S%re denom = norm2(N)**2 end select - yp(1:3,jk) = matmul(M, N) / denom - yp(4:6,jk) = matmul(M, matmul(hess_S%im, grad_S%im)) / denom + beam%yp(1:3,jk) = matmul(M, N) / denom + beam%yp(4:6,jk) = matmul(M, matmul(hess_S%im, grad_S%im)) / denom end block end do end subroutine compute_initial_conds - subroutine rkstep(params, equil, plasma, & - sox, bres, xgcn, y, yp, dgr, ddgr, igrad) - ! Runge-Kutta integrator - use gray_params, only : gray_parameters - use gray_equil, only : abstract_equil - use gray_plasma, only : abstract_plasma - - ! subroutine arguments - type(gray_parameters), intent(in) :: params - class(abstract_equil), intent(in) :: equil - class(abstract_plasma), intent(in) :: plasma - real(wp_), intent(in) :: bres, xgcn - real(wp_), intent(inout) :: y(6) - real(wp_), intent(in) :: yp(6) - real(wp_), intent(in) :: dgr(3) - real(wp_), intent(in) :: ddgr(3,3) - integer, intent(in) :: igrad,sox - - ! local variables - real(wp_), dimension(6) :: k1, k2, k3, k4 - - associate (h => params%raytracing%dst) - k1 = yp - k2 = f(y + k1*h/2) - k3 = f(y + k2*h/2) - k4 = f(y + k3*h) - y = y + h * (k1/6 + k2/3 + k3/3 + k4/6) - end associate - - contains - - function f(y) - real(wp_), intent(in) :: y(6) - real(wp_) :: f(6) - call rhs(params, equil, plasma, sox, bres, xgcn, y, dgr, ddgr, f, igrad) - end function - end subroutine rkstep - - - subroutine rhs(params, equil, plasma, & - sox, bres, xgcn, y, dgr, ddgr, dery, igrad) - ! Compute right-hand side terms of the ray equations (dery) - ! used in R-K integrator - use gray_params, only : gray_parameters - use gray_equil, only : abstract_equil - use gray_plasma, only : abstract_plasma - - ! subroutine arguments - type(gray_parameters), intent(in) :: params - class(abstract_equil), intent(in) :: equil - class(abstract_plasma), intent(in) :: plasma - real(wp_), intent(in) :: y(6) - real(wp_), intent(in) :: bres, xgcn - real(wp_), intent(in) :: dgr(3) - real(wp_), intent(in) :: ddgr(3,3) - real(wp_), intent(out) :: dery(6) - integer, intent(in) :: igrad, sox - - ! local variables - real(wp_) :: dens,btot,xg,yg - real(wp_), dimension(3) :: xv,anv,bv,derxg,deryg - real(wp_), dimension(3,3) :: derbv - - xv = y(1:3) - call plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & - bv, derbv, xg, yg, derxg, deryg) - anv = y(4:6) - call disp_deriv(params, anv, sox, xg, yg, derxg, deryg, & - bv, derbv, dgr, ddgr, igrad, dery) - end subroutine rhs - - - 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) - ! 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 - - ! subroutine rguments - type(gray_parameters), intent(in) :: params - class(abstract_equil), intent(in) :: equil - class(abstract_plasma), intent(in) :: plasma - - real(wp_), intent(in) :: xv(3), anv(3) - real(wp_), intent(in) :: dgr(3) - real(wp_), intent(in) :: ddgr(3,3) - integer, intent(in) :: sox - real(wp_), intent(in) :: bres,xgcn - real(wp_), intent(out) :: dery(6) - real(wp_), intent(out) :: psinv, dens, btot, xg, yg, anpl, anpr - real(wp_), intent(out) :: ddr, ddi, dersdst, derdnm - real(wp_), intent(out) :: bv(3) - real(wp_), intent(out) :: derxg(3) - integer, intent(in) :: igrad - - ! local variables - real(wp_):: deryg(3) - real(wp_):: derbv(3,3) - - 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) - end subroutine ywppla_upd - - - subroutine gradi_upd(params, ywrk, ak0, xc, du1, gri, ggri, error) - use gray_params, only : raytracing_parameters - use gray_errors, only : gray_error, unstable_beam, raise_error - - ! subroutine parameters - type(raytracing_parameters), intent(in) :: params - real(wp_), intent(in) :: ak0 - real(wp_), dimension(6,params%nray), intent(in) :: ywrk - real(wp_), dimension(3,params%nrayth,params%nrayr), intent(inout) :: xc, du1 - real(wp_), dimension(3,params%nray), intent(out) :: gri - real(wp_), dimension(3,3,params%nray), intent(out) :: ggri - integer(kind=gray_error), intent(inout) :: error - - ! local variables - real(wp_), dimension(3,params%nrayth,params%nrayr) :: xco, du1o - integer :: jk, j, jm, jp, k, km, kp - real(wp_) :: ux, uxx, uxy, uxz, uy, uyy, uyz, uz, uzz - real(wp_) :: dr, dfuu, dffiu, gx, gxx, gxy, gxz, gy, gyy, gyz, gz, gzz - real(wp_), dimension(3) :: dxv1, dxv2, dxv3, dgu - real(wp_), dimension(3,3) :: dgg, dff - - ! update position and du1 vectors - xco = xc - du1o = du1 - - jk = 1 - do j=1,params%nrayr - do k=1,params%nrayth - if(j>1) jk=jk+1 - xc(1:3,k,j)=ywrk(1:3,jk) - end do - end do - - ! compute grad u1 for central ray - j = 1 - jp = 2 - do k=1,params%nrayth - if(k == 1) then - km = params%nrayth - else - km = k-1 - end if - if(k == params%nrayth) then - kp = 1 - else - kp = k+1 - end if - dxv1 = xc(:,k ,jp) - xc(:,k ,j) - dxv2 = xc(:,kp,jp) - xc(:,km,jp) - dxv3 = xc(:,k ,j) - xco(:,k ,j) - call solg0(dxv1,dxv2,dxv3,dgu) - du1(:,k,j) = dgu - end do - gri(:,1) = 0 - - ! compute grad u1 and grad(S_I) for all the other rays - if (params%nrayr > 1) then - dr = params%rwmax / (params%nrayr - 1) - else - dr = 1 - end if - dfuu=2*dr**2/ak0 - jm=1 - j=2 - k=0 - dffiu = dfuu - do jk=2,params%nray - k=k+1 - if(k > params%nrayth) then - jm = j - j = j+1 - k = 1 - dffiu = dfuu*jm - end if - kp = k+1 - km = k-1 - if (k == 1) then - km=params%nrayth - else if (k == params%nrayth) then - kp=1 - end if - dxv1 = xc(:,k ,j) - xc(:,k ,jm) - dxv2 = xc(:,kp,j) - xc(:,km,j) - dxv3 = xc(:,k ,j) - xco(:,k ,j) - - call solg0(dxv1, dxv2, dxv3, dgu) - - ! Check for numerical instability in the computation of ∇S_I. - ! This can happen when the distances between the rays and the - ! integration step size differ wildly; in such cases the linear - ! system Δx⋅∇u = Δu becomes bad conditioned. - ! - ! Source: https://discourse.julialang.org/t/what-is-the-best-way-to-solve-a-ill-conditioned-linear-problem/95894/10 - block - use const_and_precisions, only : comp_eps - use utils, only : singvals - real(wp_) :: A(3,3), x(3), b(3), sigma(3) - real(wp_) :: R, normA, back_error, forw_error - A = reshape([dxv1, dxv2, dxv3], shape=[3,3], order=[2,1]) - sigma = singvals(A) - b = [1, 0, 0] - x = dgu - R = maxval(sigma) / minval(sigma) ! condition number - normA = maxval(sigma) ! L₂ norm of A - back_error = norm2(matmul(A, x) - b) / (normA * norm2(b)) - forw_error = back_error * R * norm2(x) - ! TODO: figure out how to tie this threshold - ! to the derivatives of the dispersion relation - if (forw_error / comp_eps > 5000) then - error = raise_error(error, unstable_beam) - end if - end block - - du1(:,k,j) = dgu - gri(:,jk) = dgu(:)*dffiu - end do - - ! compute derivatives of grad u and grad(S_I) for rays jk>1 - ggri(:,:,1) = 0 - jm=1 - j=2 - k=0 - dffiu = dfuu - do jk=2,params%nray - k=k+1 - if(k > params%nrayth) then - jm=j - j=j+1 - k=1 - dffiu = dfuu*jm - end if - kp=k+1 - km=k-1 - if (k == 1) then - km=params%nrayth - else if (k == params%nrayth) then - kp=1 - end if - dxv1 = xc(:,k ,j) - xc(:,k ,jm) - dxv2 = xc(:,kp,j) - xc(:,km,j) - dxv3 = xc(:,k ,j) - xco(:,k ,j) - dff(:,1) = du1(:,k ,j) - du1(:,k ,jm) - dff(:,2) = du1(:,kp,j) - du1(:,km,j) - dff(:,3) = du1(:,k ,j) - du1o(:,k ,j) - call solg3(dxv1,dxv2,dxv3,dff,dgg) - - ! derivatives of u - ux = du1(1,k,j) - uy = du1(2,k,j) - uz = du1(3,k,j) - uxx = dgg(1,1) - uyy = dgg(2,2) - uzz = dgg(3,3) - uxy = (dgg(1,2) + dgg(2,1))/2 - uxz = (dgg(1,3) + dgg(3,1))/2 - uyz = (dgg(2,3) + dgg(3,2))/2 - - ! derivatives of S_I and Grad(S_I) - gx = ux*dffiu - gy = uy*dffiu - gz = uz*dffiu - gxx = dfuu*ux*ux + dffiu*uxx - gyy = dfuu*uy*uy + dffiu*uyy - gzz = dfuu*uz*uz + dffiu*uzz - gxy = dfuu*ux*uy + dffiu*uxy - gxz = dfuu*ux*uz + dffiu*uxz - gyz = dfuu*uy*uz + dffiu*uyz - - ggri(1,1,jk)=gxx - ggri(2,1,jk)=gxy - ggri(3,1,jk)=gxz - ggri(1,2,jk)=gxy - ggri(2,2,jk)=gyy - ggri(3,2,jk)=gyz - ggri(1,3,jk)=gxz - ggri(2,3,jk)=gyz - ggri(3,3,jk)=gzz - end do - - end subroutine gradi_upd - - - subroutine solg0(dxv1,dxv2,dxv3,dgg) ! solution of the linear system of 3 eqs : dgg . dxv = dff ! input vectors : dxv1, dxv2, dxv3, dff @@ -1214,11 +855,10 @@ contains end subroutine solg3 - subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & + pure subroutine plas_deriv(equil, plasma, xv, bres, xgcn, dens, btot, & bv, derbv, xg, yg, derxg, deryg, psinv_opt) use const_and_precisions, only : cm - use gray_equil, only : abstract_equil, vacuum - use gray_plasma, only : abstract_plasma + use gray_equil, only : vacuum ! subroutine arguments class(abstract_equil), intent(in) :: equil @@ -1366,9 +1006,8 @@ contains end subroutine plas_deriv - - subroutine disp_deriv(params, anv, sox, xg, yg, derxg, deryg, bv, derbv, & - dgr, ddgr, igrad, dery, anpl_, anpr, ddr, ddi, & + pure subroutine disp_deriv(var, anv, sox, xg, yg, derxg, deryg, bv, derbv, & + dgr, ddgr, beamtracing, dery, anpl_, anpr, ddr, ddi, & dersdst, derdnm_) ! Computes the dispersion relation, derivatives and other ! related quantities @@ -1379,14 +1018,14 @@ contains ! computing the absoprtion and current drive. use const_and_precisions, only : zero, half - use gray_params, only : gray_parameters, STEP_ARCLEN, & + use gray_params, only : step_enum, STEP_ARCLEN, & STEP_TIME, STEP_PHASE ! subroutine arguments ! Inputs - type(gray_parameters), intent(in) :: params + integer(kind(step_enum)), intent(in) :: var ! refractive index N̅ vector, b̅ = B̅/B magnetic field direction real(wp_), dimension(3), intent(in) :: anv, bv ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X @@ -1401,7 +1040,7 @@ contains ! gradient of the magnetic field direction real(wp_), dimension(3, 3), intent(in) :: derbv ! ∇b̅ ! raytracing/beamtracing switch - integer, intent(in) :: igrad + logical, intent(in) :: beamtracing ! Outputs @@ -1474,7 +1113,7 @@ contains dan2sdnpl = - xg*yg2*anpl/duh & - sox*xg*anpl*(2*(1 - xg) - yg2*dnl)/(del*duh) - if(igrad > 0) then + if (beamtracing) then ! Derivatives used in the complex eikonal terms (beamtracing only) block real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel @@ -1523,7 +1162,7 @@ contains ! on the constrains Λ=0. derdom = -2*an2 + 2*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl - if (igrad > 0) then + if (beamtracing) then ! Complex eikonal terms added to the above expressions block real(wp_) :: dgr2(3) @@ -1562,7 +1201,7 @@ contains ! dx̅/dσ = + ∂Λ/∂N̅ / denom ! dN̅/dσ = - ∂Λ/∂x̅ / denom ! - select case (params%raytracing%idst) + select case (var) case (STEP_ARCLEN) ! σ=s, arclength ! denom = |∂Λ/∂N̅| ! Proof: Normalising ∂Λ/∂N̅ (∥ to the group velocity) @@ -1613,7 +1252,7 @@ contains ! if (present(dersdst)) dersdst = derdnm/denom - ! |∂Λ/∂N̅|: used in the α computation (see alpha_effj) + ! |∂Λ/∂N̅|: used in the α computation if (present(derdnm_)) derdnm_ = derdnm ! r.h.s. vector @@ -1627,179 +1266,20 @@ contains ddr = an2 - an2s end if - if (present(ddr) .and. igrad > 0) then + if (present(ddr) .and. beamtracing) then ! Dispersion relation (complex eikonal) ! ddr ← Λ = N² - N²s - |∇S_I|² + ½ (b̅⋅∇S_I)² ∂²N²s/∂N∥² ddr = ddr - gr2 - half*bdotgr**2*fdia ! real part end if - ! Note: we have to return ddi even for igrad=0 because + ! Note: we have to return ddi even for beamtracing=0 because ! it's printed to udisp unconditionally if (present(ddi)) then ! Dispersion relation (complex eikonal) ! ddi ← ∂Λ/∂N̅⋅∇S_I - ddi = merge(dot_product(derdnv, dgr), zero, igrad > 0) ! imaginary part + ddi = merge(dot_product(derdnv, dgr), zero, beamtracing) ! imaginary part end if end subroutine disp_deriv - - - subroutine alpha_effj(params, equil, plasma, psinv, X, Y, density, & - temperature, k0, Bres, derdnm, Npl, Npr_cold, & - sox, Npr, alpha, dIdp, nhmin, nhmax, iokhawa, error) - ! Computes the absorption coefficient and effective current - - use const_and_precisions, only : pi, mc2=>mc2_ - use gray_params, only : ecrh_cd_parameters - use gray_equil, only : abstract_equil - use gray_plasma, only : abstract_plasma - use dispersion, only : harmnumber, warmdisp - use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl - use gray_errors, only : gray_error, negative_absorption, raise_error - - - ! subroutine arguments - - ! Inputs - - ! ECRH & CD parameters - type(ecrh_cd_parameters), intent(in) :: params - ! MHD equilibrium, plasma object - class(abstract_equil), intent(in) :: equil - class(abstract_plasma), intent(in) :: plasma - ! poloidal flux ψ - real(wp_), intent(in) :: psinv - ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω - real(wp_), intent(in) :: X, Y - ! densityity [10¹⁹ m⁻³], temperature [keV] - real(wp_), intent(in) :: density, temperature - ! vacuum wavenumber k₀=ω/c, resonant B field - real(wp_), intent(in) :: k0, Bres - ! group velocity: |∂Λ/∂N̅| where Λ=c²/ω²⋅D_R - real(wp_), intent(in) :: derdnm - ! refractive index: N∥, N⊥ (cold dispersion) - real(wp_), intent(in) :: Npl, Npr_cold - ! sign of polarisation mode: -1 ⇒ O, +1 ⇒ X - integer, intent(in) :: sox - - ! Outputs - - ! orthogonal refractive index N⊥ (solution of the warm dispersion) - complex(wp_), intent(out) :: Npr - ! absorption coefficient, current density - real(wp_), intent(out) :: alpha, dIdp - ! lowest/highest resonant harmonic numbers - integer, intent(out) :: nhmin, nhmax - ! ECCD/overall error codes - integer(kind=gray_error), intent(inout) :: iokhawa, error - - ! local variables - real(wp_) :: rbavi, rrii, rhop - integer :: nlarmor, ithn, ierrcd - real(wp_) :: mu, rbn, rbx - real(wp_) :: zeff, cst2, bmxi, bmni, fci - real(wp_), dimension(:), allocatable :: eccdpar - real(wp_) :: effjcd, effjcdav, Btot - complex(wp_) :: e(3) - - alpha = 0 - Npr = 0 - dIdp = 0 - nhmin = 0 - nhmax = 0 - iokhawa = 0 - - ! Absorption computation - - ! Skip when the temperature is zero (no plasma) - if (temperature <= 0) return - - ! Skip when the lowest resonant harmonic number is zero - mu = mc2/temperature ! μ=(mc²/kT) - call harmnumber(Y, mu, Npl**2, params%iwarm == 1, nhmin, nhmax) - if (nhmin <= 0) return - - ! Solve the full dispersion only when needed - if (params%iwarm /= 4 .or. params%ieccd /= 0) then - nlarmor = max(params%ilarm, nhmax) - if (params%ieccd /= 0) then - ! Compute the polarisation vector only when current drive is on - call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, e, & - model=params%iwarm, nlarmor=nlarmor, & - max_iters=abs(params%imx), fallback=params%imx < 0) - else - call warmdisp(X, Y, mu, Npl, Npr_cold, sox, error, Npr, & - model=params%iwarm, nlarmor=nlarmor, & - max_iters=abs(params%imx), fallback=params%imx < 0) - end if - end if - - ! Compute α from the solution of the dispersion relation - ! The absoption coefficient is defined as - ! - ! α = 2 Im(k̅)⋅s̅ - ! - ! where s̅ = v̅_g/|v_g|, the direction of the energy flow. - ! Since v̅_g = - ∂Λ/∂N̅ / ∂Λ/∂ω, using the cold dispersion - ! relation, we have that - ! - ! s̅ = ∂Λ/∂N̅ / |∂Λ/∂N̅| - ! = [2N̅ - ∂(N²s)/∂N∥ b̅ - ! + ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅] / |∂Λ/∂N̅| - ! - ! Assuming Im(k∥)=0: - ! - ! α = 4 Im(k⊥)⋅N⊥ / |∂Λ/∂N̅| - ! - block - real(wp_) :: k_im - k_im = k0 * Npr%im ! imaginary part of k⊥ - alpha = 4 * k_im*Npr_cold / derdnm - end block - - if (alpha < 0) then - error = raise_error(error, negative_absorption) - return - end if - - ! Current drive computation - if (params%ieccd <= 0) return - - zeff = plasma%zeff(psinv) - ithn = 1 - if (nlarmor > nhmin) ithn = 2 - rhop = sqrt(psinv) - call equil%flux_average(rhop, R_J=rrii, B_avg=rbavi, & - B_min=bmni, B_max=bmxi, f_c=fci) - rbavi = rbavi / bmni - Btot = Y*Bres - rbn = Btot/bmni - rbx = Btot/bmxi - - select case(params%ieccd) - case(1) - ! Cohen model - call setcdcoeff(zeff, rbn, rbx, cst2, eccdpar) - call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & - ithn, cst2, fjch, eccdpar, effjcd, iokhawa, ierrcd) - case(2) - ! No trapping - call setcdcoeff(zeff, cst2, eccdpar) - call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & - ithn, cst2, fjch0, eccdpar, effjcd, iokhawa, ierrcd) - case default - ! Neoclassical model - call setcdcoeff(zeff, rbx, fci, mu, rhop, equil%spline_cd_eff, cst2, eccdpar) - call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, & - ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd) - end select - error = error + ierrcd - - ! current drive efficiency R* = / [A⋅m/W] - effjcdav = rbavi*effjcd - dIdp = equil%sgn_bphi * effjcdav / (2*pi*rrii) - - end subroutine alpha_effj - end module gray_core diff --git a/src/gray_params.f90 b/src/gray_params.f90 index d76dd94..00fc951 100644 --- a/src/gray_params.f90 +++ b/src/gray_params.f90 @@ -135,14 +135,14 @@ module gray_params ! Raytracing parameters type raytracing_parameters - real(wp_) :: dst ! Integration step size - real(wp_) :: rwmax = 1 ! Normalized maximum radius of beam power - integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays - integer :: igrad = 0 ! Complex eikonal switch - integer :: nstep = 12000 ! Max number of integration steps - integer(kind(step_enum)) :: idst = STEP_TIME ! Choice of the integration variable - integer :: ipass = 1 ! Number of plasma passes - logical :: ipol = .false. ! Whether to compute wave polarisation (from chi, psi) + real(wp_) :: dst ! Integration step size + real(wp_) :: rwmax = 1 ! Normalized maximum radius of beam power + integer :: nray, nrayr, nrayth ! Total/radial/angular number of rays + logical :: igrad = .false. ! Complex eikonal switch + integer :: nstep = 12000 ! Max number of integration steps + integer(kind(step_enum)) :: idst = STEP_TIME ! Choice of the integration variable + integer :: ipass = 1 ! Number of plasma passes + logical :: ipol = .false. ! Whether to compute wave polarisation (from chi, psi) end type ! EC resonant heating & Current Drive parameters @@ -341,7 +341,7 @@ contains call append(header, line) ! code parameters - write(line, '("# COD igrad idst ipass ipol:",3(1x,i4),1x,l4)') & + write(line, '("# COD igrad idst ipass ipol:",l4,2(1x,i4),1x,l4)') & params%raytracing%igrad, params%raytracing%idst, & params%raytracing%ipass, params%raytracing%ipol call append(header, line) diff --git a/src/gray_plasma.f90 b/src/gray_plasma.f90 index 30a198e..066d9b3 100644 --- a/src/gray_plasma.f90 +++ b/src/gray_plasma.f90 @@ -21,7 +21,7 @@ module gray_plasma end type abstract interface - subroutine density_sub(self, psin, dens, ddens) + pure subroutine density_sub(self, psin, dens, ddens) ! Computes the density its first derivative as a function of ! normalised poloidal flux. ! @@ -32,7 +32,7 @@ module gray_plasma real(wp_), intent(out) :: dens, ddens ! density and first derivative end subroutine density_sub - function temp_fun(self, psin) result(temp) + pure function temp_fun(self, psin) result(temp) ! Computes the temperature as a function of the ! normalised poloidal flux. ! @@ -43,7 +43,7 @@ module gray_plasma real(wp_) :: temp end function temp_fun - function zeff_fun(self, psin) result(zeff) + pure function zeff_fun(self, psin) result(zeff) ! Computes the effective charge Z_eff as a ! function of the normalised poloidal flux. import :: abstract_plasma, wp_ @@ -97,7 +97,7 @@ module gray_plasma contains - subroutine analytic_density(self, psin, dens, ddens) + pure subroutine analytic_density(self, psin, dens, ddens) ! subroutine arguments class(analytic_plasma), intent(in) :: self real(wp_), intent(in) :: psin ! normalised poloidal flux @@ -116,7 +116,7 @@ contains end subroutine analytic_density - function analytic_temp(self, psin) result(temp) + pure function analytic_temp(self, psin) result(temp) ! function arguments class(analytic_plasma), intent(in) :: self real(wp_), intent(in) :: psin @@ -139,7 +139,7 @@ contains end function analytic_temp - function analytic_zeff(self, psin) result(zeff) + pure function analytic_zeff(self, psin) result(zeff) ! function arguments class(analytic_plasma), intent(in) :: self real(wp_), intent(in) :: psin @@ -155,7 +155,7 @@ contains end function analytic_zeff - subroutine numeric_density(self, psin, dens, ddens) + pure subroutine numeric_density(self, psin, dens, ddens) use logger, only : log_error ! subroutine arguments @@ -221,16 +221,16 @@ contains end block end if - if (dens < 0) then - write (msg, '("negative density:", 2(x,a,"=",g0.3))') & - 'ne', dens, 'ψ', psin - call log_error(msg, mod='coreprofiles', proc='density') - end if + ! if (dens < 0) then + ! write (msg, '("negative density:", 2(x,a,"=",g0.3))') & + ! 'ne', dens, 'ψ', psin + ! call log_error(msg, mod='coreprofiles', proc='density') + ! end if end subroutine numeric_density - function numeric_temp(self, psin) result(temp) + pure function numeric_temp(self, psin) result(temp) ! function arguments class(numeric_plasma), intent(in) :: self real(wp_), intent(in) :: psin @@ -246,7 +246,7 @@ contains end function numeric_temp - function numeric_zeff(self, psin) result(zeff) + pure function numeric_zeff(self, psin) result(zeff) ! function arguments class(numeric_plasma), intent(in) :: self real(wp_), intent(in) :: psin diff --git a/src/main.f90 b/src/main.f90 index de992b7..5d0f673 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -77,10 +77,10 @@ program main ! Do some checks on the inputs associate (p => params%raytracing) - if (p%igrad == 1 .and. p%nrayr < 5) then - p%igrad = 0 + if (p%igrad .and. p%nrayr < 5) then + p%igrad = .false. call log_message(level=WARNING, mod='main', & - msg='igrad = 1 but nrayr < 5: disabling beamtracing') + msg='igrad = .true. but nrayr < 5: disabling beamtracing') end if if (p%nrayr == 1) p%nrayth = 1 diff --git a/src/main_convert.f90 b/src/main_convert.f90 index 46f6099..cf0aa93 100644 --- a/src/main_convert.f90 +++ b/src/main_convert.f90 @@ -156,8 +156,8 @@ contains write(u, fmt) "nrayth", params%raytracing%nrayth if (params%raytracing%rwmax /= defaults%raytracing%rwmax) & write(u, fmt) "rwmax", params%raytracing%rwmax - if (params%raytracing%igrad /= defaults%raytracing%igrad) & - write(u, fmt) "igrad", params%raytracing%igrad + if (params%raytracing%igrad .neqv. defaults%raytracing%igrad) & + write(u, fmt) "igrad", format_bool(params%raytracing%igrad) if (params%raytracing%ipass /= defaults%raytracing%ipass) & write(u, fmt) "ipass", params%raytracing%ipass if (params%raytracing%ipol .neqv. defaults%raytracing%ipol) & diff --git a/tests/06-ITER-startup/inputs/gray.ini b/tests/06-ITER-startup/inputs/gray.ini index fa988a1..98e0bd6 100644 --- a/tests/06-ITER-startup/inputs/gray.ini +++ b/tests/06-ITER-startup/inputs/gray.ini @@ -1,7 +1,7 @@ [raytracing] nrayr = 11 nrayth = 16 -igrad = 1 +igrad = true dst = 0.1 nstep = 8000