Out of memory error

Linear Algebra
1

Greetings,

I am trying to compute electrostatic calculations on a mesh containing 170135 vertices and 973248 cells, I am getting an error:

    WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.
    Solving linear variational problem.
      *** Warning: Found no facets matching domain for boundary condition.
      *** Warning: Found no facets matching domain for boundary condition.

    UMFPACK V5.7.8 (Nov 9, 2018): ERROR: out of memory

    Traceback (most recent call last):
      File "Python3.py", line 58, in <module>
        fn.solve(a == L, u, bcs)
      File "/home/tadask/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/solving.py", line 220, in solve
        _solve_varproblem(*args, **kwargs)
      File "/home/tadask/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/solving.py", line 247, 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 /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1575559212425/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.
    *** Process: 0
    ***
    *** DOLFIN version: 2019.1.0
    *** Git changeset:  77c758717da1cb8017bf95109363bb1a91a4f8e5
    *** -------------------------------------------------------------------------

Here is the code I am using:

import fenics as fn
import numpy as np

epsilon_0 = 8.85e-12

mesh = fn.Mesh('Smaller_mesh.xml')
subdomains = fn.MeshFunction("size_t", mesh, 'Smaller_mesh_physical_region.xml')
boundaries = fn.MeshFunction('size_t', mesh, 'Smaller_mesh_facet_region.xml')
dx = fn.Measure('dx', domain=mesh, subdomain_data=subdomains)
V0 = fn.FunctionSpace(mesh, 'DG', 0)
V = fn.FunctionSpace(mesh, 'P', 1)


class permittivity(fn.UserExpression):
    def __init__(self, markers, **kwargs):
        self.markers = markers
        super().__init__(**kwargs)
        
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 1:
            values[0] = 5
        elif self.markers[cell.index] == 4:
            values[0] = 6
        elif self.markers[cell.index] == 5:
            values[0] = 10
        elif self.markers[cell.index] == 5:
            values[0] = 5
        else:
            values[0] = 1

eps = permittivity(subdomains, degree=1)


top = fn.DirichletBC(V, fn.Constant(0), boundaries, 3)
bottom = fn.DirichletBC(V, fn.Constant(3000), boundaries, 1)
bcs =[top, bottom]


permitivity = fn.project(eps, V0)
epsfile = fn.File('output6/e_eps.pvd')
epsfile << permitivity

u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * eps * fn.dx
L = fn.Constant(0) * v * fn.dx
u = fn.Function(V)
fn.solve(a == L, u, bcs)

electric_field = fn.project(-fn.grad(u))

potentialFile = fn.File('output6/potential.pvd')
potentialFile << u
e_fieldfile = fn.File('output6/e_field.pvd')
e_fieldfile << electric_field

Is there a way to overcome this issue?

Best wishes,
Ted

2

Many things can be done, like using “mumps” ( solve(a==l, deformation, bcs, solver_parameters={"linear_solver": "mumps"}) as the solver or running in parallel using MPI (mpirun -n 2 python3 script.py)

3

Hi dokken,

I have updated the code:

import fenics as fn
import numpy as np

epsilon_0 = 8.85e-12
# Import the mesh, identify the subdomains and boundaries
# and then make a function space and a dx element.
mesh = fn.Mesh('Smaller_mesh.xml')
subdomains = fn.MeshFunction("size_t", mesh, 'Smaller_mesh_physical_region.xml')
boundaries = fn.MeshFunction('size_t', mesh, 'Smaller_mesh_facet_region.xml')
dx = fn.Measure('dx', domain=mesh, subdomain_data=subdomains)
V0 = fn.FunctionSpace(mesh, 'DG', 0)
V = fn.FunctionSpace(mesh, 'P', 1)


class permittivity(fn.UserExpression):
    def __init__(self, markers, **kwargs):
        self.markers = markers
        super().__init__(**kwargs)
        
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 1:
            values[0] = 5
        elif self.markers[cell.index] == 4:
            values[0] = 6
        elif self.markers[cell.index] == 5:
            values[0] = 10
        elif self.markers[cell.index] == 5:
            values[0] = 5
        else:
            values[0] = 1

eps = permittivity(subdomains, degree=1)



# Set the boundary conditions using the boundary numbers
# specified in Gmsh.

top = fn.DirichletBC(V, fn.Constant(0), boundaries, 3)
bottom = fn.DirichletBC(V, fn.Constant(3000), boundaries, 1)
bcs =[top, bottom]


permitivity = fn.project(eps, V0)
epsfile = fn.File('output6/e_eps.pvd')
epsfile << permitivity


# Solve the Poisson equation with the source set to 0
# (the Laplace equation) via the L = fn.Constant('0')... line
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * eps * fn.dx
L = fn.Constant(0) * v * fn.dx
u = fn.Function(V)
fn.solve(a == L, u, bcs, solver_parameters={"linear_solver": "mumps"})


# Find the electric field by taking the negative of the
# gradient of the electric potential. Project this onto
# the function space for evaluation.
electric_field = fn.project(-fn.grad(u))

# output the results to two separate files
potentialFile = fn.File('output6/potential.pvd')
potentialFile << u
e_fieldfile = fn.File('output6/e_field.pvd')
e_fieldfile << electric_field

  *** Warning: Found no facets matching domain for boundary condition.
  *** Warning: Found no facets matching domain for boundary condition.

UMFPACK V5.7.8 (Nov 9, 2018): ERROR: out of memory

Traceback (most recent call last):
  File "Python3.py", line 63, in <module>
    electric_field = fn.project(-fn.grad(u))
  File "/home/tadask/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/fem/projection.py", line 138, in project
    cpp.la.solve(A, function.vector(), b, solver_type, preconditioner_type)
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 /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1575559212425/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset:  77c758717da1cb8017bf95109363bb1a91a4f8e5
*** -------------------------------------------------------------------------

I am still receiving an error, although at the different code line.

It seems that fn.project does not allow to apply a different solver method. The other method crashes my ubuntu.

4

Please encapsulate your code with ``` to make it readable.
Project allows for solver options, which you can see by opening an IPython shell and run from dolfin import project and project?. Then you obtain the following:

project(
    ['v', 'V=None', 'bcs=None', 'mesh=None', 'function=None', "solver_type='lu'", "preconditioner_type='default'", 'form_compiler_parameters=None'],
)
Source:   
def project(v, V=None, bcs=None, mesh=None,
            function=None,
            solver_type="lu",
            preconditioner_type="default",
            form_compiler_parameters=None):
5

dokken,

I have changed the solver_type in fn.project from default “lu” to “cg”, honestly I have no idea why, but it works. Could you please explain why?

fn.solve(a == L, u, bcs, solver_parameters={"linear_solver": "mumps"})
electric_field = fn.project(-fn.grad(u), solver_type="cg")

Best wishes,
Ted

6

See for instance Tutorial on solvers

1 Like