Parallel

Applying algorithmic differentiation tools to parallel source code is still a major research area, and most adjoint codes that work in parallel manually adjoin the parallel communication sections of their code.

One of the major advantages of the new high-level abstraction used in dolfin-adjoint is that the problem of parallelism in adjoint codes simply disappears: indeed, there is not a single line of parallel-specific code in dolfin-adjoint or libadjoint. For more details on how this works, see the papers.

Therefore, if your forward model runs in parallel, your adjoint will also, with no modification. For example, let us take the checkpointed adjoint model used in the previous section:

$ mpiexec -n 8 python tutorial5.py
...
Process 0: Convergence orders for Taylor remainder with adjoint information (should all be 2):
  [1.9744066553464978, 1.9872606129796675, 1.9936586367818951, 1.9968385300177882]

Similarly, parallelism over OpenMP works in the same way:

from dolfin import *
from dolfin_adjoint import *

parameters["num_threads"] = 8

adj_checkpointing(strategy='multistage', steps=11,
                  snaps_on_disk=2, snaps_in_ram=2, verbose=True)

n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)

ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])")),  V)

def main(nu):
    u = Function(ic)
    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)
        adj_inc_timestep()

    return u

if __name__ == "__main__":
    nu = Constant(0.0001)
    u = main(nu)

    J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
    dJdnu = compute_gradient(J, Control(nu))

    Jnu = assemble(inner(u, u)*dx)

    parameters["adjoint"]["stop_annotating"] = True
    def Jhat(nu):
        u = main(nu)
        return assemble(inner(u, u)*dx)

    conv_rate = taylor_test(Jhat, Control(nu), Jnu, dJdnu)

Download the checkpointed threaded parallel code.

$ python tutorial6.py
...
Convergence orders for Taylor remainder with adjoint information (should all be 2):
  [1.9581779061701046, 1.9787032981536112, 1.989252750128478, 1.994601330484473]