Do not require get_form in user code (#228)

Closes #223.

We only need the form when we do `GoalOrientedMeshSeq.indicate_errors`,
so everything form-related was removed from `MeshSeq` and
`AdjointMeshSeq`. Communicating the form from user code to
`GoalOrientedMeshSeq` is now done with the `read_forms` method.

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,
 )
 

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()

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",

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",

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",

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]

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",