# Restart/interruption petsc

Hello everybody,

Currently, I want to solve the complex-valued Helmholtz equation (including Neumann- and Robin boundary conditions) for thousands of frequency points in a loop. The DOF of my model is in the order of about ~100k (3D-model). Further, I have access to an HPC with much memory. Nevertheless, I face the problem that updating of the bilinear and the linear form (a & L) in a loop (some parameters are frequency-dependent) followed by

``````problem = fem.LinearProblem(aH, L, u=pres, petsc_options={"pc_factor_mat_solver_type":"mumps", "sub_pc_type":"ilu"})
problem.solve()
``````

, works fine for about ~200 frequency points. Then, petsc interrupts with petsc4py.PETSc.Error: error code 55

Has anybody an idea what happens there?
By the way, choosing GMRES or avoiding any pre-conditioner doesn’t change much and I use the latest docker-version of dolfinx …

I would assume that this is due to the fact that you recreate your PETSc solver in a loop.
However, as you have not supplied a minimal working example, I cannot help you much further.
See: Read before posting: How do I get my question answered? - #2 on information about how to create a minimal working example.

That might be the solution. But unfortunately, I have no idea how I can split updating the variational form and the solver parameter to specify the solver before calling the do-loop. Please see the following MWE for clarification:

``````import numpy as np
from dolfinx import fem, Function, FunctionSpace
from dolfinx.generation import UnitSquareMesh
from dolfinx.mesh import MeshTags, locate_entities
from mpi4py import MPI
from ufl import inner, grad, TestFunction, TrialFunction, Measure, dx

# create model and evaluate ds
mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# define function space
V = FunctionSpace(mesh, ("Lagrange", 1))
p = TrialFunction(V)
q = TestFunction(V)
pres = fem.Function(V)

# wave number
c = 353.0
rho = 1.13
frequencies = np.arange(2e2, 5e2, .1)
k = (2*np.pi*frequencies)/float(c)

# Neumann b.c.
neu_bc=Function(V)
g_neu_bc = -1j*c*rho*(-1.0)
neu_bc.interpolate(lambda x: 0*x[0]+g_neu_bc)

# Robin b.c. (frequency dependent value)
def rf(c, rho, k, a):
ka = k*a
z_r = rho*c * ((ka**2)/(1+(ka)**2))
z_i = rho*c * (ka/(1+(ka)**2))
z_c = complex(z_r,z_i)
return z_c

# iterate solution
for ii, freq in enumerate(frequencies):
aH = (inner(grad(p), grad(q)) - k[ii]**2 * inner(p, q)) * dx
robin_f=rf(c, rho, k[ii], 30.)
aH += k[ii]*1j*rho*c/robin_f*inner(p,q)*ds(2)
L = k[ii]*inner(neu_bc, q)*ds(1)

problem = fem.LinearProblem(aH, L, u=pres, petsc_options={"pc_factor_mat_solver_type":"mumps"})
#    problem = fem.LinearProblem(aH, L, u=pres)
problem.solve()
print(f'{ii}/{np.shape(frequencies)[0]}')
``````

By the way, no specification of the solver, i.e.

``````problem = fem.LinearProblem(aH, L, u=pres)
``````

, does not make any problems so far.

I would strongly recommend recasting your problem using constants:

``````import numpy as np
from dolfinx import fem
from dolfinx.generation import UnitSquareMesh
from dolfinx.mesh import MeshTags, locate_entities
from mpi4py import MPI
from ufl import inner, grad, TestFunction, TrialFunction, Measure, dx
from petsc4py import PETSc

# create model and evaluate ds
mesh = UnitSquareMesh(MPI.COMM_WORLD, 10, 10)
boundaries = [(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], 1)),
(3, lambda x: np.isclose(x[1], 0)),
(4, lambda x: np.isclose(x[1], 1))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = MeshTags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# define function space
V = fem.FunctionSpace(mesh, ("Lagrange", 1))
p = TrialFunction(V)
q = TestFunction(V)
pres = fem.Function(V)

# wave number
c = 353.0
rho = 1.13
frequencies = np.arange(2e2, 5e2, .1)
k = (2*np.pi*frequencies)/float(c)

# Neumann b.c.
neu_bc=fem.Function(V)
g_neu_bc = -1j*c*rho*(-1.0)
neu_bc.interpolate(lambda x: 0*x[0]+g_neu_bc)

# Robin b.c. (frequency dependent value)
def rf(c, rho, k, a):
ka = k*a
z_r = rho*c * ((ka**2)/(1+(ka)**2))
z_i = rho*c * (ka/(1+(ka)**2))
z_c = complex(z_r,z_i)
return z_c

kc = fem.Constant(mesh, PETSc.ScalarType(0))
aH = (inner(grad(p), grad(q)) - kc**2 * inner(p, q)) * dx
robin_c = fem.Constant(mesh, PETSc.ScalarType(0))
aH += kc*1j*rho*c/robin_c*inner(p,q)*ds(2)
L = kc*inner(neu_bc, q)*ds(1)
problem = fem.LinearProblem(aH, L, u=pres)#, petsc_options={"pc_factor_mat_solver_type":"mumps","sub_pc_type":"ilu"})

for ii, freq in enumerate(frequencies):
robin_f=rf(c, rho, k[ii], 30.)
kc.value = k[ii]
robin_c.value = robin_f
problem.solve()
print(f'{ii}/{np.shape(frequencies)[0]}')
``````

Doing so will avoid the re-creation of matrices, recompilation of variational forms etc.
I note that using your suggested options `petsc_options={"pc_factor_mat_solver_type":"mumps","sub_pc_type":"ilu"}` throws:

``````[0] KSPSetFromOptions() at /usr/local/petsc/src/ksp/ksp/interface/itcl.c:355
[0] PCSetFromOptions() at /usr/local/petsc/src/ksp/pc/interface/pcset.c:165
[0] PCSetFromOptions_ILU() at /usr/local/petsc/src/ksp/pc/impls/factor/ilu/ilu.c:55
[0] PCSetFromOptions_Factor() at /usr/local/petsc/src/ksp/pc/impls/factor/factimpl.c:272
[0] PCFactorSetDefaultOrdering_Factor() at /usr/local/petsc/src/ksp/pc/impls/factor/factor.c:20
[0] MatGetFactor() at /usr/local/petsc/src/mat/interface/matrix.c:4780
[0] See https://petsc.org/release/overview/linear_solve_table/ for possible LU and Cholesky solvers
[0] MatSolverType mumps does not support factorization type ILU for matrix type seqaij
``````

However, you can use `"superlu"` with no issues

Perfect, this is truly an elegant solution.