From 37a1b4c0dd8761526ff00279a633f8069f1b1f47 Mon Sep 17 00:00:00 2001 From: priyanshlunia <40486607+priyanshlunia@users.noreply.github.com> Date: Mon, 20 Jul 2026 11:07:58 -0400 Subject: [PATCH 1/3] GPEC - FEATURE - Add Phi_res_crit_callen critical resonant flux output Compute the resonant flux at which the vacuum island width equals the Callen critical width (w_vac = w_v_crit), as Phi_res*(w_crit/w_vac)^2 in gpout_singfld. This is the Callen analog of the SLAYER Phi_res_crit, enabling direct model comparison (issue #240). Output alongside the other Callen quantities: netCDF (Phi_res_crit_callen), the singfld ASCII table (Phi_crit_callen), and the verbose stdout table (Phi_crit_cal). Co-Authored-By: Claude Opus 4.8 --- gpec/gpout.f | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/gpec/gpout.f b/gpec/gpout.f index f2dedb06..2a666a75 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -1556,7 +1556,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, $ a_id,pp_id,cp_id,wp_id,np_id,dp_id,dd_id,pr_id,wc_id, $ bc_id,ti_id,te_id,ni_id,ne_id,we_id,wi_id,q1_id,rh_id, $ r1_id,ssp_id,lq_id,rt_id,at_id,wmin_id,wsat_id,pm_id, - $ ps_id,astat + $ ps_id,pcc_id,astat INTEGER :: itheta,ising,icoup REAL(r8) :: respsi,lpsi,rpsi,shear,hdist,sbnosurf @@ -1588,7 +1588,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, TYPE(spline_type) :: sr REAL(r8), DIMENSION(msing) :: b_crit, ti_r, te_r, ni_r, ne_r, - $ q1_r, we_r, wi_r, rh_r, r1_r, dP_r, P_r + $ q1_r, we_r, wi_r, rh_r, r1_r, dP_r, P_r, phi_crit_callen REAL(r8) :: omega_i,omega_e,jxb,omega_sol,br_th,slayer_shear, $ pleft, pright, psileft, psiright REAL(r8) :: slayer_inpr_loc @@ -1715,6 +1715,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, hw_sat(ising) = 0.0_r8 hw_min(ising) = 0.0_r8 b_crit(ising) = 0.0_r8 + phi_crit_callen(ising) = 0.0_r8 dP_r(ising) = 0.0_r8 P_r(ising) = sq%f(2) / mu0 IF (callen_threshold_flag. OR. slayer_threshold_flag) THEN @@ -1843,6 +1844,16 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, ! convert from meters to psi_n for clear comparision to hw_isl hw_v_crit(ising) = hw_v_crit(ising) / sr%f1(1) + ! Critical resonant flux for w_vac = w_v_crit (Callen); the + ! analog of the SLAYER Phi_res_crit (b_crit) for direct + ! model comparison. Since w_vac ~ sqrt(Phi_res), the flux + ! needed to reach the critical width is Phi_res*(w_crit/w_vac)^2 + ! (hw_v and hw_v_crit are both in psi_n here). + IF (ALLOCATED(vsingfld) .AND. hw_v(ising) > 0.0_r8) THEN + phi_crit_callen(ising) = ABS(vsingfld(ising)) + $ * (hw_v_crit(ising) / hw_v(ising))**2 + ENDIF + ! If we are above hw_v_crit then we can estimate the pressure lost assuming flat ! pressure inside the island. NOTE: this per-surface dP_r/P_r ! estimate uses the unmodified equilibrium pressure (sq) for @@ -1884,13 +1895,14 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, IF (verbose) THEN IF (callen_threshold_flag .OR. slayer_threshold_flag) THEN - IF(ising == 1) WRITE(*,'(1x,11a13)') - $ "psi","q","singflx","singlfx_crit","chirikov", + IF(ising == 1) WRITE(*,'(1x,12a13)') + $ "psi","q","singflx","singlfx_crit","Phi_crit_cal", + $ "chirikov", $ "w_island", "w_v", "w_v_crit", $ "w_sat","w_min","dP/P" - WRITE(*,'(1x,es13.3,f13.3,2es13.3,f13.3,6es13.3)') + WRITE(*,'(1x,es13.3,f13.3,3es13.3,f13.3,6es13.3)') $ respsi,sq%f(4),ABS(singflx_mn(resnum(ising),ising)), - $ ABS(b_crit(ising)), + $ ABS(b_crit(ising)),ABS(phi_crit_callen(ising)), $ chirikov(ising),2*hw_isl(ising), $ 2*hw_v(ising),2*hw_v_crit(ising), $ 2*hw_sat(ising),2*hw_min(ising), @@ -1967,16 +1979,16 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, WRITE(out_unit,*) WRITE(out_unit,'(1x,a12,1x,I4)')"msing =",msing WRITE(out_unit,*) - WRITE(out_unit,'(1x,a6,18(1x,a16))')"q","psi","spot", + WRITE(out_unit,'(1x,a6,19(1x,a16))')"q","psi","spot", $ "real(singflx)","imag(singflx)", $ "real(singcur)","imag(singcur)", $ "real(singbwp)","imag(singbwp)", $ "real(Delta)","imag(Delta)", $ "half_w_isl","chirikov", - $ "half_w_isl_v_crit","singflx_crit", + $ "half_w_isl_v_crit","singflx_crit","Phi_crit_callen", $ "half_w_sat","half_w_min","dP","P" DO ising=1,msing - WRITE(out_unit,'(1x,f6.3,18(es17.8e3))') + WRITE(out_unit,'(1x,f6.3,19(es17.8e3))') $ singtype(ising)%q,singtype(ising)%psifac, $ spots(ising), $ REAL(singflx_mn(resnum(ising),ising)), @@ -1985,7 +1997,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, $ REAL(singbwp(ising)),AIMAG(singbwp(ising)), $ REAL(delta(ising)),AIMAG(delta(ising)), $ hw_isl(ising),chirikov(ising), - $ hw_v_crit(ising),b_crit(ising), + $ hw_v_crit(ising),b_crit(ising),phi_crit_callen(ising), $ hw_sat(ising),hw_min(ising), $ dP_r(ising),P_r(ising) ENDDO @@ -2077,6 +2089,12 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, CALL check( nf90_put_att(fncid, bc_id, "units", "T") ) CALL check( nf90_put_att(fncid, bc_id, "long_name", $ "Critical resonant field for island growth from SLAYER") ) + CALL check( nf90_def_var(fncid, "Phi_res_crit_callen", + $ nf90_double, (/q_id/), pcc_id) ) + CALL check( nf90_put_att(fncid, pcc_id, "units", "T") ) + CALL check( nf90_put_att(fncid, pcc_id, "long_name", + $ "Critical resonant field for island growth from Callen "// + $ "model") ) CALL check( nf90_def_var(fncid, "K_isl", nf90_double, $ (/q_id/), k_id) ) CALL check( nf90_put_att(fncid, k_id, "long_name", @@ -2165,6 +2183,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, CALL check( nf90_put_var(fncid, dp_id, dP_r) ) CALL check( nf90_put_var(fncid, pr_id, P_r) ) CALL check( nf90_put_var(fncid, bc_id, b_crit) ) + CALL check( nf90_put_var(fncid, pcc_id, phi_crit_callen) ) CALL check( nf90_put_var(fncid, k_id, chirikov) ) CALL check( nf90_put_var(fncid, ti_id, ti_r) ) CALL check( nf90_put_var(fncid, te_id, te_r) ) From d2fff0760d4b5cd98818e0800d7b573359630bb7 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Thu, 23 Jul 2026 16:13:13 -0700 Subject: [PATCH 2/3] GPEC - FIX - Consistent slayer/callen naming for singflx_crit outputs Address PR #274 review: use the locked historical file naming convention for both threshold columns instead of an ad-hoc new style. - Rename internal vars b_crit -> b_crit_slayer, phi_crit_callen -> b_crit_callen throughout gpout_singfld. - stdout table: fix the "singlfx_crit" typo and rename to singflx_crit_slayer / singflx_crit_callen; widen those two columns 13 -> 20 (header and data) to fit the names and keep alignment. - ASCII gpec_singfld file: rename the same two columns and widen them a16 -> a19 / es17.8e3 -> es20.8e3. - netCDF variable names (Phi_res_crit, Phi_res_crit_callen) left unchanged to avoid breaking OMFIT. Pure naming/formatting change; values are unchanged. Verified on DIIID_ideal_example: clean build and consistent slayer/callen values across stdout, ASCII, and netCDF. Co-Authored-By: Claude Opus 4.8 --- gpec/gpout.f | 40 +++++++++++++++++++++++----------------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/gpec/gpout.f b/gpec/gpout.f index 2a666a75..3793ff7b 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -1587,8 +1587,9 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, REAL(r8), DIMENSION(0:mthsurf) :: r_tmp TYPE(spline_type) :: sr - REAL(r8), DIMENSION(msing) :: b_crit, ti_r, te_r, ni_r, ne_r, - $ q1_r, we_r, wi_r, rh_r, r1_r, dP_r, P_r, phi_crit_callen + REAL(r8), DIMENSION(msing) :: b_crit_slayer, b_crit_callen, + $ ti_r, te_r, ni_r, ne_r, + $ q1_r, we_r, wi_r, rh_r, r1_r, dP_r, P_r REAL(r8) :: omega_i,omega_e,jxb,omega_sol,br_th,slayer_shear, $ pleft, pright, psileft, psiright REAL(r8) :: slayer_inpr_loc @@ -1714,8 +1715,8 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, hw_v_crit(ising) = 0.0_r8 hw_sat(ising) = 0.0_r8 hw_min(ising) = 0.0_r8 - b_crit(ising) = 0.0_r8 - phi_crit_callen(ising) = 0.0_r8 + b_crit_slayer(ising) = 0.0_r8 + b_crit_callen(ising) = 0.0_r8 dP_r(ising) = 0.0_r8 P_r(ising) = sq%f(2) / mu0 IF (callen_threshold_flag. OR. slayer_threshold_flag) THEN @@ -1845,12 +1846,12 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, hw_v_crit(ising) = hw_v_crit(ising) / sr%f1(1) ! Critical resonant flux for w_vac = w_v_crit (Callen); the - ! analog of the SLAYER Phi_res_crit (b_crit) for direct + ! analog of the SLAYER b_crit_slayer for direct ! model comparison. Since w_vac ~ sqrt(Phi_res), the flux ! needed to reach the critical width is Phi_res*(w_crit/w_vac)^2 ! (hw_v and hw_v_crit are both in psi_n here). IF (ALLOCATED(vsingfld) .AND. hw_v(ising) > 0.0_r8) THEN - phi_crit_callen(ising) = ABS(vsingfld(ising)) + b_crit_callen(ising) = ABS(vsingfld(ising)) $ * (hw_v_crit(ising) / hw_v(ising))**2 ENDIF @@ -1889,20 +1890,21 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, $ kin%f(9),kin%f(5),omega_e,omega_i,sq%f(4),slayer_shear, $ bt0,sr%f(1),ro,mi,slayer_inpr_loc,resm,nn,ascii_flag, $ delta_s,psi0,jxb,omega_sol,br_th) - b_crit(ising)=br_th ! Tesla. Normal resonant field comparable to singflx + b_crit_slayer(ising)=br_th ! Tesla. Normal resonant field comparable to singflx ENDIF IF (verbose) THEN IF (callen_threshold_flag .OR. slayer_threshold_flag) THEN - IF(ising == 1) WRITE(*,'(1x,12a13)') - $ "psi","q","singflx","singlfx_crit","Phi_crit_cal", + IF(ising == 1) WRITE(*,'(1x,3a13,2a20,7a13)') + $ "psi","q","singflx", + $ "singflx_crit_slayer","singflx_crit_callen", $ "chirikov", $ "w_island", "w_v", "w_v_crit", $ "w_sat","w_min","dP/P" - WRITE(*,'(1x,es13.3,f13.3,3es13.3,f13.3,6es13.3)') + WRITE(*,'(1x,es13.3,f13.3,es13.3,2es20.3,f13.3,6es13.3)') $ respsi,sq%f(4),ABS(singflx_mn(resnum(ising),ising)), - $ ABS(b_crit(ising)),ABS(phi_crit_callen(ising)), + $ ABS(b_crit_slayer(ising)),ABS(b_crit_callen(ising)), $ chirikov(ising),2*hw_isl(ising), $ 2*hw_v(ising),2*hw_v_crit(ising), $ 2*hw_sat(ising),2*hw_min(ising), @@ -1979,16 +1981,19 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, WRITE(out_unit,*) WRITE(out_unit,'(1x,a12,1x,I4)')"msing =",msing WRITE(out_unit,*) - WRITE(out_unit,'(1x,a6,19(1x,a16))')"q","psi","spot", + WRITE(out_unit,'(1x,a6,13(1x,a16),2(1x,a19),4(1x,a16))') + $ "q","psi","spot", $ "real(singflx)","imag(singflx)", $ "real(singcur)","imag(singcur)", $ "real(singbwp)","imag(singbwp)", $ "real(Delta)","imag(Delta)", $ "half_w_isl","chirikov", - $ "half_w_isl_v_crit","singflx_crit","Phi_crit_callen", + $ "half_w_isl_v_crit", + $ "singflx_crit_slayer","singflx_crit_callen", $ "half_w_sat","half_w_min","dP","P" DO ising=1,msing - WRITE(out_unit,'(1x,f6.3,19(es17.8e3))') + WRITE(out_unit, + $ '(1x,f6.3,13(es17.8e3),2(es20.8e3),4(es17.8e3))') $ singtype(ising)%q,singtype(ising)%psifac, $ spots(ising), $ REAL(singflx_mn(resnum(ising),ising)), @@ -1997,7 +2002,8 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, $ REAL(singbwp(ising)),AIMAG(singbwp(ising)), $ REAL(delta(ising)),AIMAG(delta(ising)), $ hw_isl(ising),chirikov(ising), - $ hw_v_crit(ising),b_crit(ising),phi_crit_callen(ising), + $ hw_v_crit(ising), + $ b_crit_slayer(ising),b_crit_callen(ising), $ hw_sat(ising),hw_min(ising), $ dP_r(ising),P_r(ising) ENDDO @@ -2182,8 +2188,8 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, CALL check( nf90_put_var(fncid, wsat_id, 2*hw_sat) ) CALL check( nf90_put_var(fncid, dp_id, dP_r) ) CALL check( nf90_put_var(fncid, pr_id, P_r) ) - CALL check( nf90_put_var(fncid, bc_id, b_crit) ) - CALL check( nf90_put_var(fncid, pcc_id, phi_crit_callen) ) + CALL check( nf90_put_var(fncid, bc_id, b_crit_slayer) ) + CALL check( nf90_put_var(fncid, pcc_id, b_crit_callen) ) CALL check( nf90_put_var(fncid, k_id, chirikov) ) CALL check( nf90_put_var(fncid, ti_id, ti_r) ) CALL check( nf90_put_var(fncid, te_id, te_r) ) From eb0db17debdf9c453f136b253271278bfc332b8e Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Thu, 23 Jul 2026 20:09:55 -0700 Subject: [PATCH 3/3] GPEC - FIX - Widen/clarify singfld ASCII columns Follow-up to the slayer/callen naming fix, cleaning up the ASCII gpec_singfld output columns: - Widen half_w_isl_v_crit (17 chars) from a16 -> a19 so it is no longer truncated to "half_w_isl_v_cri". - Rename the ambiguous dP -> dP_w_sat (equilibrium pressure lost across the saturated island) and P -> P_res (pressure at the rational surface), matching the netCDF dP_res/P_res intent. Pure naming/formatting; values unchanged. Verified on DIIID_ideal_example (clean build, normal run, identical crit values). Co-Authored-By: Claude Opus 4.8 --- gpec/gpout.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/gpec/gpout.f b/gpec/gpout.f index 3793ff7b..fa5ccf32 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -1981,7 +1981,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, WRITE(out_unit,*) WRITE(out_unit,'(1x,a12,1x,I4)')"msing =",msing WRITE(out_unit,*) - WRITE(out_unit,'(1x,a6,13(1x,a16),2(1x,a19),4(1x,a16))') + WRITE(out_unit,'(1x,a6,12(1x,a16),3(1x,a19),4(1x,a16))') $ "q","psi","spot", $ "real(singflx)","imag(singflx)", $ "real(singcur)","imag(singcur)", @@ -1990,10 +1990,10 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spots,interpspot,nspot, $ "half_w_isl","chirikov", $ "half_w_isl_v_crit", $ "singflx_crit_slayer","singflx_crit_callen", - $ "half_w_sat","half_w_min","dP","P" + $ "half_w_sat","half_w_min","dP_w_sat","P_res" DO ising=1,msing WRITE(out_unit, - $ '(1x,f6.3,13(es17.8e3),2(es20.8e3),4(es17.8e3))') + $ '(1x,f6.3,12(es17.8e3),3(es20.8e3),4(es17.8e3))') $ singtype(ising)%q,singtype(ising)%psifac, $ spots(ising), $ REAL(singflx_mn(resnum(ising),ising)),