From 2101e69d842324a20363ffa561848996c1057ea3 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 21 Oct 2026 15:15:27 +0200 Subject: [PATCH 1/3] PENTRC - BUG FIX - Read the perpendicular pmodb drive and Jacobian pair --- pentrc/inputs.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index e98d64a9..a3a550b1 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -434,9 +434,8 @@ subroutine read_pmodb(file,jac_in,jsurf_in,tmag_in,debug,op_powin) +xj*reshape(table(:,10),(/npsi,nm/),order=(/1,2/)) endif - ! For consistency with GPEC-0.3.0 and smaller values near rationals - if(verbose) print *," -> Using div(xi_perp) = -(dB/B+kappa.xi_perp)" - divxmni = -(lagbmni+kapxmni) + ! Keep the supplied divergence and use its perpendicular gauge. + lagbmni = -(divxmni+kapxmni) ! convert to chebyshev coordinates if(jac_in=="" .or. jac_in=="default")then From cbf20c0a6009d4ff790f56f09acdb9379d43b247 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Thu, 22 Oct 2026 07:17:11 +0200 Subject: [PATCH 2/3] GPEC - BUG FIX - Form the pmodb perpendicular displacement from xi^psi and xi^s gpout_pmodb built xi_perp as xi minus the Fourier-truncated 1/(m-nq) closure xi_par, which left large residues in Bdivxprp and Bkxprp at the resonant and window-edge harmonics. Build J*xi_perp^theta and J*xi_perp^zeta directly from xi^psi and xi^s with the eqfun(2:3) metric combinations, as in the smat..zmat kinetic matrices. Output format is unchanged. --- gpec/gpout.f | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/gpec/gpout.f b/gpec/gpout.f index ce376b52..8faae6d0 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -2955,7 +2955,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) COMPLEX(r8), DIMENSION(mpert) :: eulbpar_mn,lagbpar_mn, $ divxprp_mn,curv_mn COMPLEX(r8), DIMENSION(0:mthsurf) :: xsp_fun,xms_fun, - $ bvt_fun,bvz_fun,xmz_fun,xvt_fun,xvz_fun,xmt_fun,xsp1_fun + $ bvt_fun,bvz_fun,xmz_fun,xmt_fun,xsp1_fun COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparmns, $ lagbparmns,divxprpmns,curvmns, $ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns @@ -3014,8 +3014,6 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) CALL iscdftb(mfac,mpert,bvz_fun,mthsurf,bvz_mn) CALL iscdftb(mfac,mpert,xmt_fun,mthsurf,xmt_mn) CALL iscdftb(mfac,mpert,xmz_fun,mthsurf,xmz_mn) - CALL iscdftb(mfac,mpert,xvt_fun,mthsurf,xvt_mn) - CALL iscdftb(mfac,mpert,xvz_fun,mthsurf,xvz_mn) CALL spline_eval(sq,psifac(istep),1) @@ -3027,12 +3025,12 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) CALL bicube_eval(rzphi,psifac(istep),theta(itheta),0) jac=rzphi%f(4) - cspl%fs(itheta,1)=xmt_fun(itheta)- - $ (chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) - cspl%fs(itheta,2)=xmz_fun(itheta)- - $ sq%f(4)*(chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) +c Form J*xi_perp directly from xi^psi and xi^alpha=xi^s/chi1. + cspl%fs(itheta,1)=chi1**2* + $ (-eqfun%f(2)*xsp_fun(itheta)+ + $ eqfun%f(3)*xms_fun(itheta)/chi1) + cspl%fs(itheta,2)=sq%f(4)*cspl%fs(itheta,1)- + $ jac*xms_fun(itheta)/chi1 ENDDO CALL cspline_fit(cspl,"periodic") @@ -3062,10 +3060,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) $ divxprpfun(istep,itheta) curvfun(istep,itheta)=-xsp_fun(itheta)/eqfun%f(1)*sq%f1(2)- $ (xsp_fun(itheta)*eqfun%fx(1)+eqfun%fy(1)* - $ (xms_fun(itheta)/(chi1*sq%f(4))+ - $ xmz_fun(itheta)/(jac*sq%f(4))- - $ (chi1/(jac*eqfun%f(1)))**2* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)))) + $ cspl%f(1)/jac) equilbfun(istep,itheta) = eqfun%f(1) dphi(itheta) = rzphi%f(3) ENDDO From 4d7e2781fcafd78a5899a285cd528f4ba9767638 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 22 Oct 2026 20:39:10 -0400 Subject: [PATCH 3/3] GPEC - CLEANUP - Use the direct xi_perp form in the pmodb diagnostic gpdiag_pmodb now builds J*xi_perp from xi^psi and xi^s like gpout_pmodb, for the divergence, curvature and llagb columns, so the diagnostic agrees with the main pmodb output. Its llagb vs -(Bdivx+Bkx) residual at resonant and window-edge m drops from O(1-100) to <1e-3 on the Solov'ev examples. Drops the now-unused xmt/xmz (gpout_pmodb) and xmz/xvt/xvz (gpdiag_pmodb) inverse transforms and relabels the PENTRC pmodb debug columns. gpec_pmodb, gpec_pmodb_fun, pmodbmn diagnostic and PENTRC outputs are bit-identical to cbf20c0a. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01VMDGyfQKRZQMdk3NJXsKQE --- gpec/gpdiag.f | 29 +++++++++++------------------ gpec/gpout.f | 4 +--- pentrc/inputs.f90 | 2 +- 3 files changed, 13 insertions(+), 22 deletions(-) diff --git a/gpec/gpdiag.f b/gpec/gpdiag.f index e029e89f..9f1c7e3c 100644 --- a/gpec/gpdiag.f +++ b/gpec/gpdiag.f @@ -936,7 +936,7 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) COMPLEX(r8), DIMENSION(mpert) :: eulbpar_mn,lagbpar_mn, $ llagbpar_mn,divx_mn,curv_mn,divx2_mn,curv2_mn COMPLEX(r8), DIMENSION(0:mthsurf) :: xsp_fun,xms_fun, - $ bvt_fun,bvz_fun,xmt_fun,xmz_fun,xvt_fun,xvz_fun,xmp1_fun + $ bvt_fun,bvz_fun,xmt_fun,xmp1_fun REAL(r8), DIMENSION(:,:), POINTER :: rs,zs,eqfunx,eqfuny,eqfuns COMPLEX(r8), DIMENSION(:,:), POINTER :: eulbparmns,lagbparmns, @@ -989,9 +989,6 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) CALL iscdftb(mfac,mpert,bvt_fun,mthsurf,bvt_mn) CALL iscdftb(mfac,mpert,bvz_fun,mthsurf,bvz_mn) CALL iscdftb(mfac,mpert,xmt_fun,mthsurf,xmt_mn) - CALL iscdftb(mfac,mpert,xmz_fun,mthsurf,xmz_mn) - CALL iscdftb(mfac,mpert,xvt_fun,mthsurf,xvt_mn) - CALL iscdftb(mfac,mpert,xvz_fun,mthsurf,xvz_mn) CALL iscdftb(mfac,mpert,xmp1_fun,mthsurf,xmp1_mn) xspmns(istep,:)=xsp_mn @@ -1010,12 +1007,12 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) CALL bicube_eval(rzphi,psifac(istep),theta(itheta),0) jac=rzphi%f(4) - cspl%fs(itheta,1)=xmt_fun(itheta)- - $ (chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) - cspl%fs(itheta,2)=xmz_fun(itheta)- - $ sq%f(4)*(chi1/eqfun%f(1))**2/jac* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta)) +c Form J*xi_perp directly from xi^psi and xi^alpha=xi^s/chi1. + cspl%fs(itheta,1)=chi1**2* + $ (-eqfun%f(2)*xsp_fun(itheta)+ + $ eqfun%f(3)*xms_fun(itheta)/chi1) + cspl%fs(itheta,2)=sq%f(4)*cspl%fs(itheta,1)- + $ jac*xms_fun(itheta)/chi1 cspl%fs(itheta,3)=xmt_fun(itheta) ENDDO @@ -1060,10 +1057,9 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) $ xsp_fun(itheta)*eqfun%fx(1)+ $ xms_fun(itheta)/(chi1*sq%f(4))*eqfun%fy(1) llagbparfun(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))) + $ eulbparfun(istep,itheta)+ + $ xsp_fun(itheta)*eqfun%fx(1)+ + $ eqfun%fy(1)*cspl%f(1)/jac divxfun(istep,itheta)=xmp1_fun(itheta)+ $ (jac1/jac)*xsp_fun(itheta)+ @@ -1074,10 +1070,7 @@ SUBROUTINE gpdiag_pmodb(egnum,xspmn) curvfun(istep,itheta)=-xsp_fun(itheta)/eqfun%f(1)*sq%f1(2)- $ xsp_fun(itheta)*eqfun%fx(1)-eqfun%fy(1)* - $ (xms_fun(itheta)/(chi1*sq%f(4))+ - $ xmz_fun(itheta)/(jac*sq%f(4))- - $ (chi1/(jac*eqfun%f(1)))**2* - $ (xvt_fun(itheta)+sq%f(4)*xvz_fun(itheta))) + $ cspl%f(1)/jac curv2fun(istep,itheta)=-jac*eqfun%f(1)*curvfun(istep,itheta) ENDDO diff --git a/gpec/gpout.f b/gpec/gpout.f index 8faae6d0..069ae193 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -2955,7 +2955,7 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) COMPLEX(r8), DIMENSION(mpert) :: eulbpar_mn,lagbpar_mn, $ divxprp_mn,curv_mn COMPLEX(r8), DIMENSION(0:mthsurf) :: xsp_fun,xms_fun, - $ bvt_fun,bvz_fun,xmz_fun,xmt_fun,xsp1_fun + $ bvt_fun,bvz_fun,xsp1_fun COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparmns, $ lagbparmns,divxprpmns,curvmns, $ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns @@ -3012,8 +3012,6 @@ SUBROUTINE gpout_pmodb(egnum,xspmn) CALL iscdftb(mfac,mpert,xms_fun,mthsurf,xms_mn) CALL iscdftb(mfac,mpert,bvt_fun,mthsurf,bvt_mn) CALL iscdftb(mfac,mpert,bvz_fun,mthsurf,bvz_mn) - CALL iscdftb(mfac,mpert,xmt_fun,mthsurf,xmt_mn) - CALL iscdftb(mfac,mpert,xmz_fun,mthsurf,xmz_mn) CALL spline_eval(sq,psifac(istep),1) diff --git a/pentrc/inputs.f90 b/pentrc/inputs.f90 index a3a550b1..d764fffd 100644 --- a/pentrc/inputs.f90 +++ b/pentrc/inputs.f90 @@ -539,7 +539,7 @@ subroutine read_pmodb(file,jac_in,jsurf_in,tmag_in,debug,op_powin) write(out_unit,'(a8,es16.8e3,2/)') " chi1 = ",chi1 write(out_unit,'(1x,a16,a5,8(1x,a16))')"psi_n","m",& 'real(deltaB/B)','imag(deltaB/B)','real(divxprp)','imag(divxprp)',& - 'real(deltaB)','imag(deltaB)','real(Bdivxprp)','imag(Bdivxprp)' + 'real(deltaBprp)','imag(deltaBprp)','real(Bdivxprp)','imag(Bdivxprp)' do i=0,npsi-1 do j=1,mpert write(out_unit,'(1x,es16.8e3,i5,14(1x,es16.8e3))') &