From d244ac4a4356d8ec573ac04f7768bb628cb68ab1 Mon Sep 17 00:00:00 2001 From: ignis Date: Tue, 24 Aug 2021 10:21:22 +0900 Subject: [PATCH] all fixed Dm balance terms --- binning_balance_terms_from_raw.py | 119 ++++++++++++++++++++++++------ dnstool.py | 3 +- 2 files changed, 98 insertions(+), 24 deletions(-) diff --git a/binning_balance_terms_from_raw.py b/binning_balance_terms_from_raw.py index d5c82d6..b7dab4e 100644 --- a/binning_balance_terms_from_raw.py +++ b/binning_balance_terms_from_raw.py @@ -37,9 +37,7 @@ nx, ny, nz = case.shape n_timestep = len(case.data_files) if params["test"] : n_timestep = 10 - -u_idx = 0 -y_idx = 4 + fdebug = True x_indices_all = { 'IC1': np.array([272, 285, 294, 302, 309, 315, 321, 329, 339]), @@ -65,10 +63,13 @@ print(" * x index:", x_indices) class ConditionalMeanBinning: - def __init__ (self, name='', nbins=100, boundary=0.1): + def __init__ (self, var='phi', case='default', loc=0, nbins=100, boundary=0.1): self.eps = np.finfo(np.float).eps - self.name = name + 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( @@ -87,15 +88,23 @@ class ConditionalMeanBinning: self.v_sums = np.zeros(self.grid.size) self.v_counts = np.zeros(self.grid.size, dtype=int) - def feed_samples (self, c, v): + 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_indices == i] + binned = v[bin_masks[i-1]] self.v_sums[i] += binned.sum() self.v_counts[i] += binned.size @@ -111,32 +120,85 @@ class ConditionalMeanBinning: def sum_count (self): return { - "bin_grid" : self.grid, - "bin_sum" : self.v_sums, + "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) + + 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)) -avg_objs_ddxc = [ConditionalMeanBinning(name=str(xi)) for xi in x_indices] -avg_objs_d2dxc = [ConditionalMeanBinning(name=str(xi)) for xi in x_indices] +# Conditional Average Objects +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] + +avg_objs_ddxcddxc = [ConditionalMeanBinning(var='ddxcddxc', 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, y = case.read_single_field(fpath, y_idx) + time, u0, shape, velocity, scalar = case.read_data(fpath) + u, v, w = velocity + u = u+u0 + y = scalar[-1] c = (1. - y) - ddxc = cs.ddx(c) + d2dxc = cs.d2dx(c) + d2dyc = cs.d2dy(c) + d2dzc = cs.d2dz(c) + + ddzc = cs.ddz(c) + ddyc = cs.ddy(c) + ddxc = cs.ddx(c) + + ddxcddxc = ddxc * ddxc + + uddxc = u*ddxc + vddyc = v*ddyc + wddzc = w*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): - avg_objs_ddxc[j].feed_samples(c[:,:,xi], ddxc[:,:,xi]) - avg_objs_d2dxc[j].feed_samples(c[:,:,xi], d2dxc[:,:,xi]) + mask = avg_objs_ddxc[j].calculate_conditional_mask(c[:,:,xi]) + + 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_ddxcddxc[j].feed_samples(c[:,:,xi], ddxcddxc[:,:,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() @@ -147,12 +209,23 @@ for xi_cmean_pair in zip(x_indices, cmean_verify): print() print("Save the result") -for j, xi in enumerate(x_indices): - filename = "cavg_ddxc_c_{}_{}".format(casename, avg_objs_ddxc[j].name) - np.savez(filename, **avg_objs_ddxc[j].average()) - print("saved ", filename) -for j, xi in enumerate(x_indices): - filename = "cavg_d2dxc_c_{}_{}".format(casename, avg_objs_d2dxc[j].name) - np.savez(filename, **avg_objs_d2dxc[j].average()) - print("saved ", filename) +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_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) + +save_cmean_result("cavg_ddxcddxc_c_{}_{}", avg_objs_ddxcddxc) diff --git a/dnstool.py b/dnstool.py index 428a08c..c5d1a29 100644 --- a/dnstool.py +++ b/dnstool.py @@ -131,6 +131,7 @@ class DnsCase: nx = info_dict['nx'] ny = info_dict['ny'] nz = info_dict['nz'] + u0 = info_dict['u0'] with open(fname, 'rb') as f1 : f1.seek(offset_data) @@ -140,7 +141,7 @@ class DnsCase: V = np.fromfile(f1, dtype=np.double, count=(3*count)).reshape((3,nz,ny,nx)) s = np.fromfile(f1, dtype=np.double, count=(2*count)).reshape((2,nz,ny,nx)) - return t, (nx, ny, nz), V, s + return t, u0, (nx, ny, nz), V, s def read_single_field (self, fname, ifield):