From 98dda6d6facbda193def4f3be549bd5b63be3f47 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Mon, 16 Oct 2023 11:40:10 +0200 Subject: [PATCH] =?UTF-8?q?src/equilibrium.f90:=20improve=20handling=20of?= =?UTF-8?q?=20the=20=CF=88=20normalisation?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Hide the implementation of the re-normalisation of the ψ(R,z) spline by adjusting the spline coefficients, instead of shifting and rescaling after each evaluation. - Correct the value of `psia` = ψ(X point) - ψ(O point) after the ψ(R,z) spline has been re-normalised. This fixes another instance of decoupling between the values of X and ∇X that introduce a systematic error in the numerical integration of the raytracing equations. Here the issue is caused by the different normalisations used, specifically X=X(ψ_n') and ∇X=∇X(ψ_n), where ψ_n' is the re-normalised spline of the normalised flux and ψ_n the spline of the normalised flux. See d0a5a9f for more details on problems resulting from this error. - Change `equian` and `equinum_psi` to return the normalised values for both the flux and its derivatives, to avoid confusions. Callers that needed the unnormalised derivatives now multiply explicitly by `psia`. --- src/equilibrium.f90 | 80 ++++++++++++++++++++++++++------------------ src/gray_core.f90 | 34 +++++++++---------- src/magsurf_data.f90 | 11 +++--- 3 files changed, 69 insertions(+), 56 deletions(-) diff --git a/src/equilibrium.f90 b/src/equilibrium.f90 index 65195f9..b723b86 100644 --- a/src/equilibrium.f90 +++ b/src/equilibrium.f90 @@ -29,7 +29,7 @@ module equilibrium ! Splines type(spline_1d), save :: fpol_spline ! Poloidal current function F(ψ) - type(spline_2d), save :: psi_spline ! Poloidal flux ψ(R,z) + type(spline_2d), save :: psi_spline ! Normalised poloidal flux ψ(R,z) type(spline_simple), save :: q_spline ! Safey factor q(ψ) type(linear_1d), save :: rhop_spline, rhot_spline ! Normalised radii ρ_p(ρ_t), ρ_t(ρ_p) @@ -39,8 +39,6 @@ module equilibrium ! More parameters real(wp_), save :: fpolas ! Poloidal current at the edge (r=a), F_a real(wp_), save :: psia ! Poloidal flux at the edge (r=a), ψ_a - real(wp_), save :: psiant ! Normalised poloidal flux at the true edge - real(wp_), save :: psinop ! Normalised poloidal flux at the true O point real(wp_), save :: phitedge ! Toroidal flux at the edge real(wp_), save :: btaxis ! B₀: B_φ = B₀R₀/R real(wp_), save :: rmaxis, zmaxis ! R,z of the magnetic axis (O point) @@ -65,8 +63,7 @@ module equilibrium public kspl, psi_spline, q_spline, points_tgo, model ! Members exposed to gray_core and more - public psia, psiant, psinop - public phitedge + public psia, phitedge public btaxis, rmaxis, zmaxis, sgnbphi public btrcen, rcen public rmnm, rmxm, zmnm, zmxm @@ -416,6 +413,7 @@ contains integer :: nr, nz, nrest, nzest, npsest, nrz, npsi, nbnd, ibinf, ibsup real(wp_) :: tension, rax0, zax0, psinoptmp, psinxptmp real(wp_) :: rbmin, rbmax, rbinf, rbsup, r1, z1 + real(wp_) :: psiant, psinop real(wp_), dimension(size(data%psinr)) :: rhotn real(wp_), dimension(:), allocatable :: rv1d, zv1d, fvpsi, wf integer :: ixploc, info, i, j, ij @@ -555,12 +553,13 @@ contains fpolas = fpol_spline%eval(data%psinr(npsi)) sgnbphi = sign(one ,fpolas) - ! Re-normalize ψ after spline computation + ! Re-normalize ψ_n after the spline computation + ! Note: this ensures 0 ≤ ψ_n'(R,z) < 1 inside the plasma - ! Start with un-corrected psi + ! Start with un-corrected ψ_n psia = data%psia - psinop = 0 - psiant = 1 + psinop = 0 ! ψ_n(O point) + psiant = 1 ! ψ_n(X point) - ψ_n(O point) ! Use provided boundary to set an initial guess ! for the search of O/X points @@ -643,6 +642,22 @@ contains call log_info(msg, mod='equilibrium', proc='set_equil_spline') end if + ! Adjust all the B-spline coefficients + ! Note: since ψ_n(R,z) = Σ_ij c_ij B_i(R)B_j(z), to correct ψ_n + ! it's enough to correct (shift and scale) the c_ij + psi_spline%coeffs = (psi_spline%coeffs - psinop) / psiant + + ! Same for the derivatives + psi_spline%partial(1,0)%ptr = psi_spline%partial(1,0)%ptr / psiant + psi_spline%partial(0,1)%ptr = psi_spline%partial(0,1)%ptr / psiant + psi_spline%partial(2,0)%ptr = psi_spline%partial(2,0)%ptr / psiant + psi_spline%partial(0,2)%ptr = psi_spline%partial(0,2)%ptr / psiant + psi_spline%partial(1,1)%ptr = psi_spline%partial(1,1)%ptr / psiant + + ! Do the opposite scaling to preserve un-normalised values + ! Note: this is only used for the poloidal magnetic field + psia = psia * psiant + ! Save Bt value on axis (required in flux_average and used in Jcd def) ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def) @@ -804,12 +819,12 @@ contains if (R <= rmxm .and. R >= rmnm .and. & z <= zmxm .and. z >= zmnm) then - if (present(psi)) psi = (psi_spline%eval(R, z) - psinop) / psiant - if (present(dpsidr)) dpsidr = psi_spline%deriv(R, z, 1, 0) * psia - if (present(dpsidz)) dpsidz = psi_spline%deriv(R, z, 0, 1) * psia - if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0) * psia - if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2) * psia - if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1) * psia + 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 @@ -863,13 +878,13 @@ contains ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a² if (present(psin)) psin = rho_p**2 - ! ∂ψ/∂R, ∂ψ/∂z, note that ψ = ψ_n⋅ψ_a - if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2 * psia - if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2 * psia + ! ∂ψ/∂R, ∂ψ/∂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 * psia - if (present(ddpsidzz)) ddpsidzz = 2/model%a**2 * psia + if (present(ddpsidrr)) ddpsidrr = 2/model%a**2 + if (present(ddpsidzz)) ddpsidzz = 2/model%a**2 if (present(ddpsidrz)) ddpsidrz = 0 ! F(ψ) = B₀⋅R₀, a constant @@ -1120,8 +1135,7 @@ contains call log_error(msg, mod='equilibrium', proc='psi_raxis') end if - val=psin*psiant+psinop - call sproota(val, psi_spline%knots_x, psi_spline%nknots_x, & + call sproota(psin, psi_spline%knots_x, psi_spline%nknots_x, & czc, zeroc, mest, m, ier) r1=zeroc(1) r2=zeroc(2) @@ -1186,15 +1200,15 @@ contains select case(iflag) case(1) call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz) - fvec(1) = dpsidr/psia - fvec(2) = dpsidz/psia + fvec(1) = dpsidr + fvec(2) = dpsidz case(2) call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, & ddpsidrz=ddpsidrz) - fjac(1,1) = ddpsidrr/psia - fjac(1,2) = ddpsidrz/psia - fjac(2,1) = ddpsidrz/psia - fjac(2,2) = ddpsidzz/psia + fjac(1,1) = ddpsidrr + fjac(1,2) = ddpsidrz + fjac(2,1) = ddpsidrz + fjac(2,2) = ddpsidzz case default write (msg, '("invalid iflag: ",g0)') call log_error(msg, mod='equilibrium', proc='fcnox') @@ -1260,14 +1274,14 @@ contains case(1) call equinum_psi(x(1),x(2),psinv,dpsidr) fvec(1) = psinv-f0(1) - fvec(2) = dpsidr/psia-f0(2) + fvec(2) = dpsidr-f0(2) case(2) call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, & ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz) - fjac(1,1) = dpsidr/psia - fjac(1,2) = dpsidz/psia - fjac(2,1) = ddpsidrr/psia - fjac(2,2) = ddpsidrz/psia + fjac(1,1) = dpsidr + fjac(1,2) = dpsidz + fjac(2,1) = ddpsidrr + fjac(2,2) = ddpsidrz case default write (msg, '("invalid iflag: ",g0)') call log_error(msg, mod='equilibrium', proc='fcntgo') diff --git a/src/gray_core.f90 b/src/gray_core.f90 index e902540..e5896da 100644 --- a/src/gray_core.f90 +++ b/src/gray_core.f90 @@ -1258,7 +1258,7 @@ contains subroutine plas_deriv(xv, bres, xgcn, dens, btot, bv, derbv, & xg, yg, derxg, deryg, psinv_opt) - use const_and_precisions, only : zero + use const_and_precisions, only : zero, cm use gray_params, only : iequil use equilibrium, only : psia, equinum_fpol, equinum_psi, equian, sgnbphi use coreprofiles, only : density @@ -1281,7 +1281,7 @@ contains real(wp_), dimension(3) :: dbtot,bvc real(wp_), dimension(3,3) :: dbvcdc,dbvdc,dbv real(wp_) :: brr,bphi,bzz,dxgdpsi - real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddenspsin + real(wp_) :: dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz,fpolv,dfpv,ddensdpsi xg = zero yg = 99._wp_ @@ -1342,14 +1342,14 @@ contains end if ! compute xg and derivative - call density(psinv,dens,ddenspsin) + call density(psinv, dens, ddensdpsi) xg = xgcn*dens - dxgdpsi = xgcn*ddenspsin/psia + dxgdpsi = xgcn*ddensdpsi ! B = f(psi)/R e_phi+ grad psi x e_phi/R - bphi = fpolv/rrm - brr =-dpsidz/rrm - bzz = dpsidr/rrm + bphi = fpolv/rrm + brr = -dpsidz*psia/rrm + bzz = +dpsidr*psia/rrm ! bvc(i) = B_i in cylindrical coordinates bvc = [brr, bphi, bzz] @@ -1360,12 +1360,12 @@ contains bv(3)=bvc(3) ! dbvcdc(iv,jv) = d Bcil(iv) / dxvcil(jv) - dbvcdc(1,1) = -ddpsidrz/rrm - brr/rrm - dbvcdc(2,1) = dfpv*dpsidr/rrm - bphi/rrm - dbvcdc(3,1) = ddpsidrr/rrm - bzz/rrm - dbvcdc(1,3) = -ddpsidzz/rrm - dbvcdc(2,3) = dfpv*dpsidz/rrm - dbvcdc(3,3) = ddpsidrz/rrm + dbvcdc(1,1) = -ddpsidrz*psia/rrm - brr/rrm + dbvcdc(2,1) = dfpv*dpsidr*psia/rrm - bphi/rrm + dbvcdc(3,1) = ddpsidrr*psia/rrm - bzz/rrm + dbvcdc(1,3) = -ddpsidzz*psia/rrm + dbvcdc(2,3) = dfpv*dpsidz*psia/rrm + dbvcdc(3,3) = ddpsidrz*psia/rrm ! dbvdc(iv,jv) = d Bcart(iv) / dxvcil(jv) dbvdc(1,1) = dbvcdc(1,1)*csphi - dbvcdc(2,1)*snphi @@ -1403,12 +1403,12 @@ contains deryg = 1.0e-2_wp_*dbtot/Bres bv = bv/btot do jv=1,3 - derbv(:,jv) = 1.0e-2_wp_*(dbv(:,jv) - bv(:)*dbtot(jv))/btot + derbv(:,jv) = cm * (dbv(:,jv) - bv(:)*dbtot(jv))/btot end do - derxg(1) = 1.0e-2_wp_*drrdx*dpsidr*dxgdpsi - derxg(2) = 1.0e-2_wp_*drrdy*dpsidr*dxgdpsi - derxg(3) = 1.0e-2_wp_*dpsidz *dxgdpsi + derxg(1) = cm * drrdx*dpsidr*dxgdpsi + derxg(2) = cm * drrdy*dpsidr*dxgdpsi + derxg(3) = cm * dpsidz *dxgdpsi end subroutine plas_deriv diff --git a/src/magsurf_data.f90 b/src/magsurf_data.f90 index d7fced9..f1287d0 100644 --- a/src/magsurf_data.f90 +++ b/src/magsurf_data.f90 @@ -103,7 +103,7 @@ contains 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 + equian,equinum_psi,bfield,frhotor,fq,tor_curr,psia use dierckx, only : regrid,coeff_parder implicit none @@ -175,8 +175,8 @@ contains else call equinum_psi(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz) end if - qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz) - ajphiav=-ccj*(ddpsidrr+ddpsidzz)/rmaxis + qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz)/psia + ajphiav=-ccj*(ddpsidrr+ddpsidzz)*psia/rmaxis psicon(1)=0.0_wp_ rcon(:,1)=rmaxis @@ -448,7 +448,7 @@ contains use logger, only : log_warning use dierckx, only : profil, sproota use equilibrium, only : model, frhotor, psi_spline, & - kspl, psiant, psinop, points_tgo + kspl, points_tgo use limiter, only : rwallm implicit none @@ -510,8 +510,7 @@ contains call log_warning(msg, mod='magsurf_data', proc='contours_psi') end block end if - val=h*psiant+psinop - call sproota(val, psi_spline%knots_x, psi_spline%nknots_x, & + call sproota(h, psi_spline%knots_x, psi_spline%nknots_x, & czc, zeroc, mest, m, ier) if (zeroc(1).gt.rwallm) then rcn(ic)=zeroc(1)