From 392eb70d4332aa78ac466b67bf6a256def496d1b Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 17:32:53 +0000
Subject: [PATCH 01/13] #223: Remove form-related methods in MeshSeq

---
 goalie/mesh_seq.py | 39 ++-------------------------------------
 1 file changed, 2 insertions(+), 37 deletions(-)

diff --git a/goalie/mesh_seq.py b/goalie/mesh_seq.py
index 7e10c2b2..ed849fb7 100644
--- a/goalie/mesh_seq.py
+++ b/goalie/mesh_seq.py
@@ -39,7 +39,6 @@ def __init__(self, time_partition, initial_meshes, **kwargs):
             :meth:`~.MeshSeq.get_function_spaces`
         :kwarg get_initial_condition: a function as described in
             :meth:`~.MeshSeq.get_initial_condition`
-        :kwarg get_form: a function as described in :meth:`~.MeshSeq.get_form`
         :kwarg get_solver: a function as described in :meth:`~.MeshSeq.get_solver`
         :kwarg transfer_method: the method to use for transferring fields between
             meshes. Options are "project" (default) and "interpolate". See
@@ -58,7 +57,6 @@ def __init__(self, time_partition, initial_meshes, **kwargs):
         self._fs = None
         self._get_function_spaces = kwargs.get("get_function_spaces")
         self._get_initial_condition = kwargs.get("get_initial_condition")
-        self._get_form = kwargs.get("get_form")
         self._get_solver = kwargs.get("get_solver")
         self._transfer_method = kwargs.get("transfer_method", "project")
         self._transfer_kwargs = kwargs.get("transfer_kwargs", {})
@@ -267,29 +265,6 @@ def get_initial_condition(self):
             for field, fs in self.function_spaces.items()
         }
 
-    def get_form(self):
-        """
-        Get the function mapping a subinterval index and a solution dictionary to a
-        dictionary containing parts of the PDE weak form corresponding to each solution
-        component.
-
-        Signature for the function to be returned:
-        ```
-        :arg index: the subinterval index
-        :type index: :class:`int`
-        :arg solutions: map from fields to tuples of current and previous solutions
-        :type solutions: :class:`dict` with :class:`str` keys and :class:`tuple` values
-        :return: map from fields to the corresponding forms
-        :rtype: :class:`dict` with :class:`str` keys and :class:`ufl.form.Form` values
-        ```
-
-        :returns: the function for obtaining the form
-        :rtype: see docstring above
-        """
-        if self._get_form is None:
-            raise NotImplementedError("'get_form' needs implementing.")
-        return self._get_form(self)
-
     def get_solver(self):
         """
         Get the function mapping a subinterval index and an initial condition dictionary
@@ -339,10 +314,10 @@ def _transfer(self, source, target_space, **kwargs):
 
     def _outputs_consistent(self):
         """
-        Assert that function spaces, initial conditions, and forms are given in a
+        Assert that function spaces and initial conditions are given in a
         dictionary format with :attr:`MeshSeq.fields` as keys.
         """
-        for method in ["function_spaces", "initial_condition", "form", "solver"]:
+        for method in ["function_spaces", "initial_condition", "solver"]:
             if getattr(self, f"_get_{method}") is None:
                 continue
             method_map = getattr(self, f"get_{method}")
@@ -350,9 +325,6 @@ def _outputs_consistent(self):
                 method_map = method_map(self.meshes[0])
             elif method == "initial_condition":
                 method_map = method_map()
-            elif method == "form":
-                self._reinitialise_fields(self.get_initial_condition())
-                method_map = method_map()(0)
             elif method == "solver":
                 self._reinitialise_fields(self.get_initial_condition())
                 solver_gen = method_map()(0)
@@ -435,13 +407,6 @@ def initial_condition(self):
         """
         return AttrDict(self.get_initial_condition())
 
-    @property
-    def form(self):
-        """
-        See :meth:`~.MeshSeq.get_form`.
-        """
-        return self.get_form()
-
     @property
     def solver(self):
         """

From 1a1214d9f1b71bf19cbf7a861ef4ecb1faacbf55 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 17:34:03 +0000
Subject: [PATCH 02/13] #223: Update demos without GoalOrientedMeshSeq

---
 demos/burgers-hessian.py           | 16 +----
 demos/burgers.py                   | 62 ++++++--------------
 demos/burgers1.py                  | 18 +-----
 demos/burgers2.py                  | 16 +----
 demos/gray_scott.py                | 29 +++-------
 demos/gray_scott_split.py          | 29 ++--------
 demos/ode.py                       | 93 +++++++++++++++---------------
 demos/point_discharge2d-hessian.py | 19 +-----
 demos/solid_body_rotation.py       | 30 +++-------
 9 files changed, 95 insertions(+), 217 deletions(-)

diff --git a/demos/burgers-hessian.py b/demos/burgers-hessian.py
index ebe343c3..13429722 100644
--- a/demos/burgers-hessian.py
+++ b/demos/burgers-hessian.py
@@ -23,8 +23,8 @@ def get_function_spaces(mesh):
     return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         u, u_ = mesh_seq.fields["u"]
 
         # Define constants
@@ -39,17 +39,6 @@ def form(index):
             + inner(dot(u, nabla_grad(u)), v) * dx
             + nu * inner(grad(u), grad(v)) * dx
         )
-        return {"u": F}
-
-    return form
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        u, u_ = mesh_seq.fields["u"]
-
-        # Define form
-        F = mesh_seq.form(index)["u"]
 
         # Time integrate from t_start to t_end
         tp = mesh_seq.time_partition
@@ -91,7 +80,6 @@ def get_initial_condition(mesh_seq):
     meshes,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
 )
 
diff --git a/demos/burgers.py b/demos/burgers.py
index b48415b1..addc9576 100644
--- a/demos/burgers.py
+++ b/demos/burgers.py
@@ -40,25 +40,29 @@ def get_function_spaces(mesh):
     return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-# In order to solve PDEs using the finite element method, we
-# require a weak form. For this, Goalie requires a function
-# that maps the :class:`MeshSeq` index and a dictionary of
-# solution data to the form. The form should be
-# returned inside its own dictionary, with an entry for each equation
-# being solved for.
-#
 # The solution :class:`Function`\s are automatically built on the function spaces given
 # by the :func:`get_function_spaces` function and are accessed via the :attr:`fields`
 # attribute of the :class:`MeshSeq`. This attribute provides a dictionary of tuples
 # containing the current and lagged solutions for each field.
-# Similarly, timestepping information associated with a given subinterval
-# can be accessed via the :attr:`TimePartition` attribute of
-# the :class:`MeshSeq`. For technical reasons, we need to create a :class:`Function`
-# in the `'R'` space (of real numbers) to hold constants. ::
+#
+# In order to solve the PDE, we need to choose a time integration routine and solver
+# parameters for the underlying linear and nonlinear systems. This is achieved below by
+# using a function :func:`solver` whose input is the :class:`MeshSeq` index. The
+# function should return a generator that yields the solution at each timestep, so
+# that Goalie can efficiently track the solution history. This is done by using the
+# `yield` statement before progressing to the next timestep.
+#
+# The lagged solution is assigned the initial condition for the current subinterval
+# index. For the :math:`0^{th}` index, this will be provided by the initial conditions,
+# otherwise it will be transferred from the previous mesh in the sequence.
+# Timestepping information associated with a given subinterval can be accessed via the
+# :attr:`TimePartition` attribute of the :class:`MeshSeq`. For technical reasons, we
+# need to create a :class:`Function` in the `'R'` space (of real numbers) to hold
+# constants.::
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         # Get the current and lagged solutions
         u, u_ = mesh_seq.fields["u"]
 
@@ -74,37 +78,6 @@ def form(index):
             + inner(dot(u, nabla_grad(u)), v) * dx
             + nu * inner(grad(u), grad(v)) * dx
         )
-        return {"u": F}
-
-    return form
-
-
-# We have a weak form. The dictionary usage may seem cumbersome when applied to such a
-# simple problem, but it comes in handy when solving adjoint problems associated with
-# coupled systems of equations.
-
-# In order to solve the PDE, we need to choose
-# a time integration routine and solver parameters for the underlying
-# linear and nonlinear systems. This is achieved below by using a function
-# :func:`solver` whose input is the :class:`MeshSeq` index. As noted above, the solution
-# :class:`Function`\s are automatically initialised and accessed via the
-# :attr:`fields` attribute of the :class:`MeshSeq`. The lagged solution is assigned
-# the initial condition for the current subinterval index. For the :math:`0^{th}` index,
-# this will be provided by the initial conditions, otherwise it will be transferred
-# from the previous mesh in the sequence.
-#
-# The function should return a generator that yields the solution at each timestep, so
-# that Goalie can efficiently track the solution history. This is done by using the
-# `yield` statement before progressing to the next timestep. ::
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        # Get the current and lagged solutions
-        u, u_ = mesh_seq.fields["u"]
-
-        # Define form
-        F = mesh_seq.form(index)["u"]
 
         # Time integrate from t_start to t_end
         tp = mesh_seq.time_partition
@@ -164,7 +137,6 @@ def get_initial_condition(mesh_seq):
     meshes,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
 )
 solutions = mesh_seq.solve_forward()
diff --git a/demos/burgers1.py b/demos/burgers1.py
index c5d4020d..9710f2e7 100644
--- a/demos/burgers1.py
+++ b/demos/burgers1.py
@@ -15,7 +15,7 @@
 from goalie_adjoint import *
 
 # For ease, the list of field names and functions for obtaining the
-# function spaces, forms, solvers, and initial conditions
+# function spaces, solvers, and initial conditions
 # are redefined as in the previous demo. The only difference
 # is that now we are solving the adjoint problem, which
 # requires that the PDE solve is labelled with an
@@ -29,8 +29,8 @@ def get_function_spaces(mesh):
     return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         u, u_ = mesh_seq.fields["u"]
 
         # Define constants
@@ -45,17 +45,6 @@ def form(index):
             + inner(dot(u, nabla_grad(u)), v) * dx
             + nu * inner(grad(u), grad(v)) * dx
         )
-        return {"u": F}
-
-    return form
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        u, u_ = mesh_seq.fields["u"]
-
-        # Define form
-        F = mesh_seq.form(index)["u"]
 
         # Time integrate from t_start to t_end
         tp = mesh_seq.time_partition
@@ -125,7 +114,6 @@ def end_time_qoi():
     mesh,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
diff --git a/demos/burgers2.py b/demos/burgers2.py
index 14a99826..ebb94284 100644
--- a/demos/burgers2.py
+++ b/demos/burgers2.py
@@ -24,8 +24,8 @@ def get_function_spaces(mesh):
     return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         u, u_ = mesh_seq.fields["u"]
 
         # Define constants
@@ -40,17 +40,6 @@ def form(index):
             + inner(dot(u, nabla_grad(u)), v) * dx
             + nu * inner(grad(u), grad(v)) * dx
         )
-        return {"u": F}
-
-    return form
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        u, u_ = mesh_seq.fields["u"]
-
-        # Define form
-        F = mesh_seq.form(index)["u"]
 
         # Time integrate from t_start to t_end
         tp = mesh_seq.time_partition
@@ -105,7 +94,6 @@ def end_time_qoi():
     meshes,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
diff --git a/demos/gray_scott.py b/demos/gray_scott.py
index 9d1a8c9d..8ef350d0 100644
--- a/demos/gray_scott.py
+++ b/demos/gray_scott.py
@@ -46,16 +46,14 @@ def get_initial_condition(mesh_seq):
     return {"ab": ab_init}
 
 
-# Since we are using a mixed formulation, the forms for each component equation are
-# summed together. ::
+# For the solver, we just use the default configuration of Firedrake's
+# :class:`NonlinearVariationalSolver`. Since we are using a mixed formulation, the forms
+# for each component equation are summed together. ::
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         ab, ab_ = mesh_seq.fields["ab"]
-        a, b = split(ab)
-        a_, b_ = split(ab_)
-        psi_a, psi_b = TestFunctions(mesh_seq.function_spaces["ab"][index])
 
         # Define constants
         R = FunctionSpace(mesh_seq[index], "R", 0)
@@ -66,6 +64,9 @@ def form(index):
         kappa = Function(R).assign(0.06)
 
         # Write the two equations in variational form
+        a, b = split(ab)
+        a_, b_ = split(ab_)
+        psi_a, psi_b = TestFunctions(mesh_seq.function_spaces["ab"][index])
         F = (
             psi_a * (a - a_) * dx
             + dt * D_a * inner(grad(psi_a), grad(a)) * dx
@@ -74,21 +75,8 @@ def form(index):
             + dt * D_b * inner(grad(psi_b), grad(b)) * dx
             - dt * psi_b * (a * b**2 - (gamma + kappa) * b) * dx
         )
-        return {"ab": F}
-
-    return form
-
-
-# For the solver, we just use the default configuration of Firedrake's
-# :class:`NonlinearVariationalSolver`. ::
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        ab, ab_ = mesh_seq.fields["ab"]
 
         # Setup solver objects
-        F = mesh_seq.form(index)["ab"]
         nlvp = NonlinearVariationalProblem(F, ab)
         nlvs = NonlinearVariationalSolver(nlvp, ad_block_tag="ab")
 
@@ -150,7 +138,6 @@ def qoi():
     mesh,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
diff --git a/demos/gray_scott_split.py b/demos/gray_scott_split.py
index 5adcf89b..97f2462a 100644
--- a/demos/gray_scott_split.py
+++ b/demos/gray_scott_split.py
@@ -47,16 +47,14 @@ def get_initial_condition(mesh_seq):
     return {"a": a_init, "b": b_init}
 
 
-# We now see why :func:`get_form` needs to provide a function whose return value is a
-# dictionary: its keys correspond to the different equations being solved. ::
+# Correspondingly the solver needs to be constructed from the two parts and must
+# include two nonlinear solves at each timestep. ::
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         a, a_ = mesh_seq.fields["a"]
         b, b_ = mesh_seq.fields["b"]
-        psi_a = TestFunction(mesh_seq.function_spaces["a"][index])
-        psi_b = TestFunction(mesh_seq.function_spaces["b"][index])
 
         # Define constants
         R = FunctionSpace(mesh_seq[index], "R", 0)
@@ -67,6 +65,8 @@ def form(index):
         kappa = Function(R).assign(0.06)
 
         # Write the two equations in variational form
+        psi_a = TestFunction(mesh_seq.function_spaces["a"][index])
+        psi_b = TestFunction(mesh_seq.function_spaces["b"][index])
         F_a = (
             psi_a * (a - a_) * dx
             + dt * D_a * inner(grad(psi_a), grad(a)) * dx
@@ -77,24 +77,8 @@ def form(index):
             + dt * D_b * inner(grad(psi_b), grad(b)) * dx
             - dt * psi_b * (a * b**2 - (gamma + kappa) * b) * dx
         )
-        return {"a": F_a, "b": F_b}
-
-    return form
-
-
-# Correspondingly the solver needs to be constructed from the two parts and must
-# include two nonlinear solves at each timestep. ::
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        a, a_ = mesh_seq.fields["a"]
-        b, b_ = mesh_seq.fields["b"]
 
         # Setup solver objects
-        forms = mesh_seq.form(index)
-        F_a = forms["a"]
-        F_b = forms["b"]
         nlvp_a = NonlinearVariationalProblem(F_a, a)
         nlvs_a = NonlinearVariationalSolver(nlvp_a, ad_block_tag="a")
         nlvp_b = NonlinearVariationalProblem(F_b, b)
@@ -154,7 +138,6 @@ def qoi():
     mesh,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
diff --git a/demos/ode.py b/demos/ode.py
index ceb427a2..dd456808 100644
--- a/demos/ode.py
+++ b/demos/ode.py
@@ -101,29 +101,10 @@ def get_initial_condition(point_seq):
 # because we have a scalar :math:`R`-space it is a summation of a single value. Again,
 # this machinery may seem excessive but it becomes necessary for PDE problems.
 #
-# The Forward Euler scheme may be implemented as follows. ::
+# The Forward Euler scheme may be implemented and solved as follows. ::
 
 
-def get_form_forward_euler(point_seq):
-    def form(index):
-        R = point_seq.function_spaces["u"][index]
-        v = TestFunction(R)
-        u, u_ = point_seq.fields["u"]
-        dt = Function(R).assign(point_seq.time_partition.timesteps[index])
-
-        # Setup variational problem
-        F = (u - u_ - dt * u_) * v * dx
-        return {"u": F}
-
-    return form
-
-
-# We have a method defining the Forward Euler scheme. To put it into practice, we need
-# to define the solver. This boils down to applying the update repeatedly in a time
-# loop. ::
-
-
-def get_solver(point_seq):
+def get_solver_forward_euler(point_seq):
     def solver(index):
         tp = point_seq.time_partition
 
@@ -131,7 +112,10 @@ def solver(index):
         u, u_ = point_seq.fields["u"]
 
         # Define the (trivial) form
-        F = point_seq.form(index)["u"]
+        R = point_seq.function_spaces["u"][index]
+        dt = Function(R).assign(tp.timesteps[index])
+        v = TestFunction(R)
+        F = (u - u_ - dt * u_) * v * dx
 
         # Since the form is trivial, we can solve with a single application of a Jacobi
         # preconditioner
@@ -159,8 +143,7 @@ def solver(index):
     time_partition,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form_forward_euler,
-    get_solver=get_solver,
+    get_solver=get_solver_forward_euler,
 )
 
 # We can solve the ODE using the :meth:`~.MeshSeq.solve_forward` method and extract the
@@ -210,32 +193,40 @@ def solver(index):
 # .. math::
 #    \int_0^1 (u_{i+1} - u_{i} - \Delta t u_{i+1}) v \mathrm{d}t, \forall v\in R.
 #
-# ::
+# To apply Backward Euler we create the :class:`~.PointSeq` in the same way, just with
+# `get_solver_forward_euler` substituted for `get_solver_backward_euler`. ::
 
 
-def get_form_backward_euler(point_seq):
-    def form(index):
+def get_solver_backward_euler(point_seq):
+    def solver(index):
+        tp = point_seq.time_partition
+        u, u_ = point_seq.fields["u"]
         R = point_seq.function_spaces["u"][index]
+        dt = Function(R).assign(tp.timesteps[index])
         v = TestFunction(R)
-        u, u_ = point_seq.fields["u"]
-        dt = Function(R).assign(point_seq.time_partition.timesteps[index])
 
-        # Setup variational problem
+        # This is the only change from the Forward Euler solver
         F = (u - u_ - u * dt) * v * dx
-        return {"u": F}
 
-    return form
+        sp = {"ksp_type": "preonly", "pc_type": "jacobi"}
+        dt = tp.timesteps[index]
+        t_start, t_end = tp.subintervals[index]
+        t = t_start
+        while t < t_end - 1.0e-05:
+            solve(F == 0, u, solver_parameters=sp)
+            yield
+
+            u_.assign(u)
+            t += dt
 
+    return solver
 
-# To apply Backward Euler we create the :class:`~.PointSeq` in the same way, just with
-# `get_form_forward_euler` substituted for `get_form_backward_euler`. ::
 
 point_seq = PointSeq(
     time_partition,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form_backward_euler,
-    get_solver=get_solver,
+    get_solver=get_solver_backward_euler,
 )
 solutions = point_seq.solve_forward()["u"]["forward"]
 
@@ -269,27 +260,37 @@ def form(index):
 # where :math:`\theta\in(0,1)`. The standard choice is to take :math:`\theta=\frac12`. ::
 
 
-def get_form_crank_nicolson(point_seq):
-    def form(index):
+def get_solver_crank_nicolson(point_seq):
+    def solver(index):
+        tp = point_seq.time_partition
+        u, u_ = point_seq.fields["u"]
         R = point_seq.function_spaces["u"][index]
+        dt = Function(R).assign(tp.timesteps[index])
         v = TestFunction(R)
-        u, u_ = point_seq.fields["u"]
-        dt = Function(R).assign(point_seq.time_partition.timesteps[index])
-        theta = Function(R).assign(0.5)
 
-        # Setup variational problem
+        # This is the only change from the Forward and Backward Euler solvers
+        theta = Function(R).assign(0.5)
         F = (u - u_ - dt * (theta * u + (1 - theta) * u_)) * v * dx
-        return {"u": F}
 
-    return form
+        sp = {"ksp_type": "preonly", "pc_type": "jacobi"}
+        dt = tp.timesteps[index]
+        t_start, t_end = tp.subintervals[index]
+        t = t_start
+        while t < t_end - 1.0e-05:
+            solve(F == 0, u, solver_parameters=sp)
+            yield
+
+            u_.assign(u)
+            t += dt
+
+    return solver
 
 
 point_seq = PointSeq(
     time_partition,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form_crank_nicolson,
-    get_solver=get_solver,
+    get_solver=get_solver_crank_nicolson,
 )
 
 solutions = point_seq.solve_forward()["u"]["forward"]
diff --git a/demos/point_discharge2d-hessian.py b/demos/point_discharge2d-hessian.py
index 7f9478fd..980ef420 100644
--- a/demos/point_discharge2d-hessian.py
+++ b/demos/point_discharge2d-hessian.py
@@ -36,10 +36,10 @@ def source(mesh):
     return 100.0 * exp(-((x - x0) ** 2 + (y - y0) ** 2) / r**2)
 
 
-def get_form(mesh_seq):
-    def form(index):
-        c = mesh_seq.fields["c"]
+def get_solver(mesh_seq):
+    def solver(index):
         function_space = mesh_seq.function_spaces["c"][index]
+        c = mesh_seq.fields["c"]
         h = CellSize(mesh_seq[index])
         S = source(mesh_seq[index])
 
@@ -63,18 +63,6 @@ def form(index):
             + inner(D * grad(c), grad(psi)) * dx
             - S * psi * dx
         )
-        return {"c": F}
-
-    return form
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        function_space = mesh_seq.function_spaces["c"][index]
-        c = mesh_seq.fields["c"]
-
-        # Setup variational problem
-        F = mesh_seq.form(index)["c"]
         bc = DirichletBC(function_space, 0, 1)
 
         solve(F == 0, c, bcs=bc, ad_block_tag="c")
@@ -92,7 +80,6 @@ def solver(index):
     time_partition,
     mesh,
     get_function_spaces=get_function_spaces,
-    get_form=get_form,
     get_solver=get_solver,
 )
 
diff --git a/demos/solid_body_rotation.py b/demos/solid_body_rotation.py
index b95c0cdb..af807b50 100644
--- a/demos/solid_body_rotation.py
+++ b/demos/solid_body_rotation.py
@@ -132,19 +132,18 @@ def get_initial_condition(mesh_seq):
 #    :figwidth: 90%
 #    :align: center
 
-# Now let's set up the solver. First, we need to write the
-# :meth:`get_form` method. There is no integration by parts
+# Now let's set up the solver. First, we need to define the
+# variational problem. There is no integration by parts
 # and we apply Crank-Nicolson timestepping with implicitness
 # one half. Since we have a linear PDE, we write the variational
 # problem in terms of a left-hand side and right-hand side and
 # output both of them. ::
 
 
-def get_form(mesh_seq):
-    def form(index):
-        c, c_ = mesh_seq.fields["c"]
+def get_solver(mesh_seq):
+    def solver(index):
         V = mesh_seq.function_spaces["c"][index]
-        mesh = mesh_seq[index]
+        c, c_ = mesh_seq.fields["c"]
 
         # Define velocity field
         x, y = SpatialCoordinate(mesh)
@@ -155,28 +154,14 @@ def form(index):
         dt = Function(R).assign(mesh_seq.time_partition.timesteps[index])
         theta = Function(R).assign(0.5)
 
+        # Setup variational problem
         psi = TrialFunction(V)
         phi = TestFunction(V)
         a = psi * phi * dx + dt * theta * dot(u, grad(psi)) * phi * dx
         L = c_ * phi * dx - dt * (1 - theta) * dot(u, grad(c_)) * phi * dx
-        return {"c": (a, L)}
-
-    return form
-
-
-# The :func:`get_form` function is then used by :func:`get_solver`. ::
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        function_space = mesh_seq.function_spaces["c"][index]
-        c, c_ = mesh_seq.fields["c"]
-
-        # Setup variational problem
-        a, L = mesh_seq.form(index)["c"]
 
         # Zero Dirichlet condition on the boundary
-        bcs = DirichletBC(function_space, 0, "on_boundary")
+        bcs = DirichletBC(V, 0, "on_boundary")
 
         # Setup the solver object
         lvp = LinearVariationalProblem(a, L, c, bcs=bcs)
@@ -226,7 +211,6 @@ def qoi():
     mesh,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",

From 414a112abe0718edbf58135c2c2ea21eb9a6878f Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 17:40:56 +0000
Subject: [PATCH 03/13] #223: Update demos without GoalOrientedMeshSeq 2

---
 demos/burgers_time_integrated.py | 40 +++++++++++---------------------
 goalie/adjoint.py                |  1 -
 2 files changed, 14 insertions(+), 27 deletions(-)

diff --git a/demos/burgers_time_integrated.py b/demos/burgers_time_integrated.py
index 51a775c9..cedd8c77 100644
--- a/demos/burgers_time_integrated.py
+++ b/demos/burgers_time_integrated.py
@@ -12,35 +12,14 @@
 
 from goalie_adjoint import *
 
-# Redefine the ``get_initial_condition``, ``get_function_spaces``,
-# and ``get_form`` functions as in the first Burgers demo. ::
+# Redefine the ``get_initial_condition`` and ``get_function_spaces``,
+# functions as in the first Burgers demo. ::
 
 
 def get_function_spaces(mesh):
     return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-def get_form(mesh_seq):
-    def form(index):
-        u, u_ = mesh_seq.fields["u"]
-
-        # Define constants
-        R = FunctionSpace(mesh_seq[index], "R", 0)
-        dt = Function(R).assign(mesh_seq.time_partition.timesteps[index])
-        nu = Function(R).assign(0.0001)
-
-        # Setup variational problem
-        v = TestFunction(u.function_space())
-        F = (
-            inner((u - u_) / dt, v) * dx
-            + inner(dot(u, nabla_grad(u)), v) * dx
-            + nu * inner(grad(u), grad(v)) * dx
-        )
-        return {"u": F}
-
-    return form
-
-
 def get_initial_condition(mesh_seq):
     fs = mesh_seq.function_spaces["u"][0]
     x, y = SpatialCoordinate(mesh_seq[0])
@@ -61,8 +40,18 @@ def get_solver(mesh_seq):
     def solver(index):
         u, u_ = mesh_seq.fields["u"]
 
-        # Define form
-        F = mesh_seq.form(index)["u"]
+        # Define constants
+        R = FunctionSpace(mesh_seq[index], "R", 0)
+        dt = Function(R).assign(mesh_seq.time_partition.timesteps[index])
+        nu = Function(R).assign(0.0001)
+
+        # Setup variational problem
+        v = TestFunction(u.function_space())
+        F = (
+            inner((u - u_) / dt, v) * dx
+            + inner(dot(u, nabla_grad(u)), v) * dx
+            + nu * inner(grad(u), grad(v)) * dx
+        )
 
         # Time integrate from t_start to t_end
         t_start, t_end = mesh_seq.subintervals[index]
@@ -126,7 +115,6 @@ def time_integrated_qoi(t):
     meshes,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="time_integrated",
diff --git a/goalie/adjoint.py b/goalie/adjoint.py
index 72ba6686..95e8a4a7 100644
--- a/goalie/adjoint.py
+++ b/goalie/adjoint.py
@@ -85,7 +85,6 @@ def __init__(self, time_partition, initial_meshes, **kwargs):
             :meth:`~.MeshSeq.get_function_spaces`
         :kwarg get_initial_condition: a function as described in
             :meth:`~.MeshSeq.get_initial_condition`
-        :kwarg get_form: a function as described in :meth:`~.MeshSeq.get_form`
         :kwarg get_solver: a function as described in :meth:`~.MeshSeq.get_solver`
         :kwarg get_qoi: a function as described in :meth:`~.AdjointMeshSeq.get_qoi`
         """

From a85ff3826d2dec1dd8b3d7d745b492504b0a514c Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 17:54:39 +0000
Subject: [PATCH 04/13] #223: Handle forms in GoalOrientedMeshSeq

---
 goalie/go_mesh_seq.py | 25 +++++++++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index c3c479ad..3e23f057 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -27,6 +27,28 @@ class GoalOrientedMeshSeq(AdjointMeshSeq):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
         self.estimator_values = []
+        self._forms = None
+
+    def read_forms(self, forms_dictionary):
+        """
+        Read in the variational form corresponding to each prognostic field.
+
+        :arg forms_dictionary: dictionary where the keys are the field names and the
+            values are the UFL forms
+        :type forms_dictionary: :class:`dict`
+        """
+        self._forms = forms_dictionary
+
+    @property
+    def forms(self):
+        """
+        Get the variational form associated with each prognostic field.
+
+        :returns: dictionary where the keys are the field names and the values are the
+            UFL forms
+        :rtype: :class:`dict`
+        """
+        return self._forms
 
     @PETSc.Log.EventDecorator()
     def get_enriched_mesh_seq(self, enrichment_method="p", num_enrichments=1):
@@ -67,7 +89,6 @@ def get_enriched_mesh_seq(self, enrichment_method="p", num_enrichments=1):
             meshes,
             get_function_spaces=self._get_function_spaces,
             get_initial_condition=self._get_initial_condition,
-            get_form=self._get_form,
             get_solver=self._get_solver,
             get_qoi=self._get_qoi,
             qoi_type=self.qoi_type,
@@ -189,7 +210,7 @@ def indicate_errors(
 
             # Get forms for each equation in enriched space
             enriched_mesh_seq.fields = mapping
-            forms = enriched_mesh_seq.form(i)
+            forms = enriched_mesh_seq.forms
 
             # Loop over each timestep
             for j in range(self.time_partition.num_exports_per_subinterval[i] - 1):

From fe7d21b2f9de86520ed1b6a11c257588d7458027 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 17:54:54 +0000
Subject: [PATCH 05/13] #223: Update demos with GoalOrientedMeshSeq

---
 demos/burgers_ee.py                      | 22 +++++--------
 demos/point_discharge2d-goal_oriented.py | 23 ++++---------
 demos/point_discharge2d.py               | 42 +++++++++---------------
 3 files changed, 29 insertions(+), 58 deletions(-)

diff --git a/demos/burgers_ee.py b/demos/burgers_ee.py
index 66f96b09..3fe848a0 100644
--- a/demos/burgers_ee.py
+++ b/demos/burgers_ee.py
@@ -31,7 +31,10 @@
 set_log_level(DEBUG)
 
 # Redefine the ``field_names`` variable and the getter functions as in the first
-# adjoint Burgers demo. ::
+# adjoint Burgers demo. The only difference is the inclusion of the
+# :meth:`GoalOrientedMeshSeq.read_forms()` method in the ``get_solver`` function. The
+# method is used to communicate the variational form to the mesh sequence object so that
+# Goalie can utilise it in the error estimation process described above. ::
 
 field_names = ["u"]
 
@@ -40,8 +43,8 @@ def get_function_spaces(mesh):
     return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
         u, u_ = mesh_seq.fields["u"]
 
         # Define constants
@@ -56,17 +59,9 @@ def form(index):
             + inner(dot(u, nabla_grad(u)), v) * dx
             + nu * inner(grad(u), grad(v)) * dx
         )
-        return {"u": F}
-
-    return form
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        u, u_ = mesh_seq.fields["u"]
 
-        # Define form
-        F = mesh_seq.form(index)["u"]
+        # Communicate variational form to mesh_seq
+        mesh_seq.read_forms({"u": F})
 
         # Time integrate from t_start to t_end
         P = mesh_seq.time_partition
@@ -122,7 +117,6 @@ def end_time_qoi():
     meshes,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
diff --git a/demos/point_discharge2d-goal_oriented.py b/demos/point_discharge2d-goal_oriented.py
index 8b420320..92cf5c41 100644
--- a/demos/point_discharge2d-goal_oriented.py
+++ b/demos/point_discharge2d-goal_oriented.py
@@ -32,10 +32,10 @@ def source(mesh):
     return 100.0 * exp(-((x - x0) ** 2 + (y - y0) ** 2) / r**2)
 
 
-def get_form(mesh_seq):
-    def form(index):
-        c = mesh_seq.fields["c"]
+def get_solver(mesh_seq):
+    def solver(index):
         function_space = mesh_seq.function_spaces["c"][index]
+        c = mesh_seq.fields["c"]
         h = CellSize(mesh_seq[index])
         S = source(mesh_seq[index])
 
@@ -59,20 +59,11 @@ def form(index):
             + inner(D * grad(c), grad(psi)) * dx
             - S * psi * dx
         )
-        return {"c": F}
-
-    return form
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        function_space = mesh_seq.function_spaces["c"][index]
-        c = mesh_seq.fields["c"]
-
-        # Setup variational problem
-        F = mesh_seq.form(index)["c"]
         bc = DirichletBC(function_space, 0, 1)
 
+        # Communicate variational form to mesh_seq
+        mesh_seq.read_forms({"c": F})
+
         solve(F == 0, c, bcs=bc, ad_block_tag="c")
         yield
 
@@ -99,7 +90,6 @@ def qoi():
     time_partition,
     mesh,
     get_function_spaces=get_function_spaces,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="steady",
@@ -328,7 +318,6 @@ def adaptor(mesh_seq, solutions, indicators):
     time_partition,
     mesh,
     get_function_spaces=get_function_spaces,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="steady",
diff --git a/demos/point_discharge2d.py b/demos/point_discharge2d.py
index 3c438966..b1eee69e 100644
--- a/demos/point_discharge2d.py
+++ b/demos/point_discharge2d.py
@@ -65,13 +65,21 @@ def source(mesh):
 # .. math::
 #    \tau = \min\left(\frac{h}{2\|\mathbf u\|},\frac{h\|\mathbf u\|}{6D}\right),
 #
-# where :math:`h` measures cell size. ::
+# where :math:`h` measures cell size.
+#
+# Note that :attr:`mesh_seq.fields` now returns a single
+# :class:`~firedrake.function.Function` object since the problem is steady, so there is
+# no notion of a lagged solution, unlike in previous (time-dependent) demos.
+# With these ingredients, we can now define the :meth:`get_solver` method. Don't forget
+# to apply the corresponding `ad_block_tag` to the solve call. Additionally, we must
+# communicate the defined variational form to ``mesh_seq`` using the
+# :meth:`mesh_seq.read_form()` method for Goalie to utilise it during error indication. ::
 
 
-def get_form(mesh_seq):
-    def form(index):
-        c = mesh_seq.fields["c"]
+def get_solver(mesh_seq):
+    def solver(index):
         function_space = mesh_seq.function_spaces["c"][index]
+        c = mesh_seq.fields["c"]
         h = CellSize(mesh_seq[index])
         S = source(mesh_seq[index])
 
@@ -95,30 +103,11 @@ def form(index):
             + inner(D * grad(c), grad(psi)) * dx
             - S * psi * dx
         )
-        return {"c": F}
-
-    return form
-
-
-# Note that :attr:`mesh_seq.fields` now returns a single
-# :class:`~firedrake.function.Function` object since the problem is steady, so there is
-# no notion of a lagged solution, unlike in previous (time-dependent) demos.
-# With these ingredients, we can now define the :meth:`get_solver` method. Don't forget
-# to apply the corresponding `ad_block_tag` to the solve call. ::
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        function_space = mesh_seq.function_spaces["c"][index]
-
-        c = mesh_seq.fields["c"]
-
-        # Setup variational problem
-        F = mesh_seq.form(index)["c"]
-
-        # Strongly enforce boundary conditions on the inflow, which is indexed by 1
         bc = DirichletBC(function_space, 0, 1)
 
+        # Communicate variational form to mesh_seq
+        mesh_seq.read_forms({"c": F})
+
         solve(F == 0, c, bcs=bc, ad_block_tag="c")
         yield
 
@@ -158,7 +147,6 @@ def qoi():
     time_partition,
     mesh,
     get_function_spaces=get_function_spaces,
-    get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="steady",

From 10692188e4f5a8d2028fd308e5be103ee0f81a29 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 17:58:32 +0000
Subject: [PATCH 06/13] #223: Update demos with GoalOrientedMeshSeq 2

---
 demos/burgers_oo.py | 15 ++++-----------
 1 file changed, 4 insertions(+), 11 deletions(-)

diff --git a/demos/burgers_oo.py b/demos/burgers_oo.py
index 145a1a14..a38ec020 100644
--- a/demos/burgers_oo.py
+++ b/demos/burgers_oo.py
@@ -30,8 +30,8 @@ class BurgersMeshSeq(GoalOrientedMeshSeq):
     def get_function_spaces(mesh):
         return {"u": VectorFunctionSpace(mesh, "CG", 2)}
 
-    def get_form(self):
-        def form(index):
+    def get_solver(self):
+        def solver(index):
             u, u_ = self.fields["u"]
 
             # Define constants
@@ -46,16 +46,9 @@ def form(index):
                 + inner(dot(u, nabla_grad(u)), v) * dx
                 + nu * inner(grad(u), grad(v)) * dx
             )
-            return {"u": F}
-
-        return form
-
-    def get_solver(self):
-        def solver(index):
-            u, u_ = self.fields["u"]
 
-            # Define form
-            F = self.form(index)["u"]
+            # Communicate variational form to mesh_seq
+            self.read_forms({"u": F})
 
             # Time integrate from t_start to t_end
             t_start, t_end = self.subintervals[index]

From 3305d3d34e2c4543fe1d509624c25b37b1f83a5c Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 18:15:45 +0000
Subject: [PATCH 07/13] #223: Update tests

---
 test/test_mesh_seq.py                         | 40 +++++--------------
 test_adjoint/examples/burgers.py              | 29 ++++----------
 test_adjoint/examples/point_discharge2d.py    | 27 ++++---------
 test_adjoint/examples/point_discharge3d.py    | 28 ++++---------
 test_adjoint/examples/steady_flow_past_cyl.py | 28 ++++---------
 test_adjoint/test_adjoint.py                  |  3 --
 test_adjoint/test_fp_iteration.py             |  1 -
 7 files changed, 39 insertions(+), 117 deletions(-)

diff --git a/test/test_mesh_seq.py b/test/test_mesh_seq.py
index 0b19e953..924014fb 100644
--- a/test/test_mesh_seq.py
+++ b/test/test_mesh_seq.py
@@ -36,7 +36,7 @@ def test_inconsistent_dim(self):
         msg = "Meshes must all have the same topological dimension."
         self.assertEqual(str(cm.exception), msg)
 
-    @parameterized.expand(["get_function_spaces", "get_form", "get_solver"])
+    @parameterized.expand(["get_function_spaces", "get_solver"])
     def test_notimplemented_error(self, function_name):
         mesh_seq = MeshSeq(self.time_interval, UnitSquareMesh(1, 1))
         with self.assertRaises(NotImplementedError) as cm:
@@ -47,54 +47,35 @@ def test_notimplemented_error(self, function_name):
         msg = f"'{function_name}' needs implementing."
         self.assertEqual(str(cm.exception), msg)
 
-    @parameterized.expand(["get_function_spaces", "get_initial_condition", "get_form"])
+    @parameterized.expand(["get_function_spaces", "get_initial_condition"])
     def test_return_dict_error(self, method):
         mesh = UnitSquareMesh(1, 1)
-        methods = ["get_function_spaces", "get_initial_condition", "get_form"]
+        methods = ["get_function_spaces", "get_initial_condition"]
         funcs = [lambda _: 0, lambda _: 0, lambda _: lambda *_: 0]
-        methods_map = dict(zip(methods, funcs))
-        if method == "get_form":
-            kwargs = methods_map
-            f_space = FunctionSpace(mesh, "CG", 1)
-            kwargs["get_function_spaces"] = lambda _: {"field": f_space}
-            kwargs["get_initial_condition"] = lambda _: {"field": Function(f_space)}
-        else:
-            kwargs = {method: methods_map[method]}
+        kwargs = {method: funcs[methods.index(method)]}
         with self.assertRaises(AssertionError) as cm:
             MeshSeq(self.time_interval, mesh, **kwargs)
         msg = f"{method} should return a dict"
         self.assertEqual(str(cm.exception), msg)
 
-    @parameterized.expand(["get_function_spaces", "get_initial_condition", "get_form"])
+    @parameterized.expand(["get_function_spaces", "get_initial_condition"])
     def test_missing_field_error(self, method):
         mesh = UnitSquareMesh(1, 1)
-        methods = ["get_function_spaces", "get_initial_condition", "get_form"]
+        methods = ["get_function_spaces", "get_initial_condition"]
         funcs = [lambda _: {}, lambda _: {}, lambda _: lambda *_: {}]
-        kwargs = dict(zip(methods, funcs))
-        if method == "get_form":
-            f_space = FunctionSpace(mesh, "CG", 1)
-            kwargs["get_function_spaces"] = lambda _: {"field": f_space}
-            kwargs["get_initial_condition"] = lambda _: {"field": Function(f_space)}
-        else:
-            kwargs = {method: kwargs[method]}
+        kwargs = {method: funcs[methods.index(method)]}
         with self.assertRaises(AssertionError) as cm:
             MeshSeq(self.time_interval, mesh, **kwargs)
         msg = "missing fields {'field'} in " + f"{method}"
         self.assertEqual(str(cm.exception), msg)
 
-    @parameterized.expand(["get_function_spaces", "get_initial_condition", "get_form"])
+    @parameterized.expand(["get_function_spaces", "get_initial_condition"])
     def test_unexpected_field_error(self, method):
         mesh = UnitSquareMesh(1, 1)
-        methods = ["get_function_spaces", "get_initial_condition", "get_form"]
+        methods = ["get_function_spaces", "get_initial_condition"]
         out_dict = {"field": None, "extra_field": None}
         funcs = [lambda _: out_dict, lambda _: out_dict, lambda _: lambda *_: out_dict]
-        kwargs = dict(zip(methods, funcs))
-        if method == "get_form":
-            f_space = FunctionSpace(mesh, "CG", 1)
-            kwargs["get_function_spaces"] = lambda _: {"field": f_space}
-            kwargs["get_initial_condition"] = lambda _: {"field": Function(f_space)}
-        else:
-            kwargs = {method: kwargs[method]}
+        kwargs = {method: funcs[methods.index(method)]}
         with self.assertRaises(AssertionError) as cm:
             MeshSeq(self.time_interval, mesh, **kwargs)
         msg = "unexpected fields {'extra_field'} in " + f"{method}"
@@ -106,7 +87,6 @@ def test_solver_generator_error(self):
         kwargs = {
             "get_function_spaces": lambda _: {"field": f_space},
             "get_initial_condition": lambda _: {"field": Function(f_space)},
-            "get_form": lambda msq: lambda *_: {"field": msq.fields["field"][0]},
             "get_solver": lambda _: lambda *_: {},
         }
         with self.assertRaises(AssertionError) as cm:
diff --git a/test_adjoint/examples/burgers.py b/test_adjoint/examples/burgers.py
index 1f0aa934..6e1042ea 100644
--- a/test_adjoint/examples/burgers.py
+++ b/test_adjoint/examples/burgers.py
@@ -28,13 +28,18 @@ def get_function_spaces(mesh):
     return {"uv_2d": VectorFunctionSpace(mesh, "CG", 2)}
 
 
-def get_form(self):
+def get_solver(self):
     """
-    Burgers equation weak form.
+    Burgers equation solved using a direct method and backward Euler timestepping.
     """
 
-    def form(i):
+    def solver(i):
+        t_start, t_end = self.time_partition.subintervals[i]
+        dt = self.time_partition.timesteps[i]
+
         u, u_ = self.fields["uv_2d"]
+
+        # Setup variational problem
         dt = self.time_partition.timesteps[i]
         fs = self.function_spaces["uv_2d"][i]
         R = FunctionSpace(self[i], "R", 0)
@@ -46,24 +51,6 @@ def form(i):
             + inner(dot(u, nabla_grad(u)), v) * dx
             + nu * inner(grad(u), grad(v)) * dx
         )
-        return {"uv_2d": F}
-
-    return form
-
-
-def get_solver(self):
-    """
-    Burgers equation solved using a direct method and backward Euler timestepping.
-    """
-
-    def solver(i):
-        t_start, t_end = self.time_partition.subintervals[i]
-        dt = self.time_partition.timesteps[i]
-
-        u, u_ = self.fields["uv_2d"]
-
-        # Setup variational problem
-        F = self.form(i)["uv_2d"]
 
         # Time integrate from t_start to t_end
         t = t_start
diff --git a/test_adjoint/examples/point_discharge2d.py b/test_adjoint/examples/point_discharge2d.py
index 76e02658..e0f18ceb 100644
--- a/test_adjoint/examples/point_discharge2d.py
+++ b/test_adjoint/examples/point_discharge2d.py
@@ -51,13 +51,17 @@ def source(mesh):
     return 100.0 * exp(-((x - src_x) ** 2 + (y - src_y) ** 2) / src_r**2)
 
 
-def get_form(self):
+def get_solver(self):
     """
-    Advection-diffusion with SUPG stabilisation.
+    Advection-diffusion equation
+    solved using a direct method.
     """
 
-    def form(i):
+    def solver(i):
+        fs = self.function_spaces["tracer_2d"][i]
         c = self.fields["tracer_2d"]
+
+        # Define constants
         fs = self.function_spaces["tracer_2d"][i]
         R = FunctionSpace(self[i], "R", 0)
         D = Function(R).assign(0.1)
@@ -80,23 +84,6 @@ def form(i):
             - dot(u, grad(c)) * psi * dx
             - inner(D * grad(c), grad(psi)) * dx
         )
-        return {"tracer_2d": F}
-
-    return form
-
-
-def get_solver(self):
-    """
-    Advection-diffusion equation
-    solved using a direct method.
-    """
-
-    def solver(i):
-        fs = self.function_spaces["tracer_2d"][i]
-        c = self.fields["tracer_2d"]
-
-        # Setup variational problem
-        F = self.form(i)["tracer_2d"]
 
         # Zero Dirichlet condition on the left-hand (inlet) boundary
         bc = DirichletBC(fs, 0, 1)
diff --git a/test_adjoint/examples/point_discharge3d.py b/test_adjoint/examples/point_discharge3d.py
index 87764eed..29a4d36b 100644
--- a/test_adjoint/examples/point_discharge3d.py
+++ b/test_adjoint/examples/point_discharge3d.py
@@ -62,14 +62,17 @@ def source(mesh):
     )
 
 
-def get_form(self):
+def get_solver(self):
     """
-    Advection-diffusion with SUPG
-    stabilisation.
+    Advection-diffusion equation
+    solved using a direct method.
     """
 
-    def form(i):
+    def solver(i):
+        fs = self.function_spaces["tracer_3d"][i]
         c = self.fields["tracer_3d"]
+
+        # Define constants
         fs = self.function_spaces["tracer_3d"][i]
         R = FunctionSpace(self[i], "R", 0)
         D = Function(R).assign(0.1)
@@ -93,23 +96,6 @@ def form(i):
             - dot(u, grad(c)) * psi * dx
             - inner(D * grad(c), grad(psi)) * dx
         )
-        return {"tracer_3d": F}
-
-    return form
-
-
-def get_solver(self):
-    """
-    Advection-diffusion equation
-    solved using a direct method.
-    """
-
-    def solver(i):
-        fs = self.function_spaces["tracer_3d"][i]
-        c = self.fields["tracer_3d"]
-
-        # Setup variational problem
-        F = self.form(i)["tracer_3d"]
 
         # Zero Dirichlet condition on the left-hand (inlet) boundary
         bc = DirichletBC(fs, 0, 1)
diff --git a/test_adjoint/examples/steady_flow_past_cyl.py b/test_adjoint/examples/steady_flow_past_cyl.py
index b1e30cc4..811456e9 100644
--- a/test_adjoint/examples/steady_flow_past_cyl.py
+++ b/test_adjoint/examples/steady_flow_past_cyl.py
@@ -32,14 +32,17 @@ def get_function_spaces(mesh):
     return {"up": VectorFunctionSpace(mesh, "CG", 2) * FunctionSpace(mesh, "CG", 1)}
 
 
-def get_form(self):
+def get_solver(self):
     """
-    Weak form for Stokes equation.
+    Stokes problem solved using a
+    direct method.
     """
 
-    def form(i):
-        up = self.fields["up"]
+    def solver(i):
         W = self.function_spaces["up"][i]
+        up = self.fields["up"]
+
+        # Define variational problem
         R = FunctionSpace(self[i], "R", 0)
         nu = Function(R).assign(1.0)
         u, p = split(up)
@@ -50,23 +53,6 @@ def form(i):
             - inner(p, div(v)) * dx
             - inner(q, div(u)) * dx
         )
-        return {"up": F}
-
-    return form
-
-
-def get_solver(self):
-    """
-    Stokes problem solved using a
-    direct method.
-    """
-
-    def solver(i):
-        W = self.function_spaces["up"][i]
-        up = self.fields["up"]
-
-        # Define variational problem
-        F = self.form(i)["up"]
 
         # Define inflow and no-slip boundary conditions
         y = SpatialCoordinate(self[i])[1]
diff --git a/test_adjoint/test_adjoint.py b/test_adjoint/test_adjoint.py
index 43750cb3..b917260e 100644
--- a/test_adjoint/test_adjoint.py
+++ b/test_adjoint/test_adjoint.py
@@ -114,7 +114,6 @@ def test_adjoint_same_mesh(problem, qoi_type, debug=False):
         test_case.mesh,
         get_function_spaces=test_case.get_function_spaces,
         get_initial_condition=test_case.get_initial_condition,
-        get_form=test_case.get_form,
         get_solver=test_case.get_solver,
         get_qoi=test_case.get_qoi,
         qoi_type=qoi_type,
@@ -176,7 +175,6 @@ def test_adjoint_same_mesh(problem, qoi_type, debug=False):
             test_case.mesh,
             get_function_spaces=test_case.get_function_spaces,
             get_initial_condition=test_case.get_initial_condition,
-            get_form=test_case.get_form,
             get_solver=test_case.get_solver,
             get_qoi=test_case.get_qoi,
             qoi_type=qoi_type,
@@ -254,7 +252,6 @@ def plot_solutions(problem, qoi_type, debug=True):
         test_case.mesh,
         get_function_spaces=test_case.get_function_spaces,
         get_initial_condition=test_case.get_initial_condition,
-        get_form=test_case.get_form,
         get_solver=test_case.get_solver,
         get_qoi=test_case.get_qoi,
         qoi_type=qoi_type,
diff --git a/test_adjoint/test_fp_iteration.py b/test_adjoint/test_fp_iteration.py
index c962b7d2..08d87637 100644
--- a/test_adjoint/test_fp_iteration.py
+++ b/test_adjoint/test_fp_iteration.py
@@ -73,7 +73,6 @@ def mesh_seq(self, **kwargs):
         kw.update(kwargs)
         mesh_seq = self.seq(kw.pop("time_partition"), kw.pop("mesh"), **kw)
         mesh_seq._get_function_spaces = lambda _: {}
-        mesh_seq._get_form = lambda _: lambda *_: {}
         mesh_seq._get_solver = lambda _: lambda *_: (
             yield from (None for _ in iter(int, 1))
         )

From 848fef29db2b3c4b444dc438d153541603f5582f Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 23:16:36 +0000
Subject: [PATCH 08/13] #223: Fix form coefficient mapping

---
 goalie/go_mesh_seq.py | 13 ++++---------
 1 file changed, 4 insertions(+), 9 deletions(-)

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index 3e23f057..80f285e6 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -195,21 +195,16 @@ def indicate_errors(
             enriched_spaces = {
                 f: enriched_mesh_seq.function_spaces[f][i] for f in self.fields
             }
-            mapping = {}
             for f, fs_e in enriched_spaces.items():
-                u[f] = Function(fs_e)
-                u_[f] = Function(fs_e)
-                mapping[f] = (
-                    (u[f], u_[f])
-                    if enriched_mesh_seq.field_types[f] == "unsteady"
-                    else u[f]
-                )
+                if self.field_types[f] == "steady":
+                    u[f] = u_[f] = enriched_mesh_seq.fields[f]
+                else:
+                    u[f], u_[f] = enriched_mesh_seq.fields[f]
                 u_star[f] = Function(fs_e)
                 u_star_next[f] = Function(fs_e)
                 u_star_e[f] = Function(fs_e)
 
             # Get forms for each equation in enriched space
-            enriched_mesh_seq.fields = mapping
             forms = enriched_mesh_seq.forms
 
             # Loop over each timestep

From 75bd853fb51bd4e366958153fc0f6092463609fe Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 23:17:44 +0000
Subject: [PATCH 09/13] #223: Raise an error when read_forms wasn't used

---
 goalie/go_mesh_seq.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index 80f285e6..62cccc8f 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -48,6 +48,11 @@ def forms(self):
             UFL forms
         :rtype: :class:`dict`
         """
+        if self._forms is None:
+            raise AttributeError(
+                "Forms have not been read in."
+                " Use read_forms({'field_name': F}) in get_solver to read in the forms."
+            )
         return self._forms
 
     @PETSc.Log.EventDecorator()

From f80edbf0c75532d8646b0f76b1595eb2d47620bd Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 12 Nov 2024 23:41:41 +0000
Subject: [PATCH 10/13] #223: Minor fix

---
 goalie/go_mesh_seq.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index 62cccc8f..48657c98 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -209,9 +209,6 @@ def indicate_errors(
                 u_star_next[f] = Function(fs_e)
                 u_star_e[f] = Function(fs_e)
 
-            # Get forms for each equation in enriched space
-            forms = enriched_mesh_seq.forms
-
             # Loop over each timestep
             for j in range(self.time_partition.num_exports_per_subinterval[i] - 1):
                 # In case of having multiple solution fields that are solved for one
@@ -242,7 +239,7 @@ def indicate_errors(
                     u_star_e[f] -= u_star[f]
 
                     # Evaluate error indicator
-                    indi_e = indicator_fn(forms[f], u_star_e[f])
+                    indi_e = indicator_fn(enriched_mesh_seq.forms[f], u_star_e[f])
 
                     # Transfer back to the base space
                     indi = self._transfer(indi_e, P0_spaces[i])

From 912f54018f9ce3396884c598264e31b8e77c5e2e Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Wed, 13 Nov 2024 11:49:29 +0000
Subject: [PATCH 11/13] #223: Add field name and form type check

---
 goalie/go_mesh_seq.py         | 10 ++++++++++
 test_adjoint/test_mesh_seq.py | 26 ++++++++++++++++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index 48657c98..47a1066e 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -37,6 +37,16 @@ def read_forms(self, forms_dictionary):
             values are the UFL forms
         :type forms_dictionary: :class:`dict`
         """
+        for field, form in forms_dictionary.items():
+            if field not in self.fields:
+                raise ValueError(
+                    f"Unexpected field '{field}' in forms dictionary."
+                    f" Expected one of {self.time_partition.field_names}."
+                )
+            if not isinstance(form, ufl.Form):
+                raise TypeError(
+                    f"Expected a UFL form for field '{field}', not '{type(form)}'."
+                )
         self._forms = forms_dictionary
 
     @property
diff --git a/test_adjoint/test_mesh_seq.py b/test_adjoint/test_mesh_seq.py
index 8757253b..89ac6ce3 100644
--- a/test_adjoint/test_mesh_seq.py
+++ b/test_adjoint/test_mesh_seq.py
@@ -313,6 +313,32 @@ def go_mesh_seq(self, get_function_spaces, parameters=None):
         )
 
 
+class TestGoalOrientedMeshSeq(TrivialGoalOrientedBaseClass):
+    """
+    Unit tests for a :class:`GoalOrientedMeshSeq`.
+    """
+
+    def get_function_spaces(self, mesh):
+        return {self.field: FunctionSpace(mesh, "R", 0)}
+
+    def test_read_forms_error_field(self):
+        mesh_seq = self.go_mesh_seq(self.get_function_spaces)
+        with self.assertRaises(ValueError) as cm:
+            mesh_seq.read_forms({"field2": None})
+        msg = (
+            "Unexpected field 'field2' in forms dictionary."
+            f" Expected one of ['{self.field}']."
+        )
+        self.assertEqual(str(cm.exception), msg)
+
+    def test_read_forms_error_form(self):
+        mesh_seq = self.go_mesh_seq(self.get_function_spaces)
+        with self.assertRaises(TypeError) as cm:
+            mesh_seq.read_forms({self.field: None})
+        msg = f"Expected a UFL form for field '{self.field}', not '<class 'NoneType'>'."
+        self.assertEqual(str(cm.exception), msg)
+
+
 class TestGlobalEnrichment(TrivialGoalOrientedBaseClass):
     """
     Unit tests for global enrichment of a :class:`GoalOrientedMeshSeq`.

From 1903f96fbf64a7b3ad56dd99266de7b7fd767023 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Wed, 13 Nov 2024 13:06:04 +0000
Subject: [PATCH 12/13] #223: Tidy up tests (minor)

---
 test/test_mesh_seq.py | 13 +++----------
 1 file changed, 3 insertions(+), 10 deletions(-)

diff --git a/test/test_mesh_seq.py b/test/test_mesh_seq.py
index 924014fb..7c2287e6 100644
--- a/test/test_mesh_seq.py
+++ b/test/test_mesh_seq.py
@@ -50,9 +50,7 @@ def test_notimplemented_error(self, function_name):
     @parameterized.expand(["get_function_spaces", "get_initial_condition"])
     def test_return_dict_error(self, method):
         mesh = UnitSquareMesh(1, 1)
-        methods = ["get_function_spaces", "get_initial_condition"]
-        funcs = [lambda _: 0, lambda _: 0, lambda _: lambda *_: 0]
-        kwargs = {method: funcs[methods.index(method)]}
+        kwargs = {method: lambda _: 0}
         with self.assertRaises(AssertionError) as cm:
             MeshSeq(self.time_interval, mesh, **kwargs)
         msg = f"{method} should return a dict"
@@ -61,9 +59,7 @@ def test_return_dict_error(self, method):
     @parameterized.expand(["get_function_spaces", "get_initial_condition"])
     def test_missing_field_error(self, method):
         mesh = UnitSquareMesh(1, 1)
-        methods = ["get_function_spaces", "get_initial_condition"]
-        funcs = [lambda _: {}, lambda _: {}, lambda _: lambda *_: {}]
-        kwargs = {method: funcs[methods.index(method)]}
+        kwargs = {method: lambda _: {}}
         with self.assertRaises(AssertionError) as cm:
             MeshSeq(self.time_interval, mesh, **kwargs)
         msg = "missing fields {'field'} in " + f"{method}"
@@ -72,10 +68,7 @@ def test_missing_field_error(self, method):
     @parameterized.expand(["get_function_spaces", "get_initial_condition"])
     def test_unexpected_field_error(self, method):
         mesh = UnitSquareMesh(1, 1)
-        methods = ["get_function_spaces", "get_initial_condition"]
-        out_dict = {"field": None, "extra_field": None}
-        funcs = [lambda _: out_dict, lambda _: out_dict, lambda _: lambda *_: out_dict]
-        kwargs = {method: funcs[methods.index(method)]}
+        kwargs = {method: lambda _: {"field": None, "extra_field": None}}
         with self.assertRaises(AssertionError) as cm:
             MeshSeq(self.time_interval, mesh, **kwargs)
         msg = "unexpected fields {'extra_field'} in " + f"{method}"

From 210f2eaa60dd4070acea9b16630988ac3d58c199 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Thu, 14 Nov 2024 11:00:04 +0000
Subject: [PATCH 13/13] #223: Don't update forward_old soln if steady

---
 goalie/go_mesh_seq.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index 47a1066e..d79706c7 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -212,7 +212,7 @@ def indicate_errors(
             }
             for f, fs_e in enriched_spaces.items():
                 if self.field_types[f] == "steady":
-                    u[f] = u_[f] = enriched_mesh_seq.fields[f]
+                    u[f] = enriched_mesh_seq.fields[f]
                 else:
                     u[f], u_[f] = enriched_mesh_seq.fields[f]
                 u_star[f] = Function(fs_e)
@@ -233,7 +233,8 @@ def indicate_errors(
                 for f in self.fields:
                     # Transfer solutions associated with the current field f
                     transfer(self.solutions[f][FWD][i][j], u[f])
-                    transfer(self.solutions[f][FWD_OLD][i][j], u_[f])
+                    if self.field_types[f] == "unsteady":
+                        transfer(self.solutions[f][FWD_OLD][i][j], u_[f])
                     transfer(self.solutions[f][ADJ][i][j], u_star[f])
                     transfer(self.solutions[f][ADJ_NEXT][i][j], u_star_next[f])