Dear all,

I’ve tried to make parameters optimization with pyadjoint for a, quite heavy, bio-poromechanical system, and everything is fine until an error occurs with the project function. Depending on the version of the code, the error doesn’t occurs at the same iteration (638, 1015, 1048, 2035, …). It seems that simpler is the code, greater the number of iterations done.

I’m on Ubuntu 18.04.2 LTS and I’m using FEniCS and Pyadjoint on Docker.

The error is the following, occurs at iteration 2035 - takes 30 seconds running with the MWE - :

File “1flu_1sol_dolfin.py”, line 81, in <module>

p_P1 = project(_p, L)

File “/usr/local/lib/python3.6/dist-packages/fenics_adjoint/projection.py”, line 24, in project

block = ProjectBlock(args[0], args[1], output, bcs, **sb_kwargs)

File “/usr/local/lib/python3.6/dist-packages/fenics_adjoint/projection.py”, line 45, ininit

super(ProjectBlock, self).init(a == L, output, bcs, *args, **kwargs)

File “/usr/local/lib/python3.6/dist-packages/fenics_adjoint/solving.py”, line 76, ininit

self._init_dependencies(*args, **kwargs)

File “/usr/local/lib/python3.6/dist-packages/fenics_adjoint/solving.py”, line 145, in _init_dependencies

self.add_dependency(c.block_variable, no_duplicates=True)

File “/usr/local/lib/python3.6/dist-packages/pyadjoint/block.py”, line 50, in add_dependency

dep.will_add_as_dependency()

File “/usr/local/lib/python3.6/dist-packages/pyadjoint/block_variable.py”, line 63, in will_add_as_dependency

self.save_output(overwrite=overwrite)

File “/usr/local/lib/python3.6/dist-packages/pyadjoint/tape.py”, line 46, in wrapper

return function(*args, **kwargs)

File “/usr/local/lib/python3.6/dist-packages/pyadjoint/block_variable.py”, line 51, in save_output

self._checkpoint = self.output._ad_create_checkpoint()

File “/usr/local/lib/python3.6/dist-packages/fenics_adjoint/types/function.py”, line 142, in _ad_create_checkpoint

deepcopy=True)

File “/usr/local/lib/python3.6/dist-packages/dolfin/function/function.py”, line 526, in sub

self.cpp_object().sub(i),

RuntimeError: *** Error: Duplication of MPI communicator failed (MPI_Comm_dup

I’ve tried to make a minimal working example, but it’s 80 lines long. It’s the classical Terzaghi consolidation problem of poromechanics.

It’s a mixed formulation problem, known and benchmarked. It’s almost sure that the error isn’t a math problem. Here is the code:

```
from dolfin import *
from dolfin_adjoint import *
from mshr import *
from ufl import nabla_div
# Terzaghi's problem
T, num_steps = 100.0, 10000
dt = T / (1.0*num_steps)
rho_l, mu_l, k = 1000, 1e-2, 1.8e-15
E, nu = 5000, 0.4
lambda_, mu =(E*nu)/((1+nu)*(1-2*nu)), E/(2*(1+nu))
Kf, Ks, poro=2.2e9, 1e10, 0.2
S=(poro/Kf)+(1-poro)/Ks
# Create mesh and define expression
mesh=RectangleMesh(Point(0.0,0.0), Point(1e-5, 1e-4),8,16,'crossed')
top = CompiledSubDomain("near(x[1], side)", side = 1e-4)
bottom = CompiledSubDomain("near(x[1], side)", side = 0.0)
left_right = CompiledSubDomain("(near(x[0], 0.0) )|| (near(x[0], 1e-5) )")
boundary = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary.set_all(1)
top.mark(boundary, 2)
ds = Measure('ds', domain = mesh, subdomain_data = boundary)
# Load on the top boundary
T = Expression(('0.0', '-100'),degree=1)
#Define Mixed Space (R2,R) -> (u,p)
V = VectorElement("CG", mesh.ufl_cell(), 2)
W = FiniteElement("CG", mesh.ufl_cell(), 1)
L = dolfin.FunctionSpace(mesh,W)
MS = dolfin.FunctionSpace(mesh, MixedElement([V,W]))
# Define boundary condition
c = Expression(("0.0"), degree=1)
bcu1 = DirichletBC(MS.sub(0).sub(0), c, left_right)
bcu2 = DirichletBC(MS.sub(0).sub(1), c, bottom)
bcp = DirichletBC(MS.sub(1), c, top)
bc=[bcu1,bcu2,bcp]
# Define strain
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# Define variational problem and initial condition
X0 = Function(MS)
B = TestFunction(MS)
e_u0 = Expression(('0.0', '0.0'), degree=1)
e_p0 = Expression('100.0', degree=1)
u0 = interpolate(e_u0, MS.sub(0).collapse())
p0 = interpolate(e_p0, MS.sub(1).collapse())
Xn = Function(MS)
assign(Xn, [u0, p0])
(u,p)=split(X0)
(u_n,p_n)=split(Xn)
(v,q)=split(B)
F = (1/dt)*nabla_div(u-u_n)*q*dx + (k/(mu_l))*dot(grad(p),grad(q))*dx
F += inner(2*mu*epsilon(u),epsilon(v))*dx + lambda_*nabla_div(u)*nabla_div(v)*dx -p*nabla_div(v)*dx- dot(T,v)*ds(2)
t = 0
for n in range(num_steps):
print('n=',n)
t += dt
solve(F == 0, X0, bc)
(_u,_p)=X0.split()
# the error occurs here
p_P1 = project(_p, L)
assign(Xn,X0)
```

If the lines:

```
from dolfin import *
from dolfin_adjoint import *
```

are replaced by

`from fenics import *`

the 10,000 iterations are done.

Have you a lead to avoid this error ? Is it an identified bug ?

Thank you for your time