Restart/interruption petsc

dolfinx
1

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 …

2

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.

3

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.

4

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

5

Perfect, this is truly an elegant solution.