305 lines
6.8 KiB
Python
305 lines
6.8 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.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: # 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.)
|
|
|
|
print ("calculated pi value = ", pi8)
|
|
|
|
l_0 = 2.0
|
|
hyp=l_0*pi8/ny
|
|
hxp=hyp
|
|
hzp=hyp
|
|
|
|
|
|
cs = CompactScheme(nx, ny, nz, False, True, True, 4., 2., 2.)
|
|
|
|
|
|
Y1 = 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))[:]
|
|
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()
|
|
|