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)