Skip to content

PENTRC, GPEC - BUG FIX - Form the pmodb perpendicular displacement from xi^psi and xi^s and read its drive #341

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
logan-nc merged 3 commits into from
Oct 23, 2026
Merged

PENTRC, GPEC - BUG FIX - Form the pmodb perpendicular displacement from xi^psi and xi^s and read its drive #341

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2101e69
PENTRC - BUG FIX - Read the perpendicular pmodb drive and Jacobian pair
krystophny Oct 21, 2026
cbf20c0
GPEC - BUG FIX - Form the pmodb perpendicular displacement from xi^ps…
krystophny Oct 22, 2026
4d7e278
GPEC - CLEANUP - Use the direct xi_perp form in the pmodb diagnostic
logan-nc Oct 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 11 additions & 18 deletions gpec/gpdiag.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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)+
Expand All @@ -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

Expand Down
23 changes: 8 additions & 15 deletions gpec/gpout.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,xsp1_fun
COMPLEX(r8), DIMENSION(:,:), ALLOCATABLE :: eulbparmns,
$ lagbparmns,divxprpmns,curvmns,
$ xmp1mns,xspmns,xmsmns,xmtmns,xmzmns
Expand Down Expand Up @@ -3012,10 +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 iscdftb(mfac,mpert,xvt_fun,mthsurf,xvt_mn)
CALL iscdftb(mfac,mpert,xvz_fun,mthsurf,xvz_mn)

CALL spline_eval(sq,psifac(istep),1)

Expand All @@ -3027,12 +3023,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")
Expand Down Expand Up @@ -3062,10 +3058,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
Expand Down
7 changes: 3 additions & 4 deletions pentrc/inputs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -540,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))') &
Expand Down
Loading