Skip to content

EQUIL - BUGFIX - Sample every direct-equilibrium surface at the common straight-fieldline angles #329

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
matt-pharr merged 2 commits into from
Oct 23, 2026
Merged

EQUIL - BUGFIX - Sample every direct-equilibrium surface at the common straight-fieldline angles #329

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0576ed3
EQUIL - BUGFIX - Sample every direct-equilibrium surface at the commo…
StuartBenjamin Oct 18, 2026
98c8c10
EQUIL - MINOR - Sample each lsode step in its own subprogram, direct_…
StuartBenjamin Oct 18, 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
117 changes: 84 additions & 33 deletions equil/direct.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,10 @@
c 2. direct_get_bfield.
c 3. direct_position.
c 4. direct_fl_int.
c 5. direct_fl_der.
c 6. direct_refine.
c 7. direct_output.
c 5. direct_fl_sample.
c 6. direct_fl_der.
c 7. direct_refine.
c 8. direct_output.
c-----------------------------------------------------------------------
c subprogram 0. direct_mod.
c module declarations.
Expand Down Expand Up @@ -57,9 +58,10 @@ SUBROUTINE direct_run
REAL(r8), DIMENSION(2, mpsi+1) :: xdx
REAL(r8), DIMENSION(3,3) :: v

REAL(r8) :: xm,dx,rholow,rhohigh
REAL(r8) :: xm,dx,rholow,rhohigh,y4tot,theta
REAL(r8), DIMENSION(4) :: dy
REAL(r8), DIMENSION(:,:), ALLOCATABLE :: ys
TYPE(direct_bfield_type) :: bf
TYPE(spline_type) :: ff
c-----------------------------------------------------------------------
c warning.
c-----------------------------------------------------------------------
Expand Down Expand Up @@ -139,17 +141,6 @@ SUBROUTINE direct_run
DO ipsi=mpsi,0,-1
CALL direct_fl_int(ipsi,y_out,bf)
c-----------------------------------------------------------------------
c fit data to cubic splines.
c-----------------------------------------------------------------------
CALL spline_alloc(ff,istep,4)
ff%xs(0:istep)=y_out(0:istep,4)/y_out(istep,4)
ff%fs(0:istep,1)=y_out(0:istep,2)**2
ff%fs(0:istep,2)=y_out(0:istep,0)/twopi-ff%xs(0:istep)
ff%fs(0:istep,3)=bf%f*
$ (y_out(0:istep,3)-ff%xs(0:istep)*y_out(istep,3))
ff%fs(0:istep,4)=y_out(0:istep,1)/y_out(istep,1)-ff%xs
CALL spline_fit(ff,"periodic")
c-----------------------------------------------------------------------
c allocate space for rzphi and define grids.
c-----------------------------------------------------------------------
IF(ipsi == mpsi)THEN
Expand All @@ -164,15 +155,21 @@ SUBROUTINE direct_run
eqfun%title=(/" b0 "," "," " /)
eqfun%xs=sq%xs
eqfun%ys=(/(itheta,itheta=0,mtheta)/)/REAL(mtheta,r8)
ALLOCATE(ys(0:mtheta,0:4))
ENDIF
c-----------------------------------------------------------------------
c interpolate to uniform grid.
c sample the surface at the common straight-fieldline angles.
c-----------------------------------------------------------------------
y4tot=y_out(istep,4)
CALL direct_fl_int(ipsi,y_out,bf,y4tot,ys)
DO itheta=0,mtheta
CALL spline_eval(ff,rzphi%ys(itheta),1)
rzphi%fs(ipsi,itheta,1:3)=ff%f(1:3)
rzphi%fs(ipsi,itheta,4)=(1+ff%f1(4))
$ *y_out(istep,1)*twopi*psio
theta=rzphi%ys(itheta)
CALL direct_fl_der(4,ys(itheta,0),ys(itheta,1:4),dy)
rzphi%fs(ipsi,itheta,1)=ys(itheta,2)**2
rzphi%fs(ipsi,itheta,2)=ys(itheta,0)/twopi-theta
rzphi%fs(ipsi,itheta,3)=bf%f*
$ (ys(itheta,3)-theta*y_out(istep,3))
rzphi%fs(ipsi,itheta,4)=dy(1)/dy(4)*y4tot*twopi*psio
ENDDO
c-----------------------------------------------------------------------
c store surface quantities.
Expand All @@ -183,7 +180,6 @@ SUBROUTINE direct_run
sq%fs1(ipsi,2)=bf%p1
sq%fs(ipsi,3)=y_out(istep,1)*twopi*psio
sq%fs(ipsi,4)=y_out(istep,3)*bf%f/twopi
CALL spline_dealloc(ff)
ENDDO
c-----------------------------------------------------------------------
c close output files.
Expand Down Expand Up @@ -428,24 +424,28 @@ SUBROUTINE direct_position
END SUBROUTINE direct_position
c-----------------------------------------------------------------------
c subprogram 4. direct_fl_int.
c integrates along field line.
c integrates along field line; with y4tot, a second pass also
c samples y at the straight-fieldline angles rzphi%ys.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c declarations.
c-----------------------------------------------------------------------
SUBROUTINE direct_fl_int(ipsi,y_out,bf)
SUBROUTINE direct_fl_int(ipsi,y_out,bf,y4tot,ys)

INTEGER, INTENT(IN) :: ipsi
REAL(r8), DIMENSION(0:,0:), INTENT(OUT) :: y_out
TYPE(direct_bfield_type), INTENT(OUT) :: bf
REAL(r8), INTENT(IN), OPTIONAL :: y4tot
REAL(r8), DIMENSION(0:,0:), INTENT(OUT), OPTIONAL :: ys

CHARACTER(64) :: message

INTEGER, PARAMETER :: neq=4,liw=30,lrw=22+neq*16
INTEGER :: iopt,istate,itask,itol,jac,mf,ir
INTEGER :: iopt,istate,itask,itol,jac,mf,ir,js
INTEGER, DIMENSION(liw) :: iwork
REAL(r8), PARAMETER :: eps=1e-12
REAL(r8) :: atol,rtol,rfac,deta,r,z,eta,err,psi0,psifac,dr
REAL(r8) :: eta0,y40
REAL(r8), DIMENSION(neq) :: y
REAL(r8), DIMENSION(lrw) :: rwork
c-----------------------------------------------------------------------
Expand Down Expand Up @@ -485,6 +485,8 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf)
deta=twopi/mtheta
y=0
y(2)=SQRT((r-ro)**2+(z-zo)**2)
IF(PRESENT(ys))ys(0,:)=(/eta,y/)
js=1
c-----------------------------------------------------------------------
c set up integrator parameters.
c-----------------------------------------------------------------------
Expand All @@ -502,7 +504,7 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf)
c-----------------------------------------------------------------------
c write header.
c-----------------------------------------------------------------------
IF(out_fl)THEN
IF(out_fl .AND. .NOT. PRESENT(ys))THEN
WRITE(out_2d_unit,10)ipsi,psifac
WRITE(out_2d_unit,20)
ENDIF
Expand All @@ -520,21 +522,26 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf)
c-----------------------------------------------------------------------
c compute and print output for each step.
c-----------------------------------------------------------------------
IF(out_fl)WRITE(out_2d_unit,30)
IF(out_fl .AND. .NOT. PRESENT(ys))WRITE(out_2d_unit,30)
$ istep,eta,rwork(11),y(1:2),r,z,bf%psi,err
IF(bin_fl)WRITE(bin_2d_unit)
IF(bin_fl .AND. .NOT. PRESENT(ys))WRITE(bin_2d_unit)
$ REAL(eta,4),REAL(y(1:2),4),REAL(r,4),REAL(z,4),REAL(err,4)
c-----------------------------------------------------------------------
c advance differential equations.
c-----------------------------------------------------------------------
IF(eta >= twopi .OR. istep >= enstep .OR. istate < 0
$ .OR. ABS(err) >= 1)EXIT
istep=istep+1
eta0=eta
y40=y(4)
CALL lsode(direct_fl_der,neq,y,eta,twopi,itol,rtol,atol,
$ itask,istate,iopt,rwork,lrw,iwork,liw,jac,mf)
IF(PRESENT(ys))CALL direct_fl_sample(eta0,eta,y40,y(4),
$ y4tot,rwork(21),js,ys)
ENDDO
IF(out_fl)WRITE(out_2d_unit,20)
IF(bin_fl)WRITE(bin_2d_unit)
IF(PRESENT(ys))ys(mtheta,:)=(/eta,y/)
IF(out_fl .AND. .NOT. PRESENT(ys))WRITE(out_2d_unit,20)
IF(bin_fl .AND. .NOT. PRESENT(ys))WRITE(bin_2d_unit)
c-----------------------------------------------------------------------
c abort if istep > enstep.
c-----------------------------------------------------------------------
Expand All @@ -549,7 +556,51 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf)
RETURN
END SUBROUTINE direct_fl_int
c-----------------------------------------------------------------------
c subprogram 5. direct_fl_der.
c subprogram 5. direct_fl_sample.
c samples the last lsode step, eta0 to eta1, at the angles
c rzphi%ys it covers: Newton on lsode's interpolant for
c y(4) = theta*y4tot, then y there.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c declarations.
c-----------------------------------------------------------------------
SUBROUTINE direct_fl_sample(eta0,eta1,y40,y41,y4tot,yh,js,ys)

REAL(r8), INTENT(IN) :: eta0,eta1,y40,y41,y4tot
REAL(r8), DIMENSION(*), INTENT(IN) :: yh
INTEGER, INTENT(INOUT) :: js
REAL(r8), DIMENSION(0:,0:), INTENT(INOUT) :: ys

INTEGER, PARAMETER :: neq=4
REAL(r8), PARAMETER :: eps=1e-12
INTEGER :: it,iflag
REAL(r8) :: y4j,eta,res
REAL(r8), DIMENSION(neq) :: dky
c-----------------------------------------------------------------------
c find and sample each angle in the step.
c-----------------------------------------------------------------------
DO WHILE(js < mtheta)
y4j=rzphi%ys(js)*y4tot
IF(y4j > y41)EXIT
eta=eta0+(eta1-eta0)*(y4j-y40)/(y41-y40)
DO it=1,20
CALL dintdy(eta,0,yh,neq,dky,iflag)
res=dky(4)-y4j
CALL dintdy(eta,1,yh,neq,dky,iflag)
eta=MIN(MAX(eta-res/dky(4),eta0),eta1)
IF(ABS(res) <= eps*y4tot)EXIT
ENDDO
CALL dintdy(eta,0,yh,neq,dky,iflag)
ys(js,:)=(/eta,dky/)
js=js+1
ENDDO
c-----------------------------------------------------------------------
c terminate.
c-----------------------------------------------------------------------
RETURN
END SUBROUTINE direct_fl_sample
c-----------------------------------------------------------------------
c subprogram 6. direct_fl_der.
c contains differential equations for field line averages.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
Expand Down Expand Up @@ -590,7 +641,7 @@ SUBROUTINE direct_fl_der(neq,eta,y,dy)
RETURN
END SUBROUTINE direct_fl_der
c-----------------------------------------------------------------------
c subprogram 6. direct_refine.
c subprogram 7. direct_refine.
c moves a point orthogonally to a specified flux surface.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
Expand Down Expand Up @@ -638,7 +689,7 @@ SUBROUTINE direct_refine(rfac,eta,psi0)
RETURN
END SUBROUTINE direct_refine
c-----------------------------------------------------------------------
c subprogram 7. direct_output.
c subprogram 8. direct_output.
c diagnoses input.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
Expand Down
Loading