Error Analysis of Heat Equation

variational formulation
1

Dolfin version: 2019.1.0
How installed: manual commands from a source repository in terminal
System: Ubuntu 18.04 LTS on Windows Subsystem for Linux. Also replicated in Ubuntu 19.10.

I am solving the heat equation in 1D space, u_t - u_{xx} = f(x,t), on [x_0, X] \times [t_0, T] using space-time CGq (Continuous Galerkin method of degree q) method. I then solve a specific adjoint problem corresponding to a certain quantity of interest that I am interested in using a similar space-time CGq method. Creating the adjoint problem involves choosing right hand side data, initial conditions, and boundary conditions that in the end give me the desired error in the quantity of interest. I then want to calculate the expected error in that quantity of interest.

The way the space-time CGq method works is that we assume the solution can be written as a sum of space-time basis vector functions, which themselves can be written as a product of a scalar spatial function and a scalar time function. In other words, u(x,t) = \sum_{i=0}^q \alpha_i(x)l_i(t), where the l_i(t) are Lagrange polynomials in time. Each subinterval [t_n, t_{n+1}] there are q coupled equations that need to be solved in-tandem to find the spatial coefficient functions. I then update the “initial” solution with the t_{n+1} solution calculated prior. Repeat until final time reached.

I use Fenics to handle the spatial portion of the solve and use my own coded time functions plus Gaussian quadrature in time to handle the time portion.

I perform this method for both the primal(original) problem and the adjoint problem. However, I do increase the degree of interpolation in space and time by 1 for the adjoint problem to avoid orthogonality issues cropping up. For instance, if I use degree 1 in space and degree 2 in time for the primal problem, I will use at least degree 2 in space and degree 3 in time for the adjoint problem.

I have the CGq method working for the primal problem. I have the errors and error convergence rates for different q values to prove it. However, my effectivity ratios, the ratio of the calcuated error in the quantity of interest and the exact error in the quantity of interest, are diverging as my mesh size in space and time decreases. The expected behavior is that they converge to 1 as the mesh size decreases.

I’m not certain if it is my code for the adjoint problem, as I do not have a exact solution to check against, or my error representation code. I’ve been working on the code for several weeks now, but cannot seem to find what the issue is. I’m starting to think there is a quirk here from Fenics that I am not understanding or aware of. Any help in finding out what is wrong is greatly appreciated.

The following is a more in-depth PDF file of the derivation to get from the original heat equation to the error representation form.
https://www.dropbox.com/s/ylf6e9zes4dkocj/Error%20Estimates%20for%20Heat%20Equation%20in%201D.pdf?dl=0

The code that I am using is broken into two files (Requres Fenics 2019.1.0, Python3, and NumPy to be installed):

  1. The file that you run in the command line. It is called HeatEquation1DCGq.py

  2. The file that contains the functions that are called in the first file. This includes the primal and adjoint solvers as well as the error representation solver. It is called HeatEquation1DFunctionFile.py

The code is found here: https://www.dropbox.com/s/fqctx3vr17z9e9b/Code.zip?dl=0

2

Instead of implementing the adjoint equation and functional derivative by hand, you could use dolfin_adjoint to automatically compute the adjoint and derivative by overloading your forward problem. Note that the adjoint equation uses the same function space as the state variable in pyadjoint. However, it might point you in the correct direction.

3

I’m not exactly seeing how dolfin_adjoint works or would help me. I’m trying to a priori, create an estimate of the error in a quantity of interest. I looked at the link you provided, but can’t seem to see anything relating to solving an adjoint that I understand. Is there a better example of how dolfin replaces work I am doing? Also, I need to solve the adjoint problem in such a way that I can evaluate at any point between any two mesh nodes. Does this method allow that?

4

I believe the issue might lie with my adjoint problem. The problem I’m solving is
\begin{align} \phi_t + \phi_{xx} &= 0\\ \phi(x, T) &= \phi_0(x)\\ \phi(0, t) &= \phi(1, t) = 0 \end{align}

Where \phi_0(x) has compact support on [0,1]. Physically, I expect the solution to decay to 0 due to the boundary conditions and lack of a source term. My CGq solve does that, but it is very oscillatory in the center. When compared with the Backward Euler solution which is smooth on the entire domain at each time step.

I’ve double-checked my adjoint code again and still haven’t been able to find an issue.

I cut down the adjoint code into it’s own file. Maybe it will make it easier to spot what’s the problem.

Adjoint solved using backward euler: https://www.dropbox.com/s/9nr20v4zyts1af4/HeatEquation1DAdjointBE.py?dl=0

Adjoint solved using continuous galerkin in space-time: https://www.dropbox.com/s/9nr20v4zyts1af4/HeatEquation1DAdjointBE.py?dl=0