From 272e337b214bfbd76d067869e92c58ce14242487 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 19 Oct 2026 23:37:25 +0200 Subject: [PATCH 1/5] GPEC - NEW FEATURE - Write the xi_perp-only Lagrangian field strength b_lag_perp with pmodb pmodb_flag now also writes b_lag_perp and b_lag_perp_fun to the profile NetCDF output: the Lagrangian field strength of the perpendicular displacement alone, b_lag_perp = lagb + (xi_perp^zeta/q) d_theta B = eulb + xi_perp.grad B, the gauge-invariant input that NTV consumers in perturbed Boozer or Hamada coordinates need. It is the quantity gpdiag_pmodb already prints as llagb; lagb itself is the xi^zeta = 0 gauge and is unchanged. The ASCII gpec_pmodb file keeps its format. --- gpec/gpout.f | 43 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 36 insertions(+), 7 deletions(-) diff --git a/gpec/gpout.f b/gpec/gpout.f index a02897d0..021cc0ee 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -2934,22 +2934,22 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) REAL(r8) :: jac,psi INTEGER :: p_id,t_id,i_id,m_id,r_id,z_id,b_id,bme_id,be_id, - $ bml_id,bl_id,xm_id,x_id,km_id,k_id,rzstat + $ bml_id,bl_id,xm_id,x_id,km_id,k_id,rzstat,bmp_id,bp_id REAL(r8), DIMENSION(0:mthsurf) :: dphi REAL(r8), DIMENSION(:,:), ALLOCATABLE :: rs,zs,equilbfun COMPLEX(r8), DIMENSION(mpert) :: eulbpar_mn,lagbpar_mn, - $ divxprp_mn,curv_mn + $ divxprp_mn,curv_mn,lagbprp_mn COMPLEX(r8), DIMENSION(0:mthsurf) :: xsp_fun,xms_fun, $ bvt_fun,bvz_fun,xmz_fun,xvt_fun,xvz_fun,xmt_fun,xsp1_fun COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparmns, $ lagbparmns,divxprpmns,curvmns, $ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparmout, - $ lagbparmout,divxprpmout,curvmout + $ lagbparmout,divxprpmout,curvmout,lagbprpmout COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparfun,lagbparfun, $ divxprpfun,curvfun,eulbparfout,lagbparfout,divxprpfout, - $ curvfout + $ curvfout,lagbprpfun,lagbprpfout REAL(r8), DIMENSION(:), ALLOCATABLE :: psis,ches,chex COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE ::chea,chelagbmns @@ -2965,11 +2965,13 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ xmp1mns(mstep,mpert),xspmns(mstep,mpert),xmsmns(mstep,mpert), $ xmtmns(mstep,mpert),xmzmns(mstep,mpert)) ALLOCATE(eulbparmout(mstep,lmpert),lagbparmout(mstep,lmpert), - $ divxprpmout(mstep,lmpert),curvmout(mstep,lmpert)) + $ divxprpmout(mstep,lmpert),curvmout(mstep,lmpert), + $ lagbprpmout(mstep,lmpert)) ALLOCATE(eulbparfun(mstep,0:mthsurf),lagbparfun(mstep,0:mthsurf), $ divxprpfun(mstep,0:mthsurf),curvfun(mstep,0:mthsurf), $ eulbparfout(mstep,0:mthsurf),lagbparfout(mstep,0:mthsurf), - $ divxprpfout(mstep,0:mthsurf),curvfout(mstep,0:mthsurf)) + $ divxprpfout(mstep,0:mthsurf),curvfout(mstep,0:mthsurf), + $ lagbprpfun(mstep,0:mthsurf),lagbprpfout(mstep,0:mthsurf)) c----------------------------------------------------------------------- c compute necessary components. c----------------------------------------------------------------------- @@ -3042,6 +3044,12 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ eulbparfun(istep,itheta)+ $ xsp_fun(itheta)*eqfun%fx(1)+ $ xms_fun(itheta)/(chi1*sq%f(4))*eqfun%fy(1) +c Lagrangian dB of xi_perp alone: lagb plus the gauge term +c (xi_perp^zeta/q) d_theta B (docs/pmodb_gauge.rst). + lagbprpfun(istep,itheta)=lagbparfun(istep,itheta)+ + $ eqfun%fy(1)*(xmz_fun(itheta)/(jac*sq%f(4))- + $ (chi1/(jac*eqfun%f(1)))**2* + $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta))) divxprpfun(istep,itheta)=xsp1_fun(itheta)+ $ (jac1/jac)*xsp_fun(itheta)+ $ cspl%f1(1)/jac-(twopi*ifac*nn)*cspl%f(2)/jac @@ -3064,6 +3072,8 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ mthsurf,divxprp_mn) CALL iscdftf(mfac,mpert,curvfun(istep,:), $ mthsurf,curv_mn) + CALL iscdftf(mfac,mpert,lagbprpfun(istep,:), + $ mthsurf,lagbprp_mn) c----------------------------------------------------------------------- c Save working coordinate verions for pent files. c----------------------------------------------------------------------- @@ -3084,10 +3094,12 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) CALL gpeq_bcoordsout(lagbparmout(i,:),lagbpar_mn,psi,tout,0) CALL gpeq_bcoordsout(divxprpmout(i,:),divxprp_mn,psi,tout,0) CALL gpeq_bcoordsout(curvmout(i,:) ,curv_mn,psi,tout,0) + CALL gpeq_bcoordsout(lagbprpmout(i,:),lagbprp_mn,psi,tout,0) eulbparfout(istep,:)=eulbparfun(istep,:)*EXP(ifac*nn*dphi) lagbparfout(istep,:)=lagbparfun(istep,:)*EXP(ifac*nn*dphi) divxprpfout(istep,:)=divxprpfun(istep,:)*EXP(ifac*nn*dphi) curvfout(istep,:)=curvfun(istep,:)*EXP(ifac*nn*dphi) + lagbprpfout(istep,:)=lagbprpfun(istep,:)*EXP(ifac*nn*dphi) ENDDO CALL cspline_dealloc(cspl) @@ -3134,6 +3146,17 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ "Lagrangian perturbed field") ) CALL check( nf90_put_att(fncid,bml_id,"units","Tesla") ) CALL check( nf90_put_att(fncid,bml_id,"jacobian",jac_out) ) + CALL check( nf90_def_var(fncid, "b_lag_perp_fun", nf90_double, + $ (/p_id,t_id,i_id/),bp_id) ) + CALL check( nf90_put_att(fncid,bp_id,"long_name", + $ "Lagrangian perturbed field of xi_perp") ) + CALL check( nf90_put_att(fncid,bp_id,"units","Tesla") ) + CALL check( nf90_def_var(fncid, "b_lag_perp", nf90_double, + $ (/p_id,m_id,i_id/),bmp_id) ) + CALL check( nf90_put_att(fncid,bmp_id,"long_name", + $ "Lagrangian perturbed field of xi_perp") ) + CALL check( nf90_put_att(fncid,bmp_id,"units","Tesla") ) + CALL check( nf90_put_att(fncid,bmp_id,"jacobian",jac_out) ) CALL check( nf90_def_var(fncid, "Bdivxi_perp_fun", nf90_double, $ (/p_id,t_id,i_id/),x_id) ) CALL check( nf90_put_att(fncid,x_id,"long_name", @@ -3165,6 +3188,11 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ eulbparmout),AIMAG(eulbparmout)/),(/mstep,lmpert,2/))) ) CALL check( nf90_put_var(fncid,bml_id,RESHAPE((/REAL( $ lagbparmout),AIMAG(lagbparmout)/),(/mstep,lmpert,2/))) ) + CALL check( nf90_put_var(fncid,bmp_id,RESHAPE((/REAL( + $ lagbprpmout),AIMAG(lagbprpmout)/),(/mstep,lmpert,2/))) ) + CALL check( nf90_put_var(fncid,bp_id,RESHAPE( + $ (/REAL(lagbprpfout),-helicity*AIMAG(lagbprpfout)/), + $ (/mstep,mthsurf,2/))) ) CALL check( nf90_put_var(fncid,xm_id,RESHAPE((/REAL( $ -divxprpmout),AIMAG(-divxprpmout)/),(/mstep,lmpert,2/))) ) CALL check( nf90_put_var(fncid,km_id,RESHAPE((/REAL( @@ -3385,7 +3413,8 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) DEALLOCATE(rs,zs,equilbfun) DEALLOCATE(eulbparmns,lagbparmns,divxprpmns,curvmns, $ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns) - DEALLOCATE(eulbparmout,lagbparmout,divxprpmout,curvmout) + DEALLOCATE(eulbparmout,lagbparmout,divxprpmout,curvmout, + $ lagbprpmout,lagbprpfun,lagbprpfout) DEALLOCATE(eulbparfun,lagbparfun,divxprpfun,curvfun, $ eulbparfout,lagbparfout,divxprpfout,curvfout) c----------------------------------------------------------------------- From c377c5fdd123372dff134136683d0f032194f65f Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:26:46 +0200 Subject: [PATCH 2/5] GPEC - NEW FEATURE - Emit optional perturbed Boozer and Hamada field strength --- .github/workflows/make_test.yml | 8 +++ docs/namelists.rst | 2 + docs/outputs.rst | 2 + docs/pmodb_target_gauge.rst | 77 +++++++++++++++++++++++++ docs/site.rst | 1 + gen/pmodb_gauge/generate.py | 12 ++++ gen/pmodb_gauge/prove.py | 65 +++++++++++++++++++++ gpec/gpec.f | 7 ++- gpec/gpout.f | 73 ++++++++++++++++++----- gpec/pmodb_gauge.f90 | 31 ++++++++++ gpec/pmodb_gauge_kernel.f90 | 19 ++++++ input/gpec.in | 1 + install/DEPENDENCIES.inc | 3 +- install/RULES.inc | 2 + regression/test/pmodb_target_compare.py | 20 +++++++ regression/test/pmodb_target_gauge.f90 | 49 ++++++++++++++++ 16 files changed, 354 insertions(+), 18 deletions(-) create mode 100644 docs/pmodb_target_gauge.rst create mode 100644 gen/pmodb_gauge/generate.py create mode 100644 gen/pmodb_gauge/prove.py create mode 100644 gpec/pmodb_gauge.f90 create mode 100644 gpec/pmodb_gauge_kernel.f90 create mode 100644 regression/test/pmodb_target_compare.py create mode 100644 regression/test/pmodb_target_gauge.f90 diff --git a/.github/workflows/make_test.yml b/.github/workflows/make_test.yml index 3598bf7e..8ecc4a7a 100644 --- a/.github/workflows/make_test.yml +++ b/.github/workflows/make_test.yml @@ -221,6 +221,14 @@ jobs: ulimit -s unlimited ../../../bin/dcon && ../../../bin/gpec && ../../../bin/stride && ../../../bin/pentrc + - name: Check perturbed-coordinate gauge operator + run: | + export LD_LIBRARY_PATH=$PREFIX/lib:$LD_LIBRARY_PATH + gfortran -O0 -fopenmp -Igpec -Iequil regression/test/pmodb_target_gauge.f90 gpec/pmodb_gauge.o gpec/pmodb_gauge_kernel.o gpec/ismath.o gpec/gpglobal.o -Ldeps/lib -lequil -llapack -lblas -o /tmp/pmodb_target_gauge + /tmp/pmodb_target_gauge + if /tmp/pmodb_target_gauge resonant > /tmp/resonant.log 2>&1; then exit 1; fi + python3 -c "assert 'Nonzero resonant pmodb gauge source' in open('/tmp/resonant.log').read()" + - name: Run Soloviev kinetic example run: | export GPECHOME=$(pwd) diff --git a/docs/namelists.rst b/docs/namelists.rst index 4629dade..9695e05a 100644 --- a/docs/namelists.rst +++ b/docs/namelists.rst @@ -45,3 +45,5 @@ This namelist sets the inputs, controls the behavior, and determines the outputs .. include:: ../input/pentrc.in :literal: + +See :doc:`pmodb_target_gauge` for the optional ``pmodb_gauge_flag`` output. diff --git a/docs/outputs.rst b/docs/outputs.rst index a0d80b2a..221a86d8 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -290,3 +290,5 @@ solutions.bin, solutions.bin **Flag** sol_num (match.in, MATCH) **Info** Ideal Euler-Lagrange solution built from DCON eigenfunction sol_num (1 is the most unstable). sol_num=1 writes solutions.bin; sol_num=n>1 writes solutions.bin so several eigenfunctions can be kept side by side. + +See :doc:`pmodb_target_gauge` for ``b_lag_perp`` and optional perturbed Boozer/Hamada field strength. diff --git a/docs/pmodb_target_gauge.rst b/docs/pmodb_target_gauge.rst new file mode 100644 index 00000000..e889fb8d --- /dev/null +++ b/docs/pmodb_target_gauge.rst @@ -0,0 +1,77 @@ +Optional perturbed-coordinate field strength +============================================= + +With ``pmodb_flag=t``, NetCDF ``b_lag_perp`` supplies +``eulb+xi_perp.grad B``. It is the perpendicular-displacement value called +``llagb`` in ``gpdiag_pmodb``. Existing ASCII columns keep their format. +This scalar alone is not the perturbed-Boozer input to NEO-RT or NEO-2. + +Set ``pmodb_gauge_flag=t`` in ``GPEC_OUTPUT`` for an additional NetCDF +field ``b_lag_boozer`` or ``b_lag_hamada``. This requires ``reg_flag=f`` +and a matching ``jac_type='boozer'`` or ``'hamada'``. The gauge is selected +by the working chart; ``jac_out`` selects only the harmonic representation. +The field has units Tesla and the same complex harmonic layout as ``b_lag``. +The output's ``gauge`` attribute records the physical gauge separately from +its ``jacobian`` representation. Use Boozer for NEO-type scalar-only +consumers and Hamada for Shaing-type formulas with that gauge contract. + +Verified derivation +------------------- + +GPEC angles are in turns and modes are exp(2 pi i(m theta-n zeta)). Write +``D=partial_theta+q partial_zeta``, so ``B.grad=(chi1/J)D``. On the matching +base ``J B^a=H(psi)``, with a=2 for Boozer or a=0 for Hamada. Under the +parallel shift ``xi_target=xi_X+lambda B``, ideal induction is unchanged, +``b_target=b_X+lambda B.grad B`` and ``d_target=d_X+B.grad lambda``. +The target condition for a nonaxisymmetric mode is ``d_target+a b_target/B=0``. +The product rule gives + +.. math:: + + T=d_X+a b_X/B,\qquad + T_{target}=T+\frac{\mathbf B\cdot\nabla(\lambda B^a)}{B^a},\qquad + w=\chi_1\lambda/J,\qquad + Dw=-T,\qquad w_m=\frac{iT_m}{2\pi(m-nq)},\qquad + b_{target}=b_X+w\partial_\theta B. + +The implementation forms T in real space, truncates to the existing GPEC +poloidal spectrum, inverts each nonresonant mode, and forms the scalar in +real space before the existing output-coordinate transform. Identities are +exact; numerical output retains GPEC's finite spectral/spline resolution. +A nonzero resonant source is rejected because no periodic solution exists. +At a resonant zero source the free homogeneous shift is set to zero. +The near-resonant result can be large: this feature applies no smoothing. + +The input pair is ``b_X=b_lag_perp`` and the ordinary perpendicular +coordinate divergence, using the derivative selected by gpeq_sol. Regularized fields/derivatives are deliberately +excluded because they need not form one ideal-induction displacement. + +Verification and scope +----------------------- + +``gen/pmodb_gauge/prove.py`` proves Q1-Q8 and K5/P1 using FortSym's native +zero test, with no probe fallback (verified at ae635c80). The inverse +coefficient is separately proved and emitted by +``gen/pmodb_gauge/generate.py`` to ``gpec/pmodb_gauge_kernel.f90``. +Regenerate and check the resulting diff; ordinary builds require no CAS. + +``regression/test/pmodb_target_gauge.f90`` exercises the actual operator: +parallel displacement gives zero in both gauges, a zero target source +leaves the scalar unchanged, and a nonzero resonant source fails. It also +writes a Hamada scalar transformed on Boozer equilibrium angles. +``regression/test/pmodb_target_compare.py`` compares it with rmp_torque's +c3614334/#55 kernel (external source supplied through PYTHONPATH). + +For the Hamada scalar already relabeled onto Boozer equilibrium angles, +T=2 b_H/B; both paths use the same inverse-field-line operator. Equilibrium +relabeling must precede this conversion. PEST output is not supported as a +target; its J/R^2 constraint also involves radial geometric displacement. +This feature does not infer a new equilibrium chart from a scalar spectrum. + +Source contracts: ``gpec/gpeq.f:114`` selects the derivative; +``gpec/gpdiag.f:1061`` defines the existing perpendicular scalar; +``gpec/gpout.f`` computes it before output transformation. The scalar-only +consumer formulas are those derived in rmp_torque c3614334, +``docs/boozer_gauge_ntv_derivation.md`` Eqs. 4.3-4.8 and 4.13-4.15. +Park 2011 PoP 18 110702 Eqs. 10-11 explain why a general kinetic action +with an explicit Jacobian term has a different input contract. diff --git a/docs/site.rst b/docs/site.rst index 7905ced0..a3b2fa15 100644 --- a/docs/site.rst +++ b/docs/site.rst @@ -16,3 +16,4 @@ GPEC documentation website site map source_documentation.rst installation.rst + pmodb_target_gauge.rst diff --git a/gen/pmodb_gauge/generate.py b/gen/pmodb_gauge/generate.py new file mode 100644 index 00000000..c9412fa7 --- /dev/null +++ b/gen/pmodb_gauge/generate.py @@ -0,0 +1,12 @@ +"""Prove and emit the inverse field-line coefficient (angles in turns).""" +from pathlib import Path +import fortsym +import fortsym.sympy as fs +s = fs.Symbol("s") +weight = 1/(2*fs.pi*s) +assert (2*fs.pi*s*weight-1).is_zero is True +source = fortsym.emit_kernel("pmodb_gauge_inverse", ["s"], {"weight":weight}, + generator="gen/pmodb_gauge/generate.py") +root = Path(__file__).resolve().parents[2] +(root/"gpec/pmodb_gauge_kernel.f90").write_text( + "module pmodb_gauge_kernel\ncontains\n"+source+"\nend module pmodb_gauge_kernel\n") diff --git a/gen/pmodb_gauge/prove.py b/gen/pmodb_gauge/prove.py new file mode 100644 index 00000000..306fdfb4 --- /dev/null +++ b/gen/pmodb_gauge/prove.py @@ -0,0 +1,65 @@ +import fortsym.sympy as sp + +def gauge_coordinates(sp): + """Local Jacobian and field-strength gauge equations in angles in turns. + + B.grad = chi1/J (d_theta + q d_zeta). For a target chart J B**a=H(psi), + a=0 (Hamada) or 2 (Boozer), w=chi1 lambda B**a/H obeys D w=-T, + T=div xi + a delta B_L/B, for n!=0 (zero flux-function perturbation). + """ + B, Bt, J, chi, q, lam, Lt, Lz, d, b, c = sp.symbols( + "B Bt J chi q lam Lt Lz d b c") + m,n=sp.symbols("m n") + out=[] + for a,name in ((0,"hamada"),(2,"boozer")): + # First jet of lambda B**a; axisymmetric equilibrium, div B=0. + Ftheta = Lt*B**a + a*lam*B**(a-1)*Bt + Fzeta = Lz*B**a + direct = d + chi/J*(Lt+q*Lz) + a*(b+lam*chi/J*Bt)/B + product = d+a*b/B+chi/(J*B**a)*(Ftheta+q*Fzeta) + out.append(("Q1."+name,"T_a(xi+lambda B) = T_a(xi)+B.grad(lambda B**a)/B**a",direct-product)) + # Solving the Fourier equation, with modes exp(2 pi i(m theta-n zeta)). + T,H=sp.symbols("T H") + w=-T/(2*sp.pi*sp.I*(m-n*q)) + out.append(("Q2."+name,"D w = -T; w_m = -T_m/[2 pi i(m-nq)]",2*sp.pi*sp.I*(m-n*q)*w+T)) + # On a target base J B**a=H, correction to field strength is w d_theta B. + lambda_w = w*J/chi + out.append(("Q3."+name,"delta B_target = delta B_X + w d_theta B (J B**a flux function)", + lambda_w*chi/J*Bt-w*Bt)) + # Transform the base angles: theta_B=theta_H+iota h, phi_B=phi_H+h. + h,Dh,iota,Hh,Hb=sp.symbols("h Dh iota Hh Hb") + ht,hp=sp.symbols("ht hp") + determinant=(1+iota*ht)*(1+hp)-ht*iota*hp + out.append(("Q4","det d(theta_H+iota h,phi_H+h)/d(theta_H,phi_H)=1+D h", + determinant-1-iota*ht-hp)) + # Same physical scalar and relabeling, in the convention of rmp_torque E4.14. + bX,jX,u,v=sp.symbols("bX jX u v") + dbB=bX+u*Bt; dbH=bX+v*Bt + out.append(("Q5","Hamada->Boozer equals direct X->Boozer when w_HB=w_XB-w_XH", + dbH+(u-v)*Bt-dbB)) + # Perpendicular reconstruction is unchanged by adding a field-aligned vector. + out.append(("Q6","field-aligned input: delta B_E=0, xi_perp=0 => target delta B=0", + lam*chi/J*Bt + (-lam*chi/J)*Bt)) + out.append(("Q7","xi_perp limit: zero target source implies no gauge correction", (bX+u*Bt-bX).subs(u,0))) + # PEST J/R**2 is a flux function. Its linearized constraint has radial motion. + R,xiR,divxi,db,dr=sp.symbols("R xiR divxi db dr") + eps=sp.Symbol("eps") + ratio=J*(1+eps*divxi)/(R+eps*xiR)**2 + out.append(("Q8","delta(J/R**2)/(J/R**2)=div xi-2 xi_R/R", + ratio.diff(eps).subs(eps,0)*R**2/J-divxi+2*xiR/R)) + # Complex CGL quadratic form: real and imaginary components independently. + Dr,Di,Kr,Ki=sp.symbols("Dr Di Kr Ki") + pressure=0; code=0 + for D,K in ((Dr,Kr),(Di,Ki)): + b=-(D+K); pp=-D+b; pa=-3*D-2*b + pressure-=((pa-pp)*b+pa*D)/2 + code+=(sp.Rational(5,3)*D**2+sp.Rational(1,3)*(D+3*K)**2)/2 + out.append(("K5","Complex CGL pressure work equals [5/3 |D|**2+1/3 |D+3K|**2]/2",pressure-code)) + # Pressure term in curvature is mu0 p' xi^psi/B, not mu0 p' xi^psi. + P,X=sp.symbols("P X") + out.append(("P1","(Bk - (llagb-eulb))/B = mu0 p' xi^psi/B**2", (P*X/B)/B-P*X/B**2)) + return out + +for label,equation,residual in gauge_coordinates(sp): + assert sp.sympify(residual).is_zero is True, label + print(label,"PROVED") diff --git a/gpec/gpec.f b/gpec/gpec.f index 840b2e50..c9677912 100644 --- a/gpec/gpec.f +++ b/gpec/gpec.f @@ -33,6 +33,7 @@ PROGRAM gpec_main CHARACTER(8) :: filter_types CHARACTER(128) :: infile LOGICAL :: singcoup_flag,singfld_flag,vsingfld_flag,pmodb_flag, + $ pmodb_gauge_flag, $ xbcontra_flag,xbnormal_flag,vbnormal_flag,xbnobo_flag, $ d3_flag,xbst_flag,rzphibx_flag,jacfac_flag,dw_flag, $ radvar_flag,eigen_flag,magpot_flag,xbtangent_flag, @@ -64,7 +65,8 @@ PROGRAM gpec_main $ reg_spot,chebyshev_flag,nche,nchr,nchz,use_classic_splines, $ use_res_spot NAMELIST/gpec_output/resp_flag,singcoup_flag,nrzeq_flag,nr,nz, - $ singfld_flag,pmodb_flag,xbnormal_flag,rstep,jsurf_out, + $ singfld_flag,pmodb_flag,pmodb_gauge_flag,xbnormal_flag, + $ rstep,jsurf_out, $ jac_out,power_bout,power_rout,power_bpout,power_rcout, $ tmag_out,mlim_out,eqbrzphi_flag,brzphi_flag,xrzphi_flag, $ vbrzphi_flag,vvbrzphi_flag,divzero_flag,dw_flag,opsi1,opsi2, @@ -144,6 +146,7 @@ PROGRAM gpec_main singfld_flag=.TRUE. vsingfld_flag=.FALSE. pmodb_flag=.FALSE. + pmodb_gauge_flag=.FALSE. xbnormal_flag=.TRUE. vbnormal_flag=.FALSE. xbtangent_flag=.FALSE. @@ -622,7 +625,7 @@ PROGRAM gpec_main CALL gpout_dw_matrix(coil_flag) ENDIF IF (pmodb_flag) THEN - CALL gpout_pmodb(mode,xspmn) + CALL gpout_pmodb(mode,xspmn,pmodb_gauge_flag) ENDIF IF (xbnormal_flag) THEN CALL gpout_xbnormal(mode,xspmn,sing_spots,sing_interpspot, diff --git a/gpec/gpout.f b/gpec/gpout.f index 021cc0ee..2d3e5baf 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -37,6 +37,7 @@ c declarations. c----------------------------------------------------------------------- MODULE gpout_mod + USE pmodb_gauge_mod, ONLY: pmodb_target_gauge USE netcdf USE gpglobal_mod, ONLY : recon_out USE gpresp_mod @@ -2923,18 +2924,26 @@ END SUBROUTINE gpout_dw_matrix c subprogram 8. gpout_pmodb. c compute perturbed mod b. c----------------------------------------------------------------------- - SUBROUTINE gpout_pmodb(egnum,xspmn) + SUBROUTINE gpout_pmodb(egnum,xspmn,gauge) c----------------------------------------------------------------------- c declaration. c----------------------------------------------------------------------- INTEGER, INTENT(IN) :: egnum COMPLEX(r8), DIMENSION(mpert), INTENT(IN) :: xspmn + LOGICAL, INTENT(IN) :: gauge + INTEGER :: gauge_power,target_id + COMPLEX(r8), DIMENSION(mpert) :: target_mn + COMPLEX(r8), DIMENSION(0:mthsurf) :: source_fun,gauge_div, + $ actual_dx + REAL(r8), DIMENSION(0:mthsurf) :: btheta_fun + COMPLEX(r8), ALLOCATABLE :: target_out(:,:) + CHARACTER(16) :: gauge_name INTEGER :: i,istep,ipert,itheta,cstep,tout REAL(r8) :: jac,psi INTEGER :: p_id,t_id,i_id,m_id,r_id,z_id,b_id,bme_id,be_id, - $ bml_id,bl_id,xm_id,x_id,km_id,k_id,rzstat,bmp_id,bp_id + $ bml_id,bl_id,xm_id,x_id,km_id,k_id,rzstat,bmp_id REAL(r8), DIMENSION(0:mthsurf) :: dphi REAL(r8), DIMENSION(:,:), ALLOCATABLE :: rs,zs,equilbfun @@ -2949,7 +2958,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ lagbparmout,divxprpmout,curvmout,lagbprpmout COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparfun,lagbparfun, $ divxprpfun,curvfun,eulbparfout,lagbparfout,divxprpfout, - $ curvfout,lagbprpfun,lagbprpfout + $ curvfout,lagbprpfun REAL(r8), DIMENSION(:), ALLOCATABLE :: psis,ches,chex COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE ::chea,chelagbmns @@ -2971,10 +2980,25 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ divxprpfun(mstep,0:mthsurf),curvfun(mstep,0:mthsurf), $ eulbparfout(mstep,0:mthsurf),lagbparfout(mstep,0:mthsurf), $ divxprpfout(mstep,0:mthsurf),curvfout(mstep,0:mthsurf), - $ lagbprpfun(mstep,0:mthsurf),lagbprpfout(mstep,0:mthsurf)) + $ lagbprpfun(mstep,0:mthsurf)) c----------------------------------------------------------------------- c compute necessary components. c----------------------------------------------------------------------- + IF(gauge)THEN + IF(.NOT.netcdf_flag) ERROR STOP 'pmodb gauge needs NetCDF' + IF(reg_flag) ERROR STOP 'pmodb gauge needs reg_flag=f' + SELECT CASE(jac_type) + CASE('boozer') + gauge_power=2 + gauge_name='b_lag_boozer' + CASE('hamada') + gauge_power=0 + gauge_name='b_lag_hamada' + CASE DEFAULT + ERROR STOP 'pmodb gauge needs boozer or hamada jac_type' + END SELECT + ALLOCATE(target_out(mstep,lmpert)) + ENDIF IF(timeit) CALL gpec_timer(-2) IF(verbose) WRITE(*,*)"Computing |b| and del.x_prp" CALL cspline_alloc(cspl,mthsurf,2) @@ -3010,6 +3034,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) singfac=mfac-nn*sq%f(4) xsp1_mn=xsp1_mn*(singfac**2/(singfac**2+reg_spot**2)) CALL iscdftb(mfac,mpert,xsp1_fun,mthsurf,xsp1_mn) + IF(gauge) CALL iscdftb(mfac,mpert,actual_dx,mthsurf,xmp1_mn) DO itheta=0,mthsurf CALL bicube_eval(eqfun,psifac(istep),theta(itheta),0) @@ -3045,7 +3070,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ xsp_fun(itheta)*eqfun%fx(1)+ $ xms_fun(itheta)/(chi1*sq%f(4))*eqfun%fy(1) c Lagrangian dB of xi_perp alone: lagb plus the gauge term -c (xi_perp^zeta/q) d_theta B (docs/pmodb_gauge.rst). +c (xi_perp^zeta/q) d_theta B (docs/pmodb_target_gauge.rst). lagbprpfun(istep,itheta)=lagbparfun(istep,itheta)+ $ eqfun%fy(1)*(xmz_fun(itheta)/(jac*sq%f(4))- $ (chi1/(jac*eqfun%f(1)))**2* @@ -3053,6 +3078,8 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) divxprpfun(istep,itheta)=xsp1_fun(itheta)+ $ (jac1/jac)*xsp_fun(itheta)+ $ cspl%f1(1)/jac-(twopi*ifac*nn)*cspl%f(2)/jac + IF(gauge) gauge_div(itheta)=divxprpfun(istep,itheta)+ + $ actual_dx(itheta)-xsp1_fun(itheta) divxprpfun(istep,itheta)=-eqfun%f(1)* $ divxprpfun(istep,itheta) curvfun(istep,itheta)=-xsp_fun(itheta)/eqfun%f(1)*sq%f1(2)- @@ -3064,6 +3091,19 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) equilbfun(istep,itheta) = eqfun%f(1) dphi(itheta) = rzphi%f(3) ENDDO + IF(gauge)THEN + DO itheta=0,mthsurf + CALL bicube_eval(eqfun,psifac(istep),theta(itheta),1) + btheta_fun(itheta)=eqfun%fy(1) + ENDDO + CALL pmodb_target_gauge(mfac,nn,sq%f(4), + $ equilbfun(istep,:),btheta_fun,lagbprpfun(istep,:), + $ gauge_div,gauge_power, + $ source_fun) + CALL iscdftf(mfac,mpert,source_fun,mthsurf,target_mn) + CALL gpeq_bcoordsout(target_out(istep,:),target_mn, + $ psifac(istep),tout,0) + ENDIF CALL iscdftf(mfac,mpert,eulbparfun(istep,:), $ mthsurf,eulbpar_mn) CALL iscdftf(mfac,mpert,lagbparfun(istep,:), @@ -3099,7 +3139,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) lagbparfout(istep,:)=lagbparfun(istep,:)*EXP(ifac*nn*dphi) divxprpfout(istep,:)=divxprpfun(istep,:)*EXP(ifac*nn*dphi) curvfout(istep,:)=curvfun(istep,:)*EXP(ifac*nn*dphi) - lagbprpfout(istep,:)=lagbprpfun(istep,:)*EXP(ifac*nn*dphi) + ENDDO CALL cspline_dealloc(cspl) @@ -3146,17 +3186,19 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ "Lagrangian perturbed field") ) CALL check( nf90_put_att(fncid,bml_id,"units","Tesla") ) CALL check( nf90_put_att(fncid,bml_id,"jacobian",jac_out) ) - CALL check( nf90_def_var(fncid, "b_lag_perp_fun", nf90_double, - $ (/p_id,t_id,i_id/),bp_id) ) - CALL check( nf90_put_att(fncid,bp_id,"long_name", - $ "Lagrangian perturbed field of xi_perp") ) - CALL check( nf90_put_att(fncid,bp_id,"units","Tesla") ) CALL check( nf90_def_var(fncid, "b_lag_perp", nf90_double, $ (/p_id,m_id,i_id/),bmp_id) ) CALL check( nf90_put_att(fncid,bmp_id,"long_name", $ "Lagrangian perturbed field of xi_perp") ) CALL check( nf90_put_att(fncid,bmp_id,"units","Tesla") ) CALL check( nf90_put_att(fncid,bmp_id,"jacobian",jac_out) ) + IF(gauge)THEN + CALL check(nf90_def_var(fncid,TRIM(gauge_name),nf90_double, + $ (/p_id,m_id,i_id/),target_id)) + CALL check(nf90_put_att(fncid,target_id,'units','Tesla')) + CALL check(nf90_put_att(fncid,target_id,'gauge',jac_type)) + CALL check(nf90_put_att(fncid,target_id,'jacobian',jac_out)) + ENDIF CALL check( nf90_def_var(fncid, "Bdivxi_perp_fun", nf90_double, $ (/p_id,t_id,i_id/),x_id) ) CALL check( nf90_put_att(fncid,x_id,"long_name", @@ -3190,9 +3232,9 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ lagbparmout),AIMAG(lagbparmout)/),(/mstep,lmpert,2/))) ) CALL check( nf90_put_var(fncid,bmp_id,RESHAPE((/REAL( $ lagbprpmout),AIMAG(lagbprpmout)/),(/mstep,lmpert,2/))) ) - CALL check( nf90_put_var(fncid,bp_id,RESHAPE( - $ (/REAL(lagbprpfout),-helicity*AIMAG(lagbprpfout)/), - $ (/mstep,mthsurf,2/))) ) + IF(gauge) CALL check(nf90_put_var(fncid,target_id,RESHAPE( + $ (/REAL(target_out),AIMAG(target_out)/), + $ (/mstep,lmpert,2/)))) CALL check( nf90_put_var(fncid,xm_id,RESHAPE((/REAL( $ -divxprpmout),AIMAG(-divxprpmout)/),(/mstep,lmpert,2/))) ) CALL check( nf90_put_var(fncid,km_id,RESHAPE((/REAL( @@ -3410,11 +3452,12 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) c----------------------------------------------------------------------- c deallocation cleans memory in heap c----------------------------------------------------------------------- + IF(gauge) DEALLOCATE(target_out) DEALLOCATE(rs,zs,equilbfun) DEALLOCATE(eulbparmns,lagbparmns,divxprpmns,curvmns, $ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns) DEALLOCATE(eulbparmout,lagbparmout,divxprpmout,curvmout, - $ lagbprpmout,lagbprpfun,lagbprpfout) + $ lagbprpmout,lagbprpfun) DEALLOCATE(eulbparfun,lagbparfun,divxprpfun,curvfun, $ eulbparfout,lagbparfout,divxprpfout,curvfout) c----------------------------------------------------------------------- diff --git a/gpec/pmodb_gauge.f90 b/gpec/pmodb_gauge.f90 new file mode 100644 index 00000000..ae32c768 --- /dev/null +++ b/gpec/pmodb_gauge.f90 @@ -0,0 +1,31 @@ +module pmodb_gauge_mod + use local_mod, only: r8 + use ismath_mod, only: iscdftf, iscdftb + use pmodb_gauge_kernel, only: pmodb_gauge_inverse + implicit none +contains + subroutine pmodb_target_gauge(m, n, q, b, bt, delta, divxi, power, target) + integer, intent(in) :: m(:), n, power + real(r8), intent(in) :: q, b(0:), bt(0:) + complex(r8), intent(in) :: delta(0:), divxi(0:) + complex(r8), intent(out) :: target(0:) + complex(r8) :: source(size(m)), w(0:size(b)-1) + real(r8) :: s, weight + integer :: j, nt + nt = size(b)-1 + call iscdftf(m, size(m), divxi+power*delta/b, nt, source) + do j = 1, size(m) + s = m(j)-n*q + if (abs(s) < 1.0e-10_r8) then + if (abs(source(j)) > 1.0e-12_r8*max(maxval(abs(source)), tiny(1.0_r8))) & + error stop 'Nonzero resonant pmodb gauge source' + source(j) = 0 + else + call pmodb_gauge_inverse(s, weight) + source(j) = cmplx(0.0_r8, 1.0_r8, r8)*weight*source(j) + end if + end do + call iscdftb(m, size(m), w, nt, source) + target = delta+w*bt + end subroutine pmodb_target_gauge +end module pmodb_gauge_mod diff --git a/gpec/pmodb_gauge_kernel.f90 b/gpec/pmodb_gauge_kernel.f90 new file mode 100644 index 00000000..e71e8629 --- /dev/null +++ b/gpec/pmodb_gauge_kernel.f90 @@ -0,0 +1,19 @@ +module pmodb_gauge_kernel +contains +! Generated by fortsym. Do not edit. +! Generator: gen/pmodb_gauge/generate.py +! Regenerate with: rerun the generating script + +subroutine pmodb_gauge_inverse(s, weight) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + real(real64), intent(in) :: s + real(real64), intent(out) :: weight + real(real64) :: t4, t6 + + t4 = (s * acos(-1.0_real64) * 2.0000000000000000E+000_real64) + t6 = ((1.0000000000000000E+000_real64 / t4)) + weight = t6 +end subroutine pmodb_gauge_inverse + +end module pmodb_gauge_kernel diff --git a/input/gpec.in b/input/gpec.in index a4b43aff..71c60ff9 100644 --- a/input/gpec.in +++ b/input/gpec.in @@ -61,6 +61,7 @@ singthresh_slayer_inpr=5.0 ! Scalar Prandtl number (viscosity) used by SLAYER for all rational surfaces singthresh_slayer_inpr_prof=-1.0 ! Per-rational-surface Prandtl number profile, comma-separated and ordered by ascending q, up to 20 entries (e.g. singthresh_slayer_inpr_prof=3.0,5.0,8.0,12.0 for the q=2,3,4,5 surfaces). Any entry <= 0 falls back to the scalar singthresh_slayer_inpr + pmodb_gauge_flag=f ! Extra target-gauge NetCDF field; reg_flag=f, Boozer/Hamada pmodb_flag=f ! Outputs delta-B_Lagrangian xbnormal_flag=t ! Outputs normal displacement and field in plasma vbnormal_flag=t ! Outputs vacuum normal displacement and field in plasma diff --git a/install/DEPENDENCIES.inc b/install/DEPENDENCIES.inc index d842459f..3adebae4 100644 --- a/install/DEPENDENCIES.inc +++ b/install/DEPENDENCIES.inc @@ -197,11 +197,12 @@ ../gpec/gpec.o: ../coil/coil.o ../gpec/gpdiag.o ../gpec/gpout.o ../gpec/rdcon.o ../gpec/gpeq.o: ../gpec/idcon.o ../gpec/gpglobal.o: ../equil/bicube.o ../equil/fspline.o ../equil/local.o -../gpec/gpout.o: ../coil/field.o ../gpec/gpdiag.o ../gpec/gpglobal.o ../gpec/gpresp.o ../gpec/gpvacuum.o ../gpec/idcon.o ../gpec/version.inc ../pentrc/inputs.o ../pentrc/pentrc_interface.o ../pentrc/utilities.o ../slayer/gslayer.o +../gpec/gpout.o: ../coil/field.o ../gpec/gpdiag.o ../gpec/gpglobal.o ../gpec/gpresp.o ../gpec/gpvacuum.o ../gpec/idcon.o ../gpec/pmodb_gauge.o ../gpec/version.inc ../pentrc/inputs.o ../pentrc/pentrc_interface.o ../pentrc/utilities.o ../slayer/gslayer.o ../gpec/gpresp.o: ../gpec/gpeq.o ../gpec/gpvacuum.o: ../gpec/gpeq.o ../vacuum/vacuum_ma.o ../gpec/idcon.o: ../equil/equil.o ../gpec/gpglobal.o ../gpec/ismath.o ../vacuum/vacuum_ma.o ../gpec/ismath.o: ../equil/bicube.o ../equil/fspline.o ../equil/local.o ../gpec/gpglobal.o +../gpec/pmodb_gauge.o: ../equil/local.o ../gpec/ismath.o ../gpec/pmodb_gauge_kernel.o ../gpec/rdcon.o: ../gpec/gpglobal.o ../gpec/idcon.o ../gpec/ismath.o # XDraw diff --git a/install/RULES.inc b/install/RULES.inc index be30ba9b..0d64294e 100644 --- a/install/RULES.inc +++ b/install/RULES.inc @@ -140,6 +140,8 @@ FCGPEC = $(FC) $(FFLAGS) $(OMPFLAG) $(RECURSFLAG) $(LEGACYFLAG) -I$(MATHINC) -I$ $(FCGPEC) $(MODFLAG) ../gpec -c $< -o $@ ../gpec/%.o: ../gpec/%.F $(FCGPEC) $(MODFLAG) ../gpec -c $< -o $@ +../gpec/%.o: ../gpec/%.f90 + $(FCGPEC) $(MODFLAG) ../gpec -c $< -o $@ # XDRAW CCXDRAW = $(CC) $(X11_INC_DIR) $(CFLAGS) -DUNIX -O2 diff --git a/regression/test/pmodb_target_compare.py b/regression/test/pmodb_target_compare.py new file mode 100644 index 00000000..028c5337 --- /dev/null +++ b/regression/test/pmodb_target_compare.py @@ -0,0 +1,20 @@ +"""Cross-check the native gauge test's Hamada->Boozer result with rmp_torque. + +Run pmodb_target_gauge first. Set PYTHONPATH to rmp_torque #55 source. +""" +import numpy as np +from rmp_torque.boozer_gauge import Spectrum, hamada_to_boozer + +rows = np.loadtxt('pmodb_hamada_to_boozer.txt') +modes = np.arange(-12, 13) +values = np.zeros(25, complex) +values[13] = 0.2 +bvalues = np.zeros(25, complex) +bvalues[12] = 3 +bvalues[11] = bvalues[13] = -0.3 +# GPEC exp(i(m theta-n zeta)) maps to ITP n=-2; both angles here in radians. +expected = hamada_to_boozer(Spectrum(modes, values), Spectrum(modes,bvalues), -2, 1/1.3, 12) +actual = rows[:,2]+1j*rows[:,3] +reference = expected.evaluate(2*np.pi*rows[:,0]) +np.testing.assert_allclose(actual, reference, atol=2e-11, rtol=2e-11) +print('PASS: native Hamada->Boozer output agrees with rmp_torque #55') diff --git a/regression/test/pmodb_target_gauge.f90 b/regression/test/pmodb_target_gauge.f90 new file mode 100644 index 00000000..465209ec --- /dev/null +++ b/regression/test/pmodb_target_gauge.f90 @@ -0,0 +1,49 @@ +program test_pmodb_target_gauge + use local_mod, only: r8 + use pmodb_gauge_mod, only: pmodb_target_gauge + implicit none + integer, parameter :: nt=128, nm=25 + integer :: m(nm), i, power + real(r8) :: th, q, pi, b(0:nt), bt(0:nt) + complex(r8) :: u(0:nt), du(0:nt), delta(0:nt), d(0:nt), result(0:nt) + character(16) :: mode + call get_command_argument(1,mode) + pi=acos(-1.0_r8) + q=1.3_r8 + m=[(i,i=-12,12)] + do i=0,nt + th=2*pi*i/nt + b(i)=3.0_r8*(1-0.2_r8*cos(th)) + bt(i)=2*pi*0.6_r8*sin(th) + u(i)=exp(cmplx(0.0_r8,th,r8)) + du(i)=cmplx(0.0_r8,2*pi*(1-2*q),r8)*u(i) + end do + if (mode=='resonant') then + d=exp(cmplx(0.0_r8,4*pi,r8)*[(real(i,r8)/nt,i=0,nt)]) + delta=0 + call pmodb_target_gauge(m,2,1.0_r8,b,bt,delta,d,0,result) + error stop 'nonzero resonant source was accepted' + end if + do power=0,2,2 + ! A parallel displacement has delta=u*d_theta B and T=D u. + delta=u*bt + d=du-power*delta/b + call pmodb_target_gauge(m,2,q,b,bt,delta,d,power,result) + if (maxval(abs(result))>1.e-11_r8) error stop 'parallel displacement is not zero' + ! With T=0 the target scalar must stay unchanged. + delta=0.4_r8*u + d=-power*delta/b + call pmodb_target_gauge(m,2,q,b,bt,delta,d,power,result) + if (maxval(abs(result-delta))>1.e-12_r8) error stop 'zero-source limit changed' + end do + ! Hamada scalar already relabeled to Boozer base angles; div xi_H=0. + delta=0.2_r8*u + d=0 + call pmodb_target_gauge(m,2,q,b,bt,delta,d,2,result) + open(unit=31,file='pmodb_hamada_to_boozer.txt',status='replace') + do i=0,nt-1 + write(31,'(5es24.15)') real(i,r8)/nt,b(i),real(result(i)),aimag(result(i)),bt(i) + end do + close(31) + print '(a)', 'PASS: parallel displacement and zero-source limits in both gauges' +end program test_pmodb_target_gauge From a42df548920a64b27022b14f7ce6cd6c65a9dcb1 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:32:54 +0200 Subject: [PATCH 3/5] GPEC - NEW FEATURE - Emit gauge leaves and restrict exact labels to unregularized output --- docs/pmodb_target_gauge.rst | 9 ++++--- gen/pmodb_gauge/generate.py | 21 ++++++++++----- gpec/gpout.f | 5 +++- gpec/pmodb_gauge.f90 | 22 +++++++++++----- gpec/pmodb_gauge_kernel.f90 | 52 ++++++++++++++++++++++++++++++++----- 5 files changed, 84 insertions(+), 25 deletions(-) diff --git a/docs/pmodb_target_gauge.rst b/docs/pmodb_target_gauge.rst index e889fb8d..6a1afe43 100644 --- a/docs/pmodb_target_gauge.rst +++ b/docs/pmodb_target_gauge.rst @@ -1,9 +1,10 @@ Optional perturbed-coordinate field strength ============================================= -With ``pmodb_flag=t``, NetCDF ``b_lag_perp`` supplies +With ``pmodb_flag=t`` and ``reg_flag=f``, NetCDF ``b_lag_perp`` supplies ``eulb+xi_perp.grad B``. It is the perpendicular-displacement value called -``llagb`` in ``gpdiag_pmodb``. Existing ASCII columns keep their format. +``llagb`` in ``gpdiag_pmodb``. It is omitted for regularized solutions, whose reconstructed surface +displacement need not equal xms/chi1. Existing ASCII columns keep their format. This scalar alone is not the perturbed-Boozer input to NEO-RT or NEO-2. Set ``pmodb_gauge_flag=t`` in ``GPEC_OUTPUT`` for an additional NetCDF @@ -50,8 +51,8 @@ Verification and scope ----------------------- ``gen/pmodb_gauge/prove.py`` proves Q1-Q8 and K5/P1 using FortSym's native -zero test, with no probe fallback (verified at ae635c80). The inverse -coefficient is separately proved and emitted by +zero test, with no probe fallback (verified at ae635c80). The source, complex inverse +coefficient and scalar correction are separately proved and emitted by ``gen/pmodb_gauge/generate.py`` to ``gpec/pmodb_gauge_kernel.f90``. Regenerate and check the resulting diff; ordinary builds require no CAS. diff --git a/gen/pmodb_gauge/generate.py b/gen/pmodb_gauge/generate.py index c9412fa7..f22a6ffa 100644 --- a/gen/pmodb_gauge/generate.py +++ b/gen/pmodb_gauge/generate.py @@ -1,12 +1,21 @@ -"""Prove and emit the inverse field-line coefficient (angles in turns).""" +"""Prove and emit the real leaves; Fourier indexing stays in pmodb_gauge_mod.""" from pathlib import Path import fortsym import fortsym.sympy as fs -s = fs.Symbol("s") +s, Tr, Ti = fs.symbols('s Tr Ti') weight = 1/(2*fs.pi*s) +wr, wi = -weight*Ti, weight*Tr assert (2*fs.pi*s*weight-1).is_zero is True -source = fortsym.emit_kernel("pmodb_gauge_inverse", ["s"], {"weight":weight}, - generator="gen/pmodb_gauge/generate.py") +assert (2*fs.pi*fs.I*s*(wr+fs.I*wi)+Tr+fs.I*Ti).is_zero is True +a, b, bt, dr, di, vr, vi, ur, ui = fs.symbols('a b bt dr di vr vi ur ui') +source = {'tr': vr+a*dr/b, 'ti': vi+a*di/b} +target = {'out_r': dr+ur*bt, 'out_i': di+ui*bt} +assert (source['tr']+fs.I*source['ti']-(vr+fs.I*vi)-a*(dr+fs.I*di)/b).is_zero is True +assert (target['out_r']+fs.I*target['out_i']-(dr+fs.I*di)-(ur+fs.I*ui)*bt).is_zero is True +kernels = [ + fortsym.emit_kernel('pmodb_gauge_inverse', ['s','Tr','Ti'], {'wr':wr,'wi':wi}, generator='gen/pmodb_gauge/generate.py'), + fortsym.emit_kernel('pmodb_gauge_source', ['a','b','dr','di','vr','vi'], source, generator='gen/pmodb_gauge/generate.py'), + fortsym.emit_kernel('pmodb_gauge_apply', ['dr','di','ur','ui','bt'], target, generator='gen/pmodb_gauge/generate.py')] root = Path(__file__).resolve().parents[2] -(root/"gpec/pmodb_gauge_kernel.f90").write_text( - "module pmodb_gauge_kernel\ncontains\n"+source+"\nend module pmodb_gauge_kernel\n") +(root/'gpec/pmodb_gauge_kernel.f90').write_text('module pmodb_gauge_kernel\ncontains\n'+'\n'.join(kernels)+'\nend module pmodb_gauge_kernel\n') +print('PROVED and emitted: complex source, inverse field-line coefficient, scalar correction') diff --git a/gpec/gpout.f b/gpec/gpout.f index 2d3e5baf..4b80a051 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -3186,12 +3186,14 @@ SUBROUTINE gpout_pmodb(egnum,xspmn,gauge) $ "Lagrangian perturbed field") ) CALL check( nf90_put_att(fncid,bml_id,"units","Tesla") ) CALL check( nf90_put_att(fncid,bml_id,"jacobian",jac_out) ) + IF(.NOT.reg_flag)THEN CALL check( nf90_def_var(fncid, "b_lag_perp", nf90_double, $ (/p_id,m_id,i_id/),bmp_id) ) CALL check( nf90_put_att(fncid,bmp_id,"long_name", $ "Lagrangian perturbed field of xi_perp") ) CALL check( nf90_put_att(fncid,bmp_id,"units","Tesla") ) CALL check( nf90_put_att(fncid,bmp_id,"jacobian",jac_out) ) + ENDIF IF(gauge)THEN CALL check(nf90_def_var(fncid,TRIM(gauge_name),nf90_double, $ (/p_id,m_id,i_id/),target_id)) @@ -3230,7 +3232,8 @@ SUBROUTINE gpout_pmodb(egnum,xspmn,gauge) $ eulbparmout),AIMAG(eulbparmout)/),(/mstep,lmpert,2/))) ) CALL check( nf90_put_var(fncid,bml_id,RESHAPE((/REAL( $ lagbparmout),AIMAG(lagbparmout)/),(/mstep,lmpert,2/))) ) - CALL check( nf90_put_var(fncid,bmp_id,RESHAPE((/REAL( + IF(.NOT.reg_flag) + $ CALL check( nf90_put_var(fncid,bmp_id,RESHAPE((/REAL( $ lagbprpmout),AIMAG(lagbprpmout)/),(/mstep,lmpert,2/))) ) IF(gauge) CALL check(nf90_put_var(fncid,target_id,RESHAPE( $ (/REAL(target_out),AIMAG(target_out)/), diff --git a/gpec/pmodb_gauge.f90 b/gpec/pmodb_gauge.f90 index ae32c768..ca907de5 100644 --- a/gpec/pmodb_gauge.f90 +++ b/gpec/pmodb_gauge.f90 @@ -1,7 +1,7 @@ module pmodb_gauge_mod use local_mod, only: r8 use ismath_mod, only: iscdftf, iscdftb - use pmodb_gauge_kernel, only: pmodb_gauge_inverse + use pmodb_gauge_kernel, only: pmodb_gauge_inverse, pmodb_gauge_source, pmodb_gauge_apply implicit none contains subroutine pmodb_target_gauge(m, n, q, b, bt, delta, divxi, power, target) @@ -9,11 +9,16 @@ subroutine pmodb_target_gauge(m, n, q, b, bt, delta, divxi, power, target) real(r8), intent(in) :: q, b(0:), bt(0:) complex(r8), intent(in) :: delta(0:), divxi(0:) complex(r8), intent(out) :: target(0:) - complex(r8) :: source(size(m)), w(0:size(b)-1) - real(r8) :: s, weight + complex(r8) :: source(size(m)), w(0:size(b)-1), t(0:size(b)-1) + real(r8) :: s, out_r, out_i integer :: j, nt nt = size(b)-1 - call iscdftf(m, size(m), divxi+power*delta/b, nt, source) + do j = 0, nt + call pmodb_gauge_source(real(power,r8), b(j), real(delta(j),r8), aimag(delta(j)), & + real(divxi(j),r8), aimag(divxi(j)), out_r, out_i) + t(j) = cmplx(out_r, out_i, r8) + end do + call iscdftf(m, size(m), t, nt, source) do j = 1, size(m) s = m(j)-n*q if (abs(s) < 1.0e-10_r8) then @@ -21,11 +26,14 @@ subroutine pmodb_target_gauge(m, n, q, b, bt, delta, divxi, power, target) error stop 'Nonzero resonant pmodb gauge source' source(j) = 0 else - call pmodb_gauge_inverse(s, weight) - source(j) = cmplx(0.0_r8, 1.0_r8, r8)*weight*source(j) + call pmodb_gauge_inverse(s, real(source(j),r8), aimag(source(j)), out_r, out_i) + source(j) = cmplx(out_r, out_i, r8) end if end do call iscdftb(m, size(m), w, nt, source) - target = delta+w*bt + do j = 0, nt + call pmodb_gauge_apply(real(delta(j),r8), aimag(delta(j)), real(w(j),r8), aimag(w(j)), bt(j), out_r, out_i) + target(j) = cmplx(out_r, out_i, r8) + end do end subroutine pmodb_target_gauge end module pmodb_gauge_mod diff --git a/gpec/pmodb_gauge_kernel.f90 b/gpec/pmodb_gauge_kernel.f90 index e71e8629..8eb51784 100644 --- a/gpec/pmodb_gauge_kernel.f90 +++ b/gpec/pmodb_gauge_kernel.f90 @@ -4,16 +4,54 @@ module pmodb_gauge_kernel ! Generator: gen/pmodb_gauge/generate.py ! Regenerate with: rerun the generating script -subroutine pmodb_gauge_inverse(s, weight) +subroutine pmodb_gauge_inverse(s, Tr, Ti, wr, wi) use, intrinsic :: iso_fortran_env, only: real64 implicit none - real(real64), intent(in) :: s - real(real64), intent(out) :: weight - real(real64) :: t4, t6 + real(real64), intent(in) :: s, Tr, Ti + real(real64), intent(out) :: wr, wi + real(real64) :: t5, t7, t9, t11 - t4 = (s * acos(-1.0_real64) * 2.0000000000000000E+000_real64) - t6 = ((1.0000000000000000E+000_real64 / t4)) - weight = t6 + t5 = (s * acos(-1.0_real64) * 2.0000000000000000E+000_real64) + t7 = ((1.0000000000000000E+000_real64 / t5)) + t9 = (Ti * t7 * (-1.0000000000000000E+000_real64)) + t11 = (Tr * t7) + wr = t9 + wi = t11 end subroutine pmodb_gauge_inverse +! Generated by fortsym. Do not edit. +! Generator: gen/pmodb_gauge/generate.py +! Regenerate with: rerun the generating script + +subroutine pmodb_gauge_source(a, b, dr, di, vr, vi, tr, ti) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + real(real64), intent(in) :: a, b, dr, di, vr, vi + real(real64), intent(out) :: tr, ti + real(real64) :: t6, t8, t12 + + t6 = ((1.0000000000000000E+000_real64 / b)) + t8 = (vr + (a * dr * t6)) + t12 = (vi + (a * di * t6)) + tr = t8 + ti = t12 +end subroutine pmodb_gauge_source + +! Generated by fortsym. Do not edit. +! Generator: gen/pmodb_gauge/generate.py +! Regenerate with: rerun the generating script + +subroutine pmodb_gauge_apply(dr, di, ur, ui, bt, out_r, out_i) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + real(real64), intent(in) :: dr, di, ur, ui, bt + real(real64), intent(out) :: out_r, out_i + real(real64) :: t5, t9 + + t5 = (dr + (bt * ur)) + t9 = (di + (bt * ui)) + out_r = t5 + out_i = t9 +end subroutine pmodb_gauge_apply + end module pmodb_gauge_kernel From 27970bd69c26d29b362c8bc71dccfa9eaf46ed54 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:37:43 +0200 Subject: [PATCH 4/5] DOCS - MINOR - Pin target-gauge source references --- docs/pmodb_target_gauge.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/pmodb_target_gauge.rst b/docs/pmodb_target_gauge.rst index 6a1afe43..f20d4786 100644 --- a/docs/pmodb_target_gauge.rst +++ b/docs/pmodb_target_gauge.rst @@ -69,7 +69,7 @@ relabeling must precede this conversion. PEST output is not supported as a target; its J/R^2 constraint also involves radial geometric displacement. This feature does not infer a new equilibrium chart from a scalar spectrum. -Source contracts: ``gpec/gpeq.f:114`` selects the derivative; +Source contracts: ``gpec/gpeq.f:131`` selects the derivative; ``gpec/gpdiag.f:1061`` defines the existing perpendicular scalar; ``gpec/gpout.f`` computes it before output transformation. The scalar-only consumer formulas are those derived in rmp_torque c3614334, From ba12ee0b9e4a7b1ce80add0255fb646dd8140f60 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 20 Oct 2026 00:57:03 +0200 Subject: [PATCH 5/5] DOCS - MINOR - State the resonant-source numerical tolerances --- docs/pmodb_target_gauge.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/pmodb_target_gauge.rst b/docs/pmodb_target_gauge.rst index f20d4786..201c5612 100644 --- a/docs/pmodb_target_gauge.rst +++ b/docs/pmodb_target_gauge.rst @@ -41,6 +41,8 @@ real space before the existing output-coordinate transform. Identities are exact; numerical output retains GPEC's finite spectral/spline resolution. A nonzero resonant source is rejected because no periodic solution exists. At a resonant zero source the free homogeneous shift is set to zero. +Numerically, abs(m-nq)<1e-10 is classified as resonant; a source coefficient +below 1e-12 times the largest source coefficient is treated as zero there. The near-resonant result can be large: this feature applies no smoothing. The input pair is ``b_X=b_lag_perp`` and the ordinary perpendicular