Solver failed with "PC failed due to FACTOR_OUTMEMORY" (MUMPS)

[lrtfm]

I installed firedrake with the latest version and got the following error. The code works well on the old version of the firedrake.

Error information

  Linear ns_ solve did not converge due to DIVERGED_PC_FAILED iterations 0
                 PC failed due to FACTOR_OUTMEMORY
...
[0] SNESSolve() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/snes/interface/snes.c:4810
[0] SNESSolve_KSPONLY() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/snes/impls/ksponly/ksponly.c:51
[0] KSPSolve() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/ksp/ksp/interface/itfunc.c:1103
[0] KSPSolve_Private() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/ksp/ksp/interface/itfunc.c:863
[0] KSPSetUp() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/ksp/ksp/interface/itfunc.c:417
[0] PCSetUp() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/ksp/pc/interface/precon.c:1017
[0] PCSetUp_LU() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/ksp/pc/impls/factor/lu/lu.c:133
[0] MatLUFactorNumeric() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/mat/interface/matrix.c:3186
[0] MatFactorNumeric_MUMPS() at /home/bpduan/zzyang/workspace/firedrake-gcc/firedrake-mini-petsc/src/petsc/src/mat/impls/aij/mpi/mumps/mumps.c:1688
[0] Error in external library
[0] Error reported by MUMPS in numerical factorization phase: INFOG(1)=-9, INFO(2)=2081948

ksp_view output

The ksp_view output is:

KSP Object: (ns_) 1 MPI processes
  type: preonly
  maximum iterations=10000, initial guess is zero
  tolerances:  relative=1e-05, absolute=1e-50, divergence=10000.
  left preconditioning
  using NONE norm type for convergence test
PC Object: (ns_) 1 MPI processes
  type: lu
    out-of-place factorization
    tolerance for zero pivot 2.22045e-14
    matrix ordering: external
    factor fill ratio given 0., needed 0.
      Factored matrix follows:
        Mat Object: 1 MPI processes
          type: mumps
          rows=4564, cols=4564
          package used to perform factorization: mumps
          total: nonzeros=719272, allocated nonzeros=719272
            MUMPS run parameters:
              SYM (matrix type):                   0
              PAR (host participation):            1
              ICNTL(1) (output for error):         6
              ICNTL(2) (output of diagnostic msg): 0
              ICNTL(3) (output for global info):   0
              ICNTL(4) (level of printing):        0
              ICNTL(5) (input mat struct):         0
              ICNTL(6) (matrix prescaling):        7
              ICNTL(7) (sequential matrix ordering):7
              ICNTL(8) (scaling strategy):        77
              ICNTL(10) (max num of refinements):  0
              ICNTL(11) (error analysis):          0
              ICNTL(12) (efficiency control):                         1
              ICNTL(13) (sequential factorization of the root node):  0
              ICNTL(14) (percentage of estimated workspace increase): 200
              ICNTL(18) (input mat struct):                           0
              ICNTL(19) (Schur complement info):                      0
              ICNTL(20) (RHS sparse pattern):                         0
              ICNTL(21) (solution struct):                            0
              ICNTL(22) (in-core/out-of-core facility):               0
              ICNTL(23) (max size of memory can be allocated locally):0
              ICNTL(24) (detection of null pivot rows):               0
              ICNTL(25) (computation of a null space basis):          0
              ICNTL(26) (Schur options for RHS or solution):          0
              ICNTL(27) (blocking size for multiple RHS):             -32
              ICNTL(28) (use parallel or sequential ordering):        1
              ICNTL(29) (parallel ordering):                          0
              ICNTL(30) (user-specified set of entries in inv(A)):    0
              ICNTL(31) (factors is discarded in the solve phase):    0
              ICNTL(33) (compute determinant):                        0
              ICNTL(35) (activate BLR based factorization):           0
              ICNTL(36) (choice of BLR factorization variant):        0
              ICNTL(38) (estimated compression rate of LU factors):   333
              CNTL(1) (relative pivoting threshold):      0.01
              CNTL(2) (stopping criterion of refinement): 1.49012e-08
              CNTL(3) (absolute pivoting threshold):      0.
              CNTL(4) (value of static pivoting):         -1.
              CNTL(5) (fixation for null pivots):         0.
              CNTL(7) (dropping parameter for BLR):       0.
              RINFO(1) (local estimated flops for the elimination after analysis):
                [0] 7.86601e+07
              RINFO(2) (local estimated flops for the assembly after factorization):
                [0]  4.60763e+06
              RINFO(3) (local estimated flops for the elimination after factorization):
                [0]  8.31715e+08
              INFO(15) (estimated size of (in MB) MUMPS internal data for running numerical factorization):
              [0] 22
              INFO(16) (size of (in MB) MUMPS internal data used during numerical factorization):
                [0] 29
              INFO(23) (num of pivots eliminated on this processor after factorization):
                [0] 1958
              RINFOG(1) (global estimated flops for the elimination after analysis): 7.86601e+07
              RINFOG(2) (global estimated flops for the assembly after factorization): 4.60763e+06
              RINFOG(3) (global estimated flops for the elimination after factorization): 8.31715e+08
              (RINFOG(12) RINFOG(13))*2^INFOG(34) (determinant): (0.,0.)*(2^0)
              INFOG(3) (estimated real workspace for factors on all processors after analysis): 719272
              INFOG(4) (estimated integer workspace for factors on all processors after analysis): 52460
              INFOG(5) (estimated maximum front size in the complete tree): 226
              INFOG(6) (number of nodes in the complete tree): 401
              INFOG(7) (ordering option effectively used after analysis): 2
              INFOG(8) (structural symmetry in percent of the permuted matrix after analysis): 100
              INFOG(9) (total real/complex workspace to store the matrix factors after factorization): 1421640
              INFOG(10) (total integer space store the matrix factors after factorization): 52604
              INFOG(11) (order of largest frontal matrix after factorization): 980
              INFOG(12) (number of off-diagonal pivots): 503
              INFOG(13) (number of delayed pivots after factorization): 9558
              INFOG(14) (number of memory compress after factorization): 4
              INFOG(15) (number of steps of iterative refinement after solution): 0
              INFOG(16) (estimated size (in MB) of all MUMPS internal data for factorization after analysis: value on the most memory consuming processor): 22
              INFOG(17) (estimated size of all MUMPS internal data for factorization after analysis: sum over all processors): 22
              INFOG(18) (size of all MUMPS internal data allocated during factorization: value on the most memory consuming processor): 29
              INFOG(19) (size of all MUMPS internal data allocated during factorization: sum over all processors): 29
              INFOG(20) (estimated number of entries in the factors): 719272
              INFOG(21) (size in MB of memory effectively used during factorization - value on the most memory consuming processor): 28
              INFOG(22) (size in MB of memory effectively used during factorization - sum over all processors): 28
              INFOG(23) (after analysis: value of ICNTL(6) effectively used): 0
              INFOG(24) (after analysis: value of ICNTL(12) effectively used): 1
              INFOG(25) (after factorization: number of pivots modified by static pivoting): 0
              INFOG(28) (after factorization: number of null pivots encountered): 0
              INFOG(29) (after factorization: effective number of entries in the factors (sum over all processors)): 1421640
              INFOG(30, 31) (after solution: size in Mbytes of memory used during solution phase): 0, 0
              INFOG(32) (after analysis: type of analysis done): 1
              INFOG(33) (value used for ICNTL(8)): 7
              INFOG(34) (exponent of the determinant if determinant is requested): 0
              INFOG(35) (after factorization: number of entries taking into account BLR factor compression - sum over all processors): 1421640
              INFOG(36) (after analysis: estimated size of all MUMPS internal data for running BLR in-core - value on the most memory consuming processor): 0
              INFOG(37) (after analysis: estimated size of all MUMPS internal data for running BLR in-core - sum over all processors): 0
              INFOG(38) (after analysis: estimated size of all MUMPS internal data for running BLR out-of-core - value on the most memory consuming processor): 0
              INFOG(39) (after analysis: estimated size of all MUMPS internal data for running BLR out-of-core - sum over all processors): 0
  linear system matrix = precond matrix:
  Mat Object: (ns_) 1 MPI processes
    type: seqaij
    rows=4564, cols=4564
    total: nonzeros=161968, allocated nonzeros=161968
    total number of mallocs used during MatSetValues calls=0
      using I-node routines: found 2751 nodes, limit used is 5

Information on the versions of Firedrake and PETSc

Latest version of firedrake

My code works well on old version firedrake:

commit 5d1a8512da6391d8373e4280046bdf51c0a3efee (HEAD -> master, origin/master, origin/HEAD)
Merge: 16bbad8f 5cea71c2
Author: David A Ham <[email protected]>
Date:   Fri Dec 17 10:23:53 2021 +0000

    Merge pull request #2302 from firedrakeproject/JDBetteridge/fix_blas_config_json_mac

    Mac Bug Squash

with petsc:

commit 70a4eb5794b0eb30a077b74a20497d039985a498 (HEAD -> firedrake, tag: Firedrake_20211026.0, tag: Firedrake_20211025.0, tag: Firedrake_20210924.0, origin/firedrake, origin/HEAD)
Merge: 0ca58fa61f bacf669493
Author: David Ham <[email protected]>
Date:   Thu Sep 9 16:44:00 2021 +0100

    Merge pull request #6 from firedrakeproject/ksagiyam/merge_upstream

    Ksagiyam/merge upstream

new version firedrake and PETSc

Firedrake:

commit ac22e4c55d6fad32ddc9e936cd3674fb8a75f1da (HEAD -> master, origin/master, origin/HEAD)
Merge: 95fe0efd 019cf585
Author: Lawrence Mitchell <[email protected]>
Date:   Wed Mar 16 15:22:13 2022 +0000

    Merge pull request #2378 from firedrakeproject/fabianl1908/fix-test-variant

    Fix interpolation test for Nedelec, BDM and RT elements for variant=integral

PETSc:

commit ba68173f42851ea992f421766dabf4d5e6291972 (HEAD -> firedrake, origin/firedrake, origin/HEAD)
Merge: 70a4eb5794 5f24bd3434
Author: Lawrence Mitchell <[email protected]>
Date:   Fri Feb 25 16:26:32 2022 +0000

    Merge pull request #7 from firedrakeproject/ksagiyam/merge_upstream

    Ksagiyam/merge upstream

A minimal failing example

In this example, the pressure of the NS equation is parted into two: one for water and the other one for bubble. I copy the code from https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html to apply the boundary conditions for the two pressures.

#!/usr/bin/env python
# coding: utf-8
​
# A failing example
​
from firedrake import *
from firedrake.petsc import PETSc
import numpy as np
​
# Ref: https://www.firedrakeproject.org/demos/linear_fluid_structure_interaction.py.html
class SubManifoldBC(DirichletBC):
    """A class help to make a sub-domain or sub-manifold, such as a surface in 3D domain, a line in 2D domain"""
    def __init__(self, V, value, markers):
        """
        :param V: the :class:`.FunctionSpace`
        :param value: the boundary value
        :param markers: a scalar function of which the zeros values indicate bounday nodes. 
        """
        # Call superclass init
        # We provide a dummy subdomain id.
        super(SubManifoldBC, self).__init__(V, value, 0)
        # Override the "nodes" property which says where the boundary
        # condition is to be applied.
        self.nodes = np.unique(np.where(markers.dat.data_ro_with_halos == 0)[0])
        self._markers = markers
​
def make_domain_marker(mesh, index, name=None, dtype=np.int32):
    """Maker the cell with id=index to 1"""
    V = FunctionSpace(mesh, 'DG', 0)
    marker = Function(V, name=name, dtype=dtype)
    par_loop(('{[i] : 0 <= i < f.dofs}', 'f[i, 0] = 1'), 
             dx(index), {'f': (marker, WRITE)}, is_loopy_kernel=True)
    return marker
​
def make_dof_marker(V, index, name=None, dtype=np.int32):
    """Maker the dofs in cell with id=index to 1, and leave others zero"""
    dmarker = make_domain_marker(V.ufl_domain(), index, dtype=dtype)
    dof_marker = Function(V, name=name, dtype=dtype)
    par_loop(('{[i, j] : 0 <= i < A.dofs and 0 <= j < %d}' % V.value_size, 
              'A[i, j] = if(A[i, j] > B[0, 0], A[i, j], B[0, 0])'),
        dx, {'A' : (dof_marker, INC), 'B' : (dmarker, READ)}, is_loopy_kernel=True)
    dof_marker.dat.data_with_halos[dof_marker.dat.data_with_halos[:] > 0] = 1
    return dof_marker
​
​
################################
#  Domain
# ----------------
# |              |
# |              |
# |              |
# |     water 3  |
# |              |
# |      __      |
# |    / 4  \    |
# |   |bubble|   |
# |    \ __ /    |
# |              |
# ----------------
# Rectangle: [0, 1] X [0, 2]
# Circle: center (0.5, 0.5) with R=0.25
################################
mesh = Mesh('bubble.msh')
water_id = 3
bubble_id = 4
# interface_id = 13   # The circle
​
# space for fluid velocity, presure, curvature 
U = VectorFunctionSpace(mesh, 'CG', 2) # velocity
Q = FunctionSpace(mesh, 'CG', 1)       # p1 for water and p2 for bubble
V = MixedFunctionSpace([U, Q, Q])
​
# mark subdomains
I_water = make_domain_marker(mesh, water_id, name='I_water')
I_bubble = make_domain_marker(mesh, bubble_id, name='I_bubble')
​
I_p_water = make_dof_marker(Q, water_id, name='I_p_water')
I_p_bubble = make_dof_marker(Q, bubble_id, name='I_p_bubble')
​
# set soome parameters
rho = 1000.*I_water + 100.*I_bubble
mu = 10.*I_water + 1.*I_bubble
dt = 0.01
eps = 1e-8
​
u_, p_water_, p_bubble_ = TrialFunctions(V)
v_, q_water_, q_bubble_ = TestFunctions(V)
​
# A11 = dot(rho*u_/dt, v_)*dx \
A11 =   + 2*inner(mu*sym(grad(u_)), sym(grad(v_)))*dx       \
        - dot(p_water_, div(v_))*dx(water_id) \
        - dot(p_bubble_, div(v_))*dx(bubble_id)
A12 = dot(div(u_), q_water_)*dx(water_id) \
        + dot(div(u_),  q_bubble_)*dx(bubble_id)     \
        + eps*inner(p_water_, q_water_)*dx(water_id) \
        + eps*inner(p_bubble_, q_bubble_)*dx(bubble_id)
​
A = A11 + A12
​
BC_p_water = SubManifoldBC(V.sub(1), 0., I_p_water)
BC_p_bubble = SubManifoldBC(V.sub(2), 0., I_p_bubble)
​
BC_u = DirichletBC(V.sub(0), 0., 'on_boundary')
bcs =  [BC_p_water, BC_p_bubble, BC_u]
​
# print min and max value of the diagnoal
M = assemble(A, bcs=bcs)
J = M.petscmat
d = J.getDiagonal()
d.abs()
print('min and max value on diag:', d.min(), d.max())

sol = Function(V)
# for test, we just set L = 0
ns_problem = LinearVariationalProblem(A, 0, sol, bcs=bcs)
ns_solver = LinearVariationalSolver(ns_problem, options_prefix='ns')
ns_solver.solve()

ns_mumps_failed.zip

[knepley]

From the MUMPS manual: "If an error –9
occurs, the user should increase the value of ICNTL(14) before calling the factorization (JOB=2)
again, except if LWK USER is provided LWK USER should be increased."

[lrtfm]

@knepley Thanks for your suggestion. I will try to debug in this direction.

[knepley]

I think you just need to give that option to increase memory. If it ran before, it is likely the sparsity pattern has changed.

[lrtfm]

I try to increase memory and PC failed due to FACTOR_NUMERIC_ZEROPIVOT rise.
And then I print the max and min value of the diagonal's absolute value. The value is different from before and changes every time. It seems the matrix did not initialize properly.

Output on new Firedrake:

In [1]: %run ns_mumps_fail.py
min and max value on diag: (0, 0.0) (1720, inf)

In [2]: %run ns_mumps_fail.py
min and max value on diag: (0, 0.0) (3368, 141140.44334537265)

In [3]: %run ns_mumps_fail.py
min and max value on diag: (0, 0.0) (1811, 1.4256938990112674e+268)

In [4]: %run ns_mumps_fail.py
min and max value on diag: (0, 0.0) (1738, inf)

In [5]: %run ns_mumps_fail.py
min and max value on diag: (0, 0.0) (1630, inf)

Output of the old version:

In [1]: %run ns_mumps_fail.py

min and max value on diag: (4326, 3.8318714356058965e-12) (3505, 264.25708153824326)

[lrtfm]

After lots of trial and error, I found the main reason: I used a Function with type np.int32 in the variational form.
After changing it to np.float64, everything seems to go back to normal. See the following example and its outputs:

Using a Function with a data type different from PetscScalar should be forbidden in variational form?

An example

#!/usr/bin/env python
# coding: utf-8

from firedrake import *
from firedrake.petsc import PETSc
import numpy as np

mesh = RectangleMesh(10, 10, 1, 1)
U = FunctionSpace(mesh, 'CG', 1)
u, v = TrialFunction(U), TestFunction(U)

type_map = {'int': np.int32, 'double': np.float64}
As = {}
for name, dtype in type_map.items():
    rho = Function(U, dtype=dtype)
    A = dot(rho*u, v)*dx
    M = assemble(A)
    J = M.petscmat
    d = J.getDiagonal()
    d.abs()
    print(name, '\n\t min and max value on diag:', d.min(), d.max())

Its outputs

In [1]: %run mumps_fail.py
int
         min and max value on diag: (0, 0.0) (110, 1.8447537411794655e-302)
double
         min and max value on diag: (0, 0.0) (0, 0.0)

In [2]: %run mumps_fail.py
int
         min and max value on diag: (0, 0.0) (73, 3.191999999999991)
double
         min and max value on diag: (0, 0.0) (0, 0.0)

In [3]: %run mumps_fail.py
int
         min and max value on diag: (0, 0.0) (73, 3.191999999999991)
double
         min and max value on diag: (0, 0.0) (0, 0.0)

In [4]: %run mumps_fail.py
int
         min and max value on diag: (0, 0.0) (0, 0.0)
double
         min and max value on diag: (0, 0.0) (0, 0.0)

[lrtfm]

@connorjward, could you explain why the previous version works well? Just curious to know.

[connorjward]

Sorry no I don't know when this behaviour got introduced.

[lrtfm]

Never mind, thanks anyway!

[wence-]

I've prepared a 'fix'. OP2/PyOP2#659 should hopefully raise a somewhat helpful error for this case. Arguably we should also do some processing of the form, perhaps here, to catch this error at an earlier point.

At the point when you're compiling a form there's no guarantee that your Coefficient is a Function (so you can't rely on a dtype property). That said, assemble can inspect the form.coefficients() and barf in any of them have a non-PetscScalar dtype.

[connorjward]

Done: 464c8a6