firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
484 stars 157 forks source link

Incorrect gradients with DirichletBC with overlapping boundary #2557

Open IvanYashchuk opened 1 year ago

IvanYashchuk commented 1 year ago

I noticed that the gradient wrt DirichletBC might be incorrectly non-zero even when it doesn't have an influence on the result:

from firedrake import *
import firedrake_adjoint

n = 30
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "P", 1)

x = SpatialCoordinate(mesh)
f = x[0]

u = Function(V)
v = TestFunction(V)

c0 = Constant(0.0001)
c0.adj_value = 1.0
c1 = Constant(0.0002)
c1.adj_value = 1.0

# Note the boundary is overlapping so the second bc would overwrite the first one
bcs = [DirichletBC(V, c0, "on_boundary"), DirichletBC(V, c1, "on_boundary")]

JJ = 0.5 * inner(grad(u), grad(u)) * dx - f * u * dx
F = derivative(JJ, u, v)
solve(F == 0, u, bcs=bcs)

J = assemble(inner(u, u)*dx)
dJdc0 = firedrake_adjoint.compute_gradient(J, c0)
dJdc1 = firedrake_adjoint.compute_gradient(J, c1)

h = Constant(0.0001)  # the direction of the perturbation
for c in [c0, c1]:
    Jhat = firedrake_adjoint.ReducedFunctional(J, firedrake_adjoint.Control(c))  # the functional as a pure function of nu
    conv_rate = firedrake_adjoint.taylor_test(Jhat, c, h)
Running Taylor test
Computed residuals: [3.54177434863098e-08, 1.77088717431549e-08, 8.85443587157745e-09, 4.427217935788725e-09]
Computed convergence rates: [1.0, 1.0, 1.0]
Running Taylor test
Computed residuals: [1.000001825575511e-12, 2.500009679968148e-13, 6.250116187450287e-14, 1.5625628433530902e-14]
Computed convergence rates: [1.9999970476603297, 1.9999787666459534, 1.999968795778822]

Taylor test convergence rate wrt c0 is 1.0 indicating incorrect implementation.

wence- commented 1 year ago

Not sure what the right answer is here, the problem doesn't see c0

IvanYashchuk commented 1 year ago

The problem doesn't see c0 and so the gradient wrt c0 should be zero, right? If I'd use finite differences to compute the derivative I would get zero.

dham commented 1 year ago

I'd argue that the input is illegal here and so you deserve whatever Firedrake-adjoint gives you ;).

It's not clear to me how we would detect this situation and do something different.