Irksome breaks with Cofunctions

[rckirby]

It seems we are missing something in our UFL manipulation. The example below fails with

Traceback (most recent call last):
  File "/Users/robert_kirby/fddelta.py", line 22, in <module>
    stepper = TimeStepper(F, bt, t, dt, u)
  File "/Users/robert_kirby/Code/irksome/irksome/stepper.py", line 85, in TimeStepper
    return StageDerivativeTimeStepper(
        F, butcher_tableau, t, dt, u0, bcs, appctx=appctx,
        solver_parameters=solver_parameters, nullspace=nullspace,
        bc_type=bc_type, splitting=splitting)
  File "/Users/robert_kirby/Code/irksome/irksome/stage_derivative.py", line 180, in __init__
    super().__init__(F, t, dt, u0,
    ~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^
                     butcher_tableau.num_stages, bcs=bcs,
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    ...<2 lines>...
                     splitting=splitting, bc_type=bc_type,
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                     butcher_tableau=butcher_tableau)
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/robert_kirby/Code/irksome/irksome/base_time_stepper.py", line 125, in __init__
    self.prob = NonlinearVariationalProblem(Fbig, stages, bigBCs)
                ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/opt/homebrew/Cellar/python@3.13/3.13.2/Frameworks/Python.framework/Versions/3.13/lib/python3.13/contextlib.py", line 85, in inner
    return func(*args, **kwds)
  File "/Users/robert_kirby/Code/firedrake/firedrake/adjoint_utils/variational_solver.py", line 15, in wrapper
    init(self, *args, **kwargs)
    ~~~~^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/robert_kirby/Code/firedrake/firedrake/variational_solver.py", line 116, in __init__
    check_pde_args(self.F, self.J, self.Jp)
    ~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/robert_kirby/Code/firedrake/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

from firedrake import *
from irksome import Dt, RadauIIA, TimeStepper

N = 4
msh = UnitSquareMesh(N, N)
V = FunctionSpace(msh, "CG", 1)
u = Function(V)
v = TestFunction(V)

vom_msh = VertexOnlyMesh(msh, [(0.5, 0.5)])
Vvom = FunctionSpace(vom_msh, "DG", 0)
delta = Function(Vvom).interpolate(1.0)
vvom = TestFunction(Vvom)
d = Cofunction(V.dual()).interpolate(assemble(inner(delta, vvom) * dx))

F = inner(Dt(u), v) * dx + inner(grad(u), grad(v)) * dx - d(v)

t = Constant(0)
dt = Constant(1/N)
bt = RadauIIA(2)

stepper = TimeStepper(F, bt, t, dt, u)

stepper.advance()

[rckirby]

We'll work on this with #142. I tried something there that doesn't work.

[pbrubeck]

The example above would just work (as soon as firedrakeproject/firedrake#4197 and FEniCS/ufl#370 land) if you replace d(v) by Interpolate(v, d).