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..201c5612 --- /dev/null +++ b/docs/pmodb_target_gauge.rst @@ -0,0 +1,80 @@ +Optional perturbed-coordinate field strength +============================================= + +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``. 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 +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. +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 +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 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. + +``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: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, +``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..f22a6ffa --- /dev/null +++ b/gen/pmodb_gauge/generate.py @@ -0,0 +1,21 @@ +"""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, 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 +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'+'\n'.join(kernels)+'\nend module pmodb_gauge_kernel\n') +print('PROVED and emitted: complex source, inverse field-line coefficient, scalar correction') 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 a02897d0..4b80a051 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,33 +2924,41 @@ 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 + $ 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 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 REAL(r8), DIMENSION(:), ALLOCATABLE :: psis,ches,chex COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE ::chea,chelagbmns @@ -2965,14 +2974,31 @@ 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)) 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) @@ -3008,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) @@ -3042,9 +3069,17 @@ 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_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* + $ (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 + 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)- @@ -3056,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,:), @@ -3064,6 +3112,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 +3134,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) + ENDDO CALL cspline_dealloc(cspl) @@ -3134,6 +3186,21 @@ 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) ) + 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)) + 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", @@ -3165,6 +3232,12 @@ 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/))) ) + 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)/), + $ (/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( @@ -3382,10 +3455,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) + DEALLOCATE(eulbparmout,lagbparmout,divxprpmout,curvmout, + $ 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..ca907de5 --- /dev/null +++ b/gpec/pmodb_gauge.f90 @@ -0,0 +1,39 @@ +module pmodb_gauge_mod + use local_mod, only: r8 + use ismath_mod, only: iscdftf, iscdftb + 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) + 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), t(0:size(b)-1) + real(r8) :: s, out_r, out_i + integer :: j, nt + nt = size(b)-1 + 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 + 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, 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) + 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 new file mode 100644 index 00000000..8eb51784 --- /dev/null +++ b/gpec/pmodb_gauge_kernel.f90 @@ -0,0 +1,57 @@ +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, Tr, Ti, wr, wi) + use, intrinsic :: iso_fortran_env, only: real64 + implicit none + real(real64), intent(in) :: s, Tr, Ti + real(real64), intent(out) :: wr, wi + real(real64) :: t5, t7, t9, t11 + + 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 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