26976a31dd
This change limits the evaluation of the poloidal flux spline ψ(R,z) to a particular domain where ψ(R,z) is strictly monotonic. This choice ensures that ψ_n < 1 only inside the plasma boundary; so plasma parameters like q(ψ), Te(ψ), etc. can be mapped to the physical space without ambiguities. Before this change, anywhere ψ_n happened to decrease below 1 (either as a result of a physical current or simply from the spline extrapolation), Gray would effectively create a spurious plasma region. This behavior is seriously problematic because it completely invalidates the simulation: it can alter the ray directions, the power disribution and total power if the beam enters one such region. The domain is computed by applying a transformation to the contour of the plasma boundary: for each contour point we cast a ray from the magnetic axis to that point and extend the ray until the restriction of ψ(R,z) on the ray starts decreasing or reach a maximum scaling factor. The endpoint of the ray is then taken as the new point. If ψ(R,z) is globally monotonic, the transformation is a homotety wrt the magnetic axis, so the domain will be an enlarged boundary; otherwise the shape will be more irregular (an intersection of the enlarged boundary and several level curves of ψ). Finally, each `pol_flux(r, z)` call is now guarded behind a check `inside(psi_domain, r, z)`. For points outside the domain the subroutine returns, as usual, -1 for ψ and 0 for derivatives.
1516 lines
48 KiB
Fortran
1516 lines
48 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.)
|
||
!
|
||
! Two kinds of data are supported: analytical (suffix `_an` in the
|
||
! subroutine names) or numerical (suffix `_spline`). For the latter, the
|
||
! the data is interpolated using splines.
|
||
module equilibrium
|
||
use const_and_precisions, only : wp_
|
||
use splines, only : spline_simple, spline_1d, spline_2d, linear_1d
|
||
|
||
implicit none
|
||
|
||
! Parameters of the analytical equilibrium model
|
||
type analytic_model
|
||
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)
|
||
end type
|
||
|
||
! A 2D contour in the (R,z) plane
|
||
type contour
|
||
real(wp_), allocatable :: R(:)
|
||
real(wp_), allocatable :: z(:)
|
||
end type
|
||
|
||
! Order of the splines
|
||
integer, parameter :: kspl=3, ksplp=kspl + 1
|
||
|
||
! Global variable storing the state of the module
|
||
|
||
! Splines
|
||
type(spline_1d), save :: fpol_spline ! Poloidal current function F(ψ)
|
||
type(spline_2d), save :: psi_spline ! Normalised poloidal flux ψ(R,z)
|
||
type(spline_simple), save :: q_spline ! Safey factor q(ψ)
|
||
type(linear_1d), save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p)
|
||
|
||
! Analytic model
|
||
type(analytic_model), save :: model
|
||
|
||
! Contour of the region where the ψ(R,z) spline is monotonically
|
||
! increasing. The mapping from ψ to other plasma parameters is
|
||
! physically meaningful only inside this domain.
|
||
type(contour), save :: psi_domain
|
||
|
||
! More parameters
|
||
real(wp_), save :: fpolas ! Poloidal current at the edge (r=a), F_a
|
||
real(wp_), save :: psia ! Poloidal flux at the edge (r=a), ψ_a
|
||
real(wp_), save :: phitedge ! Toroidal flux at the edge
|
||
real(wp_), save :: btaxis ! B₀: B_φ = B₀R₀/R
|
||
real(wp_), save :: rmaxis, zmaxis ! R,z of the magnetic axis (O point)
|
||
real(wp_), save :: sgnbphi ! Sign of B_φ (>0 means CCW)
|
||
real(wp_), save :: btrcen ! used only for jcd_astra def.
|
||
real(wp_), save :: rcen ! Major radius R₀, computed as fpolas/btrcen
|
||
real(wp_), save :: rmnm, rmxm ! R interval of the equilibrium grid
|
||
real(wp_), save :: zmnm, zmxm ! z interval of the equilibrium grid
|
||
real(wp_), save :: zbinf, zbsup ! z interval of the plasma boundary
|
||
|
||
private
|
||
public read_eqdsk, read_equil_an ! Reading data files
|
||
public scale_equil, change_cocos ! Transforming data
|
||
public fq, bfield, tor_curr, tor_curr_psi ! Accessing local quantities
|
||
public pol_flux, pol_curr !
|
||
public frhotor, frhopol ! Converting between poloidal/toroidal flux
|
||
public set_equil_spline, set_equil_an ! Initialising internal state
|
||
public unset_equil_spline ! Deinitialising internal state
|
||
|
||
! Members exposed for magsurf_data
|
||
public kspl, psi_spline, points_tgo, model
|
||
|
||
! Members exposed to gray_core and more
|
||
public psia, phitedge
|
||
public btaxis, rmaxis, zmaxis, sgnbphi
|
||
public btrcen, rcen
|
||
public rmnm, rmxm, zmnm, zmxm
|
||
public zbinf, zbsup
|
||
|
||
contains
|
||
|
||
subroutine read_eqdsk(params, data, err, unit)
|
||
! Reads the MHD equilibrium `data` from a G-EQDSK file (params%filenm).
|
||
! If given, the file is opened in the `unit` number.
|
||
! For a description of the G-EQDSK, see the GRAY user manual.
|
||
use const_and_precisions, only : one
|
||
use gray_params, only : equilibrium_parameters, equilibrium_data
|
||
use utils, only : get_free_unit
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(in) :: params
|
||
type(equilibrium_data), intent(out) :: data
|
||
integer, intent(out) :: err
|
||
integer, optional, intent(in) :: unit
|
||
|
||
! local variables
|
||
integer :: u, idum, i, j, nr, nz, nbnd, nlim
|
||
character(len=48) :: string
|
||
real(wp_) :: dr, dz, dps, rleft, zmid, zleft, psiedge, psiaxis
|
||
real(wp_) :: xdum ! dummy variable, used to discard data
|
||
|
||
u = get_free_unit(unit)
|
||
|
||
! Open the G-EQDSK file
|
||
open(file=params%filenm, status='old', action='read', unit=u, iostat=err)
|
||
if (err /= 0) then
|
||
call log_error('opening eqdsk file ('//trim(params%filenm)//') failed!', &
|
||
mod='equilibrium', proc='read_eqdsk')
|
||
err = 1
|
||
return
|
||
end if
|
||
|
||
! get size of main arrays and allocate them
|
||
if (params%idesc == 1) then
|
||
read (u,'(a48,3i4)') string,idum,nr,nz
|
||
else
|
||
read (u,*) nr, nz
|
||
end if
|
||
if (allocated(data%rv)) deallocate(data%rv)
|
||
if (allocated(data%zv)) deallocate(data%zv)
|
||
if (allocated(data%psin)) deallocate(data%psin)
|
||
if (allocated(data%psinr)) deallocate(data%psinr)
|
||
if (allocated(data%fpol)) deallocate(data%fpol)
|
||
if (allocated(data%qpsi)) deallocate(data%qpsi)
|
||
allocate(data%rv(nr), data%zv(nz), &
|
||
data%psin(nr, nz), &
|
||
data%psinr(nr), &
|
||
data%fpol(nr), &
|
||
data%qpsi(nr))
|
||
|
||
! Store 0D data and main arrays
|
||
if (params%ifreefmt==1) then
|
||
read (u, *) dr, dz, data%rvac, rleft, zmid
|
||
read (u, *) data%rax, data%zax, psiaxis, psiedge, xdum
|
||
read (u, *) xdum, xdum, xdum, xdum, xdum
|
||
read (u, *) xdum, xdum, xdum, xdum, xdum
|
||
read (u, *) (data%fpol(i), i=1,nr)
|
||
read (u, *) (xdum,i=1, nr)
|
||
read (u, *) (xdum,i=1, nr)
|
||
read (u, *) (xdum,i=1, nr)
|
||
read (u, *) ((data%psin(i,j), i=1,nr), j=1,nz)
|
||
read (u, *) (data%qpsi(i), i=1,nr)
|
||
else
|
||
read (u, '(5e16.9)') dr,dz,data%rvac,rleft,zmid
|
||
read (u, '(5e16.9)') data%rax,data%zax,psiaxis,psiedge,xdum
|
||
read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
|
||
read (u, '(5e16.9)') xdum,xdum,xdum,xdum,xdum
|
||
read (u, '(5e16.9)') (data%fpol(i),i=1,nr)
|
||
read (u, '(5e16.9)') (xdum,i=1,nr)
|
||
read (u, '(5e16.9)') (xdum,i=1,nr)
|
||
read (u, '(5e16.9)') (xdum,i=1,nr)
|
||
read (u, '(5e16.9)') ((data%psin(i,j),i=1,nr),j=1,nz)
|
||
read (u, '(5e16.9)') (data%qpsi(i),i=1,nr)
|
||
end if
|
||
|
||
! Get size of boundary and limiter arrays and allocate them
|
||
read (u,*) nbnd, nlim
|
||
if (allocated(data%rbnd)) deallocate(data%rbnd)
|
||
if (allocated(data%zbnd)) deallocate(data%zbnd)
|
||
if (allocated(data%rlim)) deallocate(data%rlim)
|
||
if (allocated(data%zlim)) deallocate(data%zlim)
|
||
|
||
! Load plasma boundary data
|
||
if(nbnd > 0) then
|
||
allocate(data%rbnd(nbnd), data%zbnd(nbnd))
|
||
if (params%ifreefmt == 1) then
|
||
read(u, *) (data%rbnd(i), data%zbnd(i), i=1,nbnd)
|
||
else
|
||
read(u, '(5e16.9)') (data%rbnd(i), data%zbnd(i), i=1,nbnd)
|
||
end if
|
||
end if
|
||
|
||
! Load limiter data
|
||
if(nlim > 0) then
|
||
allocate(data%rlim(nlim), data%zlim(nlim))
|
||
if (params%ifreefmt == 1) then
|
||
read(u, *) (data%rlim(i), data%zlim(i), i=1,nlim)
|
||
else
|
||
read(u, '(5e16.9)') (data%rlim(i), data%zlim(i), i=1,nlim)
|
||
end if
|
||
end if
|
||
|
||
! End of G-EQDSK file
|
||
close(u)
|
||
|
||
! Build rv,zv,psinr arrays
|
||
zleft = zmid-0.5_wp_*dz
|
||
dr = dr/(nr-1)
|
||
dz = dz/(nz-1)
|
||
dps = one/(nr-1)
|
||
do i=1,nr
|
||
data%psinr(i) = (i-1)*dps
|
||
data%rv(i) = rleft + (i-1)*dr
|
||
end do
|
||
do i=1,nz
|
||
data%zv(i) = zleft + (i-1)*dz
|
||
end do
|
||
|
||
! Normalize psin
|
||
data%psia = psiedge - psiaxis
|
||
if(params%ipsinorm == 0) data%psin = (data%psin - psiaxis)/data%psia
|
||
|
||
end subroutine read_eqdsk
|
||
|
||
|
||
subroutine read_equil_an(filenm, ipass, data, err, unit)
|
||
! Reads the MHD equilibrium `data` in the analytical format
|
||
! from params%filenm.
|
||
! If given, the file is opened in the `unit` number.
|
||
!
|
||
! TODO: add format description
|
||
use gray_params, only : equilibrium_data
|
||
use utils, only : get_free_unit
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
character(len=*), intent(in) :: filenm
|
||
integer, intent(in) :: ipass
|
||
type(equilibrium_data), intent(out) :: data
|
||
integer, intent(out) :: err
|
||
integer, optional, intent(in) :: unit
|
||
|
||
! local variables
|
||
integer :: i, u, nlim
|
||
|
||
u = get_free_unit(unit)
|
||
|
||
open(file=filenm, status='old', action='read', unit=u, iostat=err)
|
||
if (err /= 0) then
|
||
call log_error('opening equilibrium file ('//trim(filenm)//') failed!', &
|
||
mod='equilibrium', proc='read_equil_an')
|
||
err = 1
|
||
return
|
||
end if
|
||
|
||
read(u, *) model%R0, model%z0, model%a
|
||
read(u, *) model%B0
|
||
read(u, *) model%q0, model%q1, model%alpha
|
||
|
||
if(ipass >= 2) then
|
||
! get size of boundary and limiter arrays and allocate them
|
||
read (u,*) nlim
|
||
if (allocated(data%rlim)) deallocate(data%rlim)
|
||
if (allocated(data%zlim)) deallocate(data%zlim)
|
||
|
||
! store boundary and limiter data
|
||
if(nlim > 0) then
|
||
allocate(data%rlim(nlim), data%zlim(nlim))
|
||
read(u,*) (data%rlim(i), data%zlim(i), i = 1, nlim)
|
||
end if
|
||
end if
|
||
close(u)
|
||
end subroutine read_equil_an
|
||
|
||
|
||
subroutine change_cocos(data, cocosin, cocosout, error)
|
||
! 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 gray_params, only : equilibrium_data
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_data), intent(inout) :: data
|
||
integer, intent(in) :: cocosin, cocosout
|
||
integer, intent(out), optional :: error
|
||
|
||
! 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%psia)
|
||
if (.not.psiincr) isign = -isign
|
||
bsign = sign(one, data%fpol(size(data%fpol)))
|
||
if (qpos .neqv. isign * bsign * data%qpsi(size(data%qpsi)) > zero) then
|
||
! Warning: sign inconsistency found among q, Ipla and Bref
|
||
data%qpsi = -data%qpsi
|
||
if (present(error)) error = 1
|
||
else
|
||
if (present(error)) error = 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 Bphi*Ip
|
||
if (qpos .neqv. qposout) data%qpsi = -data%qpsi
|
||
|
||
! psi and Ip signs don't change accordingly
|
||
if ((phiccw .eqv. phiccwout) .neqv. (psiincr .eqv. psiincrout)) &
|
||
data%psia = -data%psia
|
||
|
||
! Convert Wb to Wb/rad or viceversa
|
||
data%psia = data%psia * (2.0_wp_*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
|
||
implicit none
|
||
|
||
! 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 scale_equil(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, equilibrium_data
|
||
use gray_params, only : iequil
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(inout) :: params
|
||
type(equilibrium_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_φ).
|
||
|
||
if (iequil < 2) then
|
||
! 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
|
||
else
|
||
! Numeric data
|
||
|
||
! Apply signs
|
||
if (params%sgni /= 0) &
|
||
data%psia = sign(data%psia, -params%sgni*one)
|
||
if (params%sgnb /= 0) &
|
||
data%fpol = sign(data%fpol, +params%sgnb*one)
|
||
! 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
|
||
|
||
|
||
subroutine set_equil_spline(params, data, err)
|
||
! Computes splines for the MHD equilibrium data and stores them
|
||
! in their respective global variables, see the top of this file.
|
||
use const_and_precisions, only : zero, one
|
||
use gray_params, only : equilibrium_parameters, equilibrium_data
|
||
use gray_params, only : iequil
|
||
use reflections, only : inside
|
||
use utils, only : vmaxmin, vmaxmini
|
||
use logger, only : log_info
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(equilibrium_parameters), intent(in) :: params
|
||
type(equilibrium_data), intent(in) :: data
|
||
integer, intent(out) :: err
|
||
|
||
! local variables
|
||
integer :: nr, nz, nrest, nzest, npsest, nrz, npsi, nbnd, ibinf, ibsup
|
||
real(wp_) :: tension, rax0, zax0, psinoptmp, psinxptmp
|
||
real(wp_) :: rbmin, rbmax, rbinf, rbsup, r1, z1
|
||
real(wp_) :: psiant, psinop
|
||
real(wp_), dimension(size(data%psinr)) :: rhotn
|
||
real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi, wf
|
||
integer :: ixploc, info, i, j, ij
|
||
character(256) :: msg ! for log messages formatting
|
||
|
||
! compute array sizes
|
||
nr = size(data%rv)
|
||
nz = size(data%zv)
|
||
nrz = nr*nz
|
||
npsi = size(data%psinr)
|
||
nrest = nr + ksplp
|
||
nzest = nz + ksplp
|
||
npsest = npsi + ksplp
|
||
|
||
! length in m !!!
|
||
rmnm = data%rv(1)
|
||
rmxm = data%rv(nr)
|
||
zmnm = data%zv(1)
|
||
zmxm = data%zv(nz)
|
||
|
||
! Spline interpolation of ψ(R, z)
|
||
|
||
if (iequil>2) then
|
||
! 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
|
||
if (allocated(psi_spline%knots_x)) deallocate(psi_spline%knots_x)
|
||
if (allocated(psi_spline%knots_y)) deallocate(psi_spline%knots_y)
|
||
if (allocated(psi_spline%coeffs)) deallocate(psi_spline%coeffs)
|
||
allocate(psi_spline%knots_x(nrest), psi_spline%knots_y(nzest))
|
||
allocate(psi_spline%coeffs(nrz))
|
||
|
||
! determine number of valid grid points
|
||
nrz=0
|
||
do j=1,nz
|
||
do i=1,nr
|
||
if (nbnd.gt.0) then
|
||
if(.not.inside(data%rbnd,data%zbnd,nbnd,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
|
||
|
||
! store valid data
|
||
allocate(rv1d(nrz),zv1d(nrz),fvpsi(nrz),wf(nrz))
|
||
ij=0
|
||
do j=1,nz
|
||
do i=1,nr
|
||
if (nbnd.gt.0) then
|
||
if(.not.inside(data%rbnd,data%zbnd,nbnd,data%rv(i),data%zv(j))) cycle
|
||
else
|
||
if(data%psin(i,j).le.0.0d0) cycle
|
||
end if
|
||
ij=ij+1
|
||
rv1d(ij)=data%rv(i)
|
||
zv1d(ij)=data%zv(j)
|
||
fvpsi(ij)=data%psin(i,j)
|
||
wf(ij)=1.0d0
|
||
end do
|
||
end do
|
||
|
||
! Fit as a scattered set of points
|
||
! 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_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
|
||
else
|
||
! iequil==2: data are valid on the full R,z grid
|
||
|
||
! reshape 2D ψ array to 1D (transposed)
|
||
allocate(fvpsi(nrz))
|
||
fvpsi = reshape(transpose(data%psin), [nrz])
|
||
|
||
! compute spline coefficients
|
||
call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, &
|
||
range=[rmnm, rmxm, zmnm, zmxm], &
|
||
tension=params%ssplps, err=err)
|
||
|
||
! if failed, re-fit with an interpolating spline (zero tension)
|
||
if(err == -1) then
|
||
call psi_spline%init(data%rv, data%zv, fvpsi, nr, nz, &
|
||
range=[rmnm, rmxm, zmnm, zmxm], &
|
||
tension=zero)
|
||
err = 0
|
||
end if
|
||
deallocate(fvpsi)
|
||
end if
|
||
|
||
if (err /= 0) then
|
||
err = 2
|
||
return
|
||
end if
|
||
|
||
! compute spline coefficients for ψ(R,z) partial derivatives
|
||
call psi_spline%init_deriv(nr, nz, 1, 0) ! ∂ψ/∂R
|
||
call psi_spline%init_deriv(nr, nz, 0, 1) ! ∂ψ/∂z
|
||
call psi_spline%init_deriv(nr, nz, 1, 1) ! ∂²ψ/∂R∂z
|
||
call psi_spline%init_deriv(nr, nz, 2, 0) ! ∂²ψ/∂R²
|
||
call 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.
|
||
psi_domain%R = [rmnm, rmnm, rmxm, rmxm]
|
||
psi_domain%z = [zmnm, zmxm, zmxm, zmnm]
|
||
|
||
! Spline interpolation of F(ψ)
|
||
|
||
! give a small weight to the last point
|
||
allocate(wf(npsi))
|
||
wf(1:npsi-1) = 1
|
||
wf(npsi) = 1.0e2_wp_
|
||
call fpol_spline%init(data%psinr, data%fpol, npsi, range=[zero, one], &
|
||
weights=wf, tension=params%ssplf)
|
||
deallocate(wf)
|
||
|
||
! set vacuum value used outside 0 ≤ ψ ≤ 1 range
|
||
fpolas = fpol_spline%eval(data%psinr(npsi))
|
||
sgnbphi = sign(one ,fpolas)
|
||
|
||
! Re-normalize ψ_n after the spline computation
|
||
! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma
|
||
|
||
! Start with un-corrected ψ_n
|
||
psia = data%psia
|
||
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=min(size(data%rbnd), size(data%zbnd))
|
||
if (nbnd>0) then
|
||
call vmaxmini(data%zbnd,nbnd,zbinf,zbsup,ibinf,ibsup)
|
||
rbinf=data%rbnd(ibinf)
|
||
rbsup=data%rbnd(ibsup)
|
||
call vmaxmin(data%rbnd,nbnd,rbmin,rbmax)
|
||
else
|
||
zbinf=data%zv(2)
|
||
zbsup=data%zv(nz-1)
|
||
rbinf=data%rv((nr+1)/2)
|
||
rbsup=rbinf
|
||
rbmin=data%rv(2)
|
||
rbmax=data%rv(nr-1)
|
||
end if
|
||
|
||
! Search for exact location of the magnetic axis
|
||
rax0=data%rax
|
||
zax0=data%zax
|
||
call points_ox(rax0,zax0,rmaxis,zmaxis,psinoptmp,info)
|
||
|
||
write (msg, '("O-point found:", 3(x,a,"=",g0.3))') &
|
||
'r', rmaxis, 'z', zmaxis, 'ψ', psinoptmp
|
||
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
|
||
|
||
! search for X-point if params%ixp /= 0
|
||
|
||
ixploc = params%ixp
|
||
if(ixploc/=0) then
|
||
if(ixploc<0) then
|
||
call points_ox(rbinf,zbinf,r1,z1,psinxptmp,info)
|
||
if(psinxptmp/=-1.0_wp_) then
|
||
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
|
||
'r', r1, 'z', z1, 'ψ', psinxptmp
|
||
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
|
||
|
||
zbinf=z1
|
||
psinop=psinoptmp
|
||
psiant=psinxptmp-psinop
|
||
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
|
||
zbsup=z1
|
||
else
|
||
ixploc=0
|
||
end if
|
||
else
|
||
call points_ox(rbsup,zbsup,r1,z1,psinxptmp,info)
|
||
if(psinxptmp.ne.-1.0_wp_) then
|
||
write (msg, '("X-point found:", 3(x,a,"=",g0.3))') &
|
||
'r', r1, 'z', z1, 'ψ', psinxptmp
|
||
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
|
||
|
||
zbsup=z1
|
||
psinop=psinoptmp
|
||
psiant=psinxptmp-psinop
|
||
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
|
||
zbinf=z1
|
||
else
|
||
ixploc=0
|
||
end if
|
||
end if
|
||
end if
|
||
|
||
if (ixploc==0) then
|
||
psinop=psinoptmp
|
||
psiant=one-psinop
|
||
|
||
! Find upper horizontal tangent point
|
||
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbsup),r1,z1,one,info)
|
||
zbsup=z1
|
||
rbsup=r1
|
||
|
||
! Find lower horizontal tangent point
|
||
call points_tgo(rmaxis,0.5_wp_*(zmaxis+zbinf),r1,z1,one,info)
|
||
zbinf=z1
|
||
rbinf=r1
|
||
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='equilibrium', proc='set_equil_spline')
|
||
end if
|
||
|
||
! Adjust all the B-spline coefficients
|
||
! Note: since ψ_n(R,z) = Σ_ij c_ij B_i(R)B_j(z), to correct ψ_n
|
||
! it's enough to correct (shift and scale) the c_ij
|
||
psi_spline%coeffs = (psi_spline%coeffs - psinop) / psiant
|
||
|
||
! Same for the derivatives
|
||
psi_spline%partial(1,0)%ptr = psi_spline%partial(1,0)%ptr / psiant
|
||
psi_spline%partial(0,1)%ptr = psi_spline%partial(0,1)%ptr / psiant
|
||
psi_spline%partial(2,0)%ptr = psi_spline%partial(2,0)%ptr / psiant
|
||
psi_spline%partial(0,2)%ptr = psi_spline%partial(0,2)%ptr / psiant
|
||
psi_spline%partial(1,1)%ptr = psi_spline%partial(1,1)%ptr / psiant
|
||
|
||
! Do the opposite scaling to preserve un-normalised values
|
||
! Note: this is only used for the poloidal magnetic field
|
||
psia = psia * psiant
|
||
|
||
! Save Bt value on axis (required in flux_average and used in Jcd def)
|
||
! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def)
|
||
|
||
call pol_curr(zero, btaxis)
|
||
btaxis = btaxis/rmaxis
|
||
btrcen = fpolas/data%rvac
|
||
rcen = data%rvac
|
||
|
||
write (msg, '(2(a,g0.3))') 'Bt_center=', btrcen, ' Bt_axis=', btaxis
|
||
call log_info(msg, mod='equilibrium', proc='set_equil_spline')
|
||
|
||
! Compute ρ_p/ρ_t mapping based on the input q profile
|
||
call setqphi_num(data%psinr, abs(data%qpsi), abs(psia), rhotn)
|
||
call set_rho_spline(sqrt(data%psinr), rhotn)
|
||
|
||
! Compute the domain of the ψ mapping
|
||
psi_domain%R = data%rbnd
|
||
psi_domain%z = data%zbnd
|
||
call rescale_boundary(psi_domain, O=[rmaxis, zmaxis], t0=1.5_wp_)
|
||
|
||
end subroutine set_equil_spline
|
||
|
||
|
||
subroutine rescale_boundary(cont, 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.
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
type(contour), intent(inout) :: cont ! (R,z) contour
|
||
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
|
||
|
||
function s(t)
|
||
! Rescriction of ψ(R, z) on the line Q(t) = O + tN
|
||
implicit none
|
||
|
||
real(wp_), intent(in) :: t
|
||
real(wp_) :: s, Q(2)
|
||
|
||
Q = O + t * N
|
||
s = psi_spline%eval(Q(1), Q(2))
|
||
end function
|
||
|
||
end subroutine rescale_boundary
|
||
|
||
|
||
subroutine scatterspl(x,y,z,w,n,kspl,sspl,xmin,xmax,ymin,ymax, &
|
||
tx,nknt_x,ty,nknt_y,coeff,ierr)
|
||
! Computes the spline interpolation of a surface when
|
||
! the data points are irregular, i.e. not on a uniform grid
|
||
use const_and_precisions, only : comp_eps
|
||
use dierckx, only : surfit
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: n
|
||
real(wp_), dimension(n), intent(in) :: x, y, z
|
||
real(wp_), dimension(n), intent(in) :: w
|
||
integer, intent(in) :: kspl
|
||
real(wp_), intent(in) :: sspl
|
||
real(wp_), intent(in) :: xmin, xmax, ymin, ymax
|
||
real(wp_), dimension(nknt_x), intent(inout) :: tx
|
||
real(wp_), dimension(nknt_y), intent(inout) :: ty
|
||
integer, intent(inout) :: nknt_x, nknt_y
|
||
real(wp_), dimension(nknt_x*nknt_y), intent(out) :: coeff
|
||
integer, intent(out) :: ierr
|
||
|
||
! local variables
|
||
integer :: iopt
|
||
real(wp_) :: resid
|
||
integer :: u,v,km,ne,b1,b2,lwrk1,lwrk2,kwrk,nxest,nyest
|
||
real(wp_), dimension(:), allocatable :: wrk1, wrk2
|
||
integer, dimension(:), allocatable :: iwrk
|
||
|
||
nxest=nknt_x
|
||
nyest=nknt_y
|
||
ne = max(nxest,nyest)
|
||
|
||
km = kspl+1
|
||
u = nxest-km
|
||
v = nyest-km
|
||
b1 = kspl*min(u,v)+kspl+1
|
||
b2 = (kspl+1)*min(u,v)+1
|
||
lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(n+ne)+ne-2*kspl)+b2+1
|
||
lwrk2 = u*v*(b2+1)+b2
|
||
kwrk = n+(nknt_x-2*kspl-1)*(nknt_y-2*kspl-1)
|
||
allocate(wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk))
|
||
|
||
iopt=0
|
||
call surfit(iopt,n,x,y,z,w,xmin,xmax,ymin,ymax,kspl,kspl, &
|
||
sspl,nxest,nyest,ne,comp_eps,nknt_x,tx,nknt_y,ty, &
|
||
coeff,resid,wrk1,lwrk1,wrk2,lwrk2,iwrk,kwrk,ierr)
|
||
|
||
deallocate(wrk1,wrk2,iwrk)
|
||
|
||
end subroutine scatterspl
|
||
|
||
|
||
subroutine setqphi_num(psinq,q,psia,rhotn)
|
||
! Computes the spline of the safety factor q(ψ)
|
||
use const_and_precisions, only : pi
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
real(wp_), dimension(:), intent(in) :: psinq,q
|
||
real(wp_), intent(in) :: psia
|
||
real(wp_), dimension(:), intent(out), optional :: rhotn
|
||
|
||
! local variables
|
||
real(wp_), dimension(size(q)) :: phit
|
||
real(wp_) :: dx
|
||
integer :: k
|
||
|
||
call q_spline%init(psinq, q, size(q))
|
||
|
||
! Toroidal flux as Φ(ψ) = 2π ∫q(ψ)dψ
|
||
phit(1)=0
|
||
do k=1,q_spline%ndata-1
|
||
dx=q_spline%data(k+1)-q_spline%data(k)
|
||
phit(k+1)=phit(k) + dx*(q_spline%coeffs(k,1) + dx*(q_spline%coeffs(k,2)/2 + &
|
||
dx*(q_spline%coeffs(k,3)/3 + dx* q_spline%coeffs(k,4)/4) ) )
|
||
end do
|
||
phitedge=phit(q_spline%ndata)
|
||
if(present(rhotn)) rhotn(1:q_spline%ndata)=sqrt(phit/phitedge)
|
||
phitedge=2*pi*psia*phitedge
|
||
end subroutine setqphi_num
|
||
|
||
|
||
subroutine set_equil_an
|
||
! Computes the analytical equilibrium data and stores them
|
||
! in their respective global variables, see the top of this file.
|
||
use const_and_precisions, only : pi, one
|
||
|
||
implicit none
|
||
|
||
real(wp_) :: dq, gamma
|
||
|
||
btaxis = model%B0
|
||
rmaxis = model%R0
|
||
zmaxis = model%z0
|
||
btrcen = model%B0
|
||
rcen = model%R0
|
||
zbinf = model%z0 - model%a
|
||
zbsup = model%z0 + model%a
|
||
sgnbphi = sign(one, model%B0)
|
||
|
||
rmxm = model%R0 + model%a
|
||
rmnm = model%R0 - model%a
|
||
zmxm = zbsup
|
||
zmnm = zbinf
|
||
|
||
! Toroidal flux at r=a:
|
||
!
|
||
! Φ(a) = B₀πa² 2γ/(γ + 1)
|
||
!
|
||
! where γ=1/√(1-ε²),
|
||
! ε=a/R₀ is the tokamak aspect ratio
|
||
gamma = 1/sqrt(1 - (model%a/model%R0)**2)
|
||
phitedge = model%B0 * pi * model%a**2 * 2*gamma/(gamma + 1)
|
||
|
||
! 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)
|
||
dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
|
||
psia = 1/(2*pi) * phitedge / (model%q0 + dq)
|
||
end subroutine set_equil_an
|
||
|
||
|
||
subroutine set_rho_spline(rhop, rhot)
|
||
! Computes the splines for converting between the poloidal (ρ_p)
|
||
! and toroidal (ρ_t) normalised radii
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
real(wp_), dimension(:), intent(in) :: rhop, rhot
|
||
|
||
call rhop_spline%init(rhot, rhop, size(rhop))
|
||
call rhot_spline%init(rhop, rhot, size(rhot))
|
||
end subroutine set_rho_spline
|
||
|
||
|
||
subroutine pol_flux(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.
|
||
use gray_params, only : iequil
|
||
use reflections, only : inside
|
||
use const_and_precisions, only : one, pi
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
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 !
|
||
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²)²
|
||
|
||
if (iequil < 2) then
|
||
! Analytical model
|
||
!
|
||
! The normalised poloidal flux ψ_n(R, z) is computed as follows:
|
||
! 1. ψ_n = ρ_p²
|
||
! 2. ρ_p = ρ_p(ρ_t), using `frhopol`, which in turns uses q(ψ)
|
||
! 3. ρ_t = √Φ_n
|
||
! 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:
|
||
!
|
||
! Φ(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 > 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
|
||
|
||
gamma = 1 / sqrt(1 - (r_g/model%R0)**2)
|
||
phi_n = model%B0 * pi*r_g**2 * 2*gamma/(gamma + 1) / phitedge
|
||
rho_t = sqrt(phi_n)
|
||
rho_p = frhopol(rho_t)
|
||
|
||
! For ∇Φ_n and ∇∇Φ_n we also need:
|
||
!
|
||
! ∂Φ∂(r²) = B₀π γ(r)
|
||
! ∂²Φ∂(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²)
|
||
! where ∇(r²) = 2[(R-R₀), (z-z₀)]
|
||
dphidr = dphidr2 * 2*(R - model%R0)
|
||
dphidz = dphidr2 * 2*(z - model%z0)
|
||
|
||
! ∇∇Φ_n = ∇[∂Φ_n/∂(r²)] ∇(r²) + ∂Φ_n/∂(r²) ∇∇(r²)
|
||
! = ∂²Φ_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)
|
||
|
||
! ψ_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 = 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
|
||
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 = 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:
|
||
!
|
||
! ∇∇ψ_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
|
||
else
|
||
! Numerical data
|
||
if (inside(psi_domain%R, psi_domain%z, psi_domain%npoints, 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)
|
||
else
|
||
! Outside
|
||
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 if
|
||
end subroutine pol_flux
|
||
|
||
|
||
subroutine pol_curr(psi_n, fpol, dfpol)
|
||
! Computes the poloidal current function F(ψ_n)
|
||
! and (optionally) its derivative dF/dψ_n given ψ_n.
|
||
use gray_params, only : iequil
|
||
|
||
implicit none
|
||
|
||
! function arguments
|
||
real(wp_), intent(in) :: psi_n ! normalised poloidal flux
|
||
real(wp_), intent(out) :: fpol ! poloidal current
|
||
real(wp_), intent(out), optional :: dfpol ! derivative
|
||
|
||
if (iequil < 2) then
|
||
! Analytical model
|
||
! F(ψ) = B₀⋅R₀, a constant
|
||
fpol = model%B0 * model%R0
|
||
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
|
||
end if
|
||
end if
|
||
end subroutine pol_curr
|
||
|
||
|
||
function frhotor(rho_p)
|
||
! Converts from poloidal (ρ_p) to toroidal (ρ_t) normalised radius
|
||
use gray_params, only : iequil
|
||
|
||
implicit none
|
||
|
||
! function arguments
|
||
real(wp_), intent(in) :: rho_p
|
||
real(wp_) :: frhotor
|
||
|
||
if (iequil < 2) then
|
||
! Analytical model
|
||
block
|
||
! 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)
|
||
real(wp_) :: dq
|
||
dq = (model%q1 - model%q0) / (model%alpha/2 + 1)
|
||
frhotor = rho_p * sqrt((model%q0 + dq*rho_p**model%alpha) &
|
||
/ (model%q0 + dq))
|
||
end block
|
||
else
|
||
! Numerical data
|
||
frhotor = rhot_spline%eval(rho_p)
|
||
end if
|
||
end function frhotor
|
||
|
||
|
||
function frhopol(rho_t)
|
||
! Converts from toroidal (ρ_t) to poloidal (ρ_p) normalised radius
|
||
use gray_params, only : iequil
|
||
use const_and_precisions, only : comp_eps
|
||
|
||
implicit none
|
||
|
||
! function arguments
|
||
real(wp_), intent(in) :: rho_t
|
||
real(wp_) :: frhopol
|
||
|
||
if (iequil < 2) then
|
||
! 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
|
||
|
||
real(wp_) :: rho_p(1), fvec(1), fjac(1,1), wa(7)
|
||
integer :: info
|
||
|
||
rho_p = [rho_t] ! 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)
|
||
frhopol = rho_p(1)
|
||
end block
|
||
else
|
||
! Numerical data
|
||
frhopol = rhop_spline%eval(rho_t)
|
||
end if
|
||
|
||
contains
|
||
|
||
subroutine equation(n, x, f, df, ldf, flag)
|
||
! The equation to solve: f(x) = ρ_t(x) - ρ_t₀ = 0
|
||
implicit none
|
||
|
||
! 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) = frhotor(x(1)) - rho_t
|
||
else
|
||
! return f'(x), computed numerically
|
||
df(1,1) = (frhotor(x(1) + e) - frhotor(x(1) - e)) / (2*e)
|
||
end if
|
||
end subroutine
|
||
|
||
end function frhopol
|
||
|
||
|
||
function fq(psin)
|
||
! Computes the safety factor q as a function of the normalised
|
||
! poloidal flux ψ.
|
||
!
|
||
! Note: this returns the absolute value of q.
|
||
use gray_params, only : iequil
|
||
|
||
implicit none
|
||
|
||
! function arguments
|
||
real(wp_), intent(in) :: psin
|
||
real(wp_) :: fq
|
||
|
||
if (iequil < 2) then
|
||
! Analytical model
|
||
! The safety factor is a power law in ρ_p:
|
||
! q(ρ_p) = q₀ + (q₁-q₀) ρ_p^α
|
||
block
|
||
real(wp_) :: rho_p
|
||
rho_p = sqrt(psin)
|
||
fq = abs(model%q0 + (model%q1 - model%q0) * rho_p**model%alpha)
|
||
end block
|
||
else
|
||
! Numerical data
|
||
fq = q_spline%eval(psin)
|
||
end if
|
||
end function fq
|
||
|
||
|
||
subroutine bfield(R, z, B_R, B_z, B_phi)
|
||
! Computes the magnetic field as a function of
|
||
! (R, z) in cylindrical coordinates
|
||
!
|
||
! Note: all output arguments are optional.
|
||
use gray_params, only : iequil
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
real(wp_), intent(in) :: R, z
|
||
real(wp_), intent(out), optional :: B_R, B_z, B_phi
|
||
|
||
! local variables
|
||
real(wp_) :: psi_n, fpol, dpsidr, dpsidz
|
||
|
||
call pol_flux(R, z, psi_n, dpsidr, dpsidz)
|
||
call 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
|
||
!
|
||
if (present(B_R)) B_R = - 1/R * dpsidz * psia
|
||
if (present(B_z)) B_z = + 1/R * dpsidr * psia
|
||
if (present(B_phi)) B_phi = fpol / R
|
||
end subroutine bfield
|
||
|
||
|
||
function tor_curr(R, z) result(J_phi)
|
||
! Computes the toroidal current J_φ as a function of (R, z)
|
||
use const_and_precisions, only : mu0_
|
||
|
||
implicit none
|
||
|
||
! function arguments
|
||
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 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 * psia
|
||
dB_zdR = + 1/R * (ddpsidrr - 1/R * dpsidr) * psia
|
||
J_phi = 1/mu0_ * (dB_Rdz - dB_zdR)
|
||
end function tor_curr
|
||
|
||
|
||
function tor_curr_psi(psi_n) result(J_phi)
|
||
! Computes the toroidal current J_φ as a function of ψ
|
||
|
||
implicit none
|
||
|
||
! function arguments
|
||
real(wp_), intent(in) :: psi_n
|
||
real(wp_) :: J_phi
|
||
|
||
! local constants
|
||
real(wp_), parameter :: eps = 1.e-4_wp_
|
||
|
||
! local variables
|
||
real(wp_) :: R1, R2
|
||
|
||
call psi_raxis(max(eps, psi_n), R1, R2)
|
||
J_phi = tor_curr(R2, zmaxis)
|
||
end function tor_curr_psi
|
||
|
||
|
||
subroutine psi_raxis(psin,r1,r2)
|
||
use gray_params, only : iequil
|
||
use dierckx, only : profil, sproota
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
real(wp_) :: psin,r1,r2
|
||
|
||
! local constants
|
||
integer, parameter :: mest=4
|
||
|
||
! local variables
|
||
integer :: iopt,ier,m
|
||
real(wp_) :: zc,val
|
||
real(wp_), dimension(mest) :: zeroc
|
||
real(wp_), dimension(psi_spline%nknots_x) :: czc
|
||
character(64) :: msg
|
||
|
||
if (iequil < 2) then
|
||
! Analytical model
|
||
val = sqrt(psin)
|
||
r1 = model%R0 - val*model%a
|
||
r2 = model%R0 + val*model%a
|
||
else
|
||
! Numerical data
|
||
iopt=1
|
||
zc=zmaxis
|
||
call profil(iopt, psi_spline%knots_x, psi_spline%nknots_x, &
|
||
psi_spline%knots_y, psi_spline%nknots_y, &
|
||
psi_spline%coeffs, kspl, kspl, zc, &
|
||
psi_spline%nknots_x, czc, ier)
|
||
if (ier > 0) then
|
||
write (msg, '("profil failed with error ",g0)') ier
|
||
call log_error(msg, mod='equilibrium', proc='psi_raxis')
|
||
end if
|
||
|
||
call sproota(psin, psi_spline%knots_x, psi_spline%nknots_x, &
|
||
czc, zeroc, mest, m, ier)
|
||
r1=zeroc(1)
|
||
r2=zeroc(2)
|
||
end if
|
||
end subroutine psi_raxis
|
||
|
||
|
||
subroutine points_ox(rz,zz,rf,zf,psinvf,info)
|
||
! Finds the location of the O,X points
|
||
use const_and_precisions, only : comp_eps
|
||
use minpack, only : hybrj1
|
||
use logger, only : log_error, log_debug
|
||
|
||
implicit none
|
||
|
||
! local constants
|
||
integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
|
||
|
||
! arguments
|
||
real(wp_), intent(in) :: rz,zz
|
||
real(wp_), intent(out) :: rf,zf,psinvf
|
||
integer, intent(out) :: info
|
||
|
||
! local variables
|
||
real(wp_) :: tol
|
||
real(wp_), dimension(n) :: xvec,fvec
|
||
real(wp_), dimension(lwa) :: wa
|
||
real(wp_), dimension(ldfjac,n) :: fjac
|
||
character(256) :: msg
|
||
|
||
xvec(1)=rz
|
||
xvec(2)=zz
|
||
tol = sqrt(comp_eps)
|
||
call hybrj1(fcnox,n,xvec,fvec,fjac,ldfjac,tol,info,wa,lwa)
|
||
if(info /= 1) then
|
||
write (msg, '("O,X coordinates:",2(x,", ",g0.3))') xvec
|
||
call log_debug(msg, mod='equilibrium', proc='points_ox')
|
||
write (msg, '("hybrj1 failed with error ",g0)') info
|
||
call log_error(msg, mod='equilibrium', proc='points_ox')
|
||
end if
|
||
rf=xvec(1)
|
||
zf=xvec(2)
|
||
call pol_flux(rf, zf, psinvf)
|
||
end subroutine points_ox
|
||
|
||
|
||
subroutine fcnox(n,x,fvec,fjac,ldfjac,iflag)
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: n,iflag,ldfjac
|
||
real(wp_), dimension(n), intent(in) :: x
|
||
real(wp_), dimension(n), intent(inout) :: fvec
|
||
real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
|
||
|
||
! local variables
|
||
real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz
|
||
character(64) :: msg
|
||
|
||
select case(iflag)
|
||
case(1)
|
||
call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz)
|
||
fvec(1) = dpsidr
|
||
fvec(2) = dpsidz
|
||
case(2)
|
||
call pol_flux(x(1), x(2), ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz, &
|
||
ddpsidrz=ddpsidrz)
|
||
fjac(1,1) = ddpsidrr
|
||
fjac(1,2) = ddpsidrz
|
||
fjac(2,1) = ddpsidrz
|
||
fjac(2,2) = ddpsidzz
|
||
case default
|
||
write (msg, '("invalid iflag: ",g0)')
|
||
call log_error(msg, mod='equilibrium', proc='fcnox')
|
||
end select
|
||
end subroutine fcnox
|
||
|
||
|
||
subroutine points_tgo(rz,zz,rf,zf,psin0,info)
|
||
use const_and_precisions, only : comp_eps
|
||
use minpack, only : hybrj1mv
|
||
use logger, only : log_error, log_debug
|
||
|
||
implicit none
|
||
|
||
! local constants
|
||
integer, parameter :: n=2,ldfjac=n,lwa=(n*(n+13))/2
|
||
|
||
! arguments
|
||
real(wp_), intent(in) :: rz,zz,psin0
|
||
real(wp_), intent(out) :: rf,zf
|
||
integer, intent(out) :: info
|
||
character(256) :: msg
|
||
|
||
! local variables
|
||
real(wp_) :: tol
|
||
real(wp_), dimension(n) :: xvec,fvec,f0
|
||
real(wp_), dimension(lwa) :: wa
|
||
real(wp_), dimension(ldfjac,n) :: fjac
|
||
|
||
xvec(1)=rz
|
||
xvec(2)=zz
|
||
f0(1)=psin0
|
||
f0(2)=0.0_wp_
|
||
tol = sqrt(comp_eps)
|
||
call hybrj1mv(fcntgo,n,xvec,f0,fvec,fjac,ldfjac,tol,info,wa,lwa)
|
||
if(info /= 1) then
|
||
write (msg, '("R,z coordinates:",5(x,g0.3))') xvec, rz, zz, psin0
|
||
call log_debug(msg, mod='equilibrium', proc='points_tgo')
|
||
write (msg, '("hybrj1mv failed with error ",g0)') info
|
||
call log_error(msg, mod='equilibrium', proc='points_tgo')
|
||
end if
|
||
rf=xvec(1)
|
||
zf=xvec(2)
|
||
end
|
||
|
||
|
||
subroutine fcntgo(n,x,f0,fvec,fjac,ldfjac,iflag)
|
||
use logger, only : log_error
|
||
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
integer, intent(in) :: n,ldfjac,iflag
|
||
real(wp_), dimension(n), intent(in) :: x,f0
|
||
real(wp_), dimension(n), intent(inout) :: fvec
|
||
real(wp_), dimension(ldfjac,n), intent(inout) :: fjac
|
||
|
||
! local variables
|
||
real(wp_) :: psinv,dpsidr,dpsidz,ddpsidrr,ddpsidrz
|
||
character(64) :: msg
|
||
|
||
select case(iflag)
|
||
case(1)
|
||
call pol_flux(x(1), x(2), psinv, dpsidr)
|
||
fvec(1) = psinv-f0(1)
|
||
fvec(2) = dpsidr-f0(2)
|
||
case(2)
|
||
call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz, &
|
||
ddpsidrr=ddpsidrr, ddpsidrz=ddpsidrz)
|
||
fjac(1,1) = dpsidr
|
||
fjac(1,2) = dpsidz
|
||
fjac(2,1) = ddpsidrr
|
||
fjac(2,2) = ddpsidrz
|
||
case default
|
||
write (msg, '("invalid iflag: ",g0)')
|
||
call log_error(msg, mod='equilibrium', proc='fcntgo')
|
||
end select
|
||
end subroutine fcntgo
|
||
|
||
|
||
subroutine unset_equil_spline
|
||
! Unsets the splines global variables, see the top of this file.
|
||
implicit none
|
||
|
||
call fpol_spline%deinit
|
||
call psi_spline%deinit
|
||
call q_spline%deinit
|
||
call rhop_spline%deinit
|
||
call rhot_spline%deinit
|
||
|
||
if (allocated(psi_domain%R)) deallocate(psi_domain%R, psi_domain%z)
|
||
|
||
end subroutine unset_equil_spline
|
||
|
||
end module equilibrium
|