From 8a57848295431a5265320299631a7b3a787c1e26 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Mon, 12 Oct 2026 15:45:21 -0700 Subject: [PATCH 1/4] EQUIL - BUGFIX - ldp_i reads |f| (as efit does); document the i-file format With f < 0 (e.g. a TokaMaker equilibrium with B_t < 0) inverse_run built a negative q and stopped ('Invalid extrapolation near axis'). Flipping the sign of f in an i-file now gives identical results. Co-Authored-By: Claude Opus 5.5 --- equil/read_eq.f | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/equil/read_eq.f b/equil/read_eq.f index 3ad187791..e042ed87f 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -801,7 +801,13 @@ SUBROUTINE read_eq_ldp_d END SUBROUTINE read_eq_ldp_d c----------------------------------------------------------------------- c subprogram 12. read_eq_ldp_i. -c reads data from L. Don Pearlstein's inverse equilibrium code. +c reads data from L. Don Pearlstein's inverse equilibrium code, also +c written by OpenFUSIONToolkit TokaMaker (save_ifile). Sequential +c unformatted records, real*8, axis first, theta index fastest: +c mx,my (integer*4: number of surfaces, number of angles) +c psi(mx) [Wb/rad], f(mx) = R*B_t [T m], p(mx) [Pa], q(mx) +c r(my,mx), z(my,mx) [m], periodic point duplicated +c optional: f*df/dpsi(mx), dp/dpsi(mx) (not read here) c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. @@ -840,7 +846,7 @@ SUBROUTINE read_eq_ldp_i psio=psi(mx)-psi(0) CALL spline_alloc(sq_in,mx,4) sq_in%xs=(psi-psi(0))/psio - sq_in%fs(:,1)=f + sq_in%fs(:,1)=ABS(f) sq_in%fs(:,2)=p*mu0 sq_in%fs(:,3)=q if(psio.lt.0)psio=-psio From 608e0942022bd95fb1bae11ef493ec970d8bbfdc Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Fri, 9 Oct 2026 21:49:00 -0700 Subject: [PATCH 2/4] EQUIL - FEATURE - ldp_i reads an optional psi bounds record [psi_axis, psi_lcfs] TokaMaker's save_ifile pads the last i-file surface inward from the LCFS. An optional final record of length 2 holding the unpadded plasma bounds now normalizes psi: xs = (psi - psi_axis)/(psi_lcfs - psi_axis), psio = psi_lcfs - psi_axis, so psihigh is the true psi_N; psihigh above the last file surface is clamped to it with a note. The record may follow the optional FF', P' records (skipped here) and is told apart by length: a length-mx read of it fails, and the reader then rewinds and skips the records already read (the position after a failed read is processor dependent, so no BACKSPACE). Without the record results are bit-identical; older builds ignore it. Co-Authored-By: Claude Opus 5.5 --- equil/read_eq.f | 50 +++++++++++++++++++++++++++++++++++++++++-------- input/equil.in | 2 +- 2 files changed, 43 insertions(+), 9 deletions(-) diff --git a/equil/read_eq.f b/equil/read_eq.f index e042ed87f..10ba97093 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -807,15 +807,20 @@ END SUBROUTINE read_eq_ldp_d c mx,my (integer*4: number of surfaces, number of angles) c psi(mx) [Wb/rad], f(mx) = R*B_t [T m], p(mx) [Pa], q(mx) c r(my,mx), z(my,mx) [m], periodic point duplicated -c optional: f*df/dpsi(mx), dp/dpsi(mx) (not read here) +c optional, each only if present, in this order: +c f*df/dpsi(mx), dp/dpsi(mx) (skipped here); +c psi_axis, psi_lcfs (2 reals): unpadded plasma bounds, which +c normalize psi; the last surface may lie inside the lcfs and +c caps psihigh. Optional records are told apart by length (mx, 2). c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. c----------------------------------------------------------------------- SUBROUTINE read_eq_ldp_i - INTEGER :: mx,my - REAL(r8), DIMENSION(:), POINTER :: psi,f,p,q + INTEGER :: mx,my,irec,ios1,ios2,ios3 + REAL(r8), DIMENSION(2) :: psib + REAL(r8), DIMENSION(:), POINTER :: psi,f,p,q,buf REAL(r8), DIMENSION(:,:), POINTER :: r,z c----------------------------------------------------------------------- c open data file, read scalar data. @@ -829,9 +834,12 @@ SUBROUTINE read_eq_ldp_i my=my-1 mx=mx-1 ALLOCATE(psi(0:mx),f(0:mx),p(0:mx),q(0:mx), - $ r(0:my,0:mx),z(0:my,0:mx)) + $ r(0:my,0:mx),z(0:my,0:mx),buf(0:mx)) c----------------------------------------------------------------------- -c read binary data. +c read binary data. The psi bounds record follows z, ffp or pp. A +c length-mx read of it fails (short record); the file position is +c then processor dependent, so rewind and skip the 7 fixed records +c and any ffp, pp records before reading it. c----------------------------------------------------------------------- READ(in_unit)psi READ(in_unit)f @@ -839,13 +847,39 @@ SUBROUTINE read_eq_ldp_i READ(in_unit)q READ(in_unit)r READ(in_unit)z + READ(in_unit,IOSTAT=ios1)buf + ios2=ios1 + IF(ios1 == 0)READ(in_unit,IOSTAT=ios2)buf + ios3=ios2 + IF(ios2 == 0)READ(in_unit,IOSTAT=ios3)psib + IF(ios2 > 0)THEN + REWIND(in_unit) + DO irec=1,7+COUNT((/ios1,ios2/) == 0) + READ(in_unit) + ENDDO + READ(in_unit,IOSTAT=ios3)psib + ENDIF CALL bin_close(in_unit) c----------------------------------------------------------------------- c copy and revise 1D arrays. c----------------------------------------------------------------------- - psio=psi(mx)-psi(0) CALL spline_alloc(sq_in,mx,4) - sq_in%xs=(psi-psi(0))/psio + IF(ios3 == 0)THEN + psio=psib(2)-psib(1) + sq_in%xs=(psi-psib(1))/psio + IF(.NOT.(sq_in%xs(mx) > 0 .AND. sq_in%xs(mx) <= 1)) + $ CALL program_stop("ldp_i psi bounds disagree with psi") + IF(verbose)WRITE(*,'(1x,a,1p,2es14.6)') + $ "ldp_i psi bounds: psi_axis, psi_lcfs =",psib + IF(psihigh > sq_in%xs(mx))THEN + WRITE(*,'(1x,2(a,es12.5))')"Note: psihigh =",psihigh, + $ " clamped to the last i-file surface",sq_in%xs(mx) + psihigh=sq_in%xs(mx) + ENDIF + ELSE + psio=psi(mx)-psi(0) + sq_in%xs=(psi-psi(0))/psio + ENDIF sq_in%fs(:,1)=ABS(f) sq_in%fs(:,2)=p*mu0 sq_in%fs(:,3)=q @@ -861,7 +895,7 @@ SUBROUTINE read_eq_ldp_i c----------------------------------------------------------------------- c process inverse equilibrium and deallocate local arrays. c----------------------------------------------------------------------- - DEALLOCATE(psi,f,p,q,r,z) + DEALLOCATE(psi,f,p,q,r,z,buf) CALL inverse_run c----------------------------------------------------------------------- c terminate. diff --git a/input/equil.in b/input/equil.in index 98cbb3a91..7cbce739c 100644 --- a/input/equil.in +++ b/input/equil.in @@ -9,7 +9,7 @@ grid_type="ldp" ! Radial grid packing of equilibrium quantities. Accepts rho, ldp, pow1, pow2, or original. ldp packs points near the core and edge. pow* packs near the edge. psilow=1e-2 ! Minimum value of psi, normalized from 0 to 1 - psihigh=0.994 ! Maximum value of psi, normalized from 0 to 1 + psihigh=0.994 ! Maximum value of psi, normalized from 0 to 1 (ldp_i with psi bounds: capped at the last file surface) mpsi=128 ! Number of radial grid intervals for equilibrium quantities mtheta=256 ! Number of equally spaced poloidal grid intervals for all splines newq0=0 ! Grad-Shafranov solution invariant adjustment of the q profile to give the specified value of q at the axis. Default 0 uses input file value. From 82eca30ebb2a14732ab2ecf3a0a9397f376c232a Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 13 Oct 2026 01:36:42 -0700 Subject: [PATCH 3/4] EQUIL - MINOR - Plainer read_eq_ldp_i comments Co-Authored-By: Claude Opus 5.5 --- equil/read_eq.f | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/equil/read_eq.f b/equil/read_eq.f index 10ba97093..aeefd72fd 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -809,9 +809,9 @@ END SUBROUTINE read_eq_ldp_d c r(my,mx), z(my,mx) [m], periodic point duplicated c optional, each only if present, in this order: c f*df/dpsi(mx), dp/dpsi(mx) (skipped here); -c psi_axis, psi_lcfs (2 reals): unpadded plasma bounds, which -c normalize psi; the last surface may lie inside the lcfs and -c caps psihigh. Optional records are told apart by length (mx, 2). +c psi_axis, psi_lcfs (2 reals): psi is normalized to these; the +c last surface may lie inside the lcfs and then caps psihigh. +c Optional records are told apart by length (mx, 2). c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. @@ -836,10 +836,10 @@ SUBROUTINE read_eq_ldp_i ALLOCATE(psi(0:mx),f(0:mx),p(0:mx),q(0:mx), $ r(0:my,0:mx),z(0:my,0:mx),buf(0:mx)) c----------------------------------------------------------------------- -c read binary data. The psi bounds record follows z, ffp or pp. A -c length-mx read of it fails (short record); the file position is -c then processor dependent, so rewind and skip the 7 fixed records -c and any ffp, pp records before reading it. +c read binary data. The psi bounds record follows z, ffp or pp. +c Reading it into a profile array fails (short record); the file +c position is then processor dependent, so rewind and skip the 7 +c fixed records and any ffp, pp records before reading it. c----------------------------------------------------------------------- READ(in_unit)psi READ(in_unit)f From 8dc2d0daf56ae0d2f467a2872f638e54fde8852e Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 15 Oct 2026 14:31:38 +0200 Subject: [PATCH 4/4] EQUIL - MINOR - check sign of i-file psi --- equil/read_eq.f | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/equil/read_eq.f b/equil/read_eq.f index fddc3144d..7158aee63 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -876,8 +876,11 @@ SUBROUTINE read_eq_ldp_i IF(ios3 == 0)THEN psio=psib(2)-psib(1) sq_in%xs=(psi-psib(1))/psio - IF(.NOT.(sq_in%xs(mx) > 0 .AND. sq_in%xs(mx) <= 1)) - $ CALL program_stop("ldp_i psi bounds disagree with psi") + IF(.NOT.(MINVAL(sq_in%xs) >= 0 .AND. sq_in%xs(mx) <= 1))THEN + WRITE(*,'(1x,a,1p,4es14.6)')"psi_axis, psi_lcfs, psi(0), " + $ //"psi(mx) =",psib,psi(0),psi(mx) + CALL program_stop("ldp_i psi bounds disagree with psi") + ENDIF IF(verbose)WRITE(*,'(1x,a,1p,2es14.6)') $ "ldp_i psi bounds: psi_axis, psi_lcfs =",psib IF(psihigh > sq_in%xs(mx))THEN