setup_advection_diffusion_source_inversion_benchmark

pyapprox.fenics_models.advection_diffusion_wrappers.setup_advection_diffusion_source_inversion_benchmark(measurement_times=array([0.05, 0.15]), source_strength=0.5, source_width=0.1, true_sample=array([[0.25], [0.75], [4.0], [4.0], [4.0]]), noise_stdev=0.4, max_eval_concurrency=1)

Compute functionals of the following model of transient diffusion of a contaminant

$$\begin{split}\frac{\partial u}{\partial t}(x,t,\rv) + \nabla u(x,t,\rv)-\nabla\cdot\left[k(x,\rv) \nabla u(x,t,\rv)\right] &=g(x,t) \qquad (x,t,\rv)\in D\times [0,1]\times\rvdom\\ \mathcal{B}(x,t,\rv)&=0 \qquad\qquad (x,t,\rv)\in \partial D\times[0,1]\times\rvdom\\ u(x,t,\rv)&=u_0(x,\rv) \qquad (x,t,\rv)\in D\times\{t=0\}\times\rvdom\end{split}$$

Following [MNRJCP2006], [LMSISC2014] we set

$$g(x,t)=\frac{s}{2\pi h^2}\exp\left(-\frac{\lvert x-x_\mathrm{src}\rvert^2}{2h^2}\right)$$

the initial condition as $u(x,z)=0$, $B(x,t,z)$ to be zero Neumann boundary conditions, i.e.

$$\nabla u\cdot n = 0 \quad\mathrm{on} \quad\partial D$$

and we model the diffusivity $k=1$ as a constant.

The quantities of interest are point observations $u(x_l)$ taken at $P$ points in time $\{t_p\}_{p=1}^P$ at $L$ locations $\{x_l\}_{l=1}^L$. The final time $T$ is the last observation time.

These functionals can be used to define the posterior distribution

$$\pi_{\text{post}}(\rv)=\frac{\pi(\V{y}|\rv)\pi(\rv)}{\int_{\rvdom} \pi(\V{y}|\rv)\pi(\rv)d\rv}$$

where the prior is the tensor product of independent and identically distributed uniform variables on $[0,1]$ i.e. $\pi(\rv)=1$, and the likelihood is given by

$$\pi(\V{y}|\rv)=\frac{1}{(2\pi)^{d/2}\sigma}\exp\left(-\frac{1}{2}\frac{(y-f(\rv))^T(y-f(\rv))}{\sigma^2}\right)$$

and $y$ are noisy observations of the solution u at the 9 points of a uniform $3\times 3$ grid covering the physical domain $D$ at successive times $\{t_p\}_{p=1}^P$. Here the noise is indepenent and Normally distrbuted with mean zero and variance $\sigma^2$.

Parameters

measurement_timesnp.ndarray (P)

The times $\{t_p\}_{p=1}^P$ at which measurements of the contaminant concentration are taken

source_strengthfloat

The source strength $s$

source_widthfloat

The source width $h$

true_samplenp.ndarray (2)

The true location of the source used to generate the observations used in the likelihood function

noise_stdevfloat

The standard deviation $sigma$ of the observational noise

max_eval_concurrencyinteger

The maximum number of simulations that can be run in parallel. Should be no more than the maximum number of cores on the computer being used

Returns

benchmarkpya.Benchmark

Object containing the benchmark attributes documented below

funcallable

The quantity of interest $f(w)$ with signature

fun(w) -> np.ndarray

where w is a 2D np.ndarray with shape (nvars+3,nsamples) and the output is a 2D np.ndarray with shape (nsamples,1). The first nvars rows of w are realizations of the random variables. The last 3 rows are configuration variables specifying the numerical discretization of the PDE model. Specifically the first and second configuration variables specify the levels $l_{x_1}$ and $l_{x_2}$ which dictate the resolution of the FEM mesh in the directions ${x_1}$ and ${x_2}$ respectively. The number of cells in the ${x_i}$ direction is given by $2^{l_{x_i}+2}$. The third configuration variable specifies the level $l_t$ of the temporal discretization. The number of timesteps satisfies $2^{l_{t}+2}$ so the timestep size is and $T/2^{l_{t}+2}$.

variablepya.IndependentMultivariateRandomVariable

Object containing information of the joint density of the inputs z which is the tensor product of independent and identically distributed uniform variables on $[0,1]$.

Notes

The example from [MNRJCP2006] can be obtained by setting s=0.5, h=0.1, measurement_times=np.array([0.05,0.15]) and noise_stdev=0.1

The example from [LMSISC2014] can be obtained by setting s=2, h=0.05, measurement_times=np.array([0.1,0.2]) and noise_stdev=0.1

References

MNRJCP2006(1,2)

Youssef M. Marzouk, Habib N. Najm, Larry A. Rahn, Stochastic spectral methods for efficient Bayesian solution of inverse problems, Journal of Computational Physics, Volume 224, Issue 2, 2007, Pages 560-586,

LMSISC2014(1,2)

Jinglai Li and Youssef M. Marzouk. Adaptive Construction of Surrogates for the Bayesian Solution of Inverse Problems, SIAM Journal on Scientific Computing 2014 36:3, A1163-A1186

Examples

>>> from pyapprox.benchmarks.benchmarks import setup_benchmark
>>> benchmark=setup_benchmark('advection-diffusion',nvars=2)
>>> print(benchmark.keys())
dict_keys(['fun', 'variable'])