472 lines
10 KiB
Python
472 lines
10 KiB
Python
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.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):
|
|
import pprint
|
|
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 (np.linalg.norm((dx - dx_fortran) / dx_fortran))
|
|
print (((dx - dx_fortran) / dx_fortran).min())
|
|
print (((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 (np.linalg.norm((dx - exact) / exact))
|
|
|
|
print (((dx - exact) / exact).min())
|
|
print (((dx - exact) / exact).max())
|
|
|
|
print ("dx_fortran - exact")
|
|
print (np.linalg.norm((dx_fortran - exact) / exact))
|
|
|
|
print (((dx_fortran - exact) / exact).min())
|
|
print (((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 (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 (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):
|
|
|
|
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):
|
|
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():
|
|
|
|
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 (cos_fortran.min(), cos_fortran.max())
|
|
|
|
print (cos_exact.min(), cos_exact.max())
|
|
|
|
print (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)[:]
|
|
true[:] = (1.1 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:]
|
|
dY1 = cs.ddx(Y1)[:]
|
|
|
|
print (dY1.min(), dY1.max())
|
|
|
|
print (true.min(), true.max())
|
|
|
|
relerr = (dY1 - true) / true
|
|
|
|
print (np.nanmin(relerr), np.nanmax(relerr))
|
|
|
|
|
|
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[:] = true[:]
|
|
|
|
dydx = np.memmap("dphi", dtype=np.float64, mode="w+", shape=cs.shape)
|
|
dydx[:] = dY1[:]
|
|
|
|
|
|
# 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()
|
|
|