Exercise: Inverse Modeling with ADCME

Exercise: Inverse Modeling with ADCME

The starter code can be downloaded here.


The thermal diffusivity is the measure of the ease with which heat can diffuse through a material. Let $u$ be the temperature, and $\kappa$ be the thermal diffusivity. The heat transfer process is described by the Fourier law

\[\frac{\partial u(\mathbf{x}, t)}{\partial t} = \kappa\Delta u(\mathbf{x}, t) + f(\mathbf{x}, t), \quad t\in (0,T), x\in \Omega \tag{1}\]

Here $f$ is the heat source and $\Omega$ is the domain.

To make use of the heat equation, we need additional information.

Assume that we want to experiment with a piece of new material. The thermal diffusivity coefficient of the material is an unknown function of the space. Our goal of the experiment is to find out the thermal diffusivity coefficient. To this end, we place some sensors in the domain or on the boundary. The measurements are sparse in the sense that only the temperature from those sensors–-but nowhere else–-are collected. Namely, let the sensors be located at $\{\mathbf{x}_i\}_{i=1}^M$, then we can observe $\{\hat u(\mathbf{x}_i, t)\}_{i=1}^M$, i.e., the measurements of $\{ u(\mathbf{x}_i, t) \}_{i=1}^M$. We also assume that the boundary conditions, initial conditions and the source terms are known.

Problem 1: Parameter Inverse Problem in 1D

We first consider the 1D case. In this problem, the material is a rod $\Omega=[0,1]$. We consider a homogeneous (zero) fixed boundary condition on the right side, and an insulated boundary on the left side. The initial temperature is zero everywhere, i.e., $u(x, 0)=0$, $x\in [0,1]$. The source term is $f(x, t) = \exp(-10(x-0.5)^2)$, and $\kappa(x)$ is a function of space

\[\kappa(x) = a + bx\]

Our task is to estimate the coefficient $a$ and $b$ in $\kappa(x)$. To this end, we place a sensor at $x=0,$ and the sensor records the temperature as a time series $u_0(t)$, $t\in (0,1)$.

Recall that in lecture slide 35/47 of AD.pdf, we formulate the inverse modeling problem as a PDE-constrained optimization problem

\[\begin{aligned} \min_{a, b}\ & \int_{0}^t ( u(0, t)- u_0(t))^2 dt\\ \mathrm{s.t.}\ & \frac{\partial u(x, t)}{\partial t} = \kappa(x)\Delta u(x, t) + f(x, t), \quad t\in (0,T), x\in (0,1) \\ & -\kappa(0)\frac{\partial u(0,t)}{\partial x} = 0, t>0\\ & u(1, t) = 0, t>0\\ & u(x, 0) = 0, x\in [0,1]\\ & \kappa(x) = a x + b \end{aligned}\]

We consider the discretization of the above forward problem. We divide the domain $[0,1]$ into $n$ equispaced intervals. We consider the time horizon $T = 1$, and divide the time horizon $[0,T]$ into $m$ equispaced intervals. We use a finite difference scheme to solve the 1D heat equation Equations (1)–(3). Specifically, we use an implicit scheme for stability:

\[\frac{u^{k+1}_i-u^k_i}{\Delta t} = \kappa_i \frac{u^{k+1}_{i+1}+u^{k+1}_{i-1}-2u^{k+1}_i}{\Delta x^2} + f_i^{k+1},\]
\[k=1,2,\ldots,m, i=1,2,\ldots, n \tag{4}\]

where $\Delta t$ is the time interval, $\Delta x$ is the space interval, $u_i^k$ is the numerical approximation to $u((i-1)\Delta x, (k-1)\Delta t)$, $\kappa_i$ is the numerical approximation to $\kappa((i-1)\Delta x) = a + b(i-1)\Delta x$, and $f_i^{k} = f((i-1)\Delta x, (k-1)\Delta t)$.

For the insulated boundary, we introduce the ghost node $u_0^k$ at location $x=-\Delta x$, and the insulated boundary condition can be numerically discretized by

\[-\kappa_1 \frac{u_2^{k}-u_0^k}{2\Delta x} = 0\tag{5}\]


\[U^k = \begin{bmatrix}u_1^k\\ u_2^k\\ \vdots\\ \\u_n^k\end{bmatrix}\]

The index starts from 1 and ends with $n$. Using the finite difference scheme, together with eliminating the boundary values $u_0^k$, $u_{n+1}^k$, we have the following formula

\[AU^{k+1} = U^k + F^{k+1}\]
  1. Express the matrix $A\in \mathbb{R}^{n\times n}$ in terms of $\Delta t$, $\Delta x$ and $\{\kappa_i\}_{i=1}^{n}$. What is $F^{k+1}\in \mathbb{R}^n$ ?

Hint: Can you eliminate $u_0^k$ and $u_{n+1}^k$ in Eq. (4) using Eq. (5) and $u_{n+1}^k=0$ ?

The starter code Case1D/starter1.jl precomputes the force vector $F^k$ and packs it into a matrix $F\in \mathbb{R}^{(m+1)\times n}$.

  1. Use spdiag[spdiag] to construct A as a SparseTensor (see the starter code for details). spdiag is an ADCME function. See the documentation[spdiag] for the syntax. Here $\kappa$ is given by
\[\kappa(x) = 2+1.5x\]

API Reference

For debugging, check that your $A_{ij}$ is tridiagonal. You can use run(sess, A) to evaluate the SparseTensorA. You should get the following values:


The computational graph of the dynamical system can be efficiently constructed using while_loop.

For debugging, you can plot the temperature on the left side, i.e., $u(0,t)$. You should have something similar to the following plot

Now we are ready to perform inverse modeling.

Problem 2: Function Inverse Problem

Let us consider the function inverse problem, where we do not know the form of $\kappa(x)$. To this end, we substitute $\kappa(x)$ using a neural network. We will use physics constrained learning (PCL) to train $\kappa(x)$ from the temperature data $u_0(x, t)$, where $x$ is the location of a sensor.

Since we do not know the form of $\kappa(x)$, we need more data to solve the inverse problem. Therefore, we assume that we place sensors at the first 25 locations in the discretized grid. The observation data data_pcl.txt is a $(m+1)\times 25$ matrix; each column corresponds to the observation at one sensor.

In this problem, let us parametrize $\kappa(x)$ with a fully connected neural network with 3 hidden layers, 20 neurons per layer, and $\tanh$ activation functions. In ADCME, such a neural network can be constructed using

y = fc(x, [20,20,20,1])

Here x is a $n \times 1$ input, y is a $n \times 1$ output, $[20,20,20,1]$ is the number of neurons per layer (last output layer only has 1 neuron), and fc stands for "fully-connected".

Assume that the neural network is written as $\kappa_\theta(x)$, where $\theta$ is the weights and biases.

Hint: Your curve should look like the following

Hint: You can add noise using

uc = @. uc * (1 + 0.01*randn(length(uc)))
uc = @. uc * (1 + 0.1*randn(length(uc)))

Here @. is for elementwise operations.

Problem 3: Parameter Inverse Problem in 2D

In this problem, we will explore more deeply how ADCME works with numerical solvers. We will use an important technique, custom operators, for incorporating C++ codes. This will be useful when you want to accelerate a performance-critical part, or you want to reuse existing codes. To make the problem simple, the C++ kernel has been prepared for you. For this exercise, you will be working with Case2D/starter.jl.

We consider the 2D case and $T=1 $. We assume that $\Omega=[0,1]^2$. We impose zero boundary conditions on the entire boundary $\Gamma_D=\partial\Omega$. Additionally, we assume that the initial condition is zero everywhere. Two sensors are located at $(0.2,0.2)$ and $(0.8,0.8)$ and these sensors record time series of the temperature $u_1(t)$ and $u_2(t)$. The thermal diffusivity coefficient is a linear function of the space coordinates

\[\kappa(x, y) = a + bx + cy\]

where $a, b$ and $c$ are three coefficients we want to find out from the data $u_1(t)$ and $u_2(t)$.

We use the finite difference method to discretize the PDE. Consider $\Omega=[0,1]\times [0,1]$, we use a uniform grid and divide the domain into $m\times n$ squares, with length $\Delta x$. We also divide $[0,T]$ into $N_T$ intervals of equal length. The implicit scheme for the equation is

\[\frac{u_{ij}^{k+1}-u_{ij}^k}{\Delta t} = \kappa_{ij}\frac{u_{i+1,j}^{k+1}+u_{i,j+1}^{k+1}+u_{i-1,j}^{k+1}+u_{i,j-1}^{k+1}-4u_{ij}^{k+1}}{\Delta x^2} + f_{ij}^{k+1}\]

where $i=2,3,\ldots, m, j=2,3,\ldots, n, k=1,2,\ldots, N_T$.

Here $u_{ij}^k$ is an approximation to $u((i-1)h, (j-1)h, (k-1)\Delta t)$, and $f_{ij}^k = f((i-1)h, (j-1)h, (k-1)\Delta t)$.

We flatten $\{u_{ij}^k\}$ to a 1D vector $U^k$, using $i$ as the leading dimension, i.e., the ordering of the vector is $u_{11}^k, u_{12}^k, \ldots$; We also flatten $f_{ij}^{k+1}$ and $\kappa_{ij}$ as well.

In this problem, we extend the AD framework using custom operators (also known as external function support in the AD community). In the starter code Case2D/starter.jl, we provide a function, heat_equation, a differentiable heat equation solver, which is already implemented for you using C++. By using custom operators, we replace the PDE solver node in the computational graph with our own, more efficient, implementation.

Read the instructions on how to compile the custom operator, and answer the following two questions.

Hint: you should obtain something similar to

The parameters used in this problem are: $m=50$, $n=50$, $T=1$, $N_T=50$, $f(\mathbf{x},t) = e^{-t}\exp(-50((x-0.5)^2+(y-0.5)^2))$, $a = 1.5$, $b=1.0$, $c=2.0$.

The data file data.txt is a $(N_T+1)\times 2$ matrix, where the first and the second columns are $u_1(t)$ and $u_2(t)$ respectively.