A high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangian drift-diffusion models coupled with Eulerian discontinuous spectral element method
AA HIGH - ORDER SEMI -L AGRANGIAN METHOD FOR THECONSISTENT M ONTE -C ARLO SOLUTION OF STOCHASTIC L AGRANGIAN DRIFT - DIFFUSION MODELS COUPLED WITH E ULERIAN DISCONTINUOUS SPECTRAL ELEMENT METHOD
A P
REPRINT
Natarajan, H. , Popov, P.P. , Jacobs, G.B. Department of Aerospace Engineering, San Diego State University, San Diego, USA* Corresponding author email address: [email protected] 17, 2020 A BSTRACT
The explicit semi-Lagrangian method method for solution of Lagrangian transport equations asdeveloped in [Natarajan and Jacobs,
Computer and Fluids , 2020] is adopted for the solution ofstochastic differential equations that is consistent with Discontinuous Spectral Element Method(DSEM) approximations of Eulerian conservation laws. The method extends the favorable prop-erties of DSEM that include its high-order accuracy, its local and boundary fitted properties andits high performance on parallel platforms for the concurrent Monte-Carlo, semi-Lagrangian andEulerian solution of a class of time-dependent problems that can be described by coupled Eulerian-Lagrangian formulations. Such formulations include the probabilistic models used for the simulationof chemically reacting turbulent flows or particle-laden flows. Consistent with an explicit, DSEMdiscretization, the semi-Lagrangian method seeds particles at Gauss quadrature collocation nodeswithin a spectral element. The particles are integrated explicitly in time according to a drift ve-locity and a Wiener increment forcing and form the nodal basis for an advected interpolant. Thisinterpolant is mapped back in a semi-Lagrangian fashion to the Gauss quadrature points througha least squares fit using constraints for element boundary values. Stochastic Monte-Carlo samplesare averaged element-wise on the quadrature nodes. The stable explicit time step Wiener incrementis sufficiently small to prevent particles from leaving the element’s bounds. The semi-Lagrangianmethod is hence local and parallel and does not have the grid complexity, and parallelization chal-lenges of the commonly used Lagrangian particle solvers in particle-mesh methods for solution ofEulerian-Lagrangian formulations. Formal proof is presented that the semi-Lagrangian algorithmevolves the solution according to the Eulerian Fokker-Planck equation. Numerical tests in one andtwo dimensions for drift-diffusion problems show that the method converges exponentially for con-stant and non-constant advection and diffusion velocities. K eywords semi-Lagrangian, Eulerian-Lagrangian, stochastic differential equation, discontinuous spectral elementmethod Chaotic dynamics govern the behavior of a range of physics, such as turbulent flows, molecular dynamics, and plasmas.The mathematical modeling of such stochastic physics requires a formulation based on a multi-dimensional probabilitydensity function (PDF). Molecular diffusion processes, for example, are well-known to be described by the MaxwellianPDF in phase space. In turbulence modeling the Fokker-Planck (FP) equations govern the PDF of sub-grid velocity a r X i v : . [ phy s i c s . c o m p - ph ] S e p high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT fluctuations and correlations thereof [1]. Similarly, in chemically reacting turbulent flow the FP model governs theprobability density function that is dependent on the number of species involved in the chemical reaction [1, 2]. Yetanother example is the Vlasov model that describes the motion of charged particles in phase space (See for example[3]).The dimensionality of the PDF in stochastic models is usually high. In a chemical reaction, for example, the numberof species, i.e. the dimension d of the PDF, can easily be on the order of . If we use N degrees of freedom toapproximate the stochastic partial differential equations that governs the PDF, then the degrees of freedom required forsolution is on the order of N d . This can easily yield problem sizes beyond the limitations of modern day computationalinfrastructure, even for simple problems with a relatively low number of spatial dimensions.In order to overcome this so-called "curse of dimensionality", the equations that govern the model are usually notsolved directly. Rather, a Monte-Carlo approach is used that provides samples from which the PDF can be constructed.The Fokker-Planck equation for chemical species, for example, is commonly solved using an equivalent model basedon a stochastic differential equation (SDE) [4]. In this approach, fictitious Monte Carlo(MC) particle tracers that carrythe species’ information are advected in physical space according to the SDE. At any time, the spatially dependentPDF is recovered using averaging techniques on the MC realizations.In a similar fashion, it is well known that the diffusion equation can be solved with Monte-Carlo techniques basedon random walk models and stochastic Wiener processes. In grid based random walk (RW) methods [5], fictitiousparticles that represent the concentration field are seeded at equidistant grid points. The grid points are spaced by adistance of √ D ∆ t , where D is the diffusion coefficient and ∆ t , a time increment. The random walking particles jumpto a neighboring grid point with equal probabilities. For diffusion processes, this method can be proven equivalent tothe second-order central finite difference approximation of the second-order differential terms in diffusion equations[5].Strong RW methods fall into the broader class of particle methods [6] that do not depend on an underlying mesh. Inthe strong RW method the tracers are randomly initialized. They move a distance of √ D ∆ t over a time ∆ t accordingto a Wiener process, dW t , that is defined by an independent random number selected from a normal distribution with amean of zero. At any given time and point in space, a PDF can be determined from the tracers through binning and/ordistribution of the tracer’s influence using distribution functions. Because the Monte-Carlo method is well-known toconvergence according to the inverse of the square root of the number samples, a large number of samples is required.To achieve an error on the order of − , for example, one million samples are required at each point in space. Inpractical simulations reported in literature, the number of samples is usually much smaller, and the sampling error isis of engineering accuracy within a few percent.To reduce the computational cost and improve accuracy, Ref. [7] proposed the so-called "global random walk"(GRW),which is a modification of the weak RW method. In GRW a share of the tracer particles are not moved and theremaining share is scattered to the neighboring nodes according to a Bernoulli distribution. This reduces the numberof required Monte-Carlo realizations since the particles are distributed according to a single random number. TheGRW method generalizes to a finite difference method for diffusion processes and is generally limited to low-orderaccuracy in space.In many physics models the stochastic Lagrangian model couples to a system of Eulerian partial differential equationsthat governs the field dynamics for the particle tracer. In turbulence modeling, for example, the stochastic tracer thatmodels the PDF of subgrid turbulence stresses is coupled to an averaged or filtered flow model, i.e. the filtered oraveraged Navier-Stokes equations. In plasmas, the Maxwell equations govern the electric and magnetic fields thatforce the stochastic motion of charged particles and vice-versa.High-order accurate schemes like discontinuous spectral element methods (DSEM) [8, 9] are a particularly goodchoice to solve these time-dependent Eulerian equations. Because of their low dispersion and diffusion errors DSEMsare generally better at propagating waves over longer distances and they capture small scales with fewer degreesof freedom as compared to low-order methods. Moreover, DSEM approximates the governing Eulerian equations onunstructured grids of quadrilateral or hexahedral elements which allows for the simulation of complex geometry. Sincethe method is local, i.e. no overlap between elements, DSEM is highly parallel. DSEM Navier-Stokes solvers havebeen shown extensively to obtain high accuracy and convergence using unstructured grids on complex geometries[10, 11]. Moreover, both in theory [10] and in testing through benchmarks (e.g. [8, 12, 13]), DSEMs have beenshown to have superior computational efficiency and parallelism for computation of smooth flows as compared tomore traditional discretization methods.Because of the dynamic nature of the tracers, the consistent high-order coupling of (Monte-Carlo) tracer particlesto the DSEM framework is challenging and computationally expensive. Several studies report on the coupling ofthe stochastic tracers to a DSEM field solver. In Refs. [9, 14, 15, 16, 17, 18, 19], consistent interpolation methods2 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT were developed. Several high-order distribution functions were proposed to distribute the particle influence on to theEulerian grid. In Ref. [20] linear distributions functions are used that are local to an element. Refs. [21, 22] couplesan ensemble average solution determined on hexahedral domains to the unstructured Eulerian DSEM solver. In allof these approaches, either accuracy and/or the locality of the method is compromised, which is detrimental to thecomputational efficient solution of the model.As an alternative to Lagrangian particles tracers, we introduced an explicit high order semi-Lagrangian (SL) methodfor the solution of deterministic transport equations in Ref. [23]. The SL method solves the Lagrangian form of thetransport equations and can be used instead of the particle solver in Eulerian-Lagrangian formulations. By seedingparticles on the DSEM solver’s quadrature nodes within a spectral element, the connection between the particle andfield solver is direct and consistently high-order accurate. Particles are integrated one time step forward along theircharacteristic path. The time step is restricted such that particles do not cross the element boundaries. The advectedparticle solution is remapped to the collocation nodes using a least-squares fit with boundary and mass conservationconstraints. The SL method thus remains local and does not require additional attention for parallel computing otherthan at the elements interfaces when the advected solution is patched at interfaces.In this paper, we adopt and test the semi-Lagrangian method for the Monte-Carlo solution of probability densityfunction equations and diffusion equations with a stochastic differential equation. Monte-Carlo tracers are trackedstochastically with the semi-Lagrangian method developed in [23] and are sampled at the Gauss quadrature points.There is hence no need for a binning or a distribution method. Because the method is local, the SL approach ensureshigh parallel efficiency and high order accurate boundary condition implementation that has eluded and plagued SDEmethods coupled with high-order field solvers thus far. The resulting approach shares some similarities with theEulerian Monte Carlo method [24, 25], but also has distinct advantages over that approach.The paper is organized as follows. First, generic stochastic differential equations and their equivalent Eulerian formsare discussed. Next, a staggered-grid discontinuous spectral element method for the solution of an Eulerian fieldis briefly summarized. Before introducing the semi-Lagrangian algorithm for solution of the stochastic differentialequation, we review for reference some common random walk methods. Tests are conducted for one dimensional,constant and non-constant diffusion problems as well as for stochastic problems. Finally, the SL solver is coupledwith a Navier-Stokes solver for the solution of a species equation in a temporally developing shear layer. Conclusionsand future steps are reserved for the final section.
We consider the canonical stochastic differential equation in the Itô sense for transport of MC particles in the physicalspace, X t : d X t = u ( X t , t ) dt + σ ( X t , t ) d W t . (1)This SDE is central to many models for a range of problem as described in the introduction, including the filtered massdensity function model for the modeling species transport in subgrid turbulent scales, which is the broader focus ofour research [1]. It is usually complimented by a transport equation for a variable in compositional space that we donot consider here. In (1), the drift velocity, u , interpolates from a coupled field solver. For the problems we consider,these are usually the Navier-Stokes equations. In most of the test cases below, we will assume a prescribed u field.The diffusion is driven by a Wiener process W t with diffusion coefficient, σ .The probability density, P ( x , t ) can be recovered from Monte-Carlo (MC) realizations of the Lagrangian stochasticdifferential equation at any given point in space and time. This MC solution is well-known to be equivalent to solvingthe Eulerian, Fokker-Planck equation for P ( x , t ) (e.g. [4]) given by ∂ P ( x , t ) ∂t + ∂ ( u j ( x , t ) P ( x , t )) ∂x j = ∂ ( D ( x , t ) P ( x , t )) ∂x j ∂x j , (2)with D = σ / .For a constant diffusion coefficient D ( x , t ) = D c = constant , we can rewrite (2) as ∂φ ( x , t ) ∂t + u j ( x , t ) ∂φ ( x , t ) ∂x j = D c ∂ φ ( x , t ) ∂x j ∂x j − φ ( x , t ) ∂u j ( x , t ) ∂x j , (3)where we have used the notation φ ( x, t ) = P ( x, t ) to be consistent with the deterministic formulation and method asdescribed [23]. For a conservative medium with a divergence free velocity field, the equation further reduces to thegeneric convection diffusion equation ∂φ ( x , t ) ∂t + u j ( x , t ) ∂φ ( x , t ) ∂x j = D c ∂ φ ( x , t ) ∂x j ∂x j . (4)3 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
We use the many known analytical solutions for this equation to assess the accuracy of the semi-Lagrangian methodfor solution of the SDE below.The new procedure introduced in this work is capable of solving a wide class of Fokker-Planck problems, including (2).Specifically, we solve for a random field, φ ∗ ( x , t ) , whose PDF is denoted by P φ ( ψ ; x , t ) (where ψ is the sample spacevariable of φ ∗ ). The semi-Lagrangian scheme then leads to the following Fokker-Planck equation for P φ ( ψ ; x , t ) : ∂ P φ ∂t + u j ∂ P φ ∂x j = ∂∂x j (cid:18) D ∂ P φ ∂x j (cid:19) − ∂∂ψ [ S P φ ] , (5)where S ( ψ ; x , t ) is a deterministic source term defined over the sample space. Note that while (5) is given for aone-dimensional random variable φ ∗ , the extension to multi-dimensional φ ∗ is trivial. In the present work we set S ( ψ ; x , t ) ≡ − ψ ∂u j ∂x j for the purpose of recovering (3) for the case when D = D c = constant .Setting φ ( x , t ) = (cid:82) ψ P φ ( ψ ; x , t ) dψ to be the mean of φ ∗ ( x , t ) at a specific point ( x , t ) , it can be easily shown thattaking the first moment, (cid:82) ψ · · · dψ , of each term in (5) leads to (3). This allows us to compare our method withmethods such as global, strong and weak RW, whose Fokker-Planck equation is (2).While some Monte Carlo solvers aim to find a solution to (2), in many applications (5) is just as useful as a startingpoint, and there is no need to go through (2). As an example, in the large eddy simulation/filtered mass density(LES/FMDF) method of Jaberi et al.[1], in the limit as the filter size goes down to , the FMDF transport equation(eq.29 of [1]) for a one-dimensional compositional variable becomes equivalent to ∂ (cid:16) (cid:104) ρ (cid:105) (cid:102) P φ (cid:17) ∂t + ∂∂x j (cid:16) (cid:104) ρ (cid:105) (cid:101) u j (cid:102) P φ (cid:17) = ∂∂x j (cid:32) Γ ∂ (cid:102) P φ ∂x j (cid:33) − ∂∂ψ (cid:104) (cid:104) ρ (cid:105) S (cid:102) P φ (cid:105) , (6)where (cid:102) P φ is the Favre (i.e., density-weighted) PDF of φ ∗ ( x , t ) , (cid:104) ρ (cid:105) and (cid:101) u j are the mean density and Favre-averagedvelocity, Γ is the combination of turbulent and molecular diffusivity, and the source term S combines the effects ofthe mixing model and chemical reaction. In Jaberi et al., the authors perform a Monte Carlo solution of (6) using fullyLagrangian particles which evolve by an SDE. The drift term of this SDE’s spatial component contains a gradientof the diffusivity so as to make the Fokker-Planck equation for the particle system (essentially the multi-dimensional,anisotropic version of (2)) equivalent to (6). Alternatively, (6) can be recovered from (5) by setting u j = (cid:101) u j − (cid:104) ρ (cid:105) ∂ (cid:104) ρ (cid:105) ∂x j Γ and setting D = Γ (cid:104) ρ (cid:105) . With these definitions of u j and D , (5) can be multiplied through by (cid:104) ρ (cid:105) and (6) follows, providedthe density consistency condition, ∂ (cid:104) ρ (cid:105) ∂t + ∂∂x j ( (cid:104) ρ (cid:105) (cid:101) u j ) , (7)is satisfied, meaning that the definition of ρ ( ψ ) must be such that the mean density satisfies the averaged continuityequation. The need to satisfy density consistency does not make the present method any more cumbersome thanLagrangian particle methods, which have their own density consistency conditions [26, 27] that they must satisfy.Therefore, while (2) yields good test cases of the semi-Lagrangian method, the Fokker-Planck equation (5), which thesemi-Lagrangian method solves naturally is just as useful for modeling probabilistic systems. Following [23], the semi-Lagrangian (SL) method is consistently coupled to the staggered grid discontinuous spectralelement method (DSEM) as first introduced by Kopriva [8]. In this version of DSEM the solution variable is collocatedat Gauss quadrature nodes and the fluxes on Lobatto quadrature nodes. The collocation at Gauss quadrature nodes arespecifically beneficial for the simple and consistent coupling between the SL method and DSEM because it it leads topreservation of the high-order, local nature of DSEM as we showed in [23]. For completeness, we briefly summarizeessential aspects of the staggered grid DSEM method again. For a detailed description, we refer to [11, 12, 23].In DSEM, the physical domain Ω is divided into K non-overlapping elements, Ω = ∪ Kk =1 Ω k . In the context of DSEM,elements are often referred to as subdomains, a nomenclature that we follow in this paper. Each physical subdomain is4 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT then mapped onto a unit computational cube using iso-parametric transformation [9]. The governing Eulerian equationis given by, ∂ ˜ Q ∂t + ˜ ∇ · ˜ F = 0 , (8)where, ˜ Q = | J | Q , ˜ ∇ · ˜ F = ∂ ˜ f∂ξ + ∂ ˜ g∂η + ∂ ˜ h∂ζ . | J | is the determinant of the transformation from the physical to thecomputational domain.The solution and flux collocation points are chosen according to Chebyshev Gauss and Lobatto quadrature points,which along tensorial grid lines, ≤ ξ ≤ , are given by, ξ i +1 / = 12 (cid:20) − cos (cid:18) i + 1 / N + 1 (cid:19) π (cid:21) i = 0 , , ..., N − , (9)and ξ i = 12 (cid:20) − cos (cid:18) iπN (cid:19)(cid:21) i = 0 , , ..., N, (10)respectively. Here, we have used the integer subscript, i , to identify Lobatto points and i + 1 / to identify Gauss pointsthat are located in between two Lobatto points i and i + 1 . In three dimensions, the solution interpolant ˜ Q is then ˜ Q ( ξ, η, ζ ) = N − (cid:88) i =0 N − (cid:88) j =0 N − (cid:88) k =0 ˜ Q i +1 / ,j +1 / ,k +1 / h i +1 / ( ξ ) h j +1 / ( η ) h k +1 / ( ζ ) , (11)where h i +1 / ( ξ ) is the Lagrange interpolation polynomial of degree N-1 defined on the Gauss quadrature points ξ m +1 / and h i +1 / ( ξ ) = N − (cid:89) m =0 m (cid:54) = p ξ − ξ m +1 / ξ i +1 / − ξ m +1 / , i = 0 , , ..., N − , (12)is the Lagrangian polynomial of degree N -1. The fluxes, ˜ F , are collocated similarly on the Lobatto points. Throughinterpolation between the Gauss grid and the Lobatto grid, the fluxes can be determined as a function of the solution, ˜ Q . Through an approximated Riemann solver, an interface flux is determined from interface solutions on neighbouringsubdomains. The derivatives of the fluxes, ˜ ∇ · ˜ F , are determined at the Gauss points. Then, it remains to update theGauss solution in time. We typically use an explicit integrator such as a standard fourth order explicit Runge-Kuttatime stepping method. Before we present the semi-Lagrangian methods based on DSEM, we review some of the most common random walkmethods in one-dimension which we will use for comparison and reference to the SL-DSEM.
In the strong random walk method, the spatial location x pj of N p Monte Carlo tracers are advected according to theSDE in (1) using the first order Euler-Maruyama method [28] as follows: x pj ( t + ∆ t ) = x pj ( t ) + ∆ t u ( x pj ( t )) + √ DdW t j = 1 ...N p . (13)Here, ∆ t is the time step. The particle’s solution in compositional space, φ pj , is advected along its characteristic path.The particles are randomly seeded and traced within a computational domain defined on the interval [ x a , x b ] .By sampling within bins (or elements) the probability density function, P ( x, t ) can be constructed. We use N b equidis-tant bins between x i and x i +1 with the center location of each bin given by with x i +1 / = i ∆ x + x a + ∆ x i = 0 ....N b − (14)and ∆ x = ( x b − x a ) / ( N b ) .The analytical average, (cid:104) φ (cid:105) ( x i +1 / ) within a bin with center location x i +1 / , is determined using the first moment of P ( x, t ) with respect to φ as (cid:104) φ (cid:105) i +1 / = (cid:90) x i +1 x i φ ( x ) P ( x i +1 / , t ) dx (15)This is equivalent to ensemble averaging of the MC realizations within a bin.5 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
The weak random walk algorithm is grid based (e.g. [5]). Starting again from a computational domain on the interval [ x a , x b ] , we define an equidistant grid with N b equidistant elements using the nodes, x i = i ∆ x + x a i = 0 ....N b , (16)where ∆ x = ( x b − x a ) / ( N b ) is the grid spacing.Total number of particles, N p is uniformly distributed over the nodes. The number density weighted solution at a givennode is initialized as m i = 12 φ i ∆ xN p i = 1 ....N. (17)with φ i the initial condition. Each particle at a given node moves to a neighboring node with equal probability fromtime t n to t n +1 . The time step size is ∆ t = t n +1 − t n , and it must be related to the grid spacing ∆ x according to, ∆ x = √ D ∆ t. (18)to consistently capture the diffusion in (4). The updated distribution of particles at time step n + 1 can be determinedas, m n +1 i = 12 ( m ni − + m ni +1 ) . (19)At the new step, φ n +1 i , is recovered as follows φ n +1 i = 2 m n +1 i ∆ xN p i = 1 ....N. (20) The GRW method introduced in [7] is similar to the random walk method. The GRW method does not move indi-vidual particles with equal probability, however. Instead it moves particles in large groups according to a prescribedprobability density function to reduce the number of samples and thus reduce computational cost. The domain andthe initial particle distribution are the same as for the random walk method given in (16) and (17), respectively. Let δm n ( j, i ) denote the density weighted solution for a group of particles at time t n moving from node x j to x i . For agiven time step only a fraction r of the number of particles move to the neighboring nodes, the rest of them determinedby δm n ( i, i ) = (1 − r ) m ni i = 1 ....N. (21)remain at the same node. The parameter r connects ∆ x and ∆ t as follows r = 2 D ∆ t (∆ x ) (22)to ensure a consistent solution of (4). Assuming the particles are moved only to the nearest neighboring node, thedistribution of particles at t n +1 for a given node x i is determined as m n +1 i = δm n ( i, i ) + δm n ( i + 1 , i ) + δm n ( i − , i ) i = 1 ....N. (23)The second and third term on the right hand side of this equation represent contributions from the the neighbour-ing nodes. The neighbouring groups move according to a Bernoulli distribution given by b m ( α ) = 2 − m C αm . Thisdistribution is sampled by a random number generator to provide α so that δm n ( i, i + 1) = α, i = 1 ....N, (24)and δm n ( i, i −
1) = m ni − δm n ( i, i ) − α, i = 1 ....N, (25)Using m n +1 i in (23), φ n +1 i is again recovered using (20). 6 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
The EMC method, first developed by Valiño [24], and refined by Sabel’nikov and Soulard [25] is another alternativefor solving Fokker-Planck equations. In contrast to the abovementioned particle methods, EMC solves for the PDFof φ ∗ by tracking a set of fields Φ ( s ) , for s = 1 , ..., N f , which are defined on the entire domain. The fields Φ ( s ) areevolved by a stochastic partial differential equation (SPDE) such as the following: d Φ ( s ) + u j ∂ Φ ( s ) ∂x j dt + ∂ Φ ( s ) ∂x j σdW ( s ) j − ∂∂x j (cid:18) D ∂ Φ ( s ) ∂x j (cid:19) dt = S (Φ ( s ) ; x , t ) dt, (26)where D = σ / and S ( φ ; x , t ) is a general source term. In the context of reactive flow simulations, this sourceterm will be a combination of the reaction source term and the effect of molecular diffusion. We note that the Wienerincrements dW ( s ) j are spatially global, i.e., the same Wiener increment sample dW ( s ) j is used for all points x in (26).Following [25], (26) leads to the Fokker-Planck equation ∂ P φ ∂t + ∂∂x j ( u j P φ ) = ∂∂x j (cid:18) D ∂ P φ ∂x j (cid:19) − ∂∂ψ ( S ( ψ ) P φ ) , (27)where P φ ( ψ ; x , t ) is the PDF of Φ( x , t ) at specified values of x and t , and ψ is the sample space variable of Φ . It iseasily seen that (27) is equivalent to (5).Like EMC, DSEM-SL determines its ensemble of smooth fields by applying the same Wiener increment to all pointsin a given sample field. As a result, in the limit of arbitrarily high spatial resolution, DSEM-SL solutions converge tothe same SPDE solved by EMC. This is formally proven in Appendix A.A downside to EMC is the appearance of the term ∂ Φ ( s ) ∂x j σdW ( s ) j of (26) in the formulation. The combinination of aderivative approximation with the Wiener increment dW ( s ) j , is particularly challenging. DSEM-SL avoids this termand has other advantages that will be discussed in the presentation of DSEM-SL in the next section.Finally, EMC methods have so far been implemented only with low-order FV [24] and ENO [25] spatial discretiza-tions, whereas DSEM-SL exhibits spectral spatial convergence. The semi-Lagrangian algorithm for simulation of the stochastic differential equation is based on the semi-Lagrangianmethod that we developed in [23] for the Monte-Carlo simulation of deterministic Lagrangian transport equations inEulerian-Lagrangian formulations. To solve stochastic models using Monte-Carlo sampling from tracers that behaveaccording to the stochastic differential equation, multiple polynomial solutions are generated with the deterministicsemi-Lagrangian method according to a Wiener process. Similar to the strong random walk method, each polynomialrealization represents a Monte-Carlo sample and can be used to reconstruct the density function at quadrature points.Below, we discuss the semi-Lagrangian method in one-dimension and highlight the implementation of stochasticcomponents. The multi-dimensional algorithm can be formulated on a tensorial grid as discussed in [23]. For brevity,we refer for details for the multi-dimensional algorithm to that article.
To be consistent with the Eulerian DSEM solver that provides the drift velocity, u , at the tracer location in (1), weinitialize N p = N particles within a subdomain, Ω k , at an initial time t at the N Gauss quadrature points in (9).The drift velocity is directly available at these quadrature points and hence does not require computational intensiveinterpolation that is necessary for general Lagrangian particle methods. A single sample, s , of the solution, φ s , at agiven time t n is approximated by a Lagrange interpolant as follows, φ sn ( ξ ) = N − (cid:88) i =0 φ sn ( ξ sni +1 / ) h i +1 / ( ξ ) , s = 1 , .., N s (28)7 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT where N s is the total number of Monte-Carlo samples used per grid point. The Lagrange polynomials, h i +1 / ( ξ ) , ofdegree N − are defined on the Chebyshev Gauss points ξ si +1 / according to (9) h i +1 / ( ξ ) = N − (cid:89) i =0 i (cid:54) = j ξ − ξ sni +1 / ξ snj +1 / − ξ sni +1 / , j = 0 , , ..., N − . (29) The particles and its associated sample polynomial solution, φ s , are advected in the physical space along its charac-teristic path according to the stochastic differential equation (1). To integrate the SDE in Itô form, we use an explicitfirst-order Euler-Maruyama scheme[28], so that in local coordinates the time step will have the form ξ s (cid:63) i +1 / = ξ sni +1 / + (cid:32) ∆ t (cid:20) u + D ∂ ξ∂x ∂x∂ξ (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) ξ sni +1 / + √ D ∆ W st (cid:33) / ∂x∂ξ . i = 0 , , ..., N − , (30)and in physical coordinates its form will be x s ∗ i +1 / = x i +1 / + u (cid:0) x i +1 / , t (cid:1) ∆ t + (cid:113) D (cid:0) x i +1 / , t (cid:1) ∆ W st . (31)Here, ∆ W st is the Wiener increment that is obtained from a random number generator according to a normal distri-bution with a mean of zero and a variance of ∆ t . As indicated by the notation ∆ W st , the same sample of the Wienerincrement is used for the advection of all spatial points which belong to the solution φ s . This is similar to what isdone in Eulerian Monte Carlo, and in contrast to a scheme such as the Lagrangian particle method, in which eachnew point gets its own sample of the Wiener increment. The advantages of the present approach are twofold: firstly,it preserves the spatial smoothness of φ s , which is required for the correct convergence of the spatial discretizationschemes. Secondly, it reduces computational effort, as much fewer calls to the random number generator are needed.The above advection formulations both yield convergence to the same result and are, in fact, identical for non-curvilinear grids. Here we use the physical coordinate formulation, (31), which is preferable because it avoids theneed to compute the higher-order metric term in (30). We note that this avoidance of metrics is a significant advantageof the semi-Lagrangian method over Eulerian Monte Carlo, as it reduces computational effort and yields a procedurewhich is better behaved on singular or close to singular grids.To obtain high-order accuracy for the time integration of the Itô form SDE one can consider Runge-Kutta methods asdiscussed in [28]. The algorithms extends naturally from the Euler-Maruyama to high-order time-integrators as wehave shown in[23]. For SDEs, however, these high-order time-integrators are increasingly complex with increasingorder and a topic of ongoing research. We have not considered them in this work, but aim to report on this in futureinvestigations.While there is no formal stability criterion for the temporal update of the linear characteristic equation, we preventan advected particle from leaving a subdomain by restricting the time step. This has two reasons. Firstly, if theadvected particle locations within a subdomain deviates only marginally from the quadrature point locations, then theVandermonde interpolating matrix can be expect to be reasonable well-conditioned ensuring that the remapping whichrequires an inversion of this matrix is not singular. Secondly, the nodes stay local to the element, which means that themethod is local and parallel and that connection of the interfaces can be performed in relatively simple manner usingthe interpolation method described below.Thus, the time step restriction is set by the following condition, | u | max ∆ t + √ D ∆ t ≤ ∆ x min (32)where, ∆ x min is the minimum grid spacing between two particles in the physical space and | u | max is the maximumadvection speed. For a pure diffusion problem without drift, this reduces to, ∆ t ≤ (∆ x min ) D , (33)which is equivalent to a Fourier number condition, but one that can be violated without loss of stability. In contrast tothe present scheme, the Wiener increment appearing in the advection term in EMC schemes yield a Fourier number8 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method
A P
REPRINT stability condition . Thus, another advantage (albeit one which is not used here) of semi-Lagrangian schemes is theability to take larger time steps.The solution after advection is denoted by φ s (cid:63) ( ξ ) and is given as, φ s (cid:63) ( ξ ) = N − (cid:88) i =0 φ (cid:63) ( ξ s (cid:63) i +1 / ) h s (cid:63) i +1 / ( ξ ) , (34)where h s (cid:63) i +1 / ( ξ ) are the Lagrange polynomials of degree N − defined on the advected points ξ s (cid:63) i +1 / , h s (cid:63) i ( ξ ) = N − (cid:89) i =0 i (cid:54) = j ξ − ξ s (cid:63) i +1 / ξ s (cid:63) j +1 / − ξ s (cid:63) i +1 / , j = 0 , , ..., N − . (35)In general the advected polynomial’s nodal solution values, φ s (cid:63) ( ξ s (cid:63) i +1 / ) is obtained by integrating φ s ( ξ si +1 / ) in thecompositional space according to (5) and the advection of φ along the flow, φ s(cid:63) ( ξ (cid:63)i +1 / ) = φ si +1 / + ∆ t (cid:32) − φ si +1 / (cid:18) ∂u∂x (cid:19) i +1 / (cid:33) . (36)Here, (cid:16) ∂u∂ξ (cid:17) i +1 / is obtained from the DSEM field solver solver. In the final remapping stage of the algorithm, the advected polynomial is projected back onto the Gauss-Chebyshevquadrature nodes through interpolation as follows, (cid:98) φ s n +1 ( ξ sj +1 / ) = N − (cid:88) j =0 φ s (cid:63) ( ξ s (cid:63) j +1 / ) h s (cid:63) j +1 / ( ξ si +1 / ) , (37)Here, we use the hat symbol to denote the intermediate solution at t n +1 . To account for connectivity between elementsand boundary conditions, we constrain this intermediate solution following [23]. Boundary conditions and interfaceconstraints are applied using interpolation. We determine the boundary values using polynomial interpolation accord-ing to (37), ˆ φ s n +1 b = N − (cid:88) j =0 φ s (cid:63) ( ξ (cid:63)j +1 / ) h (cid:63)j +1 / ( ξ b ) b = 1 , (38)By upwinding, a unique interface value is determined from the interfaces values of two neighbouring subdomains. If u ∆ t + √ DdW t is positive at the interface, then we use the information from the left element. ˆ φ s n +1 b(cid:63) = f (cid:32) ˆ φ s n +1 b =1 (cid:12)(cid:12)(cid:12)(cid:12) Ω k , ˆ φ s n +1 b =2 (cid:12)(cid:12)(cid:12)(cid:12) Ω k − (cid:33) (39)Boundary conditions are implemented in the same way as interface condition by using a specified ghost solution atcomputational domains boundaries.To project the interpolated polynomial, ˆ φ s n +1 , combined with the boundary constraints onto the Gauss-Chebyshevquadrature we use a least-squares method to solve the overdetermined system of equations as described in [23].We note here that, for a multi-dimensional φ ∗ , the majority of the computational cost of the remapping stage does notscale up with increasing dimension of φ ∗ , since each component of φ ∗ is defined on the same advected points. Thisis a significant advantage over EMC methods, for which the spatial discretization has to be applied to each componentof the random field, and thus scales linearly with the random variable’s dimension.9 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
To determine the mean of the polynomial solution, we ensemble average at the grid points only as follows: (cid:104) φ (cid:105) n +1 i +1 / = 1 N s N s (cid:88) s =1 φ s n +1 i +1 / i = 0 , , ..., N − . (40)Because the averaging is performed on the quadrature nodes, this is equivalent to averaging the polynomial on eachelement. We can also recover the PDF of the solution on each grid point by binning the samples.When using Dirichlet boundary conditions, the samples can develop high gradients at the boundaries. By re-seedingthe samples from the averaged solution once every few time steps we can reduce the high gradient. In this paper weperform re-seeding after every 100 time steps. φ s n +1 i +1 / = (cid:104) φ (cid:105) n +1 i +1 / i = 0 , , ..., N − , s = 1 , , ..., N s . (41) It can be proven that DSEM-SL is consistent, i.e. it solves the Fokker-Planck equation (5) implicitly for the PDF of φ ∗ . This proof and the derivation of the equivalent stochastic PDE is presented in Appendix A.The accuracy of the method depends on three known approximation errors that include, (1) the number of samples, (2)the accuracy of the spatial approximation and (3) time integration accuracy. The sampling error converges accordingto the inverse of the square root of the number of samples. The spatial approximation is spectrally accurate accordingto the high-order approximation in each element. We use a Euler-Maruyama time integration method, which is firstorder and hence the time integration accuracy is of O (∆ t ) . In the test cases below, we fix the ∆ t value and study theconvergence of interpolation and the sampling errors. We confirm that the DSEM-SL converges according to theseexpected error estimates. We assess the error and the behavior of the semi-Lagrangian method for several one and two dimensional test cases.Results are compared to the analytical solutions and solutions obtained with classic random walk method and EulerianMonte Carlo method as described above.Accuracy is measured using the L norm of the solution error which is calculated by summing up local L error normsin each subdomain, k , as, (cid:107) e (cid:107) L = K (cid:88) k =1 (cid:115)(cid:90) Ω k ( φ − φ exact ) J dξ, (42)where J is the Jacobian for the transformation from the physical space to the computational space. We also assesconservation properties of the method are by inspecting the following global mass and energy norms, (cid:107) M (cid:107) = K (cid:88) k =1 (cid:82) Ω k φdξ (cid:82) Ω k φ exact dξ , (43)and (cid:107) E (cid:107) = K (cid:88) k =1 (cid:82) Ω k ( φ ) dξ (cid:82) Ω k φ exact dξ . (44)respectively. As a first test, we consider the diffusion of a sine wave with the drift velocity u set to zero and the diffusion coefficientset to D c =1 according to (4). In the domain x = [0 , the initial condition is set to φ ( x, = sin(2 πx ) + 2 . Periodicboundary conditions are specified. In order to keep the time integration error low and to satisfy the stability criterion(33), a time step of ∆ t = − is used. Simulations are carried out for 50 time steps for different number of elementsizes, H , different polynomial orders, P = N − and different number of samples, N s .10 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT (a) (b)(c) (d)Figure 1: 1D diffusion of a sine wave using DSEM-SL method using one element, N s = samples and ∆ t = − . (a)shows the average solution after 50 time steps along with 10 samples when P =8 and H =1. The L error norm, (cid:107) e (cid:107) L ,the mass norm, (cid:107) M (cid:107) , and the energy norm, (cid:107) E (cid:107) are plotted versus time, t , in subfigures (b), (c) and (d), respectively.Figure 1 compares the time evolution of the error and conservation norms of the DSEM-SL scheme for different poly-nomial orders keeping the number of elements and the number of samples fixed with H =1 and N s = respectively.The time evolution of the L error shows that the error decreases as the polynomial order is increased consistent withexponential convergence in P for even and odd order polynomials separately. The difference in error between oddand even polynomial approximation is a result of the symmetry of the sine function, which favors the even numberof interpolating points for polynomials of an odd degree. At a polynomial order P =8 the interpolation error is of thesame order as the sampling error. The method accurately conserves mass for upto six decimals for all polynomialorders and the mass conservation improves as the polynomial order increases. The energy norm evolution shows thatfor low polynomial orders there is a small loss in energy, which reduces for increasing polynomial orders.Figure 2a illustrates the effect of the number of samples on the P convergence. The sampling error, which can beexpected to be on the order of N − / s , is found to be approximately O (10 − ) for N s = and is similar to the spatialapproximation error for P = 8. With a reduced number of samples N s the error at P =8 increases consistently accordingto the sampling convergence. This is confirmed by the linear trend of the error versus log ( N − / s ) in Figure 2b.A log-log plot of the error versus the grid spacing h = 1 /H in Figure 3 using samples is linear and show thatthe methods converges in an algebraic manner according to O ( h p ) . It can be observed that the error convergenceplot when P =2 and P =4 shows anomaly for three elements H =3. The local error is plotted in Figure 4 shows thatfor three elements and even polynomial order, a quadrature point is located exactly at the center location of thesine wave. The approximation in the center element then yields overshoots at the edges of the center elements andthe erratic convergence behavior. To avoid this behaviour, the h -convergence was computed after shifting the sin function. A log-log plot of the error versus the grid spacing is shown in Figure 5 using an initial condition of φ ( x, = sin(2 π ( x − . . The plot shows the expected algebraic convergence trend without anomalies.11 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT (a) (b)Figure 2: 1D diffusion of a sine wave using DSEM-SL method using one element and ∆ t = − . (a) plots the P convergence of the L error for different number of samples, N s . (b) plots the N s convergence of the error when P =8.Figure 3: 1D diffusion of a sine wave using DSEM-SL method using samples and ∆ t = − . Plot of the H convergence of the L error for different polynomial orders.A computation that employs a Dirichlet boundary condition shows no discernible differences with one that uses peri-odic boundary conditions (Figure 6). In practical simulation of more complex problems over longer times, the computational burden to generate N s = samples is too high for current day computational resources. Typically in engineering computations fewer samples areused per point on the order of tens to hundreds, yielding sampling errors of a few percent. In Figure 7 (a), we illustratethe performance of the DSEM-SL method of the diffused sine wave generated with a hundred samples, N s =100. Overa time span of t =1e-2, the amplitude of the sine wave has reduced significantly, a measure for the diffusion. Per theexpectation and comparable to Lagrangian methods, the semi-Lagrangian solution is in good comparison with theanalytical solution within a few percentages accuracy. Figure 7 (b) shows that when Dirichlet boundary conditions areused in longer time simulations, some samples develop high gradient at the boundaries. To prevent the samples fromdeveloping high gradients at the boundaries, we need to re-seed the samples from the average solution after every fewtime steps. In this example we re-seed the samples after every 100 time steps. To assess the performance of DSEM-SL in relation to exisiting methods, we compare it to the strong random walk(RW) method, weak random walk method and the generalized random walk (GRW) method as described in Section 4.For constant diffusion of a sine function, we focus on spatial accuracy and its convergence. To do so, we ensure that12 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method
A P
REPRINT (a) (b)(c) (d)Figure 4: 1D diffusion of a sine wave using DSEM-SL method using samples and ∆ t = − . Plot of the localerror for different number of elements. (a) P =1, (b) P =2, (c) P =3 and (d) P =4.Figure 5: 1D diffusion of a sine wave using DSEM-SL method using samples and ∆ t = − . Initial conditionused φ ( x, = sin(2 π ( x − . . Plot of the H convergence of the L error for different polynomial orders.13 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
Figure 6: Comparing of the log-linear plot of the L error versus P for a one dimensional diffusion of a sine waveusing DSEM-SL method with periodic and Dirichlet boundary conditions. One element is considered and ∆ t = − (a) (b)Figure 7: a) shows the plot comparing the DSEM-SL method with the analytic solution for a one dimensional Sinefunction. DSEM-SL method is run using H =3, P =4 using N s =100 samples and the the averaged solution is plottedagainst the analytic solution for t =0 and t =1e-2. b) shows that when using Dirichlet boundary condition, some samplesdevelop high gradient at the boundaries after a long period of time.the sampling error and the time integration error is kept low by using samples and a time step of − . Becauseof the large sample rate, the simulations are computationally intensive and we compute 50 time steps only. The shortsimulation times lead to a lower limit on the number of grid points for the weak RW and GRW methods to N =31 since ∆ x ≤ √ t .After 50 time steps, the number of grid points N is plotted versus the (cid:107) e (cid:107) L error norm in Figure 8. The DSEM-SL method shows exponential convergence whereas the strong RW, weak RW and the GRW methods have algebraicconvergence only. The DSEM-SL method method can achieve an error of around − using up to five times fewernumber of points compared to the GRW method. Because the diffusion of the sine wave displayed some odd convergence behaviors that are directly related to sym-metries in the polynomial point distribution and sine function behavior, we test another pure diffusion case with thediffusion coefficient, D =1 for a different initial condition. In a domain x = [ − , , we set the initial condition as aGaussian function according to an analytical solution of (2) as φ ( x, t ) = 1 √ πt exp (cid:18) − x t (cid:19) , (45)14 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
Figure 8: 1D diffusion of a sine wave comparing DSEM-SL method, strong RW method, Weak RW method and GRWmethod using samples. (a) plots the L error vs the number of grid points, N .at t = t = . . Dirichlet boundary conditions are specified according to the analytical solution. The time step is setto ∆ t = − and simulations are carried out for 100 time steps.Figure 9 compares the time evolution of the conservation properties and the error of the DSEM-SL scheme for dif-ferent polynomial orders keeping the number of elements and the number of samples fixed with H =1 and N s = respectively. Figure 9b shows that the (cid:107) e (cid:107) L error decreases as the polynomial order increases. Figure 9c and dindicate mass and energy conservation upto four decimal places. The plot of the P convergence of the L error fordifferent number of samples, N s is shown in Figure 10a. The L error at P =12 for different N s values is plotted inFigure 10b. Similar to the 1D sine wave test case, the P convergence curves follow the expected spectral convergenceuntil the polynomial interpolation error is of the same order of the sampling error and the convergence of the samplingerror follows the expected trend with a slope of / √ N s . The difference however with the sine wave test case is that theerror behavior vs even/odd polynomial orders for H = 1 is not noticeable anymore, which confirms that this behavioris specific to the sine case. If an odd number of points is used to approximate the sine function, then the middle pointalways has the exact value. This is not the case for the Gaussian initial condition. The h -convergence in the (cid:107) e (cid:107) L error is shown in Figure 11. The plot shows an expected algebraic convergence in the error as the number of elementsis increased from H =1 to H =5 using N s = samples for polynomial orders P =4 and P =5. To test the scheme for formulations that involve both advection and diffusion physics, we consider the analyticalOrnstein–Uhlenbeck solution [29] for the probability density function, P ( x, t ) , in the Fokker-Planck equation (2).The analytical solution for the Ornstein-Uhlenbeck process with u = − αx is given by, P ( x, t ) = (cid:114) α πD (1 − exp ( − ατ )) exp (cid:20) − α D ( x − x exp ( − ατ )) (1 − exp ( − ατ )) (cid:21) , (46)where, τ = t − t , t is the initial time. In this test we take α =1, D =1, t =0.25 and x = 2 .In order to use the DSEM-SL method, we re-write the 1D version of (5) with the specific source term S ( ψ ; x, t ) = − ψ ∂u∂x , ∂ P φ ∂t + u ∂ P φ ∂x = ∂∂x (cid:18) D ∂ P φ ∂x (cid:19) + ∂∂ψ (cid:20) ψ ∂u∂x P φ (cid:21) (47)Taking the first moment (cid:82) ψ · dψ , of all terms and applying integration by parts to the last term on the RHS, we get ∂φ∂t + u ∂φ∂x = ∂∂x (cid:18) D ∂φ∂x (cid:19) − ∂u∂x φ, (48)15 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT (a) (b)(c) (d)Figure 9: 1D diffusion of a Gaussian function using DSEM-SL method using one element, N s = samples and ∆ t = − . (a) shows the averaged solution and 10 samples after 100 time steps for P =12 and H =1. A second-orderEulerian Monte Carlo solution with similar spatial resolution is shown for comparison. The L error norm, (cid:107) e (cid:107) L , themass norm, (cid:107) M (cid:107) , and the energy norm, (cid:107) E (cid:107) are plotted versus time, t , in subfigures (b), (c) and (d), respectively.(a) (b)Figure 10: 1D diffusion of a Gaussian function using DSEM-SL method using one element and ∆ t = − . (a) plotsthe P convergence of the L error for different number of samples, N s . (b) plots the N s convergence of the error when P =12. 16 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
Figure 11: 1D diffusion of a Gaussian function using DSEM-SL method using samples and ∆ t = − . Plot of the H convergence of the L error for different polynomial orders.(a) (b)Figure 12: DSEM-SL method for a one dimensional Ornstein–Uhlenbeck process. (a) plots the solution at initial time, t =0.25 and at t =1, the solution is plotted against the analytic solution and a strong RW method. DSEM-SL is run using H =1 with N s =100 samples and P =17. (b) plots the P -convergence in (cid:107) e (cid:107) L error when H =1 and N s = .which is equivalent to (3) for constant D . We solve the equivalent SDE for the particles position and the transportequation for φ ∗ along the particles’ trajectories, dX t = u ( X t , t ) dt + √ D ( X t , t ) dW t , (49) Dφ ∗ Dt = − φ ∗ ∂u∂x . (50)The simulations are initialized according to the analytical solution at t =0.25 and are carried out in a domain x ∈ [ − , using N s =100 samples and a polynomial order P =17. The samples are re-seeded every 100 time steps to prevent highgradients appearing near the boundaries.Figure 12a plots the solution of the DSEM-SL method vs the analytical solution, the strong RW method and the EMCmethod (implemented with a second-order central differencing scheme) at t =1. For the latter, we provide results at twolevels of resolution: a grid with the same spatial resolution as the DSEM-SL scheme, and a finer grid with timesmore points. The strong RW method uses the same number grid points and samples as the DSEM-SL method; N = and N s = samples. The solution using DSEM-SL method matches the analytical solution better and is smoothercompared to both the strong RW method and the EMC method with the same particle resolution. The higher resolutionEMC solution has accuracy similar to that of DSEM-SL, but requires times more points.Figure 12b shows the (cid:107) e (cid:107) L error convergence on polynomial order using H =1 and N s = . The P convergence isobserved to show spectral convergence for polynomial orders from P = to P = .17 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
Figure 13: 1D non constant diffusion of a sine function: DSEM-SL method N s = samples with P =11. Plot com-pares the DSEM-SL solution with the high order finite difference, strong RW and Eulerian Monte Carlo solutions at t = . . Next we test the DSEM-SL algorithm for a non-constant diffusion coefficient of D = x and u = 0 , for which Equation(5) can then be written as, ∂ P φ ∂t = ∂∂x (cid:18) x ∂ P φ ∂x (cid:19) . (51)To show equivalence with (2), we compare to a strong RW solution with dX t = 2 X t dt + (cid:113) X t dW t (52) dφ ∗ = 0 , (53)for which the Fokker-Planck equation is ∂ P ∂t + ∂∂x (2 x P ) = ∂ ∂x (cid:0) x P (cid:1) = ∂∂x (cid:18) x ∂ P ∂x + 2 x P (cid:19) , (54)which, after cancellation of the ∂∂x (2 x P ) term on both sides, has the same functional form (51). Note that making(51) and (54) equivalent requires using different drift terms in the DSEM-SL and strong RW procedures, due to thedifferent form of the diffusive terms in (5) and (2), respectively. Taking the first moment (cid:82) ψ · dψ , of all terms in (51)and (54), both equations yield the same PDE for φ ( x, t ) , ∂φ∂t = ∂∂x (cid:18) x ∂φ∂x (cid:19) . (55)Note that we added a conserved composition variable to the strong RW solution - this does not change the functionalform of (54), but allows us to solve for negative values of φ .The DSEM-SL method is solved using H = , P = and N s = with an initial condition, p ( x,
0) = sin(2 πx ) in adomain x ∈ [0 , . A high order finite difference method (FDM) is used as the reference solution. Figure (13) plotsthe DSEM-SL, FDM, EMC and strong RW solutions at t = . . The DSEM-SL solution matches the high-order FDMsolution. 18 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT (a) (b)(c)Figure 14: Diffusion of a two dimensional sine function: DSEM-SL method using one element, N s = samples and ∆ t = − . The L error norm, (cid:107) e (cid:107) L , the mass norm, (cid:107) M (cid:107) , and the energy norm, (cid:107) E (cid:107) are plotted versus time, t , insubfigures (a), (b) and (c), respectively. The DSEM-SL method extends naturally to multiple dimensions on tensorial grids. To test, we consider a purediffusion (drift velocity set to zero) of a tensor product of sine waves in two dimensions. The initial condition is, φ ( x, y ) = sin(2 πx ) sin(2 πy ) + 2 in a domain x = [0 , ; y = [0 , . We set the diffusion coefficient, D = . The time step isset to ∆ t = − and the simulations are carried out for 50 time steps using Dirichlet boundary conditions for differentpolynomial orders, P and different number of samples, N s . This is the two-dimensional extension of the test casedescribed in 6.1.Figure 14a plots the time evolution of (cid:107) e (cid:107) L error for H =1 and N s = samples. It shows that time evolution of theerror decreases as the polynomial order increases for P =3 till P =7 where the order of the interpolation error becomesequal to the order of the sampling error. Figure 14b and c plot the time evolution of mass and energy conservationrespectively. Conservation is satisfied up-to 3 decimal places for the cases from P =3 to P =6.Figure 15a plots the P convergence of the (cid:107) e (cid:107) L error for the number of samples ranging from, N s = to .The plot shows the expected exponential convergence in P until the (cid:107) e (cid:107) L polynomial error is of the order of thesampling error. Figure 15b shows the linear trend of the sampling error versus log ( N − / s ) which is similar to theone-dimensional case. To illustrate that the algorithm works when coupled with a Navier-Stokes solver, we consider the transport of speciesin an unstable temporally developing shear layer. To this end, we couple a two-dimensional DSEM-SL solver to a19 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method
A P
REPRINT (a) (b)Figure 15: Diffusion of a two dimensional sine wave using DSEM-SL method using one element and ∆ t = − . (a)plots the P convergence of the L error for different number of samples, N s . (b) plots the N s convergence of the errorwhen P =7.DSEM Navier-Stokes solver. The fluid flow is governed by the Navier-Stokes equations given by, ∂Q∂t + ∂F ai ∂x i − ∂F vi ∂x i = 0 , (56)where, Q = ρρu ρu ρe , F ai = ρu i ρu u i + P δ i ρu u i + P δ i ( ρe + P ) u i , F vi = σ i σ i − q i + u k σ ik . (57)The scalar transport of the species, φ is given by, ∂φ∂t + u i ∂φ∂x i = 1 Re ∂ φ∂x i (58)The initial shear layer profile for the velocity in the x -direction, u and the species φ is set according to the tangenthyperbolic function given by, u ( y ) = 12 (1 + tanh y ) + 1 , (59) φ ( y ) = 12 (1 + tanh y ) , (60)in a domain x ∈ [0 , and y ∈ [ − , . Superimposed on the initial velocity field are perturbation modes determinefrom linear-stability analysis of a free shear layer (see for example [30]). Free stream boundary conditions are appliedin the y -direction and periodic boundary conditions in the x -direction. The Reynolds number based on the velocitychange of the shear layer and the thickness of the shear layer is Re = . The parameters used for the simulations are,number of elements in the x -direction, H x = and the number of elements on the y -direction, H y = with polynomialorder, P =8 and the number of samples, N s = .Figure 16 shows the two dimensional contour plot of the transport variable φ at different time snapshots as determinedby the DSEM-SL method. Physically, the simulation is of the temporal mixing of two co-flowing species, A and B .The transport variable φ is the species concentration variable with φ = representing species A and φ = representingspecies B . The contour plot, Figure 16a at time, t = shows the initial linear instability mode developing from themixing. As the mode develops temporally, in Figures 16(b and c), we can observe the non-linear mixing and formationof coherent structures in Figure 16d. 20 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT (a) (b)(c) (d)Figure 16: Transport of species in a two-dimensional layer flow configuration. Contour plot of the transport variable, φ at different time. (a) t =1, (b) t =5, (c) t =10 and (d) t =15. A semi-Lagrangian method is developed and tested for the consistent and concurrent solution of stochastic Lagrangiandifferential equations and Eulerian governing equations approximated with discontinuous spectral element methods(DSEMs). The semi-Lagrangian Monte-Carlo approach which accounts for deterministic drift (transport) and stochas-tic diffusion through a Wiener process is proven to be equivalent to solving Eulerian Fokker-Planck type models forstochastic physics such as filtered density function models for chemically reacting flows.The semi-Lagrangian method is consistent with an explicit Eulerian solver discretized with an explicit DSEM. Byseeding tracer particles at the Gauss quadrature nodes, the Lagrangian solution is directly available at quadrature nodesof the Eulerian solver and vice-versa. In Eulerian-Lagrangian methods, this exchange of information is commonlyperformed using computationally intensive and complicated interpolation methods.Consistent with DSEM, the semi-Lagrangian method is explicit for the drift term and uses a Wiener increment foreach semi-Lagrangian Monte-Carlo sample. By choosing the explicit time step and Wiener increment appropriately,particles are prevented from leaving the element. This ensures a local and parallel method, which is natural for DSEM.Following the explicit trace, the solution is remapped to the original quadrature points using a least-squares fit. Elementbased Monte-Carlo samples are averaged after the remapping stage at quadrature points only and hence do not requirebinning and/or distribution functions with an element as is common the procedure for this type of hybrid Eulerian-Lagrangian method. For a stable method, it is necessary to update the global solution according to a single Wienerincrement per sample. Using varying Wiener increments per quadrature node and/or per element leads to instabilityin numerical tests. To prevent steepening of the solution near Dirichlet boundary conditions, the samples can beperiodically reinitialized with the average of the Monte-Carlo samples.One-dimensional and two-dimensional tests are conducted for drift-difussion in one and two dimensions, includingfor a constant and non-constant diffusion coefficient, and drift-diffusion problems. The method is shown to be ex-ponentially convergent in space if the time integration error and sampling or smaller than the spatial approximationerror. Because Monte-Carlo sampling convergence is slow, according to the inverse of the square root of the numberof samples, a significantly smaller number of samples then required for formal spatial convergence is often used. For21 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method
A P
REPRINT a low sampling rate, the semi-Lagrangian method is shown to be stable and to provide engineering accuracy on theorder of a few percent of the solution. In a final test the Lagrangian method is coupled with a DSEM based parallelNavier-Stokes solver and is shown to have optimal parallel performance and provide expected qualitative results.In current work, we are extending the coupled semi-Lagrangian/Euler solver for simulation of chemically reactingflow based a filtered density function model as introduced by Givi [1].
Acknowledgements
Funding provided by the Computational Science Research Center and AFOSR under grant number FA9550-19-1-0387is greatly appreciated.
References [1] F. A. Jaberi, P. J. Colucci, S. James, P. Givi, S. B. Pope, Filtered mass density function for large-eddy simulationof turbulent reacting flows, Journal of Fluid Mechanics 401 (1999) 85–121.[2] D. Haworth, Progress in probability density function methods for turbulent reacting flows, Progress in Energyand Combustion Science 36 (2010) 168–259.[3] C. Birdsall, A. Langdon, Plasma physics via computer simulation, The Adam Hilger series on plasma physics,McGraw-Hill, 1985.URL https://books.google.com/books?id=7TMbAQAAIAAJ [4] P. Colucci, F. Jaberi, P. Givi, S. Pope, Filtered density function for large eddy simulation of turbulent reactingflows, Physics of Fluids 10 (2) (1998) 499–515.[5] W. Ames, Numerical Methods for Partial Differential equations, Vol. 2, Academic Press, 1977.[6] A. Tompson, R. Falgout, S. Smith, W. Bosl, S. Ashby, Analysis of subsurface contaminant migration and reme-diation using high performance computing, Advances in Water Resources 22 (3) (1998) 203 – 221.[7] C. Vamo¸s, N. Suciu, H. Vereecken, Generalized random walk algorithm for the numerical modeling of complexdiffusion processes, Journal of Computational Physics 186 (2) (2003) 527 – 544.[8] D. Kopriva, A staggered-grid multidomain spectral method for the compressible Navier-Stokes equations, Journalof Computational Physics (1998).[9] G. Jacobs, D. Kopriva, F. Mashayek, Towards efficient tracking of inertial particles with high-order multidomainmethods, Journal of Computational and Applied Mathematics 206 (2007) 392–408.[10] J. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: algorithms, analysis, and applications,Springer-Verlag, Berlin, 2008.[11] D. Kopriva, Implementing Spectral Methods for Partial Differential Equations: Algorithms for Scientists andEngineers, Springer Netherlands, 2009.[12] G. Jacobs, D. Kopriva, F. Mashayek, Validation study of a multidomain spectral element code for simulation ofturbulent flows, AIAA J. 43 (6) (2004) 1256–1264.[13] Z. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. A. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert,H. Huynh, N. Kroll, G. May, P.-O. Persson, B. van Leer, M. Visbal, High-order cfd methods: Current status andperspective, International Journal for Numerical Methods in Fluids (2012) 1–42.[14] G. Jacobs, J. Hesthaven, High-order nodal discontinuous Galerkin particle-in-cell method on unstructured grids,Journal of Computational Physics 214 (2006) 96–121.[15] G. Jacobs, J. Hesthaven, Implicit-explicit time integration of a high-order particle-in-cell method with hyperbolicdivergence cleaning, Computer Physics Communications 80 (10) (2009).[16] G. Jacobs, W. Don, A high-order WENO-Z finite difference based Particle-Source-in-Cell method for computa-tion of particle-laden flows with shocks, Journal of Computational Physics. 228 (5) (2009).[17] J. Suarez, G. Jacobs, W. Don, A higher-order Dirac-delta regularization with optimal scaling in the spectralsolution of one-dimensional singular hyperbolic conservation laws, SIAM Journal of Scientific Computing 36 (4)(2014).[18] J. Suarez, G. Jacobs, Regularization of singularities in the weighted summation of Dirac-delta functions for thespectral solution of hyperbolic conservation laws, Journal of Scientific Computing 72 (3) (2017).22 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method
A P
REPRINT [19] T. Stindl, J. Neudorfer, A. Stock, M. Auweter-Kurtz, C.-D. Munz, S. Roller, R. Schneider, Comparison of cou-pling techniques in a high-order discontinuous Galerkin-based particle-in-cell solver, Journal of Physics D: Ap-plied Physics 44 (19) (2011) 194004.[20] J. Komperda, Z. Ghiasi, D. Li, A. Peyvan, F. Jaberi, F. Mashayek, A hybrid discontinuous spectral elementmethod and filtered mass density function solver for turbulent reacting flows, Numerical Heat Transfer, B 78 (1)(2020).[21] S. Sammak, M. Brazell, P. Givi, D. Mavriplis, A hybrid DG-monte carlo FDF simulator, Computers & Fluids140 (2016) 158–166. doi:https://doi.org/10.1016/j.compfluid.2016.09.003 .URL [22] S. Sammak, A. Nouri, M. Brazell, D. Mavriplis, P. Givi, Discontinuous Galerkin-Monte-Carlo solver for largeeddy simulation of compressible turbulent flows, AIAA Paper 2017-0982, American Institute of Aeronautics andAstronautics (2017).[23] H. Natarajan, G. Jacobs, An explicit semi-lagrangian, spectral method for solution of lagrangian transport equa-tions in Eulerian-Lagrangian formulations, Computers and Fluids 207 (2020).[24] L. Valiño, A field monte carlo formulation for calculating the probability density function of a single scalar in aturbulent flow, Flow, Turbul. Combust. 60 (1998) 157–172.[25] V. Sabel’nikov, O. Soulard, Rapidly decorrelating velocity-field model as a tool for solving one-point fokker-planck equations for probability density functions of turbulent reactive scalars, Phys. Rev. E 72 (2005) 016301.[26] M. Muradoglu, S. Pope, D. Caughey, The hybrid method for the pdf equations of turbulent reactive flows: con-sistency conditions and correction algorithms, Journal of Computational Physics 172 (2001) 841–878.[27] P. Popov, S. Pope, Implicit and explicit schemes for mass consistency preservation in hybrid particle/finite-volume algorithms for turbulent reactive flows, Journal of Computational Physics 257 (2014) 352–373.[28] P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer Verlag, 1992.[29] G. E. Uhlenbeck, L. S. Ornstein, On the theory of the brownian motion, Phys. Rev. 36 (1930) 823–841. doi:10.1103/PhysRev.36.823 .URL https://link.aps.org/doi/10.1103/PhysRev.36.823 [30] H. Natarajan, G. Jacobs, Study of linear and non-linear instabilities in a multiple jet flow configuration, Proceed-ings of the ASME 2017 International Mechanical Engineering Congress and Exposition (2017).[31] L. C. Evans, An Introduction to Stochastic Differential Equations, American Mathematical Society, 2013.URL
Appendix A: Proof of consistency
This Appendix presents a proof of consistency of the DSEM-SL algorithm introduced in section 5, i.e it is shown thatin the limit as N → ∞ and ∆ t ↓ , φ s ( x , t ) evolves by a stochastic PDE (5).We start by considering the s − th sample, and a particle which is initially located at x . The position of the particleafter the advection step is x s ∗ j = x j + u j ( x , t ) ∆ t + (cid:112) D ( x , t )∆ W sj . (61)Similarly, according to (36) the value of φ ∗ at x s ∗ after the advection step is φ s ∗ ( x s ∗ , t + ∆ t ) = φ s ∗ ( x , t ) + S ( φ s ∗ ( x , t ) ; x , t ) ∆ t, (62)where S ( ψ ; x , t ) is a general source term. The specific version S ( ψ ; x , t ) = − ψ ∂u j ∂x j is used in section 5.For further analysis, the variable ∆ X j representing the difference between the location before and after advection isintroduced as follows: ∆ X j = x s ∗ j − x ( i ) j = u j ( x , t ) ∆ t + (cid:112) D ( x , t )∆ W sj , (63)Because of the Wiener increment, ∆ X j ∈ O (cid:0) ∆ t / (cid:1) . 23 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
In this work, the particle is not permitted to leave the bounds of the element, and so it follows that ∆ X j ∈ O (∆ x ) with ∆ x a representative grid spacing within an element (such as the average or minimum grid spacing). The timestep thus relates to ∆ x as ∆ t ∼ ∆ x and we can Taylor expand in x from φ s ∗ ( x s ∗ , t + ∆ t ) to φ s ∗ ( x , t + ∆ t ) , asfollows φ s ∗ ( x , t + ∆ t ) = φ s ∗ ( x s ∗ , t + ∆ t ) − ∆ X j ∂φ s ∗ ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t + (64)
12 ∆ X j ∆ X j ∂ φ s ∗ ∂x s ∗ j ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t + o (∆ t ) , where we use little o notation to denote terms which converge to 0 faster than the little o argument, i.e. lim y ↓ o ( y ) y = 0 ,for any y . The Taylor expansion gives φ s ∗ ( x , t + ∆ t ) prior to the remapping step in the semi-Lagagrangian algorithm.For a sufficiently fine grid, the spectral spatial interpolation error from the remapping step is smaller than ∆ x (andtherefore ∆ t ). Because the interpolation error of the remapping is o (∆ t ) , φ s ∗ ( x , t + ∆ t ) after remapping is thus alsorepresented by (64).Substituting (62) into (64) we find that φ s ∗ ( x , t + ∆ t ) = φ s ∗ ( x , t ) + S ( φ s ∗ ( x , t ) ; x , t ) ∆ t − (65) − ∆ X j ∂φ s ∗ ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t + 12 ∆ X i ∆ X j ∂ φ s ∗ ∂x s ∗ i ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t + o (∆ t ) . To derive an Eulerian form of this expression, we now seek to express ∆ X j ∂φ s ∗ ∂x s ∗ j (cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t and ∆ X i ∆ X j ∂ φ s ∗ ∂x s ∗ i ∂x s ∗ j (cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t in terms of x only. By the chain rule, we have that ∂φ s ∗ ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t = ∂φ s ∗ ∂x k (cid:12)(cid:12)(cid:12)(cid:12) x ,t +∆ t ∂x k ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t . (66)By (62), we have that ∂φ s ∗ ∂x k (cid:12)(cid:12)(cid:12)(cid:12) x ,t +∆ t = ∂φ s ∗ ∂x k (cid:12)(cid:12)(cid:12)(cid:12) x ,t + ∆ t ∂S ( φ s ∗ ( x , t ) ; x , t ) ∂x k . (67)By (61) the inverse derivative of ∂x k ∂x s ∗ j (cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t is ∂x s ∗ j ∂x k (cid:12)(cid:12)(cid:12)(cid:12) x ,t = δ jk + ∆ t ∂u j ∂x k (cid:12)(cid:12)(cid:12)(cid:12) x ,t + ∆ W sj ∂ (cid:16) √ D (cid:17) ∂x k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x ,t , (68)Applying the matrix inverse formula ( I + A ) − = I − A + A + · · · and the identity ∆ W sj ∆ W sk = ∆ tδ jk [31], wefind that ∂x k ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t = δ kj − ∆ t (cid:32) ∂u k ∂x j − ∂ √ D∂x k ∂ √ D∂x j (cid:33) − ∆ W sk ∂ √ D∂x j + o (∆ t ) , (69)where all the terms on the right hand side are evaluated at ( x , t ) ; we will use the convention for the rest of the derivationthat a term is evaluated at ( x , t ) if there is not specific indication otherwise. Substituting (69) and (67) into (66) itfollows that ∆ X j ∂φ s ∗ ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t = ∆ X j (cid:32) ∂φ s ∗ ∂x j − ∂φ s ∗ ∂x k ∆ W sk ∂ √ D∂x j (cid:33) + o (∆ t ) . (70)24 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT
Using the same analysis as in (66-69), we get that ∆ X i ∆ X j ∂ φ s ∗ ∂x s ∗ i ∂x s ∗ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t = ∆ X i ∆ X j ∂ φ s ∗ ∂x i ∂x j + o (∆ t ) . (71)Substituting (70,71) into (65), we can simplify using the identities ∆ X j = u j ( x , t ) ∆ t + (cid:112) D ( x , t )∆ W sj and ∆ W sj ∆ W sk = ∆ tδ jk and we get φ s ∗ ( x , t + ∆ t ) − φ s ∗ ( x , t ) = ∆ t (cid:18) S ( φ s ∗ ( x , t ) ; x , t ) − u j ∂φ s ∗ ∂x j + ∂∂x j (cid:18) D ∂φ s ∗ ∂x j (cid:19)(cid:19) −− ∂φ s ∗ ∂x j √ D ∆ W sj + o (∆ t ) . (72)Finally, dividing (72) by ∆ t and taking the limit as ∆ t ↓ , we get that ∂φ s ∗ ∂t + u j ∂φ s ∗ ∂x j + ∂φ s ∗ ∂x j √ D ˙ W sj − ∂∂x j (cid:18) D ∂φ s ∗ ∂x j (cid:19) = S ( φ s ∗ ; x , t ) , (73)where ˙ W sj is the weak time derivative of the multivariate Wiener process. This is a stochastic PDE (in Ito form) whichuniquely defines the PDF P φ ( ψ ; x , t ) of the random variable φ ∗ . Note that, just as ∆ W sj is the same Wiener incrementfor all initial points x in the s − th sample (i.e., all collocation points), then so also is ˙ W sj spatially independent. Theabove derivations shed light on why this spatial independence is necessary, since the derivative ∂x k ∂x s ∗ j (cid:12)(cid:12)(cid:12) x s ∗ ,t +∆ t from(69) would not be well-defined if the Wiener increments for different points x in the s − th sample were independent.From (73), we can derive the Fokker-Planck equation for P φ ( ψ ; x , t ) . To do this, we use the following Lemma from[25]: Lemma 1 (Sabel’nikov and Soulard)
If the stochastic field φ s ∗ ( x , t ) evolves by the stochastic PDE (in Stratanovichform) ∂φ s ∗ ∂t dt + dv j ◦ ∂φ s ∗ ∂x j = S ( φ s ∗ ; x , t ) , (74) with dv j being a stochastic advection term, and S ( ψ ; x , t ) being a deterministic and Lipschitz continuous function of ψ, x and t , then the Fokker-Planck equation for the PDF, P φ ( ψ ; x , t ) , of φ s ∗ is ∂ P φ ∂t + (cid:104) dv j (cid:105) dt + 12 (cid:68) dv j ∂dv k ∂x k (cid:69) dt ∂ P φ ∂x j = ∂∂x j (cid:18) (cid:104) dv j dv k (cid:105) dt ∂ P φ ∂x k (cid:19) − ∂∂ψ [ S ( ψ ; x , t ) P φ ] . (75)To use Lemma 1, we need to cast (73) in Stratanovich form. From Evans [31] we have the following Ito to Stratanovichconversion formula: B ( W , t ) ◦ d W = B ( W , t ) d W + 12 ∂B ij ∂x j ( W , t ) dt, (76)where W is the multivariate Wiener process, and B ( y , t ) is a C matrix function, with B ij being its components.Applying (76) to B ij ( W , t ) = (cid:112) D ( x , t ) δ i ∂φ s ∗ ∂x j for a fixed x , we have that25 high-order semi-Lagrangian method for the consistent Monte-Carlo solution of stochastic Lagrangiandrift-diffusion models coupled with Eulerian discontinuous spectral element method A P
REPRINT (cid:112) D ( x , t ) dW j ◦ ∂φ s ∗ ∂x j = (cid:112) D ( x , t ) ∂φ s ∗ ∂x j ◦ dW j = (cid:112) D ( x , t ) ∂φ s ∗ ∂x j dW j + 12 ∂∂W j (cid:18)(cid:112) D ( x , t ) ∂φ s ∗ ∂x j (cid:19) dt = (cid:112) D ( x , t ) ∂φ s ∗ ∂x j dW j + (cid:112) D ( x , t )2 ∂∂x j (cid:18) ∂φ s ∗ ∂W j (cid:19) dt = (cid:112) D ( x , t ) ∂φ s ∗ ∂x j dW j − (cid:112) D ( x , t )2 ∂∂x j (cid:18)(cid:112) D ( x , t ) ∂φ s ∗ ∂x j (cid:19) dt = √ D ∂φ s ∗ ∂x j dW j − ∂∂x j (cid:18) D ∂φ s ∗ ∂x j (cid:19) dt + √ D ∂ √ D∂x j ∂φ s ∗ ∂x j dt (77)where we used that (cid:112) D ( x , t ) is independent of W for the equality on the first line and to get from the second to thethird line. To get from the first to the second line, we applied (73) (76) is finally used to get from the third to the fourthline. Multiplying (73) by dt and substituting the last line of (77) into it, we find that ∂φ s ∗ ∂t dt + u j ∂φ s ∗ ∂x j dt + √ DdW sj ◦ ∂φ s ∗ ∂x j − √ D ∂ √ D∂x j ∂φ s ∗ ∂x j dt = S ( φ s ∗ ; x , t ) dt. (78)With dv j ≡ (cid:16) u j − √ D ∂ √ D∂x j (cid:17) dt + √ DdW sj , we have that (cid:104) dv j (cid:105) dt = (cid:16) u j − √ D ∂ √ D∂x j (cid:17) dtdt = u j − √ D ∂ √ D∂x j , (79)since dW sj has zero mean. Similarly, because dv j = √ DdW sj + o ( dt / ) and ∂dv k ∂x k = ∂ √ D∂x k dW sj + o ( dt / ) , and (cid:104) dW i dW j (cid:105) = δ ij dt , the second-order mean terms in (75) are respectively (cid:68) dv j ∂dv k ∂x k (cid:69) dt = √ D ∂ √ D∂x j , (80)and (cid:104) dv j dv k (cid:105) dt = Dδ jk . (81)Substituting (79-81) into (75) and simplifying, we get the Fokker-Planck equation for P ( ψ ; x , t ) ∂ P φ ∂t + u j ∂ P φ ∂x j = ∂∂x j (cid:18) D ∂ P φ ∂x j (cid:19) − ∂∂ψ [ S ( ψ ; x , t ) P φ ] ,,