From 22e6fa615dc30ca91023ac6345d38adf3144e9ba Mon Sep 17 00:00:00 2001 From: ignis Date: Sat, 6 Jul 2019 17:33:26 +0900 Subject: [PATCH] python wrapper using f2py --- code/pycompact/compact_doc.py | 2 + code/pycompact/pycompact.py | 300 ++++++++++++++++++++++++++++++++++ 2 files changed, 302 insertions(+) create mode 100644 code/pycompact/compact_doc.py create mode 100644 code/pycompact/pycompact.py diff --git a/code/pycompact/compact_doc.py b/code/pycompact/compact_doc.py new file mode 100644 index 0000000..f25cb7c --- /dev/null +++ b/code/pycompact/compact_doc.py @@ -0,0 +1,2 @@ +from compact import compact +print compact.__doc__ diff --git a/code/pycompact/pycompact.py b/code/pycompact/pycompact.py new file mode 100644 index 0000000..1934a80 --- /dev/null +++ b/code/pycompact/pycompact.py @@ -0,0 +1,300 @@ +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] + print xsrc.shape + 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() +