Skip to content

Commit

Permalink
Cleanups of xflux code
Browse files Browse the repository at this point in the history
  • Loading branch information
Jeff Candy committed Feb 20, 2024
1 parent 09c423a commit 75b9d87
Show file tree
Hide file tree
Showing 4 changed files with 52 additions and 64 deletions.
8 changes: 5 additions & 3 deletions cgyro/src/cgyro_write_initdata.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,17 +126,19 @@ subroutine cgyro_write_initdata
is,int(z(is)),dens(is),temp(is),mass(is),dlnndr(is),dlntdr(is),nu(is)
enddo

! Profile shear
if (profile_shear_flag == 1) then
write(io,*)
write(io,'(a)') ' i s(a/Ln) (a/Ln)_L (a/Ln)_R | s(a/Lt) (a/Lt)_L (a/Lt)_R '
do is=1,n_species
dn = sdlnndr(is)*length/rho/2
dt = sdlntdr(is)*length/rho/2
dn = sdlnndr(is)*length/rho/4
dt = sdlntdr(is)*length/rho/4
write(io,'(t1,i2,3(1x,1pe9.3),2x,3(1x,1pe9.3))') &
is,sdlnndr(is),dlnndr(is)-dn,dlnndr(is)+dn,sdlntdr(is),dlntdr(is)-dt,dlntdr(is)+dt
enddo
endif


! Running from input.gacode
if (profile_model == 2) then
dn = rho/(rhos/a_meters)
kyrat = abs(q/rmin*rhos/a_meters)
Expand Down
45 changes: 20 additions & 25 deletions f2py/pygacode/cgyro/data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,10 +201,7 @@ def getxflux(self):

def xfluxave(self,w,moment,e=0.2,nscale=0):

"""
Do complicated spatial averages for xflux
RESULT: self.lky_flux_ave
"""
# Averaging functionfor global fluxes

print('INFO: (xfluxave) Computing partial-domain averages')

Expand All @@ -221,41 +218,39 @@ def xfluxave(self,w,moment,e=0.2,nscale=0):
else:
sc[:] = 1.0

# Sum moments over fields (3) and toroidal modes (4)
# NOTE: lky_flux_n[2,ng,ns,nfield,n_n,nt]
if moment == 'n':
z = np.sum(self.lky_flux_n,axis=(3,4))
elif moment == 'e':
z = np.sum(self.lky_flux_e,axis=(3,4))
elif moment == 'v':
z = np.sum(self.lky_flux_v,axis=(3,4))
else:
raise ValueError('(xfluxave) Invalid moment.')


#--------------------------------------------
# Useful arrays required outside this routine
self.lky_xr = np.zeros((ns,ng))
self.lky_xi = np.zeros((ns,ng))
self.lky_flux_ave = np.zeros((ns,2))
# Arrays required outside this routine
self.lky_xr = np.zeros([ng,ns])
self.lky_xi = np.zeros([ng,ns])
self.lky_flux_ave = np.zeros([ns,2])
#--------------------------------------------

imin,imax=time_index(self.t,w)

for ispec in range(ns):
for l in range(ng):
self.lky_xr[ispec,l] = time_average(z[0,l,ispec,:],self.t,imin,imax)*sc[ispec]
self.lky_xi[ispec,l] = time_average(z[1,l,ispec,:],self.t,imin,imax)*sc[ispec]
imin,imax = time_index(self.t,w)

# Time averages of real and imaginary parts
self.lky_xr[:,:] = time_average(z[0,:,:,:],self.t,imin,imax)
self.lky_xi[:,:] = time_average(z[1,:,:,:],self.t,imin,imax)

l = np.arange(1,ng)
u = (2*np.pi*e)*l
for ispec in range(ns):
# Flux partial average over [-e,e]
g0 = self.lky_xr[ispec,0]
g1 = g0
for l in range(1,ng):
u = 2*np.pi*l*e
g0 = g0+2*np.sin(u)*self.lky_xr[ispec,l]/u
g1 = g1+2*np.sin(u)*self.lky_xr[ispec,l]/u*(-1)**l
g0 = self.lky_xr[0,ispec]+2*np.sum(np.sin(u)*self.lky_xr[l,ispec]/u)
g1 = self.lky_xr[0,ispec]+2*np.sum(np.sin(u)*self.lky_xr[l,ispec]/u*(-1)**l)

# Average over true (positive) interval
self.lky_flux_ave[ispec,0] = g0
self.lky_flux_ave[ispec,0] = g0*sc[ispec]
# Average over negative interval
self.lky_flux_ave[ispec,1] = g1
self.lky_flux_ave[ispec,1] = g1*sc[ispec]

def getbigfield(self):

Expand Down
61 changes: 26 additions & 35 deletions f2py/pygacode/cgyro/data_plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1006,40 +1006,29 @@ def plot_xflux(self,xin):
nscale = xin['nscale']
ymin = xin['ymin']
ymax = xin['ymax']
mirror = xin['abs'])

if xin['fig'] is None:
fig = plt.figure(MYDIR,figsize=(xin['lx'],xin['ly']))

self.getxflux()

ns = self.n_species
nl = self.n_global+1
ng = self.n_global+1
t = self.t

ky = self.ky
ave = np.zeros((self.n_n,ns))

if moment == 'phi':
moment = 'e'

# NOTE: lky_flux_* -> [ 2, nl , ns , n_n , nt ]
# 0 1 2 3 4


if moment == 'n':
ntag = 'Density~flux'
mtag = '\Gamma'
z = np.sum(self.lky_flux_n,axis=3)
ftag = 'xflux_n'
elif moment == 'e':
ntag = 'Energy~flux'
mtag = 'Q'
z = np.sum(self.lky_flux_e,axis=3)
ftag = 'xflux_e'
elif moment == 'v':
ntag = 'Momentum~flux'
mtag = '\Pi'
z = np.sum(self.lky_flux_v,axis=3)
ftag = 'xflux_v'
else:
print('ERROR: (plot_xflux) Invalid moment.')
sys.exit()
Expand All @@ -1055,7 +1044,6 @@ def plot_xflux(self,xin):
else:
mnorm = ''


#============================================================
# Otherwise plot
ax = fig.add_subplot(111)
Expand All @@ -1065,51 +1053,54 @@ def plot_xflux(self,xin):

color = ['k','m','b','c','g','r']

imin,imax=time_index(t,w)
imin,imax = time_index(t,w)
mpre,mwin = wintxt(imin,imax,t)

ax.set_title(r'$\mathrm{'+ntag+'} \quad $'+mwin)

a = -np.pi+2*np.pi*np.arange(0.0,1.0,0.001)

na = 128

a = np.linspace(-np.pi,np.pi,na)
ah = np.linspace(0,-2*np.pi,na)

for ispec in range(ns):

u = specmap(self.mass[ispec],self.z[ispec])

# Flux curve
g = np.zeros(len(t))
g = self.lky_xr[ispec,0]
for l in range(1,nl):
g = g+2*(np.cos(l*a)*self.lky_xr[ispec,l]-np.sin(l*a)*self.lky_xi[ispec,l])
#---------------------------------
# Global flux versus x
g = np.zeros(na) ; g[:] = self.lky_xr[0,ispec]
for l in range(1,ng):
g[:] = g[:]+2*(np.cos(l*a)*self.lky_xr[l,ispec]-np.sin(l*a)*self.lky_xi[l,ispec])
ax.plot(a/(2*np.pi),g,color=color[ispec])

g = np.zeros(na) ; g[:] = self.lky_xr[0,ispec]
for l in range(1,ng):
g[:] = g[:]+2*(np.cos(l*ah)*self.lky_xr[l,ispec]-np.sin(l*ah)*self.lky_xi[l,ispec])
ax.plot(a/(2*np.pi),g,color=color[ispec],linestyle='--')
#---------------------------------
# Flux partial average over [-e,e]

#---------------------------------
# Flux partial average over "positive" [-e,e] interval
g0 = self.lky_flux_ave[ispec,0]
label = r'$'+mtag+mnorm+'_'+u+'/'+mtag+'_\mathrm{GB}: '+str(round(g0,3))+'$'
ax.plot([-e,e],[g0,g0],'o-',color=color[ispec],alpha=0.2,linewidth=3,label=label)
#---------------------------------

#---------------------------------
# Flux partial average over "negative" interval
# Flux partial average over "negative" [-e,e] interval
g1 = self.lky_flux_ave[ispec,1]
ax.plot([0.5-e,0.5],[g1,g1],'o--',color=color[ispec],alpha=0.2,linewidth=3)
ax.plot([-0.5,-0.5+e],[g1,g1],'o--',color=color[ispec],alpha=0.2,linewidth=3)
#---------------------------------

#---------------------------------
# Flux spectral average
gs = self.lky_xr[ispec,0]+2*np.pi/4*self.lky_xr[ispec,1]
#---------------------------------


#---------------------------------
# Flux domain average
ga = self.lky_xr[ispec,0]
ga = self.lky_xr[0,ispec]
ax.plot([-0.5,0.5],[ga,ga],color=color[ispec],alpha=0.5)
#---------------------------------

print('INFO: (plot_xflux) Ave [inner/inner_spec, outer, domain] = '
'{:.2f}/{:.2f}, {:.2f}, {:.2f}'.format(g0,gs,g1,ga))
print('INFO: (plot_xflux) Ave [inner-e, outer-e, domain] = '
'{:.2f}, {:.2f}, {:.2f}'.format(g0,g1,ga))

if ymax != 'auto':
ax.set_ylim(top=float(ymax))
Expand Down
2 changes: 1 addition & 1 deletion platform/exec/exec.MINT
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,5 +9,5 @@ numa=${5}
mpinuma=${6}

cd $simdir
mpiexec -env OMP_NUM_THREADS $nomp -n $nmpi $exec
mpiexec -env OMP_NUM_THREADS $nomp -n $nmpi $exec

0 comments on commit 75b9d87

Please sign in to comment.