From c7d0d8370ce530ec1b7de5a7711ae04cb03af8be Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Wed, 29 Mar 2023 22:13:43 +0200 Subject: [PATCH] src/gray_core.f90: make some {disp,plas}_deriv outputs optional Some of the outputs of disp_deriv and plas_deriv are only needed when updating the local plasma quantities (ywppla_upd) and not when integrating the raytracing equations (rkstep). This change save some unnecessary computations and variable definitions. Also add some comments to disp_deriv --- src/gray_core.f90 | 344 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 243 insertions(+), 101 deletions(-) diff --git a/src/gray_core.f90 b/src/gray_core.f90 index df6bb4e..320bd14 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -1004,18 +1004,14 @@ contains real(wp_), dimension(6), intent(out) :: dery integer, intent(in) :: igrad ! local variables - real(wp_) :: psinv,dens,btot,xg,yg,anpl,anpr,ajphi - real(wp_) :: ddr,ddi,dersdst,derdnm + real(wp_) :: psinv,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(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg, & - ajphi) - + call plas_deriv(xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg) anv = y(4:6) - call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & - dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) + call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad,dery) end subroutine rhs @@ -1039,16 +1035,13 @@ contains real(wp_), dimension(3), intent(out) :: derxg integer, intent(in) :: igrad ! local variables - real(wp_) :: gr2,ajphi - real(wp_), dimension(3) :: dgr2,deryg + real(wp_), dimension(3) :: deryg real(wp_), dimension(3,3) :: derbv real(wp_), parameter :: anplth1 = 0.99_wp_, anplth2 = 1.05_wp_ - gr2 = dot_product(dgr, dgr) - dgr2 = 2 * matmul(ddgr, dgr) - call plas_deriv(xv,bres,xgcn,psinv,dens,btot,bv,derbv,xg,yg,derxg,deryg,ajphi) - call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,gr2,dgr2,dgr,ddgr, & - dery,anpl,anpr,ddr,ddi,dersdst,derdnm,igrad) + call plas_deriv(xv,bres,xgcn,dens,btot,bv,derbv,xg,yg,derxg,deryg,psinv) + call disp_deriv(anv,sox,xg,yg,derxg,deryg,bv,derbv,dgr,ddgr,igrad, & + dery,anpl,anpr,ddr,ddi,dersdst,derdnm) error=0 if( abs(anpl) > anplth1) then @@ -1262,9 +1255,9 @@ contains - subroutine plas_deriv(xv, bres, xgcn, psinv, dens, btot, bv, derbv, & - xg, yg, derxg, deryg, ajphi) - use const_and_precisions, only : zero, ccj=>mu0inv + subroutine plas_deriv(xv, bres, xgcn, dens, btot, bv, derbv, & + xg, yg, derxg, deryg, psinv_opt) + use const_and_precisions, only : zero use gray_params, only : iequil use equilibrium, only : psia, equinum_fpol, equinum_psi, equian, sgnbphi use coreprofiles, only : density @@ -1272,19 +1265,21 @@ contains implicit none ! subroutine arguments - real(wp_), dimension(3), intent(in) :: xv - real(wp_), intent(in) :: xgcn,bres - real(wp_), intent(out) :: psinv,dens,btot,xg,yg - real(wp_), dimension(3), intent(out) :: bv,derxg,deryg + real(wp_), dimension(3), intent(in) :: xv + real(wp_), intent(in) :: xgcn, bres + real(wp_), intent(out) :: dens, btot, xg, yg + real(wp_), dimension(3), intent(out) :: bv, derxg, deryg real(wp_), dimension(3,3), intent(out) :: derbv + real(wp_), optional, intent(out) :: psinv_opt ! local variables integer :: jv - real(wp_) :: xx,yy,zz + real(wp_) :: xx,yy,zz + real(wp_) :: psinv real(wp_) :: csphi,drrdx,drrdy,dphidx,dphidy,rr,rr2,rrm,snphi,zzm real(wp_), dimension(3) :: dbtot,bvc real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv - real(wp_) :: brr,bphi,bzz,ajphi,dxgdpsi + real(wp_) :: brr,bphi,bzz,dxgdpsi real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddenspsin xg = zero @@ -1292,13 +1287,16 @@ contains psinv = -1._wp_ dens = zero btot = zero - ajphi = zero derxg = zero deryg = zero bv = zero derbv = zero - if(iequil==0) return + if(iequil==0) then + ! copy optional output + if (present(psinv_opt)) psinv_opt = psinv + return + end if dbtot = zero dbv = zero @@ -1331,6 +1329,9 @@ contains call equinum_fpol(psinv,fpolv,dfpv) end if + ! copy optional output + if (present(psinv_opt)) psinv_opt = psinv + ! compute yg and derivative if(psinv < zero) then bphi = fpolv/rrm @@ -1406,47 +1407,83 @@ contains derxg(1) = 1.0e-2_wp_*drrdx*dpsidr*dxgdpsi derxg(2) = 1.0e-2_wp_*drrdy*dpsidr*dxgdpsi derxg(3) = 1.0e-2_wp_*dpsidz *dxgdpsi - - ! current density computation in Ampere/m^2, ccj==1/mu_0 - ajphi = ccj*(dbvcdc(1,3) - dbvcdc(3,1)) - !ajr=ccj*(dbvcdc(3,2)/rrm-dbvcdc(2,3)) - !ajz=ccj*(bvc(2)/rrm+dbvcdc(2,1)-dbvcdc(1,2)) - end subroutine plas_deriv - subroutine disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, & - gr2, dgr2, dgr, ddgr, dery, anpl, anpr, & - ddr, ddi, dersdst, derdnm, igrad) + subroutine disp_deriv(anv, sox, xg, yg, derxg, deryg, bv, derbv, & + dgr, ddgr, igrad, dery, anpl_, anpr, ddr, ddi, & + dersdst, derdnm_) + ! Computes the dispersion relation, derivatives and other + ! related quantities + ! + ! The mandatory outputs are used for both integrating the ray + ! trajectory (`rhs` subroutine) and updating the ray state + ! (`ywppla_upd` suborutine); while the optional ones are used for + ! computing the absoprtion and current drive. + use const_and_precisions, only : zero, one, half, two use gray_params, only : idst implicit none ! subroutine arguments - real(wp_), intent(in) :: xg,yg,gr2,sox - real(wp_), intent(out) :: anpl,anpr,ddr,ddi,derdnm,dersdst - real(wp_), dimension(3), intent(in) :: anv,bv,derxg,deryg - real(wp_), dimension(3), intent(in) :: dgr2,dgr - real(wp_), dimension(3,3), intent(in) :: ddgr,derbv - real(wp_), dimension(6), intent(out) :: dery - integer, intent(in) :: igrad + + ! Inputs + + ! 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 + real(wp_), intent(in) :: sox + ! CMA diagram variables: X=(ω_pe/ω)², Y=ω_ce/ω + real(wp_), intent(in) :: xg, yg + ! gradients of X, Y + real(wp_), dimension(3), intent(in) :: derxg, deryg + ! gradients of the complex eikonal + real(wp_), dimension(3), intent(in) :: dgr ! ∇S_I + real(wp_), dimension(3, 3), intent(in) :: ddgr ! ∇∇S_I + ! gradient of the magnetic field direction + real(wp_), dimension(3, 3), intent(in) :: derbv ! ∇b̅ + ! raytracing/beamtracing switch + integer, intent(in) :: igrad + + ! Outputs + + ! the actual derivatives: (∂Λ/∂x̅, -∂Λ/∂N̅) + real(wp_), dimension(6), intent(out) :: dery + ! additional quantities: + ! refractive index + real(wp_), optional, intent(out) :: anpl_, anpr ! N∥, N⊥ + ! real and imaginary part of the dispersion + real(wp_), optional, intent(out) :: ddr, ddi ! Λ, ∂Λ/∂N̅⋅∇S_I + ! Jacobian ds/dσ, where s arclength, σ integration variable + real(wp_), optional, intent(out) :: dersdst + ! |∂Λ/∂N̅| + real(wp_), optional, intent(out) :: derdnm_ + + ! Note: assign values to missing optional arguments results in a segfault. + ! Since some are always needed anyway, we store them here and copy + ! them later, if needed: + real(wp_) :: anpl, derdnm ! local variables - real(wp_) :: yg2, anpl2, anpr2, del, dnl, duh, dan2sdnpl, an2, an2s - real(wp_) :: dan2sdxg, dan2sdyg, ddelnpl2, ddelnpl2x, ddelnpl2y, denom, derdel + real(wp_) :: gr2, yg2, anpl2, del, dnl, duh, dan2sdnpl, an2, an2s + real(wp_) :: dan2sdxg, dan2sdyg, denom real(wp_) :: derdom, dfdiadnpl, dfdiadxg, dfdiadyg, fdia, bdotgr real(wp_), dimension(3) :: derdxv, danpldxv, derdnv, dbgr - an2 = dot_product(anv, anv) - anpl = dot_product(anv, bv) + an2 = dot_product(anv, anv) ! N² + anpl = dot_product(anv, bv) ! N∥ = N̅⋅B̅ - anpl2 = anpl**2 - dnl = one - anpl2 - anpr2 = max(an2-anpl2,zero) - anpr = sqrt(anpr2) - yg2 = yg**2 + ! Shorthands used in the expressions below + yg2 = yg**2 ! Y² + anpl2 = anpl**2 ! N∥² + dnl = one - anpl2 ! 1 - N∥² + duh = one - xg - yg2 ! UH denom (duh=0 on the upper-hybrid resonance) + + ! Compute/copy optional outputs + if (present(anpr)) anpr = sqrt(max(an2 - anpl2, zero)) ! N⊥ + if (present(anpl_)) anpl_ = anpl ! N∥ an2s = one dan2sdxg = zero @@ -1458,17 +1495,16 @@ contains dfdiadxg = zero dfdiadyg = zero - duh = one - xg - yg2 if(xg > zero) then - ! Derivatives of the cold plasma dispersion relation + ! Derivatives of the cold plasma refractive index ! - ! Λ = N²⊥ - N²s(X,Y,N∥) = 0 + ! N²s = 1 - X - XY²⋅(1 + N∥² ± √Δ)/[2(1 - X - Y²)] ! - ! where N²s = 1 - X - XY²⋅(1 + N∥ ∓ √Δ)/[2(1 - X - Y²)] - ! Δ = (1 - N∥)² + 4N∥⋅(1 - X)/Y² - ! - del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) - an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh + ! where Δ = (1 - N∥²)² + 4N∥²⋅(1 - X)/Y² + ! + for the X mode, - for the O mode + + ! √Δ + del = sqrt(dnl**2 + 4.0_wp_*anpl2*(one - xg)/yg2) ! ∂(N²s)/∂X ! Note: this term is nonzero for X=0, but it multiplies terms @@ -1483,66 +1519,172 @@ contains - sox*xg*anpl*(two*(one - xg) - yg2*dnl)/(del*duh) if(igrad > 0) then - ! Derivatives of the eikonal (beamtracing only) - ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) & - - yg2*dnl**3)/yg2/del**3 - fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh - derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) & - - dnl**2*(one + 3.0_wp_*anpl2)*yg2 - derdel = 4.0_wp_*derdel/(yg*del)**5 - ddelnpl2y = two*(one - xg)*derdel - ddelnpl2x = yg*derdel - dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) & - /(yg2*del**5) - dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) & - *ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2) - dfdiadyg = - two*yg*xg*(one - xg)/duh**2 & - - sox*xg*yg*(two*(one - xg)*ddelnpl2 & - + yg*duh*ddelnpl2y)/(two*duh**2) + ! Derivatives used in the complex eikonal terms (beamtracing only) + block + real(wp_) :: ddelnpl2, ddelnpl2x, ddelnpl2y, derdel + + ! ∂²Δ/∂N∥² + ddelnpl2 = two*(two*(one - xg)*(one + 3.0_wp_*anpl2**2) & + - yg2*dnl**3)/yg2/del**3 + ! ∂²(N²s)/∂N∥² + fdia = - xg*yg2*(one + half*sox*ddelnpl2)/duh + + ! Intermediates results used right below + derdel = two*(one - xg)*anpl2*(one + 3.0_wp_*anpl2**2) & + - dnl**2*(one + 3.0_wp_*anpl2)*yg2 + derdel = 4.0_wp_*derdel/(yg*del)**5 + ddelnpl2y = two*(one - xg)*derdel + ddelnpl2x = yg*derdel + + ! ∂³(N²s)/∂N∥³ + dfdiadnpl = 24.0_wp_*sox*xg*(one - xg)*anpl*(one - anpl2**2) & + /(yg2*del**5) + ! ∂³(N²s)/∂N∥²∂X + dfdiadxg = - yg2*(one - yg2)/duh**2 - sox*yg2*((one - yg2) & + *ddelnpl2 + xg*duh*ddelnpl2x)/(two*duh**2) + ! ∂³(N²s)/∂N∥²∂Y + dfdiadyg = - two*yg*xg*(one - xg)/duh**2 & + - sox*xg*yg*(two*(one - xg)*ddelnpl2 & + + yg*duh*ddelnpl2y)/(two*duh**2) + end block end if end if - bdotgr = dot_product(bv, dgr) - dbgr = matmul(dgr, derbv) + matmul(bv, ddgr) + ! ∇(N∥) = ∇(N̅⋅b̅) = N̅⋅∇b̅ danpldxv = matmul(anv, derbv) - derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl + & - igrad*dgr2) & - + fdia*bdotgr*dbgr + half*bdotgr**2 & - *(derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl) - derdnv = two*anv + (half*bdotgr**2*dfdiadnpl - dan2sdnpl)*bv + ! ∂Λ/∂x̅ = - ∇(N²s) = -∂Λ/∂X ∇X - ∂Λ/∂Y ∇Y - ∂Λ/∂N∥ ∇(N∥) + derdxv = -(derxg*dan2sdxg + deryg*dan2sdyg + danpldxv*dan2sdnpl) - derdnm = norm2(derdnv) + ! ∂Λ/∂N̅ = 2N̅ - ∂(N²s)/∂N̅ = 2N̅ - ∂(N²s)/∂N∥ b̅ + ! Note: we used the identity ∇f(v̅⋅b̅) = f' ∇(v̅⋅b̅) = f'b̅. + derdnv = two*anv - dan2sdnpl*bv - derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl & - + two*igrad*gr2 - bdotgr**2*(fdia + xg*dfdiadxg & - + half*yg*dfdiadyg & - + half*anpl*dfdiadnpl) + ! ∂Λ/∂ω = ∂N²/∂ω - ∂N²s/∂X⋅∂X/∂ω - ∂N²s/∂Y⋅∂Y/∂ω - ∂N²s/∂N∥⋅∂N∥/∂ω + ! Notes: 1. N depends on ω: N²=c²k²/ω² ⇒ ∂N²/∂ω = -2N²/ω + ! N∥=ck∥/ω ⇒ ∂N∥/∂ω = -N∥/ω + ! 2. derdom is actually ω⋅∂Λ/∂ω, see below for the reason. + ! 3. N gains a dependency on ω because Λ(∇S, ω) is computed + ! on the constrains Λ=0. + derdom = -two*an2 + two*xg*dan2sdxg + yg*dan2sdyg + anpl*dan2sdnpl - ! integration variable σ - select case(idst) - case(0) - ! optical path: s + if (igrad > 0) then + ! Complex eikonal terms added to the above expressions + block + real(wp_) :: dgr2(3) + + bdotgr = dot_product(bv, dgr) ! b̅⋅∇S_I + gr2 = dot_product(dgr, dgr) ! |∇S_I|² + dgr2 = 2 * matmul(ddgr, dgr) ! ∇|∇S_I|² = 2 ∇∇S_I⋅∇S_I + + ! ∇(b̅⋅∇S_I) = ∇b̅ᵗ ∇S_I + b̅ᵗ ∇∇S_I + ! Notes: 1. We are using the convention ∇b̅_ij = ∂b_i/∂x_j. + ! 2. Fortran doesn't distinguish between column/row vectors, + ! so matmul(A, v̅) ≅ Av̅ and matmul(v̅, A) ≅ v̅ᵗA ≅ Aᵗv̅. + dbgr = matmul(dgr, derbv) + matmul(bv, ddgr) + + ! ∂Λ/∂x̅ += - ∇(|∇S_I|²) + ½ ∇[(b̅⋅∇S_I)² ∂²N²s/∂N∥²] + derdxv = derdxv - dgr2 + fdia*bdotgr*dbgr + half*bdotgr**2 & + * (derxg*dfdiadxg + deryg*dfdiadyg + danpldxv*dfdiadnpl) + + ! ∂Λ/∂N̅ += ½(b̅⋅∇S_I)² ∂³(N²s)/∂N∥³ b̅ + ! Note: we used again the identity ∇f(v̅⋅b̅) = f'b̅. + derdnv = derdnv + half*bdotgr**2*dfdiadnpl*bv + + ! ∂Λ/∂ω += ∂|∇S_I|²/∂ω + ½∂(b⋅∇S_I)²/∂ω + ½(b̅⋅∇S_I)² ∂/∂ω (∂²N²s/∂N∥²) + ! Note: as above ∇S_I gains a dependency on ω + derdom = derdom + two*gr2 - bdotgr**2 & + * (fdia + xg*dfdiadxg + half*yg*dfdiadyg & + + half*anpl*dfdiadnpl) + end block + end if + + derdnm = norm2(derdnv) + + ! Denominator of the r.h.s, depending on the + ! choice of the integration variable σ: + ! + ! dx̅/dσ = + ∂Λ/∂N̅ / denom + ! dN̅/dσ = - ∂Λ/∂x̅ / denom + ! + select case (idst) + case (0) ! σ=s, arclength + ! denom = |∂Λ/∂N̅| + ! Proof: Normalising ∂Λ/∂N̅ (∥ to the group velocity) + ! simply makes σ the arclength parameter s. denom = derdnm - case(1) - ! "time": c⋅t + case (1) ! σ=ct, time + ! denom = -ω⋅∂Λ/∂ω + ! Proof: The Hamilton equations are + ! + ! dx̅/dt = + ∂H/∂k̅ (H=ℏΩ, p=ℏk̅) + ! dp̅/dt = - ∂H/∂x̅ + ! + ! where Ω(x̅, k̅) is implicitely defined by the dispersion + ! D(k̅, ω, x̅) = 0. By differentiating the latter we get: + ! + ! ∂Ω/∂x̅ = - ∂D/∂x̅ / ∂D/∂ω = - ∂Λ/∂x̅ / ∂Λ/∂ω + ! ∂Ω/∂k̅ = - ∂D/∂k̅ / ∂D/∂ω = - ∂Λ/∂k̅ / ∂Λ/∂ω + ! + ! with Λ=D/k₀, k₀=ω/c. Finally, substituting k̅=k₀N̅: + ! + ! dx̅/d(ct) = + ∂Λ/∂N̅ / (-ω⋅∂Λ/∂ω) + ! dN̅/d(ct) = - ∂Λ/∂x̅ / (-ω⋅∂Λ/∂ω) + ! denom = -derdom - case(2) - ! real eikonal: S_R + case (2) ! s=S_R, phase + ! denom = N̅⋅∂Λ/∂N̅ + ! Note: This parametrises the curve by the value + ! of the (real) phase, with the wave frozen in time. + ! Proof: By definition N̅ = k̅/k₀ = ∇S. Using the gradient + ! theorem on the ray curve parametrised by the arclength + ! + ! S(s) = ∫₀^x̅(s) N̅⋅dl̅ = ∫₀^s N̅(σ)⋅dx̅/ds ds + ! + ! Differentiating gives dS = N̅⋅dx̅/ds ds, so + ! + ! dS = N̅⋅∂Λ/∂N̅ / |∂Λ/∂N̅| ds ⇒ + ! + ! dx̅/dS = + ∂Λ/∂N̅ / N̅⋅∂Λ/∂N̅ + ! dN̅/dS = - ∂Λ/∂x̅ / N̅⋅∂Λ/∂N̅ + ! denom = dot_product(anv, derdnv) end select - ! coefficient for integration in s - ! ds/dst, where st is the integration variable - dersdst = derdnm/denom + ! Jacobian ds/dσ, where s is the arclength, + ! for computing integrals in dσ, like: + ! + ! τ = ∫α(s)ds = ∫α(σ)(ds/dσ)dσ + ! + if (present(dersdst)) dersdst = derdnm/denom + + ! |∂Λ/∂N̅|: used in the α computation (see alpha_effj) + if (present(derdnm_)) derdnm_ = derdnm ! r.h.s. vector - dery(1:3) = derdnv(:)/denom ! ∂Λ/∂N̅ + dery(1:3) = derdnv(:)/denom ! +∂Λ/∂N̅ dery(4:6) = -derdxv(:)/denom ! -∂Λ/∂x̅ - ! dispersion relation - ddr = an2 - an2s - igrad*(gr2 - half*bdotgr**2*fdia) ! real part - ddi = dot_product(derdnv, dgr) ! imaginary part + if (present(ddr) .or. present(ddi)) then + ! Dispersion relation (geometric optics) + ! ddr ← Λ = N² - N²s(X,Y,N∥) = 0 + an2s = one - xg - half*xg*yg2*(one + anpl2 + sox*del)/duh + ddr = an2 - an2s + end if + + if (present(ddr) .and. igrad > 0) 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 + ! 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 + end if end subroutine disp_deriv