72eb224568
1. Local variables are automatically deallocated when they go out of scope. 2. When calling exit() during CLI processing some stuff wasn't being deallocated, but it doesnt matter because the OS does it anyway. So, get rid of it entirely.
2048 lines
72 KiB
Fortran
2048 lines
72 KiB
Fortran
! This module handles the loading, interpolation and evaluation of the
|
||
! MHD equilibrium data (poloidal current function, safety factor,
|
||
! poloidal flux, magnetic field, etc.)
|
||
!
|
||
module gray_equil
|
||
use const_and_precisions, only : wp_, comp_huge
|
||
use splines, only : spline_simple, spline_1d, spline_2d, linear_1d
|
||
use types, only : contour
|
||
|
||
implicit none
|
||
|
||
! macro for suppressing unused variable warnings
|
||
# define unused(x) associate(x => x); end associate
|
||
|
||
type, abstract :: abstract_equil
|
||
! Generic equilibrium interface
|
||
real(wp_) :: psi_a ! Poloidal flux at the edge minus flux on axis, ψ_a
|
||
real(wp_) :: phi_a ! Toroidal flux at the edge (r=a), Φ_a
|
||
real(wp_) :: b_axis ! Value of B_φ at the magnetic axis (used in J_cd def)
|
||
real(wp_) :: b_centre ! Value of B_φ at R_centre (used in Jcd_astra def)
|
||
real(wp_) :: r_centre ! Alternative reference radius for B_φ
|
||
real(wp_) :: sgn_bphi ! Sign of B_φ (>0 means counter-clockwise)
|
||
real(wp_) :: axis(2) ! Magnetic axis position (R₀, z₀)
|
||
real(wp_) :: r_range(2) ! R range of the equilibrium domain
|
||
real(wp_) :: z_range(2) ! z range of the equilibrium domain
|
||
real(wp_) :: z_boundary(2) ! z range of the plasma boundary
|
||
|
||
! Flux average splines (see `flux_average`)
|
||
type(spline_simple) :: spline_area, spline_volume
|
||
type(spline_simple) :: spline_dAdrho_t, spline_dVdrho_t
|
||
type(spline_simple) :: spline_R_J, spline_B_avg, spline_B_min, spline_B_max
|
||
type(spline_simple) :: spline_f_c, spline_q, spline_I_p, spline_J_phi_avg
|
||
type(spline_simple) :: spline_astra_tor, spline_jintrac_tor, spline_paral_tor
|
||
type(spline_2d) :: spline_cd_eff
|
||
|
||
contains
|
||
procedure(pol_flux_sub), deferred :: pol_flux
|
||
procedure(pol_curr_sub), deferred :: pol_curr
|
||
procedure(safety_fun), deferred :: safety
|
||
procedure(rho_conv_fun), deferred :: pol2tor, tor2pol
|
||
procedure(flux_contour_sub), deferred :: flux_contour
|
||
procedure :: flux_average
|
||
procedure :: b_field
|
||
procedure :: tor_curr
|
||
procedure :: init_flux_splines
|
||
end type
|
||
|
||
abstract interface
|
||
pure subroutine pol_flux_sub(self, R, z, psi_n, dpsidr, dpsidz, &
|
||
ddpsidrr, ddpsidzz, ddpsidrz)
|
||
! Computes the normalised poloidal flux ψ_n and its
|
||
! derivatives wrt (R, z) up to the second order.
|
||
!
|
||
! Note: all output arguments are optional.
|
||
import :: abstract_equil, wp_
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_), intent(out), optional :: &
|
||
psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz
|
||
end subroutine pol_flux_sub
|
||
|
||
pure subroutine pol_curr_sub(self, psi_n, fpol, dfpol)
|
||
! Computes the poloidal current function F(ψ_n)
|
||
! and (optionally) its derivative dF/dψ_n given ψ_n.
|
||
import :: abstract_equil, wp_
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n ! normalised poloidal flux
|
||
real(wp_), intent(out) :: fpol ! poloidal current
|
||
real(wp_), intent(out), optional :: dfpol ! derivative
|
||
end subroutine pol_curr_sub
|
||
|
||
pure function safety_fun(self, psi_n) result(q)
|
||
! Computes the safety factor q as a function of the
|
||
! normalised poloidal flux ψ_n.
|
||
!
|
||
! Note: this returns the absolute value of q.
|
||
import :: abstract_equil, wp_
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_) :: q
|
||
end function safety_fun
|
||
|
||
pure function rho_conv_fun(self, rho_in) result(rho_out)
|
||
! Converts between poloidal (ρ_p) and toroidal (ρ_t) normalised radius
|
||
import :: abstract_equil, wp_
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_in
|
||
real(wp_) :: rho_out
|
||
end function rho_conv_fun
|
||
|
||
subroutine flux_contour_sub(self, psi0, R, z, &
|
||
R_hi, z_hi, R_lo, z_lo)
|
||
! Computes a contour of the ψ(R,z)=ψ₀ flux surface.
|
||
! Notes:
|
||
! - R,z are the contour points ordered clockwise
|
||
! - (R,z)_hi and (R,z)_lo are a guess for the higher and lower
|
||
! horizontal-tangent points of the contour. These variables
|
||
! will be updated with the exact value on success.
|
||
import :: abstract_equil, wp_
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi0
|
||
real(wp_), intent(out) :: R(:), z(:)
|
||
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
|
||
end subroutine flux_contour_sub
|
||
end interface
|
||
|
||
|
||
type, extends(abstract_equil) :: analytic_equil
|
||
! Analytical equilibrium
|
||
private
|
||
real(wp_) :: q0 ! Safety factor at the magnetic axis
|
||
real(wp_) :: q1 ! Safety factor at the edge
|
||
real(wp_) :: alpha ! Exponent for the q(ρ_p) power law
|
||
real(wp_) :: R0 ! R of the magnetic axis (m)
|
||
real(wp_) :: z0 ! z of the magnetic axis (m)
|
||
real(wp_) :: a ! Minor radius (m)
|
||
real(wp_) :: B0 ! Magnetic field at the magnetic axis (T)
|
||
contains
|
||
procedure :: pol_flux => analytic_pol_flux
|
||
procedure :: pol_curr => analytic_pol_curr
|
||
procedure :: safety => analytic_safety
|
||
procedure :: tor2pol => analytic_tor2pol
|
||
procedure :: pol2tor => analytic_pol2tor
|
||
procedure :: flux_contour => analytic_flux_contour
|
||
end type
|
||
|
||
|
||
type, extends(abstract_equil) :: vacuum
|
||
! Vacuum
|
||
contains
|
||
procedure :: pol_flux => vacuum_pol_flux
|
||
procedure :: pol_curr => vacuum_pol_curr
|
||
procedure :: safety => vacuum_safety
|
||
procedure :: tor2pol => vacuum_conv
|
||
procedure :: pol2tor => vacuum_conv
|
||
procedure :: flux_contour => vacuum_flux_contour
|
||
end type
|
||
|
||
! Interface for custom type constructor
|
||
interface vacuum
|
||
procedure :: vacuum_init
|
||
end interface
|
||
|
||
|
||
type, extends(abstract_equil) :: numeric_equil
|
||
! Numerical equilibrium
|
||
private
|
||
real(wp_) :: fpol_a ! Poloidal current at the edge (r=a), F_a
|
||
! Splines
|
||
type(spline_2d) :: psi_spline
|
||
type(contour) :: psi_domain
|
||
type(spline_1d) :: fpol_spline
|
||
type(spline_simple) :: q_spline
|
||
type(linear_1d) :: rhop_spline, rhot_spline
|
||
contains
|
||
procedure :: pol_flux => numeric_pol_flux
|
||
procedure :: pol_curr => numeric_pol_curr
|
||
procedure :: safety => numeric_safety
|
||
procedure :: tor2pol => numeric_tor2pol
|
||
procedure :: pol2tor => numeric_pol2tor
|
||
procedure :: flux_contour => numeric_flux_contour
|
||
procedure :: init => numeric_init
|
||
procedure :: find_ox_point
|
||
procedure :: find_htg_point
|
||
end type
|
||
|
||
|
||
type eqdsk_data
|
||
! MHD equilibrium data from G-EQDSK file format
|
||
real(wp_), allocatable :: grid_r(:) ! R values of the uniform grid
|
||
real(wp_), allocatable :: grid_z(:) ! z values of the uniform grid
|
||
real(wp_), allocatable :: fpol(:) ! Poloidal current function, F(ψ_n)
|
||
real(wp_), allocatable :: q(:) ! Safety factor, q(ψ_n)
|
||
real(wp_), allocatable :: psi(:) ! Normalised poloidal flux, ψ_n(R)
|
||
real(wp_), allocatable :: psi_map(:,:) ! Normalised poloidal flux 2D map, ψ(R, z)
|
||
real(wp_) :: psi_a ! Poloidal flux at the edge minus flux on axis, ψ_a
|
||
real(wp_) :: r_ref ! Reference R₀ (B = B₀R₀/R without the plasma)
|
||
real(wp_) :: axis(2) ! Magnetic axis position (R₀, z₀)
|
||
type(contour) :: limiter ! limiter contour (wall)
|
||
type(contour) :: boundary ! boundary contour (plasma)
|
||
end type
|
||
|
||
|
||
private
|
||
public abstract_equil ! The abstract equilibrium object
|
||
public analytic_equil, numeric_equil, vacuum ! Implementations
|
||
public load_equil ! To load equilibrium from file
|
||
public eqdsk_data ! G-EQDSK data structure
|
||
public contour ! re-export contours eqdsk_data
|
||
|
||
contains
|
||
|
||
pure subroutine b_field(self, R, z, phi, B_R, B_z, B_phi, B, B_tot)
|
||
! Computes the magnetic field as a function of (R, z)
|
||
! The outputs include cylindrical coordinates, cartesian vector and norm.
|
||
!
|
||
! Note: all output arguments are optional.
|
||
|
||
! subroutine arguments
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_), intent(in), optional :: phi
|
||
real(wp_), intent(out), optional :: B_R, B_z, B_phi
|
||
real(wp_), intent(out), optional :: B(3), B_tot
|
||
|
||
! local variables
|
||
real(wp_) :: psi_n, fpol, dpsidr, dpsidz
|
||
real(wp_) :: B_R_, B_z_, B_phi_
|
||
|
||
call self%pol_flux(R, z, psi_n, dpsidr, dpsidz)
|
||
call self%pol_curr(psi_n, fpol)
|
||
|
||
! The field in cocos=3 is given by
|
||
!
|
||
! B = F(ψ)∇φ + ∇ψ×∇φ.
|
||
!
|
||
! Writing the gradient of ψ=ψ(R,z) as
|
||
!
|
||
! ∇ψ = ∂ψ/∂R ∇R + ∂ψ/∂z ∇z,
|
||
!
|
||
! and carrying out the cross products:
|
||
!
|
||
! B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R
|
||
! = F(ψ)φ^/R - ∂ψ/∂z R^/R + ∂ψ/∂R z^/R
|
||
!
|
||
B_R_ = - 1/R * dpsidz * self%psi_a
|
||
B_z_ = + 1/R * dpsidr * self%psi_a
|
||
B_phi_ = fpol / R
|
||
|
||
if (present(B_R)) B_R = B_R_
|
||
if (present(B_z)) B_z = B_z_
|
||
if (present(B_phi)) B_phi = B_phi_
|
||
|
||
! Convert to Cartesian coordinates
|
||
!
|
||
! Since the unit vectors are given by:
|
||
!
|
||
! { R^ = x^ cosφ + y^ sinφ
|
||
! { φ^ = -x^ sinφ + y^ cosφ,
|
||
!
|
||
! B = φ^ B_φ + R^ B_R + z^ B_z
|
||
! = (B_R cosφ - B_φ sinφ) x^ + (B_R sinφ + B_φ cosφ) y^ + z^ B_Z
|
||
! = x^ B_x + y^ B_y + B_Z z^
|
||
!
|
||
if (present(B)) then
|
||
B(1) = B_R_*cos(phi) - B_phi_*sin(phi)
|
||
B(2) = B_R_*sin(phi) + B_phi_*cos(phi)
|
||
B(3) = B_z_
|
||
end if
|
||
|
||
if (present(B_tot)) B_tot = norm2([B_R_, B_z_, B_phi_])
|
||
end subroutine b_field
|
||
|
||
|
||
pure function tor_curr(self, R, z) result(J_phi)
|
||
! Computes the toroidal current J_φ as a function of (R, z)
|
||
use const_and_precisions, only : mu0_
|
||
|
||
! function arguments
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_) :: J_phi
|
||
|
||
! local variables
|
||
real(wp_) :: dB_Rdz, dB_zdR ! derivatives of B_R, B_z
|
||
real(wp_) :: dpsidr, ddpsidrr, ddpsidzz ! derivatives of ψ_n
|
||
|
||
call self%pol_flux(R, z, dpsidr=dpsidr, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz)
|
||
|
||
! In the usual MHD limit we have ∇×B = μ₀J. Using the
|
||
! curl in cylindrical coords the toroidal current is
|
||
!
|
||
! J_φ = 1/μ₀ (∇×B)_φ = 1/μ₀ [∂B_R/∂z - ∂B_z/∂R].
|
||
!
|
||
! Finally, from B = F(ψ)∇φ + ∇ψ×∇φ we find:
|
||
!
|
||
! B_R = - 1/R ∂ψ/∂z,
|
||
! B_z = + 1/R ∂ψ/∂R,
|
||
!
|
||
! from which:
|
||
!
|
||
! ∂B_R/∂z = - 1/R ∂²ψ/∂z²
|
||
! ∂B_z/∂R = + 1/R ∂²ψ/∂R² - 1/R² ∂ψ/∂R.
|
||
!
|
||
dB_Rdz = - 1/R * ddpsidzz * self%psi_a
|
||
dB_zdR = + 1/R * (ddpsidrr - 1/R * dpsidr) * self%psi_a
|
||
J_phi = 1/mu0_ * (dB_Rdz - dB_zdR)
|
||
end function tor_curr
|
||
|
||
|
||
subroutine flux_average(self, rho_p, area, volume, dAdrho_t, dVdrho_t, &
|
||
R_J, B_avg, B_min, B_max, f_c, q, I_p, J_phi_avg, &
|
||
ratio_astra_tor, ratio_jintrac_tor, ratio_paral_tor)
|
||
! Computes quantities averaged over the flux surface ψ_n(R,z) = ρ_p²
|
||
|
||
! subroutine arguments
|
||
class(abstract_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_p ! normalised poloidal radius, ρ_p
|
||
real(wp_), intent(out), optional :: area, volume ! area and volume enclosed by the flux surface
|
||
real(wp_), intent(out), optional :: dAdrho_t, dVdrho_t ! derivatives of area, volume w.r.t. ρ_t
|
||
real(wp_), intent(out), optional :: R_J ! R_J = ⟨B⟩ / (F(ψ) ⋅ ⟨1/R²⟩), effective J_cd radius
|
||
real(wp_), intent(out), optional :: B_avg ! ⟨B⟩, average magnetic field
|
||
real(wp_), intent(out), optional :: B_min, B_max ! min,max |B| on the flux surface
|
||
real(wp_), intent(out), optional :: f_c ! fraction of circulating particles
|
||
real(wp_), intent(out), optional :: q ! q = 1/2π ∮ B_φ/B_p dl/R, safety factor
|
||
real(wp_), intent(out), optional :: I_p ! I_p = 1/μ₀ ∫ B_p⋅dS_φ, plasma current
|
||
real(wp_), intent(out), optional :: J_phi_avg ! ⟨J_φ⟩, average toroidal current
|
||
real(wp_), intent(out), optional :: ratio_astra_tor ! ratio of J_cd_astra = ⟨J_cd⋅B⟩/B₀ w.r.t. ⟨J_φ⟩
|
||
real(wp_), intent(out), optional :: ratio_jintrac_tor ! ratio of J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩ w.r.t. ⟨J_φ⟩
|
||
real(wp_), intent(out), optional :: ratio_paral_tor ! ratio of Jcd_∥ = ⟨J_cd⋅B/|B|⟩ w.r.t. ⟨J_φ⟩
|
||
|
||
if (present(area)) area = self%spline_area%eval(rho_p)
|
||
if (present(volume)) volume = self%spline_volume%eval(rho_p)
|
||
|
||
if (present(dAdrho_t)) dAdrho_t = self%spline_dAdrho_t%eval(rho_p)
|
||
if (present(dVdrho_t)) dVdrho_t = self%spline_dVdrho_t%eval(rho_p)
|
||
|
||
if (present(R_J)) R_J = self%spline_R_J%eval(rho_p)
|
||
if (present(B_avg)) B_avg = self%spline_B_avg%eval(rho_p)
|
||
if (present(B_min)) B_min = self%spline_B_min%eval(rho_p)
|
||
if (present(B_max)) B_max = self%spline_B_max%eval(rho_p)
|
||
|
||
if (present(f_c)) f_c = self%spline_f_c%eval(rho_p)
|
||
if (present(q)) q = self%spline_q%eval(rho_p)
|
||
if (present(I_p)) I_p = self%spline_I_p%eval(rho_p)
|
||
|
||
if (present(ratio_astra_tor)) ratio_astra_tor = self%spline_astra_tor%eval(rho_p)
|
||
if (present(ratio_jintrac_tor)) ratio_jintrac_tor = self%spline_jintrac_tor%eval(rho_p)
|
||
if (present(ratio_paral_tor)) ratio_paral_tor = self%spline_paral_tor%eval(rho_p)
|
||
if (present(J_phi_avg)) J_phi_avg = self%spline_J_phi_avg%eval(rho_p)
|
||
end subroutine flux_average
|
||
|
||
|
||
!
|
||
! Analytical model
|
||
!
|
||
|
||
pure subroutine analytic_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
|
||
ddpsidrr, ddpsidzz, ddpsidrz)
|
||
use const_and_precisions, only : pi
|
||
|
||
! function arguments
|
||
class(analytic_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_), intent(out), optional :: &
|
||
psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz
|
||
|
||
! local variables
|
||
real(wp_) :: r_g, rho_t, rho_p ! geometric radius, √Φ_n, √ψ_n
|
||
real(wp_) :: gamma ! γ = 1/√(1 - r²/R₀²)
|
||
real(wp_) :: dpsidphi ! (∂ψ_n/∂Φ_n)
|
||
real(wp_) :: ddpsidphidr, ddpsidphidz ! ∇(∂ψ_n/∂Φ_n)
|
||
real(wp_) :: phi_n ! Φ_n
|
||
real(wp_) :: dphidr, dphidz ! ∇Φ_n
|
||
real(wp_) :: ddphidrdr, ddphidzdz ! ∇∇Φ_n
|
||
real(wp_) :: ddphidrdz ! ∂²Φ_n/∂R∂z
|
||
real(wp_) :: q, dq ! q(ρ_p), Δq=(q₁-q₀)/(α/2 + 1)
|
||
real(wp_) :: dqdr, dqdz ! ∇q
|
||
real(wp_) :: dphidr2, ddphidr2dr2 ! dΦ_n/d(r²), d²Φ_n/d(r²)²
|
||
|
||
! The normalised poloidal flux ψ_n(R, z) is computed as follows:
|
||
! 1. ψ_n = ρ_p²
|
||
! 2. ρ_p = ρ_p(ρ_t), using `self%tor2pol`, which in turns uses q(ψ)
|
||
! 3. ρ_t = √Φ_n
|
||
! 4. Φ_n = Φ(r_g)/Φ(a), where Φ(r) is the flux of B_φ=B₀R₀/R
|
||
! through a circular surface
|
||
! 5. r_g = √[(R-R₀)²+(z-z₀)²] is the geometric minor radius
|
||
r_g = hypot(R - self%R0, z - self%z0)
|
||
|
||
! The exact flux of the toroidal field B_φ = B₀R₀/R is:
|
||
!
|
||
! Φ(r) = B₀πr² 2γ/(γ + 1) where γ=1/√(1 - r²/R₀²).
|
||
!
|
||
! 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 > self%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
|
||
|
||
gamma = 1 / sqrt(1 - (r_g/self%R0)**2)
|
||
phi_n = self%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / self%phi_a
|
||
rho_t = sqrt(phi_n)
|
||
rho_p = self%tor2pol(rho_t)
|
||
|
||
! For ∇Φ_n and ∇∇Φ_n we also need:
|
||
!
|
||
! ∂Φ∂(r²) = B₀π γ(r)
|
||
! ∂²Φ∂(r²)² = B₀π γ³(r) / (2 R₀²)
|
||
!
|
||
dphidr2 = self%B0 * pi * gamma / self%phi_a
|
||
ddphidr2dr2 = self%B0 * pi * gamma**3/(2 * self%R0**2) / self%phi_a
|
||
|
||
! ∇Φ_n = ∂Φ_n/∂(r²) ∇(r²)
|
||
! where ∇(r²) = 2[(R-R₀), (z-z₀)]
|
||
dphidr = dphidr2 * 2*(R - self%R0)
|
||
dphidz = dphidr2 * 2*(z - self%z0)
|
||
|
||
! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
|
||
! = ∂²Φ_n/∂(r²)² ∇(r²)∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
|
||
! where ∇∇(r²) = 2I
|
||
ddphidrdr = ddphidr2dr2 * 4*(R - self%R0)*(R - self%R0) + dphidr2*2
|
||
ddphidzdz = ddphidr2dr2 * 4*(z - self%z0)*(z - self%z0) + dphidr2*2
|
||
ddphidrdz = ddphidr2dr2 * 4*(R - self%R0)*(z - self%z0)
|
||
|
||
! ψ_n = ρ_p(ρ_t)²
|
||
if (present(psi_n)) psi_n = rho_p**2
|
||
|
||
! Using the definitions in `frhotor`:
|
||
!
|
||
! ∇ψ_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 = self%q0 + (self%q1 - self%q0) * rho_p**self%alpha
|
||
dq = (self%q1 - self%q0) / (self%alpha/2 + 1)
|
||
dpsidphi = (self%q0 + dq) / q
|
||
|
||
! Using the above, ∇ψ_n = ∂ψ_n/∂Φ_n ∇Φ_n
|
||
if (present(dpsidr)) dpsidr = dpsidphi * dphidr
|
||
if (present(dpsidz)) dpsidz = dpsidphi * dphidz
|
||
|
||
! For the second derivatives:
|
||
!
|
||
! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
|
||
!
|
||
! ∇(∂ψ_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 = self%alpha/2 * (self%q1 - self%q0)*rho_p**(self%alpha-2) * dpsidphi * dphidr
|
||
dqdz = self%alpha/2 * (self%q1 - self%q0)*rho_p**(self%alpha-2) * dpsidphi * dphidz
|
||
ddpsidphidr = - dpsidphi * dqdr/q
|
||
ddpsidphidz = - dpsidphi * dqdz/q
|
||
|
||
! Combining all of the above:
|
||
!
|
||
! ∇∇ψ_n = ∇(∂ψ_n/∂Φ_n) ∇Φ_n + (∂ψ_n/∂Φ_n) ∇∇Φ_n
|
||
!
|
||
if (present(ddpsidrr)) ddpsidrr = ddpsidphidr * dphidr + dpsidphi * ddphidrdr
|
||
if (present(ddpsidzz)) ddpsidzz = ddpsidphidz * dphidz + dpsidphi * ddphidzdz
|
||
if (present(ddpsidrz)) ddpsidrz = ddpsidphidr * dphidz + dpsidphi * ddphidrdz
|
||
end subroutine analytic_pol_flux
|
||
|
||
|
||
pure subroutine analytic_pol_curr(self, psi_n, fpol, dfpol)
|
||
! subroutine arguments
|
||
class(analytic_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_), intent(out) :: fpol ! poloidal current
|
||
real(wp_), intent(out), optional :: dfpol ! its derivative
|
||
|
||
unused(psi_n)
|
||
|
||
! The poloidal current function F(ψ) is just a constant:
|
||
!
|
||
! B_φ = B₀R₀/R φ^ ≡ F(ψ)∇φ ⇒ F(ψ)=B₀R₀
|
||
!
|
||
fpol = self%B0 * self%R0
|
||
if (present(dfpol)) dfpol = 0
|
||
end subroutine analytic_pol_curr
|
||
|
||
|
||
pure function analytic_safety(self, psi_n) result(q)
|
||
! function arguments
|
||
class(analytic_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_) :: q
|
||
|
||
! local variables
|
||
real(wp_) :: rho_p
|
||
|
||
! The safety factor is a power law in ρ_p:
|
||
!
|
||
! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
|
||
!
|
||
rho_p = sqrt(psi_n)
|
||
q = self%q0 + (self%q1 - self%q0) * rho_p**self%alpha
|
||
end function analytic_safety
|
||
|
||
|
||
pure function analytic_pol2tor(self, rho_in) result(rho_out)
|
||
! function arguments
|
||
class(analytic_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_in
|
||
real(wp_) :: rho_out
|
||
|
||
! local variables
|
||
real(wp_) :: dq
|
||
|
||
! The change of variable is obtained by integrating
|
||
!
|
||
! q(ψ) = 1/2π ∂Φ/∂ψ
|
||
!
|
||
! and defining ψ = ψ_a ρ_p², Φ = Φ_a ρ_t².
|
||
! The result is:
|
||
!
|
||
! - ψ_a = 1/2π Φ_a / [q₀ + Δq]
|
||
!
|
||
! - ρ_t = ρ_p √[(q₀ + Δq ρ_p^α)/(q₀ + Δq)]
|
||
!
|
||
! where Δq = (q₁ - q₀)/(α/2 + 1)
|
||
dq = (self%q1 - self%q0) / (self%alpha/2 + 1)
|
||
rho_out = rho_in * sqrt((self%q0 + dq*rho_in**self%alpha) / (self%q0 + dq))
|
||
end function analytic_pol2tor
|
||
|
||
|
||
pure function analytic_tor2pol(self, rho_in) result(rho_out)
|
||
! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
|
||
use const_and_precisions, only : comp_eps
|
||
use minpack, only : hybrj1
|
||
|
||
! function arguments
|
||
class(analytic_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_in
|
||
real(wp_) :: rho_out
|
||
|
||
! local variables
|
||
real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7)
|
||
integer :: info
|
||
|
||
! 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.
|
||
rho_p = [rho_in] ! first guess, ρ_p ≈ ρ_t
|
||
call hybrj1(equation, n=1, x=rho_p, fvec=fvec, fjac=fjac, &
|
||
ldfjac=1, tol=comp_eps, info=info, wa=wa, lwa=7)
|
||
rho_out = rho_p(1)
|
||
|
||
contains
|
||
|
||
pure subroutine equation(n, x, f, df, ldf, flag)
|
||
! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0
|
||
|
||
! optimal step size
|
||
real(wp_), parameter :: e = comp_eps**(1/3.0_wp_)
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: n, ldf, flag
|
||
real(wp_), intent(in) :: x(n)
|
||
real(wp_), intent(inout) :: f(n), df(ldf,n)
|
||
|
||
if (flag == 1) then
|
||
! return f(x)
|
||
f(1) = self%pol2tor(x(1)) - rho_in
|
||
else
|
||
! return f'(x), computed numerically
|
||
if (x(1) - e > 0) then
|
||
df(1,1) = (self%pol2tor(x(1) + e) - self%pol2tor(x(1) - e)) / (2*e)
|
||
else
|
||
! single-sided when close to ρ=0
|
||
df(1,1) = (self%pol2tor(x(1) + e) - self%pol2tor(x(1))) / e
|
||
end if
|
||
end if
|
||
end subroutine
|
||
|
||
end function analytic_tor2pol
|
||
|
||
|
||
pure subroutine analytic_flux_contour(self, psi0, R, z, &
|
||
R_hi, z_hi, R_lo, z_lo)
|
||
use const_and_precisions, only : pi
|
||
|
||
! subroutine arguments
|
||
class(analytic_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi0
|
||
real(wp_), intent(out) :: R(:), z(:)
|
||
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
|
||
|
||
! local variables
|
||
integer :: n, i
|
||
real(wp_) :: r_g ! geometric minor radius
|
||
real(wp_) :: theta ! geometric poloidal angle
|
||
|
||
unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo)
|
||
|
||
n = size(R)
|
||
r_g = sqrt(psi0) * self%a
|
||
theta = 2*pi / (n - 1)
|
||
|
||
do concurrent (i=1:n)
|
||
R(i) = self%R0 + r_g * cos(-theta*(i-1))
|
||
z(i) = self%z0 + r_g * sin(-theta*(i-1))
|
||
end do
|
||
|
||
end subroutine analytic_flux_contour
|
||
|
||
|
||
!
|
||
! Numerical equilibrium
|
||
!
|
||
|
||
pure subroutine numeric_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
|
||
ddpsidrr, ddpsidzz, ddpsidrz)
|
||
! function arguments
|
||
class(numeric_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_), intent(out), optional :: &
|
||
psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz
|
||
|
||
if (self%psi_domain%contains(R, z)) then
|
||
! Within the safety domain
|
||
if (present(psi_n)) psi_n = self%psi_spline%eval(R, z)
|
||
if (present(dpsidr)) dpsidr = self%psi_spline%deriv(R, z, 1, 0)
|
||
if (present(dpsidz)) dpsidz = self%psi_spline%deriv(R, z, 0, 1)
|
||
if (present(ddpsidrr)) ddpsidrr = self%psi_spline%deriv(R, z, 2, 0)
|
||
if (present(ddpsidzz)) ddpsidzz = self%psi_spline%deriv(R, z, 0, 2)
|
||
if (present(ddpsidrz)) ddpsidrz = self%psi_spline%deriv(R, z, 1, 1)
|
||
else
|
||
! Values for outside the domain
|
||
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
|
||
end if
|
||
|
||
end subroutine numeric_pol_flux
|
||
|
||
|
||
pure subroutine numeric_pol_curr(self, psi_n, fpol, dfpol)
|
||
! subroutine arguments
|
||
class(numeric_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_), intent(out) :: fpol ! poloidal current
|
||
real(wp_), intent(out), optional :: dfpol ! its derivative
|
||
|
||
if (psi_n <= 1 .and. psi_n >= 0) then
|
||
! Inside plasma
|
||
fpol = self%fpol_spline%eval(psi_n)
|
||
if (present(dfpol)) dfpol = self%fpol_spline%deriv(psi_n)
|
||
else
|
||
! Outside plasma
|
||
fpol = self%fpol_a
|
||
if (present(dfpol)) dfpol = 0
|
||
end if
|
||
end subroutine numeric_pol_curr
|
||
|
||
|
||
pure function numeric_safety(self, psi_n) result(q)
|
||
! function arguments
|
||
class(numeric_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_) :: q
|
||
|
||
if (psi_n < 1) then
|
||
! Inside plasma
|
||
q = self%q_spline%eval(psi_n)
|
||
else
|
||
! Outside plasma, q is undefined
|
||
q = 0
|
||
end if
|
||
end function numeric_safety
|
||
|
||
|
||
pure function numeric_pol2tor(self, rho_in) result(rho_out)
|
||
! function arguments
|
||
class(numeric_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_in
|
||
real(wp_) :: rho_out
|
||
|
||
rho_out = self%rhot_spline%eval(rho_in)
|
||
end function numeric_pol2tor
|
||
|
||
|
||
pure function numeric_tor2pol(self, rho_in) result(rho_out)
|
||
! function arguments
|
||
class(numeric_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_in
|
||
real(wp_) :: rho_out
|
||
|
||
rho_out = self%rhop_spline%eval(rho_in)
|
||
end function numeric_tor2pol
|
||
|
||
|
||
subroutine numeric_flux_contour(self, psi0, R, z, &
|
||
R_hi, z_hi, R_lo, z_lo)
|
||
use const_and_precisions, only : pi
|
||
use logger, only : log_warning
|
||
use dierckx, only : profil, sproota
|
||
|
||
! subroutine arguments
|
||
class(numeric_equil), intent(in) :: self
|
||
real(wp_), intent(in) :: psi0
|
||
real(wp_), intent(out) :: R(:), z(:)
|
||
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
|
||
|
||
! local variables
|
||
integer :: n, np, i
|
||
integer :: ier, iopt, m
|
||
real(wp_) :: theta
|
||
real(wp_) :: R_hi1, z_hi1, R_lo1, z_lo1, zc
|
||
real(wp_) :: czc(self%psi_spline%nknots_x), sols(3)
|
||
character(256) :: msg
|
||
|
||
! TODO: handle the n even case
|
||
n = size(R)
|
||
np = (n - 1)/2
|
||
theta = pi / np
|
||
|
||
! Find the lowest and highest point in the contour
|
||
call self%find_htg_point(R_hi, z_hi, R_hi1, z_hi1, psi0)
|
||
call self%find_htg_point(R_lo, z_lo, R_lo1, z_lo1, psi0)
|
||
|
||
R(1) = R_lo1
|
||
z(1) = z_lo1
|
||
R(n) = R_lo1
|
||
z(n) = z_lo1
|
||
R(np+1) = R_hi1
|
||
z(np+1) = z_hi1
|
||
|
||
! Slice up ψ(R,z) for each z ∈ [z_lo, z_hi]
|
||
do i = 2, np
|
||
zc = z_lo1 + (z_hi1 - z_lo1) * (1 - cos(theta*(i-1))) / 2
|
||
iopt = 1
|
||
associate (s => self%psi_spline)
|
||
! Compute the cubic spline s(R) = ψ(R,z) at fixed z
|
||
call profil(iopt, s%knots_x, s%nknots_x, s%knots_y, s%nknots_y, &
|
||
s%coeffs, 3, 3, zc, s%nknots_x, czc, ier)
|
||
if (ier > 0) then
|
||
write(msg, '(a, g0)') &
|
||
'when computing ψ(R,z) contour `profil` returned ier=', ier
|
||
call log_warning(msg, mod='gray_equil', proc='numeric_flux_contour')
|
||
end if
|
||
|
||
! Find the solutions of s(R) = ψ₀
|
||
call sproota(psi0, s%knots_x, s%nknots_x, czc, sols, 3, m, ier)
|
||
end associate
|
||
|
||
! Select the physical zeros (in the region where ψ is convex)
|
||
if (self%psi_domain%contains(sols(1), zc)) then
|
||
R(i) = sols(1)
|
||
R(n+1-i) = sols(2)
|
||
else
|
||
R(i) = sols(2)
|
||
R(n+1-i) = sols(3)
|
||
end if
|
||
z(i) = zc
|
||
z(n+1-i) = zc
|
||
end do
|
||
|
||
! Replace the initial guess with the exact values
|
||
R_hi = R_hi1
|
||
z_hi = z_hi1
|
||
R_lo = R_lo1
|
||
z_lo = z_lo1
|
||
end subroutine numeric_flux_contour
|
||
|
||
|
||
subroutine numeric_init(self, params, data, err)
|
||
! Initialises a numeric equilibrium
|
||
use const_and_precisions, only : zero, one
|
||
use gray_params, only : equilibrium_parameters
|
||
use gray_params, only : EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
|
||
use gray_params, only : X_AT_TOP, X_AT_BOTTOM, X_IS_MISSING
|
||
use utils, only : vmaxmin
|
||
use logger, only : log_debug, log_info
|
||
|
||
! subroutine arguments
|
||
class(numeric_equil), intent(out) :: self
|
||
type(equilibrium_parameters), intent(in) :: params
|
||
type(eqdsk_data), intent(in) :: data
|
||
integer, intent(out) :: err
|
||
|
||
! local variables
|
||
integer :: nr, nz, nrz, npsi, nbnd, ibinf, ibsup
|
||
real(wp_) :: psinoptmp, psinxptmp
|
||
real(wp_) :: rbinf, rbsup, zbinf, zbsup, R1, z1
|
||
real(wp_) :: psiant, psinop
|
||
real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi
|
||
integer :: i, j, ij
|
||
character(256) :: msg ! for log messages formatting
|
||
|
||
! compute array sizes
|
||
nr = size(data%grid_r)
|
||
nz = size(data%grid_z)
|
||
nrz = nr*nz
|
||
npsi = size(data%psi)
|
||
|
||
! length in m !!!
|
||
self%r_range = [data%grid_r(1), data%grid_r(nr)]
|
||
self%z_range = [data%grid_z(1), data%grid_z(nz)]
|
||
|
||
call log_debug('computing splines...', mod='gray_equil', proc='numeric_init')
|
||
|
||
! Spline interpolation of ψ(R, z)
|
||
|
||
select case (params%iequil)
|
||
|
||
case (EQ_EQDSK_PARTIAL)
|
||
! Data valid only inside boundary (data%psi=0 outside),
|
||
! presence of boundary anticipated here to filter invalid data
|
||
nbnd = size(data%boundary%R)
|
||
|
||
! allocate knots and spline coefficients arrays
|
||
allocate(self%psi_spline%knots_x(nr + 4), self%psi_spline%knots_y(nz + 4))
|
||
allocate(self%psi_spline%coeffs(nrz))
|
||
|
||
! determine number of valid grid points
|
||
nrz=0
|
||
do j=1,nz
|
||
do i=1,nr
|
||
if (nbnd > 0) then
|
||
if (.not. data%boundary%contains(data%grid_r(i), data%grid_z(j))) cycle
|
||
else
|
||
if (data%psi_map(i,j) <= 0) cycle
|
||
end if
|
||
nrz=nrz+1
|
||
end do
|
||
end do
|
||
|
||
! store valid data
|
||
allocate(rv1d(nrz), zv1d(nrz), fvpsi(nrz))
|
||
ij=0
|
||
do j=1,nz
|
||
do i=1,nr
|
||
if (nbnd > 0) then
|
||
if (.not. data%boundary%contains(data%grid_r(i), data%grid_z(j))) cycle
|
||
else
|
||
if (data%psi_map(i,j) <= 0) cycle
|
||
end if
|
||
ij = ij + 1
|
||
rv1d(ij) = data%grid_r(i)
|
||
zv1d(ij) = data%grid_z(j)
|
||
fvpsi(ij) = data%psi_map(i,j)
|
||
end do
|
||
end do
|
||
|
||
! Fit as a scattered set of points
|
||
! use reduced number of knots to limit memory comsumption ?
|
||
self%psi_spline%nknots_x=nr/4+4
|
||
self%psi_spline%nknots_y=nz/4+4
|
||
call self%psi_spline%init_nonreg( &
|
||
rv1d, zv1d, fvpsi, tension=params%ssplps, &
|
||
range=[self%r_range, self%z_range], err=err)
|
||
|
||
! if failed, re-fit with an interpolating spline (zero tension)
|
||
if (err == -1) then
|
||
err = 0
|
||
self%psi_spline%nknots_x=nr/4+4
|
||
self%psi_spline%nknots_y=nz/4+4
|
||
call self%psi_spline%init_nonreg( &
|
||
rv1d, zv1d, fvpsi, tension=zero, &
|
||
range=[self%r_range, self%z_range], err=err)
|
||
end if
|
||
! 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)
|
||
allocate(fvpsi(nrz))
|
||
fvpsi = reshape(transpose(data%psi_map), [nrz])
|
||
|
||
! compute spline coefficients
|
||
call self%psi_spline%init(data%grid_r, data%grid_z, fvpsi, nr, nz, &
|
||
range=[self%r_range, self%z_range], &
|
||
tension=params%ssplps, err=err)
|
||
|
||
! if failed, re-fit with an interpolating spline (zero tension)
|
||
if (err == -1) then
|
||
call self%psi_spline%init(data%grid_r, data%grid_z, fvpsi, nr, nz, &
|
||
range=[self%r_range, self%z_range], tension=zero)
|
||
err = 0
|
||
end if
|
||
end select
|
||
|
||
if (err /= 0) then
|
||
err = 2
|
||
return
|
||
end if
|
||
|
||
! compute spline coefficients for ψ(R,z) partial derivatives
|
||
call self%psi_spline%init_deriv(nr, nz, 1, 0) ! ∂ψ/∂R
|
||
call self%psi_spline%init_deriv(nr, nz, 0, 1) ! ∂ψ/∂z
|
||
call self%psi_spline%init_deriv(nr, nz, 1, 1) ! ∂²ψ/∂R∂z
|
||
call self%psi_spline%init_deriv(nr, nz, 2, 0) ! ∂²ψ/∂R²
|
||
call self%psi_spline%init_deriv(nr, nz, 0, 2) ! ∂²ψ/∂z²
|
||
|
||
! set the initial ψ(R,z) domain to the grid boundary
|
||
!
|
||
! Note: this is required to bootstrap `flux_pol` calls
|
||
! within this very subroutine.
|
||
self%psi_domain = contour(Rmin=self%r_range(1), Rmax=self%r_range(2), &
|
||
zmin=self%z_range(1), zmax=self%z_range(2))
|
||
|
||
! Spline interpolation of F(ψ)
|
||
|
||
! give a small weight to the last point
|
||
block
|
||
real(wp_) :: w(npsi)
|
||
w(1:npsi-1) = 1
|
||
w(npsi) = 1.0e2_wp_
|
||
call self%fpol_spline%init(data%psi, data%fpol, npsi, weights=w, &
|
||
range=[zero, one], tension=params%ssplf)
|
||
end block
|
||
|
||
! set vacuum value used outside 0 ≤ ψ ≤ 1 range
|
||
self%fpol_a = self%fpol_spline%eval(data%psi(npsi))
|
||
self%sgn_bphi = sign(one, self%fpol_a)
|
||
|
||
! Re-normalize ψ_n after the spline computation
|
||
! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma
|
||
|
||
! Start with un-corrected ψ_n
|
||
self%psi_a = data%psi_a
|
||
psinop = 0 ! ψ_n(O point)
|
||
psiant = 1 ! ψ_n(X point) - ψ_n(O point)
|
||
|
||
! Use provided boundary to set an initial guess
|
||
! for the search of O/X points
|
||
nbnd = size(data%boundary%R)
|
||
if (nbnd > 0) then
|
||
call vmaxmin(data%boundary%z, zbinf, zbsup, ibinf, ibsup)
|
||
rbinf = data%boundary%R(ibinf)
|
||
rbsup = data%boundary%R(ibsup)
|
||
else
|
||
zbinf = data%grid_z(2)
|
||
zbsup = data%grid_z(nz-1)
|
||
rbinf = data%grid_r((nr+1)/2)
|
||
rbsup = rbinf
|
||
end if
|
||
|
||
! Search for exact location of the magnetic axis
|
||
call self%find_ox_point(R0=data%axis(1), z0=data%axis(2), &
|
||
R1=self%axis(1), z1=self%axis(2), psi1=psinoptmp)
|
||
|
||
write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
|
||
'r', self%axis(1), 'z', self%axis(2), 'ψ', psinoptmp
|
||
call log_info(msg, mod='gray_equil', proc='numeric_init')
|
||
|
||
! Search for X-point
|
||
select case (params%ixp)
|
||
case (X_AT_BOTTOM)
|
||
call self%find_ox_point(R0=rbinf, z0=zbinf, R1=R1, z1=z1, psi1=psinxptmp)
|
||
rbinf = z1
|
||
if (psinxptmp /= -1) then
|
||
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
|
||
'R', R1, 'z', z1, 'ψ', psinxptmp
|
||
call log_info(msg, mod='gray_equil', proc='numeric_init')
|
||
psinop = psinoptmp
|
||
psiant = psinxptmp-psinop
|
||
call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbsup)/2, &
|
||
R1=r1, z1=zbsup, psi0=one)
|
||
end if
|
||
|
||
case (X_AT_TOP)
|
||
call self%find_ox_point(R0=rbsup, z0=zbsup, R1=R1, z1=z1, psi1=psinxptmp)
|
||
zbsup = z1
|
||
if (psinxptmp /= -1) then
|
||
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
|
||
'R', r1, 'z', zbsup, 'ψ', psinxptmp
|
||
call log_info(msg, mod='gray_equil', proc='numeric_init')
|
||
psinop = psinoptmp
|
||
psiant = psinxptmp - psinop
|
||
call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbinf)/2, &
|
||
R1=r1, z1=zbinf, psi0=one)
|
||
end if
|
||
|
||
case (X_IS_MISSING)
|
||
psinop = psinoptmp
|
||
psiant = 1 - psinop
|
||
! Find upper horizontal tangent point
|
||
call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbsup)/2, &
|
||
R1=rbsup, z1=zbsup, psi0=one)
|
||
! Find lower horizontal tangent point
|
||
call self%find_htg_point(R0=self%axis(1), z0=(self%axis(2)+zbinf)/2, &
|
||
R1=rbinf, z1=zbinf, psi0=one)
|
||
write (msg, '("X-point not found in", 2(x,a,"∈[",g0.3,",",g0.3,"]"))') &
|
||
'R', rbinf, rbsup, 'z', zbinf, zbsup
|
||
call log_info(msg, mod='gray_equil', proc='numeric_init')
|
||
end select
|
||
|
||
! Correct the spline coefficients: ψ_n → (ψ_n - psinop)/psiant
|
||
call self%psi_spline%transform(1/psiant, -psinop/psiant)
|
||
|
||
! Do the opposite scaling to preserve un-normalised values
|
||
! Note: this is only used for the poloidal magnetic field
|
||
self%psi_a = self%psi_a * psiant
|
||
|
||
! Compute other constants
|
||
call self%pol_curr(zero, self%b_axis)
|
||
self%b_axis = self%b_axis / self%axis(1)
|
||
self%b_centre = self%fpol_a / data%r_ref
|
||
self%r_centre = data%r_ref
|
||
self%z_boundary = [zbinf, zbsup]
|
||
|
||
write (msg, '(2(a,g0.3))') 'B_centre=', self%b_centre, ' B_axis=', self%b_axis
|
||
call log_info(msg, mod='gray_equil', proc='numeric_init')
|
||
|
||
! Compute ρ_p/ρ_t mapping based on the input q profile
|
||
block
|
||
use const_and_precisions, only : pi
|
||
real(wp_), dimension(npsi) :: phi, rho_p, rho_t
|
||
real(wp_) :: dx
|
||
integer :: k
|
||
|
||
call self%q_spline%init(data%psi, data%q, npsi)
|
||
|
||
! Toroidal flux as Φ(ψ) = 2π ∫ q(ψ)dψ
|
||
! Specifically, Φ(ψ_n) = 2π|ψ_a| ∫₀ψ_n q(ψ_n)dψ_n, with ψ_n ∈ [0,1]
|
||
phi(1) = 0
|
||
do k = 1, npsi-1
|
||
dx = self%q_spline%data(k+1) - self%q_spline%data(k)
|
||
phi(k+1) = phi(k) + dx*(self%q_spline%coeffs(k,1) + dx*(self%q_spline%coeffs(k,2)/2 + &
|
||
dx*(self%q_spline%coeffs(k,3)/3 + dx* self%q_spline%coeffs(k,4)/4) ) )
|
||
end do
|
||
|
||
self%phi_a = phi(npsi)
|
||
rho_p = sqrt(data%psi)
|
||
rho_t = sqrt(phi/self%phi_a)
|
||
self%phi_a = 2*pi * abs(self%psi_a) * self%phi_a
|
||
|
||
call self%rhop_spline%init(rho_t, rho_p, size(rho_p))
|
||
call self%rhot_spline%init(rho_p, rho_t, size(rho_t))
|
||
end block
|
||
|
||
call log_debug('splines computed', mod='gray_equil', proc='numeric_init')
|
||
|
||
! Compute the domain of the ψ mapping
|
||
self%psi_domain = data%boundary
|
||
call rescale_boundary(self%psi_domain, self%psi_spline, O=self%axis, t0=1.5_wp_)
|
||
|
||
end subroutine numeric_init
|
||
|
||
|
||
pure subroutine rescale_boundary(cont, psi, O, t0)
|
||
! Given the plasma boundary contour `cont`, the position of the
|
||
! magnetic axis `O`, and a scaling factor `t0`; this subroutine
|
||
! rescales the contour by `t0` about `O` while ensuring the
|
||
! psi_spline stays monotonic within the new boundary.
|
||
|
||
! subroutine arguments
|
||
type(contour), intent(inout) :: cont ! (R,z) contour
|
||
type(spline_2d), intent(inout) :: psi ! ψ(R,z) spline
|
||
real(wp_), intent(in) :: O(2) ! center point
|
||
real(wp_), intent(in) :: t0 ! scaling factor
|
||
|
||
! subroutine variables
|
||
integer :: i
|
||
real(wp_) :: t
|
||
real(wp_), parameter :: dt = 0.05
|
||
real(wp_) :: P(2), N(2)
|
||
|
||
do i = 1, size(cont%R)
|
||
! For each point on the contour compute:
|
||
P = [cont%R(i), cont%z(i)] ! point on the contour
|
||
N = P - O ! direction of the line from O to P
|
||
|
||
! Find the max t: s(t) = ψ(O + tN) is monotonic in [1, t]
|
||
t = 1
|
||
do while (t < t0)
|
||
if (s(t + dt) < s(t)) exit
|
||
t = t + dt
|
||
end do
|
||
|
||
! The final point is thus O + tN
|
||
P = O + t * N
|
||
cont%R(i) = P(1)
|
||
cont%z(i) = P(2)
|
||
end do
|
||
|
||
contains
|
||
|
||
pure function s(t)
|
||
! Rescriction of ψ(R, z) on the line Q(t) = O + tN
|
||
|
||
real(wp_), intent(in) :: t
|
||
real(wp_) :: s, Q(2)
|
||
|
||
Q = O + t * N
|
||
s = psi%eval(Q(1), Q(2))
|
||
end function
|
||
|
||
end subroutine rescale_boundary
|
||
|
||
|
||
subroutine find_ox_point(self, R0, z0, R1, z1, psi1)
|
||
! Given the point (R₀,z₀) as an initial guess, finds
|
||
! the exact location (R₁,z₁) where ∇ψ(R₁,z₁) = 0.
|
||
! It also returns ψ₁=ψ(R₁,z₁).
|
||
!
|
||
! Note: this is used to find the magnetic X and O point,
|
||
! because both are stationary points for ψ(R,z).
|
||
use const_and_precisions, only : comp_eps
|
||
use minpack, only : hybrj1
|
||
use logger, only : log_error, log_debug
|
||
|
||
! subroutine arguments
|
||
class(numeric_equil) :: self
|
||
real(wp_), intent(in) :: R0, z0
|
||
real(wp_), intent(out) :: R1, z1, psi1
|
||
|
||
! local variables
|
||
integer :: info
|
||
real(wp_) :: sol(2), f(2), df(2,2), wa(15)
|
||
character(256) :: msg
|
||
|
||
sol = [R0, z0] ! first guess
|
||
|
||
call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, &
|
||
tol=sqrt(comp_eps), info=info, wa=wa, lwa=15)
|
||
if (info /= 1) then
|
||
write (msg, '("guess: ", g0.3, ", ", g0.3)') R0, z0
|
||
call log_debug(msg, mod='equilibrium', proc='find_ox_point')
|
||
write (msg, '("solution: ", g0.3, ", ", g0.3)') sol
|
||
call log_debug(msg, mod='equilibrium', proc='find_ox_point')
|
||
write (msg, '("hybrj1 failed with error ",g0)') info
|
||
call log_error(msg, mod='equilibrium', proc='find_ox_point')
|
||
end if
|
||
|
||
R1 = sol(1) ! solution
|
||
z1 = sol(2)
|
||
call self%pol_flux(R1, z1, psi1)
|
||
|
||
contains
|
||
|
||
pure subroutine equation(n, x, f, df, ldf, flag)
|
||
! The equation to solve: f(R,z) = ∇ψ(R,z) = 0
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: n, flag, ldf
|
||
real(wp_), intent(in) :: x(n)
|
||
real(wp_), intent(inout) :: f(n), df(ldf,n)
|
||
|
||
if (flag == 1) then
|
||
! return f(R,z) = ∇ψ(R,z)
|
||
call self%pol_flux(R=x(1), z=x(2), dpsidr=f(1), dpsidz=f(2))
|
||
else
|
||
! return ∇f(R,z) = ∇∇ψ(R,z)
|
||
call self%pol_flux(R=x(1), z=x(2), ddpsidrr=df(1,1), &
|
||
ddpsidzz=df(2,2), ddpsidrz=df(1,2))
|
||
df(2,1) = df(1,2)
|
||
end if
|
||
end subroutine equation
|
||
|
||
end subroutine find_ox_point
|
||
|
||
|
||
subroutine find_htg_point(self, R0, z0, R1, z1, psi0)
|
||
! Given the point (R₀,z₀) as an initial guess, finds
|
||
! the exact location (R₁,z₁) where:
|
||
! { ψ(R₁,z₁) = ψ₀
|
||
! { ∂ψ/∂R(R₁,z₁) = 0 .
|
||
!
|
||
! Note: this is used to find the horizontal tangent
|
||
! point of the contour ψ(R,z)=ψ₀.
|
||
use const_and_precisions, only : comp_eps
|
||
use minpack, only : hybrj1
|
||
use logger, only : log_error, log_debug
|
||
|
||
! subroutine arguments
|
||
class(numeric_equil) :: self
|
||
real(wp_), intent(in) :: R0, z0, psi0
|
||
real(wp_), intent(out) :: R1, z1
|
||
|
||
! local variables
|
||
integer :: info
|
||
real(wp_) :: sol(2), f(2), df(2,2), wa(15)
|
||
character(256) :: msg
|
||
|
||
sol = [R0, z0] ! first guess
|
||
|
||
call hybrj1(equation, n=2, x=sol, fvec=f, fjac=df, ldfjac=2, &
|
||
tol=sqrt(comp_eps), info=info, wa=wa, lwa=15)
|
||
if (info /= 1) then
|
||
write (msg, '("guess: ", g0.3, ", ", g0.3)') R0, z0
|
||
call log_debug(msg, mod='equilibrium', proc='find_ox_point')
|
||
write (msg, '("solution: ", g0.3, ", ", g0.3)') sol
|
||
call log_debug(msg, mod='equilibrium', proc='find_htg_point')
|
||
write (msg, '("hybrj1 failed with error ",g0)') info
|
||
call log_error(msg, mod='equilibrium', proc='find_htg_point')
|
||
end if
|
||
|
||
R1 = sol(1) ! solution
|
||
z1 = sol(2)
|
||
contains
|
||
|
||
pure subroutine equation(n, x, f, df, ldf, flag)
|
||
! The equation to solve: f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R] = 0
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: n, ldf, flag
|
||
real(wp_), intent(in) :: x(n)
|
||
real(wp_), intent(inout) :: f(n), df(ldf, n)
|
||
|
||
if (flag == 1) then
|
||
! return f(R,z) = [ψ(R,z)-ψ₀, ∂ψ/∂R]
|
||
call self%pol_flux(R=x(1), z=x(2), psi_n=f(1), dpsidr=f(2))
|
||
f(1) = f(1) - psi0
|
||
else
|
||
! return ∇f(R,z) = [[∂ψ/∂R, ∂ψ/∂z], [∂²ψ/∂R², ∂²ψ/∂R∂z]]
|
||
call self%pol_flux(R=x(1), z=x(2), dpsidr=df(1,1), dpsidz=df(1,2), &
|
||
ddpsidrr=df(2,1), ddpsidrz=df(2,2))
|
||
end if
|
||
end subroutine equation
|
||
|
||
end subroutine find_htg_point
|
||
|
||
|
||
!
|
||
! Vacuum
|
||
!
|
||
|
||
pure function vacuum_init() result(self)
|
||
! function arguments
|
||
type(vacuum) :: self
|
||
|
||
self%psi_a = 0
|
||
self%phi_a = 0
|
||
self%b_axis = 0
|
||
self%b_centre = 0
|
||
self%r_centre = 1
|
||
self%sgn_bphi = 0
|
||
self%axis = [0, 0]
|
||
self%r_range = [-comp_huge, comp_huge]
|
||
self%z_range = [-comp_huge, comp_huge]
|
||
self%z_boundary = [0, 0]
|
||
end function vacuum_init
|
||
|
||
|
||
pure subroutine vacuum_pol_flux(self, R, z, psi_n, dpsidr, dpsidz, &
|
||
ddpsidrr, ddpsidzz, ddpsidrz)
|
||
! function arguments
|
||
class(vacuum), intent(in) :: self
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_), intent(out), optional :: &
|
||
psi_n, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz
|
||
|
||
unused(self); unused(R); unused(z)
|
||
|
||
! ψ(R,z) is undefined everywhere, return -1
|
||
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
|
||
end subroutine vacuum_pol_flux
|
||
|
||
|
||
pure subroutine vacuum_pol_curr(self, psi_n, fpol, dfpol)
|
||
! subroutine arguments
|
||
class(vacuum), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_), intent(out) :: fpol ! poloidal current
|
||
real(wp_), intent(out), optional :: dfpol ! its derivative
|
||
|
||
unused(self); unused(psi_n)
|
||
|
||
! There is no current, F(ψ)=0
|
||
fpol = 0
|
||
if (present(dfpol)) dfpol = 0
|
||
end subroutine vacuum_pol_curr
|
||
|
||
|
||
pure function vacuum_safety(self, psi_n) result(q)
|
||
! function arguments
|
||
class(vacuum), intent(in) :: self
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_) :: q
|
||
|
||
unused(self); unused(psi_n)
|
||
|
||
! q(ψ) is undefined, return 0
|
||
q = 0
|
||
end function vacuum_safety
|
||
|
||
|
||
pure function vacuum_conv(self, rho_in) result(rho_out)
|
||
! function arguments
|
||
class(vacuum), intent(in) :: self
|
||
real(wp_), intent(in) :: rho_in
|
||
real(wp_) :: rho_out
|
||
|
||
unused(self)
|
||
|
||
! Neither ρ_p nor ρ_t are defined, do nothing
|
||
rho_out = rho_in
|
||
end function vacuum_conv
|
||
|
||
|
||
pure subroutine vacuum_flux_contour(self, psi0, R, z, &
|
||
R_hi, z_hi, R_lo, z_lo)
|
||
! subroutine arguments
|
||
class(vacuum), intent(in) :: self
|
||
real(wp_), intent(in) :: psi0
|
||
real(wp_), intent(out) :: R(:), z(:)
|
||
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
|
||
|
||
unused(self); unused(psi0);
|
||
unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo)
|
||
|
||
! flux surfaces are undefined
|
||
R = 0
|
||
z = 0
|
||
end subroutine vacuum_flux_contour
|
||
|
||
|
||
!
|
||
! Helpers
|
||
!
|
||
|
||
subroutine load_equil(params, equil, limiter, err)
|
||
! Loads a generic MHD equilibrium and limiter
|
||
! contour from file (params%filenm)
|
||
use gray_params, only : equilibrium_parameters, EQ_VACUUM
|
||
use gray_params, only : EQ_ANALYTICAL, EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL
|
||
use types, only : contour
|
||
use logger, only : log_warning, log_debug
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(inout) :: params
|
||
class(abstract_equil), allocatable, intent(out) :: equil
|
||
type(contour), intent(out) :: limiter
|
||
integer, intent(out) :: err
|
||
|
||
! local variables
|
||
type(numeric_equil) :: ne
|
||
type(analytic_equil) :: ae
|
||
type(eqdsk_data) :: data
|
||
|
||
select case (params%iequil)
|
||
case (EQ_VACUUM)
|
||
call log_debug('vacuum, doing nothing', mod='gray_equil', proc='load_equil')
|
||
allocate(equil, source=vacuum())
|
||
return
|
||
|
||
case (EQ_ANALYTICAL)
|
||
call log_debug('loading analytical file', mod='gray_equil', proc='load_equil')
|
||
call load_analytical(params, ae, limiter, err)
|
||
if (err /= 0) return
|
||
allocate(equil, source=ae)
|
||
|
||
case (EQ_EQDSK_FULL, EQ_EQDSK_PARTIAL)
|
||
call log_debug('loading G-EQDSK file', mod='gray_equil', proc='load_equil')
|
||
call load_eqdsk(params, data, err)
|
||
if (err /= 0) return
|
||
|
||
call change_cocos(data, params%icocos, 3)
|
||
call scale_eqdsk(params, data)
|
||
|
||
allocate(equil, mold=ne)
|
||
select type (equil)
|
||
type is (numeric_equil)
|
||
call equil%init(params, data, err)
|
||
if (err /= 0) return
|
||
end select
|
||
|
||
limiter = data%limiter
|
||
|
||
end select
|
||
|
||
call equil%init_flux_splines
|
||
|
||
end subroutine load_equil
|
||
|
||
|
||
subroutine load_analytical(params, equil, limiter, err)
|
||
! Loads the analytical equilibrium and limiter contour from file
|
||
use const_and_precisions, only : pi, one
|
||
use gray_params, only : equilibrium_parameters
|
||
use logger, only : log_error
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(in) :: params
|
||
type(analytic_equil), intent(out) :: equil
|
||
type(contour), intent(out) :: limiter
|
||
integer, intent(out) :: err
|
||
|
||
! local variables
|
||
integer :: i, u, nlim
|
||
real(wp_), allocatable :: R(:), z(:)
|
||
|
||
open(file=params%filenm, status='old', action='read', newunit=u, iostat=err)
|
||
if (err /= 0) then
|
||
call log_error('opening equilibrium ('//trim(params%filenm)//') failed!', &
|
||
mod='gray_equil', proc='load_analytical')
|
||
err = 1
|
||
return
|
||
end if
|
||
|
||
! The analytical file format is:
|
||
!
|
||
! 1 R0 z0 a
|
||
! 2 B0
|
||
! 3 q0 q1 alpha
|
||
! 4 nlim
|
||
! 5 R z
|
||
! ...
|
||
|
||
! Read model parameters
|
||
read (u, *) equil%R0, equil%z0, equil%a
|
||
read (u, *) equil%B0
|
||
read (u, *) equil%q0, equil%q1, equil%alpha
|
||
read (u, *) nlim
|
||
|
||
! Read limiter points
|
||
if (nlim > 0) then
|
||
allocate(R(nlim), z(nlim))
|
||
read (u, *) (R(i), z(i), i = 1, nlim)
|
||
limiter = contour(R, z)
|
||
end if
|
||
close(u)
|
||
|
||
! Notes on cocos=3
|
||
! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0.
|
||
! 2. ψ_a = ψ_edge - ψ_axis < 0.
|
||
! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0.
|
||
! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ).
|
||
|
||
! Apply forced signs
|
||
if (params%sgni /= 0) then
|
||
equil%q0 = sign(equil%q0, -params%sgni*one)
|
||
equil%q1 = sign(equil%q1, -params%sgni*one)
|
||
end if
|
||
if (params%sgnb /= 0) then
|
||
equil%B0 = sign(equil%B0, +params%sgnb*one)
|
||
end if
|
||
|
||
! Apply rescaling factors
|
||
equil%B0 = equil%B0 * params%factb
|
||
|
||
! Toroidal flux at r=a:
|
||
!
|
||
! Φ(a) = B₀πa² 2γ/(γ + 1)
|
||
!
|
||
! where γ=1/√(1-ε²),
|
||
! ε=a/R₀ is the tokamak aspect ratio
|
||
block
|
||
real(wp_) :: gamma
|
||
gamma = 1/sqrt(1 - (equil%a/equil%R0)**2)
|
||
equil%phi_a = equil%B0 * pi * equil%a**2 * 2*gamma/(gamma + 1)
|
||
end block
|
||
|
||
! In cocos=3 the safety factor is
|
||
!
|
||
! q(ψ) = 1/2π ∂Φ/∂ψ.
|
||
!
|
||
! Given the power law of the model
|
||
!
|
||
! q(ψ) = q₀ + (q₁-q₀) (ψ/ψa)^(α/2),
|
||
!
|
||
! we can find ψ_a = ψ(r=a) by integrating:
|
||
!
|
||
! ∫ q(ψ)dψ = 1/2π ∫ dΦ
|
||
! ∫₀^ψ_a q(ψ)dψ = 1/2π Φ_a
|
||
! ψa [q₀ + (q₁-q₀)/(α/2+1)] = Φa/2π
|
||
!
|
||
! ⇒ ψ_a = Φ_a 1/2π 1/(q₀ + Δq)
|
||
!
|
||
! where Δq = (q₁ - q₀)/(α/2 + 1)
|
||
block
|
||
real(wp_) :: dq
|
||
dq = (equil%q1 - equil%q0) / (equil%alpha/2 + 1)
|
||
equil%psi_a = 1/(2*pi) * equil%phi_a / (equil%q0 + dq)
|
||
end block
|
||
|
||
! Compute other constants (see abstract_equil)
|
||
equil%b_axis = equil%B0
|
||
equil%b_centre = equil%B0
|
||
equil%r_centre = equil%R0
|
||
equil%sgn_bphi = sign(one, equil%B0)
|
||
equil%axis = [equil%R0, equil%z0]
|
||
equil%r_range = equil%R0 + [-equil%a, equil%a]
|
||
equil%z_range = equil%z0 + [-equil%a, equil%a]
|
||
equil%z_boundary = equil%z0 + [-equil%a, equil%a]
|
||
end subroutine load_analytical
|
||
|
||
|
||
subroutine load_eqdsk(params, data, err)
|
||
! Loads the equilibrium `data` from a G-EQDSK file (params%filenm).
|
||
! For a description of the G-EQDSK format, see the GRAY user manual.
|
||
use const_and_precisions, only : one
|
||
use gray_params, only : equilibrium_parameters
|
||
use logger, only : log_error
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(in) :: params
|
||
type(eqdsk_data), intent(out) :: data
|
||
integer, intent(out) :: err
|
||
|
||
! local variables
|
||
integer :: u, i, j, nr, nz, nlim, nbnd
|
||
character(len=48) :: string
|
||
real(wp_) :: r_dim, z_dim, r_left, z_mid, psi_edge, psi_axis
|
||
real(wp_) :: skip_r ! dummy variables, used to discard data
|
||
integer :: skip_i !
|
||
|
||
! Open the G-EQDSK file
|
||
open(file=params%filenm, status='old', action='read', newunit=u, iostat=err)
|
||
if (err /= 0) then
|
||
call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', &
|
||
mod='gray_equil', proc='load_eqdsk')
|
||
err = 1
|
||
return
|
||
end if
|
||
|
||
! Get size of main arrays and allocate them
|
||
if (params%idesc) then
|
||
read (u,'(a48,3i4)') string, skip_i, nr, nz
|
||
else
|
||
read (u,*) nr, nz
|
||
end if
|
||
allocate(data%grid_r(nr), data%grid_z(nz), data%psi_map(nr, nz))
|
||
allocate(data%psi(nr), data%fpol(nr), data%q(nr))
|
||
|
||
! Store 0D data and main arrays
|
||
if (params%ifreefmt) then
|
||
read (u, *) r_dim, z_dim, data%r_ref, r_left, z_mid
|
||
read (u, *) data%axis, psi_axis, psi_edge, skip_r
|
||
read (u, *) skip_r, skip_r, skip_r, skip_r, skip_r
|
||
read (u, *) skip_r, skip_r, skip_r, skip_r, skip_r
|
||
read (u, *) (data%fpol(i), i=1,nr)
|
||
read (u, *) (skip_r,i=1, nr)
|
||
read (u, *) (skip_r,i=1, nr)
|
||
read (u, *) (skip_r,i=1, nr)
|
||
read (u, *) ((data%psi_map(i,j), i=1,nr), j=1,nz)
|
||
read (u, *) (data%q(i), i=1,nr)
|
||
else
|
||
read (u, '(5e16.9)') r_dim, z_dim, data%r_ref, r_left, z_mid
|
||
read (u, '(5e16.9)') data%axis, psi_axis, psi_edge, skip_r
|
||
read (u, '(5e16.9)') skip_r, skip_r, skip_r, skip_r, skip_r
|
||
read (u, '(5e16.9)') skip_r, skip_r, skip_r, skip_r, skip_r
|
||
read (u, '(5e16.9)') (data%fpol(i), i=1,nr)
|
||
read (u, '(5e16.9)') (skip_r, i=1,nr)
|
||
read (u, '(5e16.9)') (skip_r, i=1,nr)
|
||
read (u, '(5e16.9)') (skip_r, i=1,nr)
|
||
read (u, '(5e16.9)') ((data%psi_map(i,j), i=1,nr), j=1,nz)
|
||
read (u, '(5e16.9)') (data%q(i), i=1,nr)
|
||
end if
|
||
|
||
! Get size of boundary and limiter arrays and allocate them
|
||
read (u, *) nbnd, nlim
|
||
|
||
! Load plasma boundary data
|
||
if (nbnd > 0) then
|
||
block
|
||
real(wp_) :: R(nbnd), z(nbnd)
|
||
if (params%ifreefmt) then
|
||
read(u, *) (R(i), z(i), i=1,nbnd)
|
||
else
|
||
read(u, '(5e16.9)') (R(i), z(i), i=1,nbnd)
|
||
end if
|
||
data%boundary = contour(R, z)
|
||
end block
|
||
end if
|
||
|
||
! Load limiter data
|
||
if (nlim > 0) then
|
||
block
|
||
real(wp_) :: R(nlim), z(nlim)
|
||
if (params%ifreefmt) then
|
||
read(u, *) (R(i), z(i), i=1,nlim)
|
||
else
|
||
read(u, '(5e16.9)') (R(i), z(i), i=1,nlim)
|
||
end if
|
||
data%limiter = contour(R, z)
|
||
end block
|
||
end if
|
||
|
||
! End of G-EQDSK file
|
||
close(u)
|
||
|
||
! Build the grid arrays
|
||
block
|
||
real(wp_) :: dpsi, dr, dz
|
||
dr = r_dim/(nr - 1)
|
||
dz = z_dim/(nz - 1)
|
||
dpsi = one/(nr - 1)
|
||
do i = 1, nr
|
||
data%psi(i) = (i-1)*dpsi
|
||
data%grid_r(i) = r_left + (i-1)*dr
|
||
end do
|
||
do i = 1, nz
|
||
data%grid_z(i) = z_mid - z_dim/2 + (i-1)*dz
|
||
end do
|
||
end block
|
||
|
||
! Normalize the poloidal flux
|
||
data%psi_a = psi_edge - psi_axis
|
||
if (.not. params%ipsinorm) data%psi_map = (data%psi_map - psi_axis)/data%psi_a
|
||
|
||
end subroutine load_eqdsk
|
||
|
||
|
||
subroutine scale_eqdsk(params, data)
|
||
! Rescale the magnetic field (B) and the plasma current (I_p)
|
||
! and/or force their signs.
|
||
!
|
||
! Notes:
|
||
! 1. signi and signb are ignored on input if equal to 0.
|
||
! They are used to assign the direction of B_φ and I_p BEFORE scaling.
|
||
! 2. cocos=3 assumed: positive toroidal direction is CCW from above
|
||
! 3. B_φ and I_p scaled by the same factor factb to keep q unchanged
|
||
! 4. factb<0 reverses the directions of Bphi and Ipla
|
||
use const_and_precisions, only : one
|
||
use gray_params, only : equilibrium_parameters
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(inout) :: params
|
||
type(eqdsk_data), intent(inout) :: data
|
||
|
||
! Notes on cocos=3
|
||
! 1. With both I_p,B_φ > 0 we have ∂ψ/∂r<0 and ∂Φ/∂r>0.
|
||
! 2. ψ_a = ψ_edge - ψ_axis < 0.
|
||
! 3. q = 1/2π ∂Φ/∂ψ ~ ∂Φ/∂r⋅∂r/∂ψ < 0.
|
||
! 4. In general, sgn(q) = -sgn(I_p)⋅sgn(B_φ).
|
||
|
||
! Apply signs
|
||
if (params%sgni /= 0) data%psi_a = sign(data%psi_a, -params%sgni*one)
|
||
if (params%sgnb /= 0) data%fpol = sign(data%fpol, +params%sgnb*one)
|
||
|
||
! Rescale
|
||
data%psi_a = data%psi_a * params%factb
|
||
data%fpol = data%fpol * params%factb
|
||
|
||
! Compute the signs to be shown in the outputs header
|
||
! Note: this is needed if sgni, sgnb = 0 in gray.ini
|
||
params%sgni = int(sign(one, -data%psi_a))
|
||
params%sgnb = int(sign(one, +data%fpol(size(data%fpol))))
|
||
end subroutine scale_eqdsk
|
||
|
||
|
||
subroutine change_cocos(data, cocosin, cocosout, err)
|
||
! Convert the MHD equilibrium data from one coordinate convention
|
||
! (COCOS) to another. These are specified by `cocosin` and
|
||
! `cocosout`, respectively.
|
||
!
|
||
! For more information, see: https://doi.org/10.1016/j.cpc.2012.09.010
|
||
use const_and_precisions, only : zero, one, pi
|
||
use logger, only : log_warning
|
||
|
||
! subroutine arguments
|
||
type(eqdsk_data), intent(inout) :: data
|
||
integer, intent(in) :: cocosin, cocosout
|
||
integer, intent(out), optional :: err
|
||
|
||
! local variables
|
||
real(wp_) :: isign, bsign
|
||
integer :: exp2pi, exp2piout
|
||
logical :: phiccw, psiincr, qpos, phiccwout, psiincrout, qposout
|
||
|
||
call decode_cocos(cocosin, exp2pi, phiccw, psiincr, qpos)
|
||
call decode_cocos(cocosout, exp2piout, phiccwout, psiincrout, qposout)
|
||
|
||
! Check sign consistency
|
||
isign = sign(one, data%psi_a)
|
||
if (.not. psiincr) isign = -isign
|
||
bsign = sign(one, data%fpol(size(data%fpol)))
|
||
if (qpos .neqv. isign * bsign * data%q(size(data%q)) > zero) then
|
||
! Warning: sign inconsistency found among q, I_p and B_ref
|
||
data%q = -data%q
|
||
call log_warning('data not consistent with cocosin', &
|
||
mod='gray_equil', proc='change_cocos')
|
||
if (present(err)) err = 1
|
||
else
|
||
if (present(err)) err = 0
|
||
end if
|
||
|
||
! Convert cocosin to cocosout
|
||
|
||
! Opposite direction of toroidal angle phi in cocosin and cocosout
|
||
if (phiccw .neqv. phiccwout) data%fpol = -data%fpol
|
||
|
||
! q has opposite sign for given sign of B_phi⋅I_p
|
||
if (qpos .neqv. qposout) data%q = -data%q
|
||
|
||
! psi and Ip signs don't change accordingly
|
||
if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) &
|
||
data%psi_a = -data%psi_a
|
||
|
||
! Convert Wb to Wb/rad or viceversa
|
||
data%psi_a = data%psi_a * (2*pi)**(exp2piout - exp2pi)
|
||
end subroutine change_cocos
|
||
|
||
|
||
subroutine decode_cocos(cocos, exp2pi, phiccw, psiincr, qpos)
|
||
! Extracts the sign and units conventions from a COCOS index
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: cocos
|
||
integer, intent(out) :: exp2pi
|
||
logical, intent(out) :: phiccw, psiincr, qpos
|
||
|
||
! local variables
|
||
integer :: cmod10, cmod4
|
||
|
||
cmod10 = mod(cocos, 10)
|
||
cmod4 = mod(cmod10, 4)
|
||
|
||
! cocos>10 ψ in Wb, cocos<10 ψ in Wb/rad
|
||
exp2pi = (cocos - cmod10)/10
|
||
|
||
! cocos mod 10 = 1,3,5,7: toroidal angle φ increasing CCW
|
||
phiccw = (mod(cmod10, 2)== 1)
|
||
|
||
! cocos mod 10 = 1,2,5,6: ψ increasing with positive Ip
|
||
psiincr = (cmod4==1 .or. cmod4==2)
|
||
|
||
! cocos mod 10 = 1,2,7,8: q positive for positive Bφ*Ip
|
||
qpos = (cmod10<3 .or. cmod10>6)
|
||
end subroutine decode_cocos
|
||
|
||
|
||
subroutine init_flux_splines(self)
|
||
! Computes the splines of the flux surface averages
|
||
|
||
use const_and_precisions, only : zero, one, comp_huge, pi, mu0inv
|
||
|
||
! subroutine arguments
|
||
class(abstract_equil), intent(inout) :: self
|
||
|
||
! local constants
|
||
integer, parameter :: nsurf=101, npoints=101, nlam=101
|
||
|
||
! local variables
|
||
integer :: i, j, l
|
||
|
||
! Arrays of quantities evaluated on a uniform ρ_p grid
|
||
real(wp_) :: rho_p(nsurf) ! ρ_p ∈ [0, 1)
|
||
real(wp_) :: psi_n(nsurf) ! ψ_n = ρ_p²
|
||
real(wp_) :: B_avg(nsurf) ! ⟨B⟩ = 1/N ∮ B dl/B_p, average magnetic field
|
||
real(wp_) :: J_phi_avg(nsurf) ! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p, average toroidal current
|
||
real(wp_) :: I_p(nsurf) ! I_p = 1/μ₀ ∫ B_p⋅dS_φ, plasma current
|
||
real(wp_) :: q(nsurf) ! q = 1/2π ∮ B_φ/B_p dl/R, safety factor
|
||
real(wp_) :: R_J(nsurf) ! R_J = ⟨B⟩ / (F(ψ) ⋅ ⟨1/R²⟩), effective J_cd radius
|
||
real(wp_) :: area(nsurf) ! poloidal area enclosed by the flux surface
|
||
real(wp_) :: volume(nsurf) ! volume V enclosed by the flux surface
|
||
real(wp_) :: dAdrho_t(nsurf) ! dA/dρ_t, where A is the poloidal area
|
||
real(wp_) :: dVdrho_t(nsurf) ! dV/dρ_t, where V is the enclosed volume
|
||
real(wp_) :: B_min(nsurf) ! min |B| on the flux surface
|
||
real(wp_) :: B_max(nsurf) ! max |B| on the flux surface
|
||
real(wp_) :: f_c(nsurf) ! fraction of circulating particles
|
||
|
||
! Ratios of different J_cd definitions w.r.t. ⟨J_φ⟩
|
||
real(wp_) :: ratio_astra_tor(nsurf) ! J_cd_astra = ⟨J_cd⋅B⟩/B₀
|
||
real(wp_) :: ratio_jintrac_tor(nsurf) ! J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩
|
||
real(wp_) :: ratio_paral_tor(nsurf) ! Jcd_∥ = ⟨J_cd⋅B/|B|⟩
|
||
|
||
! Arrays of quantities evaluated on the poloidal flux contour
|
||
real(wp_) :: R(npoints), z(npoints) ! R, z coordinates
|
||
real(wp_) :: B_tots(npoints) ! magnetic field
|
||
real(wp_) :: B_pols(npoints) ! poloidal field
|
||
real(wp_) :: dls(npoints) ! line element
|
||
|
||
! Intermediate results
|
||
real(wp_) :: norm ! N = ∮ dl/B_p, normalisation of the ⟨⋅⟩ integral
|
||
real(wp_) :: B2_avg ! ⟨B²⟩ = 1/N ∮ B² dl/B_p, average squared field
|
||
real(wp_) :: dAdpsi ! dA/dψ, where A is the poloidal area
|
||
real(wp_) :: dVdpsi ! dV/dψ, where V is the enclosed volume
|
||
real(wp_) :: R2_inv_avg ! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p
|
||
real(wp_) :: R_inv_avg ! ⟨R⁻¹⟩ = 1/N ∮ R⁻¹ dl/B_p
|
||
real(wp_) :: fpol ! F(ψ), poloidal current function
|
||
|
||
! Neoclassical CD efficiency variables
|
||
real(wp_) :: dlam, lambda(nlam) ! uniform grid λ ∈ [0,1]
|
||
real(wp_) :: H(nlam, nsurf) ! H(ρ_p, λ) sampled on a uniform ρ_p×λ grid
|
||
|
||
! Miscellanea
|
||
real(wp_) :: B_R, B_z, B_phi, R_hi, R_lo, z_hi, z_lo
|
||
|
||
! On the magnetic axis, ρ_p = 0
|
||
psi_n(1) = 0
|
||
rho_p(1) = 0
|
||
area(1) = 0
|
||
volume(1) = 0
|
||
I_p(1) = 0
|
||
J_phi_avg(1) = self%tor_curr(self%axis(1), self%axis(2))
|
||
B_max(1) = abs(self%b_axis)
|
||
B_min(1) = abs(self%b_axis)
|
||
B_avg(1) = abs(self%b_axis)
|
||
R_J(1) = self%axis(1)
|
||
f_c(1) = 1
|
||
dAdrho_t(1) = 0
|
||
dVdrho_t(1) = 0
|
||
|
||
ratio_paral_tor(1) = 1
|
||
ratio_astra_tor(1) = abs(self%b_axis / self%b_centre)
|
||
ratio_jintrac_tor(1) = 1
|
||
|
||
! The safety factor at ρ=0 is given by
|
||
! q₀ = B₀ / √(∂²ψ/∂R² ⋅ ∂²ψ/∂z²)
|
||
! Source: Spheromaks, Bellan
|
||
block
|
||
real(wp_) :: a, b
|
||
call self%pol_flux(self%axis(1), self%axis(2), ddpsidrr=a, ddpsidzz=b)
|
||
q(1) = self%b_axis / sqrt(a*b)/self%psi_a
|
||
end block
|
||
|
||
! Compute the uniform λ grid and H(0, λ)
|
||
H = 0
|
||
dlam = one/(nlam - 1)
|
||
do l = 1, nlam
|
||
lambda(l) = (l-1) * dlam ! λ
|
||
H(l,1) = sqrt(1 - lambda(l)) ! H(0, λ) = √(1-λ)
|
||
end do
|
||
|
||
! Guess for first contour reconstruction
|
||
R_hi = self%axis(1)
|
||
R_lo = self%axis(1)
|
||
z_hi = self%axis(2) + (self%z_boundary(2) - self%axis(2))/10
|
||
z_lo = self%axis(2) - (self%axis(2) - self%z_boundary(1))/10
|
||
|
||
! For each surface with ρ_p > 0
|
||
surfaces_loop: do i = 2, nsurf
|
||
! Note: keep ρ_p < 1 to avoid the X point
|
||
rho_p(i) = merge((i - 1) * one/(nsurf - 1), 0.9999_wp_, i < nsurf)
|
||
psi_n(i) = rho_p(i)**2
|
||
|
||
call self%flux_contour(psi_n(i), R, z, R_hi, z_hi, R_lo, z_lo)
|
||
|
||
block
|
||
! Variables for the previous and current contour point
|
||
real(wp_) :: B_tot0, B_tot1, B_pol0, B_pol1, J_phi0, J_phi1
|
||
real(wp_) :: dl, dir ! line element, sgn(dl⋅B_p)
|
||
|
||
! Initialise all integrals
|
||
area(i) = 0
|
||
volume(i) = 0
|
||
norm = 0
|
||
dAdpsi = 0
|
||
I_p(i) = 0
|
||
B_avg(i) = 0
|
||
B2_avg = 0
|
||
R2_inv_avg = 0
|
||
J_phi_avg(i) = 0
|
||
B_max(i) = -comp_huge
|
||
B_min(i) = +comp_huge
|
||
|
||
! Evaluate integrands at the first "previous" point
|
||
J_phi0 = self%tor_curr(R(1), z(1))
|
||
call self%b_field(R(1), z(1), B_R=B_R, B_z=B_z, &
|
||
B_phi=B_phi, B_tot=B_tot0)
|
||
fpol = B_phi*R(1)
|
||
B_pol0 = sqrt(B_R**2 + B_z**2)
|
||
B_tots(1) = B_tot0
|
||
B_pols(1) = B_pol0
|
||
|
||
contour_loop: do j = 1, npoints-1
|
||
block
|
||
! Compute the area by decomposing the contour into triangles
|
||
! formed by the magnetic axis and two adjacent points, whose
|
||
! area is given by the cross-product of the vertices, ½|v×w|.
|
||
real(wp_) :: O(2), P(2), Q(2), tri_area
|
||
O = self%axis
|
||
P = [R(j), z(j)]
|
||
Q = [R(j+1), z(j+1)]
|
||
tri_area = abs(cross(P-O, Q-O))/2
|
||
area(i) = area(i) + tri_area
|
||
|
||
! Compute the volume using the centroid theorem (V = 2πR⋅A) with
|
||
! the centroid R of the contour as the mean of the triangles
|
||
! centroids weighted by their area. Note: the total area cancels.
|
||
volume(i) = volume(i) + 2*pi * (O(1) + P(1)+ Q(1))/3 * tri_area
|
||
|
||
! The line differential is difference of the two vertices
|
||
dl = norm2(P - Q)
|
||
|
||
! We always assume dl ∥ B_p (numerically they're not),
|
||
! but compute the sign of dl⋅B_p for the case I_p < 0.
|
||
dir = sign(one, dot_product([B_R, B_z], P - Q))
|
||
end block
|
||
|
||
! Evaluate integrands at the current point for line integrals ∮ dl
|
||
call self%b_field(R(j+1), z(j+1), B_R=B_R, B_z=B_z, B_tot=B_tot1)
|
||
J_phi1 = self%tor_curr(R(j+1), z(j+1))
|
||
B_pol1 = sqrt(B_R**2 + B_z**2)
|
||
dls(j) = dl
|
||
B_tots(j+1) = B_tot1
|
||
B_pols(j+1) = B_pol1
|
||
|
||
! Do one step of the trapezoid method
|
||
! N = ∮ dl/B_p
|
||
! dA/dψ = ∮ (1/R) dl/B_p
|
||
! ⟨B⟩ = 1/N ∮ B dl/B_p
|
||
! ⟨B²⟩ = 1/N ∮ B² dl/B_p
|
||
! ⟨R⁻²⟩ = 1/N ∮ R⁻² dl/B_p
|
||
! ⟨J_φ⟩ = 1/N ∮ J_φ dl/B_p
|
||
! I_p = 1/μ₀ ∮ B_p⋅dl
|
||
norm = norm + dl*(1/B_pol1 + 1/B_pol0)/2
|
||
dAdpsi = dAdpsi + dl*(1/(B_pol1*R(j+1)) + 1/(B_pol0*R(j)))/2
|
||
B_avg(i) = B_avg(i) + dl*(B_tot1/B_pol1 + B_tot0/B_pol0)/2
|
||
B2_avg = B2_avg + dl*(B_tot1**2/B_pol1 + B_tot0**2/B_pol0)/2
|
||
R2_inv_avg = R2_inv_avg + dl*(1/(B_pol1*R(j+1)**2) + 1/(B_pol0*R(j)**2))/2
|
||
J_phi_avg(i) = J_phi_avg(i) + dl*(J_phi0/(B_pol0*R(j)) + J_phi1/(B_pol1*R(j+1)))/2
|
||
I_p(i) = I_p(i) + dl*(B_pol1 + B_pol0)/2 * dir * mu0inv
|
||
|
||
! Update the previous values
|
||
J_phi0 = J_phi1
|
||
B_pol0 = B_pol1
|
||
B_tot0 = B_tot1
|
||
|
||
! Compute max,min magnetic field
|
||
if(B_tot1 <= B_min(i)) B_min(i) = B_tot1
|
||
if(B_tot1 >= B_max(i)) B_max(i) = B_tot1
|
||
end do contour_loop
|
||
end block
|
||
|
||
! Normalise integrals to obtain average value
|
||
B_avg(i) = B_avg(i) / norm
|
||
B2_avg = B2_avg / norm
|
||
R2_inv_avg = R2_inv_avg / norm
|
||
R_inv_avg = dAdpsi / norm
|
||
J_phi_avg(i) = J_phi_avg(i)/dAdpsi
|
||
dVdpsi = 2*pi * norm
|
||
|
||
! J_cd_astra/⟨J_φ⟩ where J_cd_astra = ⟨J_cd⋅B⟩/B₀
|
||
! J_cd_jintrac/⟨J_φ⟩ where J_cd_jintrac = ⟨J_cd⋅B⟩/⟨B⟩
|
||
! J_cd_∥/⟨J_φ⟩ where Jcd_∥ = ⟨J_cd⋅B/|B|⟩
|
||
ratio_astra_tor(i) = abs(B2_avg *R_inv_avg/(fpol*R2_inv_avg*self%b_centre))
|
||
ratio_jintrac_tor(i) = abs(B2_avg *R_inv_avg/(fpol*R2_inv_avg*B_avg(i)))
|
||
ratio_paral_tor(i) = abs(B_avg(i)*R_inv_avg/(fpol*R2_inv_avg))
|
||
|
||
R_J(i) = B_avg(i)/abs(fpol*R2_inv_avg)
|
||
|
||
! By definition q(γ) = 1/2π ∮_γ dφ, where γ is a field line:
|
||
!
|
||
! { γ̅(0) = γ̅₀
|
||
! { dγ̅/dt = B̅(γ(t))
|
||
!
|
||
! Using the poloidal angle θ, we rewrite it as:
|
||
!
|
||
! q(γ) = 1/2π ∮_γ dφ/dθ dθ = 1/2π ∮_γ dφ/dt / dθ/dt dθ
|
||
!
|
||
! By definition of directional derivative, d/dt = B̅⋅∇, so
|
||
!
|
||
! q(γ) = 1/2π ∮_γ B̅⋅∇φ / B̅⋅∇θ dθ = 1/2π ∮_γ B_φ/B_p ρ/R dθ
|
||
!
|
||
! Finally, since ρdθ=dl in the poloidal plane and B_φ=F/R:
|
||
!
|
||
! q(ρ) = F/2π ∮ 1/R² dl/B_p = F/2π N ⟨R⁻²⟩
|
||
!
|
||
q(i) = fpol/(2*pi) * norm * R2_inv_avg
|
||
|
||
! Using Φ=Φ_a⋅ρ_t² and q=1/2π⋅∂Φ/∂ψ (see gray_equil.f90), we have
|
||
!
|
||
! ∂ψ/∂ρ_t = ∂ψ/∂Φ ∂Φ/∂ρ_t = 1/(2πq) Φ_a 2ρ_t
|
||
! = Φ_a ρ_t / (πq)
|
||
block
|
||
real(wp_) :: rho_t, dpsidrho_t
|
||
|
||
rho_t = self%pol2tor(rho_p(i))
|
||
dpsidrho_t = self%phi_a * rho_t / (self%safety(psi_n(i)) * pi)
|
||
dAdrho_t(i) = dAdpsi * dpsidrho_t
|
||
dVdrho_t(i) = dVdpsi * dpsidrho_t
|
||
end block
|
||
|
||
! Compute the fraction of circulating particles:
|
||
!
|
||
! f_c(ρ_p) = ¾ <B²/B_max²> ∫₀¹ dλ λ/⟨√[1-λ'⋅B(ρ_p)/B_max]⟩
|
||
!
|
||
! and the function:
|
||
!
|
||
! H(ρ_p, λ) = ½ B_max/B_min/f_c ∫_λ^1 dλ'/⟨√[1-λ'⋅B(ρ_p)/B_max]⟩
|
||
!
|
||
f_c(i) = 0
|
||
|
||
block
|
||
real(wp_) :: srl, rl0, rl1, dHdlam1, dHdlam0
|
||
|
||
do l = nlam, 1, -1
|
||
! Compute srl = ⟨√[1-λ B(ρ_p)/B_max]⟩ with the trapezoid method
|
||
srl = 0
|
||
rl0 = sqrt(max(1 - lambda(l)*B_tots(1)/B_max(i), zero))
|
||
do j = 1, npoints-1
|
||
rl1 = sqrt(max(1 - lambda(l)*B_tots(j+1)/B_max(i), zero))
|
||
srl = srl + dls(j) * (rl1/B_pols(j+1) + rl0/B_pols(j))/2
|
||
rl0 = rl1
|
||
end do
|
||
srl = srl / norm
|
||
|
||
! Do one step of the Riemann sum of f_c
|
||
f_c(i) = f_c(i) + 3/4.0_wp_ * B2_avg/B_max(i)**2 &
|
||
* dlam * lambda(l)/srl * merge(1.0_wp_, 0.5_wp_, l > 1 .and. l < nlam)
|
||
|
||
! Do one step of the trapezoid method for ∫_λ^1 dλ'/srl
|
||
dHdlam1 = 1 / srl
|
||
if (l /= nlam) H(l, i) = H(l+1, i) + dlam*(dHdlam1 + dHdlam0)/2
|
||
dHdlam0 = dHdlam1
|
||
end do
|
||
end block
|
||
|
||
! Add leading factors
|
||
H(:, i) = H(:, i)/2 * B_min(i)/B_max(i)/f_c(i)
|
||
|
||
end do surfaces_loop
|
||
|
||
! Fake ρ_p = 1 on the last surface (actually slightly < 1)
|
||
rho_p(nsurf) = 1
|
||
psi_n(nsurf) = 1
|
||
|
||
! Compute splines
|
||
call self%spline_area%init(rho_p, area, nsurf)
|
||
call self%spline_volume%init(rho_p, volume, nsurf)
|
||
|
||
call self%spline_dAdrho_t%init(rho_p, dAdrho_t, nsurf)
|
||
call self%spline_dVdrho_t%init(rho_p, dVdrho_t, nsurf)
|
||
|
||
call self%spline_R_J%init(rho_p, R_J, nsurf)
|
||
call self%spline_B_avg%init(rho_p, B_avg, nsurf)
|
||
call self%spline_B_max%init(rho_p, B_max, nsurf)
|
||
call self%spline_B_min%init(rho_p, B_min, nsurf)
|
||
|
||
call self%spline_f_c%init(rho_p, f_c, nsurf)
|
||
call self%spline_q%init(rho_p, q, nsurf)
|
||
call self%spline_I_p%init(rho_p, I_p, nsurf)
|
||
|
||
call self%spline_J_phi_avg%init(rho_p, J_phi_avg, nsurf)
|
||
call self%spline_astra_tor%init(rho_p, ratio_astra_tor, nsurf)
|
||
call self%spline_jintrac_tor%init(rho_p, ratio_jintrac_tor, nsurf)
|
||
call self%spline_paral_tor%init(rho_p, ratio_paral_tor, nsurf)
|
||
|
||
call self%spline_cd_eff%init(rho_p, lambda, reshape(H, [nsurf*nlam]), &
|
||
nsurf, nlam, range=[zero, one, zero, one])
|
||
|
||
contains
|
||
|
||
pure function cross(v, w)
|
||
! z-component of the cross product of 2D vectors
|
||
real(wp_), intent(in) :: v(2), w(2)
|
||
real(wp_) :: cross
|
||
cross = v(1)*w(2) - v(2)*w(1)
|
||
end function cross
|
||
|
||
end subroutine init_flux_splines
|
||
|
||
end module
|