diff --git a/src/eccd.f90 b/src/eccd.f90 index dd94257..d18445d 100644 --- a/src/eccd.f90 +++ b/src/eccd.f90 @@ -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 diff --git a/src/gray_core.f90 b/src/gray_core.f90 index 0d0496c..98a6e91 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -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 diff --git a/src/gray_equil.f90 b/src/gray_equil.f90 index 2ad4741..0f87f8a 100644 --- a/src/gray_equil.f90 +++ b/src/gray_equil.f90 @@ -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) = ¾ ∫₀¹ 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 diff --git a/src/gray_jetto1beam.f90 b/src/gray_jetto1beam.f90 index d69dae3..35b0c6e 100644 --- a/src/gray_jetto1beam.f90 +++ b/src/gray_jetto1beam.f90 @@ -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 diff --git a/src/gray_tables.f90 b/src/gray_tables.f90 index ec77af8..52eac3e 100644 --- a/src/gray_tables.f90 +++ b/src/gray_tables.f90 @@ -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 diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 deleted file mode 100644 index 2256b31..0000000 --- a/src/magsurf_data.f90 +++ /dev/null @@ -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) = ¾ ∫₀¹ 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 diff --git a/src/main.f90 b/src/main.f90 index 358d1c6..8c9a484 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -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) diff --git a/src/pec.f90 b/src/pec.f90 index b5f2080..6b035b0 100644 --- a/src/pec.f90 +++ b/src/pec.f90 @@ -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