import os import sys import struct import pprint import numpy as np from compact import compact class CompactScheme: """Python wrapper for the high-order compact finite difference scheme core. This wraps the compiled Fortran `compact` solver module, handling grid configurations, boundary periodicities, array allocations, and LU decompositions. It exposes differentiation methods ddx, ddy, and ddz to Python. """ def __init__ (self, nx, ny, nz, px, py, pz, lx, ly, lz): """Initializes the CompactScheme solver. Args: nx (int): Grid points in X direction. ny (int): Grid points in Y direction. nz (int): Grid points in Z direction. px (bool): Periodic boundary condition flag for X direction. py (bool): Periodic boundary condition flag for Y direction. pz (bool): Periodic boundary condition flag for Z direction. lx (float): Domain size in X direction. ly (float): Domain size in Y direction. lz (float): Domain size in Z direction. """ pi8 = np.arccos(-1., dtype=np.float64) self.shape = (nz, ny, nx) self.px = px self.py = py self.pz = pz h = pi8 * lx / nx self.hx = h self.hy = h self.hz = h # Allocate LU compact.lxf = np.zeros(nx, dtype=np.float64) compact.lxs = np.zeros(nx, dtype=np.float64) compact.wxf = np.zeros(nx, dtype=np.float64) compact.wxs = np.zeros(nx, dtype=np.float64) compact.lyf = np.zeros(ny, dtype=np.float64) compact.lys = np.zeros(ny, dtype=np.float64) compact.wyf = np.zeros(ny, dtype=np.float64) compact.wys = np.zeros(ny, dtype=np.float64) compact.lzf = np.zeros(nz, dtype=np.float64) compact.lzs = np.zeros(nz, dtype=np.float64) compact.wzf = np.zeros(nz, dtype=np.float64) compact.wzs = np.zeros(nz, dtype=np.float64) bcx = 0 if px else 1 bcy = 0 if py else 1 bcz = 0 if pz else 1 compact.ludcmp_calculate(nx, ny, nz, bcx, bcy, bcz) def test_ludcmp (self): pp = pprint.PrettyPrinter(indent=4) # First Derivative Non-periodic BC l1 = compact.test_nonp_lud1(self.shape[-1]) # Second Derivative Non-periodic BC l2 = compact.test_nonp_lud2(self.shape[-1]) print ("Test Internally Calculated Non-periodic Coefs") print (np.linalg.norm((l1 - compact.lxf)/compact.lxf)) print (np.linalg.norm((l2 - compact.lxs)/compact.lxs)) def py_rhs_1_np (self, x): dx = np.zeros(x.shape) h1 = 1./self.hx r1 = 7./3. r2 = 1./12. r3 = 3. a = -1.25 b = 1. c = 0.25 nd, n = x.shape dx[:, -2] = x[:, -1] - x[:, -3] dx[:, -1] = - (a*x[:, -1] + b*x[:, -2] + c*x[:, -3]) dx[:, 0] = (a*x[:, 0] + b*x[:, 1] + c*x[:, 2]) dx[:, 1] = x[:, 2] - x[:, 0] dx[:,-2] = dx[:,-2]*h1*r3 dx[:,-1] = dx[:,-1]*h1 dx[:,0] = dx[:,0]*h1 dx[:,1] = dx[:,1]*h1*r3 for i in range(2,n-2): t1=x[:,i+1]-x[:,i-1] t2=x[:,i+2]-x[:,i-2] dx[:,i]=h1*(r1*t1+r2*t2) return dx def py_tdslv(self, r, l): nd, n = r.shape r[:,0] = r[:,0] * l[0] for i in range(1,n): r[:,i] = l[i] * (r[:,i] - r[:,i-1]) for i in range(n-1)[::-1]: r[:,i] = r[:,i] - l[i] * r[:,i+1] def test_dfnonp (self): x = np.sin(1.1 * np.arange(512) * self.hx).reshape((1,-1)) exact = 1.1 * np.cos(1.1 * np.arange(512) * self.hx).reshape((1,-1)) print ("First Non-periodic RHS Test") dx = self.py_rhs_1_np(x) dx_fortran = compact.rhs1np(self.hx, x) print ("RelError Norm: ", np.linalg.norm((dx - dx_fortran) / dx_fortran)) print ("RelError Min : ", ((dx - dx_fortran) / dx_fortran).min()) print ("RelError Min : ", ((dx - dx_fortran) / dx_fortran).max()) print ("First Non-periodic TD SOLVE Test") l1 = compact.test_nonp_lud1(512) self.py_tdslv(dx, l1) compact.tdslv(dx_fortran,l1) print ("dx - exact") print ("RelError Norm: ", np.linalg.norm((dx - exact) / exact)) print ("RelError Min : ", ((dx - exact) / exact).min()) print ("RelError Min : ", ((dx - exact) / exact).max()) print ("dx_fortran - exact") print ("RelError Norm: ", np.linalg.norm((dx_fortran - exact) / exact)) print ("RelError Min : ", ((dx_fortran - exact) / exact).min()) print ("RelError Min : ", ((dx_fortran - exact) / exact).max()) ''' import pprint pp = pprint.PrettyPrinter(indent=4) pp.pprint ((dx - exact) / exact) pp.pprint ((zip ((dx - exact).ravel(), dx.ravel(), exact.ravel()))) ''' def verify_nonp_lud1(self): print ("Non-periodic coef first derivative") nx = 512 aa = np.ones(nx) * 3. aa[0] = 0.5 aa[1] = 4. aa[-2] = 4. aa[-1] = 0.5 coef = compact.stdlu(aa) coef_verify = self.py_stdlu(aa) print ("RelError Norm: ", np.linalg.norm((coef - coef_verify)/coef_verify)) def verify_nonp_lud2(self): print ("Non-periodic coef second derivative") nx = 512 aa = np.ones(nx) * 3. aa[0] = 2./11. aa[1] = 10. aa[-2] = 10. aa[-1] = 2./11. coef = compact.stdlu(aa) coef_verify = self.py_stdlu(aa) print ("RelError Norm: ", np.linalg.norm((coef - coef_verify)/coef_verify)) def py_stdlu(self, aa): coef = np.ones(aa.shape)/aa[0] print ("coef.size = ", coef.size) for i in range(1,coef.size): coef[i]=1.0/(aa[i]-coef[i-1]) return coef def ddx (self, src): """Computes the first-order spatial derivative in the X direction. Args: src (numpy.ndarray): 3D input field array matching (nz, ny, nx) shape. Returns: numpy.ndarray: 3D first derivative array in the X direction. """ if src.shape != self.shape: print ("error") 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 for i in range(nz): dst[i] = compact.dfp(self.hx, src[i], 1) else: for i in range(nz): dst[i] = compact.dfnonp(self.hx, src[i], 1) # return np.swapaxes(dst, 1, 2) return dst def ddy (self, src): """Computes the first-order spatial derivative in the Y direction. Args: src (numpy.ndarray): 3D input field array matching (nz, ny, nx) shape. Returns: numpy.ndarray: 3D first derivative array in the Y direction. """ if src.shape != self.shape: print ("error") 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 for i in range(nz): dst[i] = compact.dfp(self.hx, src[i].T, 2).T else: for i in range(nz): dst[i] = compact.dfnonp(self.hx, src[i].T, 2).T # return np.swapaxes(dst, 1, 2) return dst def ddz (self, src): """Computes the first-order spatial derivative in the Z direction. Args: src (numpy.ndarray): 3D input field array matching (nz, ny, nx) shape. Returns: numpy.ndarray: 3D first derivative array in the Z direction. """ if src.shape != self.shape: print ("error") 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) else: for i in range(ny): dst[:,i,:] = compact.dfnonp(self.hx, src[:,i,:], 3) # return np.swapaxes(dst, 1, 2) return dst def port_nonp_coef (self): # SUBROUTINE nonp_lud(xyz,xx) nz, ny, nx = self.shape xx = nx lxf = np.zeros(xx) lxs = np.zeros(xx) aa = np.zeros(xx) aa[:] = 3. aa[0]=0.5 aa[1]=4. aa[-2]=4. aa[-1]=0.5 # first derivative compact.stdlu(aa,lxf) aa[:] = 5.5 aa[0]=2./11. aa[1]=10. aa[-2]=10. aa[-1]=2./11. # second derivative compact.stdlu(aa,lxs) compact.lxf = lxf compact.lxs = lxs def read_old_data (fname): with open(fname, 'rb') as f1 : f1.seek(0) raw_info = f1.read(4+8*6+4)[4:-4] t = struct.unpack('d', raw_info[ 0: 8])[0] nx = struct.unpack('q', raw_info[ 8:16])[0] ny = struct.unpack('q', raw_info[16:24])[0] nz = struct.unpack('q', raw_info[24:32])[0] count = nx*ny*nz bSize = count*8 # size in bytes for a variable dummy_len = (4+8*3+4) + (4+8*2+4) + (4+8*2+4) + (4+8*2+4) + 4 dummy = f1.read(dummy_len) #dummy = f1.read(4) print (t, nx, ny, nz) #raw_field = f1.read(4+bSize*5+4)[4:-4] V = np.fromfile(f1, dtype=np.float64, count=(3*count)).reshape((3,nz,ny,nx)) s = np.fromfile(f1, dtype=np.float64, count=(2*count)).reshape((nz,ny,nx,2)) print (V.order) print (s.order) print (V.shape) print (s.shape) V.order="F" s.order="F" print (V.shape) print (s.shape) u = V[0] v = V[1] w = V[2] Y0 = s.T[0].T Y1 = s.T[1].T return t, nx, ny, nz, u, v, w, Y0, Y1 def read_data (fname): with open(fname, 'rb') as f1 : f1.seek(0) raw_info = f1.read(4+8*6+4)[4:-4] t = struct.unpack('d', raw_info[ 0: 8])[0] nx = struct.unpack('q', raw_info[ 8:16])[0] ny = struct.unpack('q', raw_info[16:24])[0] nz = struct.unpack('q', raw_info[24:32])[0] count = nx*ny*nz bSize = count*8 # size in bytes for a variable dummy_len = (4+8*3+4) + (4+8*2+4) + (4+8*2+4) + (4+8*2+4) + 4 dummy = f1.read(dummy_len) #dummy = f1.read(4) #raw_field = f1.read(4+bSize*5+4)[4:-4] V = np.fromfile(f1, dtype=np.float64, count=(3*count)).reshape((3,nz,ny,nx)) s = np.fromfile(f1, dtype=np.float64, count=(2*count)).reshape((2,nz,ny,nx)) print (V.flags) print (s.flags) print (V.shape) print (s.shape) u = V[0] v = V[1] w = V[2] Y0 = s[0] Y1 = s[1] return t, nx, ny, nz, u, v, w, Y0, Y1 def validate_trigonometric(): writeToFile = False shape = (256, 256, 512) nz, ny, nx = shape pi8 = np.arccos(-1.) l_0 = 2.0 hyp=l_0*pi8/ny hxp=hyp hzp=hyp cs = CompactScheme(nx, ny, nz, False, True, True, 4., 2., 2.) cs.test_ludcmp() cs.verify_nonp_lud1() cs.verify_nonp_lud2() cs.test_dfnonp() print ("Test ddx") Y1 = np.zeros(shape) true = np.zeros(shape) XX = np.arange(nx) * hxp YY = np.arange(ny) * hyp ZZ = np.arange(nz) * hzp print ("1-D sine test") cos_fortran = compact.dfnonp(hxp, np.sin(1.1*XX).reshape((1,-1)), 1) cos_exact = 1.1 * np.cos(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) Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] def compare_3d_result (true1, dY1): print ("Calculated Min/Max", dY1.min(), dY1.max()) print ("True Min/Max", true1.min(), true1.max()) eps = np.finfo(true1.dtype).eps relerr = (dY1 - true1) / (true1 + eps) print ("Relative Error", np.nanmin(relerr), np.nanmax(relerr)) print(" DDX Test ") 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(" 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))[:] dY1 = cs.ddy(Y1)[:] compare_3d_result(true, dY1) 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))[:] dY1 = cs.ddz(Y1)[:] compare_3d_result(true, dY1) if writeToFile: y = np.memmap("phi", dtype=np.float64, mode="w+", shape=cs.shape) y[:] = Y1[:] dydxtrue = np.memmap("dphitrue", dtype=np.float64, mode="w+", shape=cs.shape) dydxtrue[:] = true1[:] dydx = np.memmap("dphi", dtype=np.float64, mode="w+", shape=cs.shape) dydx[:] = dY1[:] # cs.verify_nonp_coef() def test_dns_data(): file_name = "./fort.1000" answer = "./fort.2000" t, nx, ny, nz, u, v, w, Y0, Y1 = read_data(file_name) with open(answer, 'rb') as ans_file: ans_file.seek(4) ddx_answer = np.fromfile(ans_file, dtype=np.float64, count=Y1.size, ).reshape(Y1.shape) cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2) ddx = cs.ddx(Y1) print (ddx.min(), ddx.max()) print (ddx_answer.min(), ddx_answer.max()) relerr = (ddx - ddx_answer) / (ddx_answer) err = (ddx - ddx_answer) print ("Absolute Error", np.nanmin(err), np.nanmax(err)) print ("Relative Error", np.nanmin(relerr), np.nanmax(relerr)) if __name__ == "__main__": # validate_trigonometric() print("DNS Field Test") test_dns_data()