Making all demos cleaned up for communicators

demo/ChannelFlow.py

@@ -1,6 +1,7 @@
+import warnings
 from warnings import WarningMessage
 from shenfun import *
-np.warnings.filterwarnings('ignore')
+warnings.filterwarnings('ignore')
 
 class KMM:
     """Navier Stokes channel flow solver

demo/Ginzburg_Landau.py

@@ -13,7 +13,8 @@
 import matplotlib.pyplot as plt
 from mpi4py_fft import generate_xdmf, fftw
 from shenfun import div, grad, TestFunction, TrialFunction, \
-    TensorProductSpace, Array, Function, ETDRK4, HDF5File, FunctionSpace, comm
+    TensorProductSpace, Array, Function, ETDRK4, HDF5File, FunctionSpace, \
+    comm, cleanup
 
 # Use sympy to set up initial condition
 x, y = symbols("x,y", real=True)
@@ -86,6 +87,4 @@ def update(self, u, u_hat, t, tstep, plot_tstep, write_tstep, file, **params):
     if comm.Get_rank() == 0:
         generate_xdmf("Ginzburg_Landau_{}.h5".format(N[0]))
     fftw.export_wisdom('GL.wisdom')
-    T.destroy()
-    for space in U_hat._padded_space.values():
-        space.destroy()
+    cleanup(vars())

demo/Ginzburg_Landau_sphere_IMEX.py

@@ -61,14 +61,7 @@ def __init__(self, N=(32, 64), dt=0.25, refineplot=False,
         phij = np.hstack([phij, phij[:, 0][:, None]])
         self.file_u = HDF5File(filename+'_U.h5', domain=[np.squeeze(d) for d in [thetaj, phij]], mode='a')
         self.file_c = HDF5File(filename+'_C.h5', domain=self.T.mesh(), mode='w')
-
-    def destroy(self):
-        spaces = [self.T, self.T2, self.T3, *self.u_hat._padded_space.values()]
-        if self.refineplot not in (False, 1):
-            spaces.append(self.ur_hat.function_space())
-        for T in spaces:
-            T.destroy()
-
+
     def prepare_step(self, rk=0):
         up = self.u_hat.backward(padding_factor=(1.5, 1.5))
         self.H_ = up.function_space().forward(up*abs(up)**2, self.H_)
@@ -176,4 +169,4 @@ def solve(self, t=0, tstep=0, end_time=100):
     print(f'Computing time {time()-t0: 2.4f}')
     if comm.Get_rank() == 0:
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5')
-    c.destroy()
+    cleanup(vars(c))

demo/Gray_Scott_fractal.py

@@ -35,7 +35,7 @@
 from mpi4py_fft import generate_xdmf
 from shenfun import inner, div, grad, TestFunction, TrialFunction, Function, \
     HDF5File, ETDRK4, TensorProductSpace, VectorSpace, FunctionSpace, Array, \
-    comm, get_simplified_tpmatrices
+    comm, get_simplified_tpmatrices, cleanup
 
 # Use sympy to set up initial condition
 x, y = symbols("x,y", real=True)
@@ -126,5 +126,5 @@ def update(self, uv, uv_hat, t, tstep, **params):
     integrator.setup(dt)
     UV_hat = integrator.solve(UV, UV_hat, dt, (0, end_time))
     generate_xdmf("Gray_Scott_{}.h5".format(N[0]))
-    for space in (T, TVp, UV_hat._padded_space[padding_factor]):
-        space.destroy()
+    cleanup(vars())
+

demo/KleinGordon.py

@@ -156,6 +156,4 @@ def update(self, fu, fu_hat, t, tstep, **params):
 
     if rank == 0:
         generate_xdmf("KleinGordon{}.h5".format(N[0]))
-
-    for space in (T, Tp):
-        space.destroy()
+    cleanup((T, Tp))

demo/Kuramato_Sivashinsky.py

@@ -90,5 +90,4 @@ def update(self, u, u_hat, t, tstep, plot_step, **params):
     #integrator = RK4(T, L=LinearRHS, N=NonlinearRHS, update=update, **par)
     integrator.setup(dt)
     U_hat = integrator.solve(U, U_hat, dt, (0, end_time))
-    for space in (T, Tp):
-        space.destroy()
+    cleanup((T, Tp))

demo/MKM.py

@@ -77,17 +77,14 @@ def plot(self, t, tstep):
             ub = self.u_.backward(self.ub)
             if comm.Get_rank() == 0:
                 X = self.X
-                self.im1.axes.clear()
                 self.im1.axes.contourf(X[1][:, :, 0], X[0][:, :, 0], ub[0, :, :, 0], 100)
                 self.im1.autoscale()
                 plt.figure(1)
                 plt.pause(1e-6)
-                self.im2.axes.clear()
                 self.im2.axes.contourf(X[1][:, :, 0], X[0][:, :, 0], ub[1, :, :, 0], 100)
                 self.im2.autoscale()
                 plt.figure(2)
                 plt.pause(1e-6)
-                self.im3.axes.clear()
                 self.im3.axes.contourf(X[2][:, 0, :], X[0][:, 0, :], ub[0, :, 0, :], 100)
                 self.im3.autoscale()
                 plt.figure(3)
@@ -129,8 +126,8 @@ def update(self, t, tstep):
         'Re': 180.,
         'dt': 0.0005,
         'filename': f'MKM_{N[0]}_{N[1]}_{N[1]}',
-        'conv': 1,
-        'modplot': -100,
+        'conv': 0,
+        'modplot': 100,
         'modsave': 1000,
         'moderror': 1,
         'family': 'C',
@@ -141,7 +138,8 @@ def update(self, t, tstep):
     c = MKM(**d)
     t, tstep = c.initialize(from_checkpoint=False)
     t0 = time()
-    c.solve(t=t, tstep=tstep, end_time=40)
+    c.solve(t=t, tstep=tstep, end_time=0.01)
     print('Computing time %2.4f'%(time()-t0))
     if comm.Get_rank() == 0:
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5')
+    cleanup(vars(c))

demo/MKM_MicroPolar.py

@@ -101,17 +101,14 @@ def plot(self, t, tstep):
             ub = self.u_.backward(self.ub)
             if comm.Get_rank() == 0:
                 X = self.X
-                self.im1.axes.clear()
                 self.im1.axes.contourf(X[1][:, :, 0], X[0][:, :, 0], ub[0, :, :, 0], 100)
                 self.im1.autoscale()
                 plt.figure(1)
                 plt.pause(1e-6)
-                self.im2.axes.clear()
                 self.im2.axes.contourf(X[1][:, :, 0], X[0][:, :, 0], ub[1, :, :, 0], 100)
                 self.im2.autoscale()
                 plt.figure(2)
                 plt.pause(1e-6)
-                self.im3.axes.clear()
                 self.im3.axes.contourf(X[2][:, 0, :], X[0][:, 0, :], ub[0, :, 0, :], 100)
                 self.im3.autoscale()
                 plt.figure(3)
@@ -365,16 +362,16 @@ def fromfile(self, filename="stats"):
     from time import time
     from mpi4py_fft import generate_xdmf
     t0 = time()
-    N = (128, 128, 64)
+    N = (64, 64, 32)
     d = {
         'N': N,
         'Re': 180.,
         'dt': 0.0005,
         'filename': f'MKM_MP_{N[0]}_{N[1]}_{N[2]}',
         'conv': 1,
-        'modplot': 100,
+        'modplot': -10,
         'modsave': 1000,
-        'moderror': 100,
+        'moderror': 1,
         'family': 'C',
         'checkpoint': 100,
         'sample_stats': 100,
@@ -383,8 +380,8 @@ def fromfile(self, filename="stats"):
         'timestepper': 'IMEXRK222', # IMEXRK222, IMEXRK443
         }
     c = MKM(**d)
-    t, tstep = c.initialize(from_checkpoint=True)
-    c.solve(t=t, tstep=tstep, end_time=50)
+    t, tstep = c.initialize(from_checkpoint=False)
+    c.solve(t=t, tstep=tstep, end_time=0.001)
     print('Computing time %2.4f'%(time()-t0))
     if comm.Get_rank() == 0:
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5')
@@ -401,3 +398,4 @@ def fromfile(self, filename="stats"):
         # multiply by utau=0.0635 to get the same scaling as fig 7
         plt.semilogx((x[32:] + 1)*180, np.sqrt(ww[1, 32:])*0.0635, 'b', (1 - x[:32])*180, np.sqrt(ww[1, :32])*0.0635, 'r')
         #plt.show()
+    cleanup(vars(c))

demo/MicroPolar.py

@@ -1,7 +1,8 @@
 from random import sample
 from shenfun import *
 from ChannelFlow import KMM
-np.warnings.filterwarnings('ignore')
+import warnings
+warnings.filterwarnings('ignore')
 
 class MicroPolar(KMM):
     """Micropolar channel flow solver
@@ -174,19 +175,19 @@ def convection(self):
         KMM.convection(self)
         HW = self.HW_
         up = self.up
-        dw0dxp = self.dw0dx().backward(padding_factor=self.padding_factor)
-        dw0dyp = self.dw0dy().backward(padding_factor=self.padding_factor)
-        dw0dzp = self.dw0dz().backward(padding_factor=self.padding_factor)
-        dw1dxp = self.dw1dx().backward(padding_factor=self.padding_factor)
-        dw1dyp = self.dw1dy().backward(padding_factor=self.padding_factor)
-        dw1dzp = self.dw1dz().backward(padding_factor=self.padding_factor)
-        dw2dxp = self.dw2dx().backward(padding_factor=self.padding_factor)
-        dw2dyp = self.dw2dy().backward(padding_factor=self.padding_factor)
-        dw2dzp = self.dw2dz().backward(padding_factor=self.padding_factor)
-        HW[0] = self.TDp.forward(up[0]*dw0dxp+up[1]*dw0dyp+up[2]*dw0dzp, HW[0])
-        HW[1] = self.TDp.forward(up[0]*dw1dxp+up[1]*dw1dyp+up[2]*dw1dzp, HW[1])
-        HW[2] = self.TDp.forward(up[0]*dw2dxp+up[1]*dw2dyp+up[2]*dw2dzp, HW[2])
-        HW.mask_nyquist(self.mask)
+        dw0dxp = self.dw0dxp = self.dw0dx().backward(padding_factor=self.padding_factor).v
+        #dw0dyp = self.dw0dy().backward(padding_factor=self.padding_factor)
+        #dw0dzp = self.dw0dz().backward(padding_factor=self.padding_factor)
+        #dw1dxp = self.dw1dx().backward(padding_factor=self.padding_factor)
+        #dw1dyp = self.dw1dy().backward(padding_factor=self.padding_factor)
+        #dw1dzp = self.dw1dz().backward(padding_factor=self.padding_factor)
+        #dw2dxp = self.dw2dx().backward(padding_factor=self.padding_factor)
+        #dw2dyp = self.dw2dy().backward(padding_factor=self.padding_factor)
+        #dw2dzp = self.dw2dz().backward(padding_factor=self.padding_factor)
+        #HW[0] = self.TDp.forward(up[0]*dw0dxp+up[1]*dw0dyp+up[2]*dw0dzp, HW[0])
+        #HW[1] = self.TDp.forward(up[0]*dw1dxp+up[1]*dw1dyp+up[2]*dw1dzp, HW[1])
+        #HW[2] = self.TDp.forward(up[0]*dw2dxp+up[1]*dw2dyp+up[2]*dw2dzp, HW[2])
+        #HW.mask_nyquist(self.mask)
 
     def tofile(self, tstep):
         self.file_u.write(tstep, {'u': [self.u_.backward(mesh='uniform')]}, as_scalar=True)

demo/NavierStokesDrivenCavity.py

@@ -135,8 +135,7 @@ def compute_rhs(up_hat, bh_hat):
 u_, p_ = up
 
 if 'pytest' in os.environ:
-    for T in [V0, V1, P, uiuj.function_space(), up_hat._padded_space[1.5]]:
-        T.destroy()
+    cleanup(vars())
     sys.exit(1)
 
 # Postprocessing
@@ -181,6 +180,4 @@ def compute_rhs(up_hat, bh_hat):
     print("%d %d " %(xi, yi) +("%+2.7e "*4) %(xmid, ymid, psi_min, err))
     count += 1
 
-# clean up
-for T in [V0, V1, P, P1, uiuj.function_space(), up_hat._padded_space[1.5]]:
-    T.destroy()
+cleanup(vars())

demo/NavierStokesDrivenCavityIRK3.py

@@ -72,13 +72,6 @@ def __init__(self, N=(10, 10), dt=0.01, Re=100., modplot=10, family='C'):
         self.uiuj_hat = Function(S1)
         IRK3.__init__(self, VQ)
 
-    def destroy(self):
-        for T in [*self.up_.function_space(),
-                  *self.uip.function_space(),
-                  *self.uiuj.function_space(),
-                  self.T1]:
-            T.destroy()
-
     def LinearRHS(self, up):
         u, p = up
         return self.nu*div(grad(u))-grad(p)
@@ -130,6 +123,7 @@ def update(self, up, up_hat, t, tstep, **par):
             plt.pause(0.01)
 
 if __name__ == '__main__':
+    import sys
     d = {'N': (32, 32),
          'dt': 0.1,
          'modplot': 100,
@@ -138,7 +132,7 @@ def update(self, up, up_hat, t, tstep, **par):
         }
     sol = NavierStokesIRK3(**d)
     sol.solve(sol.up_, sol.up_hat, sol.dt, (0, 100))
-
+    
     # Solve streamfunction
     r = TestFunction(sol.T[0][1])
     s = TrialFunction(sol.T[0][1])
@@ -147,7 +141,7 @@ def update(self, up, up_hat, t, tstep, **par):
     H = la.SolverGeneric2ND(S)
     phi_h = H(h)
     phi = phi_h.backward()
-
+    
     # Find minimal streamfunction value and position
     # by gradually zooming in on mesh
     W = 101
@@ -171,4 +165,4 @@ def update(self, up, up_hat, t, tstep, **par):
         dx = dx/4.
         print("%d %d " %(xi, yi) +("%+2.7e "*4) %(xmid, ymid, psi_min, err))
         count += 1
-    sol.destroy()
+    cleanup(vars(sol))

demo/OrrSommerfeld.py

@@ -133,3 +133,4 @@ def print_energy_and_divergence(self, t, tstep):
         plt.contourf(X[1][0, :, 0], X[0][:, 0, 0], p0)
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5')
     plt.show()
+    cleanup(vars(c))

demo/OrrSommerfeld2D.py

@@ -125,3 +125,4 @@ def print_energy_and_divergence(self, t, tstep):
     print('Computing time %2.4f'%(time()-t0))
     if comm.Get_rank() == 0:
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5')
+    cleanup(vars(c))

demo/RayleighBenard.py

@@ -212,4 +212,5 @@ def solve(self, t=0, tstep=0, end_time=1000):
     t0 = time()
     c.solve(t=t, tstep=tstep, end_time=2)
     print('Computing time %2.4f'%(time()-t0))
+    cleanup(vars(c))
 

demo/RayleighBenard2D.py

@@ -81,7 +81,7 @@ def update_bc(self, t):
         self.T_.get_dealiased_space(self.padding_factor).bases[0].bc.update(t)
 
     def prepare_step(self, rk):
-        self.convection()
+        self.convection(rk)
         Tp = self.T_.backward(padding_factor=self.padding_factor)
         self.uT_ = self.up.function_space().forward(self.up*Tp, self.uT_)
 
@@ -159,7 +159,6 @@ def plot(self, t, tstep):
                 self.im1.scale = np.linalg.norm(ub[1])
                 plt.pause(1e-6)
                 plt.figure(2)
-                self.im2.axes.clear()
                 self.im2.axes.contourf(self.X[1][:, :], self.X[0][:, :], self.Tb[:, :], 100)
                 self.im2.autoscale()
                 plt.pause(1e-6)
@@ -213,4 +212,4 @@ def solve(self, t=0, tstep=0, end_time=1000):
         generate_xdmf('_'.join((d['filename'], 'U'))+'.h5', order='visit')
         generate_xdmf('_'.join((d['filename'], 'T'))+'.h5', order='visit')
     print('Computing time %2.4f'%(time()-t0))
-
+    cleanup(vars(c))

demo/Stokes.py

@@ -110,9 +110,8 @@ def main(N, family):
         plt.show()
     else:
         assert np.all(abs(np.array(error)) < 1e-8), error
-    TD.destroy()
-    Q.destroy()
-
+    cleanup((TD, Q))
+
 if __name__ == '__main__':
     for family in 'CULQJ':
-        main(24, family)
+        main(24, family)

demo/Stokes2NP.py

@@ -131,3 +131,5 @@
     plt.figure()
     plt.contourf(X[0], X[1], u_[0], 100)
     #plt.show()
+
+cleanup(vars())

demo/Stokes3D.py

@@ -114,3 +114,4 @@
     plt.figure()
     plt.contourf(X[0][:, :, 0], X[1][:, :, 0], u_[0, :, :, 6], 100)
     plt.show()
+cleanup(vars())