diff --git a/pycompact.py b/pycompact.py index d56d746..e0d61c9 100644 --- a/pycompact.py +++ b/pycompact.py @@ -7,6 +7,8 @@ class CompactScheme: pi8 = np.arccos(-1., dtype=np.float64) + self.nb = 32 + self.shape = (nz, ny, nx) self.px = px @@ -195,34 +197,36 @@ class CompactScheme: return coef - def ddx (self, src): + def ddx (self, src, dst=None): - if src.shape != self.shape: - print ("error") + if src.shape != self.shape: print ("error: src.shape wrong!") nz, ny, nx = self.shape + nb = self.nb - dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + if dst is None or dst.shape != self.shape : + 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) + for j in range(0, ny, nb): + dst[i,j:j+nb] = compact.dfp(self.hx, src[i,j:j+nb], 1) else: for i in range(nz): - dst[i] = compact.dfnonp(self.hx, src[i], 1) + for j in range(0, ny, nb): + dst[i,j:j+nb] = compact.dfnonp(self.hx, src[i,j:j+nb], 1) - # return np.swapaxes(dst, 1, 2) return dst - def ddy (self, src): + def ddy (self, src, dst=None): - if src.shape != self.shape: - print ("error") + if src.shape != self.shape: print ("error: src.shape wrong!") nz, ny, nx = self.shape - dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + if dst is None or dst.shape != self.shape : + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) if self.py: # Periodic BC for i in range(nz): @@ -232,53 +236,60 @@ class CompactScheme: 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): + def ddz (self, src, dst=None): - if src.shape != self.shape: - print ("error") + if src.shape != self.shape: print ("error: src.shape wrong!") nz, ny, nx = self.shape - # dst = np.zeros((nx, ny, nz,), order="F") - dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + if dst is None or dst.shape != self.shape : + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + zxsrc = np.zeros((nz, nx)) if self.pz: # Periodic BC for i in range(ny): - dst[:,i,:] = compact.dfp(self.hx, src[:,i,:].T, 3).T + zxsrc[:,:] = src[:,i,:] + dst[:,i,:] = compact.dfp(self.hx, zxsrc.T, 3).T else: for i in range(ny): - dst[:,i,:] = compact.dfnonp(self.hx, src[:,i,:].T, 3).T + zxsrc[:,:] = src[:,i,:] + dst[:,i,:] = compact.dfnonp(self.hx, zxsrc.T, 3).T - # return np.swapaxes(dst, 1, 2) return dst - def d2dx (self, src): + def d2dx (self, src, dst=None): - if src.shape != self.shape: print ("error") + if src.shape != self.shape: print ("error: src.shape wrong!") nz, ny, nx = self.shape + nb = self.nb - dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + if dst is None or dst.shape != self.shape : + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) if self.px: # Periodic BC for i in range(nz): - dst[i] = compact.d2fp(self.hx, src[i], 1) + for j in range(0, ny, nb): + dst[i,j:j+nb] = compact.d2fp(self.hx, src[i,j:j+nb], 1) + else: for i in range(nz): - dst[i] = compact.d2fnonp(self.hx, src[i], 1) + for j in range(0, ny, nb): + dst[i,j:j+nb] = compact.d2fnonp(self.hx, src[i,j:j+nb], 1) return dst - def d2dy (self, src): + def d2dy (self, src, dst=None): - if src.shape != self.shape: print ("error") + if src.shape != self.shape: print ("error: src.shape wrong!") nz, ny, nx = self.shape - dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + if dst is None or dst.shape != self.shape : + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) if self.py: # Periodic BC for i in range(nz): @@ -289,20 +300,25 @@ class CompactScheme: return dst - def d2dz (self, src): + def d2dz (self, src, dst=None): - if src.shape != self.shape: print ("error") + if src.shape != self.shape: print ("error: src.shape wrong!") nz, ny, nx = self.shape - dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + if dst is None or dst.shape != self.shape : + dst = np.zeros((nz, ny, nx,), dtype=np.float64,) + + zxsrc = np.zeros((nz, nx)) if self.pz: # Periodic BC for i in range(ny): - dst[:,i,:] = compact.d2fp(self.hx, src[:,i,:].T, 3).T + zxsrc[:,:] = src[:,i,:] + dst[:,i,:] = compact.d2fp(self.hx, zxsrc.T, 3).T else: for i in range(ny): - dst[:,i,:] = compact.d2fnonp(self.hx, src[:,i,:].T, 3).T + zxsrc[:,:] = src[:,i,:] + dst[:,i,:] = compact.d2fnonp(self.hx, zxsrc.T, 3).T return dst @@ -513,7 +529,7 @@ def validate_trigonometric(): print(" DDX Test ") 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)[:] + dY1 = cs.ddx(Y1, dY1)[:] compare_3d_result(true, dY1) @@ -521,7 +537,7 @@ def validate_trigonometric(): print("D2DX Test ") Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] true[:] = (- 1.1 * 1.1 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] - dY1 = cs.d2dx(Y1)[:] + dY1 = cs.d2dx(Y1, dY1)[:] compare_3d_result(true, dY1) @@ -529,7 +545,7 @@ def validate_trigonometric(): print(" DDY Test ") Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] true[:] = (3.0 * np.sin(1.1 * xx) * np.cos(3.0 * yy) * np.sin(2.0 * zz))[:] - dY1 = cs.ddy(Y1)[:] + dY1 = cs.ddy(Y1, dY1)[:] compare_3d_result(true, dY1) @@ -537,7 +553,7 @@ def validate_trigonometric(): print(" DDZ Test ") Y1[:] = np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz)[:] true[:] = (2.0 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.cos(2.0 * zz))[:] - dY1 = cs.ddz(Y1)[:] + dY1 = cs.ddz(Y1, dY1)[:] compare_3d_result(true, dY1) @@ -576,6 +592,7 @@ def performance_benchmark(): print ("Test ddx") Y1 = np.zeros(shape) + dY1 = np.zeros(shape) true = np.zeros(shape) @@ -605,7 +622,7 @@ def performance_benchmark(): true[:] = (1.1 * np.cos(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] t0 = time.perf_counter() for i in range(100): - dY1 = cs.ddx(Y1)[:] + cs.ddx(Y1, dY1) t1 = time.perf_counter() print ("Elapsed time in sec for 100 repetition: ", t1-t0) compare_3d_result(true, dY1) @@ -614,7 +631,7 @@ def performance_benchmark(): true[:] = (- 1.1 * 1.1 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] t0 = time.perf_counter() for i in range(100): - dY1 = cs.d2dx(Y1)[:] + cs.d2dx(Y1, dY1) t1 = time.perf_counter() print ("Elapsed time in sec for 100 repetition: ", t1-t0) compare_3d_result(true, dY1) @@ -623,7 +640,7 @@ def performance_benchmark(): true[:] = (3.0 * np.sin(1.1 * xx) * np.cos(3.0 * yy) * np.sin(2.0 * zz))[:] t0 = time.perf_counter() for i in range(100): - dY1 = cs.ddy(Y1)[:] + cs.ddy(Y1, dY1) t1 = time.perf_counter() print ("Elapsed time in sec for 100 repetition: ", t1-t0) compare_3d_result(true, dY1) @@ -632,7 +649,7 @@ def performance_benchmark(): true[:] = (-3. * 3.0 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] t0 = time.perf_counter() for i in range(100): - dY1 = cs.d2dy(Y1)[:] + cs.d2dy(Y1, dY1) t1 = time.perf_counter() print ("Elapsed time in sec for 100 repetition: ", t1-t0) compare_3d_result(true, dY1) @@ -641,7 +658,7 @@ def performance_benchmark(): true[:] = (2.0 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.cos(2.0 * zz))[:] t0 = time.perf_counter() for i in range(100): - dY1 = cs.ddz(Y1)[:] + cs.ddz(Y1, dY1) t1 = time.perf_counter() print ("Elapsed time in sec for 100 repetition: ", t1-t0) compare_3d_result(true, dY1) @@ -651,7 +668,7 @@ def performance_benchmark(): true[:] = (-2. * 2.0 * np.sin(1.1 * xx) * np.sin(3.0 * yy) * np.sin(2.0 * zz))[:] t0 = time.perf_counter() for i in range(100): - dY1 = cs.d2dz(Y1)[:] + cs.d2dz(Y1, dY1) t1 = time.perf_counter() print ("Elapsed time in sec for 100 repetition: ", t1-t0) compare_3d_result(true, dY1)