PDE Constrained Optimisation with G-ADOPT - Boundary Values¶

In this tutorial, we undertake an inversion for an (unknown) initial condition, to match given time-dependent boundary values. This differs to our previous tutorial, where our goal was to match a given final state.

We start with our usual imports:

In [1]:

from gadopt import *
from gadopt.inverse import *

from gadopt import * from gadopt.inverse import *

Create synthetic twin experiment and record solution at all timesteps¶

Note that the setup is similar to our previous example, except that the velocity is now counter clockwise around the origin $(0,0)$ in the corner of the unit square domain. This implies an inflow at the bottom boundary and an outflow boundary on the left of the domain.

In [2]:

mesh = UnitSquareMesh(40, 40)
mesh.cartesian = True
boundary = get_boundary_ids(mesh)

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
T = Function(Q, name='Temperature')
T0 = Function(Q, name="Initial_Temperature")  # T Initial condition which we will invert for.
T0_ref = Function(Q, name="Reference_Initial_Temperature")

x, y = SpatialCoordinate(mesh)
u = Function(V, name="Velocity").interpolate(as_vector((-y, x)))

approximation = BoussinesqApproximation(Ra=1, kappa=5e-2)

delta_t = 0.1
energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint)

mesh = UnitSquareMesh(40, 40) mesh.cartesian = True boundary = get_boundary_ids(mesh) V = VectorFunctionSpace(mesh, "CG", 2) Q = FunctionSpace(mesh, "CG", 1) T = Function(Q, name='Temperature') T0 = Function(Q, name="Initial_Temperature") # T Initial condition which we will invert for. T0_ref = Function(Q, name="Reference_Initial_Temperature") x, y = SpatialCoordinate(mesh) u = Function(V, name="Velocity").interpolate(as_vector((-y, x))) approximation = BoussinesqApproximation(Ra=1, kappa=5e-2) delta_t = 0.1 energy_solver = EnergySolver(T, u, approximation, delta_t, ImplicitMidpoint)

The initial condition that we, again, later will invert for, is now centered in the domain.

In [3]:

x0, y0 = 0.5, 0.5
w = .2
r2 = (x-x0)**2 + (y-y0)**2
T0_ref.interpolate(exp(-r2/w**2))

x0, y0 = 0.5, 0.5 w = .2 r2 = (x-x0)**2 + (y-y0)**2 T0_ref.interpolate(exp(-r2/w**2))

Out[3]:

Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x72be12546ab0>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0)), 6)

In [4]:

active-ipynb

import matplotlib.pyplot as plt
fig, axes = plt.subplots()
collection = tripcolor(T0_ref, axes=axes, cmap='magma', vmax=0.5)
fig.colorbar(collection);

import matplotlib.pyplot as plt fig, axes = plt.subplots() collection = tripcolor(T0_ref, axes=axes, cmap='magma', vmax=0.5) fig.colorbar(collection);

After setting the initial condition for T, we run this simulation for 20 timesteps to ensure the entire Gaussian has left the domain. For this example, we checkpoint the solution at every timestep, so that we can later use it as the target boundary values.

In [5]:

num_timesteps = 20
T.project(T0_ref)
with CheckpointFile("Model_State.h5", "w") as model_checkpoint:
    model_checkpoint.save_mesh(mesh)
    for timestep in range(num_timesteps):
        model_checkpoint.save_function(T, idx=timestep)
        energy_solver.solve()
    # After saving idx=0, 19 at beginning of each timestep, we include idx=20 for the solution at
    # the end of the final timestep:
    model_checkpoint.save_function(T, idx=timestep)

num_timesteps = 20 T.project(T0_ref) with CheckpointFile("Model_State.h5", "w") as model_checkpoint: model_checkpoint.save_mesh(mesh) for timestep in range(num_timesteps): model_checkpoint.save_function(T, idx=timestep) energy_solver.solve() # After saving idx=0, 19 at beginning of each timestep, we include idx=20 for the solution at # the end of the final timestep: model_checkpoint.save_function(T, idx=timestep)

The solution has almost completely disappeared (note the different scalebar):

In [6]:

active-ipynb

fig, axes = plt.subplots()
collection = tripcolor(T, axes=axes, cmap='magma', vmax=0.1)
fig.colorbar(collection);

fig, axes = plt.subplots() collection = tripcolor(T, axes=axes, cmap='magma', vmax=0.1) fig.colorbar(collection);

Advection diffusion model with unknown initial condition¶

As with our previous example, we again set up the model with the same configuration, albeit where we do not know the initial condition. We will try to find the optimal initial condition such that we closely match the recorded outflow boundary values.

In [7]:

with CheckpointFile("Model_State.h5", "r") as model_checkpoint:
    mesh = model_checkpoint.load_mesh()
    mesh.cartesian = True

with CheckpointFile("Model_State.h5", "r") as model_checkpoint: mesh = model_checkpoint.load_mesh() mesh.cartesian = True

We now set up the model exactly as before:

In [8]:

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
T = Function(Q, name='Temperature')
T0 = Function(Q, name="Initial_Temperature")
T0_ref = Function(Q, name="Reference_Initial_Temperature")
T_wrong = Function(Q, name="Wrong_Initial_Temperature")

x, y = SpatialCoordinate(mesh)
u = Function(V, name="Velocity").interpolate(as_vector((-y, x)))

approximation = BoussinesqApproximation(Ra=1, kappa=5e-2)
delta_t = 0.1
energy_solver = EnergySolver(
    T, u, approximation, delta_t, ImplicitMidpoint, solver_parameters_extra={"ksp_converged_reason": DeleteParam}
)

V = VectorFunctionSpace(mesh, "CG", 2) Q = FunctionSpace(mesh, "CG", 1) T = Function(Q, name='Temperature') T0 = Function(Q, name="Initial_Temperature") T0_ref = Function(Q, name="Reference_Initial_Temperature") T_wrong = Function(Q, name="Wrong_Initial_Temperature") x, y = SpatialCoordinate(mesh) u = Function(V, name="Velocity").interpolate(as_vector((-y, x))) approximation = BoussinesqApproximation(Ra=1, kappa=5e-2) delta_t = 0.1 energy_solver = EnergySolver( T, u, approximation, delta_t, ImplicitMidpoint, solver_parameters_extra={"ksp_converged_reason": DeleteParam} )

As a first guess we use a Gaussian that is in the wrong place: centred around $(0.7, 0.7)$ instead of $(0.5, 0.5)$:

In [9]:

x0, y0 = 0.7, 0.7
w = .2
r2 = (x-x0)**2 + (y-y0)**2
T_wrong.interpolate(exp(-r2/w**2))

x0, y0 = 0.7, 0.7 w = .2 r2 = (x-x0)**2 + (y-y0)**2 T_wrong.interpolate(exp(-r2/w**2))

Out[9]:

Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x72be0f312b70>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 90)), 197)

As in our first example, we make sure to clear the tape before our actual model starts and specify the control at the right stage. During the model we load back in the solutions from the synthetic twin, but only use its values at the boundary to compute a mismatch with the current model as an integral over the left boundary. Note that we start calculating the functional already in the first timestep, and we keep on adding terms to it, all of which will still be automatically recorded by the pyadjoint tape.

In [10]:

tape = get_working_tape()
tape.clear_tape()

T0.project(T_wrong)

tape = get_working_tape() tape.clear_tape() T0.project(T_wrong)

Out[10]:

Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x72be0f312b70>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 90)), 193)

At this stage we need to define the control variable for the adjoint-based optimisation problem by wrapping T0 in a Control object. This registers T0 as the optimisation variable with pyadjoint's automatic differentiation framework. Furthermore, the L2 Riesz map specification is crucial for two reasons:

1. It defines the inner product structure on the control space, ensuring that gradients are computed in the $L^2$ Hilbert space with respect to the inner product.
2. It guarantees mesh-independent gradient representations, which is essential for consistent optimisation convergence across different mesh refinements in finite element discretisations.

Note that L2 is the default Riesz map if none is explicitly specified.

In [11]:

m = Control(T0, riesz_map="L2")

J = AdjFloat(0.0)  # Initialise functional
factor = AdjFloat(0.5)  # First & final boundary integral weighted by 0.5 to implement mid-point rule time-integration.

T.project(T0)
with CheckpointFile("Model_State.h5", "r") as model_checkpoint:
    for timestep in range(num_timesteps):
        T_target = model_checkpoint.load_function(mesh, 'Temperature', idx=timestep)
        J = J + factor * assemble((T-T_target)**2*ds(boundary.left))
        factor = 1.0  # Remaining timesteps weighted by 1
        energy_solver.solve()

    T_target = model_checkpoint.load_function(mesh, 'Temperature', idx=timestep)
    # Add final contribution weighted again by 0.5
    J = J + factor * assemble((T-T_target)**2*ds(boundary.left))

print(J)

m = Control(T0, riesz_map="L2") J = AdjFloat(0.0) # Initialise functional factor = AdjFloat(0.5) # First & final boundary integral weighted by 0.5 to implement mid-point rule time-integration. T.project(T0) with CheckpointFile("Model_State.h5", "r") as model_checkpoint: for timestep in range(num_timesteps): T_target = model_checkpoint.load_function(mesh, 'Temperature', idx=timestep) J = J + factor * assemble((T-T_target)**2*ds(boundary.left)) factor = 1.0 # Remaining timesteps weighted by 1 energy_solver.solve() T_target = model_checkpoint.load_function(mesh, 'Temperature', idx=timestep) # Add final contribution weighted again by 0.5 J = J + factor * assemble((T-T_target)**2*ds(boundary.left)) print(J)

0.07749995459676606

We define the reduced functional using the final value of J and the specified control. This allows us to rerun the model with an arbitrary initial condition. As with our previous example, we first try to simply rerun the model with the same "wrong" initial condition, and print the functional.

In [12]:

reduced_functional = ReducedFunctional(J, m)
print(reduced_functional(T_wrong))

reduced_functional = ReducedFunctional(J, m) print(reduced_functional(T_wrong))

0.07749995459676606

Now we re run the model with the "correct" initial condition from the twin experiment, ending up with a near-zero misfit.

In [13]:

x0, y0 = 0.5, 0.5
w = .2
r2 = (x-x0)**2 + (y-y0)**2
T0_ref.interpolate(exp(-r2/w**2))

x0, y0 = 0.5, 0.5 w = .2 r2 = (x-x0)**2 + (y-y0)**2 T0_ref.interpolate(exp(-r2/w**2))

Out[13]:

Coefficient(WithGeometry(FunctionSpace(<firedrake.mesh.MeshTopology object at 0x72be0f312b70>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 90)), 195)

In [14]:

print(reduced_functional(T0_ref))

print(reduced_functional(T0_ref))

6.645421118823317e-06

We can again look at the gradient. We evaluate the gradient around an initial guess of T=0 as the initial condition, noting that when a Function is created its associated data values are zero.

In [15]:

T_wrong.assign(0.0)
reduced_functional(T_wrong)

T_wrong.assign(0.0) reduced_functional(T_wrong)

Out[15]:

0.2512941576647654

Here we calculate the gradient by calling the derivative method of the reduced functional with the L2 Riesz map applied to ensure a grid-independent result.

In [16]:

gradJ = reduced_functional.derivative(apply_riesz=True)

gradJ = reduced_functional.derivative(apply_riesz=True)

In [17]:

active-ipynb

fig, axes = plt.subplots()
collection = tripcolor(gradJ, axes=axes, cmap='viridis', vmin=-5, vmax=5)
fig.colorbar(collection);

fig, axes = plt.subplots() collection = tripcolor(gradJ, axes=axes, cmap='viridis', vmin=-5, vmax=5) fig.colorbar(collection);

Invert for optimal initial condition using gradient-based optimisation algorithm¶

As in the previous example, we can now use ROL to invert for the inital condition. We last evaluated the reduced functional with a zero initial condition as the control value, so this will be our initial guess.

We first set lower and upper bound values for the control, which we can provide as functions in the same function space as the control:

In [18]:

T_lb = Function(Q).assign(0.0)
T_ub = Function(Q).assign(1.0)

T_lb = Function(Q).assign(0.0) T_ub = Function(Q).assign(1.0)

We next specify our minimisation problem using the LinMore algorithm. As this case is a little more challenging than our previous tutorial, we specify 50 iterations as the limit.

In [19]:

minimisation_problem = MinimizationProblem(reduced_functional, bounds=(T_lb, T_ub))
minimisation_parameters["Status Test"]["Iteration Limit"] = 50

# Define the LinMore Optimiser class:
optimiser = LinMoreOptimiser(
    minimisation_problem,
    minimisation_parameters,
)

minimisation_problem = MinimizationProblem(reduced_functional, bounds=(T_lb, T_ub)) minimisation_parameters["Status Test"]["Iteration Limit"] = 50 # Define the LinMore Optimiser class: optimiser = LinMoreOptimiser( minimisation_problem, minimisation_parameters, )

And again use our callback function to record convergence:

In [20]:

functional_values = []


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


reduced_functional.eval_cb_post = record_value

functional_values = [] def record_value(value, *args): if functional_values: functional_values.append(min(value, min(functional_values))) else: functional_values.append(value) reduced_functional.eval_cb_post = record_value

We next run the optimisation:

In [21]:

optimiser.run()

optimiser.run()

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

  1     2.512942e-01   9.730735e-01   9.970268e-01   2.492567e-01   2         1         4         7         2         0         0         

  2     2.512942e-01   9.730735e-01   2.492567e-01   6.231418e-02   3         1         21        19        2         8         3         

  3     1.398485e-01   7.049145e-01   6.231418e-02   6.231418e-02   4         2         24        23        0         1         3         

  4     1.398485e-01   7.049145e-01   2.131429e-02   5.328573e-03   5         2         39        29        2         10        1         

  5     1.322613e-01   7.120842e-01   5.328573e-03   5.328573e-03   6         3         60        42        0         10        1         

  6     1.252120e-01   7.397731e-01   5.245016e-03   5.328573e-02   7         4         75        50        0         10        1         

  7     1.103025e-01   7.817187e-01   1.080576e-02   5.328573e-01   8         5         89        56        0         10        1         

  8     7.722533e-02   8.109672e-01   2.793263e-02   5.328573e+00   9         6         103       62        0         10        1         

  9     5.846680e-02   4.082384e-01   6.416686e-02   5.328573e+01   10        7         116       67        0         10        1         

  10    3.346408e-02   3.108484e-01   1.055538e-02   5.328573e+02   11        8         129       71        0         10        1         

  11    3.219637e-02   3.217317e-01   1.948306e-03   5.328573e+03   12        9         142       76        0         10        1         

  12    3.092001e-02   3.375663e-01   1.819457e-03   5.328573e+04   13        10        155       81        0         10        1         

  13    2.770807e-02   3.500349e-01   5.920215e-03   5.328573e+05   14        11        168       86        0         10        1         

  14    2.262432e-02   3.478144e-01   1.097739e-02   5.328573e+06   15        12        181       91        0         10        1         

  15    2.016705e-02   3.871509e-01   1.896204e-02   5.328573e+07   16        13        194       96        0         10        1         

  16    1.522114e-02   2.604325e-01   7.444520e-03   5.328573e+08   17        14        207       101       0         10        1         

  17    1.430663e-02   2.299542e-01   3.830272e-03   5.328573e+09   18        15        220       106       0         10        1         

  18    1.341042e-02   2.205176e-01   2.955549e-03   5.328573e+10   19        16        234       112       0         10        1         

  19    1.184265e-02   2.202902e-01   5.342161e-03   5.328573e+11   20        17        248       118       0         10        1         

  20    1.015155e-02   2.030560e-01   6.910755e-03   5.328573e+12   21        18        262       124       0         10        1         

  21    9.472086e-03   1.772844e-01   2.795669e-03   5.328573e+13   22        19        275       129       0         10        1         

  22    8.426315e-03   1.420380e-01   6.207780e-03   5.328573e+14   23        20        288       134       0         10        1         

  23    7.656868e-03   1.374520e-01   5.279115e-03   5.328573e+15   24        21        301       139       0         10        1         

  24    7.193732e-03   1.226151e-01   2.563479e-03   5.328573e+16   25        22        314       144       0         10        1         

  25    6.712520e-03   1.322570e-01   3.433767e-03   5.328573e+17   26        23        327       149       0         10        1         

  26    6.044341e-03   1.383960e-01   5.568237e-03   5.328573e+18   27        24        340       154       0         10        1         

  27    5.460526e-03   1.383758e-01   5.130370e-03   5.328573e+19   28        25        353       159       0         10        1         

  28    5.004609e-03   1.193260e-01   3.516834e-03   1.000000e+20   29        26        366       164       0         10        1         

  29    4.632202e-03   1.225904e-01   3.731863e-03   1.000000e+20   30        27        379       169       0         10        1         

  30    4.413085e-03   1.204403e-01   2.436297e-03   1.000000e+20   31        28        392       174       0         10        1         

  31    3.995208e-03   1.193136e-01   5.093619e-03   1.000000e+20   32        29        406       180       0         10        1         

  32    3.706623e-03   1.362322e-01   4.176357e-03   1.000000e+20   33        30        420       186       0         10        1         

  33    3.525253e-03   1.223796e-01   2.470767e-03   1.000000e+20   34        31        433       191       0         10        1         

  34    3.359774e-03   1.092514e-01   2.513939e-03   1.000000e+20   35        32        447       197       0         10        1         

  35    3.199909e-03   1.106916e-01   2.743700e-03   1.000000e+20   36        33        462       204       0         10        1         

  36    2.991365e-03   1.142644e-01   4.087501e-03   1.000000e+20   37        34        476       210       0         10        1         

  37    2.803523e-03   9.802485e-02   4.350296e-03   1.000000e+20   38        35        490       216       0         10        1         

  38    2.680916e-03   1.150016e-01   4.979110e-03   1.000000e+20   39        36        503       221       0         10        1         

  39    2.626521e-03   8.937871e-02   1.127888e-03   1.000000e+20   40        37        516       226       0         10        1         

  40    2.592479e-03   9.852714e-02   1.452293e-03   1.000000e+20   41        38        529       231       0         10        1         

  41    2.541459e-03   1.064615e-01   2.313853e-03   1.000000e+20   42        39        542       236       0         10        1         

  42    2.445617e-03   1.218980e-01   3.973523e-03   1.000000e+20   43        40        555       241       0         10        1         

  43    2.323316e-03   1.256717e-01   4.376649e-03   1.000000e+20   44        41        568       246       0         10        1         

  44    2.236543e-03   1.121804e-01   1.613827e-03   1.000000e+20   45        42        581       251       0         10        1         

  45    2.125805e-03   9.717157e-02   2.210396e-03   1.000000e+20   46        43        595       257       0         10        1         

  46    2.017510e-03   1.131960e-01   2.517999e-03   1.000000e+20   47        44        610       264       0         10        1         

  47    1.887871e-03   1.109157e-01   2.135416e-03   1.000000e+20   48        45        624       270       0         10        1         

  48    1.659092e-03   9.322535e-02   3.512168e-03   1.000000e+20   49        46        639       277       0         10        1         

  49    1.430427e-03   8.231385e-02   3.378876e-03   1.000000e+20   50        47        652       282       0         10        1         

  50    1.224734e-03   8.324216e-02   4.142060e-03   1.000000e+20   51        48        666       288       0         10        1         
Optimization Terminated with Status: Last Type (Dummy)

And we'll write the functional values to a file so that we can test them.

In [22]:

with open("functional_boundary.txt", "w") as f:
    f.write("\n".join(str(x) for x in functional_values))

with open("functional_boundary.txt", "w") as f: f.write("\n".join(str(x) for x in functional_values))

Let's see how well we have done. At this point a total number of 50 iterations have been performed so let's plot convergence:

In [23]:

active-ipynb

plt.semilogy(functional_values)
plt.xlabel("Iteration #")
plt.ylabel("Functional value")
plt.title("Convergence")

plt.semilogy(functional_values) plt.xlabel("Iteration #") plt.ylabel("Functional value") plt.title("Convergence")

Out[23]:

Text(0.5, 1.0, 'Convergence')

This demonstrates that the functional value decreases by roughly three orders of magnitude over the 50 iterations considered. As with the previous tutorial, the functional value can be reduced further if more iterations are specified, or if the optimisation procedure is configured to continue until a specified tolerance is achieved. We can also visualise the optimised initial condition and compare to the true initial condition:

In [24]:

active-ipynb

fig, axes = plt.subplots(1,2,figsize=[8,4],subplot_kw={'aspect':1.0})
ax1 = tripcolor(T0.block_variable.checkpoint, axes=axes[0], cmap='magma', vmax=0.5)
ax2 = tripcolor(T0_ref, axes=axes[1], cmap='magma', vmax=0.5)
fig.subplots_adjust(right=0.82)
cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.68])
fig.colorbar(ax2,cax=cbar_ax);

fig, axes = plt.subplots(1,2,figsize=[8,4],subplot_kw={'aspect':1.0}) ax1 = tripcolor(T0.block_variable.checkpoint, axes=axes[0], cmap='magma', vmax=0.5) ax2 = tripcolor(T0_ref, axes=axes[1], cmap='magma', vmax=0.5) fig.subplots_adjust(right=0.82) cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.68]) fig.colorbar(ax2,cax=cbar_ax);

Not bad. Not bad at all! Thank you for listening! Crowd. Goes. Wild.