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
left, right, bottom, top = 1, 2, 3, 4  # Boundary IDs

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 left, right, bottom, top = 1, 2, 3, 4 # Boundary IDs 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 0x77393a582720>, 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);