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.