From bab3554771f287ffb9576985f5a6801c451ee1e5 Mon Sep 17 00:00:00 2001 From: Chris Hansen Date: Wed, 15 Apr 2026 11:57:35 -0400 Subject: [PATCH 01/13] TokaMaker: Improve match to Ip target when using Jphi and P profiles --- src/physics/grad_shaf_prof_phys.F90 | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/src/physics/grad_shaf_prof_phys.F90 b/src/physics/grad_shaf_prof_phys.F90 index 1e05924bd..4194a63ad 100644 --- a/src/physics/grad_shaf_prof_phys.F90 +++ b/src/physics/grad_shaf_prof_phys.F90 @@ -54,6 +54,7 @@ module grad_shaf_prof_phys type, extends(linterp_flux_func) :: jphi_flux_func integer(4) :: ngeom = 50 !< Number of points in psi for , <1/R> evaluation real(8) :: j0 = 0.d0 !< LCFS Jphi value + real(8) :: norm_last = 1.d0 !< Last Jphi normalization factor (for Ip target) real(8), pointer, dimension(:) :: jphi => NULL() !< Jphi(psi) profile values contains !> Needs docs @@ -438,14 +439,25 @@ subroutine jphi_update(self,gseq) CALL spline_fit(R_spline,"extrap") DEALLOCATE(ravgs,psi_q,qtmp) !---Update jphi normalization to match Ip target -ALLOCATE(qtmp(self%npsi)) -DO i=1,self%npsi - CALL spline_eval(R_spline,self%x(i),0) - qtmp(i) = R_spline%f(1)*R_spline%f(2) -END DO -CALL gs_flux_int(gseq,self%x,self%jphi/qtmp,self%npsi,jphi_norm) -jphi_norm=ABS(gseq%Itor_target)/jphi_norm -DEALLOCATE(qtmp) +IF(gseq%Itor_target>0.d0)THEN + ALLOCATE(qtmp(self%npsi)) + DO i=1,self%npsi + CALL spline_eval(R_spline,self%x(i),0) + qtmp(i) = R_spline%f(1)*R_spline%f(2) + END DO + CALL gs_flux_int(gseq,self%x,self%jphi/qtmp,self%npsi,jphi_norm) + DEALLOCATE(qtmp) + jphi_norm=ABS(gseq%Itor_target)/jphi_norm + self%norm_last=jphi_norm + ! WRITE(*,*)'Ip flux: ',jphi_norm +ELSE + CALL gs_itor_nl(gseq,jphi_norm) + jphi_norm=(ABS(gseq%Itor_target)*mu0/jphi_norm + self%norm_last)/2.d0 + self%norm_last=jphi_norm + ! WRITE(*,*)'Ip NL: ',jphi_norm +END IF +! jphi_norm=ABS(gseq%Itor_target)/jphi_norm +! WRITE(*,*)'Ip flux: ',jphi_norm !---Get pressure profile CALL gseq%P%update(gseq) ! Make sure pressure profile is up to date with EQ IF(ASSOCIATED(gseq%P_ani))CALL oft_abort('Jphi profiles do not support anistopic pressure','jphi_update',__FILE__) !CALL gseq%P_ani%update(gseq) From e8b0e5f5055e25ba445c90b284a3caa5e87a12bf Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 15 Apr 2026 17:51:20 -0400 Subject: [PATCH 02/13] Bootstrap hot fixes: edge spike detection, verbose flags, and bug fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Consolidates post-PR#195 improvements to bootstrap.py: - Add prominence filter to edge spike detection - Fix degenerate curve_fit bounds in analyze_bootstrap_edge_spike - Tighten Ip scale factor bracket from ±10% to ±1% - Add verbose flag to solve_with_bootstrap and find_optimal_scale - Move log Lambda calc inside relevant logic - Various bug fixes and diagnostics Co-Authored-By: Claude Opus 4.6 (1M context) --- .../OpenFUSIONToolkit/TokaMaker/bootstrap.py | 144 +++++++++++------- 1 file changed, 91 insertions(+), 53 deletions(-) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py index c3f07c9ea..f9027d9db 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py @@ -207,15 +207,23 @@ def analyze_bootstrap_edge_spike(psi_N, j_bootstrap, diagnostic_plots=False): psi_edge = psi_N[edge_mask] j_edge = j_bootstrap[edge_mask] - # Find peak in the edge region - peaks, properties = find_peaks(j_edge, height=0.)#0.5*numpy.max(j_edge)) + # Find peak in the edge region (prominence filter suppresses noise + # oscillations that can be misidentified for small bootstrap spikes) + min_prominence = 0.05 * numpy.max(j_edge) if numpy.max(j_edge) > 0 else 0.0 + peaks, properties = find_peaks(j_edge, height=0., prominence=min_prominence) if len(peaks) == 0: print("No clear peak found in the edge region") return None - # Choose peak closest to psi_N = 1 if multiple peaks exist - peak_idx = peaks[numpy.argmax(psi_edge[peaks])] + # Among peaks in the far edge (psi_N > 0.85), pick the tallest; + # fall back to rightmost peak if none are in that zone. + far_edge = psi_edge[peaks] > 0.85 + if numpy.any(far_edge): + far_peaks = peaks[far_edge] + peak_idx = far_peaks[numpy.argmax(j_edge[far_peaks])] + else: + peak_idx = peaks[numpy.argmax(psi_edge[peaks])] peak_psi = psi_edge[peak_idx] peak_height = j_edge[peak_idx] @@ -293,21 +301,35 @@ def analyze_bootstrap_edge_spike(psi_N, j_bootstrap, diagnostic_plots=False): # Attempt to fit bootstrap parameterization to calculated profile (deprecated) # Initial parameter guess sigma_init = fwhm/2.355 # Convert FWHM to sigma - p0 = [peak_height, peak_psi, sigma_init, lmin_j_BS, 1.0, j_bootstrap[-1], 0.05] - - lower_bounds = [0.99*peak_height, # fix to measured spike height - 0.8*peak_psi, # don't fix to measured spike location - 0.0, # width - 0.99*lmin_j_BS, - -50, - 0.0, - 0.001] - upper_bounds = [1.01*peak_height, # fix to measured spike height - 1.2*peak_psi, # don't fix to measured spike location - 0.33, # no spikes wider than 0.33 units of psi_N - 1.01*lmin_j_BS, - 50, # sk - 2*j_bootstrap[-1], + p0 = [peak_height, peak_psi, sigma_init, lmin_j_BS, 1.0, + max(0.0, j_bootstrap[-1]), 0.05] + + # Guard against degenerate bounds when peak_height, lmin_j_BS, + # or j_bootstrap[-1] are near zero or negative (common for low- + # pedestal L-mode-like profiles in an LH transition scan). + _eps = max(1e-6, 1e-3 * abs(peak_height)) if peak_height != 0 else 1e-6 + amp_lo = min(0.99 * peak_height, peak_height - _eps) + amp_hi = max(1.01 * peak_height, peak_height + _eps) + + _eps_off = max(1e-6, 1e-3 * abs(lmin_j_BS)) if lmin_j_BS != 0 else 1e-6 + off_lo = min(0.99 * lmin_j_BS, lmin_j_BS - _eps_off) + off_hi = max(1.01 * lmin_j_BS, lmin_j_BS + _eps_off) + + ysep_hi = max(2 * abs(j_bootstrap[-1]), abs(j_bootstrap[-1]) + 1e-6, 1e-6) + + lower_bounds = [amp_lo, + 0.8*peak_psi, + 0.0, # width + off_lo, + -50, # sk + 0.0, # y_sep + 0.001] # blend_width + upper_bounds = [amp_hi, + 1.2*peak_psi, + 0.33, + off_hi, + 50, # sk + ysep_hi, 0.2] # Perform the fit @@ -359,9 +381,9 @@ def solve_jphi(mygs,ffp_prof,pp_prof,Ip_target,pax_target): # Solve Grad-Shafranov mygs.solve() -def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, - Ip_target, psi_pad, spike_prof=None, find_j0=True, scale_j0=1.0, - tolerance=0.01, max_iter=5, diagnostic_plots=False): +def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, + Ip_target, psi_pad, spike_prof=None, find_j0=True, scale_j0=1.0, + tolerance=0.01, max_iter=5, diagnostic_plots=False, verbose=True): r'''! Optimize scaling to match input/output \f$j_0\f$ or \f$I_p\f$ @param mygs Grad-Shafranov solver object @@ -389,7 +411,8 @@ def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, # We want: Input_J0 - Output_J0 ~ 0 def get_j0_error(scale_val, n): - print(f"\n--- Checking scale_j0 = {scale_val:.4f} ---") + if verbose: + print(f"\n--- Checking scale_j0 = {scale_val:.4f} ---") matched_input_jphi = scale_val*j_inductive + spike_prof ffp_prof['type'] = 'jphi-linterp' @@ -426,15 +449,17 @@ def get_j0_error(scale_val, n): diff = input_j0 - output_j0 rel_err = abs(diff) / output_j0 - print(f" Input j_0: {input_j0:.4e}") - print(f" Output j_0: {output_j0:.4e}") - print(f" Mismatch: {rel_err*100:.3f}%") + if verbose: + print(f" Input j_0: {input_j0:.4e}") + print(f" Output j_0: {output_j0:.4e}") + print(f" Mismatch: {rel_err*100:.3f}%") return diff, rel_err, tmp_jphi # We want: Input_Ip - Output_Ip ~ 0 def get_Ip_error(scale_val, scale_j0, n): - print(f"\n--- Checking scale_Ip = {scale_val:.4f} ---") + if verbose: + print(f"\n--- Checking scale_Ip = {scale_val:.4f} ---") ffp_prof['type'] = 'jphi-linterp' ffp_prof['y'] = scale_j0*j_inductive + spike_prof @@ -449,12 +474,16 @@ def get_Ip_error(scale_val, scale_j0, n): # Calculate residual (difference) and relative error diff = output_Ip - Ip_target - rel_err = abs(diff) / Ip_target + # Use magnitude of Ip_target and guard against zero to avoid sign/zero issues + eps = numpy.finfo(float).eps + denom = max(abs(Ip_target), eps) + rel_err = abs(diff) / denom - print(f" Input Ip target: {Ip_target/1e6:.4e}") - print(f" Trial Ip target: {Ip_target*scale_val/1e6:.4e}") - print(f" Output Ip: {output_Ip/1e6:.4e}") - print(f" Mismatch: {rel_err*100:.4f}%") + if verbose: + print(f" Input Ip target: {Ip_target/1e6:.4e}") + print(f" Trial Ip target: {Ip_target*scale_val/1e6:.4e}") + print(f" Output Ip: {output_Ip/1e6:.4e}") + print(f" Mismatch: {rel_err*100:.4f}%") return diff, rel_err, None @@ -478,11 +507,11 @@ def get_Ip_error(scale_val, scale_j0, n): p1 = 1.2 else: p1 = 0.8 - else: + else: if err0 < 0: - p1 = 1.1 + p1 = 1.01 else: - p1 = 0.9 + p1 = 0.99 if find_j0: err1, rel_err1, res_jphi = get_j0_error(p1, 1) @@ -495,11 +524,13 @@ def get_Ip_error(scale_val, scale_j0, n): # --- Step 3: Secant Method Loop --- # Iterate to find the root where Input - Output = 0 for i in range(max_iter): - print(f"--- Optimization Iteration {i+1} ---") + if verbose: + print(f"--- Optimization Iteration {i+1} ---") # Avoid division by zero if abs(err1 - err0) < 1e-9: - print("Error difference too small, stopping.") + if verbose: + print("Error difference too small, stopping.") break # Secant formula: x_new = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) @@ -520,14 +551,16 @@ def get_Ip_error(scale_val, scale_j0, n): err_new, rel_err_new, res_jphi = get_Ip_error(p_new, scale_j0, i+2) if rel_err_new < tolerance: - print(f"Converged! Optimal scale factor: {p_new:.4f}") + if verbose: + print(f"Converged! Optimal scale factor: {p_new:.4f}") return p_new, res_jphi # Update points for next step (move window forward) p0, err0 = p1, err1 p1, err1 = p_new, err_new - print("Max iterations reached. Returning best last effort.") + if verbose: + print("Max iterations reached. Returning best last effort.") return p1, res_jphi def calculate_ln_lambda(Te, Ti, ne, ni, Zeff=1.0, @@ -799,7 +832,8 @@ def solve_with_bootstrap(mygs, iterations=3, diagnostic_plots=False, parameterize_jBS = False, - use_OMFIT_sauter = False): + use_OMFIT_sauter = False, + verbose = True): r'''! Self-consistently compute bootstrap current from H-mode profiles @param mygs Grad-Shafranov solver object @@ -972,7 +1006,8 @@ def calculate_profiles_and_bootstrap(psi_N, include_jBS): method='brentq', rtol=1e-6) alpha_opt = sol.root except ValueError: - print("WARNING: Root scalar failed to bracket. Defaulting to alpha=1.0") + if verbose: + print("WARNING: Root scalar failed to bracket. Defaulting to alpha=1.0") alpha_opt = 1.0 matched_j_inductive = alpha_opt * current_jphi_target @@ -989,7 +1024,8 @@ def calculate_profiles_and_bootstrap(psi_N, include_jBS): if inductive_jphi is not None: - print('\n >>> Matching input core j_phi with G-S solution') + if verbose: + print('\n >>> Matching input core j_phi with G-S solution') # Calculate new profiles pp_prof, ffp_prof, j_bs_curr, matched_j_inductive, spike_prof = calculate_profiles_and_bootstrap( @@ -1015,24 +1051,26 @@ def calculate_profiles_and_bootstrap(psi_N, include_jBS): # Enforce P' edge condition pp_prof['y'][-1] = 0. - print('\n >>> Finding optimal j_phi scale factor') + if verbose: + print('\n >>> Finding optimal j_phi scale factor') # Find optimal jphi scale final_scale_j0, final_jphi = find_optimal_scale(mygs, - psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, - Ip_target, psi_pad, spike_prof=spike_prof, find_j0=True, - diagnostic_plots=diagnostic_plots + psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, + Ip_target, psi_pad, spike_prof=spike_prof, find_j0=True, + diagnostic_plots=diagnostic_plots, verbose=verbose ) - # final_scale_j0 = 1.0 - print('\n >>> Finding optimal Ip scale factor') + if verbose: + print('\n >>> Finding optimal Ip scale factor') # Find optimal Ip_target scale final_scale_Ip, _ = find_optimal_scale(mygs, - psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, - Ip_target, psi_pad, spike_prof=spike_prof, find_j0=False, - scale_j0=final_scale_j0, - tolerance=0.001, diagnostic_plots=diagnostic_plots + psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, + Ip_target, psi_pad, spike_prof=spike_prof, find_j0=False, + scale_j0=final_scale_j0, + tolerance=0.001, diagnostic_plots=diagnostic_plots, verbose=verbose ) - print('\n >>> Iterating on H-mode equilibrium solution') + if verbose: + print('\n >>> Iterating on H-mode equilibrium solution') for n in range(iterations): # Calculate new profiles From 7cb8f374d8ae530243a49dba2c8621001a347175 Mon Sep 17 00:00:00 2001 From: Chris Hansen Date: Tue, 5 May 2026 21:19:04 -0400 Subject: [PATCH 03/13] Add explicit variable to skip non-linear solver Ip target - Replace `Itor_target` with `Ip_target` throughout code --- src/physics/grad_shaf.F90 | 27 +++++++++-------- src/physics/grad_shaf_fit.F90 | 46 ++++++++++++++--------------- src/physics/grad_shaf_prof_phys.F90 | 22 +++++++------- src/physics/grad_shaf_td.F90 | 4 +-- src/physics/grad_shaf_util.F90 | 4 +-- 5 files changed, 54 insertions(+), 49 deletions(-) diff --git a/src/physics/grad_shaf.F90 b/src/physics/grad_shaf.F90 index c62f014ee..3803f12ed 100644 --- a/src/physics/grad_shaf.F90 +++ b/src/physics/grad_shaf.F90 @@ -263,6 +263,7 @@ MODULE oft_gs TYPE :: gs_equil LOGICAL :: diverted = .FALSE. !< Equilibrium is diverted? LOGICAL :: has_plasma = .TRUE. !< Solve with plasma? (otherwise vacuum) + LOGICAL :: Ip_target_skip = .FALSE. !< Skip toroidal current target in non-linear solve? INTEGER(i4) :: mode = 0 !< RHS source mode (0 -> F*F', 1 -> F') INTEGER(i4) :: nx_points = 0 !< Number of X-points in current solution INTEGER(i4) :: nregularize = 0 !< Number of regularization terms @@ -277,7 +278,7 @@ MODULE oft_gs REAL(r8) :: mirror_n = -1.d0 !< Anisotropy exponent for mirror pressure profiles REAL(r8) :: mirror_bturn = 0.d0 !< Turning point for mirror pressure profiles REAL(r8) :: mirror_zthroat = 0.d0 !< Mirror peak field point - REAL(r8) :: Itor_target = -1.d0 !< Toroidal current target + REAL(r8) :: Ip_target = -1.d0 !< Toroidal current target REAL(r8) :: estore_target = -1.d0 !< Stored energy target REAL(r8) :: pax_target = -1.d0 !< On-axis pressure target REAL(r8) :: Ip_ratio_target = -1.d99 !< Ip ratio target @@ -1046,7 +1047,8 @@ subroutine copy_eq(self,source) self%mirror_n=source%mirror_n self%mirror_bturn=source%mirror_bturn self%mirror_zthroat=source%mirror_zthroat -self%Itor_target=source%Itor_target +self%Ip_target_skip=source%Ip_target_skip +self%Ip_target=source%Ip_target self%estore_target=source%estore_target self%pax_target=source%pax_target self%Ip_ratio_target=source%Ip_ratio_target @@ -1114,13 +1116,13 @@ subroutine gs_init_psi(self,equil,ierr,r0,a,kappa,delta,curr_source) CALL equil%psi%set(0.d0) CALL gs_vacuum_solve(self,equil%psi,tmp_vec) itor=equil%itor() - IF(equil%Itor_target>0.d0)CALL equil%psi%scale(equil%Itor_target/itor) + IF(equil%Ip_target>0.d0)CALL equil%psi%scale(equil%Ip_target/itor) ELSE IF(PRESENT(curr_source))THEN CALL equil%psi%set(0.d0) CALL tmp_vec%restore_local(curr_source) CALL gs_vacuum_solve(self,equil%psi,tmp_vec) itor=equil%itor() - IF(equil%Itor_target>0.d0)CALL equil%psi%scale(equil%Itor_target/itor) + IF(equil%Ip_target>0.d0)CALL equil%psi%scale(equil%Ip_target/itor) ELSE !---Setup Solver eigsolver%A=>self%dels @@ -1146,9 +1148,9 @@ subroutine gs_init_psi(self,equil,ierr,r0,a,kappa,delta,curr_source) call equil%psi%scale(1.d0/psi_vals(mind)) equil%psimax=1.d0 END IF - IF(equil%Itor_target>0.d0)THEN + IF(equil%Ip_target>0.d0)THEN itor=equil%itor() - CALL equil%psi%scale(equil%Itor_target/itor) + CALL equil%psi%scale(equil%Ip_target/itor) END IF !---Cleanup CALL eigsolver%pre%delete @@ -2019,6 +2021,7 @@ subroutine gs_solve(self,equil,ierr) !--- error_flag=0 self%nl_its=0 +self%Ip_target_skip=.FALSE. IF(TRIM(self%lu_solver%package)=='none')THEN CALL oft_abort("LU solver required for GS solve","gs_solve",__FILE__) ELSE @@ -2177,9 +2180,9 @@ subroutine gs_solve(self,equil,ierr) param_mat(1,1)=1.d0 param_rhs(1)=0.d0 ELSE - IF(equil%Itor_target>0.d0)THEN + IF((equil%Ip_target>0.d0).AND.(.NOT.equil%Ip_target_skip))THEN param_mat(1,:)=[itor_ffp,itor_press,0.d0] - param_rhs(1)=equil%Itor_target + param_rhs(1)=equil%Ip_target ELSE param_mat(1,1)=1.d0 param_rhs(1)=equil%ffp_scale @@ -2414,13 +2417,13 @@ subroutine gs_solve(self,equil,ierr) CALL equil%psi%scale(1.d0/equil%psimax) equil%ffp_scale=SQRT((equil%ffp_scale**2)/equil%psimax) equil%psimax=1.d0 - ! ELSE IF(self%Itor_target>0.d0)THEN + ! ELSE IF(self%Ip_target>0.d0)THEN ! itor=self%itor() ! IF(itor<=0.d0)THEN ! error_flag=-5 ! EXIT ! END IF - ! ffp_scale_prev=itor/self%Itor_target + ! ffp_scale_prev=itor/self%Ip_target ! CALL self%psi%scale(1.d0/ffp_scale_prev) ! self%ffp_scale=SQRT((self%ffp_scale**2)/ffp_scale_prev) END IF @@ -2602,11 +2605,11 @@ subroutine gs_lin_solve(self,equil,adjust_r0,ierr) param_mat=0.d0 param_rhs=0.d0 -IF(equil%Itor_target>0.d0)THEN +IF(equil%Ip_target>0.d0)THEN itor_ffp=equil%itor(psi_ffp) itor_press=equil%itor(psi_press) param_mat(1,:)=[itor_ffp,itor_press,0.d0] - param_rhs(1)=equil%Itor_target + param_rhs(1)=equil%Ip_target ELSE param_mat(1,1)=1.d0 param_rhs(1)=equil%ffp_scale diff --git a/src/physics/grad_shaf_fit.F90 b/src/physics/grad_shaf_fit.F90 index 8ef25a41c..1622f0a7e 100644 --- a/src/physics/grad_shaf_fit.F90 +++ b/src/physics/grad_shaf_fit.F90 @@ -267,7 +267,7 @@ SUBROUTINE fit_gs(gs,inpath,outpath,fitI,fitP,fit_p_scale,fit_ffp_scale,fitR0,fi !---Count coefficients ncofs=0 IF(gs%device%free)THEN - IF(fit_FFPscale.OR.(gs_active%Itor_target>0.d0))ncofs = ncofs+1 + IF(fit_FFPscale.OR.(gs_active%Ip_target>0.d0))ncofs = ncofs+1 ELSE ncofs=ncofs+1 IF(fit_FFPscale)CALL oft_abort('Lambda cannot be fit in fixed boundary mode.', & @@ -292,14 +292,14 @@ SUBROUTINE fit_gs(gs,inpath,outpath,fitI,fitP,fit_p_scale,fit_ffp_scale,fitR0,fi IF(fit_FFPscale)THEN offset=1 cofs(1)=gs_active%ffp_scale - ELSE IF(gs_active%Itor_target>0.d0)THEN + ELSE IF(gs_active%Ip_target>0.d0)THEN offset=1 - cofs(1)=gs_active%Itor_target + cofs(1)=gs_active%Ip_target END IF ELSE offset=1 - IF(gs_active%Itor_target>0.d0)THEN - cofs(1)=gs_active%Itor_target + IF(gs_active%Ip_target>0.d0)THEN + cofs(1)=gs_active%Ip_target ELSE cofs(1)=gs_active%psiscale cofs_scale(1)=1.d0/ABS(gs_active%psiscale) @@ -409,7 +409,7 @@ SUBROUTINE fit_gs(gs,inpath,outpath,fitI,fitP,fit_p_scale,fit_ffp_scale,fitR0,fi geval_count=0 !---Initialize CALL fit_error(ncons,ncofs,cofs,error,info) - ! gs_active%Itor_target=-1.d0 + ! gs_active%Ip_target=-1.d0 ! IF(fit_FFPscale)cofs(1)=gs_active%ffp_scale IF(gs_active%device%ierr<0)CALL oft_abort('Initial equilibrium solve failed to converge','fit_gs',__FILE__) !--- @@ -571,7 +571,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) ALLOCATE(cofs_in(n)) cofs_in=cofs ffp_scale_in=gs_active%ffp_scale - ip_target_in=gs_active%Itor_target + ip_target_in=gs_active%Ip_target p_scale_in=gs_active%p_scale estore_target_in=gs_active%estore_target ! bounds_in=gs_active%spatial_bounds @@ -605,14 +605,14 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) IF(fit_FFPscale)THEN offset=1 gs_active%ffp_scale=cofs(1) - ELSE IF(gs_active%Itor_target>0.d0)THEN + ELSE IF(gs_active%Ip_target>0.d0)THEN offset=1 - gs_active%Itor_target=cofs(1) + gs_active%Ip_target=cofs(1) END IF ELSE offset=1 - IF(gs_active%Itor_target>0.d0)THEN - gs_active%Itor_target=cofs(1) + IF(gs_active%Ip_target>0.d0)THEN + gs_active%Ip_target=cofs(1) ELSE gs_active%psiscale=cofs(1) END IF @@ -677,7 +677,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) IF(fit_F0)gs_active%I%f_offset = cofs(offset+1) ! !---Centering ! ffp_scale_in=gs_active%ffp_scale - ! ip_target_in=gs_active%Itor_target + ! ip_target_in=gs_active%Ip_target ! p_scale_in=gs_active%p_scale ! estore_target_in=gs_active%estore_target ! bounds_in=gs_active%spatial_bounds @@ -713,7 +713,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) CALL psi_best%add(0.d0,1.d0,gs_active%psi) END IF ffp_scale_in=gs_active%ffp_scale - ip_target_in=gs_active%Itor_target + ip_target_in=gs_active%Ip_target p_scale_in=gs_active%p_scale estore_target_in=gs_active%estore_target ! bounds_in=gs_active%spatial_bounds @@ -737,13 +737,13 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) IF(gs_active%device%free)THEN IF(fit_FFPscale)THEN offset=offset+1 - ELSE IF(gs_active%Itor_target>0.d0)THEN - WRITE(*,'(2A,ES11.3)')oft_indent,'Itor_target =',gs_active%Itor_target/mu0 + ELSE IF(gs_active%Ip_target>0.d0)THEN + WRITE(*,'(2A,ES11.3)')oft_indent,'Ip_target =',gs_active%Ip_target/mu0 offset=offset+1 END IF ELSE - IF(gs_active%Itor_target>0.d0)THEN - WRITE(*,'(2A,ES11.3)')oft_indent,'Itor_target =',gs_active%Itor_target/mu0 + IF(gs_active%Ip_target>0.d0)THEN + WRITE(*,'(2A,ES11.3)')oft_indent,'Ip_target =',gs_active%Ip_target/mu0 ELSE WRITE(*,'(2A,ES11.3)')oft_indent,'Psi_scale =',gs_active%psiscale END IF @@ -826,18 +826,18 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) jac_mat(:,offset+1)=(jac_mat(:,offset+1)-err)/dx gs_active%ffp_scale=cofs(offset+1) offset=1 - ELSE IF(gs_active%Itor_target>0.d0)THEN + ELSE IF(gs_active%Ip_target>0.d0)THEN CALL reset_eq dx = dxi/cofs_scale(offset+1) - gs_active%Itor_target=cofs(offset+1) + dx + gs_active%Ip_target=cofs(offset+1) + dx CALL run_err(.FALSE.,jac_mat(:,offset+1),m,ierr) jac_mat(:,offset+1)=(jac_mat(:,offset+1)-err)/dx - gs_active%Itor_target=cofs(offset+1) + gs_active%Ip_target=cofs(offset+1) offset=1 END IF ELSE - IF(gs_active%Itor_target>0.d0)THEN - gs_active%Itor_target=cofs(1) + IF(gs_active%Ip_target>0.d0)THEN + gs_active%Ip_target=cofs(1) ELSE CALL reset_eq dx = dxi/cofs_scale(offset+1) @@ -975,7 +975,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) ! SUBROUTINE reset_eq gs_active%ffp_scale=ffp_scale_in -gs_active%Itor_target=ip_target_in +gs_active%Ip_target=ip_target_in gs_active%p_scale=p_scale_in gs_active%estore_target=estore_target_in ! gs_active%spatial_bounds=bounds_in diff --git a/src/physics/grad_shaf_prof_phys.F90 b/src/physics/grad_shaf_prof_phys.F90 index 892c1ccca..c998b108c 100644 --- a/src/physics/grad_shaf_prof_phys.F90 +++ b/src/physics/grad_shaf_prof_phys.F90 @@ -501,6 +501,7 @@ subroutine jphi_copy(self,new) new%f_offset=self%f_offset new%ngeom = self%ngeom new%j0 = self%j0 + new%norm_last = self%norm_last ALLOCATE(new%jphi, SOURCE=self%jphi) END SELECT end subroutine jphi_copy @@ -527,7 +528,7 @@ subroutine jphi_update(self,gseq) type(spline_type) :: R_spline self%plasma_bounds=gseq%plasma_bounds IF(gseq%mode/=1)CALL oft_abort("Jphi profile requires (F^2)' formulation","jphi_update",__FILE__) -! IF(gseq%Itor_target<0.d0)CALL oft_abort("Jphi profile requires Ip target","jphi_update",__FILE__) +IF(gseq%Ip_target<0.d0)CALL oft_abort("Jphi profile requires Ip target","jphi_update",__FILE__) IF(gseq%pax_target<0.d0)CALL oft_abort("Jphi profile requires Pax target","jphi_update",__FILE__) !---Get updated flux surface geometry for Jphi -> F*F' mapping ALLOCATE(ravgs(self%ngeom+1,3),psi_q(self%ngeom+1),qtmp(self%ngeom+1)) @@ -551,7 +552,12 @@ subroutine jphi_update(self,gseq) CALL spline_fit(R_spline,"extrap") DEALLOCATE(ravgs,psi_q,qtmp) !---Update jphi normalization to match Ip target -IF(gseq%Itor_target>0.d0)THEN +IF(gseq%Ip_target_skip)THEN + CALL gs_itor_nl(gseq,jphi_norm) + jphi_norm=(ABS(gseq%Ip_target)*mu0/jphi_norm + self%norm_last)/2.d0 + self%norm_last=jphi_norm + ! WRITE(*,*)'Ip NL: ',jphi_norm +ELSE ALLOCATE(qtmp(self%npsi)) DO i=1,self%npsi CALL spline_eval(R_spline,self%x(i),0) @@ -559,16 +565,11 @@ subroutine jphi_update(self,gseq) END DO CALL gs_flux_int(gseq,self%x,self%jphi/qtmp,self%npsi,jphi_norm) DEALLOCATE(qtmp) - jphi_norm=ABS(gseq%Itor_target)/jphi_norm + jphi_norm=ABS(gseq%Ip_target)/jphi_norm self%norm_last=jphi_norm ! WRITE(*,*)'Ip flux: ',jphi_norm -ELSE - CALL gs_itor_nl(gseq,jphi_norm) - jphi_norm=(ABS(gseq%Itor_target)*mu0/jphi_norm + self%norm_last)/2.d0 - self%norm_last=jphi_norm - ! WRITE(*,*)'Ip NL: ',jphi_norm END IF -! jphi_norm=ABS(gseq%Itor_target)/jphi_norm +! jphi_norm=ABS(gseq%Ip_target)/jphi_norm ! WRITE(*,*)'Ip flux: ',jphi_norm !---Get pressure profile CALL gseq%P%update(gseq) ! Make sure pressure profile is up to date with EQ @@ -585,7 +586,8 @@ subroutine jphi_update(self,gseq) self%yp(i) = 2.d0*(self%jphi(i)*jphi_norm - R_spline%f(1)*pprime*pscale)/R_spline%f(2) END DO ! Disable Ip matching and fix F*F' scale (matching is done here instead) -IF(gseq%Itor_target>0.d0)gseq%Itor_target=-gseq%Itor_target +gseq%Ip_target_skip=.TRUE. +! IF(gseq%Ip_target>0.d0)gseq%Ip_target=-gseq%Ip_target gseq%ffp_scale=1.d0 !---Clean up CALL spline_dealloc(R_spline) diff --git a/src/physics/grad_shaf_td.F90 b/src/physics/grad_shaf_td.F90 index 1031ca9b8..d21224148 100644 --- a/src/physics/grad_shaf_td.F90 +++ b/src/physics/grad_shaf_td.F90 @@ -603,7 +603,7 @@ subroutine setup_mfop(self,eq_in) DEBUG_STACK_PUSH self%gs_equil=>eq_in self%gs_device=>eq_in%device -self%ip_target=self%gs_equil%Itor_target +self%ip_target=self%gs_equil%Ip_target self%ip_ratio_target=self%gs_equil%Ip_ratio_target self%f_scale=self%gs_equil%ffp_scale self%p_scale=self%gs_equil%p_scale @@ -632,7 +632,7 @@ subroutine update_mfop(self) self%F=>self%gs_equil%I self%P=>self%gs_equil%P ! Update current target and sync scale factors -self%ip_target=self%gs_equil%Itor_target +self%ip_target=self%gs_equil%Ip_target self%ip_ratio_target=self%gs_equil%Ip_ratio_target self%f_scale=self%gs_equil%ffp_scale self%p_scale=self%gs_equil%p_scale diff --git a/src/physics/grad_shaf_util.F90 b/src/physics/grad_shaf_util.F90 index b95e529fe..5e5346b72 100644 --- a/src/physics/grad_shaf_util.F90 +++ b/src/physics/grad_shaf_util.F90 @@ -389,7 +389,7 @@ subroutine gs_save_tokamaker(self,filename) CALL hdf5_write(self%lim_point,filename,'tokamaker/LIM_POINT') CALL hdf5_write(self%diverted,filename,'tokamaker/DIVERTED') !--- -IF(self%Itor_target>0.d0)CALL hdf5_write(self%Itor_target,filename,'tokamaker/IP_TARGET') +IF(self%Ip_target>0.d0)CALL hdf5_write(self%Ip_target,filename,'tokamaker/IP_TARGET') IF(self%Ip_ratio_target>-1.d98)CALL hdf5_write(self%Ip_ratio_target,filename,'tokamaker/IP_RATIO_TARGET') IF(self%R0_target>0.d0)CALL hdf5_write(self%R0_target,filename,'tokamaker/R0_TARGET') IF(self%Z0_target>-1.d98)CALL hdf5_write(self%Z0_target,filename,'tokamaker/Z0_TARGET') @@ -608,7 +608,7 @@ subroutine gs_load_tokamaker(self,filename,error_string) END IF !---Load targets and coil currents IF(hdf5_field_exist(filename,'tokamaker/IP_TARGET'))THEN - CALL hdf5_read(self%Itor_target,filename,'tokamaker/IP_TARGET',success=success) + CALL hdf5_read(self%Ip_target,filename,'tokamaker/IP_TARGET',success=success) IF(.NOT.success)THEN error_string='Failed to read Ip target.' RETURN From d340ca876df33c4e18458ab0a9742746b3231abe Mon Sep 17 00:00:00 2001 From: Chris Hansen Date: Tue, 5 May 2026 21:24:25 -0400 Subject: [PATCH 04/13] Fix typos in prior commit --- src/physics/grad_shaf.F90 | 2 +- src/python/wrappers/tokamaker_f.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/grad_shaf.F90 b/src/physics/grad_shaf.F90 index 3803f12ed..6637b2f9f 100644 --- a/src/physics/grad_shaf.F90 +++ b/src/physics/grad_shaf.F90 @@ -2021,7 +2021,7 @@ subroutine gs_solve(self,equil,ierr) !--- error_flag=0 self%nl_its=0 -self%Ip_target_skip=.FALSE. +equil%Ip_target_skip=.FALSE. IF(TRIM(self%lu_solver%package)=='none')THEN CALL oft_abort("LU solver required for GS solve","gs_solve",__FILE__) ELSE diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index c6614edb5..80237faa8 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -1414,7 +1414,7 @@ SUBROUTINE tokamaker_set_targets(tMaker_ptr,ip_target,ip_ratio_target,pax_target tMaker_obj%gs_equil%Z0_target=Z0_target tMaker_obj%gs_equil%pax_target=pax_target*mu0 tMaker_obj%gs_equil%estore_target=estore_target*mu0 -tMaker_obj%gs_equil%itor_target=ip_target*mu0 +tMaker_obj%gs_equil%Ip_target=ip_target*mu0 tMaker_obj%gs_equil%ip_ratio_target=ip_ratio_target END SUBROUTINE tokamaker_set_targets !--------------------------------------------------------------------------------- From 91c0ac5870b7952dc0f609790fc4a21f28db9c25 Mon Sep 17 00:00:00 2001 From: Chris Hansen Date: Tue, 5 May 2026 21:29:10 -0400 Subject: [PATCH 05/13] Add IPython artifacts to .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 8ecd1dacb..77ea7de1d 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ # Python artifacts *.pyc __pycache__ +.ipynb_checkpoints # Editor Items .settings From 5b06f66fc7d6da70090f8689b55978cda31a742b Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 27 May 2026 22:20:32 -0400 Subject: [PATCH 06/13] feat(experimental): jphi-linterp Ip-correction outer iteration Address the well-known ~0.5-1% Ip discretization mismatch between jphi_update's flux-averaged normalization (gs_flux_int) and gs_comp_globals' physical-Ip integration (R*P' + 0.5*FFP/R). Disagreement is structural -- the two formulas integrate different quantities, agreeing only when *<1/R> = 1 pointwise (circular plasma). For D-shaped DIII-D plasma, ~0.76% mismatch. Fix: outer iteration in tokamaker_solve wrapper, gated on jphi_flux_func profile type detection. After each inner Picard solve, compute Ip_phys via gs_comp_globals, compare to original Itor_target, scale Itor_target by a damped multiplicative correction, re-Picard. Continues until rel_err < 5e-4 or max_outer (5) iterations. Key design choices: * jphi_update remains stateless from caller perspective -- avoids breaking caller-side iterative wrappers (eg. bouquet's find_optimal_scale, _corrective_jphi_iteration) that depend on solve being a deterministic FFP-from-jphi mapping for a given equilibrium * Damped correction (^0.7) + hard clamp [0.7, 1.3] for safety * Only activates for jphi_flux_func profiles -- PP'/FF' users unaffected * Skipped for vacuum solves (no Ip target) Validation (DIII-D 204441 @ 4.4s recon): * Ip error: -0.7570% -> -0.0396% (19x improvement) * PIN_JPHI draw 1 bnd-diag post-Ip-align: 3.68 mm -> 0.00 mm * Recon's [Ip match] secant converges identically (deterministic callable behaviour preserved) * No build regressions Co-Authored-By: Claude Opus 4.7 (1M context) --- src/python/wrappers/tokamaker_f.F90 | 95 +++++++++++++++++++++++++++-- 1 file changed, 91 insertions(+), 4 deletions(-) diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index d6936ccfe..a611a876d 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -37,7 +37,8 @@ MODULE tokamaker_f sauter_fc, gs_calc_vloop, gs_save_tokamaker, gs_load_tokamaker USE oft_gs_fit, ONLY: fit_gs, fit_pm USE oft_gs_td, ONLY: oft_tmaker_td, eig_gs_td -USE grad_shaf_prof_phys, ONLY: create_dipole_b0_prof, dipole_ani_press, mirror_ani_slosh +USE grad_shaf_prof_phys, ONLY: create_dipole_b0_prof, dipole_ani_press, mirror_ani_slosh, & + jphi_flux_func USE diagnostic, ONLY: bscal_surf_int USE oft_base_f, ONLY: copy_string, copy_string_rev, oftpy_init IMPLICIT NONE @@ -478,8 +479,18 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka INTEGER(c_int), INTENT(out) :: nl_its !< Number of nonlinear iterations CHARACTER(KIND=c_char), INTENT(out) :: error_str(OFT_ERROR_SLEN) !< Error string (empty if no error) INTEGER(i4) :: ntargets,ierr -LOGICAL :: vac_save +LOGICAL :: vac_save,is_jphi_linterp TYPE(tokamaker_instance), POINTER :: tMaker_obj +! ---- jphi-linterp Ip-correction outer-iteration state ---- +INTEGER(i4) :: outer_it +INTEGER(i4), PARAMETER :: jphi_ip_max_outer = 5 +REAL(r8), PARAMETER :: jphi_ip_tol = 5.0d-4 ! relative +REAL(r8), PARAMETER :: jphi_ip_damp = 0.7d0 ! correction^damp +REAL(r8), PARAMETER :: jphi_ip_min = 0.7d0 ! correction clamp +REAL(r8), PARAMETER :: jphi_ip_max = 1.3d0 ! correction clamp +REAL(r8) :: itor_target_orig, ip_phys, ip_phys_target, correction, rel_err +REAL(r8) :: cent_dummy(2), v_dummy, pv_dummy, df_dummy, tf_dummy, bv_dummy +INTEGER(i4) :: nl_its_total IF(.NOT.tokamaker_ccast(tMaker_ptr,tMaker_obj,error_str))RETURN IF(.NOT.tokamaker_require_equil(tMaker_obj,error_str))RETURN IF(ANY(tMaker_obj%device%rcoils>0.d0).AND.(tMaker_obj%device%dt>0.d0))THEN @@ -494,10 +505,86 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka tMaker_obj%gs_equil%has_plasma=.FALSE. END IF tMaker_obj%device%timing=0.d0 -CALL tMaker_obj%device%solve(tMaker_obj%gs_equil,ierr) + +! Detect whether the active FFP profile is jphi-linterp. Only that +! profile type has the cut-cell + flux-averaging Ip discretization +! mismatch that this outer loop corrects. PP'/FF' profiles set +! directly do not need this loop. +is_jphi_linterp = .FALSE. +IF(ASSOCIATED(tMaker_obj%gs_equil%I))THEN + SELECT TYPE(I_prof => tMaker_obj%gs_equil%I) + CLASS IS(jphi_flux_func) + is_jphi_linterp = .TRUE. + END SELECT +END IF + +! Single solve when not jphi-linterp (preserves original behaviour +! exactly for PP'/FF' users). Also when vacuum=True (no Ip target). +IF(.NOT.is_jphi_linterp .OR. vacuum)THEN + CALL tMaker_obj%device%solve(tMaker_obj%gs_equil,ierr) + IF(vacuum)tMaker_obj%gs_equil%has_plasma=vac_save + IF(ierr/=0)CALL copy_string(gs_err_reason(ierr),error_str) + nl_its=tMaker_obj%device%nl_its + RETURN +END IF + +! ---- jphi-linterp outer iteration ---- +! TokaMaker's jphi-linterp profile type normalizes FFP via a flux- +! surface-averaged integral (gs_flux_int of jphi/(·<1/R>)), but +! gs_comp_globals reports Ip via the *physical* integral (R·P' + +! 0.5·FFP/R over plasma area). The two disagree by ~0.5-1% on D- +! shaped plasmas due to ·<1/R> ≠ R·(1/R) pointwise. Result: a +! call with set_targets(Ip=X) lands at get_globals returning X·(1-ε). +! +! Fix: iterate the inner Picard with progressively scaled Itor_target +! until get_globals matches the user's nominal target within +! jphi_ip_tol. jphi_update remains stateless from the caller's +! perspective -- only Itor_target is rescaled at the OFT solve level. +! +! ip_phys_target stays at the user's original ask; itor_target is +! temporarily inflated to compensate for the per-iter mismatch. At +! convergence the inflation factor stops growing and ip_phys ~ user's +! original target. + +! jphi_update flips sign of Itor_target after first call; use ABS for +! consistent target tracking. +itor_target_orig = ABS(tMaker_obj%gs_equil%Itor_target) +ip_phys_target = itor_target_orig +nl_its_total = 0 +ierr = 0 +DO outer_it = 1, jphi_ip_max_outer + CALL tMaker_obj%device%solve(tMaker_obj%gs_equil,ierr) + nl_its_total = nl_its_total + tMaker_obj%device%nl_its + IF(ierr /= 0) EXIT ! inner Picard failed; bail with whatever we have + ! Compute physical Ip from current equilibrium + CALL gs_comp_globals(tMaker_obj%gs_equil, ip_phys, cent_dummy, & + v_dummy, pv_dummy, df_dummy, tf_dummy, bv_dummy) + ! ip_phys here has same internal units as Itor_target (both mu0*Ip) + IF(ABS(ip_phys) <= 0.d0)THEN + EXIT ! degenerate, can't form ratio + END IF + rel_err = (ABS(ip_phys) - ip_phys_target) / ip_phys_target + IF(ABS(rel_err) < jphi_ip_tol) EXIT ! converged + ! Damped multiplicative correction toward target + correction = (ip_phys_target / ABS(ip_phys)) ** jphi_ip_damp + ! Hard clamp for safety + IF(correction > jphi_ip_max) correction = jphi_ip_max + IF(correction < jphi_ip_min) correction = jphi_ip_min + ! Apply correction to Itor_target (preserve sign convention -- jphi_update + ! flips Itor_target negative; ABS-based scaling preserves that). + tMaker_obj%gs_equil%Itor_target = SIGN(ABS(tMaker_obj%gs_equil%Itor_target) * correction, & + tMaker_obj%gs_equil%Itor_target) +END DO + +IF(ierr == 0 .AND. ABS(rel_err) >= jphi_ip_tol .AND. outer_it > jphi_ip_max_outer)THEN + ! Did not converge within max_outer; not an error, but worth noting + IF(oft_debug_print(1))WRITE(*,*) & + ' tokamaker_solve: jphi-linterp Ip-corr did not converge, final rel_err =', rel_err +END IF + IF(vacuum)tMaker_obj%gs_equil%has_plasma=vac_save IF(ierr/=0)CALL copy_string(gs_err_reason(ierr),error_str) -nl_its=tMaker_obj%device%nl_its +nl_its = nl_its_total END SUBROUTINE tokamaker_solve !--------------------------------------------------------------------------------- !> Perform linear solve to find vacuum solution for given BCs and current sources From 328163f417dcd32217d27bddb3db3291b466e3f8 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 27 May 2026 23:02:50 -0400 Subject: [PATCH 07/13] fix(experimental): restore Itor_target at exit of jphi-linterp outer loop User-reported regression: with the outer loop in place, recon flow showed ~2% Ip drift after many internal solves. Cause: outer loop scales gseq%Itor_target multiplicatively but never restores it. Each subsequent mygs.solve() call starts from the previous call's final (inflated) Itor_target -- cumulative inflation across the ~dozens of solves a typical recon flow performs. Fix: save Itor_target magnitude at entry, restore at exit (preserving whatever sign jphi_update last set). The FFP is calibrated for the inflated target during the final inner Picard, so the equilibrium correctly has Ip ~ user_target; only the bookkeeping target is reset. Future mygs.solve() calls start from the clean user-set target and re-run the outer loop if needed. Validation (PIN_JPHI=1, sigma=0, N=3 draws): * Ip-align all 3 draws: actual_Ip == target (no-op corrections) * bnd-diag post-Ip-align: 0.00 mm (all draws) * bnd-diag post-homotopy: 0.39 mm (just iso-update + homotopy noise) * F-drift 0.00%, VSC drift 0.02% on every draw * 3/3 IN_SPEC * No cumulative drift across draws Before fix (original pre-outer-loop OFT): PIN_JPHI sigma=0: 3.69 mm RMS, Ip err -0.76% After fix (outer loop + restore): PIN_JPHI sigma=0: 0.39 mm RMS, Ip err +0.005% -- ~10x and 150x Co-Authored-By: Claude Opus 4.7 (1M context) --- src/python/wrappers/tokamaker_f.F90 | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index a611a876d..d3a4d6b5f 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -546,12 +546,17 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka ! convergence the inflation factor stops growing and ip_phys ~ user's ! original target. -! jphi_update flips sign of Itor_target after first call; use ABS for -! consistent target tracking. +! Save user's original Itor_target so we can restore it before returning. +! jphi_update flips sign of Itor_target inside the Picard, but the magnitude +! always equals the user's set_targets(Ip=...) value when this routine is +! entered. We need to restore the same magnitude (with whatever sign +! jphi_update last left it) so subsequent solves don't see compounding +! inflation from prior outer-loop scaling. itor_target_orig = ABS(tMaker_obj%gs_equil%Itor_target) ip_phys_target = itor_target_orig nl_its_total = 0 ierr = 0 +rel_err = 0.d0 DO outer_it = 1, jphi_ip_max_outer CALL tMaker_obj%device%solve(tMaker_obj%gs_equil,ierr) nl_its_total = nl_its_total + tMaker_obj%device%nl_its @@ -576,7 +581,17 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka tMaker_obj%gs_equil%Itor_target) END DO -IF(ierr == 0 .AND. ABS(rel_err) >= jphi_ip_tol .AND. outer_it > jphi_ip_max_outer)THEN +! CRITICAL: restore Itor_target to its original magnitude (preserving +! the sign jphi_update left it with). Otherwise the cumulative outer- +! loop scaling compounds across successive mygs.solve() calls, drifting +! Itor_target far from the user's set_targets value (observed: 2% drift +! over a recon flow with dozens of internal solves). The current FFP +! is calibrated for the inflated target, so the equilibrium correctly +! has Ip ~ user_target; only the bookkeeping target is reset. +tMaker_obj%gs_equil%Itor_target = SIGN(itor_target_orig, & + tMaker_obj%gs_equil%Itor_target) + +IF(ierr == 0 .AND. ABS(rel_err) >= jphi_ip_tol)THEN ! Did not converge within max_outer; not an error, but worth noting IF(oft_debug_print(1))WRITE(*,*) & ' tokamaker_solve: jphi-linterp Ip-corr did not converge, final rel_err =', rel_err From 97436b2c7291001860074526dddc2e954b6d7c20 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 27 May 2026 23:32:40 -0400 Subject: [PATCH 08/13] diag: opt-in trace prints for jphi-linterp outer loop (oft_debug_print) --- src/python/wrappers/tokamaker_f.F90 | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index d3a4d6b5f..3b6369660 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -557,10 +557,16 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka nl_its_total = 0 ierr = 0 rel_err = 0.d0 +IF(oft_debug_print(1))WRITE(*,'(A,ES16.8)') & + ' [JPHI_IP] outer-loop entry target_internal=', ip_phys_target DO outer_it = 1, jphi_ip_max_outer CALL tMaker_obj%device%solve(tMaker_obj%gs_equil,ierr) nl_its_total = nl_its_total + tMaker_obj%device%nl_its - IF(ierr /= 0) EXIT ! inner Picard failed; bail with whatever we have + IF(ierr /= 0)THEN + IF(oft_debug_print(1))WRITE(*,'(A,I0,A,I0)') & + ' [JPHI_IP] iter ', outer_it, ' inner solve FAILED ierr=', ierr + EXIT + END IF ! Compute physical Ip from current equilibrium CALL gs_comp_globals(tMaker_obj%gs_equil, ip_phys, cent_dummy, & v_dummy, pv_dummy, df_dummy, tf_dummy, bv_dummy) @@ -569,6 +575,9 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka EXIT ! degenerate, can't form ratio END IF rel_err = (ABS(ip_phys) - ip_phys_target) / ip_phys_target + IF(oft_debug_print(1))WRITE(*,'(A,I0,A,3ES16.8)') & + ' [JPHI_IP] iter ', outer_it, & + ' ip_phys / target / rel_err = ', ABS(ip_phys), ip_phys_target, rel_err IF(ABS(rel_err) < jphi_ip_tol) EXIT ! converged ! Damped multiplicative correction toward target correction = (ip_phys_target / ABS(ip_phys)) ** jphi_ip_damp From b032e4fab38a0636832961971c68efa79e649bb1 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 28 May 2026 16:58:42 -0400 Subject: [PATCH 09/13] fix(experimental): safety-bail jphi-linterp Ip outer loop on corrupted ip_phys Guard the multiplicative Itor_target correction against a corrupted ip_phys returned by gs_comp_globals after a nominally-successful (ierr=0) inner solve. Three failure modes, all observed under tight coil-bound homotopy at low j_phi pinning: ip_phys <= 0 (degenerate), ip_phys NaN/Inf (singular flux state), or |ip_phys| outside [0.3, 3.0]x target (QP went singular, plasma collapsed to device center, integral is bogus). Applying a correction from that state inflates Itor_target to nonsense and corrupts downstream solves; instead bail without correction and let the caller's rollback recover. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/python/wrappers/tokamaker_f.F90 | 58 ++++++++++++++++++++++++++++- 1 file changed, 57 insertions(+), 1 deletion(-) diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index 3b6369660..f0c4255c3 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -488,9 +488,23 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka REAL(r8), PARAMETER :: jphi_ip_damp = 0.7d0 ! correction^damp REAL(r8), PARAMETER :: jphi_ip_min = 0.7d0 ! correction clamp REAL(r8), PARAMETER :: jphi_ip_max = 1.3d0 ! correction clamp +! Sanity bounds on ip_phys / ip_phys_target -- if gs_comp_globals +! reports an Ip ratio outside [jphi_ip_ratio_lo, jphi_ip_ratio_hi] +! after a nominally-successful inner solve, the equilibrium has +! corrupted (e.g. axis collapsed to device center because a tight- +! bound QP went singular but still returned ierr=0). Applying the +! standard multiplicative correction in that state would inflate +! Itor_target by a meaningless factor and corrupt downstream solves. +! Bail without correction; the caller's exception/rollback handler +! will recover. The window is intentionally wide (0.3-3.0x) -- a +! healthy iter 1 lands at ~0.994x for typical D-shape, so 0.3x is +! ~50 standard deviations of the normal first-iter undershoot. +REAL(r8), PARAMETER :: jphi_ip_ratio_lo = 0.3d0 +REAL(r8), PARAMETER :: jphi_ip_ratio_hi = 3.0d0 REAL(r8) :: itor_target_orig, ip_phys, ip_phys_target, correction, rel_err REAL(r8) :: cent_dummy(2), v_dummy, pv_dummy, df_dummy, tf_dummy, bv_dummy INTEGER(i4) :: nl_its_total +LOGICAL :: ip_phys_bad IF(.NOT.tokamaker_ccast(tMaker_ptr,tMaker_obj,error_str))RETURN IF(.NOT.tokamaker_require_equil(tMaker_obj,error_str))RETURN IF(ANY(tMaker_obj%device%rcoils>0.d0).AND.(tMaker_obj%device%dt>0.d0))THEN @@ -571,8 +585,50 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka CALL gs_comp_globals(tMaker_obj%gs_equil, ip_phys, cent_dummy, & v_dummy, pv_dummy, df_dummy, tf_dummy, bv_dummy) ! ip_phys here has same internal units as Itor_target (both mu0*Ip) + ! + ! Sanity-check ip_phys before computing a correction factor. Three + ! failure modes are guarded against here, all of which have been + ! observed under tight coil-bound homotopy at low j_phi-pinning: + ! + ! (1) ABS(ip_phys) <= 0 -- degenerate; can't form ratio. (Original + ! guard, preserved.) + ! (2) ip_phys is NaN/Inf -- happens when the inner solve produced + ! a singular flux state that gs_flux_int couldn't integrate + ! cleanly. The Fortran idiom (x /= x) is the portable NaN + ! test; combined with the finite-range check below it also + ! catches Inf. + ! (3) ABS(ip_phys) far outside [jphi_ip_ratio_lo, jphi_ip_ratio_hi] + ! x target -- inner solve returned ierr=0 but the equilibrium is + ! corrupted (typically: QP at tight bounds went singular, + ! plasma collapsed to mesh boundary, ip_phys reflects the + ! bogus integral over a "no plasma" region). Applying a + ! correction factor computed from this rel_err would inflate + ! Itor_target to nonsense and the next iter's solve would also + ! fail. Better to bail now and let the caller's rollback (e.g. + ! bouquet's homotopy pass restore-and-resolve) recover. + ! + ! Observed example (2026-05 PIN_JPHI sweep, homotopy pass 3 + ! F=+/-1% VSC=+/-1%): inner solve returned ierr=0 but axis was + ! at (R, Z) = (0.5, 0.0) (device-center default for "no plasma"); + ! gs_comp_globals then reported a nonsense ip_phys, and the next + ! outer iter scaled Itor_target by the clamped jphi_ip_max=1.3 + ! before finally failing with ierr=-6. + ip_phys_bad = .FALSE. IF(ABS(ip_phys) <= 0.d0)THEN - EXIT ! degenerate, can't form ratio + ip_phys_bad = .TRUE. + ELSE IF(ip_phys /= ip_phys)THEN ! Fortran NaN test (NaN != NaN) + ip_phys_bad = .TRUE. + ELSE IF(ABS(ip_phys) < jphi_ip_ratio_lo * ip_phys_target)THEN + ip_phys_bad = .TRUE. + ELSE IF(ABS(ip_phys) > jphi_ip_ratio_hi * ip_phys_target)THEN + ip_phys_bad = .TRUE. + END IF + IF(ip_phys_bad)THEN + IF(oft_debug_print(1))WRITE(*,'(A,I0,A,2ES16.8)') & + ' [JPHI_IP] iter ', outer_it, & + ' SAFETY BAIL -- ip_phys out of sane range, ip_phys/target =', & + ABS(ip_phys), ip_phys_target + EXIT ! caller will see equilibrium as solver produced it; restore handled below END IF rel_err = (ABS(ip_phys) - ip_phys_target) / ip_phys_target IF(oft_debug_print(1))WRITE(*,'(A,I0,A,3ES16.8)') & From 8f91d72f2d51e069b3b3bc6543d191db25095d3f Mon Sep 17 00:00:00 2001 From: d-burg Date: Tue, 2 Jun 2026 18:35:47 -0400 Subject: [PATCH 10/13] refactor: drop redundant Ip-scale secant; find_optimal_scale is core-j0 only The GS solve now holds Ip to target natively (jphi-linterp cut-cell fix + Ip-correction outer loop), so the find_optimal_scale(find_j0=False) Ip-scale secant in solve_with_bootstrap is redundant -> final_scale_Ip = 1.0. Strip the now-unused find_j0=False branch, get_Ip_error, and the find_j0/scale_j0 parameters from find_optimal_scale; it is now a clean core-j0 matcher. Callers updated. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../OpenFUSIONToolkit/TokaMaker/bootstrap.py | 116 +++++------------- 1 file changed, 32 insertions(+), 84 deletions(-) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py index f9027d9db..0459b21a4 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py @@ -382,21 +382,19 @@ def solve_jphi(mygs,ffp_prof,pp_prof,Ip_target,pax_target): mygs.solve() def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, - Ip_target, psi_pad, spike_prof=None, find_j0=True, scale_j0=1.0, + Ip_target, psi_pad, spike_prof=None, tolerance=0.01, max_iter=5, diagnostic_plots=False, verbose=True): - r'''! Optimize scaling to match input/output \f$j_0\f$ or \f$I_p\f$ + r'''! Optimize core-\f$j_0\f$ scaling so input/output \f$j_0\f$ match @param mygs Grad-Shafranov solver object @param psi_N Normalized poloidal flux profile \f$\hat{\psi}\f$ @param pressure Plasma pressure profile \f$p(\hat{\psi})\f$ [Pa] @param ffp_prof Toroidal current profile (\f$j_\phi(\hat{\psi})\f$) - @param pp_prof Pressure gradient profile (\f{P'}\f$) + @param pp_prof Pressure gradient profile (\f{P\'}\f$) @param j_inductive Inductive/total toroidal current profile \f$j_{ind}(\hat{\psi})\f$ [A/m^2] @param Ip_target Target plasma current \f$I_p\f$ [A] @param psi_pad Padding for flux surface calculations @param spike_prof Optional bootstrap current profile \f$j_{BS}(\hat{\psi})\f$ [A/m^2] - @param find_j0 If True, rescale core \f$j_\phi\f$; else, rescale \f$I_p\f$ - @param scale_j0 Used if find_j0=False; rescaling factor for \f$j_0\f$ @param tolerance Relative error tolerance for secant method @param max_iter Maximum number of secant iterations @param diagnostic_plots If True, plot input/output \f$j_\phi\f$ @@ -405,15 +403,15 @@ def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, import matplotlib.pyplot as plt n_psi = len(psi_N) - + if spike_prof is None: spike_prof = numpy.zeros_like(j_inductive) - + # We want: Input_J0 - Output_J0 ~ 0 def get_j0_error(scale_val, n): if verbose: print(f"\n--- Checking scale_j0 = {scale_val:.4f} ---") - + matched_input_jphi = scale_val*j_inductive + spike_prof ffp_prof['type'] = 'jphi-linterp' ffp_prof['y'] = matched_input_jphi @@ -421,7 +419,7 @@ def get_j0_error(scale_val, n): pax_target = pressure[0] solve_jphi(mygs,ffp_prof,pp_prof,Ip_target,pax_target) - + # Check Convergence _, f, fp, _, pp = mygs.get_profiles(npsi=n_psi, psi_pad=psi_pad) _, _, ravgs, _, _, _ = mygs.get_q(npsi=n_psi, psi_pad=psi_pad) @@ -438,61 +436,27 @@ def get_j0_error(scale_val, n): plt.ylabel(r'$j_\phi$ [MA/m$^2$]') plt.grid(ls=':') plt.show() - + # Input j_0 at this scale input_j0 = scale_val * j_inductive[0] + spike_prof[0] - + # Output j_0 from solver output_j0 = tmp_jphi[0] - + # Calculate residual (difference) and relative error diff = input_j0 - output_j0 rel_err = abs(diff) / output_j0 - + if verbose: print(f" Input j_0: {input_j0:.4e}") print(f" Output j_0: {output_j0:.4e}") print(f" Mismatch: {rel_err*100:.3f}%") - - return diff, rel_err, tmp_jphi - - # We want: Input_Ip - Output_Ip ~ 0 - def get_Ip_error(scale_val, scale_j0, n): - if verbose: - print(f"\n--- Checking scale_Ip = {scale_val:.4f} ---") - - ffp_prof['type'] = 'jphi-linterp' - ffp_prof['y'] = scale_j0*j_inductive + spike_prof - scaled_Ip_target = Ip_target*scale_val - pax_target = pressure[0] - - solve_jphi(mygs,ffp_prof,pp_prof,scaled_Ip_target,pax_target) - - eq_stats = mygs.get_stats(lcfs_pad=psi_pad) - output_Ip = eq_stats['Ip'] - - # Calculate residual (difference) and relative error - diff = output_Ip - Ip_target - # Use magnitude of Ip_target and guard against zero to avoid sign/zero issues - eps = numpy.finfo(float).eps - denom = max(abs(Ip_target), eps) - rel_err = abs(diff) / denom - - if verbose: - print(f" Input Ip target: {Ip_target/1e6:.4e}") - print(f" Trial Ip target: {Ip_target*scale_val/1e6:.4e}") - print(f" Output Ip: {output_Ip/1e6:.4e}") - print(f" Mismatch: {rel_err*100:.4f}%") - - return diff, rel_err, None + return diff, rel_err, tmp_jphi # --- Step 1: Initial Guess (1.0) --- p0 = 1.0 - if find_j0: - err0, rel_err0, res_jphi = get_j0_error(p0, 0) - else: - err0, rel_err0, res_jphi = get_Ip_error(p0, scale_j0, 0) + err0, rel_err0, res_jphi = get_j0_error(p0, 0) # Check if we got lucky immediately if rel_err0 < tolerance: @@ -501,22 +465,12 @@ def get_Ip_error(scale_val, scale_j0, n): # --- Step 2: Second Guess (Directional) --- # If Input < Output (err0 < 0), we need more current -> Try 1.2 # If Input > Output (err0 > 0), we have too much -> Try 0.8 - - if find_j0: - if err0 < 0: - p1 = 1.2 - else: - p1 = 0.8 - else: - if err0 < 0: - p1 = 1.01 - else: - p1 = 0.99 - - if find_j0: - err1, rel_err1, res_jphi = get_j0_error(p1, 1) + if err0 < 0: + p1 = 1.2 else: - err1, rel_err1, res_jphi = get_Ip_error(p1, scale_j0, 1) + p1 = 0.8 + + err1, rel_err1, res_jphi = get_j0_error(p1, 1) if rel_err1 < tolerance: return p1, res_jphi @@ -526,7 +480,7 @@ def get_Ip_error(scale_val, scale_j0, n): for i in range(max_iter): if verbose: print(f"--- Optimization Iteration {i+1} ---") - + # Avoid division by zero if abs(err1 - err0) < 1e-9: if verbose: @@ -536,25 +490,22 @@ def get_Ip_error(scale_val, scale_j0, n): # Secant formula: x_new = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) # This projects the line between (p0, err0) and (p1, err1) to zero. p_new = p1 - err1 * (p1 - p0) / (err1 - err0) - + # Safety clamp: If the projection is wild (negative or huge), constrain it # This acts as a loose bracket to prevent the solver from crashing - if p_new < 0.1: + if p_new < 0.1: p_new = 0.1 - if p_new > 5.0: + if p_new > 5.0: p_new = 5.0 - + # Execute Solver at new point - if find_j0: - err_new, rel_err_new, res_jphi = get_j0_error(p_new, i+2) - else: - err_new, rel_err_new, res_jphi = get_Ip_error(p_new, scale_j0, i+2) + err_new, rel_err_new, res_jphi = get_j0_error(p_new, i+2) if rel_err_new < tolerance: if verbose: print(f"Converged! Optimal scale factor: {p_new:.4f}") return p_new, res_jphi - + # Update points for next step (move window forward) p0, err0 = p1, err1 p1, err1 = p_new, err_new @@ -563,6 +514,7 @@ def get_Ip_error(scale_val, scale_j0, n): print("Max iterations reached. Returning best last effort.") return p1, res_jphi + def calculate_ln_lambda(Te, Ti, ne, ni, Zeff=1.0, electron_lnLambda_model='sauter', ion_lnLambda_model='unity'): @@ -1056,18 +1008,14 @@ def calculate_profiles_and_bootstrap(psi_N, include_jBS): # Find optimal jphi scale final_scale_j0, final_jphi = find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, - Ip_target, psi_pad, spike_prof=spike_prof, find_j0=True, + Ip_target, psi_pad, spike_prof=spike_prof, diagnostic_plots=diagnostic_plots, verbose=verbose ) - if verbose: - print('\n >>> Finding optimal Ip scale factor') - # Find optimal Ip_target scale - final_scale_Ip, _ = find_optimal_scale(mygs, - psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, - Ip_target, psi_pad, spike_prof=spike_prof, find_j0=False, - scale_j0=final_scale_j0, - tolerance=0.001, diagnostic_plots=diagnostic_plots, verbose=verbose - ) + # Ip-scale secant removed: the solver now holds Ip to target + # natively (jphi-linterp cut-cell fix + Ip outer loop in the gs + # solve), so rescaling Ip_target here is redundant. The core-j0 + # secant above is retained. + final_scale_Ip = 1.0 if verbose: print('\n >>> Iterating on H-mode equilibrium solution') From 21fcd91191e6717abf8d99d2285151ba093a66d3 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 4 Jun 2026 17:15:31 -0400 Subject: [PATCH 11/13] fix(TokaMaker): address Copilot review on #1 - Fortran (identical to #290 update c68d383): signal failure on the jphi-linterp Ip safety-bail via error_str so a corrupted equilibrium no longer returns as a nominal success; fix the bail debug label. - bootstrap.py `analyze_bootstrap_edge_spike`: guard an empty edge slice (psi_N may not reach psi_N>=0.7) so numpy.max no longer raises; compute the edge max once. --- .../OpenFUSIONToolkit/TokaMaker/bootstrap.py | 10 ++++++++- src/python/wrappers/tokamaker_f.F90 | 22 ++++++++++++++----- 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py index 0459b21a4..f3eef9060 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py @@ -207,9 +207,17 @@ def analyze_bootstrap_edge_spike(psi_N, j_bootstrap, diagnostic_plots=False): psi_edge = psi_N[edge_mask] j_edge = j_bootstrap[edge_mask] + # Guard against an empty edge slice (psi_N may not reach the edge region): + # numpy.max on a zero-size array raises, where the original code fell through + # to find_peaks returning no peaks. + if j_edge.size == 0: + print("No edge region (psi_N >= 0.7) in profile") + return None + # Find peak in the edge region (prominence filter suppresses noise # oscillations that can be misidentified for small bootstrap spikes) - min_prominence = 0.05 * numpy.max(j_edge) if numpy.max(j_edge) > 0 else 0.0 + j_edge_max = numpy.max(j_edge) + min_prominence = 0.05 * j_edge_max if j_edge_max > 0 else 0.0 peaks, properties = find_peaks(j_edge, height=0., prominence=min_prominence) if len(peaks) == 0: diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index cdf94cbfd..294afd8b9 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -602,7 +602,7 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka REAL(r8) :: itor_target_orig, ip_phys, ip_phys_target, correction, rel_err REAL(r8) :: cent_dummy(2), v_dummy, pv_dummy, df_dummy, tf_dummy, bv_dummy INTEGER(i4) :: nl_its_total -LOGICAL :: ip_phys_bad +LOGICAL :: ip_phys_bad, ip_phys_bail IF(.NOT.tokamaker_ccast(tMaker_ptr,tMaker_obj,error_str))RETURN IF(.NOT.tokamaker_require_equil(tMaker_obj,error_str))RETURN IF(ANY(tMaker_obj%device%rcoils>0.d0).AND.(tMaker_obj%device%dt>0.d0))THEN @@ -669,6 +669,7 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka nl_its_total = 0 ierr = 0 rel_err = 0.d0 +ip_phys_bail = .FALSE. IF(oft_debug_print(1))WRITE(*,'(A,ES16.8)') & ' [JPHI_IP] outer-loop entry target_internal=', ip_phys_target DO outer_it = 1, jphi_ip_max_outer @@ -722,11 +723,12 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka ip_phys_bad = .TRUE. END IF IF(ip_phys_bad)THEN + ip_phys_bail = .TRUE. IF(oft_debug_print(1))WRITE(*,'(A,I0,A,2ES16.8)') & ' [JPHI_IP] iter ', outer_it, & - ' SAFETY BAIL -- ip_phys out of sane range, ip_phys/target =', & + ' SAFETY BAIL -- ip_phys out of sane range, ip_phys, target =', & ABS(ip_phys), ip_phys_target - EXIT ! caller will see equilibrium as solver produced it; restore handled below + EXIT ! signalled as a hard error below so the caller's rollback can recover END IF rel_err = (ABS(ip_phys) - ip_phys_target) / ip_phys_target IF(oft_debug_print(1))WRITE(*,'(A,I0,A,3ES16.8)') & @@ -754,14 +756,24 @@ SUBROUTINE tokamaker_solve(tMaker_ptr,vacuum,nl_its,error_str) BIND(C,NAME="toka tMaker_obj%gs_equil%Itor_target = SIGN(itor_target_orig, & tMaker_obj%gs_equil%Itor_target) -IF(ierr == 0 .AND. ABS(rel_err) >= jphi_ip_tol)THEN +IF(ierr == 0 .AND. .NOT.ip_phys_bail .AND. ABS(rel_err) >= jphi_ip_tol)THEN ! Did not converge within max_outer; not an error, but worth noting IF(oft_debug_print(1))WRITE(*,*) & ' tokamaker_solve: jphi-linterp Ip-corr did not converge, final rel_err =', rel_err END IF IF(vacuum)tMaker_obj%gs_equil%has_plasma=vac_save -IF(ierr/=0)CALL copy_string(gs_err_reason(ierr),error_str) +IF(ip_phys_bail)THEN + ! The inner solve returned ierr==0 but gs_comp_globals reported a corrupted / + ! non-physical Ip (<=0, NaN, or far outside the sane window). Report this as + ! a hard failure so the Python wrapper raises and the caller's rollback / + ! restore-and-resolve path is triggered -- otherwise a garbage equilibrium + ! would be returned as a nominal success. (ierr is local to this routine; + ! error_str is what the wrapper inspects.) + CALL copy_string('jphi-linterp Ip-correction bailed: physical Ip out of valid range (corrupted equilibrium)',error_str) +ELSE IF(ierr/=0)THEN + CALL copy_string(gs_err_reason(ierr),error_str) +END IF nl_its = nl_its_total END SUBROUTINE tokamaker_solve !--------------------------------------------------------------------------------- From 1362bb8c58d3fdba31e7bd997ea3490b74fb7f96 Mon Sep 17 00:00:00 2001 From: Chris Hansen Date: Thu, 4 Jun 2026 20:19:14 -0400 Subject: [PATCH 12/13] Improve behavior of Ip matching --- src/physics/grad_shaf.F90 | 2 +- src/physics/grad_shaf_prof_phys.F90 | 12 +++++------- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/physics/grad_shaf.F90 b/src/physics/grad_shaf.F90 index 6637b2f9f..a1a4eefdb 100644 --- a/src/physics/grad_shaf.F90 +++ b/src/physics/grad_shaf.F90 @@ -2235,7 +2235,7 @@ subroutine gs_solve(self,equil,ierr) ELSE IF(equil%estore_target>0.d0)THEN param_rhs(2)=equil%estore_target param_mat(2,2)=estored*3.d0/2.d0 - ELSE IF(equil%pax_target>0.d0)THEN + ELSE IF((equil%pax_target>0.d0).AND.(.NOT.equil%Ip_target_skip))THEN param_mat(2,2)=equil%P%f(equil%plasma_bounds(2)) param_rhs(2)=equil%pax_target ELSE IF(equil%Ip_ratio_target>-1.d98)THEN diff --git a/src/physics/grad_shaf_prof_phys.F90 b/src/physics/grad_shaf_prof_phys.F90 index c998b108c..402516c2a 100644 --- a/src/physics/grad_shaf_prof_phys.F90 +++ b/src/physics/grad_shaf_prof_phys.F90 @@ -523,7 +523,7 @@ subroutine jphi_update(self,gseq) class(jphi_flux_func), intent(inout) :: self class(gs_equil), intent(inout) :: gseq INTEGER(i4) :: i -REAL(r8) :: jphi_norm,pscale,pprime +REAL(r8) :: jphi_norm,pscale,pprime,dnorm REAL(r8), ALLOCATABLE :: ravgs(:,:),qtmp(:),psi_q(:) type(spline_type) :: R_spline self%plasma_bounds=gseq%plasma_bounds @@ -552,11 +552,12 @@ subroutine jphi_update(self,gseq) CALL spline_fit(R_spline,"extrap") DEALLOCATE(ravgs,psi_q,qtmp) !---Update jphi normalization to match Ip target +CALL gseq%P%update(gseq) ! Make sure pressure profile is up to date with EQ IF(gseq%Ip_target_skip)THEN CALL gs_itor_nl(gseq,jphi_norm) - jphi_norm=(ABS(gseq%Ip_target)*mu0/jphi_norm + self%norm_last)/2.d0 + dnorm=gseq%Ip_target/jphi_norm + jphi_norm=(1.d0+dnorm)*self%norm_last/2.d0 self%norm_last=jphi_norm - ! WRITE(*,*)'Ip NL: ',jphi_norm ELSE ALLOCATE(qtmp(self%npsi)) DO i=1,self%npsi @@ -567,12 +568,8 @@ subroutine jphi_update(self,gseq) DEALLOCATE(qtmp) jphi_norm=ABS(gseq%Ip_target)/jphi_norm self%norm_last=jphi_norm - ! WRITE(*,*)'Ip flux: ',jphi_norm END IF -! jphi_norm=ABS(gseq%Ip_target)/jphi_norm -! WRITE(*,*)'Ip flux: ',jphi_norm !---Get pressure profile -CALL gseq%P%update(gseq) ! Make sure pressure profile is up to date with EQ IF(ASSOCIATED(gseq%P_ani))CALL oft_abort('Jphi profiles do not support anistopic pressure','jphi_update',__FILE__) !CALL gseq%P_ani%update(gseq) pscale=gseq%P%f(gseq%plasma_bounds(2)) pscale=gseq%pax_target/pscale @@ -589,6 +586,7 @@ subroutine jphi_update(self,gseq) gseq%Ip_target_skip=.TRUE. ! IF(gseq%Ip_target>0.d0)gseq%Ip_target=-gseq%Ip_target gseq%ffp_scale=1.d0 +gseq%p_scale=pscale !---Clean up CALL spline_dealloc(R_spline) i=self%set_cofs(self%yp) From cce2b95ff596b7b20549a2d970e59cfd10e29092 Mon Sep 17 00:00:00 2001 From: Chris Hansen Date: Thu, 4 Jun 2026 20:24:53 -0400 Subject: [PATCH 13/13] Rename target skipping flag --- src/physics/grad_shaf.F90 | 78 ++++++++++++++++------------- src/physics/grad_shaf_prof_phys.F90 | 4 +- 2 files changed, 46 insertions(+), 36 deletions(-) diff --git a/src/physics/grad_shaf.F90 b/src/physics/grad_shaf.F90 index a1a4eefdb..ec45badd7 100644 --- a/src/physics/grad_shaf.F90 +++ b/src/physics/grad_shaf.F90 @@ -263,7 +263,7 @@ MODULE oft_gs TYPE :: gs_equil LOGICAL :: diverted = .FALSE. !< Equilibrium is diverted? LOGICAL :: has_plasma = .TRUE. !< Solve with plasma? (otherwise vacuum) - LOGICAL :: Ip_target_skip = .FALSE. !< Skip toroidal current target in non-linear solve? + LOGICAL :: skip_targets = .FALSE. !< Skip toroidal current target in non-linear solve? INTEGER(i4) :: mode = 0 !< RHS source mode (0 -> F*F', 1 -> F') INTEGER(i4) :: nx_points = 0 !< Number of X-points in current solution INTEGER(i4) :: nregularize = 0 !< Number of regularization terms @@ -1047,7 +1047,7 @@ subroutine copy_eq(self,source) self%mirror_n=source%mirror_n self%mirror_bturn=source%mirror_bturn self%mirror_zthroat=source%mirror_zthroat -self%Ip_target_skip=source%Ip_target_skip +self%skip_targets=source%skip_targets self%Ip_target=source%Ip_target self%estore_target=source%estore_target self%pax_target=source%pax_target @@ -2021,7 +2021,7 @@ subroutine gs_solve(self,equil,ierr) !--- error_flag=0 self%nl_its=0 -equil%Ip_target_skip=.FALSE. +equil%skip_targets=.FALSE. IF(TRIM(self%lu_solver%package)=='none')THEN CALL oft_abort("LU solver required for GS solve","gs_solve",__FILE__) ELSE @@ -2180,12 +2180,17 @@ subroutine gs_solve(self,equil,ierr) param_mat(1,1)=1.d0 param_rhs(1)=0.d0 ELSE - IF((equil%Ip_target>0.d0).AND.(.NOT.equil%Ip_target_skip))THEN - param_mat(1,:)=[itor_ffp,itor_press,0.d0] - param_rhs(1)=equil%Ip_target - ELSE + IF(equil%skip_targets)THEN param_mat(1,1)=1.d0 param_rhs(1)=equil%ffp_scale + ELSE + IF(equil%Ip_target>0.d0)THEN + param_mat(1,:)=[itor_ffp,itor_press,0.d0] + param_rhs(1)=equil%Ip_target + ELSE + param_mat(1,1)=1.d0 + param_rhs(1)=equil%ffp_scale + END IF END IF END IF @@ -2215,35 +2220,40 @@ subroutine gs_solve(self,equil,ierr) param_rhs(2)=equil%p_scale END IF ELSE - IF(equil%R0_target>0.d0)THEN - ! - psi_geval%u=>psi_vac - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi0) - param_rhs(2)=-gpsi0(1) - ! - psi_geval%u=>psi_vcont - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi0) - psi_geval%u=>psi_ffp - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi1) - psi_geval%u=>psi_press - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi2) - param_mat(2,:)=[gpsi1(1),gpsi2(1),gpsi0(1)] - ELSE IF(equil%estore_target>0.d0)THEN - param_rhs(2)=equil%estore_target - param_mat(2,2)=estored*3.d0/2.d0 - ELSE IF((equil%pax_target>0.d0).AND.(.NOT.equil%Ip_target_skip))THEN - param_mat(2,2)=equil%P%f(equil%plasma_bounds(2)) - param_rhs(2)=equil%pax_target - ELSE IF(equil%Ip_ratio_target>-1.d98)THEN - param_rhs(2)=0.d0 - param_mat(2,:)=[itor_ffp,-itor_press*equil%Ip_ratio_target,0.d0] - ELSE + IF(equil%skip_targets)THEN param_mat(2,2)=1.d0 param_rhs(2)=equil%p_scale + ELSE + IF(equil%R0_target>0.d0)THEN + ! + psi_geval%u=>psi_vac + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi0) + param_rhs(2)=-gpsi0(1) + ! + psi_geval%u=>psi_vcont + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi0) + psi_geval%u=>psi_ffp + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi1) + psi_geval%u=>psi_press + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi2) + param_mat(2,:)=[gpsi1(1),gpsi2(1),gpsi0(1)] + ELSE IF(equil%estore_target>0.d0)THEN + param_rhs(2)=equil%estore_target + param_mat(2,2)=estored*3.d0/2.d0 + ELSE IF(equil%pax_target>0.d0)THEN + param_mat(2,2)=equil%P%f(equil%plasma_bounds(2)) + param_rhs(2)=equil%pax_target + ELSE IF(equil%Ip_ratio_target>-1.d98)THEN + param_rhs(2)=0.d0 + param_mat(2,:)=[itor_ffp,-itor_press*equil%Ip_ratio_target,0.d0] + ELSE + param_mat(2,2)=1.d0 + param_rhs(2)=equil%p_scale + END IF END IF END IF diff --git a/src/physics/grad_shaf_prof_phys.F90 b/src/physics/grad_shaf_prof_phys.F90 index 402516c2a..f786a14e5 100644 --- a/src/physics/grad_shaf_prof_phys.F90 +++ b/src/physics/grad_shaf_prof_phys.F90 @@ -553,7 +553,7 @@ subroutine jphi_update(self,gseq) DEALLOCATE(ravgs,psi_q,qtmp) !---Update jphi normalization to match Ip target CALL gseq%P%update(gseq) ! Make sure pressure profile is up to date with EQ -IF(gseq%Ip_target_skip)THEN +IF(gseq%skip_targets)THEN CALL gs_itor_nl(gseq,jphi_norm) dnorm=gseq%Ip_target/jphi_norm jphi_norm=(1.d0+dnorm)*self%norm_last/2.d0 @@ -583,7 +583,7 @@ subroutine jphi_update(self,gseq) self%yp(i) = 2.d0*(self%jphi(i)*jphi_norm - R_spline%f(1)*pprime*pscale)/R_spline%f(2) END DO ! Disable Ip matching and fix F*F' scale (matching is done here instead) -gseq%Ip_target_skip=.TRUE. +gseq%skip_targets=.TRUE. ! IF(gseq%Ip_target>0.d0)gseq%Ip_target=-gseq%Ip_target gseq%ffp_scale=1.d0 gseq%p_scale=pscale