This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Stochastic estimation of Green’s functions with application to diffusion and advection-diffusion-reaction problems

Russell G. Keanini Jerry Dahlberg Philip Brown Mehdi Morovati Hamidreza Moradi Donald Jacobs Peter T. Tkacik
Abstract

A stochastic method is described for estimating Green’s functions (GF’s), appropriate to linear advection-diffusion-reaction transport problems, evolving in arbitrary geometries. By allowing straightforward construction of approximate, though high-accuracy GF’s, within any geometry, the technique solves the central challenge in obtaining Green’s function solutions. In contrast to Monte Carlo solutions of individual transport problems, subject to specific sets of conditions and forcing, the proposed technique produces approximate GF’s that can be used: a) to obtain (infinite) sets of solutions, subject to any combination of (random and deterministic) boundary, initial, and internal forcing, b) as high fidelity direct models in inverse problems, and c) as high quality process models in thermal and mass transport design, optimization, and process control problems. The technique exploits an equivalence between the adjoint problem governing the transport problem Green’s function, G​(𝐱,t|𝐱′,t′),G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right), and the backward Kolmogorov problem governing the transition density, p​(𝐱,t|𝐱′,t′),p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right), of the stochastic process used in Green’s function construction. We address nonspecialists and report four contributions. First, a recipe is outlined for diagnosing when stochastic Green’s function estimation can be used, and for subsequently estimating the transition density and associated Green’s function. Second, a naive estimator for the transition density is proposed and tested. Third, Green’s function estimation error produced by random walker absorption at Dirichlet boundaries is suppressed using a simple random walker splitting technique. Last, spatial discontinuity in estimated GF’s, produced by the naive estimator, is suppressed using a simple area averaging method. The paper provides guidance on choosing key numerical parameters, and the technique is tested against two simple unsteady, linear heat conduction problems, and an unsteady groundwater dispersion problem, each having known, exact GF’s.

keywords:
stochastic estimation of Green’s functions , linear and nonlinear transport
PACS:
0000 , 1111
MSC:
0000 , 1111
††journal: Journal of Computational Physics
\affiliation

[inst1]organization=Department of Mechanical Engineering and Engineering Science,addressline=9201 University City Blvd., city=Charlotte, postcode=28223, state=North Carolina, country=United States

\affiliation

[inst2]organization=University of Tennessee Space Institute, city=Huntsville, postcode=35649, state=Alabama, country=United States

\affiliation

[inst3]organization=Department of Physics and Optical Science,addressline=9201 University City Blvd., city=Charlotte, postcode=28223, state=North Carolina, country=United States

1 Introduction

Considering a physical system whose dynamics are determined by a set of conservation laws and physical principles, finding a Green’s function that reliably models the system’s spatial and temporal response to forcing opens the door to a number of powerful capabilities. Focusing on thermal and chemical transport in complicated, irregular regions, the temporal thermal or chemical response at discrete locations, to an infinite array of space and time-dependent forcing, for example, can be readily determined [1, 2]. Similarly, Green’s functions can be incorporated in inverse problems, allowing, for example, extraction of thermal and pollutant source locations and strengths, based on limited temperature and chemical measurements; see references and examples in [3, 4, 5]. Green’s functions can also provide a basis for efficient thermal and mass transport design, optimization, and process control [6, 7, 8]. In these applications, a reliable Green’s function can replace approximate or empirical system models, improving the reliability and quality of a design or control process.

Unfortunately, building a Green’s function rests on three typically challenging tasks:

A) derive (or, when possible, look up) the adjoint equation governing the Green’s function (GF);


B) simultaneously derive the integral solution for the evolving field variable, η​(𝐱,t)\eta\left(\mathbf{x},t\right) - which we’ll call the ’magic rule’ [2] - stated in terms of convolutions of the Green’s function with initial and boundary conditions, as well as in-domain forcing terms, f​(𝐱,t);f(\mathbf{x},t); and


C) find an appropriate Green’s function, subject to the governing adjoint problem derived in steps A) and B).

Most texts treating Green’s function techniques do not illustrate the trial and error nature of the first two steps, leaving non-specialists with limited insight. We address this deficiency in Appendix C, simultaneously deriving, via trial and error, the adjoint equation and magic rule associated with a linear advection-diffusion problem.

Importantly, most of the difficulty surrounding Green’s function solutions is tied to task C). A number of texts illustrate this step, invariably focusing on solution of (given, non-derived) adjoint equations, stated in simple geometries, leading to analytical Green’s function solutions [1, 2, 9]. For linear transport problems in complex geometries, however, this step clearly requires a numerical attack. Unfortunately, little work has been reported on numerical construction of Green’s function’s, specifically, numerical solution of adjoint equations associated with transport problem-specific partial differential equations, such as Eqs. (1) and (2) below. While boundary element (BEM) [10, 11, 12], Green element (GEM) [13, 14, 15], and finite element methods (FEM) [16, 17, 18, 19] are widely characterized as Green’s function techniques, none of these explicitly solve problem-specific adjoint equations. As a consequence, these only allow whole-domain solutions, in contrast to the pointwise solutions obtained by the technique proposed here.

We focus on construction of approximate GF’s, appropriate to solution of advection-diffusion-reaction problems of the general form:

𝐃:∇T⊗∇η−𝐯⋅∇η−γ​η−∂η∂t=−f​(𝐱,t),\mathbf{D}\mathbf{:}\nabla^{T}\otimes\nabla\eta-\mathbf{v}\cdot\nabla\eta-\gamma\eta-\frac{\partial\eta}{\partial t}=-f(\mathbf{x},t), (1)

subject to specified boundary and initial conditions. Here, η=η​(𝐱,t),\eta=\eta\left(\mathbf{x},t\right), is a scalar field property (often a density of some kind), subject to: i) transport by diffusion/dispersion (𝐃(\mathbf{D} real-valued ) and/or ii) wave propagation (𝐃(\mathbf{D} complex-valued), iii) advective transport by a fluid velocity field, 𝐯=𝐯​(𝐱,t),\mathbf{v}=\mathbf{v}\left(\mathbf{x},t\right), iv) depletion by first order chemical reaction, −γ​η,-\gamma\eta, and v) augmentation (f​(𝐱,t)>0)\left(f\left(\mathbf{x},t\right)>0\right) or depletion (f​(𝐱,t)<0)\left(f\left(\mathbf{x},t\right)<0\right) by a source/sink field, f​(𝐱,t).f\left(\mathbf{x},t\right). In addition, 𝐃\mathbf{D} represents an anisotropic, often position- and time-dependent diffusion or dispersion tensor, γ\gamma is a solute-dependent chemical or radiological decay constant, where 𝐃,\mathbf{D}, 𝐯,\mathbf{v}, and γ\gamma have no dependence on η,\eta, while ff is at most, linearly dependent on η.\eta. Importantly, the dispersion/diffusivity, velocity and source/sink fields are known quantities.

In situations where the dispersion/diffusion tensor, 𝐃,\mathbf{D}, is diagonal - capturing diffusive or dispersive transport that’s dominant in one, two, or three mutually orthogonal directions, the substitutions η→exp⁡(−γ​t)​η\eta\rightarrow\exp\left(-\gamma t\right)\eta and f​(𝐱,t)→exp⁡(−γ​t)​f​(𝐱,t)f\left(\mathbf{x},t\right)\rightarrow\exp\left(-\gamma t\right)f\left(\mathbf{x},t\right), can be used to transform Eq. (1) into an equation of the form,

𝐃:∇T⊗∇η−𝐯⋅∇η−∂η∂t=−f​(𝐱,t),\mathbf{D}\mathbf{:}{\nabla^{T}}\otimes\nabla\eta-\mathbf{v}\cdot\nabla\eta-\frac{\partial\eta}{\partial t}=-f(\mathbf{x},t), (2)

As discussed in the third validation test below, equations of the form (1) arise, for example, in models of solute transport in aquifers [20]. Equation (2) appears in many areas of physics, biology, and engineering, including fluid mechanics (linearized Navier-Stokes equations) [21, 22], drift-diffusion [23] and reaction-diffusion problems [24, 25]), quantum mechanics (Schrödinger equation) [26], statistical mechanics (Fokker-Planck equation) [27] and plasma physics (Vlasov equation) [28]. Wave equations [29, 30] and Chapman-Kolmogorov equations [31] can also be expressed in this form.

Since ’backward time’ stochastic construction of GF’s provides the single-point response at a chosen ’response point’, (𝐱,t),\left(\mathbf{x},t\right), as produced by any given set of boundary and initial conditions, we focus most of the discussion on backward time Green’s function construction. For testing and validation purposes, however, the last validation case uses forward time stochastic construction. As discussed in test 3, forward time construction: i) can be used when the diffusion/dispersion tensor is diagnonal: 𝐃=𝐃~⋅𝐈,\mathbf{D}=\mathbf{\tilde{D}}\cdot\mathbf{I}, where 𝐃~\mathbf{\tilde{D}} is a tensor function of position and time, and 𝐈\mathbf{I} is the identity matrix, and ii) provides the response at (𝐱,t),\left(\mathbf{x},t\right), as produced by an impulse at a specific ’impulse point’, (𝐱′,t′).\left(\mathbf{x^{\prime}},t^{\prime}\right).

Focusing on backward time construction, the proposed approach launches a swarm of RW’s backward in time from a chosen response point, (𝐱,t),\left(\mathbf{x},t\right), calculating the random, backward trajectories of each random walker, using a transport-problem-specific stochastic differential equation. The transition (probability) density, p​(𝐱,t|𝐱′,t′),p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right), for observing the random walker at a chosen ’impulse point’, (𝐱′,t′),\left(\mathbf{x}^{\prime},t^{\prime}\right), given that it started at (𝐱,t),\left(\mathbf{x},t\right), is then estimated. Under conditions where the transport problem’s adjoint problem - governing evolution of the propagator, K​(𝐱,t|𝐱′,t′),K\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right), - has the same structure as the so-called backward Kolmogorov problem - governing evolution of p​(𝐱,t|𝐱′,t′)p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right) - the estimated transition density allows estimation of the propagator, and in turn, the Green’s function.

In overview, we first outline a step-by-step recipe for stochastically constructing GF’s, appropriate to transport problems governed by Eqs. (1) and (2). The proposed stochastic Green’s function estimation technique is then contrasted against Monte Carlo approaches for (approximately) solving initial-boundary value problems (IBVP’s) [32, 33, 34, 35, 36, 37]: The former bypasses complicated, often impossible derivation of stochastic representative solutions required for the latter, and, in addition, provides an approximate input-response (GF) function, applicable to generic IBVP’s, subject to any combination of Dirichlet, Neumann, and/or Robin BC’s. A simple naive estimator [38, 39, 40, 41, 42] for the transition density, p​(𝐱,t|x′,t′),p\left(\mathbf{x},t|\mathrm{x^{\prime}},t^{\prime}\right), is then derived and heuristically validated using a central limit argument. Numerical validation tests reveal two significant sources of error: depletion of the initial random walker swarm due to absorption at Dirichlet boundaries, and spatial discontinuity in estimated Green’s function’s, produced by our naive density estimator. To solve the first problem, we propose a simple particle splitting technique which preserves the number of simulated random walker’s launched, as well as satisfying continuum mass conservation at the splitting point. A simple area-averaging technique is proposed for tackling the second problem. In order to stochastically estimate GF’s, three numerical parameters must be chosen: the number of random walker’s launched, N𝐱,t,N_{\mathbf{x},t}, from a chosen response point, the interrogation area, Δ​A′​(𝐱′,t′),\Delta A^{\prime}\left(\mathbf{x^{\prime}},t^{\prime}\right), on which random walker’s are captured, and the time step, Δ​t′,\Delta t^{\prime}, used for estimating random walker trajectories. Strategies for choosing these parameters are outlined in S1 Appendix A. Finally, the Green’s function estimation technique is tested against three test problems, of increasing complexity, each having known GF’s.

2 Recipe: Stochastic estimation of Green’s functions

Using stochastic methods to estimate GF’s requires a number of analytical and numerical steps. In order to provide an overview of the proposed technique, we outline the steps here.

Step 1: Simultaneously determine the adjoint problem governing the Green’s function and the magic rule: Using Green’s identities, and by trial and error, simultaneously determine the adjoint equation, the magic rule [1, 2], and the boundary and initial conditions on the adjoint equation. We illustrate this step in Appendix C, using a simplified, constant coefficient advection-diffusion problem.

Step 2: Postulate that continuum transport reflects microscale stochastic dynamics: This step provides a physical connection between microscale and continuum transport, as well as the mathematical basis for Steps 4 and 5. Express η​(𝐱′,t′)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right) as a functional, η​(𝝌​(s′),s′),\eta\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right), of a Brownian stochastic process, 𝝌​(s′)\boldsymbol{\chi}\left(s^{\prime}\right) [31, 43, 44],

η=η​(χ​(s′),s′)\eta=\eta\left(\mathbf{\chi}\left(s^{\prime}\right),s^{\prime}\right) (3)

evolving in backward time, s′,s^{\prime}, according to

d​𝝌​(s′)=𝐛​(𝝌​(s′),s′)​d​s′+𝐁​(𝝌​(s′),s′)⋅d​𝐰​(s′)d\boldsymbol{\chi}\left(s^{\prime}\right)=\mathbf{b}\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right)ds^{\prime}+\mathbf{B}\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right)\cdot d\mathbf{w}\left(s^{\prime}\right) (4)

where 𝐰​(s′)\mathbf{w}\left(s^{\prime}\right) is an n-dimensional Weiner process, 𝐛​(𝐱′,s′)\mathbf{b}\left(\mathbf{x^{\prime}},s^{\prime}\right) is the reversed drift/velocity field, 𝐛​(𝐱′,s′)=−𝐯​(𝐱′,s′),\mathbf{b}\left(\mathbf{x^{\prime}},s^{\prime}\right)=-\mathbf{v}\left(\mathbf{x^{\prime}},s^{\prime}\right), and where the diffusion matrix is tied to 𝐁\mathbf{B} by

𝐃​(𝐱′,s′)=⟨𝐁​(𝝌​(s′),s′)⋅𝐁𝐓​(𝝌​(s′),s′)⟩/2\mathbf{D}\left(\mathbf{x^{\prime}},s^{\prime}\right)=\langle\mathbf{B}\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right)\cdot\mathbf{B^{T}}\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right)\rangle/2 (5)

with the average representing the conditional expectation E𝐱,s​𝐁​(𝝌​(s′),s′)⋅𝐁𝐓​(𝝌​(s′),s′).\mathrm{E_{\mathbf{x},s}}\mathbf{B}\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right)\cdot\mathbf{B^{T}}\left(\boldsymbol{\chi}\left(s^{\prime}\right),s^{\prime}\right). In addition, the relationships between backward and forward (actual) time instants are as follows: i) the chosen, fixed (actual) response time, t,t, corresponds to the fixed backward instant, s,s, from which the random walker swarm is launched: t⟷s;t\longleftrightarrow s; ii) the variable (actual) impulse time, t′,t^{\prime}, corresponds to the variable backward random walker capture time, s′:s^{\prime}: t′⟷s′;t^{\prime}\longleftrightarrow s^{\prime}; iii) the magnitude of forward time difference, t−t′≥0,t-t^{\prime}\geq 0, equals the backward time difference, s′−s≥0.s^{\prime}-s\geq 0.

Step 3: Determine the backward Kolmogorov problem governing backward evolution of the transition density: The probability of observing a backward moving random walker, 𝝌​(s′),\boldsymbol{\chi}\left(s^{\prime}\right), governed by (4), at (𝐱′,s′),\left(\mathbf{x^{\prime}},s^{\prime}\right), given that it started at (𝐱,s),\left(\mathbf{x},s\right), is given by the probability density function, p​(𝐱,s|𝐱′,s′).p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right). As discussed in Step 4, the present step is required in order to ensure that the adjoint problem governing the Green’s function coincides with the terminal-boundary value problem governing p​(𝐱,s|𝐱′,s′).p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right).

The backward time evolution of p​(𝐱,s|𝐱′,s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) is governed by the backward Kolmogorov equation [31, 43, 44]:

∂p​(𝐱,s|𝐱′,s′)∂s′+𝐛​(𝐱′,s′)⋅∇′p​(𝐱,s|𝐱′,s′)+Di​j​∂2p​(𝐱,s|𝐱′,s′)∂xi′​∂xj′=0\frac{\partial p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)}{\partial s^{\prime}}+\mathbf{b}\left(\mathbf{x^{\prime}},s^{\prime}\right)\cdot\nabla^{\prime}p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)+D_{ij}\frac{\partial^{2}p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)}{\partial x_{i}^{\prime}\partial x_{j}^{\prime}}=0 (6)

where p​(𝐱,s|𝐱′,s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) is subject to the terminal backward time condition:

p​(𝐱,s|𝐱′,s′)→δ​(𝐱−𝐱′)s′→sp\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)\rightarrow\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\hskip 19.91684pts^{\prime}\rightarrow s (7)

Boundary conditions on (6) that are relevant to, and correspond respectively to Dirichlet, Neumann, and Robin BC’s in typical advection-diffusion-reaction problems, are as follows [43, 44]:

p​(𝐱,s|𝐱′,s′)=0​on​δ​ΩD×[s,T)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=0\hskip 8.5359pt\mathrm{on}\hskip 8.5359pt\delta\Omega_{D}\times[s,T) (8)
∂p​(𝐱,s|𝐱′,s′)∂n~​(𝐱′,s)=0​on​δ​ΩN×[s,T)\frac{\partial p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)}{\partial\tilde{n}\left(\mathbf{x^{\prime}},s\right)}=0\hskip 8.5359pt\mathrm{on}\hskip 8.5359pt\delta\Omega_{N}\times[s,T) (9)

and

−∂p​(𝐱,s|𝐱′,s′)∂n~​(𝐱′,s)=κ​(𝐱′,s′)​p​(𝐱,s|𝐱′,s′)​on​δ​ΩR×[s,T)-\frac{\partial p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)}{\partial\tilde{n}\left(\mathbf{x^{\prime}},s\right)}=\kappa\left(\mathbf{x^{\prime}},s^{\prime}\right)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)\hskip 8.5359pt\mathrm{on}\hskip 8.5359pt\delta\Omega_{R}\times[s,T) (10)

In Eq. (8), δ​ΩD×[s,T)\delta\Omega_{D}\times[s,T) denotes the space-time boundary on which Dirichlet BC’s are imposed on the physical transport problem. Similar definitions for Neumann and Robin boundaries are given respectively in (9) and (10). The conormal vector is defined as [44] n~​(𝐱′,s)=𝐃​(𝐱′,s)⋅𝐧​(𝐱′),\tilde{n}\left(\mathbf{x^{\prime}},s\right)=\mathbf{D}\left(\mathbf{x^{\prime}},s\right)\cdot\mathbf{n}\left(\mathbf{x^{\prime}}\right), where 𝐧​(𝐱′)\mathbf{n}\left(\mathbf{x^{\prime}}\right) is the local unit outward normal. Finally, in mass transfer problems, κ​(𝐱′,s′)\kappa\left(\mathbf{x^{\prime}},s^{\prime}\right) is the local reaction coefficient [44], arising, for example, on reactive boundaries. In heat transfer problems, κ​(𝐱′,s′)\kappa\left(\mathbf{x^{\prime}},s^{\prime}\right) corresponds to the local convective heat transfer coefficient, divided by the local thermal conductivity.

Step 4: Determine the applicability of the Green’s function estimation method: For transport problems in which diffusive or dispersive transport is non-negligible, stochastic estimation of GF’s can only be used if the adjoint problem governing the so-called propagator, K​(𝐱,s|𝐱′,s′),K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), coincides with the backward Kolmogorov problem, governed by Eq. (6) and problem-specific combinations of the BC’s in Eqs. (8) through (10). Specifically:

a) the diffusion/dispersion matrix must be diagonal:

𝐃=𝐃~​(𝐱,s)⋅𝐈\mathbf{D}=\mathbf{\tilde{D}}\left(\mathbf{x},s\right)\mathbf{\cdot}\mathbf{I} (11)

This requirement emerges in Step 1: Non-diagonal 𝐃\mathbf{D} produces a convolution integral, over the problem domain, Ω,\Omega, in two unknowns, the Green’s function and η​(𝐱′,s′).\eta\left(\mathbf{x^{\prime}},s^{\prime}\right).

b) The spatial domain, Ω,\Omega, in which transport takes place, must remain fixed/non-deforming; this ensures that time-invariance and symmetry properties intrinsic to the propagator are enforced [2].

Under conditions a) and b), the problem governing the propagator, K​(𝐱,s|𝐱′,s′),K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), mirrors the problem governing p​(𝐱,s|𝐱′,s′),p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), with K​(𝐱,s|𝐱′,s′)K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) replacing p​(𝐱,s|𝐱′,s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) in Eqs. (6) through (10). Given an estimated propagator, the associated Green’s function is then given by [2]:

G​(𝐱,s|𝐱′,s′)=H​(s′−s)​K​(𝐱,s|𝐱′,s′)G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=H\left(s^{\prime}-s\right)K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) (12)

where H​(s′−s)H\left(s^{\prime}-s\right) is the Heaviside function.

Step 5: Choose a scheme for integrating Eq. (4) and an estimator for estimating p​(𝐱,s|𝐱′,s′):p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right): Various approaches, each characterized by a pointwise (strong) order of convergence, are available for numerically integrating stochastic differential equations, here Eq. (4) [45, 46, 47, 48]. Since boundary conditions on the adjoint equation are all-determining in constructing the Green’s function, and since interaction of stochastic processes with Dirichlet, Neumann, and Robin boundaries, as integrated by the Euler-Maruyama scheme [45, 47, 49] have been well-studied [44], this scheme is recommended. We note below, however, alternative, less expensive integration schemes.

Estimation of probability density functions, here p​(𝐱,s|𝐱′,s′),p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), is a well-developed field [38, 39, 40, 41, 42]. In this paper, as detailed below, we use a simple naive estimator [38, 40]. Limitations and alternatives to this approach are highlighted below.

3 Comparison of proposed Green’s function estimation method with Monte Carlo solutions of initial-boundary value problem; limitations

Monte Carlo (MC) techniques for solving linear partial differential equations, subject to problem-specific boundary and initial conditions, have a long history [50, 51, 52]. By contrast, stochastic construction of GF’s via transition density estimation appears to represent a new and more general approach, where the estimated Green’s function can be used to (approximately) solve a given linear PDE, subject to any combination of initial and boundary conditions. In order to contrast the proposed stochastic Green’s function estimation method against MC solutions of IBVP’s [32, 33, 34, 35, 36, 53], we outline the five-step construction of MC solutions, and then touch on the significant differences and advantages introduced by the present approach.

Monte Carlo solutions are based on representative stochastic solutions of elliptic, parabolic, and hyperbolic PDE’s, subject to specific boundary conditions, and in the latter two cases, initial conditions [54, 55, 56, 57, 58]. Focusing on parabolic problems governed, for example, by Eqs. (1) or (2), the recipe is as follows: i) Assume that macroscopic evolution of the field variable, η​(𝐱,t),\eta\left(\mathbf{x},t\right), reflects (long time scale, collective [21, 44]) microscale advective-diffusive transport of some conserved real or virtual property, e.g., microscopic number, momentum, or energy density [21, 59, 60], represented by 𝝌​(s),\boldsymbol{\chi}\left(s\right), and modeled by Eq. (4). ii) Using Ito’s formula [31, 43, 44], determine the differential, backward-time change in η,\eta, produced by the stochastic evolution of 𝝌​(s).\boldsymbol{\chi}\left(s\right). iii) In the equation derived in step ii), and referring to (1), replace the temporal derivative of η,\eta, along with advective and diffusive terms with remaining source and reaction terms, f​(𝐱,t)f\left(\mathbf{x},t\right) and γ​η​(𝐱,t),\gamma\eta\left(\mathbf{x},t\right), respectively. iv) Attempt to obtain a stochastic representative solution for η​(𝐱,t),\eta\left(\mathbf{x},t\right), by integrating, backward in time, the differential equation obtained in iii). In deriving the representative solution, account for the stochastic process sampling Dirichlet, Neumann, and/or Robin boundaries that appear in the continuum transport problem [55, 56, 57, 58]. v) Using the representative solution obtained in step iv), launch a random walker swarm from a chosen solution point, (𝐱,t),\left(\mathbf{x},t\right), and use appropriate absorption (Dirichlet boundaries), reflection (Neumann boundaries), and partial reflection (Robin boundaries) techniques [44] to construct a Monte Carlo estimate for η​(𝐱,t).\eta\left(\mathbf{x},t\right).

Comparing the MC approach with the proposed Green’s function estimation procedure, two significant advantages emerge. First, by directly estimating transition densities p​(𝐱,s|𝐱′,s′),p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), Green’s function estimation bypasses the challenging task of deriving representative stochastic solutions [55, 56, 57, 58]: Steps iii) and iv) are eliminated in favor of straightforward stochastic integration of Eq. (4), combined with well-developed boundary interaction techniques [44]. Second, as noted above, the technique produces an estimated, non-problem-specific input-response function (Green’s function) that can be used to solve a family of advection-diffusion-reaction problems, subject to any combination of Dirichlet, Neumann, and/or Robin BC’s.

3.1 Challenges and a limitation

The present work is experimental, focused on developing a stochastically-based Green’s function estimation technique. Prior to looking at details, we highlight two significant challenges that required solution, and note a limitation associated with the method in its current form.

First, following launch of a random walker swarm - consisting of N𝐱,sN_{\mathbf{x},s} integrated realizations of Eq. (4) - from a chosen solution point, (𝐱,s),\left(\mathbf{x},s\right), absorption of individual realizations at Dirichlet boundaries rapidly depletes the initial set of N𝐱,sN_{\mathbf{x},s} random walkers. Depletion, in turn, produces increasingly degraded-in-time Green’s function estimates. Our solution, based on respawning of weighted random walker’s, and carried out at random locations in the problem domain, enforces continuum conservation of mass at the particle respawning point, and significantly improves Green’s function estimation accuracy. See Section 6.3. Second, spatial variance in estimated GF’s likely arises due to use of a spatially discontinuous naive density estimator [38, 40], Eq. (13) below. We address this problem by introducing an area-averaging technique in which variance in local Green’s function estimates is minimized by iteratively altering the size of a smoothing window. See Section 6.4. Although not investigated, kernel estimators, which also eliminate discontinuity in naive estimators [38, 40], might also prove advantageous.

Regarding the present technique’s main limitation, Euler-Murayama integration of Eq. (4) has a low Δ​s1/2\Delta s^{1/2} order of strong convergence [45, 46, 47]. For a given degree of Green’s function estimation accuracy, significantly faster integration can be achieved using, for example, Milstein’s order 1 method [46], or order 2 strong Taylor expansion methods [47], or order 2 Runge-Kutta methods [48]. However, again, since boundary conditions on the adjoint problem are central to constructing the estimated Green’s function, and since treatment of random walker-boundary interactions has reached maturity in the case of Euler-Murayama integration [44], this scheme is used.

4 A simple estimator for the transition density

Transition density estimation remains a research-level task [38, 39, 40, 41, 42]. An accessible overview of theoretical results, focused on the simplest case of estimating ’free-space’ transition densities, po​(𝐱,s|𝐱′,s′),p_{o}\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), can be found in [39].

In this study, the transition density is estimated using the naive estimator [38]:

p​(𝐱,s|𝐱′,s′)=n​(𝐱′,s′)N𝐱,s​Δ​A′​(𝐱′,s′)+O​(Δ​s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=\frac{n\left(\mathbf{x^{\prime}},s^{\prime}\right)}{N_{\mathbf{x},s}\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)}+O\left({\Delta s^{\prime}}\right) (13)

where, the O​(Δ​s′)O\left({\Delta s^{\prime}}\right) estimation error is determined by the Euler-Murayama scheme’s weak convergence order [31, 46, 47], and where N𝐱,s,N_{\mathbf{x},s}, Δ​A′​(𝐱′,s′),\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), and n​(𝐱′,s′)n\left(\mathbf{x^{\prime}},s^{\prime}\right) are, respectively, the number of RW’s launched from the solution point, (𝐱,s),\left(\mathbf{x},s\right), the chosen RW interrogation area, centered on 𝐱′,\mathbf{x^{\prime}}, at backward time, s′,s^{\prime}, and the number of RW’s that reach or pass through Δ​A′​(𝐱′,s′),\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), over the time interval, [s′,s′+Δ​s′).[s^{\prime},s^{\prime}+\Delta s^{\prime}).

While the formula in Eq. (13) is stated in terms of the Lebesque measure in, e.g., [40], it can be derived heuristically as follows. Focus on two-dimensional random walks, as used in our tests, evolving within an arbitrary two-dimensional domain, Ω,\Omega, subject to any combination of Dirichlet, Neumann and/or Robin boundaries. For two-dimensional problems, it is useful to picture a three-dimensional, cylindrical hyperspace in which the downward (upward) vertical direction corresponds to the backward (forward) time axis, s′(↔t′),s^{\prime}\left(\leftrightarrow t^{\prime}\right), the top corresponds to the backward (forward) solution time slice, s(↔t),s\left(\leftrightarrow t\right), the bottom to the backward (forward) terminal (initial) time slice, with any horizontal slice through the hyperspace corresponding to an instantaneous snapshot of the two-dimensional spatial domain, Ω.\Omega.

Symbolically, recognize that the probability of random walker’s, 𝝌​(s′),\boldsymbol{\chi}\left(s^{\prime}\right), launched from (𝐱,s),\left(\mathbf{x},s\right), reaching a chosen interrogation area, Δ​A′​(𝐱′,s′),\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), at the instant s′,s^{\prime}, is given by:

P​[𝝌​(s′)∈Δ​A′​(𝐱′,s′)]=∫Δ​A′​(𝐱′,s′)p​(𝐱,s|𝐱​”,s′)​𝑑𝐱​”P\left[\boldsymbol{\chi}\left(s^{\prime}\right)\in\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)\right]=\int_{\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)}p\left(\mathbf{x},s|\mathbf{x"},s^{\prime}\right)d\mathbf{x"} (14)

Assuming that p​(𝐱,s|𝐱′,s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) is well-behaved, then in the limit as Δ​A′​(𝐱′,s′)→0,\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)\rightarrow 0, by the mean-value theorem,

∫Δ​A′​(𝐱′,s′)p​(𝐱,s|𝐱​”,s′)​𝑑𝐱​”=p​(𝐱,s|𝐱∗,s′)​Δ​A′​(𝐱′,s′)\int_{\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)}p\left(\mathbf{x},s|\mathbf{x"},s^{\prime}\right)d\mathbf{x"}=p\left(\mathbf{x},s|\mathbf{x^{*}},s^{\prime}\right)\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right) (15)

where 𝐱∗∈Δ​A′​(𝐱′,s′),\mathbf{x^{*}}\in\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), and 𝐱∗→𝐱′\mathbf{x^{*}}\rightarrow\mathbf{x^{\prime}} as Δ​A′​(𝐱′,s′)→0.\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)\rightarrow 0.

Thus, based on the numerical simulation of random walker evolution following launch, approximate the probability on the left side of Eq. (14) as

P​[𝝌​(s′)∈Δ​A′​(𝐱′,s′)]≈n​(𝐱′,s′)N𝐱,sP\left[\boldsymbol{\chi}\left(s^{\prime}\right)\in\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)\right]\approx\frac{n\left(\mathbf{x^{\prime}},s^{\prime}\right)}{N_{\mathbf{x},s}} (16)

Combining Eqs. (14) through (16) then leads to the estimator in (13). As highlighted in [38], and as we find below, estimators like Eq. (13) produce discontinuous estimates of p​(𝐱,s|𝐱′,s′).p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right).

4.1 Heuristic validation of estimator, Eq. (13)

We adapt a central limit argument in [61] to estimate the free-space transition density, po​(𝐱,s|𝐱′,s′),p_{o}\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), for the evolution of a numerically simulated swarm of N𝐱,sN_{\mathbf{x},s} RW’s, launched from (𝐱,s),\left(\mathbf{x},s\right), within an infinite, two-dimensional domain. Of the N𝐱,sN_{\mathbf{x},s} RW’s launched, focus on the MM RW’s that reach or pass through differential area, Δ​A′​(𝐱′,s′),\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), over the backward time interval, n​Δ​s′≤s′<(n+1)​Δ​s′,n\Delta s^{\prime}\leq s^{\prime}<\left(n+1\right)\Delta s^{\prime}, where 𝐱′∈Δ​A′​(𝐱′,s′),\mathbf{x^{\prime}}\in\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), and s′=s+n​Δ​s.s^{\prime}=s+n\Delta s.

Define the space-time vector displacement of the jt​hj^{th} random walker, in the set of MM RW’s captured by Δ​A′​(𝐱′,s′),\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right), as

𝐠𝐧(𝐣)​(s′)=∑k=1n𝚫​𝝌𝒌(𝒋)=𝝌(𝒋)​(s′)−𝐱j=1,2,…​M\mathbf{g_{n}^{\left(j\right)}}\left(s^{\prime}\right)=\sum_{k=1}^{n}\boldsymbol{\Delta}\boldsymbol{\chi_{k}^{\left(j\right)}}=\boldsymbol{\chi^{\left(j\right)}}\left(s^{\prime}\right)-\mathbf{x}\hskip 28.45274ptj=1,2,...M (17)

or, in component form,

gnx(j)\displaystyle g_{n_{x}}^{\left(j\right)} =\displaystyle= Δ​x1(j)+Δ​x2(j)+…+Δ​xn(j)\displaystyle\Delta x_{1}^{\left(j\right)}+\Delta x_{2}^{\left(j\right)}+...+\Delta x_{n}^{\left(j\right)}\hskip 113.81102pt
=\displaystyle= χx(j)​(s′)−x\displaystyle\chi_{x}^{\left(j\right)}\left(s^{\prime}\right)-x\hskip 190.63338pt
gny(j)\displaystyle g_{n_{y}}^{\left(j\right)} =\displaystyle= Δ​y1(j)+Δ​y2(j)+…+Δ​yn(j)\displaystyle\Delta y_{1}^{\left(j\right)}+\Delta y_{2}^{\left(j\right)}+...+\Delta y_{n}^{\left(j\right)}\hskip 113.81102pt
=\displaystyle= χy(j)​(s′)−y\displaystyle\chi_{y}^{\left(j\right)}\left(s^{\prime}\right)-y\hskip 190.63338pt (18)

where

𝚫​𝝌𝒌(𝒋)\displaystyle\boldsymbol{\Delta}\boldsymbol{\chi_{k}^{(j)}} =\displaystyle= 𝝌(𝒋)​((k+1)​Δ​s′)−𝝌(𝒋)​(k​Δ​s′)\displaystyle\boldsymbol{\chi^{(j)}}\left(\left(k+1\right)\Delta s^{\prime}\right)-\boldsymbol{\chi^{(j)}}\left(k\Delta s^{\prime}\right)\hskip 256.0748pt
=\displaystyle= 2​Do​[𝐖(𝐣)​((k+1)​Δ​s′)−𝐖(𝐣)​(k​Δ​s′)]\displaystyle\sqrt{2D_{o}}\left[\mathbf{W^{(j)}}\left(\left(k+1\right)\Delta s^{\prime}\right)-\mathbf{W^{(j)}}\left(k\Delta s^{\prime}\right)\right]\hskip 199.16928pt (19)

is the simulated Brownian displacement of random walker j,j, over the kt​hk^{th} time-interval.

Now, let N𝐱,s→∞,N_{\mathbf{x},s}\rightarrow\infty, so that, likewise, M→∞,M\rightarrow\infty, and define the random variable Y~nx\tilde{Y}_{n_{x}} as the sum of nn independent random variables, Δ​x1,\Delta x_{1}, Δ​x2,…,Δ​xn:\Delta x_{2},...,\Delta x_{n}:

Y~nx=gnxsnx=Δ​x1+Δ​x2+…+Δ​xnsnx\tilde{Y}_{n_{x}}=\frac{g_{n_{x}}}{s_{n_{x}}}=\frac{\Delta x_{1}+\Delta x_{2}+...+\Delta x_{n}}{s_{n_{x}}} (20)

where snx2=⟨Δ​x12⟩+⟨Δ​x22⟩+…+⟨Δ​x1n⟩,s_{n_{x}}^{2}=\langle\Delta x_{1}^{2}\rangle+\langle\Delta x_{2}^{2}\rangle+...+\langle\Delta x_{1}^{n}\rangle, is the sum of variances of the set of nn differential x-displacements. Importantly, recognize that the differential spread, i.e., variance, in backward time x-displacements (over the set of MM random walker trajectories captured on ΔA′(𝐱′,s′))\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)) decreases with increasing backward time: ⟨Δ​x12⟩=−2​Do​Δ​s′,\langle\Delta x_{1}^{2}\rangle=-2D_{o}\Delta s^{\prime}, or, in terms of forward time increments, ⟨Δ​x12⟩=2​Do​Δ​t′,\langle\Delta x_{1}^{2}\rangle=2D_{o}\Delta t^{\prime},

Letting Δ​s′→0,\Delta s^{\prime}\rightarrow 0, so that n→∞,n\rightarrow\infty, then by the central limit theorem [61], the probability density for gnxg_{n_{x}} tends to the Gaussian:

P​(gnx)→14​π​Do​(t−t′)​exp−(x−⟨χx​(t′)⟩)24​Do​(t−t′)P\left(g_{n_{x}}\right)\rightarrow\frac{1}{\sqrt{4\pi D_{o}\left(t-t^{\prime}\right)}}\exp{-\frac{\left(x-\langle\chi_{x}\left(t^{\prime}\right)\rangle\right)^{2}}{4D_{o}\left(t-t^{\prime}\right)}} (21)

where, since ⟨Δ​x1k⟩=2​Do​Δ​t′,\langle\Delta x_{1}^{k}\rangle=2D_{o}\Delta t^{\prime}, k=1,2,…,n,k=1,2,...,n, snx2=2​Do​(t−t′).s_{n_{x}}^{2}=2D_{o}\left(t-t^{\prime}\right). The same argument applied to the y-component of the random relative displacement, gnyg_{n_{y}} leads to:

P​(gny)→14​π​Do​(t−t′)​exp−(y−⟨χy​(t′)⟩)24​Do​(t−t′)P\left(g_{n_{y}}\right)\rightarrow\frac{1}{\sqrt{4\pi D_{o}\left(t-t^{\prime}\right)}}\exp{-\frac{\left(y-\langle\chi_{y}\left(t^{\prime}\right)\rangle\right)^{2}}{4D_{o}\left(t-t^{\prime}\right)}} (22)

where again, t↔s,t\leftrightarrow s, t′↔s′,t^{\prime}\leftrightarrow s^{\prime}, and t≥t′.t\geq t^{\prime}.

Finally, letting Δ​A′​(𝐱′,s′)→0,\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)\rightarrow 0, then ⟨χx​(s′)⟩→x′\langle\chi_{x}\left(s^{\prime}\right)\rangle\rightarrow x^{\prime} and ⟨χy​(s′)⟩→y′.\langle\chi_{y}\left(s^{\prime}\right)\rangle\rightarrow y^{\prime}. Thus, we find that the simple estimator in Eq. (13) for the probability of observing a random walker at (𝐱′,s′),\left(\mathbf{x^{\prime}},s^{\prime}\right), given that it was launched at (𝐱,s),\left(\mathbf{x},s\right), given by the product of the right sides of Eqs. (21) and (22), corresponds to the solution of the diffusive backward Kolmogorov equation, Eq. (6):

p​(𝐱,t|𝐱′,t′)=14​π​Do​(t−t′)​exp−(𝐱−𝐱′)24​Do​(t−t′)p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)=\frac{1}{\sqrt{4\pi D_{o}\left(t-t^{\prime}\right)}}\exp{-\frac{\left(\mathbf{x}-\mathbf{x^{\prime}}\right)^{2}}{4D_{o}\left(t-t^{\prime}\right)}} (23)

Equivalently, and as required [31, 44], this corresponds to the free-space solution to the (pure diffusion) Fokker-Planck equation, governing forward time evolution of p​(𝐱,t|𝐱′,t′),p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right), following delta function release of RW’s, at t=t′,t=t^{\prime}, from 𝐱=𝐱′.\mathbf{x}=\mathbf{x^{\prime}}.

4.2 Green’s function estimate and error in the estimate

Assuming that the adjoint problem governing the propagator, K​(𝐱,s|𝐱′,s′),K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), and the backward Kolmogorov problem governing the transition density have the same structure, then from Eq. (12), for s′−s≥0,s^{\prime}-s\geq 0, the stochastically estimated Green’s function follows immediately from Eq. (13):

G​(𝐱,s|𝐱′,s′)=p​(𝐱,s|𝐱′,s′)=n​(𝐱′,s′)N𝐱,s​Δ​A′​(𝐱′,s′)+O​(Δ​s′)G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=\frac{n\left(\mathbf{x^{\prime}},s^{\prime}\right)}{N_{\mathbf{x},s}\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)}+O\left({\Delta s^{\prime}}\right) (24)

where the Euler-Murayama scheme’s O​(Δ​s′)O\left({\Delta s^{\prime}}\right) weak convergence is again highlighted.

We note that an additional, localized, O​(Δ​s′)O\left(\Delta s^{\prime}\right) error arises in connection with the random walker respawning algorithm introduced below. As shown in Appendix B, this is an implicit discretization error, confined to random locations in the solution domain where individual RW’s are respawned following absorption on Dirichlet boundaries.

5 Absorbing boundaries and random walker respawning

In order to mitigate against loss of absorbed RW’s and preserve early-forward-time solution accuracy, we introduce a random walker respawning algorithm [62], and discuss below and in Appendix B its properties:


a) Initially, at the chosen backward solution instant, s,s, the group of N𝐱,s,N_{\mathbf{x},s}, RW’s to be launched from the solution point, 𝐱,\mathbf{x}, are all assigned a weight of 1, and each random walker weight is stored in an N𝐱,s×1−dimensionalN_{\mathbf{x},s}\times 1-\mathrm{dimensional} array, [W][W] [62].


b) After every backward time-step, Δ​s′,\Delta s^{\prime}, any random walker that reaches a Dirichlet boundary is removed, and the in-domain random walker, χ​(s′;Wm​a​x),\chi\left(s^{\prime};W_{max}\right), having the highest weight, Wm​a​x​(s′),W_{max}\left(s^{\prime}\right), in [W][W] is also removed and replaced by two RW’s, each having weights equal to Wm​a​x​(s′)/2.W_{max}\left(s^{\prime}\right)/2. The replacement pair are placed at the same location occupied by χ​(s′;Wm​a​x).\chi\left(s^{\prime};W_{max}\right).


c) After each set of mm time-steps, the array [W][W] is reordered, in descending order, according to the current set of random walker weights. Thus, high-mass RW’s that experience relatively few absorption events, in particular, those far from Dirichlet boundaries, tend towards Dirichlet boundaries over backward time, suppressing numerical thinning of near-boundary mass distributions.


When using the respawning algorithm, equation (24) must be modified:

G​(𝐱,t|𝐱′,t′)=G​(𝐱,s|𝐱′,s′)=1N𝐱,s​1Δ​A′​W​(𝐱′,s′)+O​(Δ​s′)G\left(\mathbf{x},t|\mathbf{x}^{\prime},t^{\prime}\right)=G\left(\mathbf{x},s|\mathbf{x}^{\prime},s^{\prime}\right)=\frac{1}{N_{\mathbf{x},s}}\frac{1}{\Delta A^{\prime}}W\left(\mathbf{x}^{\prime},s^{\prime}\right)+O\left({\Delta s^{\prime}}\right) (25)

where W​(𝐱′,s′)W(\mathbf{x}^{\prime},s^{\prime}) is the total weight (sum of weights) of the RW’s that reach Δ​A′​(𝐱′,s′)\Delta A^{\prime}(\mathbf{x}^{\prime},s^{\prime}).

5.1 Properties of the respawning algorithm

Appendix B shows that the respawning algorithm has two important properties which ensure the physical consistency of the algorithm, as well as a third that significantly enhances the statistical quality of the Green’s function estimation procedure:


a) The algorithm enforces, at each random respawning location, the continuum diffusion equation governing continuum random walker evolution.


b) Since the algorithm diffusively smooths, on short, single-time-step time-scales, instantaneous, localized mass gradients produced by respawning, the algorithm appears to be insensitive to the number of time steps, m,m, taken between reordering of the random walker weight array, [W].[W]. In this first study, we have not performed tests to ascertain suitable ranges for m;m; in all cases, m=10.m=10.


c) Crucially, the respawning algorithm conserves the number of random walkers, N𝐱,s,N_{\mathbf{x},s}, initially launched. As shown in the Results below, preserving N𝐱,sN_{\mathbf{x},s} is essential to obtaining low-error Green’s function estimates.

6 Numerical experiments

We focus on three test cases, each of increasing complexity: 1) two-dimensional diffusive transport in a square, 2) diffusion in a circle, and 3) unsteady two-dimensional groundwater solute transport, produced by combined advection, dispersion, and chemical depletion in a nonhomogeneous, anisotropic, infinite region. We use the first test to expose and solve two principle challenges: a) loss of random walker’s at Dirichlet boundaries and b) spatial variance in estimated GF’s. The second test examines the performance of the Green’s function estimation procedure in domains having non-planar boundaries. Since tests 1 and 2 allow (preferred, optimal) backward time Green’s function construction, in the next section, we illustrate, in detail, the backward time Green’s function estimation recipe. The third test requires forward time Green’s function estimation; there, we briefly outline the corresponding recipe.

6.1 Application of the backward time Green’s function estimation recipe to pure diffusion test cases 1 and 2

6.1.1 Step 1: Simultaneously derive the adjoint equation and magic rule

Appendix C illustrates trial and error derivation of these two central relations, focusing on a generic advection-diffusion transport problem in a finite domain. Similar trial and error approaches were used to derive the adjoint equations and magic rules used in all three test cases. In tests 1 and 2, the adjoint equation is thus given by:

∂G​(𝐱,s|𝐱′,s′)∂s′+Do​∇′2⁡G​(𝐱,s|𝐱′,s′)=−δ​(𝐱−𝐱′)​δ​(s−s′)\frac{\partial G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)}{\partial s^{\prime}}+D_{o}{\nabla^{\prime}}^{2}G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=-\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\delta\left(s-s^{\prime}\right) (26)

where DoD_{o} is a fixed diffusion coefficient.

Considering the case where a combination of Dirichlet and Neumann boundary conditions are imposed, the magic rule assumes the form:

η​(𝐱,t)=∫0t∫ΩG​(𝐱,t|𝐱′,t′)​f​(𝐱′,t)​𝑑𝐱′​𝑑t′−Do​∫0t∫∂ΩDg​(𝐱′,t′)​∇′G​(𝐱,t|𝐱′,t′)⋅𝐧^′​𝑑𝐱′​𝑑t′++Do​∫0t∫∂ΩNG​(𝐱,t|𝐱′,t′)​∇′h​(𝐱′,t′)⋅𝐧^′​𝑑𝐱′​𝑑t′++∫Ωϕ​(𝐱′)​G​(𝐱,t|𝐱′,0)​𝑑𝐱′.\eta(\mathbf{x},t)=\int_{0}^{t}\int_{\Omega}G(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})f(\mathbf{x^{\prime}},t)d\mathbf{x^{\prime}}dt^{\prime}-D_{o}\int_{0}^{t}\int_{\partial\Omega_{D}}g(\mathbf{x^{\prime}},t^{\prime})\nabla^{\prime}G(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\cdot\mathbf{\hat{n}^{\prime}}d\mathbf{x^{\prime}}dt^{\prime}+\\ +D_{o}\int_{0}^{t}\int_{\partial\Omega_{N}}G(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\nabla^{\prime}h(\mathbf{x^{\prime}},t^{\prime})\cdot\mathbf{\hat{n}^{\prime}}d\mathbf{x^{\prime}}dt^{\prime}+\\ +\int_{\Omega}\phi(\mathbf{x^{\prime}})G(\mathbf{x},t|\mathbf{x^{\prime}},0)d\mathbf{x^{\prime}}. (27)

where f​(𝐱′,t),f(\mathbf{x^{\prime}},t), g​(𝐱′,t′),g(\mathbf{x^{\prime}},t^{\prime}), h​(𝐱′,t′),h(\mathbf{x^{\prime}},t^{\prime}), and ϕ​(𝐱′)\phi(\mathbf{x^{\prime}}) represent, respectively, the areal/volumetric source/sink term, boundary conditions on the Neumann (∂ΩN)(\partial\Omega_{N}) and Dirichlet (∂ΩD)(\partial\Omega_{D}) boundaries, and the initial condition. See Eq. (S1 Appendix C), Appendix C, with the advective (last) term turned off. In addition, 𝐧^′\mathbf{\hat{n}^{\prime}} is the local outward unit normal to the problem domain, Ω′.\Omega^{\prime}.

6.1.2 Steps 2 and 3: Postulate existence of a microscale stochastic transport process, 𝝌​(s),\boldsymbol{\chi}\left(s\right), and write down the associated backward Kolmogorov problem

Since test cases 1 and 2 model pure diffusion, the backward-time stochastic differential equation is obtained by dropping the drift term in Eq. (4):

d​𝝌​(s)=𝐁​(𝝌​(s),s)⋅d​𝐰​(s)d\boldsymbol{\chi}\left(s\right)=\mathbf{B}\left(\boldsymbol{\chi}\left(s\right),s\right)\cdot d\mathbf{w}\left(s\right) (28)

where, from Eq. (5), ⟨𝐁⟩=2​Do​𝐈.\langle\mathbf{B}\rangle=\sqrt{2D_{o}}\mathbf{I}. Thus, the corresponding backward Kolmogorov equation, determined by dropping the drift term from Eq. (6), is given by

∂p​(𝐱,s|𝐱′,s′)∂s′+Do​∇′2⁡p​(𝐱,s|𝐱′,s′)=0\frac{\partial p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)}{\partial s^{\prime}}+D_{o}{\nabla^{\prime}}^{2}p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=0 (29)

where p​(𝐱,s|𝐱′,s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) satisfies the terminal condition in Eq. (7). Finally, in tests 1 and 2, Dirichlet boundary conditions, as stated in Eq. (8), are imposed at all points on both domain boundaries.

6.1.3 Step 4: Determine the applicability of the stochastic Green’s function estimation procedure

Comparing Eqs. (26) and (29), due to the non-zero right side in (26), it is clear that these do not have the same mathematical structure. However, the propagator, K​(𝐱,s|𝐱′,s′),K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), which is connected to the Green’s function, G​(𝐱,s|𝐱′,s′)G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) via Eq.(12), does satisfy the equation L∗​K​(𝐱,s|𝐱′,s′)=0,L^{*}K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=0, where L∗=∂/∂s′+Do​∇′2=∂/∂t′+Do​∇′2.L^{*}=\partial/\partial s^{\prime}+D_{o}{\nabla^{\prime}}^{2}=\partial/\partial t^{\prime}+D_{o}{\nabla^{\prime}}^{2}. Thus, using L∗L^{*} to operate on H​(τ)​K​(𝐱,s|𝐱′,s′),H\left(\tau\right)K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), where τ=t−t′,\tau=t-t^{\prime}, using ∂H/∂t′=−δ​(τ)=−δ​(−τ)=−δ​(s′−s)\partial H/\partial t^{\prime}=-\delta\left(\tau\right)=-\delta\left(-\tau\right)=-\delta\left(s^{\prime}-s\right) as well as the terminal condition, p​(𝐱,s|𝐱′,s′)→δ​(𝐱−𝐱′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)\rightarrow\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right) as s′→t,s^{\prime}\rightarrow t, we obtain Eq. (26).

Thus, since the problems governing the propagator, K​(𝐱,s|𝐱′,s′)K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) and transition density, p​(𝐱,s|𝐱′,s′)p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) have identical mathematical structure, including matching BC’s and IC’s, we can use a density estimation procedure to first estimate K​(𝐱,s|𝐱′,s′),K\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right), and then immediately obtain G​(𝐱,s|𝐱′,s′)G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right) via Eq. (12).

6.1.4 Step 5: Choose an integration scheme for the problem SDE and a density estimation method

As noted, in all three tests, we use the Euler-Murayama scheme [44, 46, 47] to integrate Eq. (28). Thus, discretize the backward time axis as: s′i+1=s′i+Δ​s′,{s^{\prime}}_{i+1}={s^{\prime}}_{i}+\Delta s^{\prime}, i=0,1,2,…​M,i=0,1,2,...\mathrm{M}, where s′o=s,{s^{\prime}}_{o}=s, the chosen backward solution time, corresponds to a chosen forward solution time, t,t, s′s^{\prime} is the chosen random walker interrogation time, and M​Δ​s′=s′−s.\mathrm{M}\Delta s^{\prime}=s^{\prime}-s. Defining Δ​𝐰​(s′i)=𝐰​(s′i+1)−𝐰​(s′i),\Delta\mathbf{w}\left({s^{\prime}}_{i}\right)=\mathbf{w}\left({s^{\prime}}_{i+1}\right)-\mathbf{w}\left({s^{\prime}}_{i}\right), as the psuedo-random, numerical approximation to the Weiner differential, d​𝐰​(s′i)d\mathbf{w}\left({s^{\prime}}_{i}\right) in Eq. (28), the EM scheme is given by [44, 47]:

𝐲​(s′i+1)=𝐲​(s′i)+2​Do​Δ​𝐰​(s′i)\mathbf{y}\left({s^{\prime}}_{i+1}\right)=\mathbf{y}\left({s^{\prime}}_{i}\right)+\sqrt{2D_{o}}\Delta\mathbf{w}\left({s^{\prime}}_{i}\right) (30)

where the numerical approximation to the actual stochastic jump, d​𝝌​(s′i)=𝝌​(s′i+1)−𝝌​(s′i),d\boldsymbol{\chi}\left({s^{\prime}}_{i}\right)=\boldsymbol{\chi}\left({s^{\prime}}_{i+1}\right)-\boldsymbol{\chi}\left({s^{\prime}}_{i}\right), in (28) is denoted as 𝐲​(s′i+1)−𝐲​(s′i).\mathbf{y}\left({s^{\prime}}_{i+1}\right)-\mathbf{y}\left({s^{\prime}}_{i}\right). Details on simulation of Δ​𝐰​(s′i)\Delta\mathbf{w}\left({s^{\prime}}_{i}\right) can be found, for example, in [44, 47].

6.2 Explanatory notes

A) The stochastic estimation technique requires specification of three numerical parameters: i) the backward (forward) time-step size, Δ​s′\Delta s^{\prime} (Δ​t),\left(\Delta t\right), ii) the interrogation area, Δ​A′​(𝐱′,s′)\Delta A^{\prime}\left(\mathbf{x}^{\prime},s^{\prime}\right) (Δ​A′​(𝐱′,s′)),\left(\Delta A^{\prime}\left(\mathbf{x}^{\prime},s^{\prime}\right)\right), and iii) the number of RW’s, N𝐱,sN_{\mathbf{x},s} (N𝐱′,t′),\left(N_{\mathbf{x^{\prime}},t^{\prime}}\right), launched from the chosen Green’s function solution point, (𝐱,s)\left(\mathbf{x},s\right) ((𝐱′,t′)).\left(\left(\mathbf{x^{\prime}},t^{\prime}\right)\right). Appendix A proposes procedures for specifying these parameters.


B) In all three tests, lengths, (x,y;x′,y′),\left(x,y;x^{\prime},y^{\prime}\right), forward and backward times, (t,t′;s,s′),\left(t,t^{\prime};s,s^{\prime}\right), and diffusion coefficients, Do,D_{o}, are expressed in non-dimensional form: i) (x,y;x′,y′)=(x~,y~;x~′,y~′)​L~−1,\left(x,y;x^{\prime},y^{\prime}\right)=\left(\tilde{x},\tilde{y};\tilde{x}^{\prime},\tilde{y}^{\prime}\right)\tilde{L}^{-1}, ii)(t,t′;s,s′)=(t~,t~′;s~,s~′)​t~s−1,\left(t,t^{\prime};s,s^{\prime}\right)=\left(\tilde{t},\tilde{t}^{\prime};\tilde{s},\tilde{s}^{\prime}\right)\tilde{t}_{s}^{-1}, , iii) Do=D~o​t~s​L~−2,D_{o}=\tilde{D}_{o}\tilde{t}_{s}\tilde{L}^{-2}, where L~\tilde{L} and t~s\tilde{t}_{s} are dimensional, problem-specific length and time-scales. In the first two tests, D~o\tilde{D}_{o} is a dimensional diffusion coefficient, and in the third, a dimensional dispersion coefficient. Throughout the rest of this paper, all functions, variables, and parameters without a tilde are dimensionless.

In all three tests, we choose a dimensionless solution time, t=10,t=10, and dimensionless length, L=1.L=1. In order to gain a sense of associated dimensional scales, consider diffusion of small molecular and ionic species in water, where D~o=O​(10−9​m2​s−1);\tilde{D}_{o}=\mathrm{O}\left(10^{-9}\ \mathrm{m^{2}}\mathrm{s^{-1}}\right); since diffusive spread of a specie introduced at a point grows as L~∼D~o​t~,\tilde{L}\sim\sqrt{\tilde{D}_{o}\tilde{t}}, then choosing time scales, t~s,\tilde{t}_{s}, equal to, e.g., 10−3​s10^{-3}\ \mathrm{s} and 107​s,10^{7}\ \mathrm{s}, requires computational domains having dimensions on the order of L~∼10​D~o​t~s∼3​(10−6)​m\tilde{L}\sim\sqrt{10\tilde{D}_{o}\tilde{t}_{s}}\sim 3\left(10^{-6}\right)\ \mathrm{m} and ∼3​(10−1)​m,\sim 3\left(10^{-1}\right)\ \mathrm{m}, respectively. Similar considerations can be used in the groundwater transport problem in test 3.


C) In test 1, we discretize the square solution domain, having dimensionless length, L=1,L=1, into 10410^{4} (square) interrogation areas. In test 2, a circular solution domain of dimensionless diameter, D=1,D=1, is inscribed and centered within the same unit square. Again, the square is discretized into 10410^{4} square interrogation areas. Thus, in both tests, Δ​x′=Δ​y′=10−2,\Delta x^{\prime}=\Delta y^{\prime}=10^{-2}, and Δ​A𝐱𝐣′,s′=Δ​x′​Δ​y′=10−4,\Delta A_{\mathbf{x_{j}^{{}^{\prime}}},s^{\prime}}=\Delta x^{\prime}\Delta y^{\prime}=10^{-4}, where 𝐱𝐣′\mathbf{x_{j}^{{}^{\prime}}} denotes the centroid of the jt​hj^{th} interrogation area. While test 3 considers ground water transport in an infinite domain, for finite solution time, t=10,t=10, and under the groundwater flow conditions considered, we can again define a unit square solution domain, again having 10410^{4} interrogation areas; see Section 7.

In tests 1 and 2, the dimensionless diffusion coefficient is specified as Do=0.05,D_{o}=0.05, while in test 3, the dimensionless dispersion coefficient is also set as Do=0.05.D_{o}=0.05. As described in test 1, after some trial and error, the backward time step is chosen as Δ​s′=10−4;\Delta s^{\prime}=10^{-4}; likewise, in test 3, Δ​t=10−4.\Delta t=10^{-4}. Finally, in tests 1 and 2, the number of RW’s launched from a chosen (forward-time) response point, (𝐱,s=0),\left(\mathbf{x},s=0\right), is N𝐱,s=0=106;N_{\mathbf{x},s=0}=10^{6}; in test 3, N𝐱′,t′=0=106,N_{\mathbf{x^{\prime}},t^{\prime}=0}=10^{6}, corresponding to the number of RW’s launched from a (forward-time) impulse point, (𝐱′,t′=0).\left(\mathbf{x^{\prime}},t^{\prime}=0\right).

Based on these choices, and as detailed in Appendix A, the estimated relative spatial variation in captured RW’s, σn/N,\sigma_{n}/N, over the chosen forward solution interval, 0<t≤10,0<t\leq 10, ranges between

∼1.6×10−5<σnN<∼1.6\sim 1.6\times 10^{-5}<\frac{\sigma_{n}}{N}<\hskip 5.69046pt\sim 1.6

where the maximum variation occurs at late forward (early backward) time, t=Δ​t=10−4t=\Delta t=10^{-4} (s′=Δ​s′=10−4).\left(s^{\prime}=\Delta s^{\prime}=10^{-4}\right). From Eq. (48) Appendix A, large, early-backward-time variations can be suppressed by choosing sufficiently large N=N𝐱,s=0N=N_{\mathbf{x},s=0} (or for forward time random walker integration, large enough N𝐱′,t′=0).N_{\mathbf{x^{\prime}},t^{\prime}=0}). Our choice for N𝐱,s=0N_{\mathbf{x},s=0} (N𝐱′,t′=0)\left(N_{\mathbf{x^{\prime}},t^{\prime}=0}\right) reflects a compromise between computational cost and variance reduction. In addition, our chosen Δ​x′=Δ​y′=Δ​x=Δ​y=10−2\Delta x^{\prime}=\Delta y^{\prime}=\Delta x=\Delta y=10^{-2} are approximately five times the (average) incremental diffusion distance, Do​Δ​s′≈2×10−3;\sqrt{D_{o}\Delta s^{\prime}}\approx 2\times 10^{-3}; thus, and consistent with the results below, we appear to be obtaining statistically reasonable estimates of n​(𝐱′,s′)n\left(\mathbf{x^{\prime}},s^{\prime}\right) [n​(𝐱,t)].[n\left(\mathbf{x},t\right)].

6.3 Test Case 1, Part A: Deleterious depletion of random walkers at Dirichlet boundaries, remediation via random walker respawning, and error due to overlarge time step size

We use the first test case to emphasize the experimental nature of constructing Green’s function estimators, briefly highlighting the challenges we encountered, as well as our solutions. In test 1, the solution (response) point is chosen as (𝐱,s)=(x=0.5,y=0.5,s=0).\left(\mathbf{x},s\right)=\left(x=0.5,y=0.5,s=0\right). In all figures, we show estimated and actual GF’s at various forward time instants, t′,t^{\prime}, where 0≤t′≤T=10.0\leq t^{\prime}\leq T=10. For reference, Fig. 1 shows the exact Green’s function.

The detrimental effect of random walker absorption at Dirichlet boundaries is illustrated in Fig. 2, where depletion of RW’s becomes increasingly apparent as backward time integration of random walker trajectories proceeds. Comparing Figs. 1 and 2, we observe that the maximum local error,

em​a​x=m​a​x​|Ge​s​t​i​m​a​t​e​d−Ge​x​a​c​t|/Ge​x​a​c​te_{max}=max\ |G_{estimated}-G_{exact}|/G_{exact} (31)

increases from 4.1%4.1\% at forward time, t′=9.9,t^{\prime}=9.9, to 5200%5200\% at t′=1.t^{\prime}=1. Moreover, while relatively late forward time estimates, obtained at t′=9.9,t^{\prime}=9.9, t′=9.5,t^{\prime}=9.5, and t′=9.0,t^{\prime}=9.0, shown respectively in Figs. 2a, 2b, and 2c, remain qualitatively consistent with the actual Green’s function in Fig. 1a-c, the early forward time estimate at t′=1,t^{\prime}=1, Fig. 2d, is spatially amorphous, bearing no resemblance to the actual Green’s function depicted in Fig. 1d. While the results in Fig. 2 were obtained using a time step, Δ​s′=10−1,\Delta s^{\prime}=10^{-1}, which is approximately 2 orders of magnitude larger than the recommended Δ​so′≈Δ​A′/Do=2​(10−3)\Delta s_{o}^{\prime}\approx\Delta A^{\prime}/D_{o}=2\left(10^{-3}\right), see Appendix A, similar qualitative results, not shown, are observed for 10−4≤Δ​s≤10−1.10^{-4}\leq\Delta s\leq 10^{-1}.

Refer to caption
Figure 1: Test 1: Exact Green’s function, obtained using an analytical solution [2], shown at forward time instants: (a) t′=9.9,t^{\prime}=9.9, (b) t′=9.5,t^{\prime}=9.5, (c) t′=9.0,t^{\prime}=9.0, and (d) t′=1.0.t^{\prime}=1.0. The total solution time is t=T=10.
Refer to caption
Figure 2: Test 1, Part A: Backward time estimated Green’s functions, shown at forward times (a) t′=9.9,t^{\prime}=9.9, (b) t′=9.5,t^{\prime}=9.5, (c) t′=9.0,t^{\prime}=9.0, and (d) t′=1.0.t^{\prime}=1.0. While the backward time-step, Δ​s′=10−1,\Delta s^{\prime}=10^{-1}, is approximately 2 orders of magnitude larger than the recommended Δ​so′≈Δ​A′/Do=2​(10−3)\Delta s_{o}^{\prime}\approx\Delta A^{\prime}/D_{o}=2\left(10^{-3}\right), see Appendix A, the major source of error arises from depletion of RW’s at absorbing Dirichlet boundaries. Maximum relative errors in each estimate are: (a) em​a​x=4.1%,e_{max}=4.1\%, (b) em​a​x=10%,e_{max}=10\%, (c) em​a​x=21%,e_{max}=21\%, and (d) em​a​x=5200%.e_{max}=5200\%.

Introduction of random walker respawning significantly improves both the accuracy and the qualitative consistency of Green’s function estimates, as shown in Figs. 3a and 3b. These depict, respectively, Ge​s​t​(𝐱,t=10|𝐱′,t′=1),G_{est}(\mathbf{x},t=10|\mathbf{x^{\prime}},t^{\prime}=1), 𝐱=(0.5,0.5),\mathbf{x}=\left(0.5,0.5\right), at the early forward time, t′=1t^{\prime}=1 (corresponding to s′=9).s^{\prime}=9). The two estimates shown in Figs. 3a and 3b use backward time steps, Δ​s′=10−1\Delta s^{\prime}=10^{-1} and Δ​s′=10−4,\Delta s^{\prime}=10^{-4}, respectively, which bracket the recommended Δ​s′=Δ​so′≈Δ​A′/Do=2​(10−3).\Delta s^{\prime}=\Delta s_{o}^{\prime}\approx\Delta A^{\prime}/D_{o}=2\left(10^{-3}\right). Again, similar qualitative results are observed for 10−4≤Δ​s≤10−1.10^{-4}\leq\Delta s\leq 10^{-1}.

In order to improve the accuracy of the estimate shown in Fig. 3a, the too-large time step, Δ​s′=10−1=50×Δ​A′/Do,\Delta s^{\prime}=10^{-1}=50\times\Delta A^{\prime}/D_{o}, is reduced to 10−4=0.05×Δ​A′/Do,10^{-4}=0.05\times\Delta A^{\prime}/D_{o}, reducing the maximum relative error, again observed at t′=1,t^{\prime}=1, from 570 % to 25 %. Although not shown, similar results are observed using Δ​s′=Δ​A′/Do.\Delta s^{\prime}=\Delta A^{\prime}/D_{o}.

Refer to caption
Figure 3: Test 1, Part A: Estimated Green’s function at t′=1.0t^{\prime}=1.0, s′=9,s^{\prime}=9,test 1, constructed using the respawning algorithm; compare with the exact solution in Fig. 1d and the estimate obtained in Fig. 2d, without respawning. The backward time step used in (a) is fifty times larger than the recommended time step, Δ​so′=Δ​A′/Do,\Delta s_{o}^{\prime}=\Delta A^{\prime}/D_{o}, see Appendix A, and in (b), twenty times smaller than Δ​so′.\Delta s_{o}^{\prime}. Maximum errors, relative to the exact solution in Fig. 1d, are: (a) em​a​x=570%,e_{max}=570\%, and (b) em​a​x=25%.e_{max}=25\%.

6.4 Test 1, Part B: Area averaging method for improving Green’s function estimation accuracy and smoothness

While reducing Δ​s′\Delta s^{\prime} to magnitudes on the order of Δ​so′\Delta s_{o}^{{}^{\prime}} or less significantly improves Green’s function estimation accuracy, relative error, at long backward integration times, here, s′=9,s^{\prime}=9, remains unacceptably large. Moreover, as shown in Fig. 3, plate b), long (backward) time estimated GF’s exihibit spatial graininess, an artifact of using the naive estimator in Eq. (25).

In order to address both issues, we introduce a simple area-averaging technique in which the value of m​(𝐱′,s′)m\left(\mathbf{x}^{\prime},s^{\prime}\right) in (25), obtained over an area, Δ​A′,\Delta A^{\prime}, encompassing grid point 𝐱′=(x′,y′),\mathbf{x}^{\prime}=\left(x^{\prime},y^{\prime}\right), is replaced with the average of all mm’s observed at all grid points lying within a square (x±a,y±a)(x\pm a,y\pm a), centered on (x′,y′).\left(x^{\prime},y^{\prime}\right). Here, the smoothing window size, aa, is the smaller of two values: the distance to the nearest boundary, or a maximum window size, am​a​x.a_{max}. At any given backward time instant, s′,s^{\prime}, am​a​xa_{max} is determined by iteratively altering aa in order to minimize the total root mean deviation,

σG¯=1Mi​n​c​∑i=1Mi​n​c′[G​(𝐱i′,s′)−G¯​(𝐱i′,s′)]2\overline{\sigma_{G}}=\frac{1}{M_{inc}}\sum_{i=1}^{M_{inc}^{\prime}}\sqrt{\left[G(\mathbf{x}_{i}^{\prime},s^{\prime})-\overline{G}(\mathbf{x}_{i}^{\prime},s^{\prime})\right]^{2}} (32)

between individual estimated and averaged Greens functions. Here, Mi​n​cM_{inc} is the number of grid points included in the sum, G​(𝐱i′,s′)G(\mathbf{x}_{i}^{\prime},s^{\prime}) is the estimated Green’s function at grid point ii, as given by (25), and G¯​(𝐱i′,s′)\overline{G}(\mathbf{x}_{i}^{\prime},s^{\prime}) is the averaged Green’s function. As noted, we use known, exact GF’s in place of G¯​(𝐱i′,s′).\overline{G}(\mathbf{x}_{i}^{\prime},s^{\prime}). The prime on the sum specifies that only those grid points where the summand is not equal to zero are included. We include this additional filter since, in many cases, there are a significant number of points in both the analytical and numerical Green’s function where both functions are zero. These flat areas, which can be seen in many late forward time GF’s, can skew the total deviation to spuriously low magnitudes.

Refer to caption
Figure 4: Test 1, Part B: Estimated Green’s functions, constructed using random walker respawning and area averaging. The backward time step, Δ​s′=10−4=0.05×Δ​A′/Do.\Delta s^{\prime}=10^{-4}=0.05\times\Delta A^{\prime}/D_{o}. Forward time instants, t′,t^{\prime}, maximum relative errors, em​a​x,e_{max}, and maximum smoothing window sizes, nm​a​x=ls′/(2​Δ​x′),n_{max}=l^{\prime}_{s}/\left(2\Delta x^{\prime}\right), where ls′l^{\prime}_{s} is the dimensionless length of the square smoothing window, are, respectively: (a) t′=9.9,t^{\prime}=9.9, em​a​x=0.31%,e_{max}=0.31\%, nm​a​x=1,n_{max}=1, (b) t′=9.5,t^{\prime}=9.5, em​a​x=0.8%,e_{max}=0.8\%, nm​a​x=3,n_{max}=3, (c) t′=9.0,t^{\prime}=9.0, em​a​x=1.4%,e_{max}=1.4\%, nm​a​x=4,n_{max}=4, and (d) t′=1.0,t^{\prime}=1.0, em​a​x=2.9%,e_{max}=2.9\%, nm​a​x=13.n_{max}=13.

Finally, we observe that combining random walker respawning and area averaging provides relatively accurate estimated GF’s, as shown in Fig.4. Comparing estimated and exact GF’s, in Figs. 4 and 1, respectively, it is clear that, even at large backward estimation times, s′=9,s^{\prime}=9, Fig. 4d, estimation accuracy and smoothness remain reasonable. At small backward estimation times, s′=0.1,s^{\prime}=0.1, s′=0.5,s^{\prime}=0.5, and s′=1.0,s^{\prime}=1.0, shown respectively in Figs. 4a - c, maximum local relative errors remain less than 1.5%,1.5\%, and increase only to 2.9%2.9\% at s′=9,s^{\prime}=9, in Fig. 4d. Note too that Green’s function magnitudes vary by more than five orders of magnitude over the chosen forward time interval, 1≤t′≤9.9.1\leq t^{\prime}\leq 9.9. Although not tested here, but as indicated by the O​(Δ​s′)\mathrm{O}\left(\Delta s^{\prime}\right) weak convergence of the Milstein-Euler integration scheme, we believe that improved accuracy can be obtained, at significantly higher computational cost, by further reductions in Δ​s′.\Delta s^{\prime}.

Two final remarks are made concerning area averaging.

a) Beyond improved estimation accuracy and Green’s function smoothing, area averaging reduces the number of random walkers required to achieve a given level of accuracy, significantly improving computational efficiency. For example, experiments show that achieving ∼1%\sim 1\% late (backward) time accuracy, without area averaging, requires on the order of twenty million RW’s.

b) For poor Green’s function estimates, like those shown in Figure 3, plate a), no optimum window size can be found. There, the smoothing window minimizing the total deviation spans the entire domain, producing a constant function equal to the average of the values of all grid points.


6.5 Test 2: Green’s function estimates near curved boundaries

Many evolution problems take place in regions having complicated, non-planar boundaries. As a preliminary test of the stochastic Green’s function estimation procedure in such problems, we perform a second set of numerical experiments in which GF’s are constructed within a two-dimensional circular domain. Referring to Figures 6 and 5, a circular domain (boundary not shown), having radius 0.5 and center at (x,y)=(0.5,0.5)\left(x,y\right)=\left(0.5,0.5\right), is inscribed and centered within the unit square shown. Here, any random walker that reaches or passes through the circle, (x−0.5)2+(y−0.5)2=0.25,\left(x-0.5\right)^{2}+\left(y-0.5\right)^{2}=0.25, is absorbed, triggering random walker respawning at a random location within the circle. For simplicity, the square is discretized into the same uniform set of square 10410^{4} interrogation areas, Δ​A′​(𝐱′,s′)=Δ​x′​Δ​y′,\Delta A^{\prime}\left(\mathbf{x^{\prime}},s^{\prime}\right)=\Delta x^{\prime}\Delta y^{\prime}, used in test 1, introducing an O​(Δ​x)\mathrm{O}\left(\Delta x\right) discretization error in boundary-adjacent Green’s function estimates. As in Example 1, there is no flow (𝐛=0)(\mathbf{b}=0), transport takes place by isotropic diffusion, with Do=0.05;D_{o}=0.05; we use a backward-time step, Δ​s′=10−4=0.05×Δ​A′/Do,\Delta s^{\prime}=10^{-4}=0.05\times\Delta A^{\prime}/D_{o}, and as in test 1, N𝐱,s=106N_{\mathbf{x},s}=10^{6} RW’s are launched, in this case, from a right-of-center solution point (𝐱,t)=(x=0.75,y=0.5,t=T=10).\left(\mathbf{x},t\right)=\left(x=0.75,y=0.5,t=T=10\right).

Refer to caption
Figure 5: Test 2: Exact Green’s function in a circle of radius 0.5,0.5, inscribed and centered within the unit square shown. The chosen solution (response) point is (x=0.75,y=0.5,t=10),\left(x=0.75,y=0.5,t=10\right), and the exact Green’s function is shown at forward time instants: (a) t′=9.9,t^{\prime}=9.9, (b) t′=9.5,t^{\prime}=9.5, (c) t′=9.0,t^{\prime}=9.0, and (d) t′=1.0.t^{\prime}=1.0.
Refer to caption
Figure 6: Test 2: Estimated, area-averaged Green’s functions in a circle of radius 0.5,0.5, inscribed and centered within the unit square shown. Random walker respawning is used and smoothing windows are circular. The unit square is discretized into 100×100100\times 100 random walker interrogation areas, Δ​A′=Δ​x′×Δ​y′=10−4,\Delta A^{\prime}=\Delta x^{\prime}\times\Delta y^{\prime}=10^{-4}, introducing an O​(Δ​x′)\mathrm{O}\left(\Delta x^{\prime}\right) discretization error in Green’s function estimates. The solution point is at (x=0.75,y=0.5).\left(x=0.75,y=0.5\right). Forward time instants, t′,t^{\prime}, maximum relative errors, em​a​x,e_{max}, and maximum smoothing window sizes, nm​a​x=a′/(Δ​x′),n_{max}=a^{\prime}/\left(\Delta x^{\prime}\right), where a′a^{\prime} is the radius of the circular smoothing window, are, respectively: (a) t′=9.9,t^{\prime}=9.9, em​a​x=1.57%,e_{max}=1.57\%, nm​a​x=2,n_{max}=2, (b) t′=9.5,t^{\prime}=9.5, em​a​x=0.86%,e_{max}=0.86\%, nm​a​x=4,n_{max}=4, (c) t′=9.0,t^{\prime}=9.0, em​a​x=0.89%,e_{max}=0.89\%, nm​a​x=6,n_{max}=6, and (d) t′=1.0,t^{\prime}=1.0, em​a​x=0.70%,e_{max}=0.70\%, nm​a​x=17.n_{max}=17.

In order to adapt the smoothing procedure to a circular domain, we define the smoothing window as a circle of radius aa centered on the grid point of interest. As before, aa is allowed to vary, being the smaller of two values: the distance to the boundary and the maximum window size. By this approach, the smoothing window is always contained within the circular domain. Although not explored here, beyond circular and square/rectangular smoothing window shapes - the former useful for domains having curved boundaries - polygonal shapes may also prove useful.

Exact and estimated GF’s are shown respectively in Figs. 5 and 6. In this case, and in contrast to test 1, the largest estimation errors are observed at the earliest backward time instant, Δ​s′=0.1.\Delta s^{\prime}=0.1. Comparing patterns in maximum relative error evolution in test 1 and 2, Figs. 4 and 6, as expected, in test 1, em​a​xe_{max} progressively increases with increasing s′,s^{\prime}, at least for s′=0.1,s^{\prime}=0.1, s′=1.0,s^{\prime}=1.0, and s′=9.0.s^{\prime}=9.0. By contrast, in test 2, em​a​xe_{max} progressively decreases at all four probed time instants. These results indicate the following: I) The fact that em​a​xe_{max}’s at s′=1.0s^{\prime}=1.0 and s′=9,s^{\prime}=9, in test 1, Fig. 4, plates c) and d), are significantly larger than those in test 2, Fig. 5, plates c) and d), suggests that use of circular smoothing windows in test 2, as opposed to square windows in test 1, reduces estimation error, even under conditions where near-boundary Green’s function estimates are subject to an O​(Δ​x′)\mathrm{O}\left(\Delta x^{\prime}\right) discretization error. II) In light of point a), the relatively large error at s′=0.1s^{\prime}=0.1 in test 2 appears to reflect our choice of the asymmetrically placed solution point, (x=0.75,y=0.5).(x=0.75,y=0.5).

Crucially, over the forward time span, 0.1≤t′≤9,0.1\leq t^{\prime}\leq 9, high-fidelity Green’s function estimates are obtained, here, over approximately six orders of magnitude in local Green’s function magnitudes.

7 Test 3: Forward time Green’s function estimation - advection-dispersion-reaction problems

The third example addresses two questions:

a) Can Green’s functions be stochastically estimated using forward time random walker integration?

b) How well does the stochastic Green’s function estimation procedure perform for advection-dispersion (or diffusion)-reaction problems?

As noted, Green’s function estimation using forward time random walker integration corresponds to estimating the response at a set of QQ space-time points, [𝐱𝟏,t;𝐱𝟐,t;…𝐱𝐐,t,],\left[\mathbf{x_{1}},t;\mathbf{x_{2}},t;...\mathbf{x_{Q}},t,\right], produced by a (delta function) impulse at a single point, (𝐱′,t′):\left(\mathbf{{x}^{\prime}},t^{\prime}\right): G​(𝐱,t|𝐱′,t′).G\left(\mathbf{x},t|\mathbf{{x}^{\prime}},t^{\prime}\right). By contrast, Green’s function estimated via backward time integration provides the response at the point, (𝐱,t),\left(\mathbf{x},t\right), as produced by a (discrete) set of QQ impulse points, [𝐱𝟏′,t′;𝐱𝟐′,t′;…​𝐱𝐐′,t′].\left[\mathbf{x_{1}^{\prime}},t^{\prime};\mathbf{x_{2}^{\prime}},t^{\prime};...\mathbf{x_{Q}^{\prime}},t^{\prime}\right]. Thus, for physical transport problems where property η\eta is not required over an entire problem domain, backward integration offers distinct cost advantages.

In order to obtain an analytical Green’s function, we take advantage of an exact solution obtained by [20] for a ground water advection-dispersion-reaction problem, governed by

∇T⊗∇:𝐃​η−∇⋅(𝐯​η)−γ​η−∂η∂t=−f​(𝐱,t),\nabla^{T}\otimes\nabla:\mathbf{D}\hskip 1.42271pt\eta-\nabla\cdot\left(\mathbf{v}\eta\right)-\gamma\eta-\frac{\partial\eta}{\partial t}=-f(\mathbf{x},t), (33)

corresponding to a slightly generalized version of Eq. (1), and designed to capture a wide range of ground water transport conditions [20], including, e.g., variable density flow in near-coastal aquifers [63].

Importantly, using the substitutions preceding Eq. (2), the sink term, γ​η,\gamma\eta, in Eq. (33) can be suppressed, so that the transformed version of Eq. (33) has the generic form of a Fokker-Planck (FP) equation, with p​(𝐱,t|𝐱′,t′)p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right) replacing η​(𝐱,t):\eta\left(\mathbf{x},t\right):

∇T⊗∇:𝐃​p−∇⋅(𝐯​p)−∂p∂t=0,\nabla^{T}\otimes\nabla:\mathbf{D}\hskip 1.42271ptp-\nabla\cdot\left(\mathbf{v}p\right)-\frac{\partial p}{\partial t}=0, (34)

and where the delta function source term treated in [20],

f​(𝐱,t)=δ​(𝐱−𝐱′)​δ​(t−t′)f\left(\mathbf{x},t\right)=\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\delta\left(t-t^{\prime}\right) (35)

corresponds to the same-time condition,

p​(𝐱,t|𝐱′,t′)→δ​(𝐱−𝐱′)t→t′p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\rightarrow\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\hskip 19.91684ptt\rightarrow t^{\prime} (36)

imposed on the FP equation. Thus, the analytical solution for η​(𝐱,t)\eta\left(\mathbf{x},t\right) obtained in [20] can be interpreted as a solution to the Fokker-Planck equation for p​(𝐱,t|𝐱′,t′).p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right).

We show in Appendix D that the transport equation solved by the estimated forward time Green’s function is given by

𝐃:∇T⊗∇η−𝐯𝐚⋅∇η−∂η∂t=−f​(𝐱,t),\mathbf{D}\mathbf{:}\nabla^{T}\otimes\nabla\eta-\mathbf{v_{a}}\cdot\nabla\eta-\frac{\partial\eta}{\partial t}=-f(\mathbf{x},t), (37)

where the actual advective flow field, 𝐯𝐚​(𝐱,t),\mathbf{v_{a}}\left(\mathbf{x},t\right), in (37) corresponds to a reversed flow field in Eq. (34): 𝐯​(𝐱,t)=−𝐯𝐚​(𝐱,t).\mathbf{v}\left(\mathbf{x},t\right)=-\mathbf{v_{a}}\left(\mathbf{x},t\right).

7.1 Recipe for forward time Green’s function estimation, in brief

Since the Fokker-Planck equation governs the transition density for observing members of a random walker swarm at response point (𝐱,t),\left(\mathbf{x},t\right), launched in the forward time direction from impulse point (𝐱′,t′),\left(\mathbf{x^{\prime}},t^{\prime}\right), computationally, we can use the same recipe given above for estimating p​(𝐱,t|𝐱′,t′)p\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right) via backward random walker integration. Here, we highlight necessary modifications.

The stochastic differential equation governing the forward time motion of individual RW’s is given by

d​𝝌​(t)=𝐯​(𝝌​(t),t)​d​t+𝐁​(𝝌​(t),t)⋅d​𝐰​(t)d\boldsymbol{\chi}\left(t\right)=\mathbf{v}\left(\boldsymbol{\chi}\left(t\right),t\right)dt+\mathbf{B}\left(\boldsymbol{\chi}\left(t\right),t\right)\cdot d\mathbf{w}\left(t\right) (38)

where, again, 𝐰​(t)\mathbf{w}\left(t\right) is an n-dimensional Weiner process, and 𝐯​(𝐱,t)\mathbf{v}\left(\mathbf{x},t\right) is the velocity field.

In forward time integration, a swarm of N𝐱′,t′N_{\mathbf{x^{\prime}},t^{\prime}} is launched from impulse point (𝐱′,t′),\left(\mathbf{x^{\prime}},t^{\prime}\right), where t≥t′,t\geq t^{\prime}, with integration stopped at a chosen final time, t=T.t=T. The spatial solution domain is discretized into NΔ​AN_{\Delta A} interrogation areas, Δ​A​(𝐱,t)=Δ​x​Δ​y,\Delta A\left(\mathbf{x},t\right)=\Delta x\Delta y, and the number of RW’s reaching each Δ​A​(𝐱,t),\Delta A\left(\mathbf{x},t\right), n​(𝐱,t),n\left(\mathbf{x},t\right), at a set of chosen times, t1,t2,…,tq=T,t_{1},t_{2},...,t_{q}=T, is tabulated. As noted immediately below, random walker respawning is not needed in Test 3. Thus, estimated GF’s follow from Eq. (24):

G​(𝐱,s|𝐱′,s′)=p​(𝐱,s|𝐱′,s′)=n​(𝐱,t)N𝐱′​t′​Δ​A​(𝐱,t)+O​(Δ​t)G\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=p\left(\mathbf{x},s|\mathbf{x^{\prime}},s^{\prime}\right)=\frac{n\left(\mathbf{x},t\right)}{N_{\mathbf{x^{\prime}}t^{\prime}}\Delta A\left(\mathbf{x},t\right)}+O\left({\Delta t}\right) (39)

Following [20], the components of the dispersion tensor, 𝐃,\mathbf{D}, are given by Dx​x=D0[a1x+a2)2ψ1(mt)],D_{xx}=D_{0}\left[a_{1}x+a_{2})^{2}\psi_{1}\left(mt\right)\right], Dy​y=D0[b1y+b2)2ψ2(mt)],D_{yy}=D_{0}\left[b_{1}y+b_{2})^{2}\psi_{2}\left(mt\right)\right], with off-diagonal terms, Dx​y=Dy​x=0.D_{xy}=D_{yx}=0. Note again that stochastic estimation of Green’s function is limited to diffusion and dispersion problems characterized by diffusive or dispersive transport dominant in one, two, or three mutually orthogonal directions. Following [20], the ground water flow field is given by vx=v0(a1x+a2)ψ1(mt))v_{x}=v_{0}\left(a_{1}x+a_{2})\psi_{1}\left(mt\right)\right) and vy=v0​(b1​y+b2)​ψ2​(m​t),v_{y}=v_{0}\left(b_{1}y+b_{2}\right)\psi_{2}\left(mt\right), with ψ1​(m​t)\psi_{1}\left(mt\right) and ψ2​(m​t)\psi_{2}\left(mt\right) being two arbitrary functions of time. For purposes of this test, we set D0=0.05D_{0}=0.05, a1=b1=1a_{1}=b_{1}=1, a2=b2=0a_{2}=b_{2}=0, v0=0.2v_{0}=0.2 and ψ1​(m​t)=ψ2​(m​t)=1\psi_{1}\left(mt\right)=\psi_{2}\left(mt\right)=1. Finally, in the analytical solution presented by [20], we set the depletion constant, γ=0.5.\gamma=0.5.

Refer to caption
Figure 7: Test 3: Exact Green’s functions on a semi-infinite domain, applicable to advective-dispersive-reactive ground water transport in a semi-infinite, anisotropic region. Here, the advective flow field is directed at 45o45^{o} from horizontal and accelerates linearly at a fixed rate, both from left to right and from bottom to top. Green’s functions, constructed using an exact ground water transport solution [20], are shown at: (a) t′=1,t^{\prime}=1, (b) t′=3,t^{\prime}=3, (c) t′=5,t^{\prime}=5, and (d) t′=9.t^{\prime}=9.
Refer to caption
Figure 8: Test 3: Estimated, smoothed Green’s functions, constructed using a backward time step Δ​s′=10−4.\Delta s^{\prime}=10^{-4}. Forward time instants, t′,t^{\prime}, maximum relative errors, em​a​x,e_{max}, and maximum smoothing window sizes, nm​a​x=ls′/(2​Δ​x′),n_{max}=l_{s}^{\prime}/\left(2\Delta x^{\prime}\right), where ls′l_{s}^{\prime} is the length of the square smoothing window, are, respectively: (a) t′=1.0,t^{\prime}=1.0, em​a​x=0.63%,e_{max}=0.63\%, no smoothing, (b) t′=3.0,t^{\prime}=3.0, em​a​x=1.73%,e_{max}=1.73\%, nm​a​x=1,n_{max}=1, (c) t′=5.0,t^{\prime}=5.0, em​a​x=0.52%,e_{max}=0.52\%, nm​a​x=1,n_{max}=1, and (d) t′=9.0,t^{\prime}=9.0, em​a​x=0.15%,e_{max}=0.15\%, nm​a​x=2.n_{max}=2.

In test 3, two features allow us to bypass the random walker respawning algorithm, and in addition, reduce the solution domain. First, since the physical problem domain is infinite, RW’s, in reality, never reach a domain boundary. In order to avoid specification of domain dimensions that are too small, allowing unphysical random walker absorption, various approaches can be used. For example, one can choose boundary distances, db​o​u​n​d​a​r​y,d_{boundary}, that are some greater-than-unity multiple of the characteristic advective and diffusive transport distance, vc​h​a​r​T+Dc​h​a​r​T,v_{char}T+\sqrt{D_{char}T}, where vc​h​a​rv_{char} and Dc​h​a​rD_{char} are characteristic velocity and diffusion/dispersivity magnitudes, and TT is the chosen total solution time. As described immediately below, we use a second approach, representing a refined, and likely more economical variation on this strategy.

Second, since the ground water velocity field and the position-dependent dispersivity both tend to zero as the origin, (x,y)=(0,0)\left(x,y\right)=\left(0,0\right) is approached, very few RW’s exit the first/positive quadrant of the infinite space. Thus, for this test, calculations are limited to the first quadrant and, based on numerical experiments showing that random walker respawning does not significantly improve Green’s function estimation accuracy, respawning is turned off.

In order to define the solution domain, at the last forward time instant, t=T=10,t=T=10, we determine the average x−x- and y−y- positions of the random walker swarm launched from the origin, μx​(t=T)\mu_{x}\left(t=T\right) and μy​(t=T),\mu_{y}\left(t=T\right), as well as standard deviations in these averages, σx​(t=T)\sigma_{x}\left(t=T\right) and σy​(t=T).\sigma_{y}\left(t=T\right). The non-normalized (but dimensionless) extent of the problem domain is then defined as: 0<x<μx+σx0<x<\mu_{x}+\sigma_{x} and 0<y<μy+σy0<y<\mu_{y}+\sigma_{y}. The solution domain is discretized into 10410^{4} interrogation areas, but with different fixed sizes in the x−x- and y−y-directions (since μx+σx≠μy+σy)\mu_{x}+\sigma_{x}\neq\mu_{y}+\sigma_{y}). As in tests 1 and 2, we choose Δ​t=10−4\Delta t=10^{-4} and N𝐱′,t′=106.N_{\mathbf{x^{\prime}},t^{\prime}}=10^{6}. Finally, the same smoothing technique used in tests 1 and 2 is applied here.

Stochastically estimated GF’s are plotted in Figure 8 at the four forward time instants: t=1,t=1, t=3,t=3, t=5,t=5, and t=9.t=9. For comparison, exact GF’s, determined from [20], are shown in Fig. 7. Again, quality Green’s function estimates are obtained at all four probed time instants, validating the use of forward time estimation. Again, under most circumstances, however, backward time estimation is recommended since one Green’s function estimate, obtained at a single response point, (𝐱,t),\left(\mathbf{x},t\right), can be used in the magic rule to estimate the evolution of η​(𝐱,t),\eta\left(\mathbf{x},t\right), subject to any set of forcing conditions.

8 Discussion

When they can be found, Green’s functions provide a powerful tool for solving linear - and incrementally, nonlinear - evolution problems having the generic forms given by (1) and (2). Unfortunately, finding GF’s can be extremely difficult, particularly within geometrically complicated regions.

This paper proposes a numerical method for stochastically estimating GF’s, appropriate to diffusive and advective-dispersive/diffusive-reactive transport of scalar- and vector-valued quantities, e.g., convective heat and mass transfer problems involving passive, as well as chemically-reacting species. We validate the technique against three known GF’s, two that apply to simple diffusion in finite regions, and a third applicable to ground water transport in an infinite region. As formulated here, the stochastic Green’s function estimation technique can be applied to scalar transport in independently determined turbulent flows. More fundamentally, the technique might also be adapted to generate, in time-incremental fashion, GF’s for computing turbulent flows. In this case: i) the approach is limited to conditions where fluid particle momentum and vorticity dynamics are strongly Brownian, ii) the transported quantity corresponds to either local fluid momentum or vorticity [21], and iii) the turbulent flow problem is parabolic, evolving from a known initial condition.

We identify three principle challenges and propose straightforward solutions for each:

a) Since the probability of observing RW’s at Dirichlet boundaries is zero, RW’s that reach such boundaries must be removed. Random walker absorption, in turn, quickly compromises the quality of estimated GF’s. In order to address this problem, we introduce a random walker respawning technique which conserves the number of RW’s initially launched, as well as enforcing continuum mass conservation and diffusion at each respawning location.

b) Estimated Green’s functions exhibit spatial graininess, an artifact associated with use of discrete naive transition density estimators [38], Eq. (13). In order to solve this problem, as well as improve Green’s function estimation accuracy, we propose an area-averaging technique that bins neighboring Green’s function estimates, taken over small smoothing windows encompassing any given interrogation area, Δ​A′.\Delta A^{\prime}. For most points, 𝐱′,\mathbf{x^{\prime}}, not near a boundary, the size of the smoothing window, a,a, is determined in iterative fashion by minimizing the point-wise deviation between numerical GF’s and averaged GF’s. In this first exploratory study, we use known rather than averaged GF’s in determining a.a.

c) To stochastically estimate Green’s functions, as given by Eq. (24) (no respawning, test 3) or (25) (with respawning, tests 1 and 2), three numerical parameters must be chosen: i) the number of RW’s, N𝐱,sN_{\mathbf{x},s} (N𝐱′,t′(N_{\mathbf{x^{\prime}},t^{\prime}} for forward time random walker integration), launched from the chosen solution/response point, (𝐱,s)\left(\mathbf{x},s\right) [chosen impulse point, (𝐱′,t′)]\left(\mathbf{x^{\prime}},t^{\prime}\right)] ii) the backward time-step size, Δ​s′\Delta s^{\prime} [forward time-step size, Δt],\Delta t], and iii) the random walker interrogation area, Δ​A′\Delta A^{\prime} [Δ​A].[\Delta A]. Appendix A outlines approaches for determining these parameters; to the best of our knowledge, and beyond [44], the literature does not offer much guidance on this important practical question.

Although not treated in this paper, the proposed stochastic Green’s function construction technique provides a dual framework for tackling coupled particle and continuum transport problems, important, for example, in microscale bioheat and biomass transport problems [44], and a feature absent in BEM, GEM, and FEM formulations. The Green’s function estimation recipe shows that this particle-continuum duality arises from the connection between stochastic differential equations modeling particle-like transport and associated Fokker-Planck and backward Kolmogorov equations, probabilistically describing the same transport.

Acknowledgements

This work was supported by the Office of Naval Research, grant no. N00014-18-1-2754. Technical and coding support provided by Nikita Nikulsen is gratefully acknowledged.

References

  • [1] K. Cole, J. Beck, A. Haji-Sheikh, B. Litkouhi, Heat Conduction Using Greens Functions, Series in Computational Methods and Physical Processes in Mechanics and Thermal Sciences, CRC Press, 2010.
  • [2] G. Barton, Elements of Green’s functions and propagation: potentials, diffusion, and waves, Oxford University Press, 1989.
  • [3] Y. Hon, M. Li, Y. Melnikov, Inverse source identification by green’s function, Engineering Analysis with Boundary Elements 34 (4) (2010) 352–358. doi:https://doi.org/10.1016/j.enganabound.2009.09.009.
  • [4] J. Beck, B. Blackwell, A. Haji-Sheikh, Comparison of some inverse heat conduction methods using experimental data, International Journal of Heat and Mass Transfer 39 (17) (1996) 3649–3657. doi:https://doi.org/10.1016/0017-9310(96)00034-8.
  • [5] A. P. Fernandes, M. B. dos Santos, G. Guimarães, An analytical transfer function method to solve inverse heat conduction problems, Applied Mathematical Modelling 39 (22) (2015) 6897–6914. doi:https://doi.org/10.1016/j.apm.2015.02.012.
  • [6] D. Menemenlis, I. Fukumori, T. Lee, Using green’s functions to calibrate an ocean general circulation model, Monthly Weather Review 133 (2005) 1224–1240.
  • [7] S. Boutami, S. Fan, Efficient pixel-by-pixel optimization of photonic devices utilizing the dyson’s equation in a green’s function formalism: Part i. implementation with the method of discrete dipole approximation, J. Opt. Soc. Am. B 36 (2019) 2378–2386.
  • [8] S. Boutami, S. Fan, Efficient pixel-by-pixel optimization of photonic devices utilizing the dyson’s equation in a green’s function formalism: Part ii. implementation using standard electromagnetic solvers, J. Opt. Soc. Am. B 36 (2019) 2387–2394.
  • [9] E. N. Economou, Green’s Functions in Quantum Physics, Springer Series in Solid-State Sciences, Springer, Berlin, 2006.
  • [10] D. J. David Joseph Biagioni, Numerical construction of Green’s functions in high dimensional elliptic problems with variable coefficients and analysis of renewable energy data via sparse and separable approximations, Ph.D. Dissertation, Applied Mathemetics, University of Colorado, Boulder, 2012.
  • [11] M. Costabel, Principles of boundary element methods, Computer Physics Reports 6 (1987) 243–274.
  • [12] C. Brebbia, S. Walker, Boundary element techniques in engineering, Newnes-Butterworths, 1980.
  • [13] A. E. Taigbenu, The green element method, International Journal for Numerical Methods in Engineering 38 (13) (1995) 2241–2263. doi:https://doi.org/10.1002/nme.1620381307.
  • [14] Y. Wu, L. Cheng, S. Fang, S. Huang, P. Jia, A green element method-based discrete fracture model for simulation of the transient flow in heterogeneous fractured porous media, Advances in Water Resources 136 (2020) 103489. doi:https://doi.org/10.1016/j.advwatres.2019.103489.
  • [15] P. Lorinczi, S. Harris, L. Elliott, Modelling of highly-heterogeneous media using a flux-vector-based green element method, Engineering Analysis with Boundary Elements 30 (2006) 818–833. doi:https://doi.org/10.1016/j.enganabound.2006.07.004.
  • [16] O. C. Zienkiewicz, The finite element method in engineering science, McGraw-Hill, 1971.
  • [17] K.-J. Bathe, Finite Element Method, John Wiley & Sons, Ltd, 2008. doi:https://doi.org/10.1002/9780470050118.ecse159.
  • [18] D. J N Reddy, An Introduction to the Finite Element Method, McGraw-Hill Education, 2005.
  • [19] T. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Civil and Mechanical Engineering, Dover Publications, 2012.
  • [20] A. Sanskrityayn, V. P. Singh, V. K. Bharati, N. Kumar, Analytical solution of two-dimensional advection–dispersion equation with spatio-temporal coefficients for point sources in an infinite medium using green’s function method, Environmental Fluid Mechanics 18 (3) (2018) 739–757.
  • [21] R. G. Keanini, Green’s function-stochastic methods framework for probing nonlinear evolution problems: Burger’s equation, the nonlinear schrödinger’s equation, and hydrodynamic organization of near-molecular-scale vorticity, Annals of Physics 326 (4) (2011) 1002–1031.
  • [22] M. Jovanovic, B. Bamieh, Modeling flow statistics using the linearized navier-stokes equations, in: Decision and Control, 2001. Proceedings of the 40th IEEE Conference on, Vol. 5, IEEE, 2001, pp. 4944–4949.
  • [23] V. G. Levich, Physicochemical hydrodynamics, Prentice hall, 1962.
  • [24] F. A. Williams, Combustion Theory: the fundamental theory of chemical reacting flow systems, Addison-Wesley, 1965.
  • [25] P. Grindrod, The Theory and Applications of Reaction-diffusion Equations: Patterns and Waves, Oxford applied mathematics and computing science series, Clarendon Press, 1996.
  • [26] A. Messiah, Quantum Mechanics, Dover Books on Physics, Dover Publications, 2014.
  • [27] L. D. Landau, E. Lifshitz, L. Pitaevskij, Course of theoretical physics. vol. 10: Physical kinetics, Oxford, 1981.
  • [28] G. Colonna, Boltzmann and vlasov equations in plasma physics, in: Plasma Modeling, 2053-2563, IOP Publishing, 2016, pp. 1–1 to 1–23. doi:10.1088/978-0-7503-1200-4ch1.
    URL http://dx.doi.org/10.1088/978-0-7503-1200-4ch1
  • [29] G. B. Whitham, Linear and nonlinear waves, Vol. 42, John Wiley & Sons, 2011.
  • [30] P. G. Drazin, R. S. Johnson, Solitons: an introduction, Vol. 2, Cambridge university press, 1989.
  • [31] C. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry, and the Natural Sciences, Springer complexity, Springer, 2004.
  • [32] M. E. Muller, Some continuous monte carlo methods for the dirichlet problem, The Annals of Mathematical Statistics 27 569–589.
  • [33] K. K. Sabelfeld, Monte Carlo methods in boundary value problems, Springer series in computational physics, Springer, Berlin, 1991.
  • [34] T. E. Booth, Exact monte carlo solution of elliptic partial differential equations, Journal of Computational Physics 39 (1981) 396–404. doi:10.1016/0021-9991(81)90159-5.
  • [35] K. K. Sabelfeld, N. A. Simonov, Random walks on boundary for solving PDEs, De Gruyter, 2013.
  • [36] A. Haji-Sheikh, E. M. Sparrow, The floating random walk and its application to monte carlo solutions of heat equations, SIAM Journal on Applied Mathematics 14 (2) (1966) 370–389.
  • [37] M. Mascagni, N. A. Simonov, Monte carlo methods for calculating some physical properties of large molecules, SIAM journal on scientific computing 26 (1) (2004) 339–357.
  • [38] B. W. Silverman, Density Estimation for Statistics and Data Analysis, Chapman & Hall, 1986.
  • [39] G. N. Milstein, J. G. M. Schoenmakers, V. Spokoiny, Transition density estimation for stochastic differential equations via forward-reverse representations, Bernoulli 10 (2004) 281–312.
  • [40] L. Devroye, L. Gyorfi, Nonparametric Density Estimation: The L1 View, Wiley Interscience Series in Discrete Mathematics, Wiley, 1985.
  • [41] E. Olariu, K. K. Cadwell, E. Hancock, D. Trueman, H. Chevrou-Severac, Current recommendations on the estimation of transition probabilities in markov cohort models for use in health care decision-making: a targeted literature review, Clinicoecon Outcomes Res. 9 (2017) 537–546.
  • [42] M. Rudemo, Empirical choice of histograms and kernel density estimators, Scandinavian Journal of Statistics 9 (1982) 65–78.
  • [43] Z. Schuss, Theory and Applications of Stochastic Processes: An Analytical Approach, Applied Mathematical Sciences, Springer New York, 2009.
  • [44] Z. Schuss, Brownian Dynamics at Boundaries and Interfaces: In Physics, Chemistry, and Biology, Applied Mathematical Sciences, Springer New York, 2013.
    URL https://books.google.com/books?id=eOi3BAAAQBAJ
  • [45] G. N. Milstein, Approximate integration of stochastic differential equations, Theory of Probability & Its Applications 19 (1975) 557–562.
  • [46] G. Milstein, Numerical Integration of Stochastic Differential Equations, Mathematics and Its Applications, Springer Netherlands, 1994.
  • [47] P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Stochastic Modelling and Applied Probability, Springer Berlin Heidelberg, 2011.
  • [48] J. A. Hansen, C. Penland, Efficient approximate techniques for integrating stochastic differential equations, Monthly Weather Review 134 (2006) 3006 – 3014. doi:10.1175/MWR3192.1.
  • [49] T. Sauer, Numerical solution of stochastic differential equations in finance (2012) 529–550.
  • [50] S. Ulam, R. Richtmyer, J. von Neumann, Statistical methods in neutron diffusion, Los Alamos Scientific Laboratory report LAMS–551 (1947).
  • [51] N. Metropolis, S. Ulam, The monte carlo method, J. Amer. Stat. Assn. 44 (1949) 335–341.
  • [52] J. Todd, Experiments in the solution of differential equations by monte carlo methods, J. Washington Acad. Sci. 12 (1954) 377–381.
  • [53] I. Dimnov, D. McKee, Monte Carlo Methods for Applied Scientists, World Scientific, 1998.
  • [54] R. Feynman, Principles of Least Action in Quantum Mechanics, Ph.D. Dissertation, Physics, Princeton University, 1942.
  • [55] M. Kac, On distributions of certain wiener functionals, Trans. Amer. Math. Soc. 65 (1949) 1–13.
  • [56] M. Kac, On some connections between probability theory and differential and integral equations (1951) 189–215.
  • [57] A. Friedman, Stochastic differential equations and applications, Vol. 1, Academic Press, 1975.
  • [58] R. Bass, Stochastic Processes, Cambridge Univ. Press, 2011.
  • [59] D. Forster, Hydrodynamic Fluctuations, Broken Symmetry, And Correlation Functions, CRC Press, 2018.
  • [60] J. Boon, S. Yip, Molecular Hydrodynamics, McGraw-Hill, 1980.
  • [61] M. Toda, R. Kubo, N. Saito, N. Hashitsume, Statistical Physics II: Nonequilibrium Statistical Mechanics, Springer Series in Solid-State Sciences, Springer Berlin Heidelberg, 2012.
  • [62] I. Lux, L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations (1st ed.), CRC Press, 1991.
  • [63] J. Bear, A. Cheng, Modeling Groundwater Flow and Contaminant Transport, Springer, 2010.

Supporting information

S1

Appendices to ’Stochastic estimation of Green’s functions with application to advection-diffusion-reaction problems’

S1 Appendix A

Choosing the number of RW’s launched, N𝐱′,t′,N_{\mathbf{x^{\prime}},t^{\prime}}, interrogation area, Δ​A​(𝐱,t),\Delta A\left(\mathbf{x},t\right), and the time step, Δ​t.\Delta t. The stochastic estimation technique requires specification of three numerical parameters: i) the backward (forward) time-step size, Δ​s′\Delta s^{\prime} (Δ​t),\left(\Delta t\right), ii) the interrogation area, Δ​A′​(𝐱′,s′)\Delta A^{\prime}\left(\mathbf{x}^{\prime},s^{\prime}\right) (Δ​A′​(𝐱′,s′)),\left(\Delta A^{\prime}\left(\mathbf{x}^{\prime},s^{\prime}\right)\right), and iii) the number of RW’s, N𝐱,sN_{\mathbf{x},s} (N𝐱′,t′),\left(N_{\mathbf{x^{\prime}},t^{\prime}}\right), launched from the chosen GF solution point, (𝐱,s)\left(\mathbf{x},s\right) ((𝐱′,t′)).\left(\left(\mathbf{x^{\prime}},t^{\prime}\right)\right). Since the literature appears to offer little guidance on choosing these parameters, and since the choices require some experimentation, we propose two straightforward approaches for identifying appropriate values. Here, we outline the technique used in this study. We focus on the simplest case in which the diffusion tensor is uniform and isotropic, 𝐃=Do​𝐈.\mathbf{D}=D_{o}\mathbf{I}. It is important to recognize that in both approaches, some numerical experimentation is required.

In order to address these important questions, we focus on forward-time integration of RW trajectories:

d​𝝌​(t)=𝐯​(𝝌​(t),t)​d​t+𝐁​(𝝌​(t),t)⋅d​𝐰​(t)d\boldsymbol{\chi}\left(t\right)=\mathbf{v}\left(\boldsymbol{\chi}\left(t\right),t\right)dt+\mathbf{B}\left(\boldsymbol{\chi}\left(t\right),t\right)\cdot d\mathbf{w}\left(t\right) (40)

and highlight two features: a) the connection between Δ​A​(𝐱,t)\Delta A\left(\mathbf{x},t\right) and Δ​t,\Delta t, and b) minimization of early-time variance in the number of RW’s captured on Δ​A​(𝐱,t).\Delta A\left(\mathbf{x},t\right).

For simplicity, limit attention to two-dimensional diffusive transport of η​(𝐱,t).\eta\left(\mathbf{x},t\right). Spatially resolved random walker displacements, 𝐝​𝝌​(t),\mathbf{d}\boldsymbol{\chi}\left(t\right), over Δ​t,\Delta t, require Δ​A​(𝐱,t)∼⟨d​χ​(t)⋅d​χ​(t)⟩=Do​Δ​t.\Delta A\left(\mathbf{x},t\right)\sim\langle d\mathbf{\chi}\left(t\right)\cdot d\mathbf{\chi}\left(t\right)\rangle=D_{o}\Delta t. Thus, we choose

Δ​A​(𝐱,t)=Do​Δ​t\Delta A\left(\mathbf{x},t\right)=D_{o}\Delta t (41)

where this choice represents the smallest area allowing statistically acceptable sampling of n​(𝐱,t):n\left(\mathbf{x},t\right): If Δ​A​(𝐱,t)\Delta A\left(\mathbf{x},t\right) is chosen smaller than Do​Δ​t,D_{o}\Delta t, then for any given number of RW’s launched from (𝐱′,t′),\left(\mathbf{x^{\prime}},t^{\prime}\right), N𝐱′,t′,N_{\mathbf{x^{\prime}},t^{\prime}}, a number of RW’s located near Δ​A​(𝐱,t),\Delta A\left(\mathbf{x},t\right), at forward instant t−Δ​t,t-\Delta t, will (numerically) jump over/around Δ​A​(𝐱,t)\Delta A\left(\mathbf{x},t\right) at t.t. The result is an underestimated n​(𝐱,t).n\left(\mathbf{x},t\right). Likewise, if Δ​A​(𝐱,t)\Delta A\left(\mathbf{x},t\right) is chosen larger than Do​Δ​t,D_{o}\Delta t, the number captured on Δ​A​(𝐱,t)\Delta A\left(\mathbf{x},t\right) is too large, overestimating n​(𝐱,t).n\left(\mathbf{x},t\right).

In order to determine acceptable N𝐱′,t′,N_{\mathbf{x^{\prime}},t^{\prime}}, consider the time-dependent standard deviation, σn,\sigma_{n}, in n​(𝐱,t),n\left(\mathbf{x},t\right), in the number of RW’s captured by the set of interrogation areas, {Δ​A1​(r​(τ)),Δ​A2​(r​(τ)),…​Δ​ANΔ​A​(r​(τ))},\{\Delta A_{1}\left(r\left(\tau\right)\right),\Delta A_{2}\left(r\left(\tau\right)\right),...\Delta A_{N_{\Delta A}}\left(r\left(\tau\right)\right)\}, lying on the circle, r=Do​τ,r=\sqrt{D_{o}\tau}, centered on the launch location, 𝐱′:\mathbf{x^{\prime}}:

σn​(τ)=1NΔ​A​(r)​∑i=1NΔ​A​(r)[ni​(r)−n¯​(r)]2\sigma_{n}\left(\tau\right)=\sqrt{\frac{1}{N_{\Delta A}\left(r\right)}\sum_{i=1}^{N_{\Delta A}\left(r\right)}\left[n_{i}\left(r\right)-\overline{n}\left(r\right)\right]^{2}} (42)

where r=r​(τ)r=r\left(\tau\right) is the mean penetration radius of the N𝐱′,t′N_{\mathbf{x^{\prime}},t^{\prime}} RW’s launched from 𝐱′\mathbf{x^{\prime}} at time t′,t^{\prime}, τ=t−t′,\tau=t-t^{\prime}, and n¯​(r)=NΔ​A−1​∑i=1NΔ​Ani​(r).\overline{n}\left(r\right)=N_{\Delta A}^{-1}\sum_{i=1}^{N_{\Delta A}}n_{i}\left(r\right).

The number of interrogation areas lying on r​(τ),r\left(\tau\right), NΔ​A,N_{\Delta A}, is given approximately by:

NΔ​A​Δ​A𝐱,t∼r​(τ)N_{\Delta A}\sqrt{\Delta A_{\mathbf{x},t}}\sim r\left(\tau\right) (43)

where, in the case that the the spatial domain is discretized into uniform interrogation areas, Δ​A𝐱,t=Δ​x⋅Δ​y,\Delta A_{\mathbf{x},t}=\Delta x\cdot\Delta y, with Δ​x≈Δ​y,\Delta x\approx\Delta y, Δ​A𝐱,t∼Δ​x≈Δ​y.\sqrt{\Delta A_{\mathbf{x},t}}\sim\Delta x\approx\Delta y. Since Δ​x≈Δ​y∼Do​Δ​t,\Delta x\approx\Delta y\sim\sqrt{D_{o}\Delta t}, then

NΔ​A∼Do​τΔ​x∼Do​τDo​Δ​t=τΔ​tN_{\Delta A}\sim\frac{\sqrt{D_{o}\tau}}{\Delta x}\sim\frac{\sqrt{D_{o}\tau}}{\sqrt{D_{o}\Delta t}}=\sqrt{\frac{\tau}{\Delta t}} (44)

From Eq. (42), and for small interrogation areas, Δ​A𝐱,t<<Ao,\Delta A_{\mathbf{x},t}<<A_{o}, AoA_{o} being the area of the solution domain,

σn​(τ)∼1NΔ​A∼(Δ​tτ)1/4\sigma_{n}\left(\tau\right)\sim\frac{1}{\sqrt{N_{\Delta A}}}\sim\left(\frac{\Delta t}{\tau}\right)^{1/4} (45)

From this expression, it becomes apparent that the maximum spatial variation, σn​(τ),\sigma_{n}\left(\tau\right), in random walker capture numbers, n​(𝐱,t),n\left(\mathbf{x},t\right), occurs at τm​i​n=Δ​t,\tau_{min}=\Delta t, where σn​(Δ​t)∼1.\sigma_{n}\left(\Delta t\right)\sim 1.

In order determine an acceptable N𝐱′,t′,N_{\mathbf{x^{\prime}},t^{\prime}}, we calculate the ratio of the mean number of RW’s captured at τ=Δ​t,\tau=\Delta t, ⟨n​(𝐱,Δ​t)⟩,\langle n\left(\mathbf{x},\Delta t\right)\rangle, to N𝐱′,t′:N_{\mathbf{x^{\prime}},t^{\prime}}:

⟨n​(𝐱,Δ​t)⟩N𝐱′,t′=P[𝝌(Δt)∈ΔA(𝐱,Δt)|𝝌(t′)=𝐱′,t′]\displaystyle\frac{\langle n\left(\mathbf{x},\Delta t\right)\rangle}{N_{\mathbf{x^{\prime}},t^{\prime}}}=P\big{[}\boldsymbol{\chi}\left(\Delta t\right)\in\Delta A\left(\mathbf{x},\Delta t\right)|\boldsymbol{\chi}\left(t^{\prime}\right)=\mathbf{x^{\prime}},t^{\prime}\big{]}
=∫Δ​A​(𝐱,Δ​t)14​π​Do​Δ​t​exp⁡[−(𝐱−𝐱′)24​Do​Δ​t]​𝑑𝐱\displaystyle\hskip-199.16928pt=\int_{\Delta A\left(\mathbf{x},\Delta t\right)}\frac{1}{4\pi D_{o}\Delta t}\exp\bigg{[}\frac{-\left(\mathbf{x}-\mathbf{x^{\prime}}\right)^{2}}{4D_{o}\Delta t}\bigg{]}d\mathbf{x} (46)

or

⟨n​(𝐱𝐣,Δ​t)⟩N𝐱′,t′=1π​exp⁡[−(𝐱𝐣−𝐱′)24​Do​Δ​t]​Δ​x​Δ​y4​Do​Δ​t∼Δ​x​Δ​y4​π​e1/4​Do​Δ​t\frac{\langle n\left(\mathbf{x_{j}},\Delta t\right)\rangle}{N_{\mathbf{x^{\prime}},t^{\prime}}}=\frac{1}{\pi}\exp\bigg{[}\frac{-\left(\mathbf{x_{j}}-\mathbf{x^{\prime}}\right)^{2}}{4D_{o}\Delta t}\bigg{]}\frac{\Delta x\Delta y}{4D_{o}\Delta t}\sim\frac{\Delta x\Delta y}{4\pi e^{1/4}D_{o}\Delta t} (47)

where 𝐱𝐣\mathbf{x_{j}} represents, e.g., the centroid of the jt​hj^{th} interrogation area, Δ​A𝐱𝐣,Δ​t,\Delta A_{\mathbf{x_{j}},\Delta t}, lying on r​(Δ​t)=Do​Δ​t,r\left(\Delta t\right)=D_{o}\Delta t, and where (𝐱𝐣−𝐱′)2∼Do​Δ​t.\left(\mathbf{x_{j}}-\mathbf{x^{\prime}}\right)^{2}\sim D_{o}\Delta t.

Finally, solving for ⟨n​(𝐱𝐣,Δ​t)⟩,\langle n\left(\mathbf{x_{j}},\Delta t\right)\rangle, we obtain the maximum relative spatial variation in captured RW’s, expressed as a function of N𝐱′,t′N_{\mathbf{x^{\prime}},t^{\prime}} and Δ​A𝐱,t=Δ​x​Δ​y:\Delta A_{\mathbf{x},t}=\Delta x\Delta y:

σn​(r)|m​a​x⟨n​(𝐱𝐣,Δ​t)⟩∼4​π​e1/4​Do​Δ​tN𝐱𝐣,Δ​t​Δ​x​Δ​y\frac{\sigma_{n}\left(r\right)|_{max}}{\langle n\left(\mathbf{x_{j}},\Delta t\right)\rangle}\sim\frac{4\pi e^{1/4}D_{o}\Delta t}{N_{\mathbf{x_{j}},\Delta t}\Delta x\Delta y} (48)

Thus, the following procedure can be used to determine Δ​A𝐱′,t′,\Delta A_{\mathbf{x^{\prime}},t^{\prime}}, Δ​t,\Delta t, and N𝐱′,t′:N_{\mathbf{x^{\prime}},t^{\prime}}:

a) Based on the nature of the transport problem under study, identify the spatial resolution required; this choice determines Δ​x,\Delta x, Δ​y,\Delta y, and Δ​A𝐱,t=Δ​x​Δ​y.\Delta A_{\mathbf{x},t}=\Delta x\Delta y.

b) Use Eq. (41) to determine Δ​t;\Delta t; for problems where Δ​x≈Δ​y,\Delta x\approx\Delta y, Δ​t=Δ​x/Do≈Δ​y/Do.\Delta t=\sqrt{\Delta x}/D_{o}\approx\sqrt{\Delta y}/D_{o}.

c) Choose a desired maximum relative (spatial, early-time) variation in RW’s captured on Δ​A𝐱,t,\Delta A_{\mathbf{x},t}, σn​(Δ​t)/⟨n​(𝐱𝐣,Δ​t)⟩,\sigma_{n}\left(\Delta t\right)/\langle n\left(\mathbf{x_{j}},\Delta t\right)\rangle, and solve Eq. (48) for N𝐱′,t′.N_{\mathbf{x^{\prime}},t^{\prime}}.

Note 1: Since all RW’s are independent, one can effectively increase the number of RW’s launched from a given solution point, (𝐱,s)(\mathbf{x},s) - which, for single-swarm simulations, we denote as N𝐱,sN_{\mathbf{x},s} - by relaunching another swarm of M𝐱,sM_{\mathbf{x},s} RW’s from the same point. The number used in the relaunch, M𝐱,s,M_{\mathbf{x},s}, can be larger or smaller than N𝐱,s;N_{\mathbf{x},s}; the corresponding numbers of RW’s captured in each interrogation area, Δ​A′,\Delta A^{\prime}, n​(𝐱′,s′)n\left(\mathbf{x}^{\prime},s^{\prime}\right) and m​(𝐱′,s′),m\left(\mathbf{x}^{\prime},s^{\prime}\right), respectively, can then be pooled to improve GF estimation accuracy.

Note 2: As long as dimensional lengths are scaled by the either the short-time or long-time diffusion length scales, Do​Δ​t\sqrt{D_{o}\Delta t} or Do​T,\sqrt{D_{o}T}, where TT is the chosen solution time interval, this procedure applies for GF estimation procedures cast in either dimensional or nondimensional form.


S1 Appendix B

Properties of the respawning algorithm. In order to focus on the computational challenges associated with absorbing boundaries, in this subsection we limit attention to pure diffusion problems in which no advection takes place. Thus, the backward-time velocity field, 𝐛,\mathbf{b}, is set to zero. The proposed respawning algorithm has three properties that appear to be essential to estimating low-error GF’s within bounded domains:

a) At all times and all locations in the solution domain, and with an error on the order of O​(Δ​x′2,Δ​s′),O\left(\Delta{x^{\prime}}^{2},\Delta s^{\prime}\right), the algorithm conserves random walker mass, which, in turn, enforces the adjoint equation, (26).

In order to show that the respawning algorithm enforces the adjoint equation, label the location where the random walker having maximum weight at any instant s′k>s,{s^{\prime}}_{k}>s, as 𝐱′j.{\mathbf{x}^{\prime}}_{j}. Define a local (x​”,y​”)−coordinate\left(x",y"\right)-\mathrm{coordinate} system, having origin at 𝐱′j,{\mathbf{x}^{\prime}}_{j}, such that the x​”−axisx"-\mathrm{axis} is aligned with the local gradient vector, ∇′Λ,\nabla^{\prime}\Lambda, in the random walker mass density field, Λ=Λ​(𝐱′,s′).\Lambda=\Lambda\left(\mathbf{x}^{\prime},s^{\prime}\right). The continuous field, Λ​(𝐱′,s′),\Lambda\left(\mathbf{x}^{\prime},s^{\prime}\right), can be constructed, for example, by interpolation on the discrete instantaneous distribution of weighted RW’s.

The process of removing a maximum-weight random walker, having weight Wm​a​x​(s′k),W_{max}\left({s^{\prime}}_{k}\right), located at 𝐱′j,{\mathbf{x}^{\prime}}_{j}, and replacing it with two RW’s, each of weight Wm​a​x​(s′k)/2,W_{max}\left({s^{\prime}}_{k}\right)/2, and both placed at 𝐱′j,{\mathbf{x}^{\prime}}_{j}, can be represented by the following expression:

12​[Λ​(𝐱​”j,s′k|𝐱,s)+Λ​(𝐱​”j,s′k|𝐱,s)]−Λ​(𝐱​”j,s′k|𝐱,s)=0\frac{1}{2}\left[\Lambda\left({\mathbf{x}"}_{j},{s^{\prime}}_{k}|\mathbf{x},s\right)+\Lambda\left({\mathbf{x}"}_{j},{s^{\prime}}_{k}|\mathbf{x},s\right)\right]-\Lambda\left({\mathbf{x}"}_{j},{s^{\prime}}_{k}|\mathbf{x},s\right)=0 (49)

random walker mass is clearly conserved since the left hand side of (49) is identically zero. Taylor expanding the first term in the positive x​”−x"- direction, the second term in the negative x​”−x"- direction, and rearranging, we obtain:

Λ​(𝐱​”j,s′k|𝐱,s)=12​[Λ​(𝐱​”j+1,s′k−1|𝐱,s)+Λ​(𝐱​”j−1,s′k−1|𝐱,s)]+O​(Δ​x​”2,Δ​s′=Δ​x​”2/D)\Lambda\left({\mathbf{x}"}_{j},{s^{\prime}}_{k}|\mathbf{x},s\right)=\frac{1}{2}\left[\Lambda\left({\mathbf{x}"}_{j+1},{s^{\prime}}_{k-1}|\mathbf{x},s\right)+\Lambda\left({\mathbf{x}"}_{j-1},{s^{\prime}}_{k-1}|\mathbf{x},s\right)\right]+O\left({\Delta x"}^{2},\Delta s^{\prime}=\Delta{x"}^{2}/D\right) (50)

where, in the limit Δ​s′→0,\Delta s^{\prime}\rightarrow 0, D​Δ​s′=Δ​x​”2,D\Delta s^{\prime}=\Delta{x"}^{2}, with Δ​x​”=x​”j+1−x​”j=x​”j−x​”j−1,\Delta{x"}={x"}_{j+1}-{x"}_{j}={x"}_{j}-{x"}_{j-1}, Δ​s′=s′k−s′k−1,\Delta s^{\prime}={s^{\prime}}_{k}-{s^{\prime}}_{k-1}, and where the order of discretization error is as shown.

Finally, expressing the left hand side of (50) as a temporal Taylor expansion about s′=s′k−1,s^{\prime}={s^{\prime}}_{k-1}, and the two terms on the right as spatial expansions about x​”j,{x"}_{j}, we obtain the local mass diffusion equation, evaluated at (𝐱′j,s′k−1),\left({\mathbf{x}^{\prime}}_{j},{s^{\prime}}_{k-1}\right), and expressed in the backward time direction

−∂Λ∂s′=D​∂2Λ∂x​”2+O​(Δ​x​”2,Δ​s′=Δ​x​”2/D)-\frac{\partial\Lambda}{\partial s^{\prime}}=D\frac{\partial^{2}\Lambda}{\partial{x"}^{2}}+O\left({\Delta x"}^{2},\Delta s^{\prime}=\Delta{x"}^{2}/D\right) (51)

Note, that replacing the time index kk with k+1k+1 shifts evaluation of (51) to (𝐱′j,s′k).\left({\mathbf{x}^{\prime}}_{j},{s^{\prime}}_{k}\right).

b) If we use the RW’s to model physical Brownian particles, and the respawning algorithm to model some process by which the particles can undergo binary fission, we see that it enforces the diffusion equation (51) in the direction of the local random walker mass flux, ∇Λ/|∇Λ|.\nabla\Lambda/\left|\nabla\Lambda\right|. Since physical particles, unlike numerical RW’s, cannot both occupy the same point that their parent occupied, we would introduce a spatial offset for the respawned RW’s. Thus, the fission of any given set of maximum-weight RW’s, {χm​a​x,1​(s′k),χm​a​x,2​(s′k),…,χm​a​x,M​(s′k)},\{\chi_{max,1}\left({s^{\prime}}_{k}\right),\chi_{max,2}\left({s^{\prime}}_{k}\right),...,\chi_{max,M}\left({s^{\prime}}_{k}\right)\}, produces a mass gradient at each of the M=M​(s′)M=M\left(s^{\prime}\right) respawning locations, χq​(s′),\chi_{q}\left(s^{\prime}\right), q=1,2,3,..,M(s′),q=1,2,3,..,M\left(s^{\prime}\right), and the algorithm, in effect, smooths - on short (single-time-step) time-scales, Δ​x′2/D=Δ​s′\Delta{x^{\prime}}^{2}/D=\Delta s^{\prime} - localized, and randomly-distributed (in space) removal of random walker mass. Hence, the algorithm appears to be insensitive to the method used to reorder, at the end of every nn time-steps, the weight array, {W}.\{W\}. Although this remains an open question, our results strongly support this view.

c) As noted in the article, the third property associated with the respawning algorithm centers on the fact that it conserves the number of RW’s launched from the solution point, N𝐱,s.N_{\mathbf{x},s}. As shown in the article, this property ensures that accurate GF estimates can be obtained throughout any given time interval, 0<t′≤T.0<t^{\prime}\leq T.


S1 Appendix C

Illustration of trial and error derivation of the adjoint equation and magic rule. A technical hurdle that likely limits widespread use of GF’s, particularly among non-specialists, revolves around the origins of the adjoint equation governing the Green’s function, and the so-called ’magic rule’ [barton], the integral solution for the transport variable, η​(𝐱,t).\eta\left(\mathbf{x},t\right). Since textbooks and research articles appear to be universally mute on how these essential building blocks are obtained - specifically, not showing how they can be simultaneously derived via straightforward trial and error - we illustrate trial and error determination of these the adjoint equation and magic rule, focusing for simplicity on fixed diffusivity advection-diffusion problems, evolving in incompressible flow fields:

∂η∂t+𝐯⋅∇η−D​∇𝟐η=f​(𝐱,t)\frac{\partial\eta}{\partial t}+\mathbf{v}\cdot\nabla\eta-D\boldsymbol{\nabla^{2}}\eta=f(\mathbf{x},t) (52)

The procedure incorporates four steps, the first of which is generally not highlighted: a) guess the form of the adjoint equation, b) multiply the candidate adjoint equation by η​(𝐱,t),\eta\left(\mathbf{x},t\right), and the transport equation, Eq. (52) by the Green’s function, G(𝐱,t|𝐱′,t′)),G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\right), c) integrate both equations in b) over the space-time hypersurface, δ​Ω×(0,t],\delta\Omega\times(0,t], formed by the product of the problem’s spatial boundary, δ​Ω,\delta\Omega, and the time axis, and d) taking the difference of the integral equations in c), use Green’s theorems and trial and error alteration of the guessed adjoint equation to arrive at an integral solution for η​(𝐱,t),\eta\left(\mathbf{x},t\right), stated strictly in terms of space-time surface integrals over δ​Ω×(0,t],\delta\Omega\times(0,t], the former involving G(𝐱,t|𝐱′,t′)),G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\right), derivatives of G(𝐱,t|𝐱′,t′)),G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\right), and known boundary and initial conditions.

In step a), it’s generally good strategy to guess that the adjoint equation has the same generic form as the transport equation, Eq. (52), but with one or more signs reversed, and with the source term, f​(𝐱,t),f(\mathbf{x},t), replaced by ±δ​(𝐱−𝐱′)​δ​(t−t′):\pm\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\delta\left(t-t^{\prime}\right):

∂G(𝐱,t|𝐱′,t′))∂t′+𝐯⋅∇′G(𝐱,t|𝐱′,t′))+D∇′𝟐G(𝐱,t|𝐱′,t′))=±δ(𝐱−𝐱′)δ(t−t′)\frac{\partial G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\right)}{\partial t^{\prime}}+\mathbf{v}\cdot\nabla^{\prime}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\right)+D\boldsymbol{{\nabla^{\prime}}^{2}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime})\right)=\pm\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\delta\left(t-t^{\prime}\right) (53)

where: i) we have arbitrarily guessed that all signs on the left of the adjoint equation are positive, and ii) have fixed the solution point, (𝐱,t).\left(\mathbf{x},t\right). Choice ii): a) leads, in straightforward fashion, to an isolated term, η​(𝐱,t),\eta\left(\mathbf{x},t\right), at the end of the procedure, and b) means that the space-time integrals formed in step c) are taken over 𝐱′\mathbf{x^{\prime}} and t′.t^{\prime}.

Carrying out steps b) and c), we obtain:

∫0t+ϵ∫Ω′η​(𝐱′,t′)​[∂G​(𝐱,t|𝐱′,t′)∂t′+𝐯⋅∇′G​(𝐱,t|𝐱′,t′)+D​∇′𝟐⁡G​(𝐱,t|𝐱′,t′)]​𝑑𝐱′​𝑑t′−\displaystyle\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\left[\frac{\partial G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)}{\partial t^{\prime}}+\mathbf{v}\cdot\nabla^{\prime}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)+D\boldsymbol{{\nabla^{\prime}}^{2}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\right]\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}-
−∫0t+ϵ∫Ω′G​(𝐱,t|𝐱′,t′)​[∂η​(𝐱′,t′)∂t+𝐯⋅∇′η​(𝐱′,t′)−D​∇′2⁡η​(𝐱′,t′)]​𝑑𝐱′​𝑑t′=\displaystyle-\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\left[\frac{\partial\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)}{\partial t}+\mathbf{v}\cdot\nabla^{\prime}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)-D\boldsymbol{\nabla^{\prime}}^{2}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\right]d\mathbf{x^{\prime}}dt^{\prime}=
±∫0t+ϵ∫Ω′η​(𝐱′,t′)​δ​(𝐱−𝐱′)​δ​(𝐭−𝐭′)​𝑑𝐱′​𝑑t′−∫0t+ϵ∫Ω′G​(𝐱,t|𝐱′,t′)​f​(𝐱′,t′)​𝑑𝐱′​𝑑t′\displaystyle\pm\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\mathbf{\delta\left(\mathbf{x}-\mathbf{x^{\prime}}\right)\delta\left(t-t^{\prime}\right)}\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}-\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)f\left(\mathbf{x^{\prime}},t^{\prime}\right)\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime} (54)

where 𝐯=𝐯​(𝐱′,t′),\mathbf{v}=\mathbf{v}\left(\mathbf{x^{\prime}},t^{\prime}\right), and where the upper limit on the time integrals, t+ϵ,t+\epsilon, is set just beyond t;t; as shown in Eq. (55) below, following time integration, this trick suppresses the a priori unknown integral at t′=t+ϵ.t^{\prime}=t+\epsilon. In addition, it allows exposure of η​(𝐱,t)\eta\left(\mathbf{x},t\right) from the first right-hand term in Eq. (S1 Appendix C).

Moving to step d) and focusing on Eq. (S1 Appendix C), use Liebniz rule on the fourth term on the left and combine with the first:

∫Ω′[G​(𝐱,t|𝐱′,t′)​η​(𝐱′,t′)]0t′+ϵ​𝑑𝐱′−2​∫0t+ϵ∫Ω′η​(𝐱′,t′)​∂G​(𝐱,t|𝐱′,t′)∂t′​𝑑𝐱′​𝑑t′\int_{\Omega^{\prime}}\left[G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\right]_{0}^{t^{\prime}+\epsilon}\hskip 2.84544ptd\mathbf{x^{\prime}}-2\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\frac{\partial G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)}{\partial t^{\prime}}\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime} (55)

Using the property that G​(𝐱,t|𝐱′,t′)G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right) is undefined at for t′>t,t^{\prime}>t, the first term simplifies to −∫Ω′G​(𝐱,t|𝐱′,t′=0)​η​(𝐱′,t′=0)​𝑑𝐱′,-\int_{\Omega^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}=0\right)\eta\left(\mathbf{x^{\prime}},t^{\prime}=0\right)\hskip 2.84544ptd\mathbf{x^{\prime}}, where η​(𝐱′,t′=0)\eta\left(\mathbf{x^{\prime}},t^{\prime}=0\right) is the specified initial condition. Since the line (d=2) or area (d=3) integral in the second term of Eq. (55) must be suppressed, we see that the guessed plus sign on the first term in Eq. (53) must be changed to a minus.

Using similar steps on the two advective terms on the left side of Eq. (S1 Appendix C (terms 2 and 5) yields:

∫0t+ϵ∫Ω′∇′⋅(G​(𝐱,t|𝐱′,t′)​η​(𝐱′,t′)​𝐯)​𝑑𝐱′​𝑑t′−2​∫0t+ϵ∫Ω′G​(𝐱,t|𝐱′,t′)​∇′⋅(η​(𝐱′,t′)​𝐯)​𝑑𝐱′​𝑑t′\int_{0}^{t+\epsilon}\int_{\Omega}^{\prime}{\nabla^{\prime}}\cdot\left(G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\mathbf{v}\right)\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}-2\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right){\nabla^{\prime}}\cdot\left(\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\mathbf{v}\right)\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime} (56)

where, for incompressible flow, ∇′⋅𝐯=0,{\nabla}^{\prime}\mathbf{\cdot}\mathbf{v}=0, allowing the 𝐯\mathbf{v} in term 5, Eq. (S1 Appendix C) to be brought in front of the divergence operator: 𝐯⋅∇η​(𝐱′,t′)=∇′⋅(𝐯​η).\mathbf{v}\cdot\nabla\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)={\nabla}^{\prime}\cdot\left(\mathbf{v}\eta\right). Both area (d=2) or volume (d=3) integrals in Eq. (56) are a priori unknown. The first can be restated as a known/calculable line (d=2) or area (d=3) integral via the divergence theorem. The second, however, must again be suppressed by changing the guessed plus sign on term 2, Eq. (53), to a minus. Thus, the two terms in Eq. (56) simplify to:

∫0t+ϵ∮δ​Ωa​d​v​e​c[G​(𝐱,t|𝐱′,t′)​η​(𝐱′,t′)​𝐯]⋅𝐧^′​(𝐱′,t′)​𝑑𝐱′​𝑑t′\int_{0}^{t+\epsilon}\oint_{\delta\Omega_{advec}}\left[G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\mathbf{v}\right]\mathbf{\cdot}\mathbf{\hat{n}^{\prime}}\left(\mathbf{x^{\prime}},t^{\prime}\right)\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime} (57)

where δ​Ωa​d​v​e​c\delta\Omega_{advec} denotes the advective (inflow plus outflow) boundaries of Ω,\Omega, and where η​(𝐱′,t′)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right) and 𝐯​(𝐱′,t′)\mathbf{v}\left(\mathbf{x^{\prime}},t^{\prime}\right) are assumed known, at all times and locations on δ​Ωa​d​v​e​c.\delta\Omega_{advec}. In addition, 𝐧^′​(𝐱′,t′)\mathbf{\hat{n}^{\prime}}\left(\mathbf{x^{\prime}},t^{\prime}\right) is the outward unit normal to δ​Ωa​d​v​e​c.\delta\Omega_{advec}. Finally, focusing on terms 3 and 6 in Eq. (S1 Appendix C), we are again lead to an a priori unknown surface/volume (d=2/3) integral:

2​∫0t+ϵ∫Ω′∇′G​(𝐱,t|𝐱′,t′)⋅∇′η​(𝐱′,t′)​𝑑𝐱′​𝑑t′2\int_{0}^{t+\epsilon}\int_{\Omega^{\prime}}{\nabla^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\mathbf{\cdot}{\nabla^{\prime}}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime} (58)

In order to suppress this term, we change the plus sign on the third term in Eq. (53) to a negative. Using the divergence theorem and integration by parts, introducing the three sign changes above into the guessed adjoint equation, Eq. (53), labeling the functions associated with specified Dirichlet, Neumann, and initial conditions on η​(𝐱,t),\eta\left(\mathbf{x},t\right), respectively, as g​(𝐱,t),g\left(\mathbf{x},t\right), h​(𝐱,t),h\left(\mathbf{x},t\right), and ϕ​(𝐱),\phi\left(\mathbf{x}\right), and choosing the plus sign on the first term of the right side of Eq. (S1 Appendix C), we finally obtain the magic rule:

η​(𝐱,t)=∫0t∫Ω′G​(𝐱,t|𝐱′,t′)​f​(𝐱′,t′)​𝑑𝐱′​𝑑t′+∫Ω′G​(𝐱,t|𝐱′,t′=0)​η​(𝐱′,t′=0)​𝑑𝐱′−\displaystyle\eta\left(\mathbf{x},t\right)=\int_{0}^{t}\int_{\Omega^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)f\left(\mathbf{x^{\prime}},t^{\prime}\right)\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}+\int_{\Omega^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}=0\right)\eta\left(\mathbf{x^{\prime}},t^{\prime}=0\right)\hskip 2.84544ptd\mathbf{x^{\prime}}-\hskip 5.69046pt
−D​∫0t∫ΩD′h​(𝐱′,t′)​∇′G​(𝐱,t|𝐱′,t′)⋅𝐧′​𝑑𝐱′​𝑑t′+D​∫0t∫ΩN′G​(𝐱,t|𝐱′,t′)​∇′η​(𝐱′,t′)⋅𝐧′​𝑑𝐱′​𝑑t′−\displaystyle-D\int_{0}^{t}\int_{\Omega_{D}^{\prime}}h\left(\mathbf{x^{\prime}},t^{\prime}\right)\boldsymbol{{\nabla^{\prime}}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\cdot\mathbf{n^{\prime}}\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}+D\int_{0}^{t}\int_{\Omega_{N}^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\boldsymbol{{\nabla^{\prime}}}\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\cdot\mathbf{n^{\prime}}\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}-\hskip 2.84544pt
−∫0t∫ΩF′G​(𝐱,t|𝐱′,t′)​η​(𝐱′,t′)​𝐯​(𝐱′,t′)⋅𝐧′​𝑑𝐱′​𝑑t′\displaystyle-\int_{0}^{t}\int_{\Omega_{F}^{\prime}}G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right)\mathbf{v}\left(\mathbf{x^{\prime}},t^{\prime}\right)\cdot\mathbf{n^{\prime}}\hskip 2.84544ptd\mathbf{x^{\prime}}dt^{\prime}\hskip 17.07182pt (59)

where ΩD′,\Omega^{\prime}_{D}, ΩN′,\Omega^{\prime}_{N}, and ΩF′\Omega^{\prime}_{F} correspond, respectively, to the portions of the domain boundary, δ​Ω′,\delta\Omega^{\prime}, on which Dirichlet, Neumann, and inflow/outflow boundary conditions are imposed, and where h​(𝐱′,t′)h\left(\mathbf{x^{\prime}},t^{\prime}\right) is the spatially and temporally specified value of η​(𝐱′,t′)\eta\left(\mathbf{x^{\prime}},t^{\prime}\right) on the Dirichlet boundary.


S1 Appendix D

Extraction of the transport equation, Eq. (37), associated with the Green’s function in test 3. The transition density estimated in test 3, governed by the Fokker-Planck equation, Eq. (34), also satisfies the backward Kolmogorov equation:

𝐃′:∇′T⊗∇′p+𝐛⋅∇′p+∂p∂t′=0,\mathbf{D^{\prime}}\mathbf{:}{\nabla^{\prime}}^{T}\otimes\nabla^{\prime}p+\mathbf{b}\cdot\nabla^{\prime}p+\frac{\partial p}{\partial t^{\prime}}=0, (60)

where 𝐃′=𝐃′​(𝐱′,t′),\mathbf{D^{\prime}}=\mathbf{D^{\prime}}\left(\mathbf{x^{\prime}},t^{\prime}\right), ∇′,\nabla^{\prime}, and 𝐛=𝐛​(𝐱′,t′)=−𝐯​(𝐱′,t′),\mathbf{b}=\mathbf{b}\left(\mathbf{x^{\prime}},t^{\prime}\right)=-\mathbf{v}\left(\mathbf{x^{\prime}},t^{\prime}\right), are expressed in terms of variable ’backward’ variables, 𝐱′\mathbf{x^{\prime}} and t′.t^{\prime}.

The corresponding adjoint equation is given by:

𝐃′:∇′T⊗∇′K+𝐛⋅∇′K+∂K∂t′=0,\mathbf{D^{\prime}}\mathbf{:}{\nabla^{\prime}}^{T}\otimes\nabla^{\prime}K+\mathbf{b}\cdot\nabla^{\prime}K+\frac{\partial K}{\partial t^{\prime}}=0, (61)

where, again, the Green’s function is given by:

G​(𝐱,t|𝐱′,t′)=H​(t−t′)​K​(𝐱,t|𝐱′,t′)G\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right)=H\left(t-t^{\prime}\right)K\left(\mathbf{x},t|\mathbf{x^{\prime}},t^{\prime}\right) (62)

Given Eq. (61), the associated transport equation is given by:

𝐃:∇T⊗∇η−𝐛⋅∇η−∂η∂t=−f​(𝐱,t),\mathbf{D}\mathbf{:}{\nabla}^{T}\otimes\nabla\eta-\mathbf{b}\cdot\nabla\eta-\frac{\partial\eta}{\partial t}=-f(\mathbf{x},t), (63)

where, in Eq. (37), we write 𝐛\mathbf{b} as 𝐯𝐚.\mathbf{v_{a}}.