istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

src/gray_core.f90: make some {disp,plas}_deriv outputs optional

Browse Source 
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
This commit is contained in:
Michele Guerini Rocco 2023-03-29 22:13:43 +02:00
parent 658389f586
commit c7d0d8370c
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
1 changed files with 243 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

344
src/gray_core.f90
View File

@ -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