From b683a20dec4e11f55afd27be6bad80be91e9e7cf Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 6 Feb 2026 17:01:43 -0500 Subject: [PATCH 1/2] GPEC - NEW - Adds estimate of the pressure lost due to a saturated island at each rational surface --- gpec/gpout.f | 56 ++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/gpec/gpout.f b/gpec/gpout.f index 68dcd25c..784faf7f 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -1581,7 +1581,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, COMPLEX(r8), DIMENSION(mpert), INTENT(IN) :: xspmn INTEGER :: i_id,q_id,m_id,p_id,c_id,bp_id,w_id,k_id,n_id,d_id, - $ a_id,pp_id,cp_id,wp_id,np_id,dp_id, + $ a_id,pp_id,cp_id,wp_id,np_id,dp_id,dd_id,pr_id, $ bc_id,ti_id, te_id, ni_id, ne_id, we_id, wi_id, q1_id, $ rh_id, r1_id,wc_id,wmin_id,wsat_id, $ astat @@ -1617,8 +1617,9 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,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 - REAL(r8) :: omega_i,omega_e,jxb,omega_sol,br_th + $ 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, + $ pleft, pright COMPLEX(r8) :: delta_s,psi0 c----------------------------------------------------------------------- @@ -1736,6 +1737,8 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, hw_sat(ising) = 0.0_r8 hw_min(ising) = 0.0_r8 b_crit(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 resm = mfac(resnum(ising)) CALL spline_eval(kin,respsi,1) @@ -1861,6 +1864,15 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, ENDIF ! convert from meters to psi_n for clear comparision to hw_isl hw_v_crit(ising) = hw_v_crit(ising) / sr%f1(1) + + ! If we are above hw_v_crit then we can estimate the pressure lost assuming flat + ! pressure inside the island + CALL spline_eval(sq,MAX(0.0, respsi - hw_sat(ising)),1) + pleft = sq%f(2) + CALL spline_eval(sq,MIN(1.0, respsi + hw_sat(ising)),1) + pright = sq%f(2) + dP_r(ising) = (pleft - pright) / mu0 ! Pascal + CALL spline_eval(sq,respsi,1) ! assumed to be at resonant surface latetr in this subroutine ENDIF c----------------------------------------------------------------------- c compute threshold by linear drift mhd with slayer module. @@ -1876,16 +1888,17 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, IF (verbose) THEN IF (callen_threshold_flag .OR. slayer_threshold_flag) THEN - IF(ising == 1) WRITE(*,'(1x,10a13)') + IF(ising == 1) WRITE(*,'(1x,11a13)') $ "psi","q","singflx","singlfx_crit","chirikov", $ "w_island", "w_v", "w_v_crit", - $ "w_sat","w_min" - WRITE(*,'(1x,es13.3,f13.3,2es13.3,f13.3,5es13.3)') + $ "w_sat","w_min","dP/P" + WRITE(*,'(1x,es13.3,f13.3,2es13.3,f13.3,6es13.3)') $ respsi,sq%f(4),ABS(singflx_mn(resnum(ising),ising)), $ ABS(b_crit(ising)), $ chirikov(ising),2*hw_isl(ising), $ 2*hw_v(ising),2*hw_v_crit(ising), - $ 2*hw_sat(ising),2*hw_min(ising) + $ 2*hw_sat(ising),2*hw_min(ising), + $ dP_r(ising)/P_r(ising) ELSE IF(ising == 1) WRITE(*,'(1x,a12,a12,a12,a12,a12)') "psi", @@ -1914,16 +1927,16 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, WRITE(out_unit,'(1x,a12,es17.8e3)')"sweet-spot =",spot WRITE(out_unit,'(1x,a12,1x,I4)')"msing =",msing WRITE(out_unit,*) - WRITE(out_unit,'(1x,a6,16(1x,a16))')"q","psi", + WRITE(out_unit,'(1x,a6,18(1x,a16))')"q","psi", $ "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_sat","half_w_min" + $ "half_w_sat","half_w_min","dP","P" DO ising=1,msing - WRITE(out_unit,'(1x,f6.3,15(es17.8e3))') + WRITE(out_unit,'(1x,f6.3,16(es17.8e3))') $ singtype(ising)%q,singtype(ising)%psifac, $ REAL(singflx_mn(resnum(ising),ising)), $ AIMAG(singflx_mn(resnum(ising),ising)), @@ -1932,7 +1945,8 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, $ REAL(delta(ising)),AIMAG(delta(ising)), $ hw_isl(ising),chirikov(ising), $ hw_v_crit(ising),b_crit(ising), - $ hw_sat(ising),hw_min(ising) + $ hw_sat(ising),hw_min(ising), + $ dP_r(ising),P_r(ising) ENDDO WRITE(out_unit,*) ENDIF @@ -1984,6 +1998,16 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, CALL check( nf90_put_att(fncid, wsat_id, "units", "psi_n") ) CALL check( nf90_put_att(fncid, wsat_id, "long_name", $ "Saturated island width from Callen model") ) + CALL check( nf90_def_var(fncid, "dP_res", nf90_double, + $ (/q_id/), dp_id) ) + CALL check( nf90_put_att(fncid, dp_id, "units", "Pa") ) + CALL check( nf90_put_att(fncid, dp_id, "long_name", + $ "Pressure lost due to saturated island") ) + CALL check( nf90_def_var(fncid, "P_res", nf90_double, + $ (/q_id/), pr_id) ) + CALL check( nf90_put_att(fncid, pr_id, "units", "Pa") ) + CALL check( nf90_put_att(fncid, pr_id, "long_name", + $ "Pressure at rational surface") ) CALL check( nf90_def_var(fncid, "Phi_res_crit", nf90_double, $ (/q_id/), bc_id) ) CALL check( nf90_put_att(fncid, bc_id, "units", "T") ) @@ -2059,6 +2083,8 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, CALL check( nf90_put_var(fncid, wc_id, 2*hw_v_crit) ) CALL check( nf90_put_var(fncid, wmin_id, 2*hw_min) ) 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, k_id, chirikov) ) CALL check( nf90_put_var(fncid, ti_id, ti_r) ) @@ -2172,9 +2198,9 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, CALL check( nf90_put_att(mncid, d_id, "long_name", $ "External Delta prime overlap") ) CALL check( nf90_def_var(mncid, "Delta_overlap_norm", - $ nf90_double,(/m_id/), dp_id) ) - CALL check( nf90_put_att(mncid, dp_id, "units", "untiless")) - CALL check( nf90_put_att(mncid, dp_id, "long_name", + $ nf90_double,(/m_id/), dd_id) ) + CALL check( nf90_put_att(mncid, dd_id, "units", "untiless")) + CALL check( nf90_put_att(mncid, dd_id, "long_name", $ "External Delta prime overlap percentage") ) CALL check( nf90_enddef(mncid) ) CALL check( nf90_put_var(mncid, p_id, RESHAPE((/ @@ -2191,7 +2217,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, CALL check( nf90_put_var(mncid, cp_id, op(2,:)) ) CALL check( nf90_put_var(mncid, wp_id, op(3,:)) ) CALL check( nf90_put_var(mncid, np_id, op(4,:)) ) - CALL check( nf90_put_var(mncid, dp_id, op(5,:)) ) + CALL check( nf90_put_var(mncid, dd_id, op(5,:)) ) CALL check( nf90_close(mncid) ) ENDIF From e5e7660d416d17d4c0764d8e8f6e579074d1eac3 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 6 Feb 2026 18:05:55 -0500 Subject: [PATCH 2/2] GPEC - WIP - Attempting to get a fully modified pressure profile output --- docs/examples/DIIID_ideal_example/coil.in | 8 +-- gpec/gpout.f | 80 ++++++++++++++++++++++- 2 files changed, 81 insertions(+), 7 deletions(-) diff --git a/docs/examples/DIIID_ideal_example/coil.in b/docs/examples/DIIID_ideal_example/coil.in index 946e00dc..fe0ff9c4 100644 --- a/docs/examples/DIIID_ideal_example/coil.in +++ b/docs/examples/DIIID_ideal_example/coil.in @@ -30,17 +30,17 @@ coil_cur(1,5)=-656 coil_cur(1,6)= 326 coil_name(2)="il" ! External correction coil set - coil_cur(2,1)= 1e3 ! 1kA, 30 degree phase n=1 field + coil_cur(2,1)= 1e5 ! 1kA, 30 degree phase n=1 field coil_cur(2,2)= 5e2 coil_cur(2,3)=-5e2 - coil_cur(2,4)=-1e3 + coil_cur(2,4)=-1e5 coil_cur(2,5)=-5e2 coil_cur(2,6)= 5e2 coil_name(3)="iu" ! External correction coil set - coil_cur(3,1)= 1e3 ! 1kA, 30 degree phase n=1 field + coil_cur(3,1)= 1e6 ! 1kA, 30 degree phase n=1 field coil_cur(3,2)= 5e2 coil_cur(3,3)=-5e2 - coil_cur(3,4)=-1e3 + coil_cur(3,4)=-1e6 coil_cur(3,5)=-5e2 coil_cur(3,6)= 5e2 diff --git a/gpec/gpout.f b/gpec/gpout.f index 784faf7f..d28b6eed 100644 --- a/gpec/gpout.f +++ b/gpec/gpout.f @@ -1583,7 +1583,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, INTEGER :: i_id,q_id,m_id,p_id,c_id,bp_id,w_id,k_id,n_id,d_id, $ a_id,pp_id,cp_id,wp_id,np_id,dp_id,dd_id,pr_id, $ bc_id,ti_id, te_id, ni_id, ne_id, we_id, wi_id, q1_id, - $ rh_id, r1_id,wc_id,wmin_id,wsat_id, + $ rh_id, r1_id,wc_id,wmin_id,wsat_id,pm_id,ps_id, $ astat INTEGER :: itheta,ising,icoup @@ -1618,10 +1618,15 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, 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 - REAL(r8) :: omega_i,omega_e,jxb,omega_sol,br_th, - $ pleft, pright + REAL(r8) :: omega_i,omega_e,jxb,omega_sol,br_th, + $ pleft, pright, psileft, psiright COMPLEX(r8) :: delta_s,psi0 + INTEGER, PARAMETER :: npsi_tmp = 1000 + INTEGER :: ipsi + REAL(r8), DIMENSION(0:npsi_tmp) :: psi_tmp, dp_tmp, p_tmp + REAL(r8), DIMENSION(0:mpsi) :: p_mod + c----------------------------------------------------------------------- c solve equation from the given poloidal perturbation. c----------------------------------------------------------------------- @@ -1913,6 +1918,64 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, CALL cspline_dealloc(fsp_sol) CALL gpeq_dealloc c----------------------------------------------------------------------- +c Calculate modified pressure profile from saturated islands +c----------------------------------------------------------------------- + IF (callen_threshold_flag) THEN + ! Create fine temporary grid for pressure calculation + DO ipsi=0,npsi_tmp + psi_tmp(ipsi) = REAL(ipsi,r8)/REAL(npsi_tmp,r8) + ENDDO + + ! Get pressure derivative at each point + DO ipsi=0,npsi_tmp + CALL spline_eval(sq,psi_tmp(ipsi),1) + dp_tmp(ipsi) = sq%f1(2) / mu0 + ENDDO + + ! Zero out pressure derivative in island regions + DO ising=1,msing + IF (hw_sat(ising) > 0.0_r8) THEN + psileft = MAX(0.0_r8, + $ singtype(ising)%psifac - hw_sat(ising)) + psiright = MIN(1.0_r8, + $ singtype(ising)%psifac + hw_sat(ising)) + DO ipsi=0,npsi_tmp + IF (psi_tmp(ipsi) >= psileft .AND. + $ psi_tmp(ipsi) <= psiright) THEN + dp_tmp(ipsi) = 0.0_r8 + ENDIF + ENDDO + ENDIF + ENDDO + + ! Integrate using trapezoidal rule to get modified pressure + CALL spline_eval(sq,1.0_r8,0) + p_tmp(npsi_tmp) = sq%f(2) / mu0 + DO ipsi=npsi_tmp-1,0,-1 + p_tmp(ipsi) = p_tmp(ipsi+1) + 0.5_r8 * + $ (dp_tmp(ipsi+1) + dp_tmp(ipsi)) / + $ REAL(npsi_tmp,r8) + ENDDO + + ! Interpolate to sq%xs grid (0:mpsi) using spline + CALL spline_alloc(spl,npsi_tmp,1) + spl%xs = psi_tmp + spl%fs(:, 1) = p_tmp + CALL spline_fit(spl,"extrap") + DO ipsi=0,mpsi + CALL spline_eval(spl,sq%xs(ipsi),0) + p_mod(ipsi) = spl%f(1) + ENDDO + CALL spline_dealloc(spl) + ELSE + ! No island correction - use equilibrium pressure + DO ipsi=0,mpsi + CALL spline_eval(sq,sq%xs(ipsi),0) + p_mod(ipsi) = sq%f(2) / mu0 + ENDDO + ENDIF + +c----------------------------------------------------------------------- c write results. c----------------------------------------------------------------------- IF(ascii_flag)THEN @@ -2068,6 +2131,14 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, CALL check( nf90_put_att(fncid, a_id, "units", "m^2") ) CALL check( nf90_put_att(fncid, a_id, "long_name", $ "Surface area of rational surface") ) + ENDIF + astat = nf90_inq_dimid(fncid, "psi_n", ps_id) + IF(astat==nf90_noerr)THEN + CALL check( nf90_def_var(fncid, "p_mod", nf90_double, + $ (/ps_id/), pm_id) ) + CALL check( nf90_put_att(fncid, pm_id, "units", "Pa") ) + CALL check( nf90_put_att(fncid, pm_id, "long_name", + $ "Modified pressure with saturated island correction") ) ENDIF CALL check( nf90_enddef(fncid) ) singflx = (/(singflx_mn(resnum(ising),ising), ising=1,msing)/) @@ -2099,6 +2170,9 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot, IF(astat/=nf90_noerr)THEN CALL check( nf90_put_var(fncid, a_id, area) ) ENDIF + IF(astat==nf90_noerr)THEN + CALL check( nf90_put_var(fncid, pm_id, p_mod) ) + ENDIF CALL check( nf90_close(fncid) ) ENDIF