From 44ef8da3b92caae3ad0ea522d024458e244f3829 Mon Sep 17 00:00:00 2001 From: ignis Date: Thu, 12 Aug 2021 03:41:32 +0900 Subject: [PATCH] binning 2nd order derivative d2/dx2(c) and pycompact complete interface/bugfix --- binning_ddxc_d2dxc_given_c_from_raw.py | 165 +++++++++++++++++++++++++ pycompact.py | 125 ++++++++++++++++--- 2 files changed, 273 insertions(+), 17 deletions(-) create mode 100644 binning_ddxc_d2dxc_given_c_from_raw.py diff --git a/binning_ddxc_d2dxc_given_c_from_raw.py b/binning_ddxc_d2dxc_given_c_from_raw.py new file mode 100644 index 0000000..da95073 --- /dev/null +++ b/binning_ddxc_d2dxc_given_c_from_raw.py @@ -0,0 +1,165 @@ +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 + +u_idx = 0 +y_idx = 4 + +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, name='', nbins=100, boundary=0.1): + self.eps = np.finfo(np.float).eps + + self.name = name + + # 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 feed_samples (self, c, v): + ''' 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] + + for i in range(1, self.grid.size): + binned = v[bin_indices == i] + 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 + } + +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_u = [ConditionalMeanBinning(name=str(xi)) for xi in x_indices] +avg_objs_ddxc = [ConditionalMeanBinning(name=str(xi)) for xi in x_indices] +avg_objs_uddxc = [ConditionalMeanBinning(name=str(xi)) for xi in x_indices] + +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, u = case.read_single_field(fpath, u_idx) + time, y = case.read_single_field(fpath, y_idx) + c = (1. - y) + ddxc = cs.ddx(c) + + for j, xi in enumerate(x_indices): + avg_objs_u[j].feed_samples(c[:,:,xi], u[:,:,xi]) + avg_objs_ddxc[j].feed_samples(c[:,:,xi], ddxc[:,:,xi]) + avg_objs_uddxc[j].feed_samples(c[:,:,xi], u[:,:,xi]*ddxc[:,:,xi]) + 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") +for j, xi in enumerate(x_indices): + filename = "cavg_u_c_{}_{}".format(casename, avg_objs_u[j].name) + np.savez(filename, **avg_objs_u[j].average()) + print("saved ", filename) + +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_uddxc_c_{}_{}".format(casename, avg_objs_uddxc[j].name) + np.savez(filename, **avg_objs_uddxc[j].average()) + print("saved ", filename) diff --git a/pycompact.py b/pycompact.py index 0b46dab..537def9 100644 --- a/pycompact.py +++ b/pycompact.py @@ -202,9 +202,6 @@ class CompactScheme: nz, ny, nx = self.shape - xsrc = np.zeros((ny, nx,), dtype=np.float64, order="F") - - # dst = np.zeros((nx, ny, nz,), order="F") dst = np.zeros((nz, ny, nx,), dtype=np.float64,) if self.px: # Periodic BC @@ -225,9 +222,6 @@ class CompactScheme: nz, ny, nx = self.shape - #xsrc = np.zeros((ny, nx,), dtype=np.float64, order="F") - - # dst = np.zeros((nx, ny, nz,), order="F") dst = np.zeros((nz, ny, nx,), dtype=np.float64,) if self.py: # Periodic BC @@ -248,21 +242,71 @@ class CompactScheme: nz, ny, nx = self.shape - # dst = np.zeros((nx, ny, nz,), order="F") dst = np.zeros((nz, ny, nx,), dtype=np.float64,) if self.pz: # Periodic BC for i in range(ny): - dst[:,i,:] = compact.dfp(self.hx, src[:,i,:], 3) - + dst[:,i,:] = compact.dfp(self.hx, src[:,i,:].T, 3).T else: for i in range(ny): - dst[:,i,:] = compact.dfnonp(self.hx, src[:,i,:], 3) + dst[:,i,:] = compact.dfnonp(self.hx, src[:,i,:].T, 3).T # return np.swapaxes(dst, 1, 2) return dst + def d2dx (self, src): + + if src.shape != self.shape: print ("error") + + nz, ny, nx = self.shape + + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + if self.px: # Periodic BC + for i in range(nz): + dst[i] = compact.d2fp(self.hx, src[i], 1) + else: + for i in range(nz): + dst[i] = compact.d2fnonp(self.hx, src[i], 1) + + return dst + + def d2dy (self, src): + + if src.shape != self.shape: print ("error") + + nz, ny, nx = self.shape + + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + if self.py: # Periodic BC + for i in range(nz): + dst[i] = compact.d2fp(self.hx, src[i].T, 2).T + else: + for i in range(nz): + dst[i] = compact.d2fnonp(self.hx, src[i].T, 2).T + + return dst + + def d2dz (self, src): + + if src.shape != self.shape: print ("error") + + nz, ny, nx = self.shape + + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + if self.pz: # Periodic BC + for i in range(ny): + dst[:,i,:] = compact.d2fp(self.hx, src[:,i,:].T, 3).T + else: + for i in range(ny): + dst[:,i,:] = compact.d2fnonp(self.hx, src[:,i,:].T, 3).T + + return dst + + def port_nonp_coef (self): # SUBROUTINE nonp_lud(xyz,xx) @@ -437,11 +481,23 @@ def validate_trigonometric(): + print ("1-D sine test 2") + cos_fortran = compact.d2fnonp(hxp, np.sin(1.1*XX).reshape((1,-1)), 1) + cos_exact = - 1.1 * 1.1 * np.sin(1.1*XX).reshape((1,-1)) + + print ("Compact Scheme: ", cos_fortran.min(), cos_fortran.max()) + + print ("Exact : ",cos_exact.min(), cos_exact.max()) + + print ("Norm of relative errors: ", np.linalg.norm((cos_fortran - cos_exact)/cos_exact)) + # print (((cos_fortran - cos_exact)/cos_exact)) + + + print ("3-D trigonometric test") - zz, yy, xx = np.meshgrid(ZZ, YY, XX) + zz, yy, xx = np.meshgrid(ZZ, YY, XX, indexing='ij') - Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] def compare_3d_result (true1, dY1): @@ -455,15 +511,24 @@ def validate_trigonometric(): print ("Relative Error", np.nanmin(relerr), np.nanmax(relerr)) print(" DDX Test ") + Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] true[:] = (1.1 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] dY1 = cs.ddx(Y1)[:] compare_3d_result(true, dY1) + print("D2DX Test ") + Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] + true[:] = (- 1.1 * 1.1 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] + dY1 = cs.d2dx(Y1)[:] + + compare_3d_result(true, dY1) + + print(" DDY Test ") Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] - true[:] = (-3.0 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] + true[:] = (3.0 * np.sin(1.1 * xx) * np.cos(3.0 * yy) * np.sin(2.0 * zz))[:] dY1 = cs.ddy(Y1)[:] compare_3d_result(true, dY1) @@ -471,7 +536,7 @@ def validate_trigonometric(): print(" DDZ Test ") Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] - true[:] = (-2.0 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] + true[:] = (2.0 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.cos(2.0 * zz))[:] dY1 = cs.ddz(Y1)[:] compare_3d_result(true, dY1) @@ -519,10 +584,36 @@ def test_dns_data(): if __name__ == "__main__": + import argparse - # validate_trigonometric() + program_description = '''\ + Test Python interface to Fortran Compact Scheme Module" + ''' + # argparse enumerated values + + # Commandline argument parser + parser = argparse.ArgumentParser(description=program_description) + parser.add_argument("-t", "--test", help="test to be run") + parser.add_argument("-l", "--list", help="list all tests available", action="store_true") + args = parser.parse_args() + params = vars(args) - print("DNS Field Test") - test_dns_data() + # Parameters + testname = params["test"] + list_tests = params["list"] + + if not list_tests and testname is None: + print ("test not selected!") + exit(-1) + + + if testname == "tri": + print("Trigonometric Function Test") + validate_trigonometric() + + + if testname == "dns": + print("DNS Field Test") + test_dns_data()