Adjoint inverse reconstruction¶

Introduction¶

In this tutorial, we will demonstrate how to perform an inversion to recover the initial temperature field of an idealised mantle convection simulation using G-ADOPT. This tutorial is published as the first synthetic experiment in Ghelichkhan et al. (2024). The full inversion showcased in the publication involves a total number of 80 timesteps. For the tutorial here we start with only 5 timesteps to go through the basics.

The tutorial involves a twin experiment, where we assess the performance of the inversion scheme by inverting the initial state of a synthetic reference simulation, known as the "Reference Twin". To create this reference twin, we run a forward mantle convection simulation and record all relevant fields (velocity and temperature) at each time step.

We have pre-run this simulation by running the forward case, and stored model output as a checkpoint file on our servers. These fields serve as benchmarks for evaluating our inverse problem's performance. To download the reference benchmark checkpoint file if it doesn't already exist, execute the following command:

In [1]:

active-ipynb

![ ! -f adjoint-demo-checkpoint-state.h5 ] && wget https://data.gadopt.org/demos/adjoint-demo-checkpoint-state.h5

![ ! -f adjoint-demo-checkpoint-state.h5 ] && wget https://data.gadopt.org/demos/adjoint-demo-checkpoint-state.h5

In this file, fields from the reference simulation are stored under the names "Temperature" and "Velocity". After importing g-adopt and the associated inverse module (gadopt.inverse - discussed further below), we can retrieve timestepping information from the pre-computed forward run as follows

In [2]:

from gadopt import *
from gadopt.inverse import *

# Open the checkpoint file and subsequently load the mesh:
checkpoint_filename = "adjoint-demo-checkpoint-state.h5"
checkpoint_file = CheckpointFile(checkpoint_filename, mode="r")
mesh = checkpoint_file.load_mesh("firedrake_default_extruded")
mesh.cartesian = True

# Specify boundary markers, noting that for extruded meshes the upper and lower boundaries are tagged as
# "top" and "bottom" respectively.
boundary = get_boundary_ids(mesh)

# Retrieve the timestepping information for the Velocity and Temperature functions from checkpoint file:
temperature_timestepping_info = checkpoint_file.get_timestepping_history(mesh, "Temperature")
velocity_timestepping_info = checkpoint_file.get_timestepping_history(mesh, "Velocity")

from gadopt import * from gadopt.inverse import * # Open the checkpoint file and subsequently load the mesh: checkpoint_filename = "adjoint-demo-checkpoint-state.h5" checkpoint_file = CheckpointFile(checkpoint_filename, mode="r") mesh = checkpoint_file.load_mesh("firedrake_default_extruded") mesh.cartesian = True # Specify boundary markers, noting that for extruded meshes the upper and lower boundaries are tagged as # "top" and "bottom" respectively. boundary = get_boundary_ids(mesh) # Retrieve the timestepping information for the Velocity and Temperature functions from checkpoint file: temperature_timestepping_info = checkpoint_file.get_timestepping_history(mesh, "Temperature") velocity_timestepping_info = checkpoint_file.get_timestepping_history(mesh, "Velocity")

We can check the information for each:

In [3]:

active-ipynb

print("Timestepping info for Temperature", temperature_timestepping_info)
print("Timestepping info for Velocity", velocity_timestepping_info)

print("Timestepping info for Temperature", temperature_timestepping_info) print("Timestepping info for Velocity", velocity_timestepping_info)

Timestepping info for Temperature {'index': array([ 0., 79.])}
Timestepping info for Velocity {'index': array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,
       13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,
       26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,
       39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51.,
       52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62., 63., 64.,
       65., 66., 67., 68., 69., 70., 71., 72., 73., 74., 75., 76., 77.,
       78., 79.]), 'delta_t': array([4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06,
       4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06, 4.e-06])}

The timestepping information reveals that there are 80 time-steps (from 0 to 79) in the reference simulation, with the temperature field stored only at the initial (index=0) and final (index=79) timesteps, while the velocity field is stored at all timesteps. We can visualise the benchmark fields using Firedrake's built-in VTK functionality. For example, initial and final temperature fields can be loaded:

In [4]:

# Load the final state, analagous to the present-day "observed" state:
Tobs = checkpoint_file.load_function(mesh, "Temperature", idx=int(temperature_timestepping_info["index"][-1]))
Tobs.rename("Observed Temperature")
# Load the reference initial state - i.e. the state that we wish to recover:
Tic_ref = checkpoint_file.load_function(mesh, "Temperature", idx=int(temperature_timestepping_info["index"][0]))
Tic_ref.rename("Reference Initial Temperature")
checkpoint_file.close()

# Load the final state, analagous to the present-day "observed" state: Tobs = checkpoint_file.load_function(mesh, "Temperature", idx=int(temperature_timestepping_info["index"][-1])) Tobs.rename("Observed Temperature") # Load the reference initial state - i.e. the state that we wish to recover: Tic_ref = checkpoint_file.load_function(mesh, "Temperature", idx=int(temperature_timestepping_info["index"][0])) Tic_ref.rename("Reference Initial Temperature") checkpoint_file.close()

These fields can be visualised using standard VTK software, such as Paraview or pyvista.

In [5]:

active-ipynb

import pyvista as pv
VTKFile("./visualisation_vtk.pvd").write(Tobs, Tic_ref)
dataset = pv.read('./visualisation_vtk.pvd')
# Create a plotter object
plotter = pv.Plotter()
# Add the dataset to the plotter
plotter.add_mesh(dataset, scalars='Observed Temperature', cmap='coolwarm')
# Adjust the camera position
plotter.camera_position = [(0.5, 0.5, 2.5), (0.5, 0.5, 0), (0, 1, 0)]
# Show the plot
plotter.show(jupyter_backend="static")

import pyvista as pv VTKFile("./visualisation_vtk.pvd").write(Tobs, Tic_ref) dataset = pv.read('./visualisation_vtk.pvd') # Create a plotter object plotter = pv.Plotter() # Add the dataset to the plotter plotter.add_mesh(dataset, scalars='Observed Temperature', cmap='coolwarm') # Adjust the camera position plotter.camera_position = [(0.5, 0.5, 2.5), (0.5, 0.5, 0), (0, 1, 0)] # Show the plot plotter.show(jupyter_backend="static")

The Inverse Code¶

The novelty of using the overloading approach provided by pyadjoint is that it requires minimal changes to our script to enable the inverse capabalities of G-ADOPT. To turn on the adjoint, one simply imports the inverse module (already done above) to enable all taping functionality from pyadjoint.

Doing so will turn Firedrake's objects to overloaded types, in a way that any UFL operation will be annotated and added to the tape, unless otherwise specified.

We first ensure that the tape is cleared of any previous operations, using the following code:

In [6]:

tape = get_working_tape()
tape.clear_tape()

tape = get_working_tape() tape.clear_tape()

In [7]:

active-ipynb

# To verify the tape is empty, we can print all blocks:
print(tape.get_blocks())

# To verify the tape is empty, we can print all blocks: print(tape.get_blocks())

[]

From here on, all user operations are specified with minimal differences relative to to our forward code. Under the hood, however, the tape will be populated by blocks that record their dependencies. Knowing the mesh was loaded above, we continue in a manner that is consistent with our most basic forward modelling tutorials.

In [8]:

# Set up function spaces:
V = VectorFunctionSpace(mesh, "CG", 2)  # Velocity function space (vector)
W = FunctionSpace(mesh, "CG", 1)  # Pressure function space (scalar)
Q = FunctionSpace(mesh, "DQ", 2)  # Temperature function space (DQ; scalar)
Z = MixedFunctionSpace([V, W])  # Mixed function space

# Specify test functions and functions to hold solutions:
z = Function(Z)  # A field over the mixed function space Z
u, p = split(z)  # Returns symbolic UFL expression for u and p
z.subfunctions[0].rename("Velocity")
z.subfunctions[1].rename("Pressure")
T = Function(Q, name="Temperature")

# Specify important constants for the problem, alongside the approximation:
Ra = Constant(1e6)  # Rayleigh number
approximation = BoussinesqApproximation(Ra)

# Define time-stepping parameters:
delta_t = Constant(4e-6)  # Constant time step
timesteps = int(temperature_timestepping_info["index"][-1]) + 1  # number of timesteps from forward

# Nullspaces for the problem are next defined:
Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=False)

# Followed by boundary conditions, noting that all boundaries are free slip, whilst the domain is
# heated from below (T = 1) and cooled from above (T = 0).
stokes_bcs = {
    boundary.bottom: {"uy": 0},
    boundary.top: {"uy": 0},
    boundary.left: {"ux": 0},
    boundary.right: {"ux": 0},
}
temp_bcs = {
    boundary.bottom: {"T": 1.0},
    boundary.top: {"T": 0.0},
}

# Setup Energy and Stokes solver
energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs)
stokes_solver = StokesSolver(
    z,
    approximation,
    T,
    bcs=stokes_bcs,
    nullspace=Z_nullspace,
    transpose_nullspace=Z_nullspace,
)

# Set up function spaces: V = VectorFunctionSpace(mesh, "CG", 2) # Velocity function space (vector) W = FunctionSpace(mesh, "CG", 1) # Pressure function space (scalar) Q = FunctionSpace(mesh, "DQ", 2) # Temperature function space (DQ; scalar) Z = MixedFunctionSpace([V, W]) # Mixed function space # Specify test functions and functions to hold solutions: z = Function(Z) # A field over the mixed function space Z u, p = split(z) # Returns symbolic UFL expression for u and p z.subfunctions[0].rename("Velocity") z.subfunctions[1].rename("Pressure") T = Function(Q, name="Temperature") # Specify important constants for the problem, alongside the approximation: Ra = Constant(1e6) # Rayleigh number approximation = BoussinesqApproximation(Ra) # Define time-stepping parameters: delta_t = Constant(4e-6) # Constant time step timesteps = int(temperature_timestepping_info["index"][-1]) + 1 # number of timesteps from forward # Nullspaces for the problem are next defined: Z_nullspace = create_stokes_nullspace(Z, closed=True, rotational=False) # Followed by boundary conditions, noting that all boundaries are free slip, whilst the domain is # heated from below (T = 1) and cooled from above (T = 0). stokes_bcs = { boundary.bottom: {"uy": 0}, boundary.top: {"uy": 0}, boundary.left: {"ux": 0}, boundary.right: {"ux": 0}, } temp_bcs = { boundary.bottom: {"T": 1.0}, boundary.top: {"T": 0.0}, } # Setup Energy and Stokes solver energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint, bcs=temp_bcs) stokes_solver = StokesSolver( z, approximation, T, bcs=stokes_bcs, nullspace=Z_nullspace, transpose_nullspace=Z_nullspace, )

Specify Problem Length¶

For the purpose of this demo, we only invert for a total of 10 time-steps. This makes it tractable to run this within a tutorial session.

To run for the simulation's full duration, change the initial_timestep to 0 below, rather than timesteps - 10.

In [9]:

initial_timestep = timesteps - 10

initial_timestep = timesteps - 10

Define the Control Space¶

In this section, we define the control space, which can be restricted to reduce the risk of encountering an undetermined problem. Here, we select the Q1 function space for the initial condition $T_{ic}$. We also provide an initial guess for the control value, which in this synthetic test is the temperature field of the reference simulation at the final time-step (timesteps - 1). In other words, our guess for the initial temperature is the final model state.

In [10]:

# Define control function space:
Q1 = FunctionSpace(mesh, "CG", 1)

# Create a function for the unknown initial temperature condition, which we will be inverting for. Our initial
# guess is set to the 1-D average of the forward model. We first load that, at the relevant timestep.
# Note that this layer average will later be used for the smoothing term in our objective functional.
with CheckpointFile(checkpoint_filename, mode="r") as checkpoint_file:
    Taverage = checkpoint_file.load_function(mesh, "Average_Temperature", idx=initial_timestep)
Tic = Function(Q1, name="Initial_Condition_Temperature").assign(Taverage)

# Given that Tic will be updated during the optimisation, we also create a function to store our initial guess,
# which we will later use for smoothing. Note that since smoothing is executed in the control space, we must
# specify boundary conditions on this term in that same Q1 space.
T0_bcs = [DirichletBC(Q1, 0., boundary.top), DirichletBC(Q1, 1., boundary.bottom)]
T0 = Function(Q1, name="Initial_Guess_Temperature").project(Tic, bcs=T0_bcs)

# We next make pyadjoint aware of our control problem:
control = Control(Tic)

# Take our initial guess and project to T, simultaneously applying boundary conditions in the Q2 space:
T.project(Tic, bcs=energy_solver.strong_bcs)

# We continue by integrating the solutions at each time-step.
# Notice that we cumulatively compute the misfit term with respect to the
# surface velocity observable.

# Define control function space: Q1 = FunctionSpace(mesh, "CG", 1) # Create a function for the unknown initial temperature condition, which we will be inverting for. Our initial # guess is set to the 1-D average of the forward model. We first load that, at the relevant timestep. # Note that this layer average will later be used for the smoothing term in our objective functional. with CheckpointFile(checkpoint_filename, mode="r") as checkpoint_file: Taverage = checkpoint_file.load_function(mesh, "Average_Temperature", idx=initial_timestep) Tic = Function(Q1, name="Initial_Condition_Temperature").assign(Taverage) # Given that Tic will be updated during the optimisation, we also create a function to store our initial guess, # which we will later use for smoothing. Note that since smoothing is executed in the control space, we must # specify boundary conditions on this term in that same Q1 space. T0_bcs = [DirichletBC(Q1, 0., boundary.top), DirichletBC(Q1, 1., boundary.bottom)] T0 = Function(Q1, name="Initial_Guess_Temperature").project(Tic, bcs=T0_bcs) # We next make pyadjoint aware of our control problem: control = Control(Tic) # Take our initial guess and project to T, simultaneously applying boundary conditions in the Q2 space: T.project(Tic, bcs=energy_solver.strong_bcs) # We continue by integrating the solutions at each time-step. # Notice that we cumulatively compute the misfit term with respect to the # surface velocity observable.

Out[10]:

Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.ExtrudedMeshTopology object at 0x7bd1ace25af0>, TensorProductElement(FiniteElement('Discontinuous Lagrange', interval, 2), FiniteElement('Discontinuous Lagrange', interval, 2), cell=TensorProductCell(interval, interval)), name=None), Mesh(VectorElement(TensorProductElement(FiniteElement('Lagrange', interval, 1), FiniteElement('Lagrange', interval, 1), cell=TensorProductCell(interval, interval)), dim=2), 1)), 31)

In [11]:

u_misfit = 0.0

# Next populate the tape by running the forward simulation.
for time_idx in range(initial_timestep, timesteps):
    stokes_solver.solve()
    energy_solver.solve()
    # Update the accumulated surface velocity misfit using the observed value.
    with CheckpointFile(checkpoint_filename, mode="r") as checkpoint_file:
        uobs = checkpoint_file.load_function(mesh, name="Velocity", idx=time_idx)
    u_misfit += assemble(dot(u - uobs, u - uobs) * ds_t)

u_misfit = 0.0 # Next populate the tape by running the forward simulation. for time_idx in range(initial_timestep, timesteps): stokes_solver.solve() energy_solver.solve() # Update the accumulated surface velocity misfit using the observed value. with CheckpointFile(checkpoint_filename, mode="r") as checkpoint_file: uobs = checkpoint_file.load_function(mesh, name="Velocity", idx=time_idx) u_misfit += assemble(dot(u - uobs, u - uobs) * ds_t)

Define the Objective Functional¶

Now that all calculations are in place, we must define the objective functional. The objective functional is our way of expressing our goal for this optimisation. It is composed of several terms, each representing a different aspect of the model's performance and regularisation.

Regularisation involves imposing constraints on solutions to prevent overfitting, ensuring that the model generalises well to new data. In this context, we use the one-dimensional (1-D) temperature profile derived from the reference simulation as our regularisation constraint. This profile, referred to below as Taverage, helps stabilise the inversion process by providing a benchmark that guides the solution towards physically plausible states.

We use Taverage as a part of the damping and smoothing terms in our regularisation. Consequently, the complete objective functional is defined mathematically as follows:

Reiterating that:

• $T_{ic}$ is the initial temperature condition.
• $T_{\text{average}}$ is the average temperature profile representing mantle's geotherm.
• $T_{F}$ is the the temperature field at final time-step.
• $T_{\text{obs}}$ is the observed temperature field at the final time-step.
• $u_{\text{obs}}$ is the observed velocity field at each time-step.
• $\alpha_u$, $\alpha_d$, $\alpha_s$ are the three different weighting terms for the velocity, damping and smoothing terms.

We define the objective functional as $$ \text{Objective Functional}= \int_{\Omega}(T - T_{\text{obs}}) ^ 2 \, dx \\ +\alpha_u\, \frac{D_{T_{obs}}}{N\times D_{u_{obs}}}\sum_{i}\int_{\partial \Omega_{\text{top}}}(u - u_{\text{obs}}) \cdot(u - u_{\text{obs}}) \, ds \\ +\alpha_s\, \frac{D_{T_{obs}}}{D_{\text{smoothing}}}\int_{\Omega} \nabla(T_{ic} - T_{\text{average}}) \cdot \nabla(T_{ic} - T_{\text{average}}) \, dx \\ +\alpha_d\, \frac{D_{T_{obs}}}{D_{\text{damping}}}\int_{\Omega}(T_{ic} - T_{\text{average}}) ^ 2 \, dx $$

With the three normlisation terms of:

• $D_{\text{damping}} = \int_{\Omega} T_{\text{average}}^2 \, dx$,
• $D_{\text{smoothing}} = \int_{\Omega} \nabla T_{\text{obs}} \cdot \nabla T_{\text{obs}} \, dx$,
• $D_{T_{obs}} = \int_{\Omega} T_{\text{obs}} ^ 2 \, dx$, and
• $D_{\text{damping}} = \int_{\partial \Omega_{\text{top}}} u_{\text{obs}} \cdot u_{\text{obs}} \, ds$

which we specify through the objective below:

In [12]:

# Define component terms of overall objective functional and their normalisation terms:
damping = assemble((T0 - Taverage) ** 2 * dx)
norm_damping = assemble(Taverage**2 * dx)
smoothing = assemble(dot(grad(T0 - Taverage), grad(T0 - Taverage)) * dx)
norm_smoothing = assemble(dot(grad(Tobs), grad(Tobs)) * dx)
norm_obs = assemble(Tobs**2 * dx)
norm_u_surface = assemble(dot(uobs, uobs) * ds_t)

# Define temperature misfit between final state solution and observation:
t_misfit = assemble((T - Tobs) ** 2 * dx)

# Weighting terms
alpha_u = 1e-1
alpha_d = 1e-3
alpha_s = 1e-3

# Define overall objective functional:
objective = (
    t_misfit +
    alpha_u * (norm_obs * u_misfit / timesteps / norm_u_surface) +
    alpha_d * (norm_obs * damping / norm_damping) +
    alpha_s * (norm_obs * smoothing / norm_smoothing)
)

# Define component terms of overall objective functional and their normalisation terms: damping = assemble((T0 - Taverage) ** 2 * dx) norm_damping = assemble(Taverage**2 * dx) smoothing = assemble(dot(grad(T0 - Taverage), grad(T0 - Taverage)) * dx) norm_smoothing = assemble(dot(grad(Tobs), grad(Tobs)) * dx) norm_obs = assemble(Tobs**2 * dx) norm_u_surface = assemble(dot(uobs, uobs) * ds_t) # Define temperature misfit between final state solution and observation: t_misfit = assemble((T - Tobs) ** 2 * dx) # Weighting terms alpha_u = 1e-1 alpha_d = 1e-3 alpha_s = 1e-3 # Define overall objective functional: objective = ( t_misfit + alpha_u * (norm_obs * u_misfit / timesteps / norm_u_surface) + alpha_d * (norm_obs * damping / norm_damping) + alpha_s * (norm_obs * smoothing / norm_smoothing) )

Define the Reduced Functional¶

In optimisation terminology, a reduced functional is a functional that takes a given value for the control and outputs the value of the objective functional defined for it. It does this without explicitly depending on all intermediary state variables, hence the name "reduced".

To define the reduced functional, we provide the class with an objective (which is an overloaded UFL object) and the control.

In [13]:

reduced_functional = ReducedFunctional(objective, control)

reduced_functional = ReducedFunctional(objective, control)

At this point, we have completed annotating the tape with the necessary information from running the forward simulation. To prevent further annotations during subsequent operations, we stop the annotation process. This ensures that no additional solves are unnecessarily recorded, keeping the tape focused only on the essential steps.

In [14]:

pause_annotation()

pause_annotation()

We can print the contents of the tape at this stage to verify that it is not empty.

In [15]:

active-ipynb

print(tape.get_blocks())

print(tape.get_blocks())

[<firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a31bdcd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a31bcef0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a31be030>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a31bf6b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a31bfe90>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1941c8d70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1941c8e90>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a31be210>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd1940edd90>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1940eef90>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd1940eddf0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd1940ee750>, <firedrake.adjoint_utils.blocks.solving.ProjectBlock object at 0x7bd1940ee030>, <firedrake.adjoint_utils.blocks.solving.ProjectBlock object at 0x7bd1940ed4f0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd1a6a61f10>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd194267c50>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd1a6b5a870>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd1a6b5b050>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd1a6a62030>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1a6cd4890>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422ca10>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422cef0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422f170>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422d9d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422f7d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422eff0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd19422ed50>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd19422e630>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1940bbcb0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1940b9f70>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd1940bb4d0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c50e330>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd1940b9d30>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61ee10>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61e090>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61e5d0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61dcd0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd1940bb830>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1940bb710>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61d5b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61d910>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61f410>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61d730>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61d4f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61e570>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61dbb0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c61ff50>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61cc50>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61ce90>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c216f30>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c215af0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd13c2170b0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c216e70>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c214230>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c216450>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c217590>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c214d70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c217230>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c216ed0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c2157f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c214b30>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c214950>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c217cb0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c214650>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c215790>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c2159d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd1ac5aabd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cddf0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c61d3d0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c61eed0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd13c61fd70>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61f350>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61fef0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61cb30>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c61d670>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c61cd10>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61f8f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61fd10>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61db50>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61d0d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61ce30>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61e8d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61fa70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c61d970>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c61e6f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd18c355fd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd18c355d30>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd18c354410>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c7cf170>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd18c356030>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7cc170>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7cef90>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7cc470>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7ceb10>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd18c3574d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd18c3578f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7ce9f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cccb0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cf590>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cd8b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cf350>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cdcd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7ce7b0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c7cc530>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7cd6d0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e0b30>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c7e2f30>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c7e2630>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd13c7e22d0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7e0830>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7e3e30>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7e0ef0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd13c7e0170>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c7e1970>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e0bf0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e3890>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e03b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e2bd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e3a10>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e0fb0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e3530>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd13c7e0650>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd13c7e0c50>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f83102f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8311a90>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8311790>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8310650>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f83116d0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f8310710>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f83111f0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f8310b90>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f8310b30>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f8311af0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8311130>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f83115b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8310dd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8310c50>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8311070>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8310d10>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8310d70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f83110d0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f8310e30>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8313a70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8313d10>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f85ec1d0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f85ec770>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8313e30>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85ec710>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85ec950>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85ec830>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85ec8f0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f8313e90>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ec9b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ec050>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ecd70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ece30>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ecad0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ecef0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ecdd0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85eca70>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f85eccb0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ee270>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ee510>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f85ee930>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f85eef30>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f85ee990>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85eeed0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85ef110>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85eeff0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85ef0b0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f85ee570>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef170>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85eeb10>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef530>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef5f0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef290>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef6b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef590>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85ef230>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f85ef470>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8510a70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8510d10>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f85110d0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8511730>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8511130>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85116d0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f8511910>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85117f0>, <firedrake.adjoint_utils.blocks.dirichlet_bc.DirichletBCBlock object at 0x7bd0f85118b0>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f8510d70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511970>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85112b0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511d30>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511df0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511a90>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511eb0>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511d90>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8511a30>, <firedrake.adjoint_utils.blocks.solving.NonlinearVariationalSolveBlock object at 0x7bd0f8511c70>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f8513230>, <firedrake.adjoint_utils.blocks.function.FunctionAssignBlock object at 0x7bd0f85134d0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8513890>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8513f50>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f85139b0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd11811a090>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd118118d70>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f82634d0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd13c61f9b0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd1187755b0>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8512f90>, <firedrake.adjoint_utils.blocks.assembly.AssembleBlock object at 0x7bd0f8510830>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8512270>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8512150>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8510770>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f85108f0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f85106b0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8512ab0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f85123f0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8510890>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8512210>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8511e50>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8510650>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f85109b0>, <pyadjoint.adjfloat.AdjFloatExprBlock object at 0x7bd0f8510a10>]

Verification of Gradients: Taylor Remainder Convergence Test¶

A fundamental tool for verifying gradients is the Taylor remainder convergence test. This test helps ensure that the gradients computed by our optimisation algorithm are accurate. For the reduced functional, $J(T_{ic})$, and its derivative, $\frac{\mathrm{d} J}{\mathrm{d} T_{ic}}$, the Taylor remainder convergence test can be expressed as:

$$ \left| J(T_{ic} + h \,\delta T_{ic}) - J(T_{ic}) - h\,\frac{\mathrm{d} J}{\mathrm{d} T_{ic}} \cdot \delta T_{ic} \right| \longrightarrow 0 \text{ at } O(h^2). $$

The expression on the left-hand side is termed the second-order Taylor remainder. This term's convergence rate of $O(h^2)$ is a robust indicator for verifying the computational implementation of the gradient calculation. Essentially, if you halve the value of $h$, the magnitude of the second-order Taylor remainder should decrease by a factor of 4.

We employ these so-called Taylor tests to confirm the accuracy of the determined gradients. The theoretical convergence rate is $O(2.0)$, and achieving this rate indicates that the gradient information is accurate down to floating-point precision.

Performing Taylor Tests¶

In our implementation, we perform a second-order Taylor remainder test for each term of the objective functional. The test involves computing the functional and the associated gradient when randomly perturbing the initial temperature field, $T_{ic}$, and subsequently halving the perturbations at each level.

Here is how you can perform a Taylor test in the code:

In [16]:

active-ipynb

# Define the perturbation in the initial temperature field
import numpy as np
Delta_temp = Function(Tic.function_space(), name="Delta_Temperature")
Delta_temp.dat.data[:] = np.random.random(Delta_temp.dat.data.shape)

# Perform the Taylor test to verify the gradients
minconv = taylor_test(reduced_functional, Tic, Delta_temp)

# Define the perturbation in the initial temperature field import numpy as np Delta_temp = Function(Tic.function_space(), name="Delta_Temperature") Delta_temp.dat.data[:] = np.random.random(Delta_temp.dat.data.shape) # Perform the Taylor test to verify the gradients minconv = taylor_test(reduced_functional, Tic, Delta_temp)

Running Taylor test

Computed residuals: [4.981239984062177e-05, 1.2453100034136091e-05, 3.113275018399343e-06, 7.783187561480504e-07]
Computed convergence rates: [1.9999999914293214, 1.9999999954283996, 1.9999999971302225]

The taylor_test function computes the Taylor remainder and verifies that the convergence rate is close to the theoretical value of $O(2.0)$. This ensures that our gradients are accurate and reliable for optimisation.

Running the inversion¶

In the final section of this tutorial, we run the optimisation method. First, we define lower and upper bounds for the optimisation problem to guide the optimisation method towards a more constrained solution.

For this simple problem, we perform a bounded nonlinear optimisation where the temperature is only permitted to lie in the range [0, 1]. This means that the optimisation problem should not search for solutions beyond these values.

In [17]:

# Define lower and upper bounds for the temperature
T_lb = Function(Tic.function_space(), name="Lower Bound Temperature")
T_ub = Function(Tic.function_space(), name="Upper Bound Temperature")

# Assign the bounds
T_lb.assign(0.0)
T_ub.assign(1.0)

# Define the minimisation problem, with the goal to minimise the reduced functional
# Note: in some scenarios, the goal might be to maximise (rather than minimise) the functional.
minimisation_problem = MinimizationProblem(reduced_functional, bounds=(T_lb, T_ub))

# Define lower and upper bounds for the temperature T_lb = Function(Tic.function_space(), name="Lower Bound Temperature") T_ub = Function(Tic.function_space(), name="Upper Bound Temperature") # Assign the bounds T_lb.assign(0.0) T_ub.assign(1.0) # Define the minimisation problem, with the goal to minimise the reduced functional # Note: in some scenarios, the goal might be to maximise (rather than minimise) the functional. minimisation_problem = MinimizationProblem(reduced_functional, bounds=(T_lb, T_ub))

Using the Lin-Moré optimiser¶

In this tutorial, we employ the trust region method of Lin and Moré (1999) implemented in ROL (Rapid Optimization Library). Lin-Moré is a truncated Newton method, which involves the repeated application of an iterative algorithm to approximately solve Newton’s equations (Dembo and Steihaug, 1983).

Lin-Moré effectively handles provided bound constraints by ensuring that variables remain within their specified bounds. During each iteration, variables are classified into "active" and "inactive" sets. Variables at their bounds that do not allow descent are considered active and are fixed during the iteration. The remaining variables, which can change without violating the bounds, are inactive. These properties make the algorithm robust and efficient for solving bound-constrained optimisation problems.

For our solution of the optimisation problem we use the pre-defined paramters set in gadopt by using minimsation_parameters. Here, we set the number of iterations to only 5, as opposed to the default 100. We also adjust the step-length for this problem, by setting it to a lower value than our default.

In [18]:

minimisation_parameters["Status Test"]["Iteration Limit"] = 5
minimisation_parameters["Step"]["Trust Region"]["Initial Radius"] = 1e-2

minimisation_parameters["Status Test"]["Iteration Limit"] = 5 minimisation_parameters["Step"]["Trust Region"]["Initial Radius"] = 1e-2

A notable feature of this optimisation approach in ROL is its checkpointing capability. For every iteration, all information necessary to restart the optimisation from that iteration is saved in the specified checkpoint_dir.

In [19]:

# Define the LinMore Optimiser class with checkpointing capability:
optimiser = LinMoreOptimiser(
    minimisation_problem,
    minimisation_parameters,
    checkpoint_dir="optimisation_checkpoint",
)

# Define the LinMore Optimiser class with checkpointing capability: optimiser = LinMoreOptimiser( minimisation_problem, minimisation_parameters, checkpoint_dir="optimisation_checkpoint", )

For sake of book-keeping the simulation, we have also implemented a user-defined way of recording information that might be used to check the optimisation performance. This callback function will be executed at the end of each iteration. Here, we write out the control field, i.e., the reconstructed intial temperature field, at the end of each iteration. To access the last value of an overloaded object we should access the .block_variable.checkpoint method as below.

For the sake of this demo, we also record the values of the reduced functional directly in order to produce a plot of the convergence.

In [20]:

solutions_vtk = VTKFile("solutions.pvd")
solution_container = Function(Tic.function_space(), name="Solutions")
functional_values = []


def callback():
    solution_container.assign(Tic.block_variable.checkpoint)
    solutions_vtk.write(solution_container)
    final_temperature_misfit = assemble(
        (T.block_variable.checkpoint - Tobs) ** 2 * dx
    )
    log(f"Terminal Temperature Misfit: {final_temperature_misfit}")


def record_value(value, *args):
    functional_values.append(value)


optimiser.add_callback(callback)
reduced_functional.eval_cb_post = record_value

# If it existed, we could restore the optimisation from last checkpoint:
# optimiser.restore()

# Run the optimisation
optimiser.run()

# Write the functional values to a file
with open("functional.txt", "w") as f:
    f.write("\n".join(str(x) for x in functional_values))

solutions_vtk = VTKFile("solutions.pvd") solution_container = Function(Tic.function_space(), name="Solutions") functional_values = [] def callback(): solution_container.assign(Tic.block_variable.checkpoint) solutions_vtk.write(solution_container) final_temperature_misfit = assemble( (T.block_variable.checkpoint - Tobs) ** 2 * dx ) log(f"Terminal Temperature Misfit: {final_temperature_misfit}") def record_value(value, *args): functional_values.append(value) optimiser.add_callback(callback) reduced_functional.eval_cb_post = record_value # If it existed, we could restore the optimisation from last checkpoint: # optimiser.restore() # Run the optimisation optimiser.run() # Write the functional values to a file with open("functional.txt", "w") as f: f.write("\n".join(str(x) for x in functional_values))

Lin-More Trust-Region Method (Type B, Bound Constraints)
  iter  value          gnorm          snorm          delta          #fval     #grad     #hess     #proj     tr_flag   iterCG    flagCG    
  0     4.357736e-02   4.000463e-01   ---            1.000000e-02   1         1         0         2         ---       ---       ---       

Terminal Temperature Misfit: 0.04066260840368385

  1     3.781235e-02   3.809691e-01   1.000000e-02   1.000000e-01   2         2         3         7         0         1         3         

Terminal Temperature Misfit: 0.03532784476773855

  2     3.484459e-03   1.247641e-01   1.000000e-01   1.000000e+00   3         3         25        22        0         10        1         

Terminal Temperature Misfit: 0.003322779307480237

  3     1.056904e-03   1.759806e-01   3.867398e-02   1.000000e+00   4         4         30        27        0         2         0         

Terminal Temperature Misfit: 0.0005033178057787317

  4     1.056904e-03   1.759806e-01   2.326458e-02   5.816145e-03   5         4         42        31        2         8         0         

Terminal Temperature Misfit: 0.0005033178057787317

  5     5.864116e-04   7.684261e-02   5.816145e-03   5.816145e-03   6         5         63        45        0         10        1         

Terminal Temperature Misfit: 0.0004449546419173902
Optimization Terminated with Status: Last Type (Dummy)

At this point a total number of 5 iterations are performed. For the example case here with 10 timesteps this should result an adequete reduction in the objective functional. Now we can look at the solution visually. For the actual simulation with 80 time-steps, this solution could be compared to Tic_ref as the "true solution".

In [21]:

active-ipynb

import pyvista as pv
VTKFile("./solution.pvd").write(optimiser.rol_solver.rolvector.dat[0])
dataset = pv.read('./solution.pvd')
# Create a plotter object
plotter = pv.Plotter()
# Add the dataset to the plotter
plotter.add_mesh(dataset, scalars=dataset[0].array_names[0], cmap='coolwarm')
plotter.add_text("Solution after 5 iterations", font_size=10)
# Adjust the camera position
plotter.camera_position = [(0.5, 0.5, 2.5), (0.5, 0.5, 0), (0, 1, 0)]
# Show the plot
plotter.show(jupyter_backend="static")

import pyvista as pv VTKFile("./solution.pvd").write(optimiser.rol_solver.rolvector.dat[0]) dataset = pv.read('./solution.pvd') # Create a plotter object plotter = pv.Plotter() # Add the dataset to the plotter plotter.add_mesh(dataset, scalars=dataset[0].array_names[0], cmap='coolwarm') plotter.add_text("Solution after 5 iterations", font_size=10) # Adjust the camera position plotter.camera_position = [(0.5, 0.5, 2.5), (0.5, 0.5, 0), (0, 1, 0)] # Show the plot plotter.show(jupyter_backend="static")

In [22]:

active-ipynb

import matplotlib.pyplot as plt
plt.plot(functional_values)
plt.xlabel("Optimisation iteration")
plt.ylabel("Reduced functional")
plt.title("Optimisation convergence")

import matplotlib.pyplot as plt plt.plot(functional_values) plt.xlabel("Optimisation iteration") plt.ylabel("Reduced functional") plt.title("Optimisation convergence")

Out[22]:

Text(0.5, 1.0, 'Optimisation convergence')