Large Mesh Error PETSc Error 76 - Adjoint Calculation

Good morning,

I have attached a minimal running example for which I get the following error:

*** -------------------------------------------------------------------------
*** 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 /Users/travis/miniconda3/conda-bld/fenics_1520063526556/work/dolfin-2017.2.0.post0/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 2017.2.0
*** Git changeset:
*** ————————————————————————————————————

The error only appears when I increase the number of cells above 30,000.

The error specifically occurs during the calculation of the gradient dJdX5 (line 84).

Could you please help me with this error? I would be very grateful.

Thanks

from dolfin import *
from dolfin_adjoint import *

eps = 1e-4
E, nu = 209*10**9, 0.3                        # Young Modulus
mu = E/(2.0*(1.0 + nu))                        # Shear Modulus
lmbda = E*nu/((1.0 + nu)*(1.0 -2.0*nu))       # Lambda

def epsilon(u):
	return as_tensor(0.5*(nabla_grad(u) + nabla_grad(u).T))

def sigma(u):
    return lmbda*tr(nabla_grad(u))*Identity(3) + 2*mu*epsilon(u) ## Stress

def simp(X5):
	return eps + (1 - eps)*X5**pp

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(3, 1, 1),90,30,30)

A = FunctionSpace(mesh, "CG", 1)
P = VectorFunctionSpace(mesh, "CG", 1)
global pp
pp = 3

class PosX(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],3)
    
class NegX(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0)
    
class PosY(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1],1)
    
class NegY(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1],0)

class PosZ(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2],1)
    
class NegZ(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[2],0)

PosX = PosX()
NegX = NegX()
PosY = PosY()
NegY = NegY()
PosZ = PosZ()
NegZ = NegZ()

materials = MeshFunction('size_t', mesh, 2)
materials.set_all(0)
PosX.mark(materials, 1)
NegX.mark(materials, 2)
PosY.mark(materials, 3)
NegY.mark(materials, 4)
PosZ.mark(materials, 5)
NegZ.mark(materials, 6)
ds = Measure("ds")(subdomain_data=materials)

bc1 = DirichletBC(P.sub(0), 0, materials, 2)
bc2 = DirichletBC(P.sub(1), 0, materials, 2)
bc3 = DirichletBC(P.sub(2), 0, materials, 2)
bc = [bc1,bc2,bc3]
T = Constant((0,-1e3,0))

X5 = interpolate(Constant(0.2),A)

u = TrialFunction(P)
v = TestFunction(P)
a = simp(X5)*inner(sigma(u), grad(v))*dx
L = dot(T,v)*ds(3)
u = Function(P)
solve(a == L, u, bc, solver_parameters={"linear_solver": "mumps"})

J = assemble(inner(T, u)*ds(3))
X5_control = Control(X5)
dJdX5 = compute_gradient(J,X5_control)

I sometimes get PETSc Error 76 when my machine runs out of memory.
Have you checked if this is the case?
I would make sense if it only happens for large number of cells.

Dear @Dave ,

This is a known dolfin-adjoint issue. The adjoint system does not use mumps, and therefore run out of memory. We are in the process of fixing this, and a branch and corresponding PullRequest has been created: https://bitbucket.org/dolfin-adjoint/pyadjoint/pull-requests/85 .
Could you try to use this branch and see if your problem is resolved?