Hi,
I’m trying to find a solution on polygon mesh. Everything works as expected when using:
from dolfin import *
However, once I follow this with:
from dolfin_adjoint import *
the code will no run and I get the following error:
AttributeError: 'dolfin.cpp.mesh.Mesh' object has no attribute '_ad_will_add_as_dependency'
If I set “annotate=False” in the solve call, it works again. From the dolfin-adjoitn API reference
This solve routine wraps the real Dolfin solve call. Its purpose is to annotate the model, recording what solves occur and what forms are involved, so that the adjoint and tangent linear models may be constructed automatically by pyadjoint.
To disable the annotation, just pass annotate=False to this routine, and it acts exactly like the Dolfin solve call.
Based on this, I assume the issue is coming from the automatic derivation of the adjoint forms in pyadjoint. I need dolfin_adjoint and its annotated version of solve as the solution u will later be used in an objective functional which is not shown in the MWE below. Any insight into why this is happening and/or potential fixes would be greatly appreciated.
I’m running FEniCs 2019.1.0 and dolfin-adjoint 2019.1.0.
Below is a MWE:
from dolfin import *
# from dolfin_adjoint import * # When this is uncommented, the code will not run.
from mshr import Polygon, generate_mesh
import matplotlib.pyplot as plt
domain_vertices = [
Point(0, 0),
Point(1, 0),
Point(1, 1),
Point(0.75, 1.1),
Point(0.7, 1.2),
Point(0.6, 1.2),
Point(0.5, 1.15),
Point(0.2, 1.1),
Point(0, 1),
Point(0, 0),
]
domain = Polygon(domain_vertices)
mesh = generate_mesh(domain,20)
V = FunctionSpace(mesh, 'P', 1)
eps = 1e-14
def boundary_top(x, on_boundary):
return on_boundary and x[0]>1e-14 and x[0]<1-eps and x[1]>1
def boundary_bottom(x, on_boundary):
return on_boundary and near(x[1], 0)
u_top = Constant(2.0)
u_bottom = Constant(1.0)
bc_top = DirichletBC(V,u_top,boundary_top)
bc_bottom = DirichletBC(V,u_bottom,boundary_bottom)
bcs = [bc_top, bc_bottom]
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0) # no source for the Laplace equation
a = dot(grad(u), grad(v))*dx
L = f*v*dx
u = Function(V)
solve(a == L, u, bcs)
# solve(a == L, u, bcs, annotate=False) # works even with from dolfin_adjoint import *
plt.figure()
plot(u)
plot(mesh)
It may be of interest that the same behaviour occurs for the Stokes BC Control example from the dolfin-adjoint website (http://www.dolfin-adjoint.org/en/latest/documentation/stokes-bc-control/stokes-bc-control.html):
from dolfin import *
# from dolfin_adjoint import * # With this commented the following will run.
set_log_level(LogLevel.ERROR)
from mshr import *
rect = Rectangle(Point(0, 0), Point(30, 10))
circ = Circle(Point(10, 5), 2.5)
domain = rect - circ
N = 50
my_mesh = generate_mesh(domain, N)
V_h = VectorElement("CG", my_mesh.ufl_cell(), 2)
Q_h = FiniteElement("CG", my_mesh.ufl_cell(), 1)
W = FunctionSpace(my_mesh, V_h * Q_h)
V, Q = W.split()
v, q = TestFunctions(W)
x = TrialFunction(W)
u, p = split(x)
s = Function(W, name="State")
V_collapse = V.collapse()
g = Function(V_collapse, name="Control")
nu = Constant(1)
# Define the circle boundary
class Circle(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (x[0]-10)**2 + (x[1]-5)**2 < 3**2
facet_marker = MeshFunction("size_t", my_mesh, my_mesh.topology().dim() - 1)
facet_marker.set_all(10)
Circle().mark(facet_marker, 2)
# Define a boundary measure with circle boundary tagged.
ds = ds(subdomain_data=facet_marker)
# Define boundary conditions
u_inflow = Expression(("x[1]*(10-x[1])/25", "0"), degree=1)
noslip = DirichletBC(W.sub(0), (0, 0),
"on_boundary && (x[1] >= 9.9 || x[1] < 0.1)")
inflow = DirichletBC(W.sub(0), u_inflow, "on_boundary && x[0] <= 0.1")
circle = DirichletBC(W.sub(0), g, facet_marker, 2)
bcs = [inflow, noslip, circle]
a = (nu*inner(grad(u), grad(v))*dx
- inner(p, div(v))*dx
- inner(q, div(u))*dx
)
L = inner(Constant((0, 0)), v)*dx
A, b = assemble_system(a, L, bcs)
solve(A, s.vector(), b)
u, p = split(s)
plot(u)
plot(p)