# Third order residual not obtained

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 *
import numpy as np

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

output = Vector(MPI.comm_world, self.V_sph_dim)

return output

def evaluate_tlm_component(self, inputs, tlm_inputs, block_variable, idx, prepared=None):

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

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

# 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

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