-
Notifications
You must be signed in to change notification settings - Fork 21
/
extractTargetFilesNonDim.py
executable file
·168 lines (149 loc) · 6.69 KB
/
extractTargetFilesNonDim.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
#!/usr/bin/env python3
import re, argparse, numpy as np, glob, os
from os import path
#nBins = 24 * 16//2 - 1
'''
def tkeFit(nu, eps): return 2.87657077 * np.power(eps, 2/3.0)
def relFit(nu, eps):
tke = tkeFit(nu,eps)
uprime = np.sqrt(2.0/3.0 * tke);
lambd = np.sqrt(15 * nu / eps) * uprime;
return uprime * lambd / nu;
'''
def epsNuFromRe(Re, uEta = 1.0):
C = 3 # np.sqrt(20.0/3)
K = 2/3.0 * C * np.sqrt(15)
eps = np.power(uEta*uEta * Re / K, 3.0/2.0)
nu = np.power(uEta, 4) / eps
return eps, nu
def findAllParams(path):
REs = set()
alldirs = glob.glob(path+'*')
for dirn in alldirs: REs.add(re.findall('RE\d\d\d', dirn)[0][2:])
REs = list( REs); REs.sort()
for i in range(len(REs)): REs[i] = float(REs[i])
return REs
def computeIntTimeScale(tau_integral):
tau_int = 0.0
N = len(tau_integral)
for i in range(N):
if(tau_integral[i]<1e16): tau_int += tau_integral[i]/N
else: tau_int += 1e16/N
return tau_int
def readAnalysisFile(dirn, analysis, size):
fname = dirn + '/' + analysis
rl_fname = dirn + '/simulation_000_00000/run_00000000/' + analysis
if path.exists( fname):
f = np.fromfile(os.path.realpath( fname), dtype=np.float64)
elif path.exists(rl_fname):
f = np.fromfile(os.path.realpath(rl_fname), dtype=np.float64)
else :
f = np.zeros([size])
print('File %s does not exist. Simulation did not start?' % fname)
nFullRows = f.size // size
nFull = nFullRows * size
assert(nFull == f.size)
#return f.reshape([nFullRows, size])
return f[:nFull].reshape([nFullRows, size])
def getAllData(dirn, eps, nu, nBins, fSkip=1, readFlux=False):
f = readAnalysisFile(dirn, 'spectralAnalysis.raw', nBins + 13)
nSamples = f.shape[0]
tIntegral = computeIntTimeScale(f[:,8])
tAnalysis = np.sqrt(nu / eps) # time space between data files
#ind0 = int(10 * tIntegral / tAnalysis) # skip initial integral times
ind0 = int(5 * tIntegral / tAnalysis) # skip initial integral times
if ind0 == 0 or ind0 > nSamples: ind0 = nSamples
#print(ind0)
data = {
'dt' : f[ind0:, 1], 'tke' : f[ind0:, 3],
'tke_filtered' : f[ind0:, 4], 'dissip_visc' : f[ind0:, 5],
'dissip_tot' : f[ind0:, 6], 'l_integral' : f[ind0:, 7],
't_integral' : f[ind0:, 8], 'grad_mean' : f[ind0:,11],
'grad_std' : f[ind0:,12], 'spectra' : f[ind0:,13:]
}
#print('n spectra = %d %f %f' % (nFullRows, tIntegral, tAnalysis))
if not readFlux: return data
f = readAnalysisFile(dirn, 'fluxAnalysis.raw', nBins)
data['flux'] = f[ind0:, :]
return data
def gatherAllData(path, re, eps, nu, nBins, fSkip):
dirn = '%s_RE%03d' % (path, re)
data = getAllData(dirn, eps, nu, nBins, fSkip)
# REMEMBER: THIS IS NEEDED IF SCANNING MULTIPLE DNS FOLDERS:
#data = None
#for run in range(40):
# dirn = '%sRE%04d_RUN%d' % (path, re, run)
# runData = getAllData(dirn, eps, nu, nBins, fSkip)
# if data is None: data = runData
# else:
# for key in runData:
# data[key] = np.append(data[key], runData[key], 0)
if data['dt'].size < 2: data = None
return data
def main(path, fSkip, nBlocks, nBlocksRL):
nBinsTgt = nBlocksRL * 16 // 2 - 1
nBins = nBlocks * 16//2 - 1
REs = findAllParams(path)
#REs = [60, 70, 82, 95, 110, 130, 150, 176, 206, 240, 280, 325, 380]
EPSs, NUs = len(REs) * [0], len(REs) * [0] # will be overwritten
for j in range(len(REs)):
EPSs[j], NUs[j] = epsNuFromRe(REs[j])
print('Re %e nu %e eps %e' % (REs[j], NUs[j], EPSs[j]))
data = gatherAllData(path, REs[j], EPSs[j], NUs[j], nBins, fSkip=fSkip)
if data == None:
print('skipped eps:%f nu:%f' % (eps, nu))
continue
logE = np.log(data['spectra'])
print(logE.shape)
avgLogSpec = np.mean(logE, axis=0)
logE = np.log(data['spectra'][:, 0:nBinsTgt])
stdLogSpec = np.std(logE, axis=0)
covLogSpec = np.cov(logE, rowvar=False)
#print(covLogSpec.shape)
modes = np.arange(1, nBins+1) # assumes box is 2 pi
avgTke, avgDissip = np.mean(data['tke']), np.mean(data['dissip_tot'])
#reLambda = np.sqrt(20/3) * avgTke / np.sqrt(NUs[j] * avgDissip)
reLambda = np.sqrt(20/3) * avgTke / np.sqrt(NUs[j] * EPSs[j])
logCov2pi = np.power( np.linalg.det(2 * np.pi * covLogSpec), 0.5/nBinsTgt)
print(-np.log(logCov2pi))
fout = open('scalars_RE%03d' % int(REs[j]), "w")
fout.write("eps %e\n" % (EPSs[j]) )
fout.write("nu %e\n" % (NUs[j]) )
fout.write("dt %e %e\n" % (np.mean(data['dt']),
np.std( data['dt']) ) )
fout.write("tKinEn %e %e\n" % (np.mean(data['tke']),
np.std( data['tke']) ) )
fout.write("epsVis %e %e\n" % (np.mean(data['dissip_visc']),
np.std( data['dissip_visc']) ) )
fout.write("epsTot %e %e\n" % (np.mean(data['dissip_tot']),
np.std( data['dissip_tot']) ) )
fout.write("lInteg %e %e\n" % (np.mean(data['l_integral']),
np.std( data['l_integral']) ) )
fout.write("tInteg %e %e\n" % (np.mean(data['t_integral']),
np.std( data['t_integral']) ) )
fout.write("avg_Du %e %e\n" % (np.mean(data['grad_mean']),
np.std( data['grad_mean']) ) )
fout.write("std_Du %e %e\n" % (np.mean(data['grad_std']),
np.std( data['grad_std']) ) )
fout.write("ReLamd %e\n" % reLambda)
fout.write("logPdenom %e\n" % -np.log(logCov2pi) )
#print(modes.shape, avgLogSpec.shape)
ary = np.append( modes.reshape([nBins,1]),
avgLogSpec.reshape([nBins,1]), 1)
np.savetxt('spectrumLogE_RE%03d' % int(REs[j]), ary, delimiter=', ')
invCovLogSpec = np.linalg.inv(covLogSpec)
np.savetxt('invCovLogE_RE%03d' % int(REs[j]), invCovLogSpec, delimiter=', ')
np.savetxt('stdevLogE_RE%03d' % int(REs[j]), stdLogSpec, delimiter=', ')
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description = "Compute a target file for RL agent from DNS data.")
parser.add_argument('simdir',
help="Simulation directory containing the 'Analysis' folder")
parser.add_argument('--fSkip', type=int, default=1,
help="Sampling frequency for analysis files. If 1, take all. If 2, take 1 skip 1, If 3, take 1, skip 2, and so on.")
parser.add_argument('--nBlocks', type=int, default=32,
help="Number of CubismUP 3D blocks in the target runs.")
parser.add_argument('--nBlocksRL', type=int, default=2,
help="Number of CubismUP 3D blocks in the training runs.")
args = parser.parse_args()
main(args.simdir, args.fSkip, args.nBlocks, args.nBlocksRL)