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 2 commits into
base: develop
Choose a base branch
from
Open

Pressure lost due to a saturated island #239

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 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
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
8 changes: 4 additions & 4 deletions docs/examples/DIIID_ideal_example/coil.in
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
132 changes: 116 additions & 16 deletions gpec/gpout.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1581,9 +1581,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 @@ -1617,10 +1617,16 @@ 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, 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-----------------------------------------------------------------------
Expand Down Expand Up @@ -1736,6 +1742,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 @@ -1861,6 +1869,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 All @@ -1876,16 +1893,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 @@ -1900,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
Expand All @@ -1914,16 +1990,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 @@ -1932,7 +2008,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 @@ -1984,6 +2061,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 @@ -2044,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)/)
Expand All @@ -2059,6 +2154,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 @@ -2073,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

Expand Down Expand Up @@ -2172,9 +2272,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 @@ -2191,7 +2291,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
Loading