Skip to content

Pressure lost due to a saturated island #239

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
logan-nc wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

Pressure lost due to a saturated island #239

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b683a20
GPEC - NEW - Adds estimate of the pressure lost due to a saturated is…
logan-nc Feb 6, 2026
e5e7660
GPEC - WIP - Attempting to get a fully modified pressure profile output
logan-nc Feb 6, 2026
c0ecd0a
Merge branch 'develop' into island_pressure_modifications
matt-pharr Jun 8, 2026
8c987d0
GPEC - MINOR - Fix merge declarations
matt-pharr Jun 8, 2026
5af4681
GPEC - FIX - Correct saturated-island modified pressure profile output
priyanshlunia Jul 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 112 additions & 16 deletions gpec/gpout.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1549,9 +1549,9 @@ 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,
$ rh_id, r1_id,wc_id,wmin_id,wsat_id,pm_id,ps_id,
$ astat

INTEGER :: itheta,ising,icoup
Expand Down Expand Up @@ -1585,11 +1585,17 @@ 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,slayer_shear
$ 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
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(mstep) :: p_mod

c-----------------------------------------------------------------------
c solve equation from the given poloidal perturbation.
c-----------------------------------------------------------------------
Expand Down Expand Up @@ -1705,6 +1711,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)
Expand Down Expand Up @@ -1830,6 +1838,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.
Expand Down Expand Up @@ -1858,16 +1875,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",
Expand All @@ -1882,6 +1900,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 the psi_n profile grid (psifac, 1:mstep)
CALL spline_alloc(spl,npsi_tmp,1)
spl%xs = psi_tmp
spl%fs(:, 1) = p_tmp
CALL spline_fit(spl,"extrap")
DO ipsi=1,mstep
CALL spline_eval(spl,psifac(ipsi),0)
p_mod(ipsi) = spl%f(1)
ENDDO
CALL spline_dealloc(spl)
ELSE
! No island correction - use equilibrium pressure
DO ipsi=1,mstep
CALL spline_eval(sq,psifac(ipsi),0)
p_mod(ipsi) = sq%f(2) / mu0
ENDDO
ENDIF

c-----------------------------------------------------------------------
c write results.
c-----------------------------------------------------------------------
IF(ascii_flag)THEN
Expand All @@ -1896,16 +1972,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)),
Expand All @@ -1914,7 +1990,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
Expand Down Expand Up @@ -1966,6 +2043,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") )
Expand Down Expand Up @@ -2027,6 +2114,12 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot,
CALL check( nf90_put_att(fncid, a_id, "long_name",
$ "Surface area of rational surface") )
ENDIF
CALL check( nf90_inq_dimid(fncid, "psi_n", ps_id) )
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") )
CALL check( nf90_enddef(fncid) )
singflx = (/(singflx_mn(resnum(ising),ising), ising=1,msing)/)
CALL check( nf90_put_var(fncid, p_id,
Expand All @@ -2041,6 +2134,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) )
Expand All @@ -2055,6 +2150,7 @@ SUBROUTINE gpout_singfld(egnum,xspmn,spot,nspot,
IF(astat/=nf90_noerr)THEN
CALL check( nf90_put_var(fncid, a_id, area) )
ENDIF
CALL check( nf90_put_var(fncid, pm_id, p_mod) )
CALL check( nf90_close(fncid) )
ENDIF

Expand Down Expand Up @@ -2143,9 +2239,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((/
Expand All @@ -2162,7 +2258,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

Expand Down