From 31a7b95fe1616eb5d3350f723408f4cc3aded4f2 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Sun, 22 Oct 2023 10:32:27 +0200 Subject: [PATCH] combine equian and equinum_* subroutines MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This change combines `equian` and `equinum_psi` into a new `pol_flux` subroutine that computes ψ(R,z) and derivatives for either numerical or analytical equilibria. Similarly, `equian` (the fpol and dfpol outputs) and `equinum_fpol` are combined intro `pol_curr` that computes F(ψ) and its derivative. Callers of these subroutines do not select a specific version based on the value of `iequil` anymore. --- src/coreprofiles.f90 | 4 +- src/equilibrium.f90 | 271 ++++++++++++++++++++++--------------------- src/gray_core.f90 | 28 ++--- src/magsurf_data.f90 | 13 +-- src/multipass.f90 | 4 +- 5 files changed, 157 insertions(+), 163 deletions(-) diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index b7926f6..633c982 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -451,7 +451,7 @@ contains ! invalid as they are given assuming a vacuum (N=1) if (present(launch_pos)) then block - use equilibrium, only : equinum_psi + use equilibrium, only : pol_flux use const_and_precisions, only : cm real(wp_) :: R, Z, psi @@ -461,7 +461,7 @@ contains ! Get the poloidal flux at the launcher ! Note: this returns -1 when the data is not available - call equinum_psi(R, z, psi) + call pol_flux(R, z, psi) if (psi > tail%start .and. psi < tail%end) then ! Fall back to the midpoint of ψ₀ and the launcher ψ diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index b723b86..fe2bdff 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -52,9 +52,8 @@ module equilibrium private public read_eqdsk, read_equil_an ! Reading data files public scale_equil, change_cocos ! Transforming data - public equian ! Analytical model - public equinum_psi, equinum_fpol ! Interpolated 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 @@ -661,7 +660,7 @@ contains ! 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 equinum_fpol(zero, btaxis) + call pol_curr(zero, btaxis) btaxis = btaxis/rmaxis btrcen = fpolas/data%rvac rcen = data%rvac @@ -803,94 +802,94 @@ contains end subroutine set_rho_spline - subroutine equinum_psi(R, z, psi, dpsidr, dpsidz, & - ddpsidrr, ddpsidzz, ddpsidrz) - ! Computes the normalised poloidal flux ψ (and its derivatives) - ! as a function of (R, z) using the numerical equilibrium + 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 implicit none - ! subroutine arguments - real(wp_), intent(in) :: R, z - real(wp_), intent(out), optional :: psi, dpsidr, dpsidz - real(wp_), intent(out), optional :: ddpsidrr, ddpsidzz, ddpsidrz - - ! Note: here lengths are measured in meters - if (R <= rmxm .and. R >= rmnm .and. & - z <= zmxm .and. z >= zmnm) then - - if (present(psi)) psi = 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 - if (present(psi)) psi = -1 - if (present(dpsidr)) dpsidr = 0 - if (present(dpsidz)) dpsidz = 0 - if (present(ddpsidrr)) ddpsidrr = 0 - if (present(ddpsidzz)) ddpsidzz = 0 - if (present(ddpsidrz)) ddpsidrz = 0 - end if - end subroutine equinum_psi - - - subroutine equinum_fpol(psi, fpol, dfpol) - ! Computes the poloidal current function F(ψ) - ! (and its derivative) using the numerical equilibrium - - implicit none - - ! subroutine arguments - real(wp_), intent(in) :: psi - real(wp_), intent(out) :: fpol - real(wp_), intent(out), optional :: dfpol - - if(psi <= 1 .and. psi >= 0) then - fpol = fpol_spline%eval(psi) - if (present(dfpol)) dfpol = fpol_spline%deriv(psi) / psia - else - fpol = fpolas - if (present(dfpol)) dfpol = 0 - end if - end subroutine equinum_fpol - - - subroutine equian(R, z, psin, fpolv, dfpv, dpsidr, dpsidz, & - ddpsidrr, ddpsidzz, ddpsidrz) - ! Computes all MHD equilibrium parameters from a simple analytical model - implicit none - ! subroutine arguments real(wp_), intent(in) :: R, z - real(wp_), intent(out), optional :: psin, fpolv, dfpv, dpsidr, dpsidz, & + real(wp_), intent(out), optional :: psi_n, dpsidr, dpsidz, & ddpsidrr, ddpsidzz, ddpsidrz - ! variables + ! local variables real(wp_) :: r_p ! poloidal radius real(wp_) :: rho_p ! poloidal radius normalised - ! ρ_p is just the geometrical poloidal radius divided by a - r_p = hypot(R - model%R0, z - model%z0) - rho_p = r_p / model%a + if (iequil < 2) then + ! Analytical model + ! ρ_p is just the geometrical poloidal radius divided by a + r_p = hypot(R - model%R0, z - model%z0) + rho_p = r_p / model%a - ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a² - if (present(psin)) psin = rho_p**2 + ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a² + if (present(psi_n)) psi_n = rho_p**2 - ! ∂ψ/∂R, ∂ψ/∂z - if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2 - if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2 + ! ∂ψ_n/∂R, ∂ψ_n/∂z + if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2 + if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2 - ! ∂²ψ/∂R², ∂²ψ/∂z², ∂²ψ/∂R∂z - if (present(ddpsidrr)) ddpsidrr = 2/model%a**2 - if (present(ddpsidzz)) ddpsidzz = 2/model%a**2 - if (present(ddpsidrz)) ddpsidrz = 0 + ! ∂²ψ_n/∂R², ∂²ψ_n/∂z², ∂²ψ_n/∂R∂z + if (present(ddpsidrr)) ddpsidrr = 2/model%a**2 + if (present(ddpsidzz)) ddpsidzz = 2/model%a**2 + if (present(ddpsidrz)) ddpsidrz = 0 + else + ! Numerical data + if (R <= rmxm .and. R >= rmnm .and. & + z <= zmxm .and. z >= zmnm) 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 - ! F(ψ) = B₀⋅R₀, a constant - if (present(fpolv)) fpolv = model%B0 * model%R0 - if (present(dfpv)) dfpv = 0 - end subroutine equian + + subroutine pol_curr(psi_n, fpol, dfpol) + ! Computes the poloidal current function F(ψ_n) + ! and (optionally) its derivative dF/dψ 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) / psia + else + fpol = fpolas + if (present(dfpol)) dfpol = 0 + end if + end if + end subroutine pol_curr function frhotor(rho_p) @@ -1017,83 +1016,93 @@ contains end function fq - subroutine bfield(rpsim, zpsim, bphi, br, bz) + 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) :: rpsim, zpsim - real(wp_), intent(out), optional :: bphi,br,bz + real(wp_), intent(in) :: R, z + real(wp_), intent(out), optional :: B_R, B_z, B_phi ! local variables - real(wp_) :: psin,fpol + real(wp_) :: psi_n, fpol, dpsidr, dpsidz - if (iequil < 2) then - ! Analytical model - call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br) - if (present(bphi)) bphi=bphi/rpsim - else - ! Numerical data - call equinum_psi(rpsim,zpsim,psi=bphi,dpsidr=bz,dpsidz=br) - if (present(bphi)) then - psin=bphi - call equinum_fpol(psin,fpol) - bphi=fpol/rpsim - end if - end if + call pol_flux(R, z, psi_n, dpsidr, dpsidz) + call pol_curr(psi_n, fpol) - if (present(br)) br=-br/rpsim - if (present(bz)) bz= bz/rpsim + ! 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(rpsim,zpsim) result(jphi) - ! Computes the toroidal current Jφ as a function of (R, z) - use const_and_precisions, only : ccj=>mu0inv - use gray_params, only : iequil + 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) :: rpsim, zpsim - real(wp_) :: jphi + real(wp_), intent(in) :: R, z + real(wp_) :: J_phi ! local variables - real(wp_) :: bzz,dbvcdc13,dbvcdc31 - real(wp_) :: dpsidr,ddpsidrr,ddpsidzz + real(wp_) :: dB_Rdz, dB_zdR ! derivatives of B_R, B_z + real(wp_) :: dpsidr, ddpsidrr, ddpsidzz ! derivatives of ψ_n - if(iequil < 2) then - ! Analytical model - call equian(rpsim,zpsim,dpsidr=dpsidr, & - ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) - else - ! Numerical data - call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, & - ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) - end if - bzz= dpsidr/rpsim - dbvcdc13=-ddpsidzz/rpsim - dbvcdc31= ddpsidrr/rpsim-bzz/rpsim - jphi=ccj*(dbvcdc13-dbvcdc31) + 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(psin) result(jphi) - ! Computes the toroidal current Jφ as a function of ψ + 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) :: psin - real(wp_) :: jphi + real(wp_), intent(in) :: psi_n + real(wp_) :: J_phi ! local variables - real(wp_) :: r1, r2 - call psi_raxis(psin, r1, r2) - jphi = tor_curr(r2, zmaxis) + real(wp_) :: R1, R2 + call psi_raxis(psi_n, R1, R2) + J_phi = tor_curr(R2, zmaxis) end function tor_curr_psi @@ -1178,7 +1187,7 @@ contains end if rf=xvec(1) zf=xvec(2) - call equinum_psi(rf,zf,psinvf) + call pol_flux(rf, zf, psinvf) end subroutine points_ox @@ -1199,12 +1208,12 @@ contains select case(iflag) case(1) - call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) + call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz) fvec(1) = dpsidr fvec(2) = dpsidz case(2) - call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & - ddpsidrz=ddpsidrz) + 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 @@ -1272,12 +1281,12 @@ contains select case(iflag) case(1) - call equinum_psi(x(1),x(2),psinv,dpsidr) + call pol_flux(x(1), x(2), psinv, dpsidr) fvec(1) = psinv-f0(1) fvec(2) = dpsidr-f0(2) case(2) - call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & - ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) + 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 diff --git a/src/gray_core.f90 b/src/gray_core.f90 index e5896da..61b1f5c 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -1260,7 +1260,7 @@ contains xg, yg, derxg, deryg, psinv_opt) use const_and_precisions, only : zero, cm use gray_params, only : iequil - use equilibrium, only : psia, equinum_fpol, equinum_psi, equian, sgnbphi + use equilibrium, only : psia, pol_flux, pol_curr, sgnbphi use coreprofiles, only : density implicit none @@ -1322,13 +1322,8 @@ contains zzm = 1.0e-2_wp_*zz rrm = 1.0e-2_wp_*rr - if(iequil==1) then - call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, & - ddpsidrr,ddpsidzz,ddpsidrz) - else - call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz) - call equinum_fpol(psinv,fpolv,dfpv) - end if + call pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz) + call pol_curr(psinv, fpolv, dfpv) ! copy optional output if (present(psinv_opt)) psinv_opt = psinv @@ -1886,7 +1881,7 @@ contains rm=sqrt(xmv(1)**2+xmv(2)**2) csphi=xmv(1)/rm snphi=xmv(2)/rm - call bfield(rm,xmv(3),bphi,br,bz) + call bfield(rm,xmv(3),br,bz,bphi) ! bv(i) = B_i in cartesian coordinates bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi @@ -2266,7 +2261,7 @@ bb: do zzk = zv(k) do j = 1, q_spline%ndata rrj = rv(j) - call bfield(rrj, zzk, bbphi, bbr, bbz) + call bfield(rrj, zzk, bbr, bbz, bbphi) btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2) if(btotal(j,k) >= btmx) btmx = btotal(j,k) if(btotal(j,k) <= btmn) btmn = btotal(j,k) @@ -2295,9 +2290,8 @@ bb: do subroutine print_maps(bres, xgcn, r0, anpl0) ! Prints several input quantities on a regular grid (unit 72) - use gray_params, only : iequil - use equilibrium, only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, & - equinum_fpol, q_spline + use equilibrium, only : rmnm, rmxm, zmnm, zmxm, & + pol_flux, pol_curr, q_spline use coreprofiles, only : density, temp use units, only : umaps, unit_active @@ -2327,12 +2321,8 @@ bb: do anpl = anpl0 * r0/rj do k = 1, q_spline%ndata zk = z(k) - if (iequil < 2) then - call equian(rj, zk, psin=psin, fpolv=bphi, dpsidr=bz, dpsidz=br) - else - call equinum_psi(rj, zk, psi=psin, dpsidr=bz, dpsidz=br) - call equinum_fpol(psin, fpol=bphi) - end if + call pol_flux(rj, zk, psin, dpsidr=bz, dpsidz=br) + call pol_curr(psin, fpol=bphi) br = -br/rj bphi = bphi/rj bz = bz/rj diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index f1287d0..f4584db 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -101,9 +101,8 @@ contains subroutine flux_average use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv - use gray_params, only : iequil use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, & - equian,equinum_psi,bfield,frhotor,fq,tor_curr,psia + bfield,frhotor,fq,tor_curr,psia,pol_flux use dierckx, only : regrid,coeff_parder implicit none @@ -170,11 +169,7 @@ contains ratio_cdbtor=1.0_wp_ ratio_pltor=1.0_wp_ fc=1.0_wp_ - if(iequil < 2) then - call equian(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) - else - call equinum_psi(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) - end if + call pol_flux(rmaxis, zmaxis, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz) qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz)/psia ajphiav=-ccj*(ddpsidrr+ddpsidzz)*psia/rmaxis @@ -230,7 +225,7 @@ contains bmmn=1.0e+30_wp_ ajphi0 = tor_curr(rctemp(1),zctemp(1)) - call bfield(rctemp(1),zctemp(1),bphi,br=brr,bz=bzz) + call bfield(rctemp(1), zctemp(1), brr, bzz, bphi) fpolv=bphi*rctemp(1) btot0=sqrt(bphi**2+brr**2+bzz**2) bpoloid0=sqrt(brr**2+bzz**2) @@ -256,7 +251,7 @@ contains ! compute line integrals on the contour psi=psinjp=height^2 rpsim=rctemp(inc1) zpsim=zctemp(inc1) - call bfield(rpsim,zpsim,br=brr,bz=bzz) + call bfield(rpsim, zpsim, brr, bzz) ajphi = tor_curr(rpsim,zpsim) bphi=fpolv/rpsim btot=sqrt(bphi**2+brr**2+bzz**2) diff --git a/src/multipass.f90 b/src/multipass.f90 index 7df3349..ea80603 100644 --- a/src/multipass.f90 +++ b/src/multipass.f90 @@ -37,7 +37,7 @@ contains rm=sqrt(xmv(1)**2+xmv(2)**2) csphi=xmv(1)/rm snphi=xmv(2)/rm - call bfield(rm,xmv(3),bphi,br,bz) + call bfield(rm,xmv(3),br,bz,bphi) bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi bv(3)=bz @@ -82,7 +82,7 @@ contains rm=sqrt(xmv(1)**2+xmv(2)**2) csphi=xmv(1)/rm snphi=xmv(2)/rm - call bfield(rm,xmv(3),bphi,br,bz) + call bfield(rm,xmv(3),br,bz,bphi) bv(1)=br*csphi-bphi*snphi bv(2)=br*snphi+bphi*csphi bv(3)=bz