From ca3920c8f9339d95ef3aad9d534edcca436397af Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Wed, 3 Mar 2021 05:29:33 +0900
Subject: [PATCH] added ddy ddz. need some fix

---
 code/pycompact/pycompact.py | 85 +++++++++++++++++++++++++++++++++----
 1 file changed, 76 insertions(+), 9 deletions(-)

diff --git a/code/pycompact/pycompact.py b/code/pycompact/pycompact.py
index 73e704a..59e04c6 100644
--- a/code/pycompact/pycompact.py
+++ b/code/pycompact/pycompact.py
@@ -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[:]