Error setting PETSc options repeatedly

Errors
1

After a repeated evaluation of a problem PETSc.Options.setValue throws error

Traceback (most recent call last):
  File "/root/shared/./test.py", line 45, in <module>
    solution = dolfinx.fem.LinearProblem(a_p, L_p,
  File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/problem
.py", line 83, in __init__
    opts[k] = v
  File "PETSc/Options.pyx", line 23, in petsc4py.PETSc.Options.__setitem__
  File "PETSc/Options.pyx", line 91, in petsc4py.PETSc.Options.setValue
petsc4py.PETSc.Error: error code 55

In my case error shows after 475 iterations of this MWE (less for a more involved code):

import dolfinx
import dolfinx.io
from mpi4py import MPI
import ufl
from ufl import inner, grad, dx
import numpy as np

def set_bc(V):
    u_bc = dolfinx.Function(V)

    bc_facets = np.where(
        np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
    bc_dofs = dolfinx.fem.locate_dofs_topological(V,
                                                  mesh.topology.dim-1,
                                                  bc_facets)

    with u_bc.vector.localForm() as loc:
        loc.setValues(bc_dofs, np.full(len(bc_dofs), 0))
    return dolfinx.DirichletBC(u_bc, bc_dofs)


mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 100, 100)
mesh.topology.create_connectivity(mesh.topology.dim-1,
                                  mesh.topology.dim)

Hcurl = ufl.FiniteElement("Nedelec 1st kind H(curl)", mesh.ufl_cell(), 1)
H1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)

N = 1000
for i in range(N):
    print(f"{i+1}/{N}")

    V = dolfinx.FunctionSpace(mesh, H1)
    u, v = ufl.TrialFunctions(V), ufl.TestFunctions(V)
    f = dolfinx.Function(V)
    f.interpolate(lambda x: (i+1) * x[0]**2 + 0*x[1])
    f.x.scatter_forward()

    V = dolfinx.FunctionSpace(mesh, Hcurl)
    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

    a_p = inner((i+1)*u, v) * dx
    L_p = inner(v, grad(f)) * dx

    solution = dolfinx.fem.LinearProblem(a_p, L_p,
                                         petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    u = solution.solve()

Running using the command

docker run -ti -e DISPLAY=$DISPLAY -v /tmp/.X11-unix:/tmp/.X11-unix -v $(pwd):/root/shared -w /root/shared --rm dolfinx/dolfinx python3 ./test.py
2

You don’t need to instantiate all these objects every iteration…

import dolfinx
import dolfinx.io
from mpi4py import MPI
import ufl
from ufl import inner, grad, dx
import numpy as np

def set_bc(V):
    u_bc = dolfinx.Function(V)

    bc_facets = np.where(
        np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
    bc_dofs = dolfinx.fem.locate_dofs_topological(V,
                                                  mesh.topology.dim-1,
                                                  bc_facets)

    with u_bc.vector.localForm() as loc:
        loc.setValues(bc_dofs, np.full(len(bc_dofs), 0))
    return dolfinx.DirichletBC(u_bc, bc_dofs)


mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 4, 4)
mesh.topology.create_connectivity(mesh.topology.dim-1,
                                  mesh.topology.dim)

Hcurl = ufl.FiniteElement("Nedelec 1st kind H(curl)", mesh.ufl_cell(), 1)
H1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)

N = 1000


V = dolfinx.FunctionSpace(mesh, H1)
f = dolfinx.Function(V)

V = dolfinx.FunctionSpace(mesh, Hcurl)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

i_plus_one = dolfinx.Constant(mesh, 1.0)
a_p = inner(i_plus_one*u, v) * dx
L_p = inner(v, grad(f)) * dx

solution = dolfinx.fem.LinearProblem(a_p, L_p,
                                     petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

for i in range(N):
    print(f"{i+1}/{N}")

    i_plus_one.value = float(i + 1)
    f.interpolate(lambda x: (i+1) * x[0]**2 + 0*x[1])
    f.x.scatter_forward()
    u = solution.solve()
3 Likes
3

It is most likely because you initialize a new LinearSolver (and petsc ksp solver) 1000 times. If you are solving the same problem on every time step (i guess your application is slightly different), you should reuse the PETSc solver.

1 Like
4

Is there a way to avoid instantiation of the solution if my forms contain functions from a MixedElement space, which I get from call phi_re, phi_im = ufl.split(solution.solve()) (like in this question) ?
This call returns new objects, so I have to create new forms and instantiate new LinearProblem.

5

BTW, forcing garbage collection for each iteration does not remove the issue.

6

I would strongly suggest to use the petsc4py interface directly, as for instance shown in:
https://jorgensd.github.io/dolfinx-tutorial/chapter2/diffusion_code.html
https://jorgensd.github.io/dolfinx-tutorial/chapter2/ns_code2.html
That way, you can assemble in to the same matrix multiple times, and you can write your code such that you avoid having

in a loop

1 Like
7

Could you, please, elaborate on not needing to ufl.split function in a loop?

V = dolfinx.FunctionSpace(mesh, ufl.MixedElement(H1, H1))
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
foo = dolfinx.Function(V)
...
for i in range(N):
    solver.solve(self._b, self.solution.fun.vector)
    foo.x.scatter_forward()
    v1, v2 = ufl.split(foo)

My solver has multiple steps and one of the steps requires the split solution. How can I manually split the solution and copy the underlying PETSc vectors without allocations in a loop?

8

Here is a minimal example of two sets of projections, where the second one is dependent on the first one.
The key is to define the function and split it before the loop:

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


mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
V0 = V.sub(0).collapse()

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
uh = dolfinx.Function(V)
u0, u1 = ufl.split(uh)


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


# First PDE
g = dolfinx.Function(V)
g.interpolate(g_expr)

a = ufl.inner(u, v) * ufl.dx
L = ufl.inner(g, v) * ufl.dx

# Second PDE
u2 = ufl.TrialFunction(V0)
v2 = ufl.TestFunction(V0)
a2 = u2 * v2 * ufl.dx
L2 = u0 * u1 * v2 * ufl.dx

# Solver for system 1
A = dolfinx.fem.create_matrix(a)
b = dolfinx.fem.create_vector(L)
solver1 = PETSc.KSP().create(mesh.mpi_comm())
solver1.setOperators(A)
solver1.setType(PETSc.KSP.Type.PREONLY)
solver1.getPC().setType(PETSc.PC.Type.LU)


# Solver for system 2
A2 = dolfinx.fem.create_matrix(a2)
b2 = dolfinx.fem.create_vector(L2)
uh2 = dolfinx.Function(V0)
solver2 = PETSc.KSP().create(mesh.mpi_comm())
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.PREONLY)
solver2.getPC().setType(PETSc.PC.Type.LU)

N = 2
for i in range(N):

    # Assemble and solve first system
    dolfinx.fem.assemble_matrix(A, a)
    A.assemble()
    dolfinx.fem.assemble_vector(b, L)
    solver1.solve(b, uh.vector)
    uh.x.scatter_forward()

    # Assemble and solve second system
    dolfinx.fem.assemble_matrix(A2, a2)
    A2.assemble()
    dolfinx.fem.assemble_vector(b2, L2)
    A2.assemble()
    solver2.solve(b2, uh2.vector)

    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "uh2.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(uh2)
2 Likes
9

Thanks for the example! It helps a lot.
One more thing. I can reuse the same solver by issuing solver.reset(), right?

10

Ive not used that command, so I cannot tell you if you should or should not use it, it Depends on your application/problem.

11

Hello, thanks for this example.

Would it be possible to have the updated version for the latest dolfinx? I’m running into issues with the subspace and with the ‘assemble’ command.

Many thanks

12

I had to wrap assemble arguments in dolfinx.fem.form to make it work with recent DolfinX version.
Maybe you’ll find something useful here.

13

Thanks, that did the trick!