From f8bd9a943abbdaf7be3ee05809c7f5715e071f91 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 21 Jun 2026 21:12:44 -0700 Subject: [PATCH 01/22] Zeff profile read-in and spline evaluation --- rdcon/dcon.f | 3 ++- rdcon/dcon_mod.f | 4 +++- rdcon/mercier.f | 15 ++++++++++++--- 3 files changed, 17 insertions(+), 5 deletions(-) diff --git a/rdcon/dcon.f b/rdcon/dcon.f index a8e056246..32fdf0125 100644 --- a/rdcon/dcon.f +++ b/rdcon/dcon.f @@ -226,7 +226,8 @@ SUBROUTINE dcon_run $ tol_nr,tol_r,crossover,ucrit,singfac_min,singfac_max, $ cyl_flag,dmlim,lim_flag,sas_flag,sing_order,sort_type, $ gal_flag,regrid_flag,sing1_flag,qlow,qhigh, - $ sing_order_ceiling,degen_tol,coil,Zeff,reform_eq_with_psilim + $ sing_order_ceiling,degen_tol,coil,reform_eq_with_psilim, + $ Zeff, psi_N_Zeff NAMELIST/rdcon_output/interp,crit_break,out_bal1, $ bin_bal1,out_bal2,bin_bal2,out_metric,bin_metric,out_fmat, $ bin_fmat,out_gmat,bin_gmat,out_kmat,bin_kmat,out_sol, diff --git a/rdcon/dcon_mod.f b/rdcon/dcon_mod.f index 99bf7444c..02c9ce502 100644 --- a/rdcon/dcon_mod.f +++ b/rdcon/dcon_mod.f @@ -80,7 +80,9 @@ MODULE rdcon_mod $ sing_start=0,nzero=0 REAL(r8) :: thmax0=1,ucrit=1e4,tol_r=1e-5,tol_nr=1e-5, $ crossover=1e-2,mthsurf0=1,prefac=1.0, - $ plasma1=0.0,vacuum1=0.0,total1=0.0,Zeff=1.52 + $ plasma1=0.0,vacuum1=0.0,total1=0.0 + + REAL(r8), DIMENSION(:), ALLOCATABLE :: Zeff, psi_N_Zeff REAL(r8) :: dpsi_intvl=0.1,dpsi1_intvl=0.1 TYPE(spline_type) :: locstab, mreterms diff --git a/rdcon/mercier.f b/rdcon/mercier.f index f86faf262..4bcf0b109 100644 --- a/rdcon/mercier.f +++ b/rdcon/mercier.f @@ -37,12 +37,20 @@ SUBROUTINE mercier_scan $ eps_eff REAL(r8), DIMENSION(:), POINTER :: avg REAL(r8), DIMENSION(2) :: rmax_loc, rmin_loc, zmax_loc, zmin_loc - TYPE(spline_type), TARGET :: ff + TYPE(spline_type), TARGET :: ff, Zeff_spline c----------------------------------------------------------------------- c prepare spline types. c----------------------------------------------------------------------- IF(MRE_flag)THEN CALL spline_alloc(ff,mtheta,22) + IF(ALLOCATED(Zeff) .AND. ALLOCATED(psi_N_Zeff))THEN + CALL spline_alloc(Zeff_spline,SIZE(psi_N_),1) + Zeff_spline%xs=psi_N_Zeff + Zeff_spline%fs(:,1)=Zeff + CALL spline_fit(Zeff_spline,"extrap") + ELSE + CALL program_stop("Provide Zeff profile values") + ENDIF ELSE CALL spline_alloc(ff,mtheta,5) ENDIF @@ -207,8 +215,9 @@ SUBROUTINE mercier_scan c----------------------------------------------------------------------- ! mu_e_on_nu_e=(ftr/(1.d0-ftr))*(1.d0+0.533d0/Zeff) ! ^ Taking banana limit of eq. B17 (& B14) in Callen - mufrac=ftr*(1.d0+0.533d0/Zeff)/ - $ ((1.d0-ftr)+ftr*(1.d0+0.533d0/Zeff)) + CALL spline_eval(Zeff_spline,psifac,0) + mufrac=ftr*(1.d0+0.533d0/Zeff_spline%f(1))/ + $ ((1.d0-ftr)+ftr*(1.d0+0.533d0/Zeff_spline%f(1))) c----------------------------------------------------------------------- c evaluate geometric prefactors of MRE stability terms from c Hegna 1999 https://doi.org/10.1063/1.873661 From ccb52e6d74552a86d05e0b13e07fd13127f4bdbf Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 21 Jun 2026 21:12:52 -0700 Subject: [PATCH 02/22] Zeff profile printing --- rdcon/rdcon_netcdf.f | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/rdcon/rdcon_netcdf.f b/rdcon/rdcon_netcdf.f index 1b2c52d34..0534e7a10 100644 --- a/rdcon/rdcon_netcdf.f +++ b/rdcon/rdcon_netcdf.f @@ -130,7 +130,11 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et) CALL check( nf90_put_att(ncid,nf90_global,'betat', betat)) CALL check( nf90_put_att(ncid,nf90_global,'betan', betan)) CALL check( nf90_put_att(ncid,nf90_global,'volume', volume)) - CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff)) + IF(ALLOCATED(Zeff) .AND. ALLOCATED(psi_N_Zeff))THEN + CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff)) + CALL check( nf90_put_att(ncid,nf90_global,'psi_N_Zeff', + $ psi_N_Zeff)) + ENDIF CALL check( nf90_put_att(ncid,nf90_global,'bwall', bwall)) CALL check( nf90_put_att(ncid,nf90_global,'bt0', bt0)) CALL check( nf90_put_att(ncid,nf90_global,'q0', q0)) From 182190e5ceb1b0f017944d3c08a55523a5d04732 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 21 Jun 2026 21:20:11 -0700 Subject: [PATCH 03/22] bugfix --- rdcon/mercier.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rdcon/mercier.f b/rdcon/mercier.f index 4bcf0b109..bd4713a42 100644 --- a/rdcon/mercier.f +++ b/rdcon/mercier.f @@ -44,7 +44,7 @@ SUBROUTINE mercier_scan IF(MRE_flag)THEN CALL spline_alloc(ff,mtheta,22) IF(ALLOCATED(Zeff) .AND. ALLOCATED(psi_N_Zeff))THEN - CALL spline_alloc(Zeff_spline,SIZE(psi_N_),1) + CALL spline_alloc(Zeff_spline,SIZE(psi_N_Zeff),1) Zeff_spline%xs=psi_N_Zeff Zeff_spline%fs(:,1)=Zeff CALL spline_fit(Zeff_spline,"extrap") From 9fe031d4ea96758c6cf7e0d042043b37ffbbc30a Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 21 Jun 2026 22:53:42 -0700 Subject: [PATCH 04/22] Adding 'read_var_len' subroutine to read in variable length arrays --- equil/local.f | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/equil/local.f b/equil/local.f index b4e25380d..fb7aba804 100644 --- a/equil/local.f +++ b/equil/local.f @@ -13,6 +13,7 @@ c 5. ascii_close. c 6. program_stop. c 7. floored_log. +c 8. read_var_len. c----------------------------------------------------------------------- c subprogram 0. local. c module declarations. @@ -259,4 +260,49 @@ FUNCTION floored_log(u) RESULT(ulog) c----------------------------------------------------------------------- RETURN END FUNCTION floored_log +c----------------------------------------------------------------------- +c subprogram 8. read_var_len. +c allocates variable length arrays from set-length read-in arrays. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE read_var_len(x_in, y_in, x_out, y_out) + + IMPLICIT NONE + + REAL(r8), DIMENSION(:), INTENT(IN) :: x_in, y_in + REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: x_out, y_out + INTEGER :: n, i + REAL(r8) :: default_value_x,default_value_y +c----------------------------------------------------------------------- +c assume last values of x_in, y_in are sentinals +c----------------------------------------------------------------------- + default_value_x = x_in(SIZE(x_in)) + default_value_y = y_in(SIZE(y_in)) +c----------------------------------------------------------------------- +c count leading non-sentinel values in x_in. +c----------------------------------------------------------------------- + n = 0 + DO i = 1, SIZE(x_in) - 1 + IF (x_in(i) /= default_value_x) THEN + n = n + 1 + ELSE + EXIT + END IF + END DO +c----------------------------------------------------------------------- +c allocate output arrays and copy if values exist. +c----------------------------------------------------------------------- + IF (n > 0) THEN + ALLOCATE(x_out(n)) + ALLOCATE(y_out(n)) + x_out(1:n) = x_in(1:n) + y_out(1:n) = y_in(1:n) + END IF +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE read_var_len END MODULE local_mod From 4e1c8e23b1e58137d7b921477e81f5dd3e79b2ce Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 21 Jun 2026 22:54:13 -0700 Subject: [PATCH 05/22] Using read_var_len to allocate internal variables Zeff_, psi_N_Zeff_ --- rdcon/dcon.f | 5 ++++- rdcon/dcon_mod.f | 3 ++- rdcon/mercier.f | 8 ++++---- rdcon/rdcon_netcdf.f | 6 +++--- 4 files changed, 13 insertions(+), 9 deletions(-) diff --git a/rdcon/dcon.f b/rdcon/dcon.f index 32fdf0125..8852772bb 100644 --- a/rdcon/dcon.f +++ b/rdcon/dcon.f @@ -261,7 +261,10 @@ SUBROUTINE dcon_run CALL gal_read_input ENDIF CALL ascii_close(in_unit) - +c----------------------------------------------------------------------- +c Allocate variable-length Zeff input profile. +c----------------------------------------------------------------------- + CALL read_var_len(Zeff,psi_N_Zeff,Zeff_,psi_N_Zeff_) c----------------------------------------------------------------------- c open output files, read, process, and diagnose equilibrium. c----------------------------------------------------------------------- diff --git a/rdcon/dcon_mod.f b/rdcon/dcon_mod.f index 02c9ce502..47604f187 100644 --- a/rdcon/dcon_mod.f +++ b/rdcon/dcon_mod.f @@ -82,7 +82,8 @@ MODULE rdcon_mod $ crossover=1e-2,mthsurf0=1,prefac=1.0, $ plasma1=0.0,vacuum1=0.0,total1=0.0 - REAL(r8), DIMENSION(:), ALLOCATABLE :: Zeff, psi_N_Zeff + REAL(r8), DIMENSION(:), ALLOCATABLE :: Zeff_, psi_N_Zeff_ + REAL(r8), DIMENSION(1000) :: Zeff=-1.0, psi_N_Zeff=-1.0 ! read-in REAL(r8) :: dpsi_intvl=0.1,dpsi1_intvl=0.1 TYPE(spline_type) :: locstab, mreterms diff --git a/rdcon/mercier.f b/rdcon/mercier.f index bd4713a42..cff17260b 100644 --- a/rdcon/mercier.f +++ b/rdcon/mercier.f @@ -43,10 +43,10 @@ SUBROUTINE mercier_scan c----------------------------------------------------------------------- IF(MRE_flag)THEN CALL spline_alloc(ff,mtheta,22) - IF(ALLOCATED(Zeff) .AND. ALLOCATED(psi_N_Zeff))THEN - CALL spline_alloc(Zeff_spline,SIZE(psi_N_Zeff),1) - Zeff_spline%xs=psi_N_Zeff - Zeff_spline%fs(:,1)=Zeff + IF(ALLOCATED(Zeff_) .AND. ALLOCATED(psi_N_Zeff_))THEN + CALL spline_alloc(Zeff_spline,SIZE(psi_N_Zeff_),1) + Zeff_spline%xs=psi_N_Zeff_ + Zeff_spline%fs(:,1)=Zeff_ CALL spline_fit(Zeff_spline,"extrap") ELSE CALL program_stop("Provide Zeff profile values") diff --git a/rdcon/rdcon_netcdf.f b/rdcon/rdcon_netcdf.f index 0534e7a10..0733d3f77 100644 --- a/rdcon/rdcon_netcdf.f +++ b/rdcon/rdcon_netcdf.f @@ -130,10 +130,10 @@ SUBROUTINE rdcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et) CALL check( nf90_put_att(ncid,nf90_global,'betat', betat)) CALL check( nf90_put_att(ncid,nf90_global,'betan', betan)) CALL check( nf90_put_att(ncid,nf90_global,'volume', volume)) - IF(ALLOCATED(Zeff) .AND. ALLOCATED(psi_N_Zeff))THEN - CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff)) + IF(ALLOCATED(Zeff_) .AND. ALLOCATED(psi_N_Zeff_))THEN + CALL check( nf90_put_att(ncid,nf90_global,'Zeff', Zeff_)) CALL check( nf90_put_att(ncid,nf90_global,'psi_N_Zeff', - $ psi_N_Zeff)) + $ psi_N_Zeff_)) ENDIF CALL check( nf90_put_att(ncid,nf90_global,'bwall', bwall)) CALL check( nf90_put_att(ncid,nf90_global,'bt0', bt0)) From 4584a7affc94d47940082467a20e50c9c25c8af2 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Fri, 3 Jul 2026 11:26:36 -0700 Subject: [PATCH 06/22] Making variable array read-in to treat each array separately. Adding explitcit indexing control to match Zeff spline assignment --- equil/local.f | 19 ++++++++----------- rdcon/dcon.f | 6 ++---- rdcon/dcon_mod.f | 2 +- rdcon/mercier.f | 2 +- 4 files changed, 12 insertions(+), 17 deletions(-) diff --git a/equil/local.f b/equil/local.f index fb7aba804..0cd9c436c 100644 --- a/equil/local.f +++ b/equil/local.f @@ -267,38 +267,35 @@ END FUNCTION floored_log c----------------------------------------------------------------------- c declarations. c----------------------------------------------------------------------- - SUBROUTINE read_var_len(x_in, y_in, x_out, y_out) + SUBROUTINE read_var_len(x_in, x_out) IMPLICIT NONE - REAL(r8), DIMENSION(:), INTENT(IN) :: x_in, y_in - REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: x_out, y_out + REAL(r8), DIMENSION(:), INTENT(IN) :: x_in ! Indexes from 1... + REAL(r8), DIMENSION(:), ALLOCATABLE, INTENT(OUT) :: x_out INTEGER :: n, i - REAL(r8) :: default_value_x,default_value_y + REAL(r8) :: default_value c----------------------------------------------------------------------- -c assume last values of x_in, y_in are sentinals +c assume last value of x_in is a sentinal c----------------------------------------------------------------------- - default_value_x = x_in(SIZE(x_in)) - default_value_y = y_in(SIZE(y_in)) + default_value = x_in(SIZE(x_in)) c----------------------------------------------------------------------- c count leading non-sentinel values in x_in. c----------------------------------------------------------------------- n = 0 DO i = 1, SIZE(x_in) - 1 - IF (x_in(i) /= default_value_x) THEN + IF (x_in(i) /= default_value) THEN n = n + 1 ELSE EXIT END IF END DO c----------------------------------------------------------------------- -c allocate output arrays and copy if values exist. +c allocate output array and copy if values exist. c----------------------------------------------------------------------- IF (n > 0) THEN ALLOCATE(x_out(n)) - ALLOCATE(y_out(n)) x_out(1:n) = x_in(1:n) - y_out(1:n) = y_in(1:n) END IF c----------------------------------------------------------------------- c terminate. diff --git a/rdcon/dcon.f b/rdcon/dcon.f index 8852772bb..ec13c5d87 100644 --- a/rdcon/dcon.f +++ b/rdcon/dcon.f @@ -262,10 +262,6 @@ SUBROUTINE dcon_run ENDIF CALL ascii_close(in_unit) c----------------------------------------------------------------------- -c Allocate variable-length Zeff input profile. -c----------------------------------------------------------------------- - CALL read_var_len(Zeff,psi_N_Zeff,Zeff_,psi_N_Zeff_) -c----------------------------------------------------------------------- c open output files, read, process, and diagnose equilibrium. c----------------------------------------------------------------------- CALL ascii_open(out_unit,"dcon.out","REPLACE") @@ -343,6 +339,8 @@ SUBROUTINE dcon_run c optionally compute modified Rutherford equation (MRE) terms. c----------------------------------------------------------------------- IF(MRE_flag)THEN + CALL read_var_len(Zeff,Zeff_) + CALL read_var_len(psi_N_Zeff,psi_N_Zeff_) CALL spline_alloc(mreterms,mpsi,30) mreterms%xs=sq%xs mreterms%fs=0 diff --git a/rdcon/dcon_mod.f b/rdcon/dcon_mod.f index 47604f187..74c161d1c 100644 --- a/rdcon/dcon_mod.f +++ b/rdcon/dcon_mod.f @@ -83,7 +83,7 @@ MODULE rdcon_mod $ plasma1=0.0,vacuum1=0.0,total1=0.0 REAL(r8), DIMENSION(:), ALLOCATABLE :: Zeff_, psi_N_Zeff_ - REAL(r8), DIMENSION(1000) :: Zeff=-1.0, psi_N_Zeff=-1.0 ! read-in + REAL(r8), DIMENSION(1:1000) :: Zeff=-1.0, psi_N_Zeff=-1.0 REAL(r8) :: dpsi_intvl=0.1,dpsi1_intvl=0.1 TYPE(spline_type) :: locstab, mreterms diff --git a/rdcon/mercier.f b/rdcon/mercier.f index cff17260b..2f2f7e4df 100644 --- a/rdcon/mercier.f +++ b/rdcon/mercier.f @@ -44,7 +44,7 @@ SUBROUTINE mercier_scan IF(MRE_flag)THEN CALL spline_alloc(ff,mtheta,22) IF(ALLOCATED(Zeff_) .AND. ALLOCATED(psi_N_Zeff_))THEN - CALL spline_alloc(Zeff_spline,SIZE(psi_N_Zeff_),1) + CALL spline_alloc(Zeff_spline,SIZE(psi_N_Zeff_)-1,1) Zeff_spline%xs=psi_N_Zeff_ Zeff_spline%fs(:,1)=Zeff_ CALL spline_fit(Zeff_spline,"extrap") From b3d713405e38a9791c4c02dffbf8a6cad6da164e Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Fri, 3 Jul 2026 11:26:57 -0700 Subject: [PATCH 07/22] Adding spline cleanup --- rdcon/dcon.f | 5 +++++ rdcon/mercier.f | 1 + 2 files changed, 6 insertions(+) diff --git a/rdcon/dcon.f b/rdcon/dcon.f index ec13c5d87..ee948fbca 100644 --- a/rdcon/dcon.f +++ b/rdcon/dcon.f @@ -77,6 +77,11 @@ SUBROUTINE dcon_dealloc(opt) DEALLOCATE(cellinfos%icell,cellinfos%iintvl, $ cellinfos%etypes,cellinfos%etypes_int,cellinfos%x1,cellinfos%x2) ENDIF + IF(MRE_flag)THEN + CALL spline_dealloc(mreterms) + DEALLOCATE(Zeff_) + DEALLOCATE(psi_N_Zeff_) + ENDIF c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- diff --git a/rdcon/mercier.f b/rdcon/mercier.f index 2f2f7e4df..b3e03ed23 100644 --- a/rdcon/mercier.f +++ b/rdcon/mercier.f @@ -288,6 +288,7 @@ SUBROUTINE mercier_scan 120 FORMAT(19(E30.15,1X)) ENDDO CALL spline_dealloc(ff) + IF(MRE_flag) CALL spline_dealloc(Zeff_spline) c----------------------------------------------------------------------- c terminate. c----------------------------------------------------------------------- From 572462737a2686cdb1e95778678d1b2a4a5b791e Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Fri, 3 Jul 2026 11:27:19 -0700 Subject: [PATCH 08/22] turning off MRE_flag, geom_flag defaults --- docs/examples/DIIID_resistive_example/rdcon.in | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/examples/DIIID_resistive_example/rdcon.in b/docs/examples/DIIID_resistive_example/rdcon.in index eef0ddbf1..36aebd085 100644 --- a/docs/examples/DIIID_resistive_example/rdcon.in +++ b/docs/examples/DIIID_resistive_example/rdcon.in @@ -75,8 +75,8 @@ bin_bal1=f ! Binary output for bal_flag poloidal functions out_bal2=f ! Ascii output for bal_flag functions bin_bal2=f ! Binary output for bal_flag functions - MRE_flag=t ! Netcdf output for various modified Rutherford equation terms (defined in mercier.f) - geom_flag=t ! Netcdf output for more flux-surface averaged terms (defined in mercier.f). To use, must set MRE_flag = t + MRE_flag=f ! Netcdf output for various modified Rutherford equation terms (defined in mercier.f) + geom_flag=f ! Netcdf output for more flux-surface averaged terms (defined in mercier.f). To use, must set MRE_flag = t / &UA_DIAGNOSE_LIST From 225d2d1fc1bf79ad913304f542e950a7a56ae50c Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 12:32:58 -0700 Subject: [PATCH 09/22] EQUIL - FEATURE - spline_fit_hermite: cubic spline fit that keeps caller-supplied knot derivatives Co-Authored-By: Claude Opus 5.5 --- equil/spline.f | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/equil/spline.f b/equil/spline.f index f7aca7896..4256cfc25 100644 --- a/equil/spline.f +++ b/equil/spline.f @@ -32,6 +32,7 @@ c 21. spline_thomas. c 22. spline_roots. c 23. spline_refine_root. +c 24. spline_fit_hermite. c----------------------------------------------------------------------- c subprogram 0. spline_type definition. c defines spline_type. @@ -1849,5 +1850,39 @@ SUBROUTINE spline_refine_root(spl, iqty, x, op_tol) c----------------------------------------------------------------------- RETURN END SUBROUTINE spline_refine_root +c----------------------------------------------------------------------- +c subprogram 24. spline_fit_hermite. +c fits cubic splines, keeping the caller's knot derivatives fs1 +c for the quantities flagged in keep (cubic Hermite interpolation). +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE spline_fit_hermite(spl,endmode,keep) + + TYPE(spline_type), INTENT(INOUT) :: spl + CHARACTER(*), INTENT(IN) :: endmode + LOGICAL, DIMENSION(spl%nqty), INTENT(IN) :: keep + + INTEGER :: iqty + REAL(r8), DIMENSION(0:spl%mx,spl%nqty) :: fs1 +c----------------------------------------------------------------------- +c fit all quantities, then restore the supplied derivatives. +c----------------------------------------------------------------------- + DO iqty=1,spl%nqty + IF(keep(iqty) .AND. ANY(spl%xpower(:,iqty) /= 0)) + $ CALL program_stop("spline_fit_hermite: "// + $ "supplied derivatives incompatible with xpower") + ENDDO + fs1=spl%fs1 + CALL spline_fit(spl,endmode) + DO iqty=1,spl%nqty + IF(keep(iqty))spl%fs1(:,iqty)=fs1(:,iqty) + ENDDO +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE spline_fit_hermite END MODULE spline_mod From a3e12fbd2b8edf180dc452e84bab5009ba0c5c93 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 12:38:39 -0700 Subject: [PATCH 10/22] EQUIL - FEATURE - profile_source: build f and p from the tabulated FF' and p' (efit, ldp_i) New equil.in variable profile_source: values - tabulated f, p; spline-fitted derivatives (default, unchanged) hermite - tabulated f, p with the tabulated derivatives as knot slopes integrate - f**2/2 and p integrated inward from the boundary values (as in Julia GPEC PR #506), tabulated derivatives as slopes read_eq_efit no longer discards FFPRIM and PPRIME, and fixes an out-of-bounds index when the g-file has no q profile. read_eq_ldp_i reads optional trailing FF' and p' records. Falls back to values with a warning when the derivatives are non-finite or disagree in sign. Co-Authored-By: Claude Opus 5.5 --- equil/equil.f | 3 +- equil/global.f | 2 + equil/read_eq.f | 141 +++++++++++++++++++++++++++++++++++++++++++++--- 3 files changed, 137 insertions(+), 9 deletions(-) diff --git a/equil/equil.f b/equil/equil.f index be7456d46..fbc12fe8f 100644 --- a/equil/equil.f +++ b/equil/equil.f @@ -39,7 +39,7 @@ SUBROUTINE equil_loadnamelists(op_psihigh, op_psilow) $ newq0,psihigh,psilow,input_only,jac_type,power_bp,power_r, $ power_b,jac_method,convert_type,power_flag,ns1, $ sp_pfac,sp_nx,sp_dx1,sp_dx2,use_galgrid,etol,enstep, - $ use_classic_splines + $ use_classic_splines,profile_source NAMELIST/equil_output/bin_2d,bin_eq_1d,bin_eq_2d,out_2d,out_eq_1d, $ out_eq_2d,bin_fl,out_fl,interp,gse_flag,dump_flag,verbose, $ out_ahg2msc @@ -150,6 +150,7 @@ SUBROUTINE equil_read(unit, op_psihigh, op_psilow) c----------------------------------------------------------------------- c read equilibrium data and diagnose. c----------------------------------------------------------------------- + sq_in_slopes=.FALSE. SELECT CASE(TRIM(eq_type)) CASE("fluxgrid") CALL read_eq_fluxgrid diff --git a/equil/global.f b/equil/global.f index b57137f02..53594f553 100644 --- a/equil/global.f +++ b/equil/global.f @@ -16,6 +16,8 @@ MODULE global_mod CHARACTER(16) :: eq_type="fluxgrid" CHARACTER(128) :: eq_filename="fluxgrid.dat" + CHARACTER(16) :: profile_source="values" + LOGICAL, DIMENSION(4) :: sq_in_slopes=.FALSE. LOGICAL :: direct_flag=.FALSE. LOGICAL :: out_eq_1d=.FALSE.,bin_eq_1d=.FALSE. diff --git a/equil/read_eq.f b/equil/read_eq.f index 3ad187791..ad14747bc 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -37,6 +37,8 @@ c 29. read_eq_chease4. c 30. read_eq_marklin_direct. c 31. read_eq_marklin_inverse. +c 32. read_eq_profiles_from_derivatives. +c 33. read_eq_profiles_warn. c----------------------------------------------------------------------- c subprogram 0. read_eq_mod. c module declarations. @@ -625,6 +627,7 @@ SUBROUTINE read_eq_efit REAL(r8) :: bcentr,cpasma,rgrid,rmaxis,rzero,ssibry1, $ ssibry2,ssimag1,ssimag2,xdim,xdum,zdim,zmaxis,zmid REAL(r8) :: tmp + REAL(r8), DIMENSION(:), ALLOCATABLE :: ffprim,pprime CHARACTER(LEN=7) :: shotstr, timestr c----------------------------------------------------------------------- c read equilibrium data. @@ -636,14 +639,15 @@ SUBROUTINE read_eq_efit READ(in_unit,'(5e16.9)')cpasma,ssimag2,xdum,rmaxis,xdum READ(in_unit,'(5e16.9)')zmaxis,xdum,ssibry2,xdum,xdum CALL spline_alloc(sq_in,nw-1,4) + ALLOCATE(ffprim(0:nw-1),pprime(0:nw-1)) READ(in_unit,'(5e16.9)')(sq_in%fs(i,1),i=0,nw-1) READ(in_unit,'(5e16.9)')(sq_in%fs(i,2),i=0,nw-1) - READ(in_unit,'(5e16.9)')(sq_in%fs(i,3),i=0,nw-1) - READ(in_unit,'(5e16.9)')(sq_in%fs(i,3),i=0,nw-1) + READ(in_unit,'(5e16.9)')(ffprim(i),i=0,nw-1) + READ(in_unit,'(5e16.9)')(pprime(i),i=0,nw-1) CALL bicube_alloc(psi_in,nw-1,nh-1,1) READ(in_unit,'(5e16.9)')((psi_in%fs(i,j,1),i=0,nw-1),j=0,nh-1) READ(in_unit,'(5e16.9)',iostat=ios)(sq_in%fs(i,3),i=0,nw-1) - IF(ios /= 0)sq_in%fs(i,3)=0 + IF(ios /= 0)sq_in%fs(:,3)=0 CALL ascii_close(in_unit) c----------------------------------------------------------------------- c translate to internal quantities. @@ -667,6 +671,9 @@ SUBROUTINE read_eq_efit sq_in%xs=(/(ia,ia=0,ma)/)/dfloat(ma) sq_in%fs(:,1)=ABS(sq_in%fs(:,1)) sq_in%fs(:,2)=MAX(sq_in%fs(:,2)*mu0,zero) + CALL read_eq_profiles_from_derivatives(ffprim*psio, + $ pprime*mu0*psio) + DEALLOCATE(ffprim,pprime) c----------------------------------------------------------------------- c copy and convert 2D quantities. c----------------------------------------------------------------------- @@ -808,8 +815,8 @@ END SUBROUTINE read_eq_ldp_d c----------------------------------------------------------------------- SUBROUTINE read_eq_ldp_i - INTEGER :: mx,my - REAL(r8), DIMENSION(:), POINTER :: psi,f,p,q + INTEGER :: mx,my,ios1,ios2 + REAL(r8), DIMENSION(:), POINTER :: psi,f,p,q,ffp,pp REAL(r8), DIMENSION(:,:), POINTER :: r,z c----------------------------------------------------------------------- c open data file, read scalar data. @@ -823,9 +830,9 @@ 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),ffp(0:mx),pp(0:mx)) c----------------------------------------------------------------------- -c read binary data. +c read binary data; f*df/dpsi and dp/dpsi records are optional. c----------------------------------------------------------------------- READ(in_unit)psi READ(in_unit)f @@ -833,6 +840,8 @@ SUBROUTINE read_eq_ldp_i READ(in_unit)q READ(in_unit)r READ(in_unit)z + READ(in_unit,IOSTAT=ios1)ffp + READ(in_unit,IOSTAT=ios2)pp CALL bin_close(in_unit) c----------------------------------------------------------------------- c copy and revise 1D arrays. @@ -843,6 +852,8 @@ SUBROUTINE read_eq_ldp_i sq_in%fs(:,1)=f sq_in%fs(:,2)=p*mu0 sq_in%fs(:,3)=q + IF(ios1 == 0 .AND. ios2 == 0) + $ CALL read_eq_profiles_from_derivatives(ffp*psio,pp*mu0*psio) if(psio.lt.0)psio=-psio c----------------------------------------------------------------------- c copy and revise 2D arrays. @@ -855,7 +866,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,ffp,pp) CALL inverse_run c----------------------------------------------------------------------- c terminate. @@ -2378,4 +2389,118 @@ SUBROUTINE read_eq_marklin_inverse c----------------------------------------------------------------------- RETURN END SUBROUTINE read_eq_marklin_inverse +c----------------------------------------------------------------------- +c subprogram 32. read_eq_profiles_from_derivatives. +c uses tabulated derivatives ffp = f*df/dpsifac and +c pp = d(mu0 p)/dpsifac for the f and mu0 p columns of sq_in, +c according to profile_source: +c "values": tabulated values, spline-fitted derivatives. +c "hermite": tabulated values and tabulated derivatives. +c "integrate": values integrated inward from the boundary +c values of f**2/2 and mu0 p, and tabulated derivatives. +c falls back to "values" if the derivatives are unusable. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) + + REAL(r8), DIMENSION(0:), INTENT(IN) :: ffp,pp + + INTEGER :: mx + REAL(r8) :: sgn_f,sgn_p,dfmax,dpmax + REAL(r8), DIMENSION(0:SIZE(ffp)-1) :: f,p,f2,fint,pint + TYPE(spline_type) :: spl +c----------------------------------------------------------------------- +c select method. +c----------------------------------------------------------------------- + SELECT CASE(TRIM(profile_source)) + CASE("values") + RETURN + CASE("hermite","integrate") + CASE DEFAULT + CALL program_stop("Cannot recognize profile_source = " + $ //TRIM(profile_source)) + END SELECT + mx=sq_in%mx + f=sq_in%fs(:,1) + p=sq_in%fs(:,2) +c----------------------------------------------------------------------- +c check derivatives are finite and agree in sign with the values. +c----------------------------------------------------------------------- + IF(.NOT.(ALL(ABS(ffp) <= HUGE(one)) .AND. + $ ALL(ABS(pp) <= HUGE(one)) .AND. ALL(f /= 0)))THEN + CALL read_eq_profiles_warn("non-finite derivatives or f = 0") + RETURN + ENDIF + sgn_f=SIGN(one,SUM((f(1:mx)-f(0:mx-1)) + $ *(ffp(1:mx)/f(1:mx)+ffp(0:mx-1)/f(0:mx-1)))) + sgn_p=SIGN(one,SUM((p(1:mx)-p(0:mx-1))*(pp(1:mx)+pp(0:mx-1)))) + IF(sgn_f /= sgn_p)THEN + CALL read_eq_profiles_warn("derivative signs disagree with" + $ //" the values for only one of f and p") + RETURN + ENDIF +c----------------------------------------------------------------------- +c integrate derivatives inward from the boundary. +c----------------------------------------------------------------------- + CALL spline_alloc(spl,mx,2) + spl%xs=sq_in%xs + spl%fs(:,1)=sgn_f*ffp + spl%fs(:,2)=sgn_f*pp + CALL spline_fit(spl,"extrap") + CALL spline_int(spl) + f2=f(mx)**2-2*(spl%fsi(mx,1)-spl%fsi(:,1)) + pint=p(mx)-(spl%fsi(mx,2)-spl%fsi(:,2)) + CALL spline_dealloc(spl) + IF(ANY(f2 <= 0))THEN + CALL read_eq_profiles_warn("integrated f**2 <= 0") + RETURN + ENDIF + fint=SIGN(SQRT(f2),f(mx)) +c----------------------------------------------------------------------- +c report consistency of tabulated values and derivatives. +c----------------------------------------------------------------------- + dfmax=MAXVAL(ABS(fint-f))/MAXVAL(ABS(f)) + dpmax=MAXVAL(ABS(pint-p))/MAX(MAXVAL(ABS(p)),TINY(one)) + IF(verbose)WRITE(*,'(1x,2a,1p,2(a,es9.2))')"profile_source = ", + $ TRIM(profile_source),", max |df|/|f| =",dfmax, + $ ", max |dp|/|p| =",dpmax + IF(dfmax > 1e-4 .OR. dpmax > 5e-2)WRITE(*,'(1x,a)') + $ "WARNING: tabulated values and derivatives of f, p disagree" +c----------------------------------------------------------------------- +c store values and derivatives. +c----------------------------------------------------------------------- + IF(profile_source == "integrate")THEN + sq_in%fs(:,1)=fint + sq_in%fs(:,2)=pint + ENDIF + sq_in%fs1(:,1)=sgn_f*ffp/sq_in%fs(:,1) + sq_in%fs1(:,2)=sgn_f*pp + sq_in_slopes(1:2)=.TRUE. +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE read_eq_profiles_from_derivatives +c----------------------------------------------------------------------- +c subprogram 33. read_eq_profiles_warn. +c warns that profile derivatives are not used. +c----------------------------------------------------------------------- +c----------------------------------------------------------------------- +c declarations. +c----------------------------------------------------------------------- + SUBROUTINE read_eq_profiles_warn(reason) + + CHARACTER(*), INTENT(IN) :: reason +c----------------------------------------------------------------------- +c write warning. +c----------------------------------------------------------------------- + WRITE(*,'(1x,3a)')"WARNING: profile_source = values used: ", + $ reason,"." +c----------------------------------------------------------------------- +c terminate. +c----------------------------------------------------------------------- + RETURN + END SUBROUTINE read_eq_profiles_warn END MODULE read_eq_mod From 49a09a9e04c859e82ddb7f6912ddda3e79317858 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 12:38:40 -0700 Subject: [PATCH 11/22] EQUIL - FEATURE - keep supplied f, p derivatives through direct and inverse sq construction sq_in and sq are fit with spline_fit_hermite so the profile slopes from profile_source are not re-derived by spline differentiation. newq0 rescaling updates the f slope consistently (f' -> f'/ffac). With profile_source = values, output is bit-for-bit unchanged. Co-Authored-By: Claude Opus 5.5 --- equil/direct.f | 9 ++++++--- equil/inverse.f | 10 +++++++--- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/equil/direct.f b/equil/direct.f index b58e3de29..2a0377a6c 100644 --- a/equil/direct.f +++ b/equil/direct.f @@ -74,7 +74,7 @@ SUBROUTINE direct_run sq_in%fs(:,4)=SQRT(sq_in%xs) sq_in%name=" sq " sq_in%title=(/"psifac"," f ","mu0 p "," q "," rho "/) - CALL spline_fit(sq_in,"extrap") + CALL spline_fit_hermite(sq_in,"extrap",sq_in_slopes) psi_in%xs=rmin+(/(ir,ir=0,psi_in%mx)/)*(rmax-rmin)/psi_in%mx psi_in%ys=zmin+(/(iz,iz=0,psi_in%my)/)*(zmax-zmin)/psi_in%my CALL bicube_fit(psi_in,"extrap","extrap") @@ -179,6 +179,8 @@ SUBROUTINE direct_run c----------------------------------------------------------------------- sq%fs(ipsi,1)=bf%f*twopi sq%fs(ipsi,2)=bf%p + sq%fs1(ipsi,1)=bf%f1*twopi + 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) @@ -191,7 +193,7 @@ SUBROUTINE direct_run c----------------------------------------------------------------------- c fit surface quantities to cubic splines. c----------------------------------------------------------------------- - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) sq%name=" sq " sq%title=(/" psi "," f "," p "," q "/) q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) @@ -205,10 +207,11 @@ SUBROUTINE direct_run DO ipsi=0,mpsi ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac + sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(1-sq%xs(mpsi)) c----------------------------------------------------------------------- diff --git a/equil/inverse.f b/equil/inverse.f index 922b027e6..06eceab6f 100644 --- a/equil/inverse.f +++ b/equil/inverse.f @@ -54,7 +54,7 @@ SUBROUTINE inverse_run sq_in%fs(:,4)=SQRT(sq_in%xs) sq_in%name=" sq " sq_in%title=(/"psifac"," f ","mu0 p "," q "," rho "/) - CALL spline_fit(sq_in,"extrap") + CALL spline_fit_hermite(sq_in,"extrap",sq_in_slopes) rz_in%xs=sq_in%xs rz_in%ys=(/(itheta,itheta=0,rz_in%my)/)/REAL(rz_in%my,r8) rz_in%name=" rz " @@ -230,12 +230,15 @@ SUBROUTINE inverse_run c----------------------------------------------------------------------- c compute surface quantities. c----------------------------------------------------------------------- + CALL spline_eval(sq_in,psifac,1) sq%fs(ipsi,1)=sq_in%f(1)*twopi sq%fs(ipsi,2)=sq_in%f(2) + sq%fs1(ipsi,1)=sq_in%f1(1)*twopi + sq%fs1(ipsi,2)=sq_in%f1(2) sq%fs(ipsi,3)=spl%fsi(mtheta,3)*twopi*pi sq%fs(ipsi,4)=spl%fsi(mtheta,4)*sq%fs(ipsi,1)/(2*twopi*psio) ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) c----------------------------------------------------------------------- c revise q profile. @@ -247,10 +250,11 @@ SUBROUTINE inverse_run DO ipsi=0,mpsi ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac + sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(one-sq%xs(mpsi)) c----------------------------------------------------------------------- From 5a8f0b9cda74b18fe6c7040e8cc155729d50a548 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 12:52:10 -0700 Subject: [PATCH 12/22] EQUIL - FIX - sq_out diagnostic keeps the supplied f, p derivatives Co-Authored-By: Claude Opus 5.5 --- equil/equil_out.f | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/equil/equil_out.f b/equil/equil_out.f index 1069179be..43c0e8f6c 100644 --- a/equil/equil_out.f +++ b/equil/equil_out.f @@ -97,7 +97,10 @@ SUBROUTINE equil_out_diagnose(sq_flag,unit) sq_out%fs(:,2)=sq%fs(:,2) sq_out%fs(:,3)=sq%fs(:,4) sq_out%fs(:,4)=SQRT(sq_out%xs) - CALL spline_fit(sq_out,"extrap") + sq_out%fs1(:,1)=sq%fs1(:,1)/twopi + sq_out%fs1(:,2)=sq%fs1(:,2) + CALL spline_fit_hermite(sq_out,"extrap", + $ (/sq_in_slopes(1:2),.FALSE.,.FALSE./)) sq_out%name=" sq " sq_out%title=(/"psifac"," f ","mu0 p "," q "," rho "/) IF(bin_eq_1d)CALL bin_open(bin_2d_unit,"sq_out.bin","UNKNOWN", From 9846ddb6acbdc65059dbeb9a57a5ac521d5c9d5c Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 12:52:10 -0700 Subject: [PATCH 13/22] EQUIL - FEATURE - profile_source defaults to integrate Against a TokaMaker equilibrium with known profiles, integrate gives F'' errors of 0.1-0.8% (0.1 < psi_n < 0.85) that fall with resolution, where values and hermite give 2-47% that grow with mpsi (single-precision f amplified by differentiation). STRIDE Delta' is correspondingly more consistent across g-file resolutions. Set profile_source="values" to recover the previous behaviour. Co-Authored-By: Claude Opus 5.5 --- equil/global.f | 2 +- input/equil.in | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/equil/global.f b/equil/global.f index 53594f553..45ee91050 100644 --- a/equil/global.f +++ b/equil/global.f @@ -16,7 +16,7 @@ MODULE global_mod CHARACTER(16) :: eq_type="fluxgrid" CHARACTER(128) :: eq_filename="fluxgrid.dat" - CHARACTER(16) :: profile_source="values" + CHARACTER(16) :: profile_source="integrate" LOGICAL, DIMENSION(4) :: sq_in_slopes=.FALSE. LOGICAL :: direct_flag=.FALSE. diff --git a/input/equil.in b/input/equil.in index 98cbb3a91..01685bf07 100644 --- a/input/equil.in +++ b/input/equil.in @@ -15,6 +15,7 @@ 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. etol = 1e-7 ! Tolerance for direct equilibrium reconstruction of flux surfaces/coordinates use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions + profile_source="integrate" ! How f and p are built when the file has FF' and p' (efit, ldp_i). integrate: integrate FF', p' inward from the boundary values and use them as spline derivatives (recommended). hermite: tabulated f, p with FF', p' as derivatives. values: tabulated f, p only (previous behaviour) input_only=f ! Generate information about the input and then quit with no further calculation use_galgrid=t ! Use the same grid as galerkin method when regrid=.True. From 414877b28a63967b25375e11cf32275866809e87 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 12:53:40 -0700 Subject: [PATCH 14/22] TESTS - NEW - spline unit tests: Hermite with supplied slopes, extrap end-slope order Co-Authored-By: Claude Opus 5.5 --- regression/spline_tests/makefile | 13 +++ .../spline_tests/test_spline_hermite.f90 | 86 +++++++++++++++++++ 2 files changed, 99 insertions(+) create mode 100644 regression/spline_tests/makefile create mode 100644 regression/spline_tests/test_spline_hermite.f90 diff --git a/regression/spline_tests/makefile b/regression/spline_tests/makefile new file mode 100644 index 000000000..2582f4ace --- /dev/null +++ b/regression/spline_tests/makefile @@ -0,0 +1,13 @@ +# Unit tests for equil/spline.f. Build GPEC first (install/make), then: make run +FC ?= gfortran +GPEC = ../.. +LIBS ?= -L$(LAPACKHOME)/lib -llapack -lblas + +test_spline_hermite: test_spline_hermite.f90 $(GPEC)/deps/lib/libequil.a + $(FC) $(FFLAGS) -I$(GPEC)/equil -o $@ $< $(GPEC)/deps/lib/libequil.a $(LIBS) + +run: test_spline_hermite + ./test_spline_hermite + +clean: + rm -f test_spline_hermite *.o *.mod diff --git a/regression/spline_tests/test_spline_hermite.f90 b/regression/spline_tests/test_spline_hermite.f90 new file mode 100644 index 000000000..3f2f23ffc --- /dev/null +++ b/regression/spline_tests/test_spline_hermite.f90 @@ -0,0 +1,86 @@ +!----------------------------------------------------------------------- +! Unit tests for spline_fit_hermite and the "extrap" end slopes. +! 1. cubic Hermite with exact knot slopes converges at 4th order. +! 2. "extrap" end slopes converge at 3rd order or better. +! 3. columns not kept are bit-for-bit those of spline_fit. +! Exits with status 1 on failure. +!----------------------------------------------------------------------- +PROGRAM test_spline_hermite + USE spline_mod + IMPLICIT NONE + + INTEGER, PARAMETER :: nres=4, mfine=2001 + INTEGER, DIMENSION(nres), PARAMETER :: mxs=(/16,32,64,128/) + INTEGER :: ires,ix,mx + REAL(r8), DIMENSION(nres) :: err_val,err_end + REAL(r8) :: x,order_val,order_end + REAL(r8), DIMENSION(:,:), ALLOCATABLE :: fs1_ref + TYPE(spline_type) :: spl + LOGICAL :: ok=.TRUE. + + DO ires=1,nres + mx=mxs(ires) + CALL spline_alloc(spl,mx,2) + spl%xs=(/(ix,ix=0,mx)/)/REAL(mx,r8) + spl%fs(:,1)=f(spl%xs) + spl%fs(:,2)=f(spl%xs) + spl%fs1(:,1)=df(spl%xs) + CALL spline_fit_hermite(spl,"extrap",(/.TRUE.,.FALSE./)) + err_val(ires)=0 + DO ix=0,mfine + x=ix/REAL(mfine,r8) + CALL spline_eval(spl,x,0) + err_val(ires)=MAX(err_val(ires),ABS(spl%f(1)-f1(x))) + ENDDO + err_end(ires)=MAX(ABS(spl%fs1(0,2)-df1(0._r8)),ABS(spl%fs1(mx,2)-df1(1._r8))) + ALLOCATE(fs1_ref(0:mx,2)) + fs1_ref=spl%fs1 + spl%fs(:,2)=f(spl%xs) + CALL spline_fit(spl,"extrap") + IF(ANY(spl%fs1(:,2) /= fs1_ref(:,2)))THEN + WRITE(*,*)"FAIL: column not kept differs from spline_fit, mx = ",mx + ok=.FALSE. + ENDIF + DEALLOCATE(fs1_ref) + CALL spline_dealloc(spl) + ENDDO + order_val=LOG(err_val(nres-1)/err_val(nres))/LOG(2._r8) + order_end=LOG(err_end(nres-1)/err_end(nres))/LOG(2._r8) + WRITE(*,'(a,4es10.2,a,f5.2)')" hermite max error:",err_val,", order",order_val + WRITE(*,'(a,4es10.2,a,f5.2)')" extrap end slope error:",err_end,", order",order_end + IF(order_val < 3.8_r8)THEN + WRITE(*,*)"FAIL: Hermite with exact slopes is not 4th order" + ok=.FALSE. + ENDIF + IF(order_end < 2.8_r8)THEN + WRITE(*,*)"FAIL: extrap end slopes are not 3rd order" + ok=.FALSE. + ENDIF + IF(.NOT.ok)STOP 1 + WRITE(*,*)"PASS" + +CONTAINS + + ELEMENTAL REAL(r8) FUNCTION f1(x) + REAL(r8), INTENT(IN) :: x + f1=SIN(3*x)+EXP(-x)*x**2 + END FUNCTION f1 + + ELEMENTAL REAL(r8) FUNCTION df1(x) + REAL(r8), INTENT(IN) :: x + df1=3*COS(3*x)+EXP(-x)*(2*x-x**2) + END FUNCTION df1 + + FUNCTION f(x) + REAL(r8), DIMENSION(:), INTENT(IN) :: x + REAL(r8), DIMENSION(SIZE(x)) :: f + f=f1(x) + END FUNCTION f + + FUNCTION df(x) + REAL(r8), DIMENSION(:), INTENT(IN) :: x + REAL(r8), DIMENSION(SIZE(x)) :: df + df=df1(x) + END FUNCTION df + +END PROGRAM test_spline_hermite From 822e8b598f658ac93b2e60d0e9bb89b4d54a4504 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 13:40:47 -0700 Subject: [PATCH 15/22] EQUIL - BUGFIX - Keep the 5-entry sq title in direct_run A second assignment of 4 titles reallocated sq%title, so the dump record was 6 bytes short and read_eq_dump failed (I/O past end of record). Co-Authored-By: Claude Opus 5.5 --- equil/direct.f | 1 - 1 file changed, 1 deletion(-) diff --git a/equil/direct.f b/equil/direct.f index 2a0377a6c..675b8f534 100644 --- a/equil/direct.f +++ b/equil/direct.f @@ -195,7 +195,6 @@ SUBROUTINE direct_run c----------------------------------------------------------------------- CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) sq%name=" sq " - sq%title=(/" psi "," f "," p "," q "/) q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) c----------------------------------------------------------------------- c revise q profile. From 8717d83240805cf4b6609680c8d359975744d6d6 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 13:40:48 -0700 Subject: [PATCH 16/22] EQUIL - FIX - eq_type dump keeps the f, p slopes and eqfun equil_out_dump appends the supplied-derivative flags with sq%fs1(:,1:2), then eqfun; read_eq_dump reads them when present (older dumps read as before) and no longer fails when the equilibrium is re-read (STRIDE psilim reform). DCON efit -> dump -> dcon reproduces f, mu0p, q, D_I, D_R exactly. Co-Authored-By: Claude Opus 5.5 --- equil/equil_out.f | 5 +++++ equil/read_eq.f | 28 +++++++++++++++++++++++++--- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/equil/equil_out.f b/equil/equil_out.f index 43c0e8f6c..18e6c9d22 100644 --- a/equil/equil_out.f +++ b/equil/equil_out.f @@ -848,6 +848,11 @@ SUBROUTINE equil_out_dump WRITE(dump_unit)rzphi%name WRITE(dump_unit)rzphi%periodic c----------------------------------------------------------------------- +c write supplied f, p derivatives and eqfun (read only if present). +c----------------------------------------------------------------------- + WRITE(dump_unit)sq_in_slopes(1:2),sq%fs1(:,1:2) + WRITE(dump_unit)eqfun%fs +c----------------------------------------------------------------------- c close dump file. c----------------------------------------------------------------------- CALL bin_close(dump_unit) diff --git a/equil/read_eq.f b/equil/read_eq.f index ad14747bc..3b2f1887a 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -1436,7 +1436,8 @@ END SUBROUTINE read_eq_rtaylor c----------------------------------------------------------------------- SUBROUTINE read_eq_dump - INTEGER :: nqty + INTEGER :: nqty,ios + LOGICAL, DIMENSION(2) :: slopes LOGICAL, PARAMETER :: diagnose=.TRUE. INTEGER :: ix,iy,iqty,iside,ipsi @@ -1468,14 +1469,15 @@ SUBROUTINE read_eq_dump c write 1D splines. c----------------------------------------------------------------------- READ(dump_unit)mpsi,nqty + IF(sq%allocated)CALL spline_dealloc(sq) CALL spline_alloc(sq,mpsi,nqty) READ(dump_unit)sq%xs,sq%fs, $ sq%xpower,sq%x0,sq%title,sq%name,sq%periodic - CALL spline_fit(sq,"extrap") c----------------------------------------------------------------------- c read 2D splines. c----------------------------------------------------------------------- READ(dump_unit)mpsi,mtheta,nqty + IF(rzphi%allocated)CALL bicube_dealloc(rzphi) CALL bicube_alloc(rzphi,mpsi,mtheta,nqty) READ(dump_unit)rzphi%xs READ(dump_unit)rzphi%ys @@ -1490,6 +1492,24 @@ SUBROUTINE read_eq_dump READ(dump_unit)rzphi%name READ(dump_unit)rzphi%periodic c----------------------------------------------------------------------- +c read supplied f, p derivatives and eqfun if present; fit sq. +c----------------------------------------------------------------------- + READ(dump_unit,IOSTAT=ios)slopes,sq%fs1(:,1:2) + IF(ios == 0)sq_in_slopes(1:2)=slopes + CALL spline_fit_hermite(sq,"extrap", + $ (/sq_in_slopes(1:2),.FALSE.,.FALSE./)) + IF(eqfun%allocated)CALL bicube_dealloc(eqfun) + CALL bicube_alloc(eqfun,mpsi,mtheta,3) + READ(dump_unit,IOSTAT=ios)eqfun%fs + IF(ios == 0)THEN + eqfun%xs=sq%xs + eqfun%ys=rzphi%ys + eqfun%title=(/" b0 "," "," " /) + CALL bicube_fit(eqfun,"extrap","periodic") + ELSE + CALL bicube_dealloc(eqfun) + ENDIF +c----------------------------------------------------------------------- c restore x powers and fit. c----------------------------------------------------------------------- ALLOCATE(xfac(0:rzphi%mx)) @@ -1516,10 +1536,12 @@ SUBROUTINE read_eq_dump DO ipsi=0,mpsi ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac + sq%fs1(ipsi,1)=sq%fs1(ipsi,1)/ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit(sq,"extrap") + CALL spline_fit_hermite(sq,"extrap", + $ (/sq_in_slopes(1:2),.FALSE.,.FALSE./)) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(one-sq%xs(mpsi)) c----------------------------------------------------------------------- From 9154d006a1352df300084d41d3d1532bf63a1038 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 14:24:11 -0700 Subject: [PATCH 17/22] EQUIL - CLEANUP - review fixes for profile_source and the dump slopes Derivative sign taken from f and p together, so a flat p (or f) no longer forces the values fallback; the dump writes eqfun only when it exists; shorter equil.in entry; unit test drops redundant wrappers. Co-Authored-By: Claude Opus 5.5 --- equil/equil_out.f | 2 +- equil/read_eq.f | 19 ++++++++++--------- input/equil.in | 2 +- regression/spline_tests/.gitignore | 1 + .../spline_tests/test_spline_hermite.f90 | 19 +++---------------- 5 files changed, 16 insertions(+), 27 deletions(-) create mode 100644 regression/spline_tests/.gitignore diff --git a/equil/equil_out.f b/equil/equil_out.f index 18e6c9d22..88d67fcf9 100644 --- a/equil/equil_out.f +++ b/equil/equil_out.f @@ -851,7 +851,7 @@ SUBROUTINE equil_out_dump c write supplied f, p derivatives and eqfun (read only if present). c----------------------------------------------------------------------- WRITE(dump_unit)sq_in_slopes(1:2),sq%fs1(:,1:2) - WRITE(dump_unit)eqfun%fs + IF(eqfun%allocated)WRITE(dump_unit)eqfun%fs c----------------------------------------------------------------------- c close dump file. c----------------------------------------------------------------------- diff --git a/equil/read_eq.f b/equil/read_eq.f index 3b2f1887a..8401b812e 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -2430,7 +2430,7 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) REAL(r8), DIMENSION(0:), INTENT(IN) :: ffp,pp INTEGER :: mx - REAL(r8) :: sgn_f,sgn_p,dfmax,dpmax + REAL(r8) :: sf,sp,sgn,dfmax,dpmax REAL(r8), DIMENSION(0:SIZE(ffp)-1) :: f,p,f2,fint,pint TYPE(spline_type) :: spl c----------------------------------------------------------------------- @@ -2455,10 +2455,11 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) CALL read_eq_profiles_warn("non-finite derivatives or f = 0") RETURN ENDIF - sgn_f=SIGN(one,SUM((f(1:mx)-f(0:mx-1)) - $ *(ffp(1:mx)/f(1:mx)+ffp(0:mx-1)/f(0:mx-1)))) - sgn_p=SIGN(one,SUM((p(1:mx)-p(0:mx-1))*(pp(1:mx)+pp(0:mx-1)))) - IF(sgn_f /= sgn_p)THEN + sf=SUM((f(1:mx)-f(0:mx-1)) + $ *(ffp(1:mx)/f(1:mx)+ffp(0:mx-1)/f(0:mx-1))) + sp=SUM((p(1:mx)-p(0:mx-1))*(pp(1:mx)+pp(0:mx-1))) + sgn=SIGN(one,sf+sp) + IF(sf*sp < 0)THEN CALL read_eq_profiles_warn("derivative signs disagree with" $ //" the values for only one of f and p") RETURN @@ -2468,8 +2469,8 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) c----------------------------------------------------------------------- CALL spline_alloc(spl,mx,2) spl%xs=sq_in%xs - spl%fs(:,1)=sgn_f*ffp - spl%fs(:,2)=sgn_f*pp + spl%fs(:,1)=sgn*ffp + spl%fs(:,2)=sgn*pp CALL spline_fit(spl,"extrap") CALL spline_int(spl) f2=f(mx)**2-2*(spl%fsi(mx,1)-spl%fsi(:,1)) @@ -2497,8 +2498,8 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) sq_in%fs(:,1)=fint sq_in%fs(:,2)=pint ENDIF - sq_in%fs1(:,1)=sgn_f*ffp/sq_in%fs(:,1) - sq_in%fs1(:,2)=sgn_f*pp + sq_in%fs1(:,1)=sgn*ffp/sq_in%fs(:,1) + sq_in%fs1(:,2)=sgn*pp sq_in_slopes(1:2)=.TRUE. c----------------------------------------------------------------------- c terminate. diff --git a/input/equil.in b/input/equil.in index 01685bf07..9ba3d176f 100644 --- a/input/equil.in +++ b/input/equil.in @@ -15,7 +15,7 @@ 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. etol = 1e-7 ! Tolerance for direct equilibrium reconstruction of flux surfaces/coordinates use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions - profile_source="integrate" ! How f and p are built when the file has FF' and p' (efit, ldp_i). integrate: integrate FF', p' inward from the boundary values and use them as spline derivatives (recommended). hermite: tabulated f, p with FF', p' as derivatives. values: tabulated f, p only (previous behaviour) + profile_source="integrate" ! f, p from the file's FF', p' (efit, ldp_i). integrate: integrated from the boundary, FF', p' as spline derivatives. hermite: file f, p with FF', p' as derivatives. values: file f, p only (previous behaviour) input_only=f ! Generate information about the input and then quit with no further calculation use_galgrid=t ! Use the same grid as galerkin method when regrid=.True. diff --git a/regression/spline_tests/.gitignore b/regression/spline_tests/.gitignore new file mode 100644 index 000000000..ca5a5db51 --- /dev/null +++ b/regression/spline_tests/.gitignore @@ -0,0 +1 @@ +test_spline_hermite diff --git a/regression/spline_tests/test_spline_hermite.f90 b/regression/spline_tests/test_spline_hermite.f90 index 3f2f23ffc..d13694e40 100644 --- a/regression/spline_tests/test_spline_hermite.f90 +++ b/regression/spline_tests/test_spline_hermite.f90 @@ -22,9 +22,9 @@ PROGRAM test_spline_hermite mx=mxs(ires) CALL spline_alloc(spl,mx,2) spl%xs=(/(ix,ix=0,mx)/)/REAL(mx,r8) - spl%fs(:,1)=f(spl%xs) - spl%fs(:,2)=f(spl%xs) - spl%fs1(:,1)=df(spl%xs) + spl%fs(:,1)=f1(spl%xs) + spl%fs(:,2)=f1(spl%xs) + spl%fs1(:,1)=df1(spl%xs) CALL spline_fit_hermite(spl,"extrap",(/.TRUE.,.FALSE./)) err_val(ires)=0 DO ix=0,mfine @@ -35,7 +35,6 @@ PROGRAM test_spline_hermite err_end(ires)=MAX(ABS(spl%fs1(0,2)-df1(0._r8)),ABS(spl%fs1(mx,2)-df1(1._r8))) ALLOCATE(fs1_ref(0:mx,2)) fs1_ref=spl%fs1 - spl%fs(:,2)=f(spl%xs) CALL spline_fit(spl,"extrap") IF(ANY(spl%fs1(:,2) /= fs1_ref(:,2)))THEN WRITE(*,*)"FAIL: column not kept differs from spline_fit, mx = ",mx @@ -71,16 +70,4 @@ ELEMENTAL REAL(r8) FUNCTION df1(x) df1=3*COS(3*x)+EXP(-x)*(2*x-x**2) END FUNCTION df1 - FUNCTION f(x) - REAL(r8), DIMENSION(:), INTENT(IN) :: x - REAL(r8), DIMENSION(SIZE(x)) :: f - f=f1(x) - END FUNCTION f - - FUNCTION df(x) - REAL(r8), DIMENSION(:), INTENT(IN) :: x - REAL(r8), DIMENSION(SIZE(x)) :: df - df=df1(x) - END FUNCTION df - END PROGRAM test_spline_hermite From ef882527aee107bede190896486cc15fe0fc4131 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 6 Oct 2026 15:00:46 -0700 Subject: [PATCH 18/22] EQUIL - CLEANUP - drop profile_source = hermite Tabulated values paired with tabulated derivatives put the values' rounding noise into f'' (the worst option in the comparison). spline_fit_hermite stays, used by integrate to keep the derivatives. Co-Authored-By: Claude Opus 5.5 --- equil/read_eq.f | 26 ++++++++++++-------------- input/equil.in | 2 +- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/equil/read_eq.f b/equil/read_eq.f index 8401b812e..acdf5fff6 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -2413,14 +2413,14 @@ SUBROUTINE read_eq_marklin_inverse END SUBROUTINE read_eq_marklin_inverse c----------------------------------------------------------------------- c subprogram 32. read_eq_profiles_from_derivatives. -c uses tabulated derivatives ffp = f*df/dpsifac and -c pp = d(mu0 p)/dpsifac for the f and mu0 p columns of sq_in, -c according to profile_source: -c "values": tabulated values, spline-fitted derivatives. -c "hermite": tabulated values and tabulated derivatives. -c "integrate": values integrated inward from the boundary -c values of f**2/2 and mu0 p, and tabulated derivatives. -c falls back to "values" if the derivatives are unusable. +c profile_source = "integrate": builds the f and mu0 p columns of +c sq_in from ffp = f*df/dpsifac and pp = d(mu0 p)/dpsifac. Values +c are integrated inward from the boundary values of f**2/2 and +c mu0 p; knot derivatives ffp/f and pp are kept by +c spline_fit_hermite. Tabulated values are not paired with the +c tabulated derivatives: their rounding noise would set f''. +c "values": tabulated values only. Falls back to "values" if the +c derivatives are unusable. c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. @@ -2439,7 +2439,7 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) SELECT CASE(TRIM(profile_source)) CASE("values") RETURN - CASE("hermite","integrate") + CASE("integrate") CASE DEFAULT CALL program_stop("Cannot recognize profile_source = " $ //TRIM(profile_source)) @@ -2494,11 +2494,9 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) c----------------------------------------------------------------------- c store values and derivatives. c----------------------------------------------------------------------- - IF(profile_source == "integrate")THEN - sq_in%fs(:,1)=fint - sq_in%fs(:,2)=pint - ENDIF - sq_in%fs1(:,1)=sgn*ffp/sq_in%fs(:,1) + sq_in%fs(:,1)=fint + sq_in%fs(:,2)=pint + sq_in%fs1(:,1)=sgn*ffp/fint sq_in%fs1(:,2)=sgn*pp sq_in_slopes(1:2)=.TRUE. c----------------------------------------------------------------------- diff --git a/input/equil.in b/input/equil.in index 9ba3d176f..0a80b7d4f 100644 --- a/input/equil.in +++ b/input/equil.in @@ -15,7 +15,7 @@ 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. etol = 1e-7 ! Tolerance for direct equilibrium reconstruction of flux surfaces/coordinates use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions - profile_source="integrate" ! f, p from the file's FF', p' (efit, ldp_i). integrate: integrated from the boundary, FF', p' as spline derivatives. hermite: file f, p with FF', p' as derivatives. values: file f, p only (previous behaviour) + profile_source="integrate" ! f, p from the file's FF', p' (efit, ldp_i). integrate: integrated from the boundary, FF', p' as spline derivatives. values: file f, p only (previous behaviour) input_only=f ! Generate information about the input and then quit with no further calculation use_galgrid=t ! Use the same grid as galerkin method when regrid=.True. From a47fc0205ac53f7b71fded73eec634764808d4de Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Mon, 12 Oct 2026 15:40:26 -0700 Subject: [PATCH 19/22] EQUIL - FEATURE - profile_source defaults to values; integrate is opt-in Review (d-burg): with double-precision g-files the tabulated f, p are clean, and integrate can move outer surfaces because the quadrature assumes a profile shape between knots near pedestal kinks. Keep the previous behaviour (values) as the default and offer integrate for single-precision g-files, where differentiating rounded f, p sets f''. With the default, the DIIID_ideal (DCON) and DIIID_resistive (RDCON) examples are bit-identical to the base branch. Co-Authored-By: Claude Opus 5.5 --- equil/global.f | 2 +- equil/read_eq.f | 18 ++++++++++-------- input/equil.in | 2 +- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/equil/global.f b/equil/global.f index 45ee91050..53594f553 100644 --- a/equil/global.f +++ b/equil/global.f @@ -16,7 +16,7 @@ MODULE global_mod CHARACTER(16) :: eq_type="fluxgrid" CHARACTER(128) :: eq_filename="fluxgrid.dat" - CHARACTER(16) :: profile_source="integrate" + CHARACTER(16) :: profile_source="values" LOGICAL, DIMENSION(4) :: sq_in_slopes=.FALSE. LOGICAL :: direct_flag=.FALSE. diff --git a/equil/read_eq.f b/equil/read_eq.f index acdf5fff6..fbea12931 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -2413,14 +2413,16 @@ SUBROUTINE read_eq_marklin_inverse END SUBROUTINE read_eq_marklin_inverse c----------------------------------------------------------------------- c subprogram 32. read_eq_profiles_from_derivatives. -c profile_source = "integrate": builds the f and mu0 p columns of -c sq_in from ffp = f*df/dpsifac and pp = d(mu0 p)/dpsifac. Values -c are integrated inward from the boundary values of f**2/2 and -c mu0 p; knot derivatives ffp/f and pp are kept by -c spline_fit_hermite. Tabulated values are not paired with the -c tabulated derivatives: their rounding noise would set f''. -c "values": tabulated values only. Falls back to "values" if the -c derivatives are unusable. +c profile_source = "values" (default): tabulated values only. +c "integrate" (opt-in, for single-precision g-files): builds the f +c and mu0 p columns of sq_in from ffp = f*df/dpsifac and +c pp = d(mu0 p)/dpsifac. Values are integrated inward from the +c boundary values of f**2/2 and mu0 p; knot derivatives ffp/f and +c pp are kept by spline_fit_hermite. Tabulated values are not +c paired with the tabulated derivatives: their rounding noise +c would set f''. The quadrature between knots can move outer +c surfaces near sharp profile features. Falls back to "values" if +c the derivatives are unusable. c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. diff --git a/input/equil.in b/input/equil.in index 0a80b7d4f..761ec70a5 100644 --- a/input/equil.in +++ b/input/equil.in @@ -15,7 +15,7 @@ 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. etol = 1e-7 ! Tolerance for direct equilibrium reconstruction of flux surfaces/coordinates use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions - profile_source="integrate" ! f, p from the file's FF', p' (efit, ldp_i). integrate: integrated from the boundary, FF', p' as spline derivatives. values: file f, p only (previous behaviour) + profile_source="values" ! How f, p are built for files with FF', p' (efit, ldp_i). values: tabulated f, p (default). integrate: integrated from the boundary with FF', p' as spline derivatives; opt-in for single-precision g-files, can move outer surfaces input_only=f ! Generate information about the input and then quit with no further calculation use_galgrid=t ! Use the same grid as galerkin method when regrid=.True. From 78977ab94141f9a1dfa26fd9cdcb4383e9aa95c7 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Mon, 12 Oct 2026 20:21:59 -0700 Subject: [PATCH 20/22] EQUIL - CLEANUP - one sq_in evaluation per surface, uniform sq_in_slopes, ldp_i profile_source warning - inverse_run evaluates sq_in once per surface with mode 1 instead of a third evaluation at the same psifac. - Every spline_fit_hermite call on sq, sq_in and sq_out passes sq_in_slopes; entries 3:4 are always false. - ldp_i warns when profile_source is not "values" and the file has no ffp, pp records, as the other fallbacks do. Default (values) DCON DIIID_ideal, RDCON DIIID_resistive and ldp_i outputs are bit-identical to a4b144d5, as are sq_out.bin, the dump and dump re-reads (newq0 = 0 and 1.1), and ldp_i with integrate. Co-Authored-By: Claude Opus 5.5 --- equil/equil_out.f | 3 +-- equil/inverse.f | 3 +-- equil/read_eq.f | 13 +++++++------ 3 files changed, 9 insertions(+), 10 deletions(-) diff --git a/equil/equil_out.f b/equil/equil_out.f index 88d67fcf9..3428b754f 100644 --- a/equil/equil_out.f +++ b/equil/equil_out.f @@ -99,8 +99,7 @@ SUBROUTINE equil_out_diagnose(sq_flag,unit) sq_out%fs(:,4)=SQRT(sq_out%xs) sq_out%fs1(:,1)=sq%fs1(:,1)/twopi sq_out%fs1(:,2)=sq%fs1(:,2) - CALL spline_fit_hermite(sq_out,"extrap", - $ (/sq_in_slopes(1:2),.FALSE.,.FALSE./)) + CALL spline_fit_hermite(sq_out,"extrap",sq_in_slopes) sq_out%name=" sq " sq_out%title=(/"psifac"," f ","mu0 p "," q "," rho "/) IF(bin_eq_1d)CALL bin_open(bin_2d_unit,"sq_out.bin","UNKNOWN", diff --git a/equil/inverse.f b/equil/inverse.f index 06eceab6f..ac19224bb 100644 --- a/equil/inverse.f +++ b/equil/inverse.f @@ -180,7 +180,7 @@ SUBROUTINE inverse_run CALL spline_alloc(spl,mtheta,5) DO ipsi=0,mpsi psifac=rzphi%xs(ipsi) - CALL spline_eval(sq_in,psifac,0) + CALL spline_eval(sq_in,psifac,1) spl%xs=rzphi%ys c----------------------------------------------------------------------- c interpolate to new surface. @@ -230,7 +230,6 @@ SUBROUTINE inverse_run c----------------------------------------------------------------------- c compute surface quantities. c----------------------------------------------------------------------- - CALL spline_eval(sq_in,psifac,1) sq%fs(ipsi,1)=sq_in%f(1)*twopi sq%fs(ipsi,2)=sq_in%f(2) sq%fs1(ipsi,1)=sq_in%f1(1)*twopi diff --git a/equil/read_eq.f b/equil/read_eq.f index fbea12931..0c13bfdba 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -852,8 +852,11 @@ SUBROUTINE read_eq_ldp_i sq_in%fs(:,1)=f sq_in%fs(:,2)=p*mu0 sq_in%fs(:,3)=q - IF(ios1 == 0 .AND. ios2 == 0) - $ CALL read_eq_profiles_from_derivatives(ffp*psio,pp*mu0*psio) + IF(ios1 == 0 .AND. ios2 == 0)THEN + CALL read_eq_profiles_from_derivatives(ffp*psio,pp*mu0*psio) + ELSEIF(TRIM(profile_source) /= "values")THEN + CALL read_eq_profiles_warn("no ffp, pp records in the file") + ENDIF if(psio.lt.0)psio=-psio c----------------------------------------------------------------------- c copy and revise 2D arrays. @@ -1496,8 +1499,7 @@ SUBROUTINE read_eq_dump c----------------------------------------------------------------------- READ(dump_unit,IOSTAT=ios)slopes,sq%fs1(:,1:2) IF(ios == 0)sq_in_slopes(1:2)=slopes - CALL spline_fit_hermite(sq,"extrap", - $ (/sq_in_slopes(1:2),.FALSE.,.FALSE./)) + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) IF(eqfun%allocated)CALL bicube_dealloc(eqfun) CALL bicube_alloc(eqfun,mpsi,mtheta,3) READ(dump_unit,IOSTAT=ios)eqfun%fs @@ -1540,8 +1542,7 @@ SUBROUTINE read_eq_dump sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac ENDDO - CALL spline_fit_hermite(sq,"extrap", - $ (/sq_in_slopes(1:2),.FALSE.,.FALSE./)) + CALL spline_fit_hermite(sq,"extrap",sq_in_slopes) ENDIF qa=sq%fs(mpsi,4)+sq%fs1(mpsi,4)*(one-sq%xs(mpsi)) c----------------------------------------------------------------------- From 0960ae5723e298e298948c7dcfd1316a42d96870 Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Tue, 13 Oct 2026 01:37:36 -0700 Subject: [PATCH 21/22] EQUIL - MINOR - Unit test checks the kept slopes exactly; shorter profile_source comments Co-Authored-By: Claude Opus 5.5 --- equil/equil_out.f | 2 +- equil/read_eq.f | 19 ++--- input/equil.in | 2 +- .../spline_tests/test_spline_hermite.f90 | 75 +++++-------------- 4 files changed, 30 insertions(+), 68 deletions(-) diff --git a/equil/equil_out.f b/equil/equil_out.f index 3428b754f..bb4fb90d8 100644 --- a/equil/equil_out.f +++ b/equil/equil_out.f @@ -847,7 +847,7 @@ SUBROUTINE equil_out_dump WRITE(dump_unit)rzphi%name WRITE(dump_unit)rzphi%periodic c----------------------------------------------------------------------- -c write supplied f, p derivatives and eqfun (read only if present). +c write f, p knot derivatives and eqfun; optional for read_eq_dump. c----------------------------------------------------------------------- WRITE(dump_unit)sq_in_slopes(1:2),sq%fs1(:,1:2) IF(eqfun%allocated)WRITE(dump_unit)eqfun%fs diff --git a/equil/read_eq.f b/equil/read_eq.f index 0c13bfdba..943973687 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -2414,16 +2414,12 @@ SUBROUTINE read_eq_marklin_inverse END SUBROUTINE read_eq_marklin_inverse c----------------------------------------------------------------------- c subprogram 32. read_eq_profiles_from_derivatives. -c profile_source = "values" (default): tabulated values only. -c "integrate" (opt-in, for single-precision g-files): builds the f -c and mu0 p columns of sq_in from ffp = f*df/dpsifac and -c pp = d(mu0 p)/dpsifac. Values are integrated inward from the -c boundary values of f**2/2 and mu0 p; knot derivatives ffp/f and -c pp are kept by spline_fit_hermite. Tabulated values are not -c paired with the tabulated derivatives: their rounding noise -c would set f''. The quadrature between knots can move outer -c surfaces near sharp profile features. Falls back to "values" if -c the derivatives are unusable. +c profile_source = "integrate": replaces the f and mu0 p columns +c of sq_in by integrals of ffp = f*df/dpsifac and +c pp = d(mu0 p)/dpsifac inward from the boundary values, and sets +c ffp/f and pp as knot derivatives for spline_fit_hermite. Leaves +c sq_in unchanged for "values", or, with a warning, if the +c derivatives are unusable. c----------------------------------------------------------------------- c----------------------------------------------------------------------- c declarations. @@ -2451,7 +2447,8 @@ SUBROUTINE read_eq_profiles_from_derivatives(ffp,pp) f=sq_in%fs(:,1) p=sq_in%fs(:,2) c----------------------------------------------------------------------- -c check derivatives are finite and agree in sign with the values. +c check derivatives are finite and that f and p agree on their sign +c relative to the values (that common sign is applied below). c----------------------------------------------------------------------- IF(.NOT.(ALL(ABS(ffp) <= HUGE(one)) .AND. $ ALL(ABS(pp) <= HUGE(one)) .AND. ALL(f /= 0)))THEN diff --git a/input/equil.in b/input/equil.in index 761ec70a5..d922f0ed6 100644 --- a/input/equil.in +++ b/input/equil.in @@ -15,7 +15,7 @@ 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. etol = 1e-7 ! Tolerance for direct equilibrium reconstruction of flux surfaces/coordinates use_classic_splines = f ! Use a classical cubic spline instead of tri-diagonal solution for splines with extrapolation boundary conditions - profile_source="values" ! How f, p are built for files with FF', p' (efit, ldp_i). values: tabulated f, p (default). integrate: integrated from the boundary with FF', p' as spline derivatives; opt-in for single-precision g-files, can move outer surfaces + profile_source="values" ! efit, and ldp_i files with FF', p' records. "values" (default): spline the tabulated f, p. "integrate": integrate FF', p' inward from the boundary and keep them as knot derivatives; for g-files with single-precision values; can shift pedestal surfaces input_only=f ! Generate information about the input and then quit with no further calculation use_galgrid=t ! Use the same grid as galerkin method when regrid=.True. diff --git a/regression/spline_tests/test_spline_hermite.f90 b/regression/spline_tests/test_spline_hermite.f90 index d13694e40..d57f1d415 100644 --- a/regression/spline_tests/test_spline_hermite.f90 +++ b/regression/spline_tests/test_spline_hermite.f90 @@ -1,73 +1,38 @@ !----------------------------------------------------------------------- -! Unit tests for spline_fit_hermite and the "extrap" end slopes. -! 1. cubic Hermite with exact knot slopes converges at 4th order. -! 2. "extrap" end slopes converge at 3rd order or better. -! 3. columns not kept are bit-for-bit those of spline_fit. +! Unit test for spline_fit_hermite: flagged columns keep the supplied +! knot derivatives; other columns are those of spline_fit. ! Exits with status 1 on failure. !----------------------------------------------------------------------- PROGRAM test_spline_hermite USE spline_mod IMPLICIT NONE - INTEGER, PARAMETER :: nres=4, mfine=2001 - INTEGER, DIMENSION(nres), PARAMETER :: mxs=(/16,32,64,128/) - INTEGER :: ires,ix,mx - REAL(r8), DIMENSION(nres) :: err_val,err_end - REAL(r8) :: x,order_val,order_end - REAL(r8), DIMENSION(:,:), ALLOCATABLE :: fs1_ref + INTEGER, PARAMETER :: mx=32 + INTEGER :: ix + REAL(r8), DIMENSION(0:mx) :: slopes + REAL(r8), DIMENSION(0:mx,2) :: fs1_ref TYPE(spline_type) :: spl LOGICAL :: ok=.TRUE. - DO ires=1,nres - mx=mxs(ires) - CALL spline_alloc(spl,mx,2) - spl%xs=(/(ix,ix=0,mx)/)/REAL(mx,r8) - spl%fs(:,1)=f1(spl%xs) - spl%fs(:,2)=f1(spl%xs) - spl%fs1(:,1)=df1(spl%xs) - CALL spline_fit_hermite(spl,"extrap",(/.TRUE.,.FALSE./)) - err_val(ires)=0 - DO ix=0,mfine - x=ix/REAL(mfine,r8) - CALL spline_eval(spl,x,0) - err_val(ires)=MAX(err_val(ires),ABS(spl%f(1)-f1(x))) - ENDDO - err_end(ires)=MAX(ABS(spl%fs1(0,2)-df1(0._r8)),ABS(spl%fs1(mx,2)-df1(1._r8))) - ALLOCATE(fs1_ref(0:mx,2)) - fs1_ref=spl%fs1 - CALL spline_fit(spl,"extrap") - IF(ANY(spl%fs1(:,2) /= fs1_ref(:,2)))THEN - WRITE(*,*)"FAIL: column not kept differs from spline_fit, mx = ",mx - ok=.FALSE. - ENDIF - DEALLOCATE(fs1_ref) - CALL spline_dealloc(spl) - ENDDO - order_val=LOG(err_val(nres-1)/err_val(nres))/LOG(2._r8) - order_end=LOG(err_end(nres-1)/err_end(nres))/LOG(2._r8) - WRITE(*,'(a,4es10.2,a,f5.2)')" hermite max error:",err_val,", order",order_val - WRITE(*,'(a,4es10.2,a,f5.2)')" extrap end slope error:",err_end,", order",order_end - IF(order_val < 3.8_r8)THEN - WRITE(*,*)"FAIL: Hermite with exact slopes is not 4th order" + CALL spline_alloc(spl,mx,2) + spl%xs=(/(ix,ix=0,mx)/)/REAL(mx,r8) + spl%fs(:,1)=SIN(3*spl%xs) + spl%fs(:,2)=spl%fs(:,1) + CALL spline_fit(spl,"extrap") + fs1_ref=spl%fs1 + slopes=3*COS(3*spl%xs) + spl%fs1(:,1)=slopes + CALL spline_fit_hermite(spl,"extrap",(/.TRUE.,.FALSE./)) + IF(ANY(spl%fs1(:,1) /= slopes))THEN + WRITE(*,*)"FAIL: flagged column lost the supplied derivatives" ok=.FALSE. ENDIF - IF(order_end < 2.8_r8)THEN - WRITE(*,*)"FAIL: extrap end slopes are not 3rd order" + IF(ANY(spl%fs1(:,2) /= fs1_ref(:,2)))THEN + WRITE(*,*)"FAIL: column not flagged differs from spline_fit" ok=.FALSE. ENDIF + CALL spline_dealloc(spl) IF(.NOT.ok)STOP 1 WRITE(*,*)"PASS" -CONTAINS - - ELEMENTAL REAL(r8) FUNCTION f1(x) - REAL(r8), INTENT(IN) :: x - f1=SIN(3*x)+EXP(-x)*x**2 - END FUNCTION f1 - - ELEMENTAL REAL(r8) FUNCTION df1(x) - REAL(r8), INTENT(IN) :: x - df1=3*COS(3*x)+EXP(-x)*(2*x-x**2) - END FUNCTION df1 - END PROGRAM test_spline_hermite From 823c713f37fc6e8bf1c5e564ea35322d0936c3ec Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 15 Oct 2026 12:24:27 +0200 Subject: [PATCH 22/22] EQUIL - MINOR - get rid of unused regression folder additions --- regression/spline_tests/.gitignore | 1 - regression/spline_tests/makefile | 13 ------- .../spline_tests/test_spline_hermite.f90 | 38 ------------------- 3 files changed, 52 deletions(-) delete mode 100644 regression/spline_tests/.gitignore delete mode 100644 regression/spline_tests/makefile delete mode 100644 regression/spline_tests/test_spline_hermite.f90 diff --git a/regression/spline_tests/.gitignore b/regression/spline_tests/.gitignore deleted file mode 100644 index ca5a5db51..000000000 --- a/regression/spline_tests/.gitignore +++ /dev/null @@ -1 +0,0 @@ -test_spline_hermite diff --git a/regression/spline_tests/makefile b/regression/spline_tests/makefile deleted file mode 100644 index 2582f4ace..000000000 --- a/regression/spline_tests/makefile +++ /dev/null @@ -1,13 +0,0 @@ -# Unit tests for equil/spline.f. Build GPEC first (install/make), then: make run -FC ?= gfortran -GPEC = ../.. -LIBS ?= -L$(LAPACKHOME)/lib -llapack -lblas - -test_spline_hermite: test_spline_hermite.f90 $(GPEC)/deps/lib/libequil.a - $(FC) $(FFLAGS) -I$(GPEC)/equil -o $@ $< $(GPEC)/deps/lib/libequil.a $(LIBS) - -run: test_spline_hermite - ./test_spline_hermite - -clean: - rm -f test_spline_hermite *.o *.mod diff --git a/regression/spline_tests/test_spline_hermite.f90 b/regression/spline_tests/test_spline_hermite.f90 deleted file mode 100644 index d57f1d415..000000000 --- a/regression/spline_tests/test_spline_hermite.f90 +++ /dev/null @@ -1,38 +0,0 @@ -!----------------------------------------------------------------------- -! Unit test for spline_fit_hermite: flagged columns keep the supplied -! knot derivatives; other columns are those of spline_fit. -! Exits with status 1 on failure. -!----------------------------------------------------------------------- -PROGRAM test_spline_hermite - USE spline_mod - IMPLICIT NONE - - INTEGER, PARAMETER :: mx=32 - INTEGER :: ix - REAL(r8), DIMENSION(0:mx) :: slopes - REAL(r8), DIMENSION(0:mx,2) :: fs1_ref - TYPE(spline_type) :: spl - LOGICAL :: ok=.TRUE. - - CALL spline_alloc(spl,mx,2) - spl%xs=(/(ix,ix=0,mx)/)/REAL(mx,r8) - spl%fs(:,1)=SIN(3*spl%xs) - spl%fs(:,2)=spl%fs(:,1) - CALL spline_fit(spl,"extrap") - fs1_ref=spl%fs1 - slopes=3*COS(3*spl%xs) - spl%fs1(:,1)=slopes - CALL spline_fit_hermite(spl,"extrap",(/.TRUE.,.FALSE./)) - IF(ANY(spl%fs1(:,1) /= slopes))THEN - WRITE(*,*)"FAIL: flagged column lost the supplied derivatives" - ok=.FALSE. - ENDIF - IF(ANY(spl%fs1(:,2) /= fs1_ref(:,2)))THEN - WRITE(*,*)"FAIL: column not flagged differs from spline_fit" - ok=.FALSE. - ENDIF - CALL spline_dealloc(spl) - IF(.NOT.ok)STOP 1 - WRITE(*,*)"PASS" - -END PROGRAM test_spline_hermite