I'm solving Poisson's equation on a 3d BoxMesh, using Lagrangian polynoms of order 1 and then computing the derivative of the solution. Using a grid size of (25, 25, 25) works fine however when increasing to (50, 50, 50) it crashes with the following error:
Calling DOLFIN just-in-time (JIT) compiler, this may take some time.
--- Instant: compiling ---
Solving linear variational problem.
UMFPACK V5.7.1 (Oct 10, 2014): ERROR: out of memory
Traceback (most recent call last):
File "test_out_of_memory.py", line 49, in <module>
derivative = project(phi.dx(0), function_space)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/projection.py", line 147, in project
cpp.la_solve(A, function.vector(), b, solver_type, preconditioner_type)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py", line 5298, in la_solve
return _la.la_solve(*args)
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 /build/dolfin-4SStI2/dolfin-2016.2.0/dolfin/la/PETScLUSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2016.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
Aborted (core dumped)
So it seems to be able to solve the diff.eq. however when computing the derivative it crashes. I monitored memory usage with top and this shows there was still plenty left the moment it crashed:
top - 15:16:18 up 4:27, 1 user, load average: 1,35, 0,87, 0,75
Tasks: 1 total, 1 running, 0 sleeping, 0 stopped, 0 zombie
%Cpu(s): 28,2 us, 1,8 sy, 0,0 ni, 69,6 id, 0,3 wa, 0,0 hi, 0,0 si, 0,0 st
KiB Mem : 16341316 total, 7292176 free, 6335008 used, 2714132 buff/cache
KiB Swap: 16685052 total, 16685052 free, 0 used. 9409472 avail Mem
PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND
11205 dominik 20 0 3191376 2,322g 54244 R 100,0 14,9 1:41.47 python
This is the last entry before it crashed. I'm running on Ubuntu 16.04 (Linux MyPC 4.4.0-66-generic #87-Ubuntu SMP Fri Mar 3 15:29:05 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux) if this is of any importance.
This is the code which I used:
from dolfin.common.constants import DOLFIN_EPS
from dolfin.cpp.mesh import BoxMesh, Point
from dolfin.fem.bcs import DirichletBC
from dolfin.fem.solving import solve
from dolfin.functions.constant import Constant
from dolfin.functions.expression import Expression
from dolfin.functions.function import Function, TestFunction, TrialFunction
from dolfin.functions.functionspace import FunctionSpace
from dolfin import dx, grad, inner, project
x_min, x_max = 0., 1.
y_min, y_max = 0., 1.
z_min, z_max = 0., 1.
grid_size = (50, 50, 50)
min_corner = Point(x_min, y_min, z_min)
max_corner = Point(x_max, y_max, z_max)
mesh = BoxMesh(min_corner, max_corner, *grid_size)
function_space = FunctionSpace(mesh, str('Lagrange'), 1)
def boundary_domain(x):
return (x[0] < x_min + DOLFIN_EPS or x[0] > x_max - DOLFIN_EPS or
x[1] < y_min + DOLFIN_EPS or x[1] > y_max - DOLFIN_EPS or
x[2] < z_min + DOLFIN_EPS or x[2] > z_max - DOLFIN_EPS)
boundary_value = Constant(0.)
boundary_condition = DirichletBC(function_space, boundary_value, boundary_domain)
trial_function = TrialFunction(function_space)
test_function = TestFunction(function_space)
source_term = 'exp(- pow(x[0] / 0.2, 2) / 2.0 ' \
'- pow(x[1] / 0.2, 2) / 2.0 ' \
'- pow(x[2] / 0.2, 2) / 2.0)'
source_term = Expression(str(source_term), degree=3)
bilinear_form = inner(grad(trial_function), grad(test_function)) * dx
linear_form = source_term * test_function * dx
phi = Function(function_space)
solve(
bilinear_form == linear_form,
phi,
boundary_condition,
solver_parameters={
str('linear_solver'): str('mumps'),
}
)
derivative = project(phi.dx(0), function_space)
