test added to function

py2d/SGSModel.py

@@ -131,14 +131,41 @@ def dleithlocal_method(self):#, Psi_hat, Omega_hat, Kx, Ky, Ksq, Delta):
         eddy_viscosity = eddy_viscosity_leith(Cl, Delta, characteristic_Omega)
 
         # Calculate the PI term for local PI = ∇.(ν_e ∇ω )
-        Grad_Omega_hat_dirx = Kx * np.fft.fft2( eddy_viscosity * np.fft.ifft2(Kx * Omega_hat ) )
-        Grad_Omega_hat_diry = Ky * np.fft.fft2( eddy_viscosity * np.fft.ifft2(Ky * Omega_hat) )
+        Grad_Omega_hat_dirx = Kx*np.fft.fft2( eddy_viscosity * np.fft.ifft2(Kx*Omega_hat) )
+        Grad_Omega_hat_diry = Ky*np.fft.fft2( eddy_viscosity * np.fft.ifft2(Ky*Omega_hat) )
         Grad_Omega_hat = Grad_Omega_hat_dirx + Grad_Omega_hat_diry
 
+        #''' test: difference between local and not
+        c_dynamic_old = coefficient_dleith_PsiOmega(Psi_hat, Omega_hat, characteristic_Omega, Kx, Ky, Ksq, Delta)
+        Cl_old = c_dynamic_old ** (1/3)
+        eddy_viscosity_old = eddy_viscosity_leith(Cl_old, Delta, characteristic_Omega)
+        Grad_Omega_hat_old = eddy_viscosity_old *(Ksq*Omega_hat)
+
+        import matplotlib.pyplot as plt
+        VMIN, VMAX = -2, 2
+        plt.figure(figsize=(10,3))
+        plt.subplot(1,3,1)
+        plt.pcolor(np.fft.ifft2(Grad_Omega_hat_old).real,vmin=VMIN,vmax=VMAX,cmap='bwr');plt.colorbar()
+        plt.subplot(1,3,2)
+        plt.pcolor(np.fft.ifft2(Grad_Omega_hat).real,vmin=VMIN,vmax=VMAX,cmap='bwr');plt.colorbar()
+        plt.subplot(1,3,3)
+        plt.pcolor(np.fft.ifft2(Grad_Omega_hat_old-Grad_Omega_hat).real,vmin=VMIN,vmax=VMAX,cmap='bwr');plt.colorbar()
+        plt.show()
+
+        plt.figure(figsize=(7,3))
+        plt.subplot(1,2,1)
+        plt.title(r'$C_L$')
+        plt.pcolor(Cl,vmin=VMIN,vmax=VMAX,cmap='bwr');plt.colorbar()
+        plt.subplot(1,2,2)
+        plt.title(r'$\nu_e$')
+        plt.pcolor(eddy_viscosity,cmap='gray_r');plt.colorbar()
+        plt.show()
+        stop_test
+        #'''
+
         #PiOmega_hat is instead replaced
         eddy_viscosity = 0
         Cl = 0
-
         self.PiOmega_hat, self.eddy_viscosity, self.C_MODEL = PiOmega_hat, eddy_viscosity, Cl
 
         return PiOmega_hat, eddy_viscosity, Cl