python wrapper using f2py
This commit is contained in:
ignis
parent
b7a16a467e
commit
22e6fa615d
2 changed files with 302 additions and 0 deletions
2
code/pycompact/compact_doc.py
Normal file
2
code/pycompact/compact_doc.py
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
from compact import compact
|
||||
print compact.__doc__
|
||||
300
code/pycompact/pycompact.py
Normal file
300
code/pycompact/pycompact.py
Normal file
|
|
@ -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()
|
||||
|
||||
Loading…
Add table
Reference in a new issue