Incorrect boundary conditions for linear variational problems whose RHS depends on the solution

[jrmaddison]

For example:

solve(inner(trial, test) * dx
      == inner(w, test) * dx + inner(grad(w), grad(test)) * dx,
      w, bcs)

seems to be equivalent to

for bc in bcs:
    bc.apply(w)
solve(inner(trial, test) * dx
      == inner(w, test) * dx + inner(grad(w), grad(test)) * dx,
      w, bcs)

This then leads to unexpected behavior e.g. in the following example

from firedrake import *

mesh = UnitSquareMesh(10, 10)
X = SpatialCoordinate(mesh)
space = FunctionSpace(mesh, "Lagrange", 1)
test, trial = TestFunction(space), TrialFunction(space)

bcs = [DirichletBC(space, 0.0, "on_boundary")]

u = Function(space, name="u")
u.interpolate(X[0])

v = Function(space, name="v")
solve(inner(trial, test) * dx
      == inner(u, test) * dx + inner(grad(u), grad(test)) * dx,
      v, bcs)

w = Function(space, name="w")
w.assign(u)

solve(inner(trial, test) * dx
      == inner(w, test) * dx + inner(grad(w), grad(test)) * dx,
      w, bcs)

print(sqrt(abs(assemble(inner(v - w, v - w) * dx))))

[jrmaddison]

It's not clear to me if this is a bug or incorrect usage, but perhaps an exception could be raised if this usage is invalid?

[ksagiyam]

The solution w cannot appear in lhs or rhs in linear problems (in Firedrake sense); I think you are right that an exception could be raised.

In the first solve you are solving a linear system $(v, test) = (u, test) + (\nabla u, \nabla test)$ for $v$ given $u$, and if you intend to solve this linear system, you always get the form used in the first solve (and never the one used in the second solve).