How to get adjoint solution from compute_gradient()

Hello, I would like to be able to access the adjoint variable (solution of the adjoint system).

I calculate the gradient \frac{dJ}{dm} by using the command compute_gradient(J, m). I’m assuming that the gradient calculation is done using the adjoint approach, so following the notation in the dolfin-adjoint mathematical background of adjoints, the adjoint equation is solved (\frac{\partial F(u,m)}{\partial u})^* \lambda=\frac{\partial J}{\partial u}^* and the result is used to calculate \frac{dJ}{dm}^* = -\frac{\partial F}{\partial m}^* \lambda + \frac{\partial J}{\partial m}^*. I’d like to be able to get the variable \lambda.

I’ve tried using J.block_variable.adj_value, but I’m not sure that I’m using it correctly because I get back a value of 1, when it should be a vector.

In this simple example, I’m using dolfin-adjoint 2019.1.0 and solving a linear elasticity pde.

from fenics import *
from fenics_adjoint import *
import numpy as np
from ufl import nabla_div

# Define parameters
nelx = 30
nely = 15
nelz = 15
lx = float(nelx)
ly = float(nely)
lz = float(nelz)
E = 1
v = 0.3
mu = Constant(E / (2*(1+v)))
lmbda = Constant(E*v / ((1+v)*(1-2*v)))

# Create a mesh
mesh = BoxMesh.create(
    [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)], # define opposing corners
    [nelx, nely, nelz], # number of elements in each direction
mesh = Mesh(mesh)

# FE function space and functions
V = VectorFunctionSpace(mesh, "CG", 1)
D = FunctionSpace(mesh, "DG", 0)
u = TrialFunction(V)
v = TestFunction(V)
u_sol = Function(V) # function for storing the solution
density = Function(D) # density field

# Bounday condition
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)
left = Left()
bc = DirichletBC(V, Constant((0., 0., 0.)), left)

# Stress and strain functions
def sigma(u):
    n = len(u) # size of u
    return lmbda*nabla_div(u)*Identity(n) + 2*mu*epsilon(u)

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

# Define density field
density.vector()[:] = np.ones(nelx*nely*nelz)

# Define forms
f = Constant((0., 0., -1.))
L = dot(f, v) * dx
a = inner(density*sigma(u), nabla_grad(v)) * dx

# Solve
A, b = assemble_system(a, L, bc)
solver = KrylovSolver("cg", "ilu")
solver.solve(A, u_sol.vector(), b)

# Get derivative
J = assemble(action(L, u_sol)) # define compliance
m = Control(density) # define density as the design variable
dJ_dm = compute_gradient(J, m)

# Get adjoint solution
adjoint_sol = J.block_variable.adj_value
print('adjoint value:', adjoint_sol)

As each block has an adjoint (not just the solve block), you need to access the correct one.
I think you want:


Thank you so much. I do think that is the adjoint I want, however, I am a bit surprised by the result.

The function J that I’m using is actually self-adjoint, so if you calculate \frac{dJ}{dm} using the adjoint approach you’d find that \lambda=-u where u is the solution to the pde. However, the values that I get from u_sol.block_variable.adj_value[:] are quite different than the solution values u_sol.vector()[:]. Do you have any idea why this discrepancy would be occuring?

Shouldn’t you actually get that \lambda=u
as J=\int_\Omega f \cdot u \mathrm{d} x , thus \frac{\partial J}{\partial u}[v]=\int_{\Omega} f\cdot v~\mathrm{dx} which is equal to the LHS of your original problem?
To get this, you can add the following code:

tape = get_working_tape()
solve_block = tape.get_blocks()[1]

Yes, you’re right, it should be that \lambda = u. The negative sign was my mistake.

Thank you so much for your help, that last bit of code did the trick.