Force clear PETSc solver from memory?

Hi everyone.

I’m running into the same issue that is brought up in these two topics:

However, in my case I need to re-mesh the domain at every time step, since the geometry/topology changes, and hence necessarily create a new solver at each time step (even if it is expensive). Is there really no option to clear the existing "LinearProblem"s from memory?

As mentioned in the topics above, forcing garbage collection or setting “solver = ” at the end of each iteration doesn’t seem to work.

Many thanks,
Alex

2 Likes

I ended up with the following: [1], [2]. You might need to update it for the current DolfinX though.

Does both your topology and geometry change at every time step?

  • If it is the topology that changes, you need re-meshing. Might I ask how you have implemented the re-meshing and re initialization of all the objects?
  • If only the geometry changes, there is no need for re-meshing and re-initialization of the Linear solver.

So if I understand correctly, you create an object outside of the loop that contains the PETSc solver, and the on each loop iteration you reset it and feed it the current matrix and vector?

I see how this sidesteps the problem of creating a new solver on each iteration, but it still doesn’t answer the original question, why do old solvers linger in memory.

PETSc4py documents a “destroy” method but even applying that to the PETSc objects themselves doesn’t solve the problem.

i.e. why does even something like this not work:

import gc
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()

problem.solver.destroy()
problem.A.destroy()
problem.b.destroy()
del(problem)
gc.collect()

The domain consists of a 2D rectangle with a circular cutout moving through from one side to the other. The position of the cutout is given externally, and not determined by the solution on each loop.

So the topology only changes when the cutout touches the outer boundary. I would like to maintain a controlled resolution on the cutout boundary (similar to the example here) which is why I opted for remeshing rather than deforming the mesh (the time it takes is negligible).

I have written separate python functions for the data import (including boundary data), meshing and the fem solution like this

for counter in range(N):
    print(counter)
    data, boundary_coords, cutout_coords = load_data(counter,N,flow_data,cutout_data)
    mesh, ct, ft, markers = meshing(boundary_coords, cutout_coords)
    u, p, W = femsolve(mesh,markers,data,ft,ct) #includes setting BCs using 'data' and 'ft'

Without a minimal example reproducing the error is it hard to give any more pointer to what could be wrong. Note that the following minimal example works with no error thrown:

from tqdm import tqdm
import dolfinx
import dolfinx.io
import ufl
from mpi4py import MPI
from petsc4py import PETSc


def solve():
    mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, Nx, Nx)
    V = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 1))

    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    uh = dolfinx.fem.Function(V)

    def g_expr(x):
        return (x[1], 3 * x[0])

    g = dolfinx.fem.Function(V)
    g.interpolate(g_expr)

    a = ufl.inner(u, v) * ufl.dx
    L = ufl.inner(g, v) * ufl.dx
    problem = dolfinx.fem.petsc.LinearProblem(a, L, [], uh)
    problem.solve()


N = 10000
Nx = 50
for i in tqdm(range(N)):
    solve()

running with dolfinx based on: this commit.

Hi dokken,

First of all thank you for looking into this and sorry for not making things clear.

I figured out in the mean time that the issue seems to be related to the options passed to the PETSc in the LinearProblem class. Taking your MWE and replacing

problem = dolfinx.fem.petsc.LinearProblem(a, L, [], uh)

with

problem = dolfinx.fem.petsc.LinearProblem(a, L, [], uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})

the code crashes eventually:

  2%|▏         | 214/10000 [00:08<06:09, 26.52it/s]
---------------------------------------------------------------------------
Error                                     Traceback (most recent call last)
Input In [3], in <cell line: 31>()
     30 Nx = 50
     31 for i in tqdm(range(N)):
---> 32     solve()

Input In [3], in solve()
     23 a = ufl.inner(u, v) * ufl.dx
     24 L = ufl.inner(g, v) * ufl.dx
---> 25 problem = dolfinx.fem.petsc.LinearProblem(a, L, [], uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
     26 problem.solve()

File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/petsc.py:589, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_params, jit_params)
    587 opts.prefixPush(problem_prefix)
    588 for k, v in petsc_options.items():
--> 589     opts[k] = v
    590 opts.prefixPop()
    591 self._solver.setFromOptions()

File PETSc/Options.pyx:23, in petsc4py.PETSc.Options.__setitem__()

File PETSc/Options.pyx:91, in petsc4py.PETSc.Options.setValue()

Error: error code 55
[0] PetscOptionsSetValue() at /usr/local/petsc/src/sys/objects/options.c:1165
[0] PetscOptionsSetValue_Private() at /usr/local/petsc/src/sys/objects/options.c:1215
[0] Out of memory. Allocated: 0, Used by process: 154685440
[0] Number of options 512 < max number of options 512, can not allocate enough space

By copy-pasting the source code of LinearProblem() and inserting

...
588 for k, v in petsc_options.items():
589     opts[k] = v
590 print(opts.getAll())
591 opts.prefixPop()

I can confirm that the options somehow accumulate over iterations, and eventually exceed 512 which seems to be the limit on the number of options that PETSc allows.

Having identified this bug I can figure out a workaround for myself, but perhaps this is worth fixing in a future update?

All the best,
Alex

The problem is a pure petsc issue:
MWE:

from petsc4py import PETSc
import gc
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps",
                 "petsc_options_left": None}


def add_options(prefix, petsc_options):
    opts = PETSc.Options()
    opts.prefixPush(prefix)
    for k, v in petsc_options.items():
        opts[k] = v
    opts.prefixPop()
    opts.clear()
    opts.destroy()
    del opts
    gc.collect()


for i in range(500):
    add_options(f"opts_{i}", petsc_options)

There seems to be no way of destroying or clearing the PETSc options once they have been initiated.

EDIT: I’ve made an issue at: PETSc.Options not cleared or destroyed in Python (#1201) · Issues · PETSc / petsc · GitLab

1 Like