From a289d8dfead06f38c8d371668ff848c171f2bd80 Mon Sep 17 00:00:00 2001 From: Jeff Candy Date: Thu, 19 Dec 2024 21:15:40 -0800 Subject: [PATCH] cgyro precision control fixes --- cgyro/bin/cgyrodmd | 39 ++++++++++++++++-------- cgyro/src/cgyro_globals.F90 | 17 ++++++----- cgyro/src/cgyro_init_manager.F90 | 10 +++++-- cgyro/src/cgyro_write_initdata.f90 | 44 +++++++++++++++++++-------- cgyro/src/cgyro_write_timedata.f90 | 11 +++++-- f2py/pygacode/cgyro/data.py | 48 ++++++++++++++---------------- 6 files changed, 106 insertions(+), 63 deletions(-) diff --git a/cgyro/bin/cgyrodmd b/cgyro/bin/cgyrodmd index 176ef3fdc..14709552c 100755 --- a/cgyro/bin/cgyrodmd +++ b/cgyro/bin/cgyrodmd @@ -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( @@ -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 = {} @@ -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 @@ -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): diff --git a/cgyro/src/cgyro_globals.F90 b/cgyro/src/cgyro_globals.F90 index ee759452e..c14bbacf8 100644 --- a/cgyro/src/cgyro_globals.F90 +++ b/cgyro/src/cgyro_globals.F90 @@ -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: @@ -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)' !---------------------------------------------------- !--------------------------------------------------------------- diff --git a/cgyro/src/cgyro_init_manager.F90 b/cgyro/src/cgyro_init_manager.F90 index 956fb3372..f557048c6 100644 --- a/cgyro/src/cgyro_init_manager.F90 +++ b/cgyro/src/cgyro_init_manager.F90 @@ -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 diff --git a/cgyro/src/cgyro_write_initdata.f90 b/cgyro/src/cgyro_write_initdata.f90 index 14ab40b6b..9285aa5b4 100644 --- a/cgyro/src/cgyro_write_initdata.f90 +++ b/cgyro/src/cgyro_write_initdata.f90 @@ -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 @@ -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 @@ -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 diff --git a/cgyro/src/cgyro_write_timedata.f90 b/cgyro/src/cgyro_write_timedata.f90 index e81cea301..135bfc412 100644 --- a/cgyro/src/cgyro_write_timedata.f90 +++ b/cgyro/src/cgyro_write_timedata.f90 @@ -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) @@ -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) @@ -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) diff --git a/f2py/pygacode/cgyro/data.py b/f2py/pygacode/cgyro/data.py index 5fde54b8d..739177c3a 100644 --- a/f2py/pygacode/cgyro/data.py +++ b/f2py/pygacode/cgyro/data.py @@ -6,8 +6,6 @@ from gacodefuncs import * except: from ..gacodefuncs import * - -BYTE='float' # class to read cgyro output data class cgyrodata: @@ -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() @@ -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=' ') @@ -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.') @@ -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') @@ -536,7 +533,6 @@ def getnorm(self,norm): except: pass - self.rhonorm = self.rho self.rhoi = r'\rho_s'