Tania Bakhos
Fast Krylov solvers for parabolic PDE-based inverse problems

Institute for Computational and Mathematical Engineering (ICME)
475 Via Ortega
Suite 060
Stanford
CA 94305
taniab@stanford.edu
Arvind K. Saibaba
Peter K. Kitanidis

We propose a fast iterative solver for large-scale PDE-based nonlinear inverse problems, where measurements are used to reconstruct the spatial variation of parameters. Our motivation is transient hydraulic tomography, which is a method to estimate hydraulic parameters related to the subsurface from pressure measurements obtained by a series of pumping tests. Standard approaches for solving the inverse problem require repeated construction of the Jacobian, which represents the sensitivity of the measurements to the unknown parameters. This is prohibitively expensive because it requires repeated solutions of time-dependent parabolic partial differential equations corresponding to multiple sources and receivers. Additionally, solving the differential equations by time-stepping methods requires computing and storing the entire time history. We use a Laplace Transform-based exponential time integrator to independently compute the solution at multiple times, thus reducing the computational cost.

The governing equations for the parabolic PDEs are

$$M(p) \partial_t u(t) + K(p) u(t) = f(t) b \\$$

$$u(0) = u_0 ,$$

where $M$ and $K$ are the mass and stiffness matrices respectively and are the result of discretization by standard linear finite elements. We denote by $p$ the spatially varying parameters we wish to estimate, given discrete pressure head measurements $u$ . The solution can be written in terms of solutions to shifted linear systems

$$(K(p)+s_k M(p) )x_k = M u_0 + \hat{f}(s_k) b, \qquad u_N(t) = \sum_k {w}_k x_k,$$ (1)

where $\hat{f}(s_k)$ denotes the Laplace Transform of $f(t)$ , $u_0$ is the initial condition and $s_k$ , $w_k$ are the $N$ quadrature nodes and weights respectively, chosen corresponding to the modified Talbot contour. The approximate solution at a given time can thus be computed without computing the entire time history. However, solving the shifted sytem of equations can be computationally expensive for large-scale problems.

Krylov subspace methods are particularly attractive to solve these shifted systems of equations (1) because the shift-invariance property of Krylov subspaces allows us to build a single solution space across all shifts and an approximate solution for each shift is then obtained by projecting into a smaller subspace. A single preconditioner of the form $K+\tau M$ , for an appropriately chosen $\tau$ , does not successfully precondition all the systems when the range of values of $s_k$ is large. We consider a flexible Arnoldi algorithm for shifted systems that employs multiple preconditioners of the form $(K+\tau_j M)$ for $j = 1,\dots,m$ . Changing the preconditioner at each iteration allows for better preconditioning of the shifted system across the entire range of shifts.

The performance of our algorithm will be demonstrated on some challenging synthetic examples of large-scale inverse problems arising from transient hydraulic tomography. We demonstrate that the Krylov subspace accelerated Laplace Transform solver provides a significantly cheaper alternative to time-stepping based solvers.