all fixed Dm balance terms

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
-
+    fdebug = True
-u_idx = 0
-y_idx = 4
 
 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_grid" : self.grid,
-            "bin_sum" : self.v_sums, 
+            "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]
+# Conditional Average Objects
-avg_objs_d2dxc = [ConditionalMeanBinning(name=str(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]
+
+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])
+        mask = avg_objs_ddxc[j].calculate_conditional_mask(c[:,:,xi])
-        avg_objs_d2dxc[j].feed_samples(c[:,:,xi], d2dxc[:,:,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):
+np.save("avg_uvwc_dc_d2c_{}".format(casename), mean_array)
-    filename = "cavg_d2dxc_c_{}_{}".format(casename, avg_objs_d2dxc[j].name)
+
-    np.savez(filename, **avg_objs_d2dxc[j].average())
+def save_cmean_result (fmt, obj_list):
-    print("saved ", filename)
+    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)

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):