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

move flux averages inside gray_equil

Browse Source
This commit is contained in:
Michele Guerini Rocco 2024-10-23 15:15:37 +02:00 committed by rnhmjoj
parent 8c144c3892
commit 7477fffb43
Signed by: rnhmjoj
GPG Key ID: BFBAF4C975F76450
8 changed files with 527 additions and 440 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
src/eccd.f90
View File

@ -138,12 +138,13 @@ contains
eccdpar(5)=fp0s
end subroutine setcdcoeff_cohen
subroutine setcdcoeff_ncl(zeff,rbx,fc,amu,rhop,cst2,eccdpar)
use magsurf_data, only : cd_eff
subroutine setcdcoeff_ncl(zeff,rbx,fc,amu,rhop,cd_eff,cst2,eccdpar)
use splines, only : spline_2d
use dierckx, only : profil
use logger, only : log_warning
real(wp_), intent(in) :: zeff,rbx,fc,amu,rhop
type(spline_2d), intent(in) :: cd_eff
real(wp_), intent(out) :: cst2
real(wp_), dimension(:), allocatable, intent(out) :: eccdpar
real(wp_), dimension(cd_eff%nknots_y) :: chlm

23
src/gray_core.f90
View File

@ -15,13 +15,13 @@ contains
use gray_equil, only : abstract_equil
use gray_plasma, only : abstract_plasma
use gray_tables, only : init_tables, store_ec_profiles, store_ray_data, &
store_beam_shape, find_flux_surfaces, &
kinetic_profiles, ec_resonance, inputs_maps
store_beam_shape, find_flux_surfaces, &
flux_surfaces, kinetic_profiles, flux_averages, &
ec_resonance, inputs_maps
use gray_errors, only : gray_error, is_critical, has_new_errors, &
print_err_raytracing, print_err_ecrh_cd
use beams, only : xgygcoeff, launchangles2n
use beamdata, only : pweight
use magsurf_data, only : compute_flux_averages
use pec, only : pec_init, spec, postproc_profiles, dealloc_pec, &
rhop_tab, rhot_tab
use utils, only : vmaxmin
@ -119,9 +119,6 @@ contains
call launchangles2n(params%antenna, anv0)
if (params%equilibrium%iequil /= EQ_VACUUM) then
! Initialise the magsurf_data module
call compute_flux_averages(params, equil, results%tables)
! Initialise the output profiles
call pec_init(params%output, equil, rhout)
end if
@ -139,13 +136,14 @@ contains
! Pre-determinted tables
results%tables%kinetic_profiles = kinetic_profiles(params, equil, plasma)
results%tables%flux_averages = flux_averages(params, equil)
results%tables%flux_surfaces = flux_surfaces(params, equil)
results%tables%ec_resonance = ec_resonance(params, equil, bres)
results%tables%inputs_maps = inputs_maps(params, equil, plasma, bres, xgcn)
! print ψ surface for q=3/2 and q=2/1
! print ψ rational surfaces
call find_flux_surfaces(qvals=[1.0_wp_, 1.5_wp_, 2.0_wp_], &
params=params, equil=equil, &
tbl=results%tables%flux_surfaces)
equil=equil, tbl=results%tables%flux_surfaces)
! print initial position
write (msg, '("initial position:",3(x,g0.3))') params%antenna%pos
@ -1659,7 +1657,6 @@ contains
use dispersion, only : harmnumber, warmdisp
use eccd, only : setcdcoeff, eccdeff, fjch0, fjch, fjncl
use gray_errors, only : gray_error, negative_absorption, raise_error
use magsurf_data, only : fluxval
! subroutine arguments
@ -1773,7 +1770,9 @@ contains
ithn = 1
if (nlarmor > nhmin) ithn = 2
rhop = sqrt(psinv)
call fluxval(rhop, rri=rrii, rbav=rbavi, bmn=bmni, bmx=bmxi, fc=fci)
call equil%flux_average(rhop, R_J=rrii, B_avg=rbavi, &
B_min=bmni, B_max=bmxi, f_c=fci)
rbavi = rbavi / bmni
Btot = Y*Bres
rbn = Btot/bmni
rbx = Btot/bmxi
@ -1791,7 +1790,7 @@ contains
ithn, cst2, fjch0, eccdpar, effjcd, iokhawa, ierrcd)
case default
! Neoclassical model
call setcdcoeff(zeff, rbx, fci, mu, rhop, cst2, eccdpar)
call setcdcoeff(zeff, rbx, fci, mu, rhop, equil%spline_cd_eff, cst2, eccdpar)
call eccdeff(Y, Npl, Npr%re, density, mu, e, nhmin, nhmax, &
ithn, cst2, fjncl, eccdpar, effjcd, iokhawa, ierrcd)
end select

424
src/gray_equil.f90
View File

@ -24,14 +24,25 @@ module gray_equil
real(wp_) :: r_range(2) = [-comp_huge, comp_huge] ! R range of the equilibrium domain
real(wp_) :: z_range(2) = [-comp_huge, comp_huge] ! z range of the equilibrium domain
real(wp_) :: z_boundary(2) = [0, 0] ! 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
@ -77,19 +88,17 @@ module gray_equil
real(wp_) :: rho_out
end function rho_conv_fun
subroutine flux_contour_sub(self, psi0, R_min, R, z, &
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_min is a value such that R>R_min for any contour point
! - (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(in) :: R_min
real(wp_), intent(out) :: R(:), z(:)
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
end subroutine flux_contour_sub
@ -274,6 +283,49 @@ contains
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
!
@ -517,14 +569,13 @@ contains
end function analytic_tor2pol
pure subroutine analytic_flux_contour(self, psi0, R_min, R, z, &
R_hi, z_hi, R_lo, z_lo)
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(in) :: R_min
real(wp_), intent(out) :: R(:), z(:)
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
@ -533,7 +584,6 @@ contains
real(wp_) :: r_g ! geometric minor radius
real(wp_) :: theta ! geometric poloidal angle
unused(R_min)
unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo)
n = size(R)
@ -636,8 +686,8 @@ contains
end function numeric_tor2pol
subroutine numeric_flux_contour(self, psi0, R_min, R, z, &
R_hi, z_hi, R_lo, z_lo)
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
@ -645,7 +695,6 @@ contains
! subroutine arguments
class(numeric_equil), intent(in) :: self
real(wp_), intent(in) :: psi0
real(wp_), intent(in) :: R_min
real(wp_), intent(out) :: R(:), z(:)
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
@ -654,13 +703,15 @@ contains
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), zeroc(4)
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)
@ -671,26 +722,31 @@ contains
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, a, g0)') &
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
call sproota(psi0, s%knots_x, s%nknots_x, czc, zeroc, 4, m, ier)
! Find the solutions of s(R) = ψ₀
call sproota(psi0, s%knots_x, s%nknots_x, czc, sols, 3, m, ier)
end associate
if (zeroc(1) > R_min) then
R(i) = zeroc(1)
R(n+1-i) = zeroc(2)
! 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) = zeroc(2)
R(n+1-i) = zeroc(3)
R(i) = sols(2)
R(n+1-i) = sols(3)
end if
z(i) = zc
z(n+1-i) = zc
@ -1223,16 +1279,15 @@ contains
end function vacuum_conv
pure subroutine vacuum_flux_contour(self, psi0, R_min, R, z, &
R_hi, z_hi, R_lo, z_lo)
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(in) :: R_min
real(wp_), intent(out) :: R(:), z(:)
real(wp_), intent(inout) :: R_hi, z_hi, R_lo, z_lo
unused(self); unused(psi0); unused(R_min)
unused(self); unused(psi0);
unused(R_hi); unused(z_hi); unused(R_lo); unused(z_lo)
! flux surfaces are undefined
@ -1295,6 +1350,8 @@ contains
end select
call equil%init_flux_splines
end subroutine load_equil
@ -1644,4 +1701,327 @@ contains
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

2
src/gray_jetto1beam.f90
View File

@ -88,6 +88,8 @@ subroutine gray_jetto1beam(ijetto, mr, mz, r, z, psin, psia, rax, zax, &
! Compute splines
call equil%init(params%equilibrium, data, err)
if (err /= 0) return
call equil%init_flux_splines
end block init_equilibrium
! Compute splines of kinetic profiles

120
src/gray_tables.f90
View File

@ -11,7 +11,8 @@ module gray_tables
public find_flux_surfaces, store_flux_surface ! Filling the main tables
public store_ray_data, store_ec_profiles !
public store_beam_shape !
public kinetic_profiles, ec_resonance, inputs_maps ! Extra tables
public kinetic_profiles, flux_averages ! Extra tables
public flux_surfaces, ec_resonance, inputs_maps !
contains
@ -43,16 +44,6 @@ contains
type(gray_parameters), intent(in) :: params
type(gray_tables), target, intent(out) :: tbls
call tbls%flux_averages%init( &
title='flux-averages', id=56, active=is_active(params, 56), &
labels=[character(64) :: 'ρ_p', 'ρ_t', 'B_avg', &
'B_max', 'B_min', 'area', 'vol', 'I_pl', &
'J_φ_avg', 'fc', 'ratio_Ja', 'ratio_Jb', 'q'])
call tbls%flux_surfaces%init( &
title='flux-surfaces', id=71, active=is_active(params, 71), &
labels=[character(64) :: 'i', 'ψ_n', 'R', 'z'])
call tbls%ec_profiles%init( &
title='ec-profiles', id=48, active=is_active(params, 48), &
labels=[character(64) :: 'ρ_p', 'ρ_t', 'J_φ', &
@ -146,8 +137,7 @@ contains
rho_p = i * drho_p
rho_t = equil%pol2tor(rho_p)
psi_n = rho_p**2
!call equil%flux_average(rho_p, J_phi=J_phi)
J_phi = 0
call equil%flux_average(rho_p, J_phi_avg=J_phi)
call plasma%density(psi_n, dens, ddens)
call tbl%append([psi_n, rho_t, dens, plasma%temp(psi_n), &
equil%safety(psi_n), J_phi])
@ -156,6 +146,104 @@ contains
end function kinetic_profiles
function flux_averages(params, equil) result(tbl)
! Generates the flux surface averages table
use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
! function arguments
type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(table) :: tbl
! local variables
integer, parameter :: n_grid = 100
integer :: i
real(wp_) :: rho_t, rho_p, drho_p
real(wp_) :: B, B_max, B_min, area, vol
real(wp_) :: J_phi, I_p, f_c, q, astra, jintrac
call tbl%init( &
title='flux-averages', id=56, active=is_active(params, 56), &
labels=[character(64) :: 'ρ_p', 'ρ_t', 'B_avg', &
'B_max', 'B_min', 'area', 'vol', 'I_pl', &
'J_φ_avg', 'f_c', 'J_astra/J_φ', 'J_jintrac/J_φ', &
'q'])
if (.not. tbl%active) return
if (params%equilibrium%iequil == EQ_VACUUM) return
! Δρ_p for the uniform ρ_p grid
drho_p = 1.0_wp_ / (n_grid - 1)
do i = 0, n_grid
rho_p = i * drho_p
rho_t = equil%pol2tor(rho_p)
call equil%flux_average(rho_p, &
B_avg=B, B_max=B_max, B_min=B_min, &
area=area, volume=vol, I_p=I_p, &
J_phi_avg=J_phi, f_c=f_c, q=q, &
ratio_astra_tor=astra, &
ratio_jintrac_tor=jintrac)
call tbl%append([rho_p, rho_t, B, B_max, B_min, area, &
vol, J_phi, f_c, astra, jintrac, q])
end do
end function flux_averages
function flux_surfaces(params, equil) result(tbl)
! Generates the flux surface contours
use gray_params, only : gray_parameters, EQ_VACUUM
use gray_equil, only : abstract_equil
! function arguments
type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(table) :: tbl
! local variables
integer, parameter :: n_surf = 10 , n_points = 101
real(wp_) :: rho_p, psi_n, drho_p
real(wp_) :: R(n_points), z(n_points)
real(wp_) :: R_hi, R_lo, z_hi, z_lo
integer :: i
call tbl%init( &
title='flux-surfaces', id=71, active=is_active(params, 71), &
labels=[character(64) :: 'i', 'ψ_n', 'R', 'z'])
if (.not. tbl%active) return
if (params%equilibrium%iequil == EQ_VACUUM) return
! Δρ_p for the uniform ρ_p grid
drho_p = 1.0_wp_ / (n_surf - 1)
! Guess for first contour reconstruction
R_hi = equil%axis(1)
R_lo = equil%axis(1)
z_hi = equil%axis(2) + (equil%z_boundary(2) - equil%axis(2))/10
z_lo = equil%axis(2) - (equil%axis(2) - equil%z_boundary(1))/10
do i = 0, n_surf - 1
rho_p = i * drho_p
psi_n = rho_p**2
call equil%flux_contour(psi_n, R, z, R_hi, z_hi, R_lo, z_lo)
call store_flux_surface(tbl, psi_n, R, z)
end do
! plasma boundary
psi_n = 0.9999_wp_
call equil%flux_contour(psi_n, R, z, R_hi, z_hi, R_lo, z_lo)
call store_flux_surface(tbl, psi_n, R, z)
end function flux_surfaces
function ec_resonance(params, equil, B_res) result(tbl)
! Generates the EC resonance surface table
@ -283,7 +371,7 @@ contains
end function inputs_maps
subroutine find_flux_surfaces(qvals, params, equil, tbl)
subroutine find_flux_surfaces(qvals, equil, tbl)
! Finds the ψ for a set of values of q and stores the
! associated surfaces to the flux surface table
@ -295,7 +383,6 @@ contains
! subroutine arguments
real(wp_), intent(in) :: qvals(:)
type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(table), intent(inout) :: tbl
@ -337,8 +424,7 @@ contains
R_lo = equil%axis(1)
z_lo = (equil%z_boundary(1) + equil%axis(2))/2
call equil%flux_contour(psi_n, params%misc%rwall, &
R_cont, z_cont, R_hi, z_hi, R_lo, z_lo)
call equil%flux_contour(psi_n, R_cont, z_cont, R_hi, z_hi, R_lo, z_lo)
call store_flux_surface(tbl, psi_n, R_cont, z_cont)
end block

376
src/magsurf_data.f90
View File

@ -1,376 +0,0 @@
module magsurf_data
use const_and_precisions, only : wp_
use splines, only : spline_simple, spline_2d
implicit none
type(spline_simple), save :: spline_volume, spline_R_J, spline_B_avg_min
type(spline_simple), save :: spline_B_max, spline_B_min, spline_area, spline_f_c
type(spline_simple), save :: spline_J_astra, spline_J_jintrac, spline_J_paral
type(spline_simple), save :: spline_dAdrho_t, spline_dVdrho_t
type(spline_2d), save :: cd_eff
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
subroutine compute_flux_averages(params, equil, tables)
! Computes the splines of the flux surface averages
use const_and_precisions, only : zero, one, comp_huge, pi, mu0inv
use types, only : table, wrap
use gray_params, only : gray_parameters, gray_tables
use gray_tables, only : store_flux_surface
use gray_equil, only : abstract_equil
! subroutine arguments
type(gray_parameters), intent(in) :: params
class(abstract_equil), intent(in) :: equil
type(gray_tables), intent(inout), optional :: tables
! 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) = equil%tor_curr(equil%axis(1), equil%axis(2))
B_max(1) = abs(equil%b_axis)
B_min(1) = abs(equil%b_axis)
B_avg(1) = abs(equil%b_axis)
R_J(1) = equil%axis(1)
f_c(1) = 1
dAdrho_t(1) = 0
dVdrho_t(1) = 0
ratio_paral_tor(1) = 1
ratio_astra_tor(1) = abs(equil%b_axis / equil%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 equil%pol_flux(equil%axis(1), equil%axis(2), ddpsidrr=a, ddpsidzz=b)
q(1) = equil%b_axis / sqrt(a*b)/equil%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 = equil%axis(1)
R_lo = equil%axis(1)
z_hi = equil%axis(2) + (equil%z_boundary(2) - equil%axis(2))/10
z_lo = equil%axis(2) - (equil%axis(2) - equil%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 equil%flux_contour(psi_n(i), params%misc%rwall, &
R, z, R_hi, z_hi, R_lo, z_lo)
if ((mod(i, 10) == 0 .or. i == nsurf) .and. present(tables)) then
call store_flux_surface(tables%flux_surfaces, psi_n(i), R, z)
end if
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 = equil%tor_curr(R(1), z(1))
call equil%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 = equil%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 equil%b_field(R(j+1), z(j+1), B_R=B_R, B_z=B_z, B_tot=B_tot1)
J_phi1 = equil%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*equil%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 = equil%pol2tor(rho_p(i))
dpsidrho_t = equil%phi_a * rho_t / (equil%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 spline_volume%init(rho_p, volume, nsurf)
call spline_B_avg_min%init(rho_p, B_avg/B_min, nsurf)
call spline_R_J%init(rho_p, R_J, nsurf)
call spline_B_max%init(rho_p, B_max, nsurf)
call spline_B_min%init(rho_p, B_min, nsurf)
call spline_J_astra%init(rho_p, ratio_astra_tor, nsurf)
call spline_J_jintrac%init(rho_p, ratio_jintrac_tor, nsurf)
call spline_J_paral%init(rho_p, ratio_paral_tor, nsurf)
call spline_area%init(rho_p, area, nsurf)
call spline_f_c%init(rho_p, f_c, nsurf)
call spline_dAdrho_t%init(rho_p, dAdrho_t, nsurf)
call spline_dVdrho_t%init(rho_p, dVdrho_t, nsurf)
call cd_eff%init(rho_p, lambda, reshape(H, [nsurf*nlam]), &
nsurf, nlam, range=[zero, one, zero, one])
if (.not. present(tables)) return
if (.not. tables%flux_averages%active) return
do i = 1, nsurf
call tables%flux_averages%append([ &
rho_p(i), equil%pol2tor(rho_p(i)), B_avg(i), B_max(i), &
B_min(i), area(i), volume(i), I_p(i), J_phi_avg(i), &
f_c(i), ratio_astra_tor(i), ratio_jintrac_tor(i), q(i)])
end do
end subroutine compute_flux_averages
subroutine fluxval(rho_p, area, vol, dervol, dadrhot, dvdrhot, &
rri, rbav, bmn, bmx, fc, ratja, ratjb, ratjpl)
use const_and_precisions, only : wp_
! subroutine arguments
real(wp_), intent(in) :: rho_p
real(wp_), intent(out), optional :: &
vol, area, rri, rbav, dervol, bmn, bmx, fc, &
ratja, ratjb, ratjpl, dadrhot, dvdrhot
if (present(area)) area = spline_area%eval(rho_p)
if (present(vol)) vol = spline_volume%eval(rho_p)
if (present(dervol)) dervol = spline_volume%deriv(rho_p)
if (present(dadrhot)) dAdrhot = spline_dAdrho_t%eval(rho_p)
if (present(dvdrhot)) dVdrhot = spline_dVdrho_t%eval(rho_p)
if (present(rri)) rri = spline_R_J%eval(rho_p)
if (present(rbav)) rbav = spline_B_avg_min%eval(rho_p)
if (present(bmn)) bmn = spline_B_min%eval(rho_p)
if (present(bmx)) bmx = spline_B_max%eval(rho_p)
if (present(fc)) fc = spline_f_c%eval(rho_p)
if (present(ratja)) ratja = spline_J_astra%eval(rho_p)
if (present(ratjb)) ratjb = spline_J_jintrac%eval(rho_p)
if (present(ratjpl)) ratjpl = spline_J_paral%eval(rho_p)
end subroutine fluxval
end module magsurf_data

4
src/main.f90
View File

@ -238,7 +238,6 @@ contains
! (of different beams) and recomputes the summary table
use gray_params, only : gray_parameters
use gray_equil, only : abstract_equil
use magsurf_data, only : compute_flux_averages
use pec, only : pec_init, postproc_profiles, dealloc_pec, &
rhot_tab
@ -263,9 +262,6 @@ contains
integer :: list, prof, summ
integer, allocatable :: beams(:)
! Initialise the magsurf_data module
call compute_flux_averages(params, equil)
! Initialise the output profiles
call pec_init(params%output, equil)

13
src/pec.f90
View File

@ -13,7 +13,6 @@ contains
subroutine pec_init(params, equil, rt_in)
use gray_params, only : output_parameters, RHO_POL
use gray_equil, only : abstract_equil
use magsurf_data, only : fluxval
! subroutine arguments
type(output_parameters), intent(in) :: params
@ -72,13 +71,14 @@ contains
! psi grid at mid points, size n+1, for use in pec_tab
rtabpsi1(it) = rhop1**2
call fluxval(rhop1,area=areai1,vol=voli1)
call equil%flux_average(rhop1, area=areai1, volume=voli1)
dvol(it) = abs(voli1 - voli0)
darea(it) = abs(areai1 - areai0)
voli0 = voli1
areai0 = areai1
call fluxval(rhop_tab(it),ratja=ratjai,ratjb=ratjbi,ratjpl=ratjpli)
call equil%flux_average(rhop_tab(it), ratio_astra_tor=ratjai, &
ratio_jintrac_tor=ratjbi, ratio_paral_tor=ratjpli)
ratjav(it) = ratjai
ratjbv(it) = ratjbi
ratjplv(it) = ratjpli
@ -251,7 +251,6 @@ contains
! Radial average values over power and current density profile
use const_and_precisions, only : pi, comp_eps
use gray_equil, only : abstract_equil
use magsurf_data, only : fluxval
class(abstract_equil), intent(in) :: equil
real(wp_), intent(in) :: pabs,currt
@ -294,7 +293,7 @@ contains
rhopjava = equil%tor2pol(rhotjava)
if (pabs > 0) then
call fluxval(rhoppav,dvdrhot=dvdrhotav)
call equil%flux_average(rhoppav, dVdrho_t=dvdrhotav)
dpdvp = pabs*2.0_wp_/(sqrt(pi)*drhotpav*dvdrhotav)
call profwidth(rhot_tab,dpdv,rhotp,dpdvmx,drhotp)
else
@ -305,8 +304,8 @@ contains
end if
if (sccsa > 0) then
call fluxval(rhopjava,dadrhot=dadrhotava,ratja=ratjamx,ratjb=ratjbmx, &
ratjpl=ratjplmx)
call equil%flux_average(rhopjava, dAdrho_t=dadrhotava, ratio_astra_tor=ratjamx, &
ratio_jintrac_tor=ratjbmx, ratio_paral_tor=ratjplmx)
ajphip = currt*2.0_wp_/(sqrt(pi)*drhotjava*dadrhotava)
call profwidth(rhot_tab,ajphiv,rhotjfi,ajmxfi,drhotjfi)
else