UMFPACK error: out of memory despite system having free memory

1

Hello all,

When I use a large enough mesh (100K+ DOF), I am able to create the following error during a call to solve()

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 5.560e+02 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory

Traceback (most recent call last):
  File "bug_demonstration.py", line 82, in <module>
    solve(Fboth == 0, u, bcs, J=J)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 266, in _solve_varproblem
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside /tmp/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2018.1.0
*** Git changeset:  948dc42cc4e06ed9227d0201ad50f94ac94cbf9f
*** -------------------------------------------------------------------------

Here is the minimum script that produces this error

import sys
from dolfin import *


def get_markers(subdomains):
    """
    Automatically determines and returns which subdomain marker belongs to gel subdomain and which belongs to cell
    subdomain based on # of elements in subdomain
    :param subdomains:
    :return:
    """
    sub_arr = subdomains.array()
    markers = list(set(sub_arr))

    print(list(sub_arr).count(markers[0]))
    print(list(sub_arr).count(markers[1]))
    if list(sub_arr).count(markers[0]) > list(sub_arr).count(markers[1]):
        gel_marker = markers[0]
        cell_marker = markers[1]
    else:
        gel_marker = markers[1]
        cell_marker = markers[0]

    return int(gel_marker), int(cell_marker)

## Solver Parameters
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

## Load Mesh
file_in = sys.argv[1]
mesh = Mesh('ProcessedXMLs/' + file_in + '.xml')
subdomains = MeshFunction("size_t", mesh, 'ProcessedXMLs/' + file_in + '_geometrical.xml')
dx = Measure('dx',domain=mesh, subdomain_data=subdomains)
gmark, cmark = get_markers(subdomains)
print("gmark: ", gmark)
print("cmark: ", cmark)

## Functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
du = TrialFunction(V)
v  = TestFunction(V)
u  = Function(V)

## Boundary conditions
xy_bound = 141.7
z_bound = 120
zero = Constant((0.0, 0.0, 0.0))
bcs = []
sbd = []
sbd.append(CompiledSubDomain("near(x[0], side)", side = 0))
sbd.append(CompiledSubDomain("near(x[1], side)", side = 0))
sbd.append(CompiledSubDomain("near(x[2], side)", side = 0))
sbd.append(CompiledSubDomain("near(x[0], side)", side = xy_bound))
sbd.append(CompiledSubDomain("near(x[1], side)", side = xy_bound))
sbd.append(CompiledSubDomain("near(x[2], side)", side = z_bound))
[bcs.append((DirichletBC(V, zero, sub))) for sub in sbd]

## Parameters and equations
d = len(u)
I = Identity(d)
F = I + grad(u)
B = Constant((0.0, 0.0, 0.0))
T = Constant((0.0, 0.0, 0.0))
E0, shr0, nu0 = 10.0, 1.0, 0.45
E1, shr1, nu1 = 10.0, 0.9, 0.45
mu0, lmbda0 = Constant(E0 / (2 * (1 + nu0))), Constant(E0 * nu0 / ((1 + nu0) * (1 - 2 * nu0)))
mu1, lmbda1 = Constant(E1 / (2 * (1 + nu1))), Constant(E1 * nu1 / ((1 + nu1) * (1 - 2 * nu1)))
Fe0 = variable(F * (1 / shr0) * I)
Fe1 = variable(F * (1 / shr1) * I)
psi0 = 1 / 2 * mu0 * (inner(Fe0, Fe0) - 3 - 2 * ln(det(Fe0))) + 1 / 2 * lmbda0 * (
        1 / 2 * (det(Fe0) ** 2 - 1) - ln(det(Fe0)))
psi1 = 1 / 2 * mu1 * (inner(Fe1, Fe1) - 3 - 2 * ln(det(Fe1))) + 1 / 2 * lmbda1 * (
        1 / 2 * (det(Fe1) ** 2 - 1) - ln(det(Fe1)))
f_int = derivative(psi0 * dx(gmark) + psi1 * dx(cmark), u, v)
f_ext = derivative(dot(B, u) * dx('everywhere') + dot(T, u) * ds, u, v)
Fboth = f_int - f_ext
J = derivative(Fboth, u, du)

## Solve
solve(Fboth == 0, u, bcs, J=J)

I am running fenics in docker on ubunutu. I monitor the container memory usage throughout the simulation and it doesn’t appear to exceed a couple GB. The container is limited to the full system memory of 62 GB. I’ve also run fenics natively on ubuntu. I recieve the same error. I also monitor the system-wide memory usage during the simulation and it doesn’t appear to exceed 20%.

Despite all this, the error indicates that I’m running out of memory somehwere. Is this a bug? This thread discusses a similar problem and it is determined to be a bug but I am not using dolfin-adjoint.

Any feedback would be greatly appreciated, thanks.

2 Likes
2

Hey,

this could indeed be a bug of umfpack. As far as I know this is because fenics is using a 32 bit version of this (someone please correct me if I’m wrong).
What you can do, however, is to use mumps, which is 64 bit, and can support more memory. For a bilinear form a, linear form L and boundary conditions given in bcs, the solve command for this looks like

solve(a==L, u, bcs, solver_parameters={'linear_solver' : 'mumps'})

If you want more freedom in configuring mumps, I’ll refer you to my answer here: Preconditioner not working or bug?

1 Like
3

Hello and thanks for the response.

Whenever I try something like

solve(... solver_parameters={'linear_solver' : 'mumps'})

I receive this error

Traceback (most recent call last):
  File "bug_demonstration.py", line 82, in <module>
    solve(Fboth == 0, u, bcs, J=J, solver_parameters={'linear_solver' : 'mumps'})
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 265, in _solve_varproblem
    solver.parameters.update(solver_parameters)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/parameter/__init__.py", line 40, in update
    self[key] = params[key]
RuntimeError: Parameter not found in Parameters object

Any idea why that might be?

4

It seems like mumps is not installed on your system. Maybe send me the output of

list_lu_solver_methods()

Alternatively, you can directly try to use one of the methods given in there. Not that the default is umfpack, as you have seen.

5

When I call list_lu_solver_methods(), I receive this output:

LU method  |  Description                                                 
--------------------------------------------------------------------------
default    |  default LU solver                                           
mumps      |  MUMPS (MUltifrontal Massively Parallel Sparse direct Solver)
petsc      |  PETSc built in LU solver                                    
superlu    |  SuperLU                                                     
umfpack    |  UMFPACK (Unsymmetric MultiFrontal sparse LU factorization)
6

That is not the problem here.

Since nwest calls solve with an expression like F==0 a NonlinearVariationalSolver is used in the background. The parameters are nested in this case. The linear solver is a parameter of the NewtonSolver parameter dict:

solve(Fboth == 0, u, bcs, J=J, solver_parameters = {"newton_solver":{ "linear_solver" : "mumps"}})

Using this line should at least not give the error from above.

2 Likes
7

Hey,

sorry, its my fault, I didnt look close enough. My code works, if you are solving a linear problem. For a nonlinear problem, like you have, the following should do the trick:

solve(Fboth == 0, u, bcs, J=J, solver_parameters={'newton_solver' : {'linear_solver' : 'mumps'}})

as in klunkean’s answer above.
But then again, I think it’d be easier for you to create a NonlinearVariationalProblem and specify the parameters there.

8

Thanks for the help. I no longer get the “Parameter not found in Parameters object”

However, I still receive similar error. With the solve line:

solve(Fboth == 0, u, bcs, J=J, solver_parameters = {"newton_solver":{ "linear_solver" : "mumps"}})

I receive this error:

Traceback (most recent call last):
  File "bug_demonstration.py", line 82, in <module>
    solve(Fboth == 0, u, bcs, J=J, solver_parameters = {"newton_solver":{ "linear_solver" : "mumps"}})
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 266, in _solve_varproblem
    solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PCSETUP_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
*** Git changeset:  948dc42cc4e06ed9227d0201ad50f94ac94cbf9f

And for what it’s worth, the error pops up more quickly than the original error does.

10

Hello,

Is there any update for the memory limitation of umfpack solver?
If umfpack update has not published, what is the alternative of umfpack solver for a large non-linear problem?
I’ve tried the gmres, but it’s not always convergent.

Thank you in advance.

11

Hi.

Were you able to solve this problem? I am also facing the same issue.

Thanks in advance.

12

Hi. I have used this solve to fix it the same problem, and its work fine to calculate de Displacement solution in my model. But, i still need to calculate the stress tensor and it give me a error, i think, for the model size.
The code is the following and if i activate the lines to calculate the stress it don’t work.:

from dolfin import*
from mshr import*
from fenics import*
from ufl import nabla_div
import matplotlib.pyplot as plt
import numpy as np

mesh = Mesh("Caseron_actual.xml")
#plot(mesh);
#Scaled Variables
E = Constant(20e3) #Módulo de Young
nu = Constant(0.3) #Radio de Poisson
mu = E/(2*(1+nu)) #Módulo de cizalle
lambda_ = E*nu/((1+nu)*(1-2*nu)) #Lamé
rho = Constant(0.027) #densidad

# Strain function
def epsilon(u):
    return sym(nabla_grad(u))

# Stress function
def sigma(u):
    return lambda_*nabla_div(u)*Identity(3) + 2*mu*epsilon(u)
# lambda is a reserved python keyword, naming convention recommends using a single trailing underscore for such cases

#Load
g_z = 0.0
b_z = -rho
g = Constant((0.0,g_z,0.0))
b = Constant((0.0,0.0,b_z))

#Function Spaces
V = VectorFunctionSpace(mesh, "P", 2)
u_tr = TrialFunction(V)
u_test = TestFunction(V)

a = inner(sigma(u_tr),epsilon(u_test))*dx #bilineal form
l = dot(b,u_test)*dx + dot(g, u_test)*ds #lineal form

x_max = 1216.0
y_max = 929.
tol=1e-8
def bottom(x, on_boundary):
    return (on_boundary and near(x[2],0,tol))
bc1 =DirichletBC(V.sub(2), Constant(0.),bottom)

def left(x, on_boundary):
    return (on_boundary and near(x[0], 0,tol))
bc2 =DirichletBC(V.sub(0), Constant(0.),left)

def right(x, on_boundary):
    return (on_boundary and near(x[0], x_max,tol))
bc3 =DirichletBC(V.sub(0), Constant(0.),right)

def north(x, on_boundary):
    return (on_boundary and near(x[1], y_max,tol))
bc4 =DirichletBC(V.sub(1), Constant(0.),north)

def south(x, on_boundary):
    return (on_boundary and near(x[1], 0,tol))
bc5 =DirichletBC(V.sub(1), Constant(0.),south)
bc=[bc1,bc2,bc3,bc4,bc5]

u = Function(V)
#problem = LinearVariationalProblem(a, l, u, bc)
#solver = LinearVariationalSolver(problem)
#solver.solve();
solve(a==l, u, bc, solver_parameters={'linear_solver' : 'mumps'})

#Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
#sig = Function(Vsig, name="Cauchy_Stress")
#sig.assign(project(-sigma(u), Vsig));

file_results = XDMFFile("Verde_CA.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
#file_results.write(sig, 0.)

I hope you can help me

13

I guess you mean that when you uncomment the line

#sig.assign(project(-sigma(u), Vsig));

your code will not work, right? This is due FEniCS using the default linear solver (i.e. umfpack) to do the projection. You can, however, implement this easily manually.

For some ufl expression expr, calling

project(expr, V)

for some FunctionSpace object V is equivalent to calling

trial = TrialFunction(V)
test = TestFunction(V)
a = trial*test*dx  # or inner(trial,test)*dx in case of a VectorFunctionSpace
L = expr*test*dx # or inner(expr, test)*dx for a VectorFunctionSpace
solution = Function(V)
solve(a==L, solution)

and here you can again specifiy the solver parameters for the linear solve.
In particular, for your problem, you can use

trial = TrialFunction(Vsig)
test = TestFunction(Vsig)
a = inner(trial, test) * dx
L = inner(-sigma(u), test) * dx
sig = Function(Vsig, name="Cauchy Stress")
solve(a == L, sig, solver_parameters={"linear_solver": "mumps"})

to compute sigma with the help of the mumps solver.

1 Like
14

Tkanks, it works perfect!

15

Hi friend, im using a more dense 3d mesh, but the solvers that i have used dont works. My solver is:

u = Function(V,name="Displacement")
solve(a==l, u, bc, solver_parameters={"linear_solver" : "gmres", "preconditioner" : "ilu"},form_compiler_parameters={"optimize":True})

And the error is:

imagen
imagen1466×515 78.3 KB

Do you know any solver to fix it?? I hope you can help me

16

Sorry, but I cannot give you a hint for this. Finding a linear solver and preconditioner suitable for a problem is a challenging task. If you are solving linear elasticity problems, you have a symmetric and coercive bilinear form, so I guess that the CG method should be the optimal choice, and you can try to use a symmetric preconditioner (ICC for example).

1 Like