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)
```