Merge develop into main

.github/CODEOWNERS

@@ -0,0 +1,15 @@
+# Each line is a file pattern followed by one or more owners.
+
+# These owners will be the default owners for everything in
+# the repo. Unless a later match takes precedence,
+# will be requested for review when someone opens a pull request.
+*      @jwallwork23 @stephankramer
+
+# Order is important; the last matching pattern takes the most
+# precedence. When someone opens a pull request that only
+# modifies specified file types, only that owner will be
+# requested for a review, not the global owner(s).
+
+# You can also use email addresses if you prefer. They'll be
+# used to look up users just like we do for commit author
+# emails.

demos/burgers-goal_oriented.py

@@ -8,21 +8,18 @@
 # multiple meshes to adapt. We also chose to apply a QoI which integrates in time as
 # well as space.
 #
-# We copy over the setup as before. The only difference is that we import from
-# ``goalie_adjoint`` rather than ``goalie``. ::
+# We copy over the setup as before. ::
 
 import matplotlib.pyplot as plt
 from animate.adapt import adapt
 from animate.metric import RiemannianMetric
 from firedrake import *
 
-from goalie_adjoint import *
+from goalie import *
 
-fields = ["u"]
-
-
-def get_function_spaces(mesh):
-    return {"u": VectorFunctionSpace(mesh, "CG", 2)}
+n = 32
+meshes = [UnitSquareMesh(n, n), UnitSquareMesh(n, n)]
+fields = [Field("u", family="Lagrange", degree=2, vector=True)]
 
 
 def get_initial_condition(mesh_seq):
@@ -37,7 +34,7 @@ def get_initial_condition(mesh_seq):
 
 def get_solver(mesh_seq):
     def solver(index):
-        u, u_ = mesh_seq.fields["u"]
+        u, u_ = mesh_seq.field_functions["u"]
 
         # Define constants
         R = FunctionSpace(mesh_seq[index], "R", 0)
@@ -77,7 +74,7 @@ def get_qoi(mesh_seq, i):
     dt = Function(R).assign(mesh_seq.time_partition.timesteps[i])
 
     def time_integrated_qoi(t):
-        u = mesh_seq.fields["u"][0]
+        u = mesh_seq.field_functions["u"][0]
         return dt * inner(u, u) * ds(2)
 
     return time_integrated_qoi
@@ -87,11 +84,8 @@ def time_integrated_qoi(t):
 # as in a `previous demo <./burgers2.py.html>`__, except that we export
 # every timestep rather than every other timestep. ::
 
-n = 32
-meshes = [UnitSquareMesh(n, n, diagonal="left"), UnitSquareMesh(n, n, diagonal="left")]
 end_time = 0.5
 dt = 1 / n
-
 num_subintervals = len(meshes)
 time_partition = TimePartition(
     end_time,
@@ -107,7 +101,6 @@ def time_integrated_qoi(t):
 mesh_seq = GoalOrientedMeshSeq(
     time_partition,
     meshes,
-    get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
     get_solver=get_solver,
     get_qoi=get_qoi,
@@ -178,7 +171,7 @@ def adaptor(mesh_seq, solutions=None, indicators=None):
     num_elem = mesh_seq.count_elements()
     pyrint(f"fixed point iteration {iteration + 1}:")
     for i, (complexity, ndofs, nelem) in enumerate(
-        zip(complexities, num_dofs, num_elem)
+        zip(complexities, num_dofs, num_elem, strict=True)
     ):
         pyrint(
             f"  subinterval {i}, complexity: {complexity:4.0f}"
@@ -339,7 +332,7 @@ def adaptor(mesh_seq, solutions=None, indicators=None):
     num_elem = mesh_seq.count_elements()
     pyrint(f"fixed point iteration {iteration + 1}:")
     for i, (complexity, ndofs, nelem) in enumerate(
-        zip(complexities, num_dofs, num_elem)
+        zip(complexities, num_dofs, num_elem, strict=False)
     ):
         pyrint(
             f"  subinterval {i}, complexity: {complexity:4.0f}"
@@ -365,7 +358,6 @@ def adaptor(mesh_seq, solutions=None, indicators=None):
 mesh_seq = GoalOrientedMeshSeq(
     time_partition,
     meshes,
-    get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
     get_solver=get_solver,
     get_qoi=get_qoi,

demos/burgers-hessian.py

@@ -7,7 +7,12 @@
 # we consider the time-dependent case. Moreover, we consider a :class:`MeshSeq` with
 # multiple subintervals and hence multiple meshes to adapt.
 #
-# As before, we copy over what is now effectively boiler plate to set up our problem. ::
+# As before, we copy over what is now effectively boiler plate to set up our problem.
+#
+# The only difference is that we need to specifically define the initial mesh for each
+# subinterval and pass them as a list. When a single mesh is passed to the
+# :class:`~.MeshSeq` constructor, it is shallowed copied, which is insufficient for mesh
+# adaptation. ::
 
 import matplotlib.pyplot as plt
 from animate.adapt import adapt
@@ -16,16 +21,14 @@
 
 from goalie import *
 
-field_names = ["u"]
-
-
-def get_function_spaces(mesh):
-    return {"u": VectorFunctionSpace(mesh, "CG", 2)}
+n = 32
+meshes = [UnitSquareMesh(n, n), UnitSquareMesh(n, n)]
+fields = [Field("u", family="Lagrange", degree=2, vector=True)]
 
 
 def get_solver(mesh_seq):
     def solver(index):
-        u, u_ = mesh_seq.fields["u"]
+        u, u_ = mesh_seq.field_functions["u"]
 
         # Define constants
         R = FunctionSpace(mesh_seq[index], "R", 0)
@@ -61,24 +64,20 @@ def get_initial_condition(mesh_seq):
     return {"u": Function(fs).interpolate(as_vector([sin(pi * x), 0]))}
 
 
-n = 32
-meshes = [UnitSquareMesh(n, n, diagonal="left"), UnitSquareMesh(n, n, diagonal="left")]
 end_time = 0.5
 dt = 1 / n
-
 num_subintervals = len(meshes)
 time_partition = TimePartition(
     end_time,
     num_subintervals,
     dt,
-    field_names,
+    fields,
     num_timesteps_per_export=2,
 )
 
 mesh_seq = MeshSeq(
     time_partition,
     meshes,
-    get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
     get_solver=get_solver,
 )
@@ -148,7 +147,7 @@ def adaptor(mesh_seq, solutions):
     num_elem = mesh_seq.count_elements()
     pyrint(f"fixed point iteration {iteration + 1}:")
     for i, (complexity, ndofs, nelem) in enumerate(
-        zip(complexities, num_dofs, num_elem)
+        zip(complexities, num_dofs, num_elem, strict=True)
     ):
         pyrint(
             f"  subinterval {i}, complexity: {complexity:4.0f}"

demos/burgers.py

@@ -24,26 +24,28 @@
 
 from goalie import *
 
-# In this problem, we have a single prognostic variable,
-# :math:`\mathbf u`. Its name is recorded in a list of
-# field names. ::
-
-field_names = ["u"]
-
-# For each such field, we need to be able to specify how to
-# build a :class:`FunctionSpace`, given some mesh. Since there
-# could be more than one field, function spaces are given as a
-# dictionary, indexed by the prognostic solution field names. ::
+# We begin by defining the two meshes of the unit sequare that we'd like to solve over.
+# For simplicity, we just use the same mesh twice: a :math:`32\times32` grid of the unit
+# square, with each grid-box divided into right-angled triangles. ::
 
+n = 32
+mesh = UnitSquareMesh(n, n)
 
-def get_function_spaces(mesh):
-    return {"u": VectorFunctionSpace(mesh, "CG", 2)}
+# In the Burgers problem, we have a single prognostic variable, :math:`\mathbf u`. Its
+# name and other metadata are recorded in a :class:`~.Field` object. One important piece
+# of metadata is the finite element used to define function spaces for the field (given
+# some mesh). This can be defined either using the :class:`finat.ufl.FiniteElement`
+# class, or using the same arguments as can be passed to
+# :class:`firedrake.functionspace.FunctionSpace` (e.g., `mesh`, `family`, `degree`). In
+# this case, we use a :math:`\mathbb{P}2` space so specify `family="Lagrange"` and
+# `degree=2`.Since Burgers is a vector equation, we need to specify `vector=True`. ::
 
+fields = [Field("u", family="Lagrange", degree=2, vector=True)]
 
 # 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.
+# by the :func:`get_function_spaces` function and are accessed via the
+# :attr:`field_functions` attribute of the :class:`MeshSeq`. This attribute provides a
+# dictionary of tuples containing the current and lagged solutions for each field.
 #
 # 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
@@ -64,7 +66,7 @@ def get_function_spaces(mesh):
 def get_solver(mesh_seq):
     def solver(index):
         # Get the current and lagged solutions
-        u, u_ = mesh_seq.fields["u"]
+        u, u_ = mesh_seq.field_functions["u"]
 
         # Define constants
         R = FunctionSpace(mesh_seq[index], "R", 0)
@@ -105,37 +107,29 @@ def get_initial_condition(mesh_seq):
     return {"u": Function(fs).interpolate(as_vector([sin(pi * x), 0]))}
 
 
-# Now that we have the above functions defined, we move onto the
-# concrete parts of the solver. To begin with, we require a
-# sequence of meshes, simulation end time and a timestep. ::
+# Now that we have the above functions defined, we need to define the time
+# discretisation used for the solver. To do this, we create a :class:`TimePartition` for
+# the problem with two subintervals. ::
 
-n = 32
-meshes = [
-    UnitSquareMesh(n, n),
-    UnitSquareMesh(n, n),
-]
 end_time = 0.5
 dt = 1 / n
-
-# These can be used to create a :class:`TimePartition` for the
-# problem with two subintervals. ::
-
-num_subintervals = len(meshes)
+num_subintervals = 2
 time_partition = TimePartition(
     end_time,
     num_subintervals,
     dt,
-    field_names,
+    fields,
     num_timesteps_per_export=2,
 )
 
-# Finally, we are able to construct a :class:`MeshSeq` and
-# solve Burgers equation over the meshes in sequence. ::
+# Finally, we are able to construct a :class:`~.MeshSeq` and solve Burgers equation over
+# the meshes in sequence. Note that the second argument can be either a list of meshes
+# or just a single mesh. If a single mesh is passed then this will be used for all
+# subintervals. ::
 
 mesh_seq = MeshSeq(
     time_partition,
-    meshes,
-    get_function_spaces=get_function_spaces,
+    mesh,
     get_initial_condition=get_initial_condition,
     get_solver=get_solver,
 )

demos/burgers1.py

@@ -6,32 +6,25 @@
 # automatic differentiation functionality in order to
 # automatically form and solve discrete adjoint problems.
 #
-# We always begin by importing Goalie. Adjoint mode is used
-# so that we have access to the :class:`AdjointMeshSeq` class.
-# ::
+# We always begin by importing Goalie. ::
 
 from firedrake import *
 
-from goalie_adjoint import *
+from goalie import *
 
-# For ease, the list of field names and functions for obtaining the
-# 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
-# ``ad_block_tag`` that matches the corresponding prognostic
-# variable name. ::
+# For ease, the list of fields and functions for obtaining the 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
+# ``ad_block_tag`` that matches the corresponding prognostic variable name. ::
 
-field_names = ["u"]
-
-
-def get_function_spaces(mesh):
-    return {"u": VectorFunctionSpace(mesh, "CG", 2)}
+n = 32
+mesh = UnitSquareMesh(n, n)
+fields = [Field("u", family="Lagrange", degree=2, vector=True)]
 
 
 def get_solver(mesh_seq):
     def solver(index):
-        u, u_ = mesh_seq.fields["u"]
+        u, u_ = mesh_seq.field_functions["u"]
 
         # Define constants
         R = FunctionSpace(mesh_seq[index], "R", 0)
@@ -83,26 +76,19 @@ def get_initial_condition(mesh_seq):
 
 def get_qoi(mesh_seq, i):
     def end_time_qoi():
-        u = mesh_seq.fields["u"][0]
+        u = mesh_seq.field_functions["u"][0]
         return inner(u, u) * ds(2)
 
     return end_time_qoi
 
 
-# Now that we have the above functions defined, we move onto the
-# concrete parts of the solver, which mimic the original demo. ::
+# Next, we define the :class:`~.TimePartition`. In cases where we only solve over a
+# single time subinterval (as in this demo), the partition is trivial and we can use the
+# :class:`~.TimeInterval` constructor, which requires fewer arguments. ::
 
-n = 32
-mesh = UnitSquareMesh(n, n)
 end_time = 0.5
 dt = 1 / n
-
-# Another requirement to solve the adjoint problem using
-# Goalie is a :class:`TimePartition`. In our case, there is a
-# single mesh, so the partition is trivial and we can use the
-# :class:`TimeInterval` constructor. ::
-
-time_partition = TimeInterval(end_time, dt, field_names, num_timesteps_per_export=2)
+time_partition = TimeInterval(end_time, dt, fields, num_timesteps_per_export=2)
 
 # Finally, we are able to construct an :class:`AdjointMeshSeq` and
 # thereby call its :meth:`solve_adjoint` method. This computes the QoI
@@ -112,7 +98,6 @@ def end_time_qoi():
 mesh_seq = AdjointMeshSeq(
     time_partition,
     mesh,
-    get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
     get_solver=get_solver,
     get_qoi=get_qoi,