From 04a42d1c9124424558f099f937650d0a6acee9f3 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Mon, 15 Apr 2024 19:01:40 +0000
Subject: [PATCH 01/17] Transfer mixed function spaces in indicate_errors

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

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index b9a32f56..e09adde7 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -8,6 +8,7 @@
 from .log import pyrint
 from firedrake import Function, FunctionSpace, MeshHierarchy, TransferManager
 from firedrake.petsc import PETSc
+from animate.interpolation import interpolate
 from collections.abc import Iterable
 import numpy as np
 import ufl
@@ -102,7 +103,7 @@ def _get_transfer_function(enrichment_method):
         if enrichment_method == "h":
             return TransferManager().prolong
         else:
-            return lambda source, target: target.interpolate(source)
+            return lambda source, target: interpolate(source, target)
 
     def _create_indicators(self):
         """

From fc2cf230c0e499bbb0011402a8fb44c8439bec82 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Mon, 15 Apr 2024 19:48:07 +0000
Subject: [PATCH 02/17] #152: Add mantle convection reference

---
 demos/demo_references.bib | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/demos/demo_references.bib b/demos/demo_references.bib
index 31fe2797..3aa488e3 100644
--- a/demos/demo_references.bib
+++ b/demos/demo_references.bib
@@ -16,3 +16,14 @@ @Book{Hundsdorfer:2003
   year={2003},
   publisher={Springer}
 }
+
+@article{Davies:2022,
+  title={Towards automatic finite-element methods for geodynamics via Firedrake},
+  author={Davies, D Rhodri and Kramer, Stephan C and Ghelichkhan, Sia and Gibson, Angus},
+  journal={Geoscientific Model Development},
+  volume={15},
+  number={13},
+  pages={5127--5166},
+  year={2022},
+  publisher={Copernicus Publications G{\"o}ttingen, Germany}
+}

From 8e45ed4e8fcba5b27cd0618f7b667186a213780b Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Mon, 15 Apr 2024 20:58:29 +0000
Subject: [PATCH 03/17] #153: Extract lagged fields in mixed cases

---
 goalie/adjoint.py  | 49 ++++++++++++++++++++++++++++++++++++++--------
 goalie/mesh_seq.py | 27 +++++++++++++++++++++----
 2 files changed, 64 insertions(+), 12 deletions(-)

diff --git a/goalie/adjoint.py b/goalie/adjoint.py
index b3c5a756..b94fee85 100644
--- a/goalie/adjoint.py
+++ b/goalie/adjoint.py
@@ -393,14 +393,47 @@ def wrapped_solver(subinterval, initial_condition_map, **kwargs):
                     if block.adj_sol is not None:
                         solutions.adjoint[i][j].assign(block.adj_sol)
 
-                    # Lagged forward solution comes from dependencies
-                    dep = self._dependency(field, i, block)
-                    if not self.steady and dep is not None:
-                        solutions.forward_old[i][j].assign(dep.saved_output)
-
-                    # Adjoint action also comes from dependencies
-                    if get_adj_values and dep is not None:
-                        solutions.adj_value[i][j].assign(dep.adj_value)
+                    if not self.steady:
+                        # Lagged solution comes from dependencies for unsteady fields
+                        if self.field_types[field] == "unsteady":
+                            dep = self._dependency(field, i, block)
+                            solutions.forward_old[i][j].assign(dep.saved_output)
+                            # Adjoint action also comes from dependencies
+                            if get_adj_values:
+                                solutions.adj_value[i][j].assign(dep.adj_value)
+                        # Lagged solution comes from previous block for steady fields
+                        if self.field_types[field] == "steady":
+                            if stride == 1:
+                                if j == 0:
+                                    if i == 0:
+                                        forward_old = self.initial_condition[field]
+                                    else:
+                                        forward_old = self._transfer(
+                                            solutions.forward[i - 1][-1], fs[i]
+                                        )
+                                else:
+                                    forward_old = solutions.forward[i][j - 1]
+                                if get_adj_values:
+                                    # TODO this does not consider the j==0 case
+                                    if j == num_solve_blocks - 1:
+                                        if i != num_subintervals - 1:
+                                            adj_value = self._transfer(
+                                                out.adj_value, fs[i + 1]
+                                            )
+                                            solutions.adj_value[i + 1][0].assign(
+                                                adj_value
+                                            )
+                                    else:
+                                        solutions.adj_value[i][j + 1].assign(
+                                            out.adj_value
+                                        )
+                            else:
+                                old_block = solve_blocks[solve_blocks.index(block) - 1]
+                                old_out = self._output(field, i, old_block)
+                                forward_old = old_out.saved_output
+                                if get_adj_values:
+                                    solutions.adj_value[i][j].assign(old_out.adj_value)
+                            solutions.forward_old[i][j].assign(forward_old)
 
                     # The adjoint solution at the 'next' timestep is determined from the
                     # adj_sol attribute of the next solve block
diff --git a/goalie/mesh_seq.py b/goalie/mesh_seq.py
index 8ad5087d..4703d8a5 100644
--- a/goalie/mesh_seq.py
+++ b/goalie/mesh_seq.py
@@ -760,10 +760,29 @@ def solve_forward(self, solver_kwargs={}):
                     if out is not None:
                         solutions.forward[i][j].assign(out.saved_output)
 
-                    # Lagged solution comes from dependencies
-                    dep = self._dependency(field, i, block)
-                    if not self.steady and dep is not None:
-                        solutions.forward_old[i][j].assign(dep.saved_output)
+                    if not self.steady:
+                        # Lagged solution comes from dependencies for unsteady fields
+                        if self.field_types[field] == "unsteady":
+                            dep = self._dependency(field, i, block)
+                            solutions.forward_old[i][j].assign(dep.saved_output)
+                        # Lagged solution comes from previous block for steady fields
+                        elif self.field_types[field] == "steady":
+                            if stride == 1:
+                                if j == 0:
+                                    if i == 0:
+                                        forward_old = self.initial_condition[field]
+                                    else:
+                                        forward_old = self._transfer(
+                                            solutions.forward[i - 1][-1], fs[i]
+                                        )
+                                else:
+                                    forward_old = solutions.forward[i][j - 1]
+                            else:
+                                old_block = solve_blocks[solve_blocks.index(block) - 1]
+                                old_out = self._output(field, i, old_block)
+                                if out is not None:
+                                    forward_old = old_out.saved_output
+                            solutions.forward_old[i][j].assign(forward_old)
 
             # Transfer the checkpoint between subintervals
             if i < num_subintervals - 1:

From ba70fdfb985c9a9eb8e1c95056cc9febd6c08dbe Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Mon, 15 Apr 2024 21:03:47 +0000
Subject: [PATCH 04/17] #153: Remove consideration of field_types in solution
 objects

---
 goalie/function_data.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/goalie/function_data.py b/goalie/function_data.py
index 9ad428ea..ce5b46f8 100644
--- a/goalie/function_data.py
+++ b/goalie/function_data.py
@@ -47,10 +47,10 @@ def _create_data(self):
                             ]
                             for i, fs in enumerate(self.function_spaces[field])
                         ]
-                        for label in self._label_dict[field_type]
+                        for label in self.labels
                     }
                 )
-                for field, field_type in zip(tp.fields, tp.field_types)
+                for field in tp.fields
             }
         )
 
@@ -201,7 +201,7 @@ def __init__(self, time_partition, meshes):
         :arg meshes: the list of meshes used to discretise the problem in space
         """
         self._label_dict = {
-            field_type: ("error_indicator",) for field_type in ("steady", "unsteady")
+            time_dep: ("error_indicator",) for time_dep in ("steady", "unsteady")
         }
         super().__init__(
             time_partition,

From 12232ad5e0fc9ee2e2675c7d6b0afd89e5ab81e8 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Mon, 15 Apr 2024 21:07:16 +0000
Subject: [PATCH 05/17] #152: Mantle convection demo

---
 demos/mantle_convection.py | 243 +++++++++++++++++++++++++++++++++++++
 1 file changed, 243 insertions(+)
 create mode 100644 demos/mantle_convection.py

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
new file mode 100644
index 00000000..e3a1cb82
--- /dev/null
+++ b/demos/mantle_convection.py
@@ -0,0 +1,243 @@
+# Mantle convection modelling
+#############################
+
+# In all demos that we have considered so far, the equations that we have solved all
+# involve a time derivative. However, in some cases, we may have a time-dependent
+# problem where an equation might not involve a time derivative, but is still
+# time-dependent in the sense that it depends on other fields that are time-dependent
+# and that have been previously solved for. An example of where this might happen is
+# in a continuous pressure projection approach, ice sheet and glaciological modelling,
+# mantle convection modelling, etc. In this demo, we illustrate how Goalie can be used
+# to solve such problems.
+
+# We consider the problem of a mantle convection in a 2D unit square domain. The
+# governing equations and Firedrake implementation are based on the 2D Cartesian
+# incompressible isoviscous case from :cite:`Davies:2022`. We refer the reader to the
+# paper for a detailed description of the problem and implementation. Here we
+# immediately present the governing equations involving a Stokes system and an energy
+# equation, which we solve for the velocity :math:`\mathbf{u}`, pressure :math:`p`, and
+# temperature :math:`T`:
+#
+# .. math::
+#    \begin{align}
+#        \nabla \cdot \mu \left[\nabla \mathbf{u} + (\nabla \mathbf{u})^T \right] - \nabla p + \mathrm{Ra}\,T\,\mathbf{k} &= 0, \\
+#        \frac{\partial T}{\partial t} \cdot \mathbf{u}\cdot\nabla T - \nabla \cdot (\kappa\nabla T) &= 0,
+#    \end{align}
+#
+# where :math:`\mu`, :math:`\kappa`, and :math:`\mathrm{Ra}` are constant viscosity,
+# thermal diffusivity, and Rayleigh number, respectively, and :math:`\mathbf{k}` is
+# the unit vector in the direction opposite to gravity. Note that while the Stokes
+# equations do not involve a time derivative, they are still time-dependent as they
+# depend on the temperature field :math:`T`, which is time-dependent.
+#
+# As always, we begin by importing Firedrake and Goalie, and defining constants.
+
+from firedrake import *
+from goalie_adjoint import *
+
+Ra, mu, kappa = Constant(1e4), Constant(1.0), Constant(1.0)
+k = Constant((0, 1))
+
+# The problem is solved simultaneously for the velocity :math:`\mathbf{u}` and pressure
+# :math:`p` using a *mixed* formulation, which was introduced in a `previous demo on
+# advection-diffusion reaction <./gray_scott.py.html>`__.
+
+fields = ["up", "T"]
+
+
+def get_function_spaces(mesh):
+    V = VectorFunctionSpace(mesh, "CG", 2)  # Velocity function space
+    W = FunctionSpace(mesh, "CG", 1)  # Pressure function space
+    Z = MixedFunctionSpace([V, W])  # Mixed function space for velocity and pressure
+    Q = FunctionSpace(mesh, "CG", 1)  # Temperature function space
+    return {"up": Z, "T": Q}
+
+
+# We must prescribe boundary and initial conditions to solve the problem. Note that
+# Goalie requires defining the initial condition for the mixed field ``up`` despite
+# the equations not involving a time derivative. In this case, the prescribed initial
+# condition should be understood as the *initial guess* for the solver.
+
+
+def get_bcs(mesh_seq):
+    def bcs(index):
+        Z = mesh_seq.function_spaces["up"][index]
+        Q = mesh_seq.function_spaces["T"][index]
+
+        # Boundary IDs
+        left, right, bottom, top = 1, 2, 3, 4
+
+        # Boundary conditions for velocity
+        bcvx = DirichletBC(Z.sub(0).sub(0), 0, sub_domain=(left, right))
+        bcvy = DirichletBC(Z.sub(0).sub(1), 0, sub_domain=(bottom, top))
+        # Boundary conditions for temperature
+        bctb = DirichletBC(Q, 1.0, sub_domain=bottom)
+        bctt = DirichletBC(Q, 0.0, sub_domain=top)
+
+        return {"up": [bcvx, bcvy], "T": [bctb, bctt]}
+
+    return bcs
+
+
+def get_initial_condition(mesh_seq):
+    x, y = SpatialCoordinate(mesh_seq[0])
+    T_init = Function(mesh_seq.function_spaces["T"][0])
+    T_init.interpolate(1.0 - y + 0.05 * cos(pi * x) * sin(pi * y))
+    up_init = Function(mesh_seq.function_spaces["up"][0])  # Zero by default
+    return {"up": up_init, "T": T_init}
+
+
+# The weak forms for the Stokes and energy equations are defined as follows. Note that
+# the lagged solution ``up_`` is not used in compiling the form, but is still by
+# internal Goalie machinery.
+
+
+def get_form(mesh_seq):
+    def form(index, sols):
+        dt = mesh_seq.time_partition.timesteps[index]
+
+        up, up_ = sols["up"]
+        u, p = split(up)
+        v, w = TestFunctions(mesh_seq.function_spaces["up"][index])
+
+        T, T_ = sols["T"]
+        q = TestFunction(mesh_seq.function_spaces["T"][index])
+
+        # Crank-Nicolson time discretisation for temperature
+        Ttheta = 0.5 * (T + T_)
+
+        # Stokes equations
+        stress = 2 * mu * sym(grad(u))
+        F_up = (
+            inner(grad(v), stress) * dx
+            - div(v) * p * dx
+            - (dot(v, k) * Ra * Ttheta) * dx
+        )
+        F_up += -w * div(u) * dx  # Continuity equation
+
+        # Energy equation
+        F_T = (
+            q * (T - T_) / dt * dx
+            + q * dot(u, grad(Ttheta)) * dx
+            + dot(grad(q), kappa * grad(Ttheta)) * dx
+        )
+
+        return {"up": F_up, "T": F_T}
+
+    return form
+
+
+def get_solver(mesh_seq):
+    def solver(index, ic):
+        Z = mesh_seq.function_spaces["up"][index]
+        up = Function(Z, name="up")
+        up_ = Function(Z, name="up_old")
+
+        Q = mesh_seq.function_spaces["T"][index]
+        T = Function(Q, name="T")
+        T_ = Function(Q, name="T_old")
+        T_.assign(ic["T"])
+
+        # Dictionary of weak forms and boundary conditions for both fields
+        F = mesh_seq.form(index, {"up": (up, up_), "T": (T, T_)})
+        bcs = mesh_seq.bcs(index)
+
+        # Solver parameters dictionary
+        solver_parameters = {
+            "mat_type": "aij",
+            "snes_type": "ksponly",
+            "ksp_type": "preonly",
+            "pc_type": "lu",
+            "pc_factor_mat_solver_type": "mumps",
+        }
+
+        nlvp_up = NonlinearVariationalProblem(F["up"], up, bcs=bcs["up"])
+        nlvs_up = NonlinearVariationalSolver(
+            nlvp_up,
+            solver_parameters=solver_parameters,
+            ad_block_tag="up",
+        )
+        nlvp_T = NonlinearVariationalProblem(F["T"], T, bcs=bcs["T"])
+        nlvs_T = NonlinearVariationalSolver(
+            nlvp_T,
+            solver_parameters=solver_parameters,
+            ad_block_tag="T",
+        )
+
+        # Time integrate from t_start to t_end
+        P = mesh_seq.time_partition
+        num_timesteps = P.num_timesteps_per_subinterval[index]
+        for _ in range(num_timesteps):
+            nlvs_up.solve()
+            nlvs_T.solve()
+            T_.assign(T)
+        return {"up": up, "T": T}
+
+    return solver
+
+
+# Finally, we define the quantity of interest (QoI) as the square of the velocity field
+# :math:`\mathbf{u}`. The QoI will be evaluated at the final time step.
+
+
+def get_qoi(mesh_seq, sols, i):
+    def qoi():
+        up = sols["up"]
+        u, _ = split(up)
+        return dot(u, u) * dx
+
+    return qoi
+
+
+# We can now create a GoalOrientedMeshSeq object and compute the error indicators.
+# For demonstration purposes, we only consider two subintervals and a low number of
+# timesteps.
+
+num_subintervals = 2
+meshes = [UnitSquareMesh(32, 32, quadrilateral=True) for _ in range(num_subintervals)]
+
+dt = 1e-3
+num_timesteps = 40
+end_time = dt * num_timesteps
+dt_per_export = [10 for _ in range(num_subintervals)]
+
+# To account for the lack of time derivative in the Stokes equations, we use the
+# ``field_types`` argument of the ``TimePartition`` object to specify that the ``up``
+# field is *steady* (i.e. without a time derivative) and that the ``T`` field is
+# *unsteady* (i.e. involves a time derivative). The order in ``field_types`` must
+# match the order of the fields in the ``fields`` list (see above).
+
+time_partition = TimePartition(
+    end_time,
+    num_subintervals,
+    dt,
+    fields,
+    num_timesteps_per_export=dt_per_export,
+    field_types=["steady", "unsteady"],
+)
+
+mesh_seq = GoalOrientedMeshSeq(
+    time_partition,
+    meshes,
+    get_function_spaces=get_function_spaces,
+    get_initial_condition=get_initial_condition,
+    get_bcs=get_bcs,
+    get_form=get_form,
+    get_solver=get_solver,
+    get_qoi=get_qoi,
+    qoi_type="end_time",
+)
+
+solutions, indicators = mesh_seq.indicate_errors()
+
+# We can plot the error indicator fields for exported timesteps using the in-built
+# plotting function ``plot_indicator_snapshots``.
+
+fig, axes, tcs = plot_indicator_snapshots(indicators, time_partition, "T", levels=50)
+fig.savefig("mantle_convection.jpg")
+
+# .. figure:: burgers-ee.jpg
+#    :figwidth: 90%
+#    :align: center
+#
+# This demo can also be accessed as a `Python script <mantle_convection.py>`__.

From 196f0c341673a1527d25d14d59aef553902fe8a2 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 16 Apr 2024 22:53:57 +0000
Subject: [PATCH 06/17] #55: Update all variables in indicate_errors

---
 .flake8                    |  2 ++
 demos/mantle_convection.py |  7 +++++++
 goalie/go_mesh_seq.py      | 19 +++++++++++--------
 3 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/.flake8 b/.flake8
index 2ea09c6e..125cc371 100644
--- a/.flake8
+++ b/.flake8
@@ -8,6 +8,8 @@ ignore =
 	E402,
 	# Do not use variables named "I", "O", or "l"
 	E741,
+	# whitespace before ':'
+	E203,
 	# Unable to detect undefined names
 	F403,
 	# Name may be undefined, or defined from star imports
diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index e3a1cb82..8724bbe2 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -127,6 +127,13 @@ def form(index, sols):
     return form
 
 
+# In the solver, it is important to solve for the fields in the order in which they are
+# defined in the ``fields`` list. This is to ensure that the error indicators are
+# computed correctly. Therefore, in the time integration loop, we first solve the
+# Stokes equations for the velocity and pressure fields, and then solve the energy
+# equation for the temperature field.
+
+
 def get_solver(mesh_seq):
     def solver(index, ic):
         Z = mesh_seq.function_spaces["up"][index]
diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index e09adde7..54fcc024 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -185,17 +185,20 @@ def indicate_errors(
                     f", not type '{type(forms)}'."
                 )
 
-            # Loop over each strongly coupled field
-            for f in self.fields:
-                # Loop over each timestep
-                solutions = self.solutions.extract(layout="field")
-                indicators = self.indicators.extract(layout="field")
-                for j in range(len(solutions[f]["forward"][i])):
-                    # Update fields
+            # Loop over each timestep
+            for j in range(self.time_partition.num_exports_per_subinterval[i] - 1):
+                # Loop over each strongly coupled field
+                for f in self.fields:
+                    # Transfer solutions associated to the current field
                     transfer(self.solutions[f][FWD][i][j], u[f])
                     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])
+                    # Transfer solutions associated to the fields solved next
+                    f_idx = self.fields.index(f)
+                    if len(self.fields) > 1 and f_idx < len(self.fields) - 1:
+                        for f_next in self.fields[f_idx + 1 :]:
+                            transfer(self.solutions[f_next][FWD_OLD][i][j], u[f_next])
 
                     # Combine adjoint solutions as appropriate
                     u_star[f].assign(0.5 * (u_star[f] + u_star_next[f]))
@@ -214,7 +217,7 @@ def indicate_errors(
                     # Transfer back to the base space
                     indi = self._transfer(indi_e, P0_spaces[i])
                     indi.interpolate(abs(indi))
-                    indicators[f][i][j].interpolate(ufl.max_value(indi, 1.0e-16))
+                    self.indicators[f][i][j].interpolate(ufl.max_value(indi, 1.0e-16))
 
         return self.solutions, self.indicators
 

From 89861bcdd0160238225148199039cea0e003504a Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 23 Jun 2024 06:21:09 +0000
Subject: [PATCH 07/17] Merge branch 'main' into 152_mixed_demo

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

diff --git a/goalie/go_mesh_seq.py b/goalie/go_mesh_seq.py
index b7939c4e..28c25772 100644
--- a/goalie/go_mesh_seq.py
+++ b/goalie/go_mesh_seq.py
@@ -204,11 +204,6 @@ def indicate_errors(
                     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])
-                    # Transfer solutions associated to the fields solved next
-                    f_idx = self.fields.index(f)
-                    if len(self.fields) > 1 and f_idx < len(self.fields) - 1:
-                        for f_next in self.fields[f_idx + 1 :]:
-                            transfer(self.solutions[f_next][FWD_OLD][i][j], u[f_next])
 
                     # Combine adjoint solutions as appropriate
                     u_star[f].assign(0.5 * (u_star[f] + u_star_next[f]))

From 203ed279b61422a6e8efbef3a3805ab59ca7291c Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 23 Jun 2024 07:22:58 +0000
Subject: [PATCH 08/17] #195: Account for mixed case

---
 goalie/mesh_seq.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/goalie/mesh_seq.py b/goalie/mesh_seq.py
index 8ea3d4d6..ed18c005 100644
--- a/goalie/mesh_seq.py
+++ b/goalie/mesh_seq.py
@@ -551,7 +551,12 @@ def _solve_forward(self, update_solutions=True, solver_kwargs=None):
             if i < num_subintervals - 1:
                 checkpoint = AttrDict(
                     {
-                        field: self._transfer(self.fields[field][0], fs[i + 1])
+                        field: self._transfer(
+                            self.fields[field]
+                            if self.field_types[field] == "steady"
+                            else self.fields[field][0],
+                            fs[i + 1],
+                        )
                         for field, fs in self._fs.items()
                     }
                 )

From 4b8f7d36fe1dc293d45e2d7bc407ffdaace87846 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 23 Jun 2024 07:25:52 +0000
Subject: [PATCH 09/17] Revert "#195: Account for mixed case"

This reverts commit 203ed279b61422a6e8efbef3a3805ab59ca7291c.
---
 goalie/mesh_seq.py | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/goalie/mesh_seq.py b/goalie/mesh_seq.py
index ed18c005..8ea3d4d6 100644
--- a/goalie/mesh_seq.py
+++ b/goalie/mesh_seq.py
@@ -551,12 +551,7 @@ def _solve_forward(self, update_solutions=True, solver_kwargs=None):
             if i < num_subintervals - 1:
                 checkpoint = AttrDict(
                     {
-                        field: self._transfer(
-                            self.fields[field]
-                            if self.field_types[field] == "steady"
-                            else self.fields[field][0],
-                            fs[i + 1],
-                        )
+                        field: self._transfer(self.fields[field][0], fs[i + 1])
                         for field, fs in self._fs.items()
                     }
                 )

From 32638d2e5aea15edc7b7d547d5f3c7933383987c Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 25 Jun 2024 19:19:15 +0000
Subject: [PATCH 10/17] #152: Update demo

---
 demos/mantle_convection.py | 83 ++++++++++++++++----------------------
 1 file changed, 34 insertions(+), 49 deletions(-)

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index 8724bbe2..65c6776a 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -33,6 +33,7 @@
 # As always, we begin by importing Firedrake and Goalie, and defining constants.
 
 from firedrake import *
+
 from goalie_adjoint import *
 
 Ra, mu, kappa = Constant(1e4), Constant(1.0), Constant(1.0)
@@ -53,30 +54,10 @@ def get_function_spaces(mesh):
     return {"up": Z, "T": Q}
 
 
-# We must prescribe boundary and initial conditions to solve the problem. Note that
-# Goalie requires defining the initial condition for the mixed field ``up`` despite
-# the equations not involving a time derivative. In this case, the prescribed initial
-# condition should be understood as the *initial guess* for the solver.
-
-
-def get_bcs(mesh_seq):
-    def bcs(index):
-        Z = mesh_seq.function_spaces["up"][index]
-        Q = mesh_seq.function_spaces["T"][index]
-
-        # Boundary IDs
-        left, right, bottom, top = 1, 2, 3, 4
-
-        # Boundary conditions for velocity
-        bcvx = DirichletBC(Z.sub(0).sub(0), 0, sub_domain=(left, right))
-        bcvy = DirichletBC(Z.sub(0).sub(1), 0, sub_domain=(bottom, top))
-        # Boundary conditions for temperature
-        bctb = DirichletBC(Q, 1.0, sub_domain=bottom)
-        bctt = DirichletBC(Q, 0.0, sub_domain=top)
-
-        return {"up": [bcvx, bcvy], "T": [bctb, bctt]}
-
-    return bcs
+# We must initial conditions to solve the problem. Note that we define the initial
+# condition for the mixed field ``up`` despite the equations not involving a time
+# derivative. In this case, the prescribed initial condition should be understood as the
+# *initial guess* for the solver.
 
 
 def get_initial_condition(mesh_seq):
@@ -88,19 +69,16 @@ def get_initial_condition(mesh_seq):
 
 
 # The weak forms for the Stokes and energy equations are defined as follows. Note that
-# the lagged solution ``up_`` is not used in compiling the form, but is still by
-# internal Goalie machinery.
+# the mixed field ``up`` does not have a lagged term.
 
 
 def get_form(mesh_seq):
-    def form(index, sols):
-        dt = mesh_seq.time_partition.timesteps[index]
-
-        up, up_ = sols["up"]
+    def form(index):
+        up = mesh_seq.fields["up"]
         u, p = split(up)
         v, w = TestFunctions(mesh_seq.function_spaces["up"][index])
 
-        T, T_ = sols["T"]
+        T, T_ = mesh_seq.fields["T"]
         q = TestFunction(mesh_seq.function_spaces["T"][index])
 
         # Crank-Nicolson time discretisation for temperature
@@ -116,6 +94,7 @@ def form(index, sols):
         F_up += -w * div(u) * dx  # Continuity equation
 
         # Energy equation
+        dt = mesh_seq.time_partition.timesteps[index]
         F_T = (
             q * (T - T_) / dt * dx
             + q * dot(u, grad(Ttheta)) * dx
@@ -135,19 +114,26 @@ def form(index, sols):
 
 
 def get_solver(mesh_seq):
-    def solver(index, ic):
+    def solver(index):
         Z = mesh_seq.function_spaces["up"][index]
-        up = Function(Z, name="up")
-        up_ = Function(Z, name="up_old")
-
         Q = mesh_seq.function_spaces["T"][index]
-        T = Function(Q, name="T")
-        T_ = Function(Q, name="T_old")
-        T_.assign(ic["T"])
+
+        up = mesh_seq.fields["up"]
+        T, T_ = mesh_seq.fields["T"]
 
         # Dictionary of weak forms and boundary conditions for both fields
-        F = mesh_seq.form(index, {"up": (up, up_), "T": (T, T_)})
-        bcs = mesh_seq.bcs(index)
+        F = mesh_seq.form(index)
+
+        # Boundary IDs
+        left, right, bottom, top = 1, 2, 3, 4
+        # Boundary conditions for velocity
+        bcux = DirichletBC(Z.sub(0).sub(0), 0, sub_domain=(left, right))
+        bcuy = DirichletBC(Z.sub(0).sub(1), 0, sub_domain=(bottom, top))
+        bcs_up = [bcux, bcuy]
+        # Boundary conditions for temperature
+        bctb = DirichletBC(Q, 1.0, sub_domain=bottom)
+        bctt = DirichletBC(Q, 0.0, sub_domain=top)
+        bcs_T = [bctb, bctt]
 
         # Solver parameters dictionary
         solver_parameters = {
@@ -158,27 +144,26 @@ def solver(index, ic):
             "pc_factor_mat_solver_type": "mumps",
         }
 
-        nlvp_up = NonlinearVariationalProblem(F["up"], up, bcs=bcs["up"])
+        nlvp_up = NonlinearVariationalProblem(F["up"], up, bcs=bcs_up)
         nlvs_up = NonlinearVariationalSolver(
             nlvp_up,
             solver_parameters=solver_parameters,
             ad_block_tag="up",
         )
-        nlvp_T = NonlinearVariationalProblem(F["T"], T, bcs=bcs["T"])
+        nlvp_T = NonlinearVariationalProblem(F["T"], T, bcs=bcs_T)
         nlvs_T = NonlinearVariationalSolver(
             nlvp_T,
             solver_parameters=solver_parameters,
             ad_block_tag="T",
         )
 
-        # Time integrate from t_start to t_end
-        P = mesh_seq.time_partition
-        num_timesteps = P.num_timesteps_per_subinterval[index]
+        # Time integrate over the subinterval
+        num_timesteps = mesh_seq.time_partition.num_timesteps_per_subinterval[index]
         for _ in range(num_timesteps):
             nlvs_up.solve()
             nlvs_T.solve()
+            yield
             T_.assign(T)
-        return {"up": up, "T": T}
 
     return solver
 
@@ -187,9 +172,9 @@ def solver(index, ic):
 # :math:`\mathbf{u}`. The QoI will be evaluated at the final time step.
 
 
-def get_qoi(mesh_seq, sols, i):
+def get_qoi(mesh_seq, i):
     def qoi():
-        up = sols["up"]
+        up = mesh_seq.fields["up"]
         u, _ = split(up)
         return dot(u, u) * dx
 
@@ -228,11 +213,11 @@ def qoi():
     meshes,
     get_function_spaces=get_function_spaces,
     get_initial_condition=get_initial_condition,
-    get_bcs=get_bcs,
     get_form=get_form,
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
+    transfer_method="interpolate",
 )
 
 solutions, indicators = mesh_seq.indicate_errors()

From 9b8f0dca98f6273d82b2d4d8a77519f028c63ce5 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Tue, 25 Jun 2024 19:19:45 +0000
Subject: [PATCH 11/17] #152: Project in demo

---
 demos/mantle_convection.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index 65c6776a..ab354329 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -217,7 +217,6 @@ def qoi():
     get_solver=get_solver,
     get_qoi=get_qoi,
     qoi_type="end_time",
-    transfer_method="interpolate",
 )
 
 solutions, indicators = mesh_seq.indicate_errors()

From 85944580bcabdf633d2f0f12a9d9ec3b34dc5961 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 24 Nov 2024 11:42:48 +0000
Subject: [PATCH 12/17] Revert changes in adjoint.py

---
 goalie/adjoint.py | 49 ++++++++---------------------------------------
 1 file changed, 8 insertions(+), 41 deletions(-)

diff --git a/goalie/adjoint.py b/goalie/adjoint.py
index 15ca6982..95e8a4a7 100644
--- a/goalie/adjoint.py
+++ b/goalie/adjoint.py
@@ -565,47 +565,14 @@ def wrapped_solver(subinterval, initial_condition_map, **kwargs):
                     if block.adj_sol is not None:
                         solutions.adjoint[i][j].assign(block.adj_sol)
 
-                    if not self.steady:
-                        # Lagged solution comes from dependencies for unsteady fields
-                        if self.field_types[field] == "unsteady":
-                            dep = self._dependency(field, i, block)
-                            solutions.forward_old[i][j].assign(dep.saved_output)
-                            # Adjoint action also comes from dependencies
-                            if get_adj_values:
-                                solutions.adj_value[i][j].assign(dep.adj_value)
-                        # Lagged solution comes from previous block for steady fields
-                        if self.field_types[field] == "steady":
-                            if stride == 1:
-                                if j == 0:
-                                    if i == 0:
-                                        forward_old = self.initial_condition[field]
-                                    else:
-                                        forward_old = self._transfer(
-                                            solutions.forward[i - 1][-1], fs[i]
-                                        )
-                                else:
-                                    forward_old = solutions.forward[i][j - 1]
-                                if get_adj_values:
-                                    # TODO this does not consider the j==0 case
-                                    if j == num_solve_blocks - 1:
-                                        if i != num_subintervals - 1:
-                                            adj_value = self._transfer(
-                                                out.adj_value, fs[i + 1]
-                                            )
-                                            solutions.adj_value[i + 1][0].assign(
-                                                adj_value
-                                            )
-                                    else:
-                                        solutions.adj_value[i][j + 1].assign(
-                                            out.adj_value
-                                        )
-                            else:
-                                old_block = solve_blocks[solve_blocks.index(block) - 1]
-                                old_out = self._output(field, i, old_block)
-                                forward_old = old_out.saved_output
-                                if get_adj_values:
-                                    solutions.adj_value[i][j].assign(old_out.adj_value)
-                            solutions.forward_old[i][j].assign(forward_old)
+                    # Lagged forward solution comes from dependencies
+                    dep = self._dependency(field, i, block)
+                    if not self.steady and dep is not None:
+                        solutions.forward_old[i][j].assign(dep.saved_output)
+
+                    # Adjoint action also comes from dependencies
+                    if get_adj_values and dep is not None:
+                        solutions.adj_value[i][j].assign(dep.adj_value)
 
                     # The adjoint solution at the 'next' timestep is determined from the
                     # adj_sol attribute of the next solve block

From d4d5a3392c44b54bf86a629d0449aa5923976a31 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 24 Nov 2024 11:56:22 +0000
Subject: [PATCH 13/17] #152: Do not use get_form

---
 demos/mantle_convection.py | 79 ++++++++++----------------------------
 1 file changed, 21 insertions(+), 58 deletions(-)

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index ab354329..34d697bb 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -68,23 +68,24 @@ def get_initial_condition(mesh_seq):
     return {"up": up_init, "T": T_init}
 
 
-# The weak forms for the Stokes and energy equations are defined as follows. Note that
-# the mixed field ``up`` does not have a lagged term.
+# In the solver, weak forms for the Stokes and energy equations are defined as follows.
+# Note that the mixed field ``up`` does not have a lagged term.
 
 
-def get_form(mesh_seq):
-    def form(index):
+def get_solver(mesh_seq):
+    def solver(index):
+        Z = mesh_seq.function_spaces["up"][index]
+        Q = mesh_seq.function_spaces["T"][index]
+
         up = mesh_seq.fields["up"]
         u, p = split(up)
-        v, w = TestFunctions(mesh_seq.function_spaces["up"][index])
-
         T, T_ = mesh_seq.fields["T"]
-        q = TestFunction(mesh_seq.function_spaces["T"][index])
 
         # Crank-Nicolson time discretisation for temperature
         Ttheta = 0.5 * (T + T_)
 
-        # Stokes equations
+        # Variational problem for the Stokes equations
+        v, w = TestFunctions(mesh_seq.function_spaces["up"][index])
         stress = 2 * mu * sym(grad(u))
         F_up = (
             inner(grad(v), stress) * dx
@@ -93,37 +94,14 @@ def form(index):
         )
         F_up += -w * div(u) * dx  # Continuity equation
 
-        # Energy equation
-        dt = mesh_seq.time_partition.timesteps[index]
+        # Variational problem for the energy equation
+        q = TestFunction(mesh_seq.function_spaces["T"][index])
         F_T = (
             q * (T - T_) / dt * dx
             + q * dot(u, grad(Ttheta)) * dx
             + dot(grad(q), kappa * grad(Ttheta)) * dx
         )
 
-        return {"up": F_up, "T": F_T}
-
-    return form
-
-
-# In the solver, it is important to solve for the fields in the order in which they are
-# defined in the ``fields`` list. This is to ensure that the error indicators are
-# computed correctly. Therefore, in the time integration loop, we first solve the
-# Stokes equations for the velocity and pressure fields, and then solve the energy
-# equation for the temperature field.
-
-
-def get_solver(mesh_seq):
-    def solver(index):
-        Z = mesh_seq.function_spaces["up"][index]
-        Q = mesh_seq.function_spaces["T"][index]
-
-        up = mesh_seq.fields["up"]
-        T, T_ = mesh_seq.fields["T"]
-
-        # Dictionary of weak forms and boundary conditions for both fields
-        F = mesh_seq.form(index)
-
         # Boundary IDs
         left, right, bottom, top = 1, 2, 3, 4
         # Boundary conditions for velocity
@@ -144,13 +122,13 @@ def solver(index):
             "pc_factor_mat_solver_type": "mumps",
         }
 
-        nlvp_up = NonlinearVariationalProblem(F["up"], up, bcs=bcs_up)
+        nlvp_up = NonlinearVariationalProblem(F_up, up, bcs=bcs_up)
         nlvs_up = NonlinearVariationalSolver(
             nlvp_up,
             solver_parameters=solver_parameters,
             ad_block_tag="up",
         )
-        nlvp_T = NonlinearVariationalProblem(F["T"], T, bcs=bcs_T)
+        nlvp_T = NonlinearVariationalProblem(F_T, T, bcs=bcs_T)
         nlvs_T = NonlinearVariationalSolver(
             nlvp_T,
             solver_parameters=solver_parameters,
@@ -168,20 +146,7 @@ def solver(index):
     return solver
 
 
-# Finally, we define the quantity of interest (QoI) as the square of the velocity field
-# :math:`\mathbf{u}`. The QoI will be evaluated at the final time step.
-
-
-def get_qoi(mesh_seq, i):
-    def qoi():
-        up = mesh_seq.fields["up"]
-        u, _ = split(up)
-        return dot(u, u) * dx
-
-    return qoi
-
-
-# We can now create a GoalOrientedMeshSeq object and compute the error indicators.
+# We can now create a MeshSeq object and solve the forward problem.
 # For demonstration purposes, we only consider two subintervals and a low number of
 # timesteps.
 
@@ -197,7 +162,7 @@ def qoi():
 # ``field_types`` argument of the ``TimePartition`` object to specify that the ``up``
 # field is *steady* (i.e. without a time derivative) and that the ``T`` field is
 # *unsteady* (i.e. involves a time derivative). The order in ``field_types`` must
-# match the order of the fields in the ``fields`` list (see above).
+# match the order of the fields in the ``fields`` list above.
 
 time_partition = TimePartition(
     end_time,
@@ -208,26 +173,24 @@ def qoi():
     field_types=["steady", "unsteady"],
 )
 
-mesh_seq = GoalOrientedMeshSeq(
+mesh_seq = MeshSeq(
     time_partition,
     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",
+    transfer_method="interpolate",
 )
 
-solutions, indicators = mesh_seq.indicate_errors()
+solutions = mesh_seq.solve_forward()
 
-# We can plot the error indicator fields for exported timesteps using the in-built
+# We can plot the temperature fields for exported timesteps using the in-built
 # plotting function ``plot_indicator_snapshots``.
 
-fig, axes, tcs = plot_indicator_snapshots(indicators, time_partition, "T", levels=50)
+fig, axes, tcs = plot_snapshots(solutions, time_partition, "T", "forward", levels=25)
 fig.savefig("mantle_convection.jpg")
 
-# .. figure:: burgers-ee.jpg
+# .. figure:: mantle_convection.jpg
 #    :figwidth: 90%
 #    :align: center
 #

From 453416ae8f48be0720b15468ab182b26a0d56ec9 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 24 Nov 2024 11:56:47 +0000
Subject: [PATCH 14/17] #152: Consider steadiness per field

---
 goalie/mesh_seq.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/goalie/mesh_seq.py b/goalie/mesh_seq.py
index 1a89eb46..16496639 100644
--- a/goalie/mesh_seq.py
+++ b/goalie/mesh_seq.py
@@ -493,7 +493,7 @@ def _solve_forward(self, update_solutions=True, solver_kwargs=None):
                         next(solver_gen)
                     # Update the solution data
                     for field, sol in self.fields.items():
-                        if not self.steady:
+                        if not self.field_types[field] == "steady":
                             assert isinstance(sol, tuple)
                             solutions[field].forward[i][j].assign(sol[0])
                             solutions[field].forward_old[i][j].assign(sol[1])

From a335f1f25350848ff153abaa3ed4ecd93c208992 Mon Sep 17 00:00:00 2001
From: ddundo <davord.no@gmail.com>
Date: Sun, 24 Nov 2024 12:19:10 +0000
Subject: [PATCH 15/17] #152: Minor fixes

---
 demos/mantle_convection.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index 34d697bb..6c5a74ab 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -34,7 +34,7 @@
 
 from firedrake import *
 
-from goalie_adjoint import *
+from goalie import *
 
 Ra, mu, kappa = Constant(1e4), Constant(1.0), Constant(1.0)
 k = Constant((0, 1))
@@ -54,7 +54,7 @@ def get_function_spaces(mesh):
     return {"up": Z, "T": Q}
 
 
-# We must initial conditions to solve the problem. Note that we define the initial
+# We must set initial conditions to solve the problem. Note that we define the initial
 # condition for the mixed field ``up`` despite the equations not involving a time
 # derivative. In this case, the prescribed initial condition should be understood as the
 # *initial guess* for the solver.
@@ -185,7 +185,7 @@ def solver(index):
 solutions = mesh_seq.solve_forward()
 
 # We can plot the temperature fields for exported timesteps using the in-built
-# plotting function ``plot_indicator_snapshots``.
+# plotting function ``plot_snapshots``.
 
 fig, axes, tcs = plot_snapshots(solutions, time_partition, "T", "forward", levels=25)
 fig.savefig("mantle_convection.jpg")

From f753e0d3cfc2b5c631ad413d5d7158b6b1705fc0 Mon Sep 17 00:00:00 2001
From: Davor Dundovic <33790330+ddundo@users.noreply.github.com>
Date: Sun, 24 Nov 2024 21:14:36 +0100
Subject: [PATCH 16/17] #164 Apply suggestions from code review

Co-authored-by: Joe Wallwork <22053413+jwallwork23@users.noreply.github.com>
---
 demos/mantle_convection.py | 29 ++++++++++++++---------------
 1 file changed, 14 insertions(+), 15 deletions(-)

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index 6c5a74ab..ed56fe9d 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -1,20 +1,19 @@
 # Mantle convection modelling
-#############################
+============================
 
 # In all demos that we have considered so far, the equations that we have solved all
-# involve a time derivative. However, in some cases, we may have a time-dependent
-# problem where an equation might not involve a time derivative, but is still
-# time-dependent in the sense that it depends on other fields that are time-dependent
-# and that have been previously solved for. An example of where this might happen is
-# in a continuous pressure projection approach, ice sheet and glaciological modelling,
-# mantle convection modelling, etc. In this demo, we illustrate how Goalie can be used
+# involve a time derivative. Those are clearly *time-dependent* equations. However,
+# time-dependent equations need not involve a time derivative. For example, they might
+# include fields that vary with time. Examples of where this might happen are
+# in continuous pressure projection approaches, ice sheet and glaciological modelling, and
+# mantle convection modelling. In this demo, we illustrate how Goalie can be used
 # to solve such problems.
 
 # We consider the problem of a mantle convection in a 2D unit square domain. The
 # governing equations and Firedrake implementation are based on the 2D Cartesian
 # incompressible isoviscous case from :cite:`Davies:2022`. We refer the reader to the
 # paper for a detailed description of the problem and implementation. Here we
-# immediately present the governing equations involving a Stokes system and an energy
+# present the governing equations involving a Stokes system and an energy
 # equation, which we solve for the velocity :math:`\mathbf{u}`, pressure :math:`p`, and
 # temperature :math:`T`:
 #
@@ -47,15 +46,15 @@
 
 
 def get_function_spaces(mesh):
-    V = VectorFunctionSpace(mesh, "CG", 2)  # Velocity function space
-    W = FunctionSpace(mesh, "CG", 1)  # Pressure function space
-    Z = MixedFunctionSpace([V, W])  # Mixed function space for velocity and pressure
-    Q = FunctionSpace(mesh, "CG", 1)  # Temperature function space
+    V = VectorFunctionSpace(mesh, "CG", 2, name="velocity")
+    W = FunctionSpace(mesh, "CG", 1, name="pressure")
+    Z = MixedFunctionSpace([V, W], name="velocity-pressure")
+    Q = FunctionSpace(mesh, "CG", 1, name="temperature")
     return {"up": Z, "T": Q}
 
 
 # We must set initial conditions to solve the problem. Note that we define the initial
-# condition for the mixed field ``up`` despite the equations not involving a time
+# condition for the mixed field ``"up"`` despite the equations not involving a time
 # derivative. In this case, the prescribed initial condition should be understood as the
 # *initial guess* for the solver.
 
@@ -69,7 +68,7 @@ def get_initial_condition(mesh_seq):
 
 
 # In the solver, weak forms for the Stokes and energy equations are defined as follows.
-# Note that the mixed field ``up`` does not have a lagged term.
+# Note that the mixed field ``"up"`` does not have a lagged term.
 
 
 def get_solver(mesh_seq):
@@ -159,7 +158,7 @@ def solver(index):
 dt_per_export = [10 for _ in range(num_subintervals)]
 
 # To account for the lack of time derivative in the Stokes equations, we use the
-# ``field_types`` argument of the ``TimePartition`` object to specify that the ``up``
+# ``field_types`` argument of the ``TimePartition`` object to specify that the ``"up"``
 # field is *steady* (i.e. without a time derivative) and that the ``T`` field is
 # *unsteady* (i.e. involves a time derivative). The order in ``field_types`` must
 # match the order of the fields in the ``fields`` list above.

From 0a8f0f8d4dd5945ac32613cea49f3c52b0430b58 Mon Sep 17 00:00:00 2001
From: Davor Dundovic <33790330+ddundo@users.noreply.github.com>
Date: Sun, 24 Nov 2024 21:22:05 +0100
Subject: [PATCH 17/17] #152: Minor formatting fix

---
 demos/mantle_convection.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/demos/mantle_convection.py b/demos/mantle_convection.py
index ed56fe9d..114a8115 100644
--- a/demos/mantle_convection.py
+++ b/demos/mantle_convection.py
@@ -1,5 +1,5 @@
 # Mantle convection modelling
-============================
+# ===========================
 
 # In all demos that we have considered so far, the equations that we have solved all
 # involve a time derivative. Those are clearly *time-dependent* equations. However,