import numpy as np from compact import compact class CompactScheme: def __init__ (self, nx, ny, nz, px, py, pz, lx, ly, lz): pi8 = np.arccos(-1., dtype=np.float64) self.shape = (nz, ny, nx) self.coefx1 = np.zeros((nx)) self.coefy1 = np.zeros((ny)) self.coefz1 = np.zeros((nz)) self.coefx2 = np.zeros((nx)) self.coefy2 = np.zeros((ny)) self.coefz2 = np.zeros((nz)) self.coefx3 = np.zeros((nx)) self.coefy3 = np.zeros((ny)) self.coefz3 = np.zeros((nz)) self.coefx4 = np.zeros((nx)) self.coefy4 = np.zeros((ny)) self.coefz4 = np.zeros((nz)) self.px = px self.py = py self.pz = pz h = pi8 * lx / nx self.hx = h self.hy = h self.hz = h 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 ddx (self, src): 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: for i in range(nz): xsrc[:] = src[i] dst[i] = compact.dfp(self.hx, xsrc, 1) else: for i in range(nz): xsrc[:] = src[i] dst[i] = compact.dfnonp(self.hx, xsrc, 1) # return np.swapaxes(dst, 1, 2) return dst def ddy (self): return def ddz (self): return 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): import struct import sys import os 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): import struct import sys import os 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.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[0] Y1 = s[1] return t, nx, ny, nz, u, v, w, Y0, Y1 def validate_trigonometric(): shape = (256, 256, 512) nz, ny, nx = shape pi8 = np.arccos(-1.) print (pi8) l_0 = 2.0 hyp=l_0*pi8/ny hxp=hyp hzp=hyp Y1 = np.zeros(shape) dY = np.zeros(shape) true = np.zeros(shape) XX = np.arange(nx) * hxp YY = np.arange(ny) * hyp ZZ = np.arange(nz) * hzp zz, yy, xx = np.meshgrid(ZZ, YY, XX) 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))[:] cs = CompactScheme(nx, ny, nz, False, True, True, 4., 2., 2.) 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[:] = true[:] dydx = np.memmap("dphi", dtype=np.float64, mode="w+", shape=cs.shape) dydx[:] = cs.ddx(Y1)[:] print (dydx.min(), dydx.max()) print (dydxtrue.min(), dydxtrue.max()) relerr = (dydx - dydxtrue) / dydxtrue print (np.nanmin(relerr), np.nanmax(relerr)) # cs.verify_nonp_coef() def test_dns_data(): import sys file_name = sys.argv[1] t, nx, ny, nz, u, v, w, Y0, Y1 = read_data(file_name) cs = CompactScheme(nx, ny, nz, False, True, True, 4, 2, 2) y = np.memmap("yr", dtype=np.float64, mode="w+", shape=cs.shape) y[:] = Y1[:] print (y.min(), y.max()) dydx = np.memmap("dyr", dtype=np.float64, mode="w+", shape=cs.shape) dydx[:] = cs.ddx(Y1)[:] dydx.flush() print (dydx.min(), dydx.max()) if __name__ == "__main__": validate_trigonometric() # test_dns_data()