Fixing all yaksa leaks in demos

demo/NavierStokesDrivenCavity.py

@@ -135,6 +135,8 @@ 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()
     sys.exit(1)
 
 # Postprocessing

demo/Stokes.py

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

demo/dirichlet_dirichlet_poisson2D.py

@@ -62,6 +62,9 @@ def main(N, family):
     error = np.sqrt(inner(1, (uj-uq)**2))
     print(f'dirichlet_dirichlet_poisson2D {family:14s} L2 error {error:2.6e}')
 
+    if family == 'legendre':
+        H.destroy()
+
     if 'pytest' not in os.environ:
         import matplotlib.pyplot as plt
         plt.figure()

demo/hermite_poisson2D.py

@@ -79,4 +79,4 @@
         plt.plot((np.squeeze(X[0])[0], np.squeeze(X[0])[-1]), (y, y), 'k')
 
     #plt.show()
-    T.destroy()
+T.destroy()

demo/sphere_helmholtz.py

@@ -80,3 +80,5 @@
     mlab.show()
 else:
     assert error < 1e-6
+
+T.destroy()

shenfun/forms/inner.py

@@ -507,8 +507,7 @@ def inner(expr0, expr1, output_array=None, assemble=None, kind=None, fixed_resol
         output_array += wh
         wh.fill(0)
     if orthogonal:
-        for trans in orthogonal.transfer:
-            trans.destroy()
+        orthogonal.destroy_transfer()
     return output_array
 
 work = CachedArrayDict()

shenfun/legendre/la.py

@@ -15,6 +15,8 @@ class Helmholtz_2dirichlet:
 
     Somewhat experimental.
 
+    The solver needs to be destroyed after use is finished.
+
     """
 
     def __init__(self, matrices):
@@ -83,6 +85,9 @@ def solve_eigen_problem(self, A, B, solver):
             a[2, :-2] = B[2]
             self.lmbda[s], self.V[s, s] = scipy_la.eig_banded(a, lower=True)
 
+    def destroy(self):
+        self.transAB.destroy() 
+
     def __call__(self, b, u, solver=1):
         from shenfun import FunctionSpace, TensorProductSpace, TPMatrix, Identity, la
 
@@ -155,5 +160,6 @@ def __call__(self, b, u, solver=1):
             self.u_B = Helmy(self.rhs_B, self.u_B, fast=False)
             self.transAB.backward(self.u_B, u)
             u[:] = self.V.dot(u)
+            FT.destroy()
 
         return u

shenfun/tensorproductspace.py

@@ -1323,6 +1323,14 @@ def shape(self, forward_output=True):
         if forward_output is not True:
             return self.forward.input_pencil.subshape
         return self.forward.output_array.shape
+
+    def destroy_transfer(self):
+        """Destroy all MPI communicators used for transfer objects of current space
+        
+        Does not destroy the main communicator.
+        """
+        for trans in self.transfer:
+            trans.destroy()
 
     def _get_ndiag_cum_dofs(self):
         """Return the cumulative sum of degrees of freedom along nondiagonal axes"""

shenfun/utilities/__init__.py

@@ -396,6 +396,10 @@ def measure(g):
         else:
             raise NotImplementedError
 
+    if Top.is_padded:
+        Vup.destroy_transfer()
+        Vvp.destroy_transfer() 
+
     if forward_output is True:
 
         if isinstance(output_array, Function):
@@ -404,6 +408,7 @@ def measure(g):
 
         uv_hat = Function(To)
         uv_hat = uv.forward(uv_hat)
+        Top.destroy_transfer()
         return uv_hat
 
     return uv

tests/test_curl.py

@@ -169,6 +169,7 @@ def test_curl_cc():
     d = curl(grad(u_hat))
     assert d.terms() == [[],[],[]]
     T.destroy()
+    T1.destroy()
 
 if __name__ == '__main__':
     #test_curl('d')

tests/test_curvilinear.py

@@ -36,6 +36,7 @@ def test_cylinder():
     u = TrialFunction(V)
     du = div(grad(u))
     assert du.tolatex() == '\\left( \\frac{\\partial^2 u^{x}}{\\partial x^2 }+\\frac{1}{x}\\frac{\\partial  u^{x}}{\\partial x  }+\\frac{1}{x^{2}}\\frac{\\partial^2 u^{x}}{\\partial y^2 }- \\frac{2}{x}\\frac{\\partial  u^{y}}{\\partial y  }- \\frac{1}{x^{2}}u^{x}+\\frac{\\partial^2 u^{x}}{\\partial z^2 }\\right) \\mathbf{b}_{x} \\\\+\\left( \\frac{\\partial^2 u^{y}}{\\partial x^2 }+\\frac{3}{x}\\frac{\\partial  u^{y}}{\\partial x  }+\\frac{2}{x^{3}}\\frac{\\partial  u^{x}}{\\partial y  }+\\frac{1}{x^{2}}\\frac{\\partial^2 u^{y}}{\\partial y^2 }+\\frac{\\partial^2 u^{y}}{\\partial z^2 }\\right) \\mathbf{b}_{y} \\\\+\\left( \\frac{\\partial^2 u^{z}}{\\partial x^2 }+\\frac{1}{x}\\frac{\\partial  u^{z}}{\\partial x  }+\\frac{1}{x^{2}}\\frac{\\partial^2 u^{z}}{\\partial y^2 }+\\frac{\\partial^2 u^{z}}{\\partial z^2 }\\right) \\mathbf{b}_{z} \\\\'
+    T.destroy()
 
 @pytest.mark.parametrize('space', ('cylinder', 'sphere'))
 def test_vector_laplace(space):
@@ -62,6 +63,7 @@ def test_vector_laplace(space):
     b0 = a0.matvec(u_hat, b0)
     b1 = a1.matvec(u_hat, b1)
     assert np.linalg.norm(b0-b1) < 1e-8
+    T.destroy()
 
 if __name__ == '__main__':
     #test_cylinder()

tests/test_dot.py

@@ -15,17 +15,22 @@ def test_dot():
     u = Function(V, buffer=(x, -y))
     bu = dot(u, u)
     assert np.linalg.norm(bu-Function(T, buffer=x**2+y**2)) < 1e-8
+    bu.function_space().destroy_transfer()
     gradu = project(grad(u), S)
     du = dot(gradu, u, forward_output=True)
     xy = Function(V, buffer=(x, y))
     assert np.linalg.norm(du - xy) < 1e-12
+    du.function_space().destroy_transfer() 
     du = dot(u, gradu, forward_output=True)
     assert np.linalg.norm(du - xy) < 1e-12
+    du.function_space().destroy_transfer() 
     gu = dot(gradu, gradu, forward_output=True)
     dd = Function(S, buffer=(1, 0, 0, 1))
     assert np.linalg.norm(gu - dd) < 1e-12
+    gu.function_space().destroy()
     T.destroy()
-
+
+
 #@pytest.mark.skip('skipping')
 def test_dot_curvilinear():
     # Define spherical coordinates without the poles
@@ -45,17 +50,21 @@ def test_dot_curvilinear():
     gij = T.coors.get_metric_tensor(config['basisvectors'])
     ue = Function(T, buffer=sp.sin(theta)**2*gij[1, 1]+sp.cos(phi)**2*gij[2, 2])
     assert np.linalg.norm(ue-u2) < 1e-6
+    u2.function_space().destroy_transfer()
     u = Function(V, buffer=(0, r, 0))
     gu = Function(S, buffer=(0, -1, 0, 1, 0, 0, 0, 0, 1/sp.tan(theta)))
     bu = dot(gu, u)
     but = Function(V, buffer=(-r, 0, 0))
     assert np.linalg.norm(bu-but) < 1e-6
+    bu.function_space().destroy_transfer()
     ub = dot(u, gu)
     ubt = Function(V, buffer=(r, 0, 0))
     assert np.linalg.norm(ub-ubt) < 1e-6
+    ub.function_space().destroy_transfer()
     gg = dot(gu, gu)
     ggt = Function(S, buffer=(-1, 0, 0, 0, -1, 0, 0, 0, 1/sp.tan(theta)**2))
     assert np.linalg.norm(gg-ggt) < 1e-4
+    gg.function_space().destroy_transfer()
     T.destroy()
     config['basisvectors'] = basisvectors
 

tests/test_forms.py

@@ -96,6 +96,7 @@ def test_neg(basis):
             s += sai+sbi
     assert sp.simplify(s) == 0
 
+T.destroy()
 K0 = shenfun.FunctionSpace(N, 'F', dtype='D')
 K1 = shenfun.FunctionSpace(N, 'F', dtype='D')
 K2 = shenfun.FunctionSpace(N, 'F', dtype='d')
@@ -178,6 +179,10 @@ def test_tensor2():
     #assert np.allclose(g.v[1], 0)
     #assert np.allclose(g.v[2], 0)
     #assert np.allclose(g.v[3], 2)
+    T.destroy()
+
+T.destroy()
+C.destroy()
 
 if __name__ == '__main__':
     # test_mul(u2)

tests/test_tensorproductspace.py

@@ -532,7 +532,6 @@ def test_assign(fam):
         VT = VectorSpace(T)
         u_hat = Function(VT)
         u_hat[:, 1:4, 1:4] = 1
-        Tp = T.get_refined((2*N[0], 2*N[1]))
         VTp = VectorSpace(Tp)
         ub_hat = Function(VTp)
         u_hat.assign(ub_hat)
@@ -568,6 +567,7 @@ def test_refine():
     u3 = u_hat.refine(2*np.array(N))
     T.destroy()
     Tp.destroy()
+    Vp.destroy()
 
 def test_eval_expression():
     import sympy as sp

tests/test_transforms.py

@@ -579,6 +579,8 @@ def test_dlt():
     assert np.linalg.norm(c3-1) < 1e-8
     c3 = shenfun.legendre.dlt.leg2cheb(cl, c3, axis=1, transpose=False)
     assert np.linalg.norm(c3(xx)-cl(xx)) < 1e-8
+    TL.destroy()
+    TT.destroy()
     TL = shenfun.TensorProductSpace(shenfun.comm, (F1, L, F), dtype='d')
     TT = shenfun.TensorProductSpace(shenfun.comm, (F1, T, F), dtype='d')
     cl = shenfun.Function(TL, val=1)
@@ -592,6 +594,8 @@ def test_dlt():
     cb1 = cl.backward(kind={'legendre': 'recursive'})
     cb2 = cl.backward(kind={'legendre': 'fast'})
     assert np.linalg.norm(cb1-cb2) < 1e-7, np.linalg.norm(cb1-cb2)
+    TL.destroy()
+    TT.destroy()
 
 def test_leg2cheb():
     for N in (100, 701, 1200):