Miscellaneous notes¶
Naming Functions
and Constants
¶
Consider the example presented in the tutorial again:
from dolfin import *
from dolfin_adjoint import *
n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)
ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])"), degree=2), V)
u = ic.copy(deepcopy=True)
u_next = Function(V)
v = TestFunction(V)
nu = Constant(0.0001)
timestep = Constant(0.01)
F = (inner((u_next - u)/timestep, v)
+ inner(grad(u_next)*u_next, v)
+ nu*inner(grad(u_next), grad(v)))*dx
bc = DirichletBC(V, (0.0, 0.0), "on_boundary")
t = 0.0
end = 0.1
while (t <= end):
solve(F == 0, u_next, bc)
u.assign(u_next)
t += float(timestep)
J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
dJdic = compute_gradient(J, Control(u), forget=False)
dJdnu = compute_gradient(J, Control(nu))
Note that the Constant nu
passed to Control(nu)
had to be exactly the same variable
as was used in the forward form. This could be quite inconvenient, if the form creation occurs in one file, and the
adjoint is driven from another. To facilitate such cases, it is possible to give the Constant
a name
:
nu = Constant(0.0001, name="nu")
which may then be used to drive the adjoint. However, the Control
class does not have
enough information to determine what kind of control it is: therefore in this case the ConstantControl
class must be used:
dJdnu = compute_gradient(J, ConstantControl("nu"))
A full example is given in the following code.
from dolfin import *
from dolfin_adjoint import *
n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)
ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])"), degree=2), V)
def main(nu):
u = ic.copy(deepcopy=True)
u_next = Function(V)
v = TestFunction(V)
timestep = Constant(0.01)
F = (inner((u_next - u)/timestep, v)
+ inner(grad(u_next)*u_next, v)
+ nu*inner(grad(u_next), grad(v)))*dx
bc = DirichletBC(V, (0.0, 0.0), "on_boundary")
t = 0.0
end = 0.1
while (t <= end):
solve(F == 0, u_next, bc)
u.assign(u_next)
t += float(timestep)
return u
if __name__ == "__main__":
nu = Constant(0.0001, name="nu")
u = main(nu)
J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
dJdnu = compute_gradient(J, ConstantControl("nu"))
Jnu = assemble(inner(u, u)*dx) # current value
parameters["adjoint"]["stop_annotating"] = True # stop registering equations
def Jhat(nu): # the functional as a pure function of nu
u = main(nu)
return assemble(inner(u, u)*dx)
conv_rate = taylor_test(Jhat, ConstantControl("nu"), Jnu, dJdnu)
Similarly, it is possible to give Functions
names:
u = Function(ic, name="Velocity")
u_next = Function(ic, name="VelocityNext")
which can then be used in other places with FunctionControl
:
dJdic = compute_gradient(J, FunctionControl("Velocity"))
A full example is given in the following code.
from dolfin import *
from dolfin_adjoint import *
n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)
ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])"), degree=2), V)
nu = Constant(0.0001, name="nu")
def main(ic):
u = ic.copy(deepcopy=True, name="Velocity")
u_next = Function(V, name="VelocityNext")
v = TestFunction(V)
timestep = Constant(0.01)
F = (inner((u_next - u)/timestep, v)
+ inner(grad(u_next)*u_next, v)
+ nu*inner(grad(u_next), grad(v)))*dx
bc = DirichletBC(V, (0.0, 0.0), "on_boundary")
t = 0.0
end = 0.1
while (t <= end):
solve(F == 0, u_next, bc)
u.assign(u_next)
t += float(timestep)
return u
if __name__ == "__main__":
u = main(ic)
J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
dJdic = compute_gradient(J, FunctionControl("Velocity"))
Jic = assemble(inner(u, u)*dx) # current value
parameters["adjoint"]["stop_annotating"] = True # stop registering equations
def Jhat(ic):
u = main(ic)
return assemble(inner(u, u)*dx)
conv_rate = taylor_test(Jhat, FunctionControl("Velocity"), Jic, dJdic, value=ic)
Lower-level interfaces¶
A lower-level interface is available to loop over the adjoint solutions (backwards in time) for some manual calculation.
This uses the compute_adjoint
function:
J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
for (variable, solution) in compute_adjoint(J):
...
Here, solution
is a dolfin.Function
storing a particular adjoint solution,,
and variable
is a libadjoint.Variable
that describes
which forward solution it corresponds to.
Similarly, it is possible to loop over the tangent linear solutions with the compute_tlm
function:
param = Control(nu)
for (variable, solution) in compute_tlm(param):
...