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

src/equilibrium.f90: use enums

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-01-30 10:14:24 +01:00
parent fac0c6ded8
commit 7c5b443847
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
4 changed files with 380 additions and 329 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

621
src/equilibrium.f90
View File

@ -8,6 +8,8 @@
module equilibrium module equilibrium
use const_and_precisions, only : wp_ use const_and_precisions, only : wp_
use splines, only : spline_simple, spline_1d, spline_2d, linear_1d use splines, only : spline_simple, spline_1d, spline_2d, linear_1d
use gray_params, only : EQ_VACUUM, EQ_ANALYTICAL, &
EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
implicit none implicit none
@ -360,36 +362,34 @@ contains
! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0. ! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0.
! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ). ! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ).
if (iequil < 2) then select case(iequil)
! Analytical model case (EQ_ANALYTICAL) ! Analytical model
! Apply signs
if (params%sgni /= 0) then
model%q0 = sign(model%q0, -params%sgni*one)
model%q1 = sign(model%q1, -params%sgni*one)
end if
if (params%sgnb /= 0) then
model%B0 = sign(model%B0, +params%sgnb*one)
end if
! Rescale
model%B0 = model%B0 * params%factb
! Apply signs case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) ! Numeric data
if (params%sgni /= 0) then ! Apply signs
model%q0 = sign(model%q0, -params%sgni*one) if (params%sgni /= 0) &
model%q1 = sign(model%q1, -params%sgni*one) data%psia = sign(data%psia, -params%sgni*one)
end if if (params%sgnb /= 0) &
if (params%sgnb /= 0) then data%fpol = sign(data%fpol, +params%sgnb*one)
model%B0 = sign(model%B0, +params%sgnb*one) ! Rescale
end if data%psia = data%psia * params%factb
! Rescale data%fpol = data%fpol * params%factb
model%B0 = model%B0 * params%factb
else
! Numeric data
! Apply signs ! Compute the signs to be shown in the outputs header when cocos≠0,10.
if (params%sgni /= 0) & ! Note: In these cases the values sgni,sgnb from gray.ini are unused.
data%psia = sign(data%psia, -params%sgni*one) params%sgni = int(sign(one, -data%psia))
if (params%sgnb /= 0) & params%sgnb = int(sign(one, +data%fpol(size(data%fpol))))
data%fpol = sign(data%fpol, +params%sgnb*one) end select
! Rescale
data%psia = data%psia * params%factb
data%fpol = data%fpol * params%factb
! Compute the signs to be shown in the outputs header when cocos≠0,10.
! Note: In these cases the values sgni,sgnb from gray.ini are unused.
params%sgni = int(sign(one, -data%psia))
params%sgnb = int(sign(one, +data%fpol(size(data%fpol))))
end if
end subroutine scale_equil end subroutine scale_equil
@ -398,7 +398,7 @@ contains
! in their respective global variables, see the top of this file. ! in their respective global variables, see the top of this file.
use const_and_precisions, only : zero, one use const_and_precisions, only : zero, one
use gray_params, only : equilibrium_parameters, equilibrium_data use gray_params, only : equilibrium_parameters, equilibrium_data
use gray_params, only : iequil use gray_params, only : iequil, X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
use utils, only : vmaxmin, vmaxmini, inside use utils, only : vmaxmin, vmaxmini, inside
use logger, only : log_info use logger, only : log_info
@ -434,77 +434,80 @@ contains
! Spline interpolation of ψ(R, z) ! Spline interpolation of ψ(R, z)
if (iequil>2) then select case (iequil)
! data valid only inside boundary (data%psin=0 outside), e.g. source==ESCO
! presence of boundary anticipated here to filter invalid data
nbnd = min(size(data%rbnd), size(data%zbnd))
! allocate knots and spline coefficients arrays case (EQ_EQDSK_PARTIAL)
if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x) ! Data valid only inside boundary (data%psin=0 outside),
if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y) ! presence of boundary anticipated here to filter invalid data
if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs) nbnd = min(size(data%rbnd), size(data%zbnd))
allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest))
allocate(psi_spline%coeffs(nrz))
! determine number of valid grid points ! allocate knots and spline coefficients arrays
nrz=0 if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x)
do j=1,nz if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y)
do i=1,nr if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs)
if (nbnd.gt.0) then allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest))
if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle allocate(psi_spline%coeffs(nrz))
else
if(data%psin(i,j).le.0.0d0) cycle ! determine number of valid grid points
end if nrz=0
nrz=nrz+1 do j=1,nz
do i=1,nr
if (nbnd.gt.0) then
if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle
else
if(data%psin(i,j).le.0.0d0) cycle
end if
nrz=nrz+1
end do
end do end do
end do
! store valid data ! store valid data
allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz)) allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz))
ij=0 ij=0
do j=1,nz do j=1,nz
do i=1,nr do i=1,nr
if (nbnd.gt.0) then if (nbnd.gt.0) then
if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle if(.not.inside(data%rbnd,data%zbnd,data%rv(i),data%zv(j))) cycle
else else
if(data%psin(i,j).le.0.0d0) cycle if(data%psin(i,j).le.0.0d0) cycle
end if end if
ij=ij+1 ij=ij+1
rv1d(ij)=data%rv(i) rv1d(ij)=data%rv(i)
zv1d(ij)=data%zv(j) zv1d(ij)=data%zv(j)
fvpsi(ij)=data%psin(i,j) fvpsi(ij)=data%psin(i,j)
wf(ij)=1.0d0 wf(ij)=1.0d0
end do
end do end do
end do
! Fit as a scattered set of points ! Fit as a scattered set of points
! use reduced number of knots to limit memory comsumption ? ! use reduced number of knots to limit memory comsumption ?
psi_spline%nknots_x=nr/4+4
psi_spline%nknots_y=nz/4+4
tension = params%ssplps
call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
rmnm, rmxm, zmnm, zmxm, &
psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, err)
! if failed, re-fit with an interpolating spline (zero tension)
if(err == -1) then
err = 0
tension = 0
psi_spline%nknots_x=nr/4+4 psi_spline%nknots_x=nr/4+4
psi_spline%nknots_y=nz/4+4 psi_spline%nknots_y=nz/4+4
tension = params%ssplps
call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, & call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
rmnm, rmxm, zmnm, zmxm, & rmnm, rmxm, zmnm, zmxm, &
psi_spline%knots_x, psi_spline%nknots_x, & psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, & psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, err) psi_spline%coeffs, err)
end if
deallocate(rv1d, zv1d, wf, fvpsi) ! if failed, re-fit with an interpolating spline (zero tension)
! reset nrz to the total number of grid points for next allocations if(err == -1) then
nrz = nr*nz err = 0
else tension = 0
! iequil==2: data are valid on the full R,z grid psi_spline%nknots_x=nr/4+4
psi_spline%nknots_y=nz/4+4
call scatterspl(rv1d, zv1d, fvpsi, wf, nrz, kspl, tension, &
rmnm, rmxm, zmnm, zmxm, &
psi_spline%knots_x, psi_spline%nknots_x, &
psi_spline%knots_y, psi_spline%nknots_y, &
psi_spline%coeffs, err)
end if
deallocate(rv1d, zv1d, wf, fvpsi)
! reset nrz to the total number of grid points for next allocations
nrz = nr*nz
case (EQ_EQDSK_FULL)
! Data are valid on the full R,z grid
! reshape 2D ψ array to 1D (transposed) ! reshape 2D ψ array to 1D (transposed)
allocate(fvpsi(nrz)) allocate(fvpsi(nrz))
@ -523,7 +526,7 @@ contains
err = 0 err = 0
end if end if
deallocate(fvpsi) deallocate(fvpsi)
end if end select
if (err /= 0) then if (err /= 0) then
err = 2 err = 2
@ -592,11 +595,11 @@ contains
'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp 'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
call log_info(msg, mod='equilibrium', proc='set_equil_spline') call log_info(msg, mod='equilibrium', proc='set_equil_spline')
! search for X-point if params%ixp /= 0 ! Search for X-point
ixploc = params%ixp ixploc = params%ixp
if(ixploc/=0) then
if(ixploc<0) then select case (ixploc)
case (X_AT_BOTTOM)
call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info) call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
if(psinxptmp/=-1.0_wp_) then if(psinxptmp/=-1.0_wp_) then
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
@ -611,7 +614,8 @@ contains
else else
ixploc=0 ixploc=0
end if end if
else
case (X_AT_TOP)
call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info) call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info)
if(psinxptmp.ne.-1.0_wp_) then if(psinxptmp.ne.-1.0_wp_) then
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') & write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
@ -626,26 +630,22 @@ contains
else else
ixploc=0 ixploc=0
end if end if
end if
end if
if (ixploc==0) then case (X_IS_MISSING)
psinop=psinoptmp psinop=psinoptmp
psiant=one-psinop psiant=one-psinop
! Find upper horizontal tangent point
! Find upper horizontal tangent point call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info) zbsup=z1
zbsup=z1 rbsup=r1
rbsup=r1 ! Find lower horizontal tangent point
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
! Find lower horizontal tangent point zbinf=z1
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info) rbinf=r1
zbinf=z1 write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
rbinf=r1 'r', rbinf, rbsup, 'z', zbinf, zbsup
write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') & call log_info(msg, mod='equilibrium', proc='set_equil_spline')
'r', rbinf, rbsup, 'z', zbinf, zbsup end select
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
end if
! Adjust all the B-spline coefficients ! Adjust all the B-spline coefficients
! Note: since ψ_n(R,z) = Σ_ij c_ij B_i(R)B_j(z), to correct ψ_n ! Note: since ψ_n(R,z) = Σ_ij c_ij B_i(R)B_j(z), to correct ψ_n
@ -908,128 +908,131 @@ contains
real(wp_) :: dqdr, dqdz ! ∇q real(wp_) :: dqdr, dqdz ! ∇q
real(wp_) :: dphidr2, ddphidr2dr2 ! dΦ_n/d(r²), d²Φ_n/d(r²)² real(wp_) :: dphidr2, ddphidr2dr2 ! dΦ_n/d(r²), d²Φ_n/d(r²)²
if (iequil < 2) then ! Values for vacuum/outside the domain
! Analytical model if (present(psi_n)) psi_n = -1
! if (present(dpsidr)) dpsidr = 0
! The normalised poloidal flux ψ_n(R, z) is computed as follows: if (present(dpsidz)) dpsidz = 0
! 1. ψ_n = ρ_p² if (present(ddpsidrr)) ddpsidrr = 0
! 2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ) if (present(ddpsidzz)) ddpsidzz = 0
! 3. ρ_t = √Φ_n if (present(ddpsidrz)) ddpsidrz = 0
! 4. Φ_n = Φ(r)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R
! through a circular surface
! 5. r = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius
r_g = hypot(R - model%R0, z - model%z0)
! The exact flux of the toroidal field B_φ = B₀R₀/R is: select case (iequil)
! case (EQ_ANALYTICAL)
! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²). ! Analytical model
! !
! Notes: ! The normalised poloidal flux ψ_n(R, z) is computed as follows:
! 1. the function Φ(r) is defined for r≤R₀ only. ! 1. ψ_n = ρ_p²
! 2. as r → 0, γ → 1, so Φ ~ B₀πr². ! 2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ)
! 3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞. ! 3. ρ_t = √Φ_n
! 4. |B_R|, |B_z| → +-∞. ! 4. Φ_n = Φ(r)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R
! ! through a circular surface
if (r_g > model%R0) then ! 5. r = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius
if (present(psi_n)) psi_n = -1 r_g = hypot(R - model%R0, z - model%z0)
if (present(dpsidr)) dpsidr = 0
if (present(dpsidz)) dpsidz = 0
if (present(ddpsidrr)) ddpsidrr = 0
if (present(ddpsidzz)) ddpsidzz = 0
if (present(ddpsidrz)) ddpsidrz = 0
return
end if
gamma = 1 / sqrt(1 - (r_g/model%R0)**2) ! The exact flux of the toroidal field B_φ = B₀R₀/R is:
phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge !
rho_t = sqrt(phi_n) ! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²).
rho_p = frhopol(rho_t) !
! Notes:
! 1. the function Φ(r) is defined for r≤R₀ only.
! 2. as r → 0, γ → 1, so Φ ~ B₀πr².
! 3. as r → 1⁻, Φ → 2B₀πr² but dΦ/dr → -∞.
! 4. |B_R|, |B_z| → +-∞.
!
if (r_g > model%R0) then
if (present(psi_n)) psi_n = -1
if (present(dpsidr)) dpsidr = 0
if (present(dpsidz)) dpsidz = 0
if (present(ddpsidrr)) ddpsidrr = 0
if (present(ddpsidzz)) ddpsidzz = 0
if (present(ddpsidrz)) ddpsidrz = 0
return
end if
! For ∇Φ_n and ∇∇Φ_n we also need: gamma = 1 / sqrt(1 - (r_g/model%R0)**2)
! phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge
! ∂Φ∂(r²) = B₀π γ(r) rho_t = sqrt(phi_n)
! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²) rho_p = frhopol(rho_t)
!
dphidr2 = model%B0 * pi * gamma / phitedge
ddphidr2dr2 = model%B0 * pi * gamma**3/(2 * model%R0**2) / phitedge
! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²) ! For ∇Φ_n and ∇∇Φ_n we also need:
! where ∇(r²) = 2[(R-R₀), (z-z₀)] !
dphidr = dphidr2 * 2*(R - model%R0) ! ∂Φ∂(r²) = B₀π γ(r)
dphidz = dphidr2 * 2*(z - model%z0) ! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²)
!
dphidr2 = model%B0 * pi * gamma / phitedge
ddphidr2dr2 = model%B0 * pi * gamma**3/(2 * model%R0**2) / phitedge
! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) ! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²)
! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²) ! where ∇(r²) = 2[(R-R₀), (z-z₀)]
! where ∇∇(r²) = 2I dphidr = dphidr2 * 2*(R - model%R0)
ddphidrdr = ddphidr2dr2 * 4*(R - model%R0)*(R - model%R0) + dphidr2*2 dphidz = dphidr2 * 2*(z - model%z0)
ddphidzdz = ddphidr2dr2 * 4*(z - model%z0)*(z - model%z0) + dphidr2*2
ddphidrdz = ddphidr2dr2 * 4*(R - model%R0)*(z - model%z0)
! ψ_n = ρ_p(ρ_t)² ! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
if (present(psi_n)) psi_n = rho_p**2 ! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
! where ∇∇(r²) = 2I
ddphidrdr = ddphidr2dr2 * 4*(R - model%R0)*(R - model%R0) + dphidr2*2
ddphidzdz = ddphidr2dr2 * 4*(z - model%z0)*(z - model%z0) + dphidr2*2
ddphidrdz = ddphidr2dr2 * 4*(R - model%R0)*(z - model%z0)
! Using the definitions in `frhotor`: ! ψ_n = ρ_p(ρ_t)²
! if (present(psi_n)) psi_n = rho_p**2
! ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
!
! ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ
! = Φ_a/ψ_a 1/2πq
!
! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then:
!
! ∂ψ_n/∂Φ_n = (q₀ + Δq)/q
!
q = model%q0 + (model%q1 - model%q0) * rho_p**model%alpha
dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
dpsidphi = (model%q0 + dq) / q
! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n ! Using the definitions in `frhotor`:
if (present(dpsidr)) dpsidr = dpsidphi * dphidr !
if (present(dpsidz)) dpsidz = dpsidphi * dphidz ! ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
!
! ∂ψ_n/∂Φ_n = Φ_a/ψ_a ∂ψ/∂Φ
! = Φ_a/ψ_a 1/2πq
!
! Using ψ_a = 1/2π Φ_a / (q₀ + Δq), then:
!
! ∂ψ_n/∂Φ_n = (q₀ + Δq)/q
!
q = model%q0 + (model%q1 - model%q0) * rho_p**model%alpha
dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
dpsidphi = (model%q0 + dq) / q
! For the second derivatives: ! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
! if (present(dpsidr)) dpsidr = dpsidphi * dphidr
! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n if (present(dpsidz)) dpsidz = dpsidphi * dphidz
!
! ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q
!
! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have:
!
! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n
! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n.
!
dqdr = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidr
dqdz = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidz
ddpsidphidr = - dpsidphi * dqdr/q
ddpsidphidz = - dpsidphi * dqdz/q
! Combining all of the above: ! For the second derivatives:
! !
! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
! !
if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr ! ∇(∂ψ_n/∂Φ_n) = - (∂ψ_n/∂Φ_n) ∇q/q
if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz !
if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz ! From q(ψ) = q₀ + (q₁-q₀) ψ_n^α/2, we have:
else !
! Numerical data ! ∇q = α/2 (q-q₀) ∇ψ_n/ψ_n
if (inside(psi_domain%R, psi_domain%z, R, z)) then ! = α/2 (q-q₀)/ψ_n (∂ψ_n/∂Φ_n) ∇Φ_n.
! Within the interpolation range !
if (present(psi_n)) psi_n = psi_spline%eval(R, z) dqdr = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidr
if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) dqdz = model%alpha/2 * (model%q1 - model%q0)*rho_p**(model%alpha-2) * dpsidphi * dphidz
if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) ddpsidphidr = - dpsidphi * dqdr/q
if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) ddpsidphidz = - dpsidphi * dqdz/q
if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2)
if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) ! Combining all of the above:
else !
! Outside ! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
if (present(psi_n)) psi_n = -1 !
if (present(dpsidr)) dpsidr = 0 if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr
if (present(dpsidz)) dpsidz = 0 if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz
if (present(ddpsidrr)) ddpsidrr = 0 if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz
if (present(ddpsidzz)) ddpsidzz = 0
if (present(ddpsidrz)) ddpsidrz = 0 case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
end if ! Numerical data
end if if (inside(psi_domain%R, psi_domain%z, R, z)) then
! Within the interpolation range
if (present(psi_n)) psi_n = psi_spline%eval(R, z)
if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0)
if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1)
if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0)
if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2)
if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1)
end if
end select
end subroutine pol_flux end subroutine pol_flux
@ -1043,21 +1046,28 @@ contains
real(wp_), intent(out) :: fpol ! poloidal current real(wp_), intent(out) :: fpol ! poloidal current
real(wp_), intent(out), optional :: dfpol ! derivative real(wp_), intent(out), optional :: dfpol ! derivative
if (iequil < 2) then select case (iequil)
! Analytical model case (EQ_VACUUM)
! F(ψ) = B₀⋅R₀, a constant ! Vacuum, no plasma
fpol = model%B0 * model%R0 fpol = 0
if (present(dfpol)) dfpol = 0
else
! Numerical data
if(psi_n <= 1 .and. psi_n >= 0) then
fpol = fpol_spline%eval(psi_n)
if (present(dfpol)) dfpol = fpol_spline%deriv(psi_n)
else
fpol = fpolas
if (present(dfpol)) dfpol = 0 if (present(dfpol)) dfpol = 0
end if
end if case (EQ_ANALYTICAL)
! Analytical model
! F(ψ) = B₀⋅R₀, a constant
fpol = model%B0 * model%R0
if (present(dfpol)) dfpol = 0
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
! Numerical data
if(psi_n <= 1 .and. psi_n >= 0) then
fpol = fpol_spline%eval(psi_n)
if (present(dfpol)) dfpol = fpol_spline%deriv(psi_n)
else
fpol = fpolas
if (present(dfpol)) dfpol = 0
end if
end select
end subroutine pol_curr end subroutine pol_curr
@ -1069,30 +1079,34 @@ contains
real(wp_), intent(in) :: rho_p real(wp_), intent(in) :: rho_p
real(wp_) :: frhotor real(wp_) :: frhotor
if (iequil < 2) then select case (iequil)
! Analytical model
block case (EQ_ANALYTICAL)
! The change of variable is obtained by integrating ! Analytical model
! block
! q(ψ) = 1/2π ∂Φ/∂ψ ! The change of variable is obtained by integrating
! !
! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t². ! q(ψ) = 1/2π ∂Φ/∂ψ
! The result is: !
! ! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t².
! - ψ_a = 1/2π Φ_a / [q₀ + Δq] ! The result is:
! !
! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)] ! - ψ_a = 1/2π Φ_a / [q₀ + Δq]
! !
! where Δq = (q₁ - q₀)/(α/2 + 1) ! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)]
real(wp_) :: dq !
dq = (model%q1 - model%q0) / (model%alpha/2 + 1) ! where Δq = (q₁ - q₀)/(α/2 + 1)
frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) & real(wp_) :: dq
/ (model%q0 + dq)) dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
end block frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) &
else / (model%q0 + dq))
! Numerical data end block
frhotor = rhot_spline%eval(rho_p)
end if case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
! Numerical data
frhotor = rhot_spline%eval(rho_p)
end select
end function frhotor end function frhotor
@ -1105,26 +1119,32 @@ contains
real(wp_), intent(in) :: rho_t real(wp_), intent(in) :: rho_t
real(wp_) :: frhopol real(wp_) :: frhopol
if (iequil < 2) then select case (iequil)
! Analytical model case (EQ_VACUUM)
block ! Vacuum, no plasma
! In general there is no closed form for ρ_p(ρ_t) in the frhopol = 0
! analytical model, we thus solve numerically the equation
! ρ_t(ρ_p) = ρ_t₀ for ρ_p.
use minpack, only : hybrj1
real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7) case (EQ_ANALYTICAL)
integer :: info ! Analytical model
block
! In general there is no closed form for ρ_p(ρ_t) in the
! analytical model, we thus solve numerically the equation
! ρ_t(ρ_p) = ρ_t₀ for ρ_p.
use minpack, only : hybrj1
rho_p = [rho_t] ! first guess, ρ_p ≈ ρ_t real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7)
call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, & integer :: info
ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7)
frhopol = rho_p(1) rho_p = [rho_t] ! first guess, ρ_p ≈ ρ_t
end block call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, &
else ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7)
! Numerical data frhopol = rho_p(1)
frhopol = rhop_spline%eval(rho_t) end block
end if
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
! Numerical data
frhopol = rhop_spline%eval(rho_t)
end select
contains contains
@ -1162,19 +1182,25 @@ contains
real(wp_), intent(in) :: psin real(wp_), intent(in) :: psin
real(wp_) :: fq real(wp_) :: fq
if (iequil < 2) then select case(iequil)
! Analytical model case (EQ_VACUUM)
! The safety factor is a power law in ρ_p: ! Vacuum, q is undefined
! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α fq = 0
block
real(wp_) :: rho_p case (EQ_ANALYTICAL)
rho_p = sqrt(psin) ! Analytical model
fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha) ! The safety factor is a power law in ρ_p:
end block ! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
else block
! Numerical data real(wp_) :: rho_p
fq = q_spline%eval(psin) rho_p = sqrt(psin)
end if fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha)
end block
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
! Numerical data
fq = q_spline%eval(psin)
end select
end function fq end function fq
@ -1183,6 +1209,7 @@ contains
! (R, z) in cylindrical coordinates ! (R, z) in cylindrical coordinates
! !
! Note: all output arguments are optional. ! Note: all output arguments are optional.
use gray_params, only : iequil
! subroutine arguments ! subroutine arguments
real(wp_), intent(in) :: R, z real(wp_), intent(in) :: R, z
@ -1191,6 +1218,14 @@ contains
! local variables ! local variables
real(wp_) :: psi_n, fpol, dpsidr, dpsidz real(wp_) :: psi_n, fpol, dpsidr, dpsidz
if (iequil == EQ_VACUUM) then
! Vacuum, no plasma nor field
if (present(B_R)) B_R = 0
if (present(B_z)) B_z = 0
if (present(B_phi)) B_phi = 0
return
end if
call pol_flux(R, z, psi_n, dpsidr, dpsidz) call pol_flux(R, z, psi_n, dpsidr, dpsidz)
call pol_curr(psi_n, fpol) call pol_curr(psi_n, fpol)

74
src/gray_core.f90
View File

@ -12,7 +12,7 @@ contains
use coreprofiles, only : temp, fzeff use coreprofiles, only : temp, fzeff
use dispersion, only : expinit use dispersion, only : expinit
use gray_params, only : gray_parameters, gray_data, gray_results, & use gray_params, only : gray_parameters, gray_data, gray_results, &
print_parameters print_parameters, EQ_VACUUM
use beams, only : xgygcoeff, launchangles2n use beams, only : xgygcoeff, launchangles2n
use beamdata, only : pweight, rayi2jk use beamdata, only : pweight, rayi2jk
use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd use gray_errors, only : is_critical, print_err_raytracing, print_err_ecrh_cd
@ -100,12 +100,14 @@ contains
! Initialise the dispersion module ! Initialise the dispersion module
if(params%ecrh_cd%iwarm > 1) call expinit if(params%ecrh_cd%iwarm > 1) call expinit
! Initialise the magsurf_data module if (params%equilibrium%iequil /= EQ_VACUUM) then
call flux_average ! requires frhotor for dadrhot,dvdrhot ! Initialise the magsurf_data module
call flux_average ! requires frhotor for dadrhot,dvdrhot
! Initialise the output profiles ! Initialise the output profiles
call pec_init(params%output%ipec, rhout) call pec_init(params%output%ipec, rhout)
nnd = size(rhop_tab) ! number of radial profile points nnd = size(rhop_tab) ! number of radial profile points
end if
call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, & call alloc_multipass(nnd, iwait, iroff, iop, iow, yynext, yypnext, yw0, ypw0, stnext, &
stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, & stv, p0ray, taus, tau1, etau1, cpls, cpl1, lgcpl1, jphi_beam, &
@ -138,11 +140,13 @@ contains
! print Btot=Bres ! print Btot=Bres
! print ne, Te, q, Jphi versus psi, rhop, rhot ! print ne, Te, q, Jphi versus psi, rhop, rhot
call print_bres(bres) if (params%equilibrium%iequil /= EQ_VACUUM) then
call print_prof(params%profiles) call print_bres(bres)
call print_maps(bres, xgcn, & call print_prof(params%profiles)
norm2(params%antenna%pos(1:2)) * 0.01_wp_, & call print_maps(bres, xgcn, &
sin(params%antenna%beta*degree)) norm2(params%antenna%pos(1:2)) * 0.01_wp_, &
sin(params%antenna%beta*degree))
end if
! ========= pre-proc prints END ========= ! ========= pre-proc prints END =========
! =========== main loop BEGIN =========== ! =========== main loop BEGIN ===========
@ -518,9 +522,11 @@ contains
icd_beam = sum(ccci(:,i)) icd_beam = sum(ccci(:,i))
call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print call vmaxmin(tau0,params%raytracing%nray,taumn,taumx) ! taumn,taumx for print
! compute power and current density profiles for all rays if (params%equilibrium%iequil /= EQ_VACUUM) then
call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,dpdv_beam,jphi_beam,jcd_beam, & ! compute power and current density profiles for all rays
pins_beam,currins_beam) call spec(psjki,ppabs,ccci,iiv,pabs_beam,icd_beam,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 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 icd_pass(iox) = icd_pass(iox) + icd_beam
@ -563,17 +569,20 @@ contains
end if end if
end if end if
call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam,dpdv_beam,currins_beam, & if (params%equilibrium%iequil /= EQ_VACUUM) then
pins_beam,ip) ! *print power and current density profiles for current beam call print_pec(rhop_tab,rhot_tab,jphi_beam,jcd_beam, &
dpdv_beam,currins_beam,pins_beam,ip) ! *print power and current density profiles for current beam
call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam,jphi_beam, & call postproc_profiles(pabs_beam,icd_beam,rhot_tab,dpdv_beam, &
rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip,rhotp,drhotp,rhotj, & jphi_beam, rhotpav,drhotpav,rhotjava,drhotjava,dpdvp,jphip, &
drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam rhotp,drhotp,rhotj,drhotj,dpdvmx,jphimx,ratjamx,ratjbmx) ! *compute profiles width for current beam
call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav, &
rhotjava,drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj, &
drhotp,drhotj,ratjamx,ratjbmx,stv(1),psipv(index_rt), &
chipv(index_rt),index_rt,sum(p0ray),cpl_beam1,cpl_beam2) ! *print 0D results for current beam
end if
call print_finals(pabs_beam,icd_beam,dpdvp,jphip,rhotpav,rhotjava, &
drhotpav,drhotjava,dpdvmx,jphimx,rhotp,rhotj,drhotp,drhotj,ratjamx, &
ratjbmx,stv(1),psipv(index_rt),chipv(index_rt),index_rt,sum(p0ray), &
cpl_beam1,cpl_beam2) ! *print 0D results for current beam
! ============ post-proc END ============ ! ============ post-proc END ============
end do beam_loop end do beam_loop
@ -1872,13 +1881,20 @@ contains
bv(1)=br*csphi-bphi*snphi bv(1)=br*csphi-bphi*snphi
bv(2)=br*snphi+bphi*csphi bv(2)=br*snphi+bphi*csphi
bv(3)=bz bv(3)=bz
call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk))
if (jk == 1) then if (norm2(bv) > 0) then
call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv) call pol_limit(anv,bv,bres,sox,ext0(jk),eyt0(jk))
call polellipse(qq,uu,vv,psipol0,chipol0)
psipol0=psipol0/degree ! convert from rad to degree if (jk == 1) then
chipol0=chipol0/degree call stokes_ce(ext0(jk),eyt0(jk),qq,uu,vv)
call polellipse(qq,uu,vv,psipol0,chipol0)
psipol0=psipol0/degree ! convert from rad to degree
chipol0=chipol0/degree
end if
else
! X/O mode are undefined if B=0
psipol0 = 0
chipol0 = 0
end if end if
else else
call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv) call stokes_ell(chipol0*degree,psipol0*degree,qq,uu,vv)

10
src/main.f90
View File

@ -242,7 +242,7 @@ contains
! Set global variables ! Set global variables
select case (params%equilibrium%iequil) select case (params%equilibrium%iequil)
case (EQ_VACUUM, EQ_ANALYTICAL) case (EQ_ANALYTICAL)
call set_equil_an call set_equil_an
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL) case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
@ -387,7 +387,7 @@ contains
EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
use utils, only : range2rect use utils, only : range2rect
use limiter, only : limiter_set_globals=>set_globals use limiter, only : limiter_set_globals=>set_globals
use const_and_precisions, only : cm use const_and_precisions, only : cm, comp_huge
! subroutine arguments ! subroutine arguments
type(gray_parameters), intent(inout) :: params type(gray_parameters), intent(inout) :: params
@ -415,7 +415,11 @@ contains
! Max radius, either due to the plasma extent or equilibrium grid ! Max radius, either due to the plasma extent or equilibrium grid
select case (params%equilibrium%iequil) select case (params%equilibrium%iequil)
case (EQ_VACUUM, EQ_ANALYTICAL) case (EQ_VACUUM)
! Use a very large R, ~ unbounded
R_max = comp_huge
case (EQ_ANALYTICAL)
! use R₀+a ! use R₀+a
block block
use equilibrium, only : model use equilibrium, only : model

4
tests/11-vacuum/inputs/equilibrium.txt
View File

@ -1,4 +0,0 @@
0.0 0.0 0.0 : rr0m,zr0m,rpam ! rhot[m] = min(sqrt((r-rr0m)**2+(z-zr0m)**2), rpam); tor flux phi = pi*b0*rhot**2
0.0 : b0 ! Bphi[T] @ rr0m[m]
1 3 2 : q0, qa, alq ! q = q0 + (qa-q0)*sqrt(psin)**alq
0 : nlim ! number of points in first wall (rlim,zlim) polygon