DirichletBC Function as control, taylor test fails

Hi, I am solving an optimization problem on the unit square, with the control being a scalar function V \in \mathbb{R}. For this problem, V needs to be constant in the y direction: V= V(x). (This can be interpreted as the thickness of a sheet of material, which varies only with x.) It has become apparent I need to discard the redundant degrees of freedom. To enforce this, I tried solving a simple PDE (prior to the main PDE) of the form \int \frac{\partial{u}}{\partial{y}}v=0, with boundary conditions specified only on the bottom portion of the domain y=0. (This seemed simpler than trying to enforce a variational equality.) I cannot get the taylor test to pass. Iā€™d be grateful for any suggestions on how to get this to work or alternative ways to enforce this type of constraint, thanks!

from dolfin import *
from dolfin_adjoint import *

mesh = UnitSquareMesh(50, 50)

V = FunctionSpace(mesh, "CG", 2)
bd_thickness = Function(V, name="boundary_thickness")
bd_thickness.assign(interpolate(Expression('x[0]', degree=1), V))

bd = CompiledSubDomain("near(x[1],0) && on_boundary")
bc = DirichletBC(V, bd_thickness, bd)

u_ = TrialFunction(V)
v = TestFunction(V)
j_hat = Constant(('0', '1'), name="j_hat")
a = inner(inner(grad(u_), j_hat), v)*dx
L = v*Constant(0, name="zero")*dx 

u = Function(V, name="u")
solve(a == L, u, bc)

This seems to be correct: u varies with x only, and is described only by the values of bd_thickness along the bottom boundary y=0.

J = assemble(u*dx)
dJdnu = compute_gradient(J, Control(bd_thickness))

The gradient appears constant along y.

perturbation = interpolate(Expression("0.01*sin(x[0]*pi)", pi=pi, degree=2), V)
Jhat = ReducedFunctional(J, Control(bd_thickness))
conv_rate = taylor_test(Jhat, bd_thickness, perturbation)
tape = get_working_tape()

The taylor test fails. However, the tape looks correct.

This is expected as it is a linear problem. Thus, higher order derivatives are zero and so the residual should be zero. If you look at the taylor test the residuals are machine precision 0.

If you change your functional to u**2*dx you should see convergence rate of 2 and non-zero residuals.

1 Like

That seems so obvious in retrospect - thank you Sebastian!