Third order residual not obtained

dolfin-adjoint
1

Hi, in verifying a very simple “shape optimization” code I have written with dolfin-adjoint, I was performing a test with Hessians and I noted I am not obtaining order 3 residuals. I am posting an example. The overloaded block is just a matrix-vector product (so, 0 Hessian), that enables a transfer from functions on a unit interval to a square.

from dolfin import *
from dolfin_adjoint import *
import numpy as np
from pyadjoint import Block
from pyadjoint.overloaded_function import overload_function

# Overloading

def backend_displacement(g, M2, VD):
    gv = g.vector()[:]
    W = Function(VD)
    W.vector()[:] = M2 @ gv
    return W

class DisplacementBlock(Block):
    def __init__(self, g, M2, VD, **kwargs):
        super(DisplacementBlock, self).__init__()
        self.kwargs = kwargs
        self.add_dependency(g)

        self.M2 = M2  # matrix bringing us from g.function_space() to VD
        self.VD = VD  # the volumetric function space, vector fields version
        self.VD_dim = VD.dim()  # the volumetric function space, vector fields version
        self.V_sph_dim = g.function_space().dim()

    def __str__(self):
        return 'DisplacementBlock'

    def recompute_component(self, inputs, block_variable, idx, prepared):
        return backend_displacement(inputs[0], self.M2, self.VD)

    def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None):
        adj_input = adj_inputs[0]  # this is a vector, with which we test the adjoint jacobian
        adj_action = self.M2.T @ adj_input
        output = Vector(MPI.comm_world, self.V_sph_dim)
        output[:] = adj_action

        return output

    def evaluate_tlm_component(self, inputs, tlm_inputs, block_variable, idx, prepared=None):
        # see http://www.dolfin-adjoint.org/en/latest/documentation/pyadjoint_docs.html

        return self.recompute_component(tlm_inputs, block_variable, idx, prepared)

    def evaluate_hessian_component(self, inputs, hessian_inputs, adj_inputs, block_variable, idx,
                                   relevant_dependencies, prepared=None):
        output = Vector(MPI.comm_world, self.V_sph_dim)
        output[:] = 0
        return output



radial_displacement = overload_function(backend_displacement, DisplacementBlock)

domain = UnitSquareMesh(5,5)
interval = UnitIntervalMesh(10)

V_int = FunctionSpace(interval, "CG", 1)
V_def = VectorFunctionSpace(domain, "CG", 1)

M2 = np.random.rand(V_def.dim(), V_int.dim())

q = Function(V_int)  # null displacement (but doesn't need to be so)
W = radial_displacement(q, M2, V_def)
ALE.move(domain, W)

# PDE

V = FunctionSpace(domain, "CG", 1)
u = Function(V)

w = TrialFunction(V)
v = TestFunction(V)
bcs = DirichletBC(V, 0, "on_boundary")

dx = Measure("dx", domain=domain)

solve(inner(grad(w), grad(v)) * dx == v * dx, u, bcs)

# Cost functional

J = assemble(sin(u) * dx)

Jhat = ReducedFunctional(J, [Control(q)])
d = Function(V_int)
d.vector()[:] = (np.random.rand(len(d.vector())) - .5) * 2 * 1e-2
td = taylor_to_dict(Jhat, q, d)

print(td['R2']['Rate'])
print(td['R2'])

The result for the Hessian residual:

{'Residual': [6.073355259626869e-09, 1.5186571271107686e-09, 3.7970434087315546e-10, 9.493119657780799e-11], 'Rate': [1.9996975778376476, 1.9998477865513609, 1.9999223191972153]}

which is not of third order as expected. The deformation that is being performed only minimally perturbs the mesh, so, no weird compenetrations are happening. What could be wrong then? Thanks a lot in advance

2

The problem is that the evaluate_hessian_component method should not return 0. This is not documented very well (or at all, really) in the pyadjoint docs.

There is a nonzero contribution from the derivative wrt the adjoint value.
Have a look here for the formula used in the hessian computations in pyadjoint (the notation can be a bit difficult but v^{(2)}_{(1)i} corresponds to the hessian input/output): The Art of Differentiating Computer Programs: An Introduction to Algorithmic ... - Uwe Naumann - Google Bøker

Since the second order derivatives are 0 the hessian becomes the same as the adjoint just on the hessian inputs:

    def evaluate_hessian_component(self, inputs, hessian_inputs, adj_inputs, block_variable, idx,
                                   relevant_dependencies, prepared=None):
        return self.evaluate_adj_component(inputs, hessian_inputs, block_variable, idx)
1 Like