Skip to content

Commit 49ece65

Browse files
authored
TokaMaker: Add <1/R^2> to fields computed by get_q() (OpenFUSIONToolkit#313)
- Convert ravgs return value for `get_q()` and `sauter_fc()` from array to dict for clarity
1 parent 0cb54b0 commit 49ece65

12 files changed

Lines changed: 227 additions & 214 deletions

File tree

src/examples/TokaMaker/ARC/ARCV3A_baseline_ex.ipynb

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -571,7 +571,7 @@
571571
],
572572
"source": [
573573
"psi,f,fp,p,pp = mygs.get_profiles()\n",
574-
"psi_q,qvals,ravgs,dl,rbounds,zbounds = mygs.get_q(psi_pad=0.01)\n",
574+
"psi_q,qvals, _,_,_,_ = mygs.get_q(psi_pad=0.01)\n",
575575
"fig, ax = plt.subplots(3,1,sharex=True)\n",
576576
"psi_eqdsk = np.linspace(0.0,1.0,eqdsk['nr'])\n",
577577
"psi_sample = np.linspace(0.0,1.0,10)\n",

src/examples/TokaMaker/DIIID/DIIID_baseline_ex.ipynb

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -622,7 +622,7 @@
622622
],
623623
"source": [
624624
"psi,f,fp,p,pp = mygs.get_profiles()\n",
625-
"psi_q,qvals,ravgs,dl,rbounds,zbounds = mygs.get_q(psi_pad=0.01)\n",
625+
"psi_q,qvals, _,_,_,_ = mygs.get_q(psi_pad=0.01)\n",
626626
"fig, ax = plt.subplots(3,1,sharex=True)\n",
627627
"psi_eqdsk = np.linspace(0.0,1.0,eqdsk['nr'])\n",
628628
"psi_sample = np.linspace(0.0,1.0,10)\n",
@@ -1063,7 +1063,7 @@
10631063
],
10641064
"source": [
10651065
"psi,f,fp,p,pp = mygs.get_profiles()\n",
1066-
"psi_q,qvals,ravgs,dl,rbounds,zbounds = mygs.get_q(psi_pad=0.01)\n",
1066+
"psi_q,qvals, _,_,_,_ = mygs.get_q(psi_pad=0.01)\n",
10671067
"fig, ax = plt.subplots(3,1,sharex=True)\n",
10681068
"psi_eqdsk = np.linspace(0.0,1.0,eqdsk['nr'])\n",
10691069
"psi_sample = np.linspace(0.0,1.0,10)\n",

src/examples/TokaMaker/ITER/ITER_Hmode_ex.ipynb

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -634,7 +634,7 @@
634634
"outputs": [],
635635
"source": [
636636
"psi,f,fp,p,pp = mygs.get_profiles(npsi=257,psi_pad=1e-3)\n",
637-
"_, q, ravgs, _, _, _ = mygs.get_q(npsi=257,psi_pad=1e-3) # get flux averaged R and 1/R from equilibrium solution\n",
637+
"_, q, _,_,_,_ = mygs.get_q(npsi=257,psi_pad=1e-3) # get flux averaged R and 1/R from equilibrium solution\n",
638638
"\n",
639639
"# Retrieve desired j_phi profile components\n",
640640
"final_jphi = results['total_j_phi']\n",
@@ -729,7 +729,7 @@
729729
],
730730
"source": [
731731
"fig, axs = plt.subplots(2,1,figsize=(6,8),sharex=True)\n",
732-
"_,q,_,_,_,_ = mygs.get_q(npsi=n_sample)\n",
732+
"_, q, _,_,_,_ = mygs.get_q(npsi=n_sample)\n",
733733
"\n",
734734
"# Plot j_total\n",
735735
"axs[0].plot(psi,final_jphi/1e+6,label=r'Output $j_{\\phi}$') # in MA/m^2\n",
@@ -973,7 +973,7 @@
973973
"source": [
974974
"fig, axs = plt.subplots(2,1,figsize=(6,8),sharex=True)\n",
975975
"\n",
976-
"_,q,_,_,_,_ = mygs.get_q(npsi=n_sample)\n",
976+
"_, q, _,_,_,_ = mygs.get_q(npsi=n_sample)\n",
977977
"\n",
978978
"# Plot j_total\n",
979979
"axs[0].plot(psi,results['total_j_phi']/1e+6,label=r'Output $j_{\\phi}$') # in MA/m^2\n",

src/examples/TokaMaker/fixed_boundary/fixed_boundary_ex1.ipynb

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -719,7 +719,7 @@
719719
],
720720
"source": [
721721
"psi,f,fp,p,pp = mygs.get_profiles()\n",
722-
"psi_q,qvals,ravgs,dl,rbounds,zbounds = mygs.get_q()\n",
722+
"psi_q,qvals, _,_,_,_ = mygs.get_q()\n",
723723
"fig, ax = plt.subplots(3,1,sharex=True)\n",
724724
"psi_eqdsk = np.linspace(0.0,1.0,EQ_in['nr'])\n",
725725
"# Plot F*F'\n",
@@ -891,7 +891,7 @@
891891
],
892892
"source": [
893893
"psi,f,fp,p,pp = mygs.get_profiles()\n",
894-
"psi_q,qvals,ravgs,dl,rbounds,zbounds = mygs.get_q(psi_pad=0.005)\n",
894+
"psi_q,qvals,ravgs, _,_,_ = mygs.get_q(psi_pad=0.005)\n",
895895
"fig, ax = plt.subplots(4,1,sharex=True)\n",
896896
"psi_eqdsk = np.linspace(0.0,1.0,EQ_in['nr'])\n",
897897
"psi_sample = np.linspace(0.025,1.0,10)\n",
@@ -911,7 +911,7 @@
911911
"ax[2].plot(psi_eqdsk,EQ_in['qpsi'],'--')\n",
912912
"ax[2].set_ylabel(\"q\")\n",
913913
"# Plot < J_phi >\n",
914-
"ax[3].plot(psi_q,np.interp(psi_q,psi,pp)*ravgs[0,:]*mu0+np.interp(psi_q,psi,f*fp)*ravgs[1,:])\n",
914+
"ax[3].plot(psi_q,np.interp(psi_q,psi,pp)*ravgs['<R>']*mu0+np.interp(psi_q,psi,f*fp)*ravgs['<1/R>'])\n",
915915
"ax[3].set_ylabel(r\"< $J_{\\phi}$ >\")\n",
916916
"_ = ax[-1].set_xlabel(r\"$\\hat{\\psi}$\")"
917917
]

src/physics/grad_shaf.F90

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -4339,7 +4339,7 @@ subroutine gs_get_qprof(gseq,nr,psi_q,prof,dl,rbounds,zbounds,ravgs)
43394339
real(8), optional, intent(out) :: dl !< Arc length of surface `psi_q(1)` (should be LCFS)
43404340
real(8), optional, intent(out) :: rbounds(2,2) !< Radial bounds of surface `psi_q(1)` (should be LCFS)
43414341
real(8), optional, intent(out) :: zbounds(2,2) !< Vertical bounds of surface `psi_q(1)` (should be LCFS)
4342-
real(8), optional, intent(out) :: ravgs(nr,3) !< Flux surface averages <R>, <1/R>, and dV/dPsi
4342+
real(8), optional, intent(out) :: ravgs(nr,4) !< Flux surface averages <R>, <1/R>, <1/R^2>, and dV/dPsi
43434343
real(8) :: psi_surf,rmax,x1,x2,raxis,zaxis,fpol,qpsi
43444344
real(8) :: pt(3),pt_last(3),pt_proj(3),f(3),psi_tmp(1),gop(3,3)
43454345
type(oft_lag_brinterp), target :: psi_int
@@ -4402,7 +4402,7 @@ subroutine gs_get_qprof(gseq,nr,psi_q,prof,dl,rbounds,zbounds,ravgs)
44024402
CALL field%setup(gseq%device%fe_rep)
44034403
IF(PRESENT(ravgs))THEN
44044404
field%compute_geom=.TRUE.
4405-
active_tracer%neq=5
4405+
active_tracer%neq=6
44064406
ELSE
44074407
field%compute_geom=.FALSE.
44084408
active_tracer%neq=3
@@ -4486,7 +4486,8 @@ subroutine gs_get_qprof(gseq,nr,psi_q,prof,dl,rbounds,zbounds,ravgs)
44864486
IF(PRESENT(ravgs))THEN
44874487
ravgs(j,1)=active_tracer%v(4)/active_tracer%v(2)
44884488
ravgs(j,2)=active_tracer%v(5)/active_tracer%v(2)
4489-
ravgs(j,3)=-2.d0*pi*active_tracer%v(2) ! First derivative of FS volume (V')
4489+
ravgs(j,3)=active_tracer%v(6)/active_tracer%v(2)
4490+
ravgs(j,4)=-2.d0*pi*active_tracer%v(2) ! First derivative of FS volume (V')
44904491
END IF
44914492
end do
44924493
CALL active_tracer%delete
@@ -4849,6 +4850,7 @@ subroutine gsinv_apply(self,cell,f,gop,val)
48494850
IF(self%compute_geom)THEN
48504851
val(4)=pt(1)*val(2)
48514852
val(5)=val(2)/pt(1)
4853+
val(6)=val(2)/(pt(1)**2)
48524854
END IF
48534855
! val(3:8)=val(3:8)
48544856
deallocate(j)

src/physics/grad_shaf_prof_phys.F90

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -531,7 +531,7 @@ subroutine jphi_update(self,gseq)
531531
IF(gseq%Ip_target<0.d0)CALL oft_abort("Jphi profile requires Ip target","jphi_update",__FILE__)
532532
IF(gseq%pax_target<0.d0)CALL oft_abort("Jphi profile requires Pax target","jphi_update",__FILE__)
533533
!---Get updated flux surface geometry for Jphi -> F*F' mapping
534-
ALLOCATE(ravgs(self%ngeom+1,3),psi_q(self%ngeom+1),qtmp(self%ngeom+1))
534+
ALLOCATE(ravgs(self%ngeom+1,4),psi_q(self%ngeom+1),qtmp(self%ngeom+1))
535535
IF(gseq%diverted)THEN
536536
psi_q=[(REAL(i-1,8)/REAL(self%ngeom,8),i=1,self%ngeom)]
537537
psi_q(1)=psi_q(2)

src/physics/grad_shaf_util.F90

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ MODULE oft_gs_util
2727
USE tracing_2d, ONLY: active_tracer, tracinginv_fs, set_tracer
2828
USE mhd_utils, ONLY: mu0
2929
USE oft_gs, ONLY: gs_factory, flux_func, gs_dflux, gs_itor_nl, gs_test_bounds, gs_b_interp, &
30-
gs_get_qprof, gsinv_interp, gs_psi2r, gs_psi2pt, gs_epsilon, gs_update_bounds
30+
gsinv_interp, gs_psi2r, gs_psi2pt, gs_epsilon, gs_update_bounds
3131
USE oft_gs_profiles
3232
USE grad_shaf_prof_phys, ONLY: create_jphi_ff, jphi_flux_func
3333
IMPLICIT NONE
@@ -286,7 +286,7 @@ subroutine gs_calc_vloop(self,vloop)
286286
real(8) :: itor_loc ! local toroidal current in integration
287287
real(8) :: itor ! toroidal current
288288
real(8) :: I_NI ! non-inductive F*F'
289-
real(8) :: eta_jsq ! eta*j_NI**2
289+
real(8) :: eta_jsq ! eta*j_NI**2
290290
real(8) :: goptmp(3,3)
291291
real(8) :: v ! volume
292292
real(8) :: pt(3) ! radial coordinate

src/python/OpenFUSIONToolkit/TokaMaker/_core.py

Lines changed: 23 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1386,7 +1386,7 @@ def get_q(self,psi=None,psi_pad=0.02,npsi=50,compute_geo=False):
13861386
@param psi_pad End padding (axis and edge) for uniform sampling (ignored if `psi` is not None)
13871387
@param npsi Number of points for uniform sampling (ignored if `psi` is not None)
13881388
@param compute_geo Compute geometric values for LCFS
1389-
@result \f$\hat{\psi}\f$, \f$q(\hat{\psi})\f$, \f$[<R>,<1/R>,dV/dPsi]\f$, length of last surface,
1389+
@result \f$\hat{\psi}\f$, \f$q(\hat{\psi})\f$, \f${<R>,<1/R>,dV/dPsi}\f$, length of last surface,
13901390
[r(R_min),r(R_max)], [r(z_min),r(z_max)]
13911391
'''
13921392
if self._tMaker_equil is None:
@@ -1399,7 +1399,7 @@ def sauter_fc(self,psi=None,psi_pad=0.02,npsi=50):
13991399
@param psi Explicit sampling locations in \f$\hat{\psi}\f$
14001400
@param psi_pad End padding (axis and edge) for uniform sampling (ignored if `psi` is not None)
14011401
@param npsi Number of points for uniform sampling (ignored if `psi` is not None)
1402-
@result \f$ f_c \f$, [\f$<R>,<1/R>,<a>\f$], [\f$<|B|>,<|B|^2>\f$]
1402+
@result \f$ f_c \f$, \f${<R>,<1/R>,<a>}\f$, \f$[<|B|>,<|B|^2>]\f$
14031403
'''
14041404
if self._tMaker_equil is None:
14051405
raise ValueError("Equilibrium object is `None`")
@@ -2723,7 +2723,7 @@ def get_q(self,psi=None,psi_pad=0.02,npsi=50,compute_geo=False):
27232723
@param psi_pad End padding (axis and edge) for uniform sampling (ignored if `psi` is not None)
27242724
@param npsi Number of points for uniform sampling (ignored if `psi` is not None)
27252725
@param compute_geo Compute geometric values for LCFS
2726-
@result \f$\hat{\psi}\f$, \f$q(\hat{\psi})\f$, \f$[<R>,<1/R>,dV/dPsi]\f$, length of last surface,
2726+
@result \f$\hat{\psi}\f$, \f$q(\hat{\psi})\f$, \f${<R>,<1/R>,<1/R^2>,dV/dPsi}\f$, length of last surface,
27272727
[r(R_min),r(R_max)], [r(z_min),r(z_max)]
27282728
'''
27292729
if psi is None:
@@ -2736,7 +2736,7 @@ def get_q(self,psi=None,psi_pad=0.02,npsi=50,compute_geo=False):
27362736
psi_save = numpy.copy(psi)
27372737
psi = numpy.ascontiguousarray(1.0-psi, dtype=numpy.float64)
27382738
qvals = numpy.zeros((psi.shape[0],), dtype=numpy.float64)
2739-
ravgs = numpy.zeros((3,psi.shape[0]), dtype=numpy.float64)
2739+
ravgs = numpy.zeros((4,psi.shape[0]), dtype=numpy.float64)
27402740
if compute_geo:
27412741
dl = c_double(1.0)
27422742
else:
@@ -2747,16 +2747,22 @@ def get_q(self,psi=None,psi_pad=0.02,npsi=50,compute_geo=False):
27472747
tokamaker_get_q(self._equil_ptr,psi.shape[0],psi,qvals,ravgs,ctypes.byref(dl),rbounds,zbounds,error_string)
27482748
if error_string.value != b'':
27492749
raise Exception(error_string.value)
2750+
ravg_dict = {
2751+
'<R>': ravgs[0,:],
2752+
'<1/R>': ravgs[1,:],
2753+
'<1/R^2>': ravgs[2,:],
2754+
'dV/dPsi': ravgs[3,:]
2755+
}
27502756
if self.psi_convention == 0:
27512757
if compute_geo:
2752-
return psi_save,qvals,ravgs,dl.value,rbounds,zbounds
2758+
return psi_save,qvals,ravg_dict,dl.value,rbounds,zbounds
27532759
else:
2754-
return psi_save,qvals,ravgs,None,None,None
2760+
return psi_save,qvals,ravg_dict,None,None,None
27552761
else:
27562762
if compute_geo:
2757-
return psi,qvals,ravgs,dl.value,rbounds,zbounds
2763+
return psi,qvals,ravg_dict,dl.value,rbounds,zbounds
27582764
else:
2759-
return psi,qvals,ravgs,None,None,None
2765+
return psi,qvals,ravg_dict,None,None,None
27602766

27612767
def trace_surf(self,psi):
27622768
r'''! Trace surface for a given poloidal flux
@@ -2849,7 +2855,7 @@ def get_stats(self,lcfs_pad=None,axis_pad=0.02,li_normalization='std',geom_type=
28492855
lcfs_pad = 0.0
28502856
if self.diverted or (not self.free_boundary):
28512857
lcfs_pad = 0.01
2852-
_,qvals,_,dl,rbounds,zbounds = self.get_q(numpy.r_[1.0-lcfs_pad,0.95,axis_pad],compute_geo=True) # Given backward so last point is LCFS (for dl)
2858+
_, qvals, _, dl,rbounds,zbounds = self.get_q(numpy.r_[1.0-lcfs_pad,0.95,axis_pad],compute_geo=True) # Given backward so last point is LCFS (for dl)
28532859
# Get diverted topology information
28542860
if self.diverted:
28552861
x_points, _ = self.get_xpoints()
@@ -2978,7 +2984,7 @@ def calc_sauter_fc(self,psi=None,psi_pad=0.02,npsi=50):
29782984
@param psi Explicit sampling locations in \f$\hat{\psi}\f$
29792985
@param psi_pad End padding (axis and edge) for uniform sampling (ignored if `psi` is not None)
29802986
@param npsi Number of points for uniform sampling (ignored if `psi` is not None)
2981-
@result \f$ f_c \f$, [\f$<R>,<1/R>,<a>\f$], [\f$<|B|>,<|B|^2>\f$]
2987+
@result \f$ f_c \f$, \f${<R>,<1/R>,<a>}\f$, \f$[<|B|>,<|B|^2>]\f$
29822988
'''
29832989
if psi is None:
29842990
psi = numpy.linspace(psi_pad,1.0-psi_pad,npsi,dtype=numpy.float64)
@@ -2996,10 +3002,15 @@ def calc_sauter_fc(self,psi=None,psi_pad=0.02,npsi=50):
29963002
tokamaker_sauter_fc(self._equil_ptr,psi.shape[0],psi,fc,r_avgs,modb_avgs,error_string)
29973003
if error_string.value != b'':
29983004
raise Exception(error_string.value)
3005+
ravg_dict = {
3006+
'<R>': r_avgs[0,:],
3007+
'<1/R>': r_avgs[1,:],
3008+
'<a>': r_avgs[2,:]
3009+
}
29993010
if self.psi_convention == 0:
3000-
return psi_save,fc,r_avgs,modb_avgs
3011+
return psi_save,fc,ravg_dict,modb_avgs
30013012
else:
3002-
return psi,fc,r_avgs,modb_avgs
3013+
return psi,fc,ravg_dict,modb_avgs
30033014

30043015
def calc_delstar_curr(self,psi):
30053016
r'''! Get toroidal current density from \f$ \psi \f$ through \f$ \Delta^{*} \f$ operator

0 commit comments

Comments
 (0)