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)[source]¶ 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 firstnvars
rows ofw
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
Examples
>>> from pyapprox.benchmarks.benchmarks import setup_benchmark >>> benchmark=setup_benchmark('advection-diffusion',nvars=2) >>> print(benchmark.keys()) dict_keys(['fun', 'variable'])