255 lines
9.1 KiB
Python
255 lines
9.1 KiB
Python
import argparse
|
|
import datetime
|
|
import sys
|
|
import os
|
|
import numpy as np
|
|
import dnstool
|
|
from pycompact import CompactScheme
|
|
|
|
|
|
program_description = '''\
|
|
Read all DNS data files and collect c and ddx(c) at selected x coordinates.
|
|
Compute conditional mean of ddx(c) given condition c by data binning
|
|
'''
|
|
|
|
# Commandline argument parser
|
|
parser = argparse.ArgumentParser(description=program_description)
|
|
parser.add_argument("-c", "--case", help="target case name, required")
|
|
parser.add_argument("-l", "--list", help="list all case names", action='store_true')
|
|
parser.add_argument("-t", "--test", help="test run, use only 10 timesteps", action='store_true')
|
|
args = parser.parse_args()
|
|
params = vars(args)
|
|
|
|
cases = dnstool.case_library()
|
|
|
|
if params["list"] :
|
|
for c in cases:
|
|
print(c, cases[c].case_root)
|
|
sys.exit(0)
|
|
else:
|
|
if params["case"] is None :
|
|
sys.exit("--case CASE not given")
|
|
|
|
casename = params["case"]
|
|
case = cases[casename]
|
|
nx, ny, nz = case.shape
|
|
|
|
n_timestep = len(case.data_files)
|
|
if params["test"] :
|
|
n_timestep = 10
|
|
fdebug = True
|
|
|
|
x_indices_all = {
|
|
'IC1': np.array([272, 285, 294, 302, 309, 315, 321, 329, 339]),
|
|
'IC2': np.array([251, 273, 286, 298, 309, 320, 330, 342, 359]),
|
|
'IC3': np.array([246, 270, 285, 299, 311, 322, 333, 346, 363]),
|
|
'IC4': np.array([262, 279, 291, 300, 309, 317, 326, 337, 349]),
|
|
'IC8': np.array([254, 277, 291, 302, 312, 321, 330, 340, 354]),
|
|
'IC10': np.array([285, 294, 300, 304, 308, 312, 316, 320, 326]),
|
|
'IC11': np.array([264, 281, 292, 301, 309, 317, 325, 334, 347]),
|
|
'IC12': np.array([253, 276, 290, 301, 311, 320, 330, 340, 355]),
|
|
'IC13': np.array([243, 269, 287, 300, 312, 323, 334, 347, 365]),
|
|
'VIC1': np.array([286, 294, 299, 304, 308, 312, 315, 320, 326]),
|
|
'VIC2': np.array([264, 280, 291, 300, 309, 317, 325, 335, 348]),
|
|
'VIC4': np.array([273, 287, 295, 302, 309, 315, 321, 328, 338]),
|
|
'VIC10': np.array([292, 298, 302, 305, 308, 310, 313, 316, 320]),
|
|
'VIC13': np.array([253, 276, 291, 302, 312, 321, 331, 341, 355]),}
|
|
|
|
x_indices = x_indices_all[casename]
|
|
|
|
print("Calculating Conditional Mean of ddx(c) given c")
|
|
print(" * Case:", casename)
|
|
print(" * x index:", x_indices)
|
|
|
|
|
|
class ConditionalMeanBinning:
|
|
def __init__ (self, var='phi', case='default', loc=0, nbins=100, boundary=0.1):
|
|
self.eps = np.finfo(np.float).eps
|
|
|
|
self.var = var
|
|
self.case = case
|
|
self.loc = loc
|
|
self.name = "cavg_{}_c_{}_{}".format(var, case, loc)
|
|
|
|
# Point grid containing interval boundaries
|
|
self.grid = np.hstack(
|
|
(np.linspace(0, boundary, nbins),
|
|
np.linspace(boundary, 1-boundary, nbins)[1:],
|
|
np.linspace(1-boundary, 1, nbins)[1:],
|
|
)
|
|
)
|
|
self.grid[0] = -self.eps # to include c = 0 samples to the first bin
|
|
|
|
# interval center points
|
|
self.cstar = np.zeros(self.grid.size + 1)
|
|
self.cstar[1:-1] = (self.grid[1:]+self.grid[:-1])/2.
|
|
self.cstar[-1] = 1.0
|
|
|
|
self.v_sums = np.zeros(self.grid.size)
|
|
self.v_counts = np.zeros(self.grid.size, dtype=int)
|
|
|
|
def calculate_conditional_mask (self, c):
|
|
''' Feed arrays c, v of samples and conditions '''
|
|
|
|
bin_indices = np.digitize(c, self.grid, right=True)
|
|
# assign bin indices. bin intervals are right inclusive
|
|
# i-th bin: c_bin[i-1] < c <= c_bin[i]
|
|
|
|
return np.vstack([bin_indices == i for i in range(1, self.grid.size)]).reshape((-1, *bin_indices.shape))
|
|
|
|
def feed_samples (self, c, v, bin_masks=None):
|
|
''' Feed arrays c, v of samples and conditions '''
|
|
|
|
if bin_masks is None:
|
|
bin_masks = self.calculate_conditional_mask(c)
|
|
|
|
for i in range(1, self.grid.size):
|
|
binned = v[bin_masks[i-1]]
|
|
self.v_sums[i] += binned.sum()
|
|
self.v_counts[i] += binned.size
|
|
|
|
def average (self):
|
|
v_avg = np.zeros(self.v_sums.size+1)
|
|
v_avg[:-1] = self.v_sums / (self.v_counts + self.eps)
|
|
# add eps to prevent nan values from divide by zero
|
|
|
|
d = self.sum_count()
|
|
d["bin_cstar"] = self.cstar
|
|
d["bin_avg"] = v_avg
|
|
return d
|
|
|
|
def sum_count (self):
|
|
return {
|
|
"bin_grid" : self.grid,
|
|
"bin_sum" : self.v_sums,
|
|
"bin_count" : self.v_counts
|
|
}
|
|
|
|
def save_to_file (self):
|
|
np.savez(self.name, **self.average())
|
|
print("saved ", self.name)
|
|
|
|
def sutherland(c, dm0=0.02*4./3., rvis=5., pre=2.10e4, ac=26.7, bc=3.):
|
|
t0 = 1
|
|
t1 = 1 + bc
|
|
|
|
Ts = (rvis * t1 - (t1**(3. / 2.))) / (t1**(3. / 2.) - rvis)
|
|
As = t0 + Ts
|
|
|
|
theta = (1 + bc * c)
|
|
|
|
return dm0 * As * np.sqrt(theta) / (1. + Ts / theta)
|
|
|
|
|
|
cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2)
|
|
|
|
# c mean at sampling x coordinates (verification purpose)
|
|
cmean_verify = np.zeros(len(x_indices))
|
|
|
|
# Conditional Average Objects
|
|
avg_objs_ddxdmddxc = [ConditionalMeanBinning(var='ddxdmddxc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_ddydmddyc = [ConditionalMeanBinning(var='ddydmddyc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_ddzdmddzc = [ConditionalMeanBinning(var='ddzdmddzc', case=casename, loc=xi) for xi in x_indices]
|
|
|
|
avg_objs_ddxc = [ConditionalMeanBinning(var='ddxc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_ddyc = [ConditionalMeanBinning(var='ddyc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_ddzc = [ConditionalMeanBinning(var='ddzc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_d2dxc = [ConditionalMeanBinning(var='d2dxc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_d2dyc = [ConditionalMeanBinning(var='d2dyc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_d2dzc = [ConditionalMeanBinning(var='d2dzc', case=casename, loc=xi) for xi in x_indices]
|
|
|
|
avg_objs_uddxc = [ConditionalMeanBinning(var='uddxc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_vddyc = [ConditionalMeanBinning(var='vddyc', case=casename, loc=xi) for xi in x_indices]
|
|
avg_objs_wddzc = [ConditionalMeanBinning(var='wddzc', case=casename, loc=xi) for xi in x_indices]
|
|
|
|
targets=[ 'ddxc', 'ddyc', 'ddzc', 'd2dxc', 'd2dyc', 'd2dzc', 'uddxc', 'vddyc', 'wddzc', ]
|
|
mean_array = np.zeros((10, nx))
|
|
|
|
# Loop over dns data files
|
|
for i, fname in enumerate(case.data_files[:n_timestep]):
|
|
print(datetime.datetime.now(), '{:5d}/{}'.format(i+1, n_timestep), fname)
|
|
|
|
fpath = os.path.join(case.case_root, fname)
|
|
time, u0, shape, velocity, scalar = case.read_data(fpath)
|
|
u, v, w = velocity
|
|
u = u+u0
|
|
y = scalar[-1]
|
|
c = (1. - y)
|
|
dm = sutherland(c)
|
|
|
|
d2dxc = cs.d2dx(c)
|
|
d2dyc = cs.d2dy(c)
|
|
d2dzc = cs.d2dz(c)
|
|
|
|
ddxdm = cs.ddx(dm)
|
|
ddydm = cs.ddy(dm)
|
|
ddzdm = cs.ddz(dm)
|
|
|
|
ddzc = cs.ddz(c)
|
|
ddyc = cs.ddy(c)
|
|
ddxc = cs.ddx(c)
|
|
|
|
uddxc = u*ddxc
|
|
vddyc = v*ddyc
|
|
wddzc = w*ddzc
|
|
|
|
ddxdmddxc = ddxdm*ddxc
|
|
ddydmddyc = ddydm*ddyc
|
|
ddzdmddzc = ddzdm*ddzc
|
|
|
|
for idx, field in enumerate([u, v, w, c, ddxc, ddyc, ddzc, d2dxc, d2dyc, d2dzc]):
|
|
mean_array[idx] += field.sum(axis=(0,1))
|
|
|
|
for j, xi in enumerate(x_indices):
|
|
mask = avg_objs_ddxc[j].calculate_conditional_mask(c[:,:,xi])
|
|
|
|
avg_objs_ddxdmddxc[j].feed_samples(c[:,:,xi], ddxdmddxc[:,:,xi], mask)
|
|
avg_objs_ddydmddyc[j].feed_samples(c[:,:,xi], ddydmddyc[:,:,xi], mask)
|
|
avg_objs_ddzdmddzc[j].feed_samples(c[:,:,xi], ddzdmddzc[:,:,xi], mask)
|
|
|
|
avg_objs_ddxc[j].feed_samples(c[:,:,xi], ddxc[:,:,xi], mask)
|
|
avg_objs_ddyc[j].feed_samples(c[:,:,xi], ddyc[:,:,xi], mask)
|
|
avg_objs_ddzc[j].feed_samples(c[:,:,xi], ddzc[:,:,xi], mask)
|
|
|
|
avg_objs_d2dxc[j].feed_samples(c[:,:,xi], d2dxc[:,:,xi], mask)
|
|
avg_objs_d2dyc[j].feed_samples(c[:,:,xi], d2dyc[:,:,xi], mask)
|
|
avg_objs_d2dzc[j].feed_samples(c[:,:,xi], d2dzc[:,:,xi], mask)
|
|
|
|
avg_objs_uddxc[j].feed_samples(c[:,:,xi], uddxc[:,:,xi], mask)
|
|
avg_objs_vddyc[j].feed_samples(c[:,:,xi], vddyc[:,:,xi], mask)
|
|
avg_objs_wddzc[j].feed_samples(c[:,:,xi], wddzc[:,:,xi], mask)
|
|
|
|
cmean_verify[j] += np.sum(c[:,:,xi])
|
|
|
|
print(datetime.datetime.now(), 'Finished')
|
|
|
|
print()
|
|
print("Verify mean(c) at sampling x coordinates")
|
|
cmean_verify /= (n_timestep * ny * nz)
|
|
for xi_cmean_pair in zip(x_indices, cmean_verify):
|
|
print (xi_cmean_pair)
|
|
|
|
print()
|
|
print("Save the result")
|
|
|
|
np.save("avg_uvwc_dc_d2c_{}".format(casename), mean_array)
|
|
|
|
def save_cmean_result (fmt, obj_list):
|
|
for j, xi in enumerate(x_indices):
|
|
obj_list[j].save_to_file()
|
|
|
|
save_cmean_result("cavg_ddxdmddxc_c_{}_{}", avg_objs_ddxdmddxc)
|
|
save_cmean_result("cavg_ddydmddyc_c_{}_{}", avg_objs_ddydmddyc)
|
|
save_cmean_result("cavg_ddzdmddzc_c_{}_{}", avg_objs_ddzdmddzc)
|
|
|
|
save_cmean_result("cavg_ddxc_c_{}_{}", avg_objs_ddxc)
|
|
save_cmean_result("cavg_ddyc_c_{}_{}", avg_objs_ddyc)
|
|
save_cmean_result("cavg_ddzc_c_{}_{}", avg_objs_ddzc)
|
|
|
|
save_cmean_result("cavg_uddxc_c_{}_{}", avg_objs_uddxc)
|
|
save_cmean_result("cavg_vddyc_c_{}_{}", avg_objs_vddyc)
|
|
save_cmean_result("cavg_wddzc_c_{}_{}", avg_objs_wddzc)
|
|
|
|
save_cmean_result("cavg_d2dxc_c_{}_{}", avg_objs_d2dxc)
|
|
save_cmean_result("cavg_d2dyc_c_{}_{}", avg_objs_d2dyc)
|
|
save_cmean_result("cavg_d2dzc_c_{}_{}", avg_objs_d2dzc)
|