research/incomp-flame-post
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

added ddy ddz. need some fix

Browse source
This commit is contained in:
ignis 2021-03-03 05:29:33 +09:00
parent 2cd32e62c3
commit ca3920c8f9
1 changed files with 76 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

85
code/pycompact/pycompact.py
View file

@ -222,11 +222,55 @@ class CompactScheme:
# return np.swapaxes(dst, 1, 2)
return dst
def ddy (self):
return
def ddy (self, src):
def ddz (self):
return
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):
xysrc = src[i].reshape(-1).reshape((nx,ny), order='F')
dst[i] = compact.dfp(self.hx, xysrc, 2).reshape(-1).reshape((ny,nx), order='C')
else:
for i in range(nz):
xysrc = src[i].reshape(-1).reshape((nx,ny), order='F')
dst[i] = compact.dfnonp(self.hx, xysrc, 2).reshape(-1).reshape((ny,nx), order='C')
# return np.swapaxes(dst, 1, 2)
return dst
def ddz (self, src):
if src.shape != self.shape:
print ("error")
nz, ny, nx = self.shape
xsrc = np.zeros((nz, 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.pz: # Periodic BC
for i in range(ny):
xsrc[:] = src[:,i,:]
dst[i] = compact.dfp(self.hx, xsrc, 3)
else:
for i in range(ny):
xsrc[:] = src[:,i,:]
dst[i] = compact.dfnonp(self.hx, xsrc, 3)
# return np.swapaxes(dst, 1, 2)
return dst
def port_nonp_coef (self):
@ -413,16 +457,39 @@ def validate_trigonometric():
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)[:]
print (dY1.min(), dY1.max())
compare_3d_result(true, dY1)
print (true.min(), true.max())
relerr = (dY1 - true) / true
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)[:]
print (np.nanmin(relerr), np.nanmax(relerr))
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:
@ -430,7 +497,7 @@ def validate_trigonometric():
y[:] = Y1[:]
dydxtrue = np.memmap("dphitrue", dtype=np.float64, mode="w+", shape=cs.shape)
dydxtrue[:] = true[:]
dydxtrue[:] = true1[:]
dydx = np.memmap("dphi", dtype=np.float64, mode="w+", shape=cs.shape)
dydx[:] = dY1[:]