Readme file for "Information density-based mesh refinement"

@note This program implements the ideas and algorithms described in the paper "Estimating and using information in inverse problems" by Wolfgang Bangerth, Chris R. Johnson, Dennis K. Njeru, and Bart can Bloemen Waanders, 2022. See there for more information.

Motivation

Inverse problems are problems where we would like to infer properties of a system from measurements of the system's state or response to external stimuli. The specific example this program addresses is that we want to identify the source term (i.e., right hand side function) in an advection-diffusion equation from point measurements of the solution of the equation. A typical application is that we would like to find out the locations and strengths of pollution sources based on measuring the concentration of the polluting substance at a number of points.

It is clear that in order to solve such problems, one needs to "know" something about the system's state (here: the pollution concentration) through measurements. Intuitively, it is also clear that we know "more" about the pollution sources by (i) measuring at more points, and (ii) by measuring downstream from the sources than we would if we had measured upstream. Intuitive concepts such as this motivate wondering whether we can define an "information density" function whose value at a point $\mathbf x$ describes how much we know about potential sources located at $\mathbf x$.

The paper which this code accompanies explores the concept of information in inverse problems. It defines an "information density" by solving auxiliary problems for each measurement, and then outlines possible applications for these information densities in three vignettes: spatially variable regularization; mesh refinement; and optimal experimental design. It then considers one of these in detail through numerical experiments, namely mesh refinement. This program implements the algorithms shown there and produces the numerical results.

To run the code

After running cmake and compiling via make (or, if you have used the -G ... option of cmake, compiling the program via your favorite integrated development environment), you can run the executable by either just saying make run or using ./mesh_refinement on the command line. The default is to compile in "debug mode"; you can switch to "release mode" by saying make release and then compiling everything again.

The program contains a switch that decides which mesh refinement algorithm to use. By default, it refines the mesh based on the information criterion discussed in the paper; it runs a sequence of 7 mesh refinement cycles. In debug mode, running the program as is takes about 50 CPU minutes on a reasonably modern laptop. (The program takes about five and a half minutes in release mode.) It parallelizes certain operations, so the actual run time may be shorter depending on how many cores are available.

For each cycle, it outputs the solution as a VTU file, along with the $A$, $B$, $C$, and $M$ matrices discussed in the paper. These matrices can then be used to compute the eigenvalues of the $H$ matrix defined by $H = B^T A^{-T} C A^{-1} B + \beta M$ where $\beta$ is the regularization parameters.

Some of the pictures shown in the paper are also reproduced as part of this code gallery program. See the paper for captions and more information.