Skip to content

Commit

Permalink
cgyro precision control fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@jcandy
jcandy committed Dec 20, 2024
1 parent 6c2029a commit a289d8d
Show file tree
Hide file tree
Showing 6 changed files with 106 additions and 63 deletions.
39 changes: 27 additions & 12 deletions cgyro/bin/cgyrodmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,9 @@ def opts():

mytext = '''\
output:
NEED DOCUMENTATION HERE
This tool will analyze a linear CGYRO run based on the time
history of various field. Subdoninant and even unstable modes
can be computed.
'''

parser = argparse.ArgumentParser(
Expand Down Expand Up @@ -109,7 +111,9 @@ fdict = {'phi' :r'$\phi$',
'ni':'$n_i$',
'ne':'$n_e$',
'vi':'$v_i$',
've':'$v_e$'}
've':'$v_e$',
'ei':'$E_i$',
'ee':'$E_e$'}

# symbols
args = {}
Expand All @@ -122,6 +126,10 @@ args[ostr[2]] = {'color':'k','marker':'+'}
sim = cgyrodata(mydir,silent=True)
sim.getbigfield()

if sim.hiprec_flag == 0:
print('ERROR: (cgyrodmd) Need to run cgyro with HIPREC_FLAG=1')
sys.exit()

t = sim.t
n_radial = sim.n_radial
n_theta = sim.theta_plot
Expand All @@ -148,35 +156,42 @@ y = sim.kxky_v[0,:,:,0,0,:imax]+1j*sim.kxky_v[1,:,:,0,0,:imax]
ovec['vi'],a0 = map1d(y,sim.q)
y = sim.kxky_v[0,:,:,1,0,:imax]+1j*sim.kxky_v[1,:,:,1,0,:imax]
ovec['ve'],a0 = map1d(y,sim.q)
y = sim.kxky_e[0,:,:,0,0,:imax]+1j*sim.kxky_e[1,:,:,0,0,:imax]
ovec['ei'],a0 = map1d(y,sim.q)
y = sim.kxky_e[0,:,:,1,0,:imax]+1j*sim.kxky_e[1,:,:,1,0,:imax]
ovec['ee'],a0 = map1d(y,sim.q)

# step for DMD
dt = k*(t[1]-t[0])

#---------------------------------------------------------------------------
# RUN DMD
# Initialize (classic) DMD object
dmd = DMD(svd_rank=0,exact=True,sorted_eigs='abs')

evec = {}
modes = {}
# dictionary of eigenvalues
edict = {}
# dictionary of eigenmodes
mdict = {}
for x in ostr:
down = downsample(ovec[x],k)
dmd.fit(down)
evec[x] = 1j*np.log(dmd.eigs)/dt
modes[x] = dmd.modes
edict[x] = 1j*np.log(dmd.eigs)/dt
mdict[x] = dmd.modes

# determine most unstable modes (zvec)

zvec = np.zeros([3,nmode],dtype=complex)
# vector of sorted eigenvalues
evec = np.zeros([3,nmode],dtype=complex)
# vector of sorted eigenmodes
mvec = np.zeros([3,n_radial*n_theta,nmode],dtype=complex)
for i,x in enumerate(ostr):
z = evec[x]
z = edict[x]
zi = z.imag
k = np.flip(np.argsort(zi))
zvec[i,:] = z[k[:nmode]]
mvec[i,:,:] = modes[x][:,k[:nmode]]
evec[i,:] = z[k[:nmode]]
mvec[i,:,:] = mdict[x][:,k[:nmode]]

# find true eigenmodes and compute errors

m01 = []
m02 = []
for i in range(nmode):
Expand Down
17 changes: 10 additions & 7 deletions cgyro/src/cgyro_globals.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,6 @@ module cgyro_globals

use, intrinsic :: iso_fortran_env

! Data output precision setting
integer, parameter :: BYTE=8 ! Change to 8 for double precision

!---------------------------------------------------------------
! Input parameters:
Expand Down Expand Up @@ -271,11 +269,16 @@ module cgyro_globals
character(len=2) :: mpiio_small_stripe_str
character(len=3) :: mpiio_stripe_str
!
! Standard precision for IO (there are optionally reset to higher precision later)
character(len=8) :: fmtstr ='(es11.4)'
integer :: fmtstr_len = 12
character(len=15) :: fmtstrn ='(10(es11.4,1x))'
character(len=9) :: fmtstr_hi ='(es18.12)'
! Precision for IO (these are set in cgyro_init_manager)
!
! BYTE=4 (standard single precision)
! BYTE=8 (double precision for DMD analysis)
!
integer :: BYTE
character(len=8) :: fmtstr
integer :: fmtstr_len
character(len=15) :: fmtstrn
character(len=9) :: fmtstr_prec ='(es18.12)'
!----------------------------------------------------

!---------------------------------------------------------------
Expand Down
10 changes: 8 additions & 2 deletions cgyro/src/cgyro_init_manager.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,16 @@ subroutine cgyro_init_manager
integer :: ie,ix

if (hiprec_flag == 1) then
fmtstr = '(es16.9)'
BYTE = 8
fmtstr = '(es16.9)'
fmtstr_len = 17
fmtstrn = '(10(es16.9,1x))'
endif
else
BYTE = 4
fmtstr ='(es11.4)'
fmtstr_len = 12
fmtstrn ='(10(es11.4,1x))'
endif

!------------------------------------------------------
! Initialize startup timers
Expand Down
44 changes: 31 additions & 13 deletions cgyro/src/cgyro_write_initdata.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,6 +231,7 @@ subroutine cgyro_write_initdata
enddo
! Added 17 Dec 2024
write (io,fmtstr) z_eff
write (io,'(i0)') hiprec_flag
close(io)

endif
Expand All @@ -242,19 +243,35 @@ subroutine cgyro_write_initdata
!
if (silent_flag == 0 .and. i_proc == 0) then
open(unit=io,file=trim(path)//'bin.cgyro.geo',status='replace',access='stream')
write(io) real(theta,kind=4)
write(io) real(g_theta_geo,kind=4)
write(io) real(bmag,kind=4)
write(io) real(omega_stream(:,1,nt1),kind=4)
write(io) real(omega_trap(:,1,nt1),kind=4)
write(io) real(omega_rdrift(:,1),kind=4)
write(io) real(omega_adrift(:,1),kind=4)
write(io) real(omega_aprdrift(:,1),kind=4)
write(io) real(omega_cdrift(:,1),kind=4)
write(io) real(omega_cdrift_r(:,1),kind=4)
write(io) real(omega_gammap(:),kind=4)
write(io) real(k_perp(ic_c(n_radial/2+1,:),nt1),kind=4)
write(io) real(captheta,kind=4)
if (hiprec_flag == 0) then
write(io) real(theta,kind=4)
write(io) real(g_theta_geo,kind=4)
write(io) real(bmag,kind=4)
write(io) real(omega_stream(:,1,nt1),kind=4)
write(io) real(omega_trap(:,1,nt1),kind=4)
write(io) real(omega_rdrift(:,1),kind=4)
write(io) real(omega_adrift(:,1),kind=4)
write(io) real(omega_aprdrift(:,1),kind=4)
write(io) real(omega_cdrift(:,1),kind=4)
write(io) real(omega_cdrift_r(:,1),kind=4)
write(io) real(omega_gammap(:),kind=4)
write(io) real(k_perp(ic_c(n_radial/2+1,:),nt1),kind=4)
write(io) real(captheta,kind=4)
else
write(io) real(theta)
write(io) real(g_theta_geo)
write(io) real(bmag)
write(io) real(omega_stream(:,1,nt1))
write(io) real(omega_trap(:,1,nt1))
write(io) real(omega_rdrift(:,1))
write(io) real(omega_adrift(:,1))
write(io) real(omega_aprdrift(:,1))
write(io) real(omega_cdrift(:,1))
write(io) real(omega_cdrift_r(:,1))
write(io) real(omega_gammap(:))
write(io) real(k_perp(ic_c(n_radial/2+1,:),nt1))
write(io) real(captheta)
endif
close(io)

endif
Expand Down Expand Up @@ -308,6 +325,7 @@ subroutine cgyro_write_initdata
endif
write(io,'(1pe13.6)') (spectraldiss((pi/n_toroidal)*in,nup_alpha),in=0,n_toroidal-1)
write(io,'(1pe13.6)') (spectraldiss((2*pi/n_radial)*p,nup_radial),p=-n_radial/2,n_radial/2-1)
write(io,'(i0)') hiprec_flag
close(io)

endif
Expand Down
11 changes: 8 additions & 3 deletions cgyro/src/cgyro_write_timedata.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -499,7 +499,7 @@ subroutine write_precision(datafile,fn)
! Append

open(unit=io,file=datafile,status='old',position='append')
write(io,fmtstr_hi) fn_sum
write(io,fmtstr_prec) fn_sum
close(io)

case(3)
Expand All @@ -508,7 +508,7 @@ subroutine write_precision(datafile,fn)

open(unit=io,file=datafile,status='old')
do i=1,i_current
read(io,fmtstr_hi) fn_sum
read(io,fmtstr_prec) fn_sum
enddo
endfile(io)
close(io)
Expand Down Expand Up @@ -856,7 +856,12 @@ subroutine write_binary(datafile,fn,n_fn)

open(unit=io,file=datafile,status='old',position='append',access='stream')

write(io) cmplx(fn,kind=BYTE)
if (BYTE == 4) then
write(io) cmplx(fn,kind=4)
else
write(io) cmplx(fn,kind=8)
endif

close(io)

case (3)
Expand Down
48 changes: 22 additions & 26 deletions f2py/pygacode/cgyro/data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,6 @@
from gacodefuncs import *
except:
from ..gacodefuncs import *

BYTE='float'

# class to read cgyro output data
class cgyrodata:
Expand All @@ -17,8 +15,25 @@ def __init__(self,sim_directory,silent=False,fast=False):

self.silent = silent
self.dir = sim_directory

# read time vector
hastime = self.gettime()

# get mesh/input data
self.getgrid()

# set data resolution
if self.hiprec_flag:
# double
self.BYTE = 'float64'
else:
# single
self.BYTE = 'float32'

if not self.silent:
print('INFO: (data.py) Detected precision '+self.BYTE)


if hastime and not fast:
self.getdata()

Expand All @@ -28,7 +43,7 @@ def extract(self,f):
start = time.time()
if os.path.isfile(self.dir+'bin'+f):
fmt = 'bin'
data = np.fromfile(self.dir+'bin'+f,dtype=BYTE)
data = np.fromfile(self.dir+'bin'+f,dtype=self.BYTE)
elif os.path.isfile(self.dir+'out'+f):
fmt = 'out'
data = np.fromfile(self.dir+'out'+f,dtype='float',sep=' ')
Expand All @@ -47,34 +62,15 @@ def gettime(self):
return False

#-----------------------------------------------------------------
# Read time vector -- autodetect number of columns (ncol)
#

# START: 3-column output (delete this code eventually)
tfile = self.dir+'out.cgyro.time'
with open(tfile,'r') as f:
ncol = len(f.readline().split())

if ncol == 3:
newlines = []
with open(tfile,'r') as f:
for line in f.readlines():
newlines.append(line.strip()+' 0.0')
with open(tfile,'w') as outfile:
outfile.write("\n".join(newlines))
print('INFO: (data.py) Updated to 4-column format for out.cgyro.time')

# Read time vector
#
data = np.fromfile(self.dir+'out.cgyro.time',dtype='float',sep=' ')
nt = len(data)//4
data = np.reshape(data,(4,nt),'F')
# END: 3-column output (delete this code eventually)

self.t = data[0,:]
self.err1 = data[1,:]
try:
self.err2 = data[2,:]
except:
self.err2 = data[1,:]
self.err2 = data[2,:]
self.n_time = nt
if not self.silent:
print('INFO: (data.py) Read time vector in out.cgyro.time.')
Expand Down Expand Up @@ -512,6 +508,7 @@ def getgrid(self):
self.sbeta[i],p = self.eget(data,p)
# Added 17 Dec 2024
self.z_eff,p = self.eget(data,p)
self.hiprec_flag,p = self.eget(data,p)

if p == -1:
print('WARNING: (getgrid) Data format outdated. Please run cgyro -t')
Expand All @@ -536,7 +533,6 @@ def getnorm(self,norm):
except:
pass


self.rhonorm = self.rho
self.rhoi = r'\rho_s'

Expand Down

0 comments on commit a289d8d

Please sign in to comment.