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 0x77195f8ac260>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 1)), 7)

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 0x77195cfb8740>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 187)), 402)

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 0x77195cfb8740>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 187)), 398)

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]:

riesz_map = RieszMap(Q, "L2", constant_jacobian=True)
m = Control(T0, riesz_map=riesz_map)

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)

riesz_map = RieszMap(Q, "L2", constant_jacobian=True) m = Control(T0, riesz_map=riesz_map) 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.07749995459676613

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.07749995459676613

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 0x77195cfb8740>, FiniteElement('Lagrange', triangle, 1), name=None), Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 187)), 400)

In [14]:

print(reduced_functional(T0_ref))

print(reduced_functional(T0_ref))

6.645421118823348e-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.25129415766476565

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         23        20        2         9         3         

  3     1.398485e-01   7.049145e-01   6.231418e-02   6.231418e-02   4         2         28        25        0         2         3         

  4     1.398485e-01   7.049145e-01   2.131429e-02   5.328573e-03   5         2         43        31        2         10        1         

  5     1.322613e-01   7.120842e-01   5.328573e-03   5.328573e-03   6         3         64        44        0         10        1         

  6     1.252120e-01   7.397731e-01   5.245016e-03   5.328573e-02   7         4         79        52        0         10        1         

  7     1.103025e-01   7.817187e-01   1.080576e-02   5.328573e-01   8         5         93        58        0         10        1         

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

  9     5.846680e-02   4.082384e-01   6.416686e-02   5.328573e+01   10        7         120       69        0         10        1         

  10    3.346408e-02   3.108484e-01   1.055538e-02   5.328573e+02   11        8         133       73        0         10        1         

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  29    4.632202e-03   1.225904e-01   3.731862e-03   1.000000e+20   30        27        383       171       0         10        1         

  30    4.413085e-03   1.204403e-01   2.436295e-03   1.000000e+20   31        28        396       176       0         10        1         

  31    3.995209e-03   1.193136e-01   5.093613e-03   1.000000e+20   32        29        410       182       0         10        1         

  32    3.706623e-03   1.362327e-01   4.176371e-03   1.000000e+20   33        30        424       188       0         10        1         

  33    3.525256e-03   1.223798e-01   2.470716e-03   1.000000e+20   34        31        437       193       0         10        1         

  34    3.359781e-03   1.092506e-01   2.513848e-03   1.000000e+20   35        32        451       199       0         10        1         

  35    3.199928e-03   1.106901e-01   2.743464e-03   1.000000e+20   36        33        466       206       0         10        1         

  36    2.991417e-03   1.142631e-01   4.086568e-03   1.000000e+20   37        34        480       212       0         10        1         

  37    2.803616e-03   9.802220e-02   4.348351e-03   1.000000e+20   38        35        494       218       0         10        1         

  38    2.680971e-03   1.149856e-01   4.973379e-03   1.000000e+20   39        36        507       223       0         10        1         

  39    2.626607e-03   8.933875e-02   1.127027e-03   1.000000e+20   40        37        520       228       0         10        1         

  40    2.592494e-03   9.852078e-02   1.452015e-03   1.000000e+20   41        38        533       233       0         10        1         

  41    2.541470e-03   1.063980e-01   2.301562e-03   1.000000e+20   42        39        546       238       0         10        1         

  42    2.445022e-03   1.217209e-01   3.956142e-03   1.000000e+20   43        40        559       243       0         10        1         

  43    2.319965e-03   1.258086e-01   4.309204e-03   1.000000e+20   44        41        572       248       0         10        1         

  44    2.232511e-03   1.108882e-01   1.573250e-03   1.000000e+20   45        42        585       253       0         10        1         

  45    2.082652e-03   1.039674e-01   3.388146e-03   1.000000e+20   46        43        599       259       0         10        1         

  46    2.012684e-03   1.082656e-01   1.258520e-03   1.000000e+20   47        44        613       265       0         10        1         

  47    1.812424e-03   1.073113e-01   3.296413e-03   1.000000e+20   48        45        627       271       0         10        1         

  48    1.478647e-03   1.040070e-01   6.600093e-03   1.000000e+20   49        46        641       277       0         10        1         

  49    1.298538e-03   8.828648e-02   2.784035e-03   1.000000e+20   50        47        654       282       0         10        1         

  50    1.230997e-03   7.933676e-02   1.544047e-03   1.000000e+20   51        48        668       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.