566 lines
14 KiB
Python
566 lines
14 KiB
Python
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()
|
|
|