Hi, I’m trying to implement the interpolation of a coefficient as part of an adjoint optimisation problem. Say, for example a coefficient z that is a function of two variables, x_1 and x_2
where these two variables are control inputs and the coefficient is to be used in the bigger problem. There is no analytical form of the coefficient calculation.
Evaluation of the function at a point does not seems to be annotated for the adjoint derivation. So this does not work:
J = z(x1,x2)
I tried to write a custom function and overload it to work with the adjoint, as in the MWE below. However, the Taylor test returns a convergence around 1 instead of 2, so this seems to be incorrect.
How would I be able to fix the code below to yield the correct behaviour? Or is there perhaps a more effective method to include the point evaluation of a function in the adjoint?
from fenics import * from fenics_adjoint import * import numpy as np from pyadjoint import Block from pyadjoint.overloaded_function import overload_function def get_coefficient(func, coord1, coord2): return func(coord1, coord2) backend_get_coefficient = get_coefficient class CoefficientBlock(Block): def __init__(self, func, coord1, coord2, **kwargs): super(CoefficientBlock, self).__init__() self.kwargs = kwargs self.func = func self.add_dependency(coord1) self.add_dependency(coord2) def __str__(self): return "CoefficientBlock" def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, prepared=None): grad_idx = project(self.func.dx(idx), self.func.function_space()) return grad_idx(inputs, inputs) * adj_inputs def recompute_component(self, inputs, block_variable, idx, prepared): return backend_get_coefficient(self.func, inputs, inputs) get_coefficient = overload_function(get_coefficient, CoefficientBlock) mesh = UnitSquareMesh(10, 10) V1 = FunctionSpace(mesh, "Lagrange", 1) u = Function(V1) x = SpatialCoordinate(u) z = project(x*x, V1) x1 = Constant(np.random.rand()) x2 = Constant(np.random.rand()) J = get_coefficient(z, x1, x2) controls = [x1, x2] m = [Control(c) for c in controls] h = [Constant(0.01*np.random.rand()) for c in controls] Jhat = ReducedFunctional(J, m) Jhat.derivative() taylor_test(Jhat, controls, h)
Results from Taylor test:
Computed residuals: [7.3138620463933775e-06, 3.6569310231966887e-06, 1.8284655115428332e-06, 9.142327556603943e-07] Computed convergence rates: [1.0, 1.0000000000437994, 1.0000000001751974]
Any help is much appreciated!