BUG: Using a 0 RHS doesn't work for `solve(a==L)` when `a` is a Linear form

[jonnybolton16]

Describe the bug
When solving a linear problem with 0 RHS, using L = 0 works with the solve method of LinearVariationalSolver. However, it does not work with solve function. In this case, L = fd.Form([]) (or some other zero-form) is required.

Steps to Reproduce

import firedrake as fd
import matplotlib.pyplot as plt

n = 10
mesh = fd.UnitSquareMesh(n, n)
x, y = fd.SpatialCoordinate(mesh)

V = fd.FunctionSpace(mesh, 'CG', 1)
u = fd.TrialFunction(V)
v = fd.TestFunction(V)
u_num = fd.Function(V)

a = fd.dot(fd.grad(v), fd.grad(u)) * fd.dx

bcs = [
    fd.DirichletBC(V,  3*fd.sin(fd.pi*y),     1 ),
    fd.DirichletBC(V,    fd.sin(fd.pi*y),     2 ),
    fd.DirichletBC(V, -2*fd.sin(fd.pi*x), (3, 4))
]

L = 0
fd.LinearVariationalSolver(
    fd.LinearVariationalProblem(a, L, u_num, bcs=bcs)
).solve()
fd.trisurf(u_num)
plt.title('Solution from solver.solve()')
plt.show()

L = fd.Form([])
fd.solve(a==L, u_num, bcs=bcs)
fd.trisurf(u_num)
plt.title('Solution from solve(a==fd.Form([]), ...)')
plt.show()

L = 0
fd.solve(a==L, u_num, bcs=bcs)

Expected behavior
I would like solve(a == 0) to work, the same way it does for a LinearVariationalSolver

Error message

Traceback (most recent call last):
  File "/mnt/c/Users/jonny/OneDrive - University of Leeds/Desktop/waveEnergyPhD/temp.py", line 36, in <module>
    fd.solve(a==L, u_num, bcs=bcs)
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jonny/firedrake/firedrake_venv/lib/python3.12/site-packages/firedrake/adjoint_utils/solving.py", line 57, in wrapper
    output = solve(*args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jonny/firedrake/firedrake_venv/lib/python3.12/site-packages/firedrake/solving.py", line 143, in solve
    _solve_varproblem(*args, **kwargs)
  File "/home/jonny/firedrake/firedrake_venv/lib/python3.12/site-packages/firedrake/solving.py", line 177, in _solve_varproblem
    problem = vs.NonlinearVariationalProblem(eq.lhs, u, bcs, J, Jp,
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/usr/lib/python3.12/contextlib.py", line 81, in inner
    return func(*args, **kwds)
           ^^^^^^^^^^^^^^^^^^^
  File "/home/jonny/firedrake/firedrake_venv/lib/python3.12/site-packages/firedrake/adjoint_utils/variational_solver.py", line 15, in wrapper
    init(self, *args, **kwargs)
  File "/home/jonny/firedrake/firedrake_venv/lib/python3.12/site-packages/firedrake/variational_solver.py", line 121, in __init__
    check_pde_args(self.F, self.J, self.Jp)
  File "/home/jonny/firedrake/firedrake_venv/lib/python3.12/site-packages/firedrake/variational_solver.py", line 29, in check_pde_args
    raise ValueError("Provided residual is not a linear form")
ValueError: Provided residual is not a linear form

Environment:

• OS: WSL
• Python version: 3.12.3
• Output of pip list

Package              Version
-------------------- -----------
cachetools           6.2.4
certifi              2025.11.12
cgen                 2025.1
charset-normalizer   3.4.4
checkpoint_schedules 1.0.4
codepy               2025.1
colorama             0.4.6
constantdict         2025.3
contourpy            1.3.3
cycler               0.12.1
decorator            4.4.2
fenics-ufl           2025.2.0
firedrake            2025.10.2
firedrake-fiat       2025.10.1
fonttools            4.61.1
genpy                2022.1
h5py                 3.15.1
idna                 3.11
immutabledict        4.2.2
iniconfig            2.3.0
islpy                2025.2.5
kiwisolver           1.4.9
libsupermesh         2025.4
loopy                2025.2
Mako                 1.3.10
MarkupSafe           3.0.3
matplotlib           3.10.8
matspy               1.0.0
mpi-pytest           2025.7
mpi4py               4.1.1
mpmath               1.3.0
numpy                2.4.0
packaging            25.0
pandas               2.3.3
petsc4py             3.24.0
petsctools           2025.3
pillow               12.0.0
pip                  24.0
pkgconfig            1.5.5
platformdirs         4.5.1
pluggy               1.6.0
progress             1.6.1
pyadjoint-ad         2025.10.1
pycparser            2.23
Pygments             2.19.2
pymbolic             2025.1
pyparsing            3.3.1
pytest               9.0.2
python-dateutil      2.9.0.post0
pytools              2025.2.5
pytz                 2025.2
recursivenodes       0.2.0
requests             2.32.5
rtree                1.4.1
scipy                1.16.3
siphash24            1.8
six                  1.17.0
symengine            0.14.1
sympy                1.14.0
typing_extensions    4.15.0
tzdata               2025.3
urllib3              2.6.2
vtk                  9.5.2

[jonnybolton16]

It's interpretting a==0 as a nonlinear problem, but then expects to see a Function not a TrialFunction

[dham]

This is a longstanding limitation of limitation of the solve interface. LinearVariationalSolver is preferred, but if you want to use solve then you need a zero form on the RHS: Constant(0)*v*dx.

[jonnybolton16]

Is there a way to tell solving._solve_varproblem to go and create a LinearVariationalProblem when eq.lhs is a Linear form and eq.rhs is 0?

[dham]

Not currently. You could modify solve to check for the case where the LHS is a bilinear form and the RHS is zero and replace the RHS in that case by the one above.