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):
    ...