# Stochastic Discontinuous Galerkin Methods (SDGM) Based on Fluctuation-Dissipation Balance

\lettrine[lines=3,lhang=0.3,nindent=0em] We introduce a general framework for approximating parabolic Stochastic Partial Differential Equations (SPDEs) based on fluctuation-dissipation balance. Using this approach we formulate Stochastic Discontinuous Galerkin Methods (SDGM). We show how methods with linear-time computational complexity can be developed for handling domains with general geometry and generating stochastic terms, handling both Dirichlet and Neumann boundary conditions. We demonstrate our approach on example systems and contrast with alternative approaches using direct stochastic discretizations based on random fluxes. We show how our Fluctuation-Dissipation Discretizations (FDD) framework allows for compensating for differences in dissipative properties of discrete numerical operators relative to their continuum counter-parts. This allows us to handle general heterogeneous discretizations, accurately capturing statistical relations. Our FDD framework provides a general approach for formulating SDGM discretizations and other numerical methods for robust approximation of stochastic differential equations.

Keywords: stochastic partial differential equations; fluctuation-dissipation balance; discontinuous Galerkin methods

## 1 Introduction

Stochastic Partial Differential Equations (SPDEs) [33, 54, 47, 62] arise in many settings including statistical inference [17, 48, 41], reduced descriptions of dynamical systems [58, 31, 37, 63, 49], and applications in the natural sciences and engineering [8, 10, 31, 29, 59, 43]. Applications include formulations of stochastic phase-field equations [26, 45, 14, 57], stochastic concentration equations [54, 7, 39], and fluctuating hydrodynamics [43, 29, 27, 6, 10, 46, 12] with fluid-structure interactions [10, 42, 8, 61]. The development of effective computational methods for SPDEs poses significant challenges often not encountered in the deterministic setting. The solutions of SPDEs take on the form of a measure on a space of functions or more often on distributions or generalized functions (separable Banach spaces) [44, 33, 54]. Stochastic equations are often driven by Gaussian noise in time and space with statistical structures ranging from -correlations (which are often not well-posed) [54, 62] to colored noise with prescribed spatial-temporal correlation functions [47, 8]. Such stochastic fields, when they do exist, are often non-differentiable in time, or do not have well-defined pointwise values, which requires interpretations in terms of a measure on an appropriate space of functions or distributions [47, 33]. These issues are further compounded in numerical methods when attempting to formulate approximations of the stochastic terms in the presence of truncation errors that arise from discretizations of the underlying differential operators [10, 7].

Many applications require preservation of inherent statistical structures to accurately capture stationary quantities or thermodynamic properties. We shall develop numerical methods to do this by driving them stochastically in a manner that compensates for discretization errors. We develop our approach utilizing considerations of fluctuation-dissipation balance. The fluctuation-dissipation principle is a well-established result in equilibrium statistical mechanics having its historic origins in the investigation of noise in electrical circuits (Johnson-Nyquist noise) [38]. For circuits, it was observed that the amplitude of fluctuations is closely related to electrical impedance [50]. These observations subsequently led to the formulation of general principles for physical systems relating the amplitude of fluctuations to the linear responses to perturbations [16, 52, 56]. We abstract these notions to the numerical setting by considering how fluctuations should relate to numerically captured dissipative relaxations of the system. This allows us to obtain discretization-dependent relations, which in turn can be used to develop strategies for the prescription of stochastic driving fields in numerical methods that preserve inherent statistical structures. These relations allow for taking into account the discretization errors and resulting particularities of the dissipative properties of the utilized discrete numerical operators. We have successfully used related ideas in our prior works when formulating stochastic numerical methods in [10, 7, 53, 9].

In this work, we abstract these ideas to formulate a general mathematical framework for stochastic numerical methods for SPDEs which we refer to as Fluctuation Dissipation Discretizations (FDD). For spatial-temporal numerical discretizations, we derive conditions relating the covariances of the stochastic driving terms to the corresponding covariances of the fluctuations in the statistical steady-state. We show how this FDD framework can be used to develop general strategies for prescribing spatial or temporal discretizations for stochastic equations that control discretization artifacts or attain discrete statistical structures amenable to efficient computational methods.

In practice, a central challenge is efficiently generating the prescribed stochastic driving fields. For finite differences with uniform and staggered mesh discretizations, the FDD framework shows covariances can be factored in terms of divergence and gradient mesh operators recovering methods for efficient stochastic driving terms related to random fluxes [10, 11, 60]. For non-uniform discretizations, as arise in finite volume and finite element methods, FDD gives much more complicated covariances and work has been done to address this by employing factorizations and multigrid techniques to obtain efficient Gibbs samplers to generate stochastic driving fields [7, 53]. Here, we take a different approach for non-uniform discretizations and general geometries avoiding the need for multigrid-based Gibbs samplers. We show how a judicious choice of spatial discretization can be used to achieve efficient stochastic methods. We develop a family of Stochastic Discontinuous Galerkin Methods (SDGM) based on FDD, which are applicable to general unstructured curvilinear grids and boundary conditions. A key feature of discretization a natural block-structure of the prescribed FDD covariances, which can be readily factored. This allows us to develop stochastic methods that have linear-time computational complexity under -refinement, which are also amenable also to parallelization.

We remark that our approaches based on fluctuation-dissipation balance are similar in spirit to other recent developments in numerical analysis following a long trend motivated by mimicking key inherent features of the mathematical problem being approximated [3, 5]. This includes mimetic methods (MM) [36], discrete exterior calculus (DEC) [28], and finite element exterior calculus (FEEC) [5, 32]. Each of these, in their own way, aim to preserve in the numerics inherent geometric features, relations in mechanics, or other properties important to the application domain. In the deterministic setting, this is often observed to be essential to obtaining efficient convergent methods, especially for finite element methods [3, 5, 32]. Here, our introduced FDD framework provides a way to preserve inherent statistical structures important when approximating stochastic systems.

We organize our paper by first discussing some background on stochastic partial differential equations (SPDEs) in Section 2. We then formulate our Fluctuation Dissipation Discretizations (FDD) framework and show how the FDD framework may be used to develop general strategies for prescribing spatial or temporal discretizations for approximating stochastic equations in Section 3. We then show how these ideas can be used to formulate Stochastic Discontinuous Galerkin Methods (SDGM) in Section 4. We develop computational methods having linear-time computational complexity for generating stochastic driving fields obeying our FDD framework in Section 4.1 and 4.2. We then make comparisons between our FDD-SDGM methods with the alternative approach of directly discretizing the stochastic terms using random fluxes in Section 5. These results highlight the utility of our FDD approach in providing a means to control discretizations artifacts and spurious long-range correlations. We then present results showing how our FDD-SDGM methods can handle general geometries and periodic, Neumann, and Dirichlet boundary conditions in Sections 5.1–5.3. The results show how our FDD approach can be used to devise general strategies for mitigating artifacts from the underlying spatial-temporal discretizations to develop robust numerical methods for stochastic differential equations.

## 2 Stochastic Partial Differential Equations

We consider parabolic stochastic partial differential equations (SPDEs) of the form

(1) |

where is a uniform elliptic operator, source term that is square-integrable, and a stochastic force [54, 47]. We take to be a Gaussian process that is -correlated in time with mean zero and prescribed covariance structure

(2) |

As a specific example, we shall consider the stochastic diffusion equation capturing linearized concentration fluctuations [7], which can be expressed in this form as

(3) |

where

(4) |

This corresponds to and spatial correlations . Here, denotes the diffusivity, and the reference concentration around which the system was linearized. This has -correlations in time and correlations in space [7]. The statistical steady-state solution has equilibrium fluctuations . As a consequence, the solutions of such stochastic differential equations are irregular, being non-differentiable in time and not having point-wise values but rather solutions in the sense of a measure on a space of distributions [33, 54, 47, 51]. This poses significant challenges both for analysis and numerical approximations. Many approaches for SPDEs rely on spectral approximations of the solutions [34, 62, 54]. We take a different approach, using instead finite differences and finite elements to discretize the differential operators, temporal dynamics, and stochastic driving fields [10, 7, 53, 9, 30, 1, 27].

For the linear stochastic differential equations 1, the solutions are Gaussian processes and have significant statistical structures that we can utilize in constructing our numerical approximations. The Gaussianity ensures that the process is determined if we know the mean values and spatial-temporal correlation functions of the process [33, 47]. The solution can be expressed formally as

(5) |

Here, denotes the semigroup of solution operators for the parabolic PDE with operator and is a formal space-time Brownian motion with integral given the Ito interpretation [33, 47, 51]. The mean of this process is given formally by

(6) |

When , the stationary distribution has . The spatial-temporal correlation are then

(7) |

where , and is the steady-state covariance of fluctuations of . This can be generalized to non-zero by subtracting off the mean behavior to consider with . Since the process is Gaussian we can formally express the probability measure in terms of a formal density as

(8) |

where is the normalization factor. This should be interpreted as defining the statistics of finite dimensional marginals [47].

### 2.1 Spatial Discretization

We approximate solutions of the stochastic partial differential equation 1 using a discrete stochastic dynamical system of the form

(9) |

We denote by the state of the system with , by the system forcing with , and by the stochastic driving force with . We take to be a Gaussian process with mean zero and a yet-to-be-determined covariance

(10) |

where . The solution is a Gaussian process completely determined by its mean and spatial-temporal covariances. Similar to equations 5–8, we have that the system has the mean

(11) |

When we have and the temporal correlations are

(12) |

where . For the covariance , we assume throughout that we have ergodicity so that as with the steady-state covariance . This can be generalized when is non-zero by subtracting off the mean with .

### 2.2 Temporal Discretization

We shall also consider the role of temporal discretization errors. We consider for illustration the case of an Euler discretization

(13) |

where is the Brownian increment . This is a discrete-time process (multi-variate Gaussian for finite number of steps) which is determined by its mean and covariance. The mean is

(14) |

In the case when we have . The covariance is

(15) |

where , and is the steady-state covariance. This can be generalized for non-zero by subtracting the mean with .

For these spatial-temporal discretization cases, we show how we can introduce stochastic driving terms or that compensate in numerical methods for spatial truncation errors or those arising from the choice of temporal approximations and numerical integrator.

## 3 Fluctuation-Dissipation Discretizations (FDD) Framework

We take the general approach for discretizations of the stochastic terms in equation 1 of using variations of fluctuation-dissipation balance adapted to the setting of numerical discretizations in space or in time. We express equation 9 representing a semi-discretization with continuous time as

(16) |

We give this the interpretation of an Ito Stochastic Process [51]. Letting , we have . We have from Ito’s Lemma [51] and equation 16 that

(17) |

where . This gives in the statistical steady-state with as the relationship

(18) |

This establishes a fluctuation-dissipation relationship for the discrete system given in equations 9 and 16. Throughout, we assume that is negative-definite and that the stochastic driving fields yield ergodic behaviors for the stochastic process.

In the case that , we have

(19) |

If we make a specific choice for for the stochastic driving terms for the system, this establishes what equilibrium fluctuations will be obtained.

We can also obtain analogous relations when time has been discretized as in equation 13. In this case, we have for that

(20) |

This was obtained by substituting for using equation 13 with . In the statistical steady-state, we have as . This gives the fluctuation-dissipation relation taking temporal discretization into account

(21) |

We see the temporal discretization results in a higher-order correction for the stochastic driving fields involving the time-step size to compensate for the additional truncation errors of the numerical methods.

In the case that , we have the further useful relation

(22) |

For a specific choice for for the stochastic driving terms for the system, this establishes what equilibrium fluctuations will be obtained. We see how the temporal discretization augments the obtained equilibrium fluctuations relative to equation 19.

The fluctuation-dissipation relations in equation 18 and 21 provide a few strategies for discretizing stochastic equations. For many stochastic equations arising from applications in physics the equilibrium fluctuations are known from the energy of the system . From equation 8, when is quadratic in , or when the system is linearized, this yields the equilibrium covariance . When discretizing the energy this yields a covariance for the discretized system. One strategy is to discretize the stochastic driving terms by using relations in equation 18 or 21 to obtain a covariance for the stochastic driving terms in the discrete system. This ensures the stochastic driving terms interact with the discrete numerical operator L to yield equilibrium fluctuations with covariance C, even in the presence of discretization errors.

An alternative strategy is to use equation 19 or 22 to formulate that yields an equilibrium covariance with acceptable properties such as localization in space with rapid decay to approximate -correlations or reduce artifacts at coarse-refined interfaces within structured discretizations [7]. Without these considerations, arbitrarily approximating directly can result in ringing phenomena or other artifacts in the numerics [7]. We also see from equations 12 and 15 the relations can be used to help ensure accurate autocorrelation functions for the discretized stochastic system.

We have used related ideas in our prior works [10, 8, 58] to derive stochastic discretizations for finite difference methods on uniform meshes, staggered meshes, and spatially adaptive meshes [10, 60, 7] and for finite element methods [53]. In the subsequent sections, we utilize the FDD framework and related ideas to develop Stochastic Discontinuous Galerkin Methods (SDGM) to handle domains with general geometries and Neumann and Dirichlet boundary conditions.

## 4 Stochastic Discontinuous Galerkin Methods (SDGM)

We develop Stochastic Discontinuous Galerkin Methods (SDGM) by introducing stochastic driving terms for approximating solutions of the SPDE in equation 1. The Discontinuous Galerkin (DG) method was first introduced in 1973 by Reed and Hill for solving the neutron transport equation [55]. Subsequently, Cockburn and Shu extended the DG method to general systems of nonlinear hyperbolic conservation laws [22, 21, 20, 19, 24]. Several extensions have been developed for elliptic and parabolic problems [2, 13, 15, 23]. These extensions subsequently have been presented in the context of a unified framework in [4].

### 4.1 Discontinuous Galerkin (DG) Methods

The DG method is a finite element method, suitable for use on unstructured meshes, that makes use of a discontinuous, piecewise polynomial function space. To begin, we describe the DG discretization of the the deterministic heat equation on a spatial domain ,

(23) | |||||

(24) | |||||

(25) |

where Dirichlet and Neumann boundary conditions are imposed on the boundary . We discretize the spatial domain by introducing a mesh

(26) |

One of the main advantages of the DG method is the great amount of flexibility it affords in the construction of the mesh. The mesh can be be entirely unstructured, non-conforming, and can consist of elements of different shapes. In this work, we restrict the elements of the mesh to be straight-sided quadrilaterals. The generalization to curved elements is possible using a standard isoparametric mapping [35]. We now fix a polynomial degree . On each element , we define the space of bivariate polynomials of degree at most in each variable,

(27) |

where the sum is taken over all multi-indices such that , and the notation is used to mean . We now introduce the discontinuous finite element space defined by

(28) |

obtained by requiring that each function , when restricted to an element , lies in the corresponding local polynomial space . We note that no continuity is enforced between the mesh elements. We consider two neighboring elements, and , that share a face . A function is not well-defined on , and thus we define to be the trace of on from within , and, analogously, to be the trace of on from within . Similarly, refers to a vector normal to , facing outward from , and is the normal facing outward from . At this point, it will be useful to introduce the average and jump of a scalar function , which we define by

(29) |

Similarly, for a vector-valued function , we define

(30) |

We remark that the jump of a scalar is a vector parallel to the normal direction, and the jump of a vector is a scalar.

To obtain the DG discretization of equation 23–25, we transform the equation into a system of first order equations by introducing the gradient , thus obtaining the equivalent system of equations

(31) | ||||

(32) |

We look for an approximate solution , . To obtain the variational form, we multiply equation 32 by a test function , and equation 31 by a test function . We integrate the resulting equations over the spatial domain . We then integrate the divergence and gradient terms by parts element-by-element, giving rise to the following weak form,

(33) | |||

(34) |

where denotes all interior and exterior edges in the triangulation , and and are yet-to-be-defined numerical flux functions. In this work, we consider the local discontinuous Galerkin (LDG) method [23], which proceeds by choosing the fluxes on interior edges by

(35) | ||||

(36) |

for parameters and . On the Dirichlet boundary, we choose

(37) | ||||

(38) |

and on the Neumann boundary,

(39) | ||||

(40) |

The parameter is necessary to stabilize the method, and can be thought of as introducing an artificial viscosity [25]. Of particular interest to this work is the so-called minimal dissipation LDG method [18], which allows for to be taken as identically zero on all interior edges by a careful choice of the parameter . Because the flux is independent of , it is possible to solve for in an element-by-element fashion.

We choose a nodal basis for the space , and write to represent the vector of degrees of freedom defining . Likewise, represents the vector of degrees of freedom defining . We then rewrite equations 33 and 34 as

(41) | ||||

(42) |

or, equivalently,

(43) |

where is the is standard mass matrix, and and are the divergence and gradient operators defined by the respective bilinear forms. The mass matrix times a vector-valued quantity is understood to be applied component-wise. We note that a simple integration by parts shows that . The matrix corresponds to the stabilization terms caused by . This allows us to define the Laplacian operator

(44) |

which is a symmetric, negative-definite linear operator. In the particular cases of Neumann and periodic boundary conditions, we choose to be identically zero, and thus , so we can write .

### 4.2 Stochastic Driving Fields for DG

We now consider discretization of the stochastic diffusion equation 3 where throughout we take diffusivity and reference concentration . This corresponds to a stochastic version of equations 23–25. In equation 23, we account for the fluctuations by taking the forcing term to be a Gaussian stochastic field that is -correlated in time with mean zero and with spatial covariance . The DG discretization gives a linear stochastic dynamical system of the form 9, where the linear operator is given by

(45) |

The covariance structure of the fluctuations is related to the linear operator and equilibrium covariance of the system by the FDD framework through equation 18 by

(46) |

At equilibrium, we have from equation 3 that . In the discrete setting, the associated energy motivates taking the discretized covariance to be for the DG solution . From equation 18 and 45, we use stochastic driving terms with covariance

(47) |

To obtain this expression, we have used that and are symmetric linear operators. By equation 19, this prescribes for a given choice of mesh a covariance for the stochastic driving fields that takes into account the dissipative properties of the DG Laplacian operator to ensure at statistical steady-state the covariance is achieved.

In order for this to be useful in practice, we must develop computational methods to efficiently generate the random forcing terms with the prescribed covariance given by equation 47. In principle, taking a Cholesky factorization would provide such random variates through with a standard Gaussian random variable. However, this technique does not fully utilize the sparse block structure of the DG discretizations, and thus is often prohibitively expensive. Furthermore, such direct solvers can prove problematic to parallelize.

In the following sections, we shall develop more efficient methods for generating the stochastic variates, while also appropriately treating the boundary conditions, by utilizing the block structure of the matrices arising in DG. Our methods scale as , where the total number of degrees of freedom is given by , is the number of mesh elements, the polynomial degree, and the spatial dimension. We focus primarily on -refinement where is held fixed and develop methods with linear computational complexity as . Our methods have the additional advantage of being computed element-wise, avoiding communication between mesh elements, facilitating straightforward parallelization.

#### 4.2.1 Neumann and periodic boundary conditions with fluctuations

We now consider the case of periodic or pure Neumann boundary conditions. In this case, we modify the finite element space to consist only of mean-zero functions,

(48) |

in order to ensure that the Laplacian operator is nonsingular. Other than this modification, the method is as described in the preceding section. Using the minimal dissipation LDG method, we can, in this case, set the LDG stabilization parameter to be identically zero on all edges. As a result, the stabilization matrix is also zero, and thus the covariance of is given by

(49) |

We recall that the DG mass matrix is given by

(50) |

where the functions form a basis for the finite element space . This matrix is symmetric positive-definite, and due to the discontinuous nature of the basis functions , has a natural block-diagonal structure, with blocks corresponding to each element in the triangulation . Thus, each block of the mass matrix is also symmetric positive-definite, allowing for efficient block-wise computation of the Cholesky factorization of the inverse mass matrix.

These operations can be performed element-by-element and in parallel. Given the block structure, the Cholesky factorization requires operations. The number of elements in the mesh is , the polynomial order used, and the spatial dimension. When performing -refinement (fixing degree ), we have linear computational complexity with scaling where . For -refinement we see the complexity would be dimension dependent.

We can use this approach by defining the matrix by

(51) |

so that we obtain

(52) |

We can generate variates using where is a Gaussian with -correlation in time and spatial components with mean zero and variance one. This has the covariance

(53) |

thus giving a linear-time complexity method under -refinement for generating the random variates prescribed by the FDD framework in equation 18 and 47.

#### 4.2.2 Dirichlet boundary conditions with fluctuations

We now consider the case of Dirichlet boundary conditions, which are enforced by penalizing the difference between the approximate solution and the prescribed boundary value by means of a penalty parameter , as in equation 38. In order for the problem to be well-posed, we must choose on all edges in . As before, we set on all interior edges. This results in a nonzero penalty matrix , whose entries are given by

(54) |

This matrix is symmetric and negative-semidefinite, with an element-wise block diagonal structure. The nonzero blocks of correspond to elements of the mesh with at least one edge on the domain boundary.

We remark that the discontinuous Galerkin method enforces the Dirichlet boundary conditions weakly through the penalization procedure described above. This formulation can permit fluctuations at the domain boundary. In this case, we consider the target covariance to be given by . The FDD framework based on fluctuation-dissipation balance in equation 18 gives the prescribed covariance in equation 47. To efficiently sample a random variable with this covariance, we require the eigendecomposition of the penalty matrix,

(55) |

This factorization can be computed efficiently in linear-time complexity under -refinement by taking advantage of the block structure of the matrix . We then let and be independent, identically distributed Gaussian random variables with mean zero and variance one. We define the matrix by

(56) |

where as in the Neumann case, is the Cholesky factor of the inverse mass matrix. We also define by

(57) |

Then, by setting , we obtain

(58) |

This formulation leads to fluctuations on the domain boundary, which may be considered undesirable. If we wish not to introduce fluctuations at the boundary, we may modify the target covariance as follows. We recall that our basis functions are defined to be nodal interpolants. We refer to an index as a boundary index if the nodal interpolation point corresponding to basis function lies on the domain boundary. We then define the modified target covariance matrix by

(59) |

In this case, the fluctuation covariance is given by

(60) |

Here, we note that is no longer symmetric. In order to ensure symmetry, we modify the Laplacian operator as follows. Rather than weakly enforce the Dirichlet conditions using a penalty term, we strongly enforce these conditions by fixing the degrees of freedom at the boundary to be equal to their specified Dirichlet values, . The resulting linear operator has the form

(61) |

where the matrix is a diagonal matrix defined by

(62) |

We note that , and thus the boundary penalty matrix is not needed in this formulation. Additionally, we have that and . Thus, we have

(63) |

Thus, in this case, we can define the matrix by

(64) |

where as before is the Cholesky factor of the inverse mass matrix. We define the fluctuation forcing term by , and thus obtain

(65) |

Given the block structure each required factorization operation only costs linear-time computational complexity under -refinement. In this way, we can efficiently generate in practice the needed stochastic driving terms with the prescribed covariance structure given by the FDD framework in equation 18 and 47.

#### 4.2.3 Temporal discretization

We use an Euler-Maruyama time discretization [40], which results in the fully-discrete system given by

(66) |

with time step . The term is used to denote a Gaussian random variable with mean zero and covariance

(67) |

This random variable is generated at each time step by sampling the random fluctuation forcing term according to the desired boundary conditions, using the methodology described in Sections 4.2.1 and 4.2.2.

We mention that while we could in principle use the relation in equation 21 to make further corrections to the covariance of the stochastic driving terms . In practice, the time steps we take in our numerical calculations are sufficiently small that these higher-order corrections do not play a significant role. Therefore, for the sake of simplicity, we do not incorporate these terms into our discretization in this work. We also remark that other stochastic time-step integrators could also be developed with relations derived like in equation 21 to make further corrections to the covariance of the stochastic driving terms.

## 5 Results

We investigate how these methods perform in practice by considering problems involving different geometries and boundary conditions. We show how the stochastic discretizations perform when accounting for the Neumann and Dirichlet boundary conditions. We also consider how the stochastic methods perform on different types of meshes ranging from structured to unstructured, and when transforming from straight-sided meshes to curved geometries. We investigate how the methods behave using both low-order and high-order polynomial spaces. These investigations illustrate some of the key features of our stochastic DG methods.

We compare the results of our fluctuation-dissipation based discretizations for our stochastic DG methods described in Section 4.2 with those that would be obtained from the intuitive approach of directly using random fluxes. The random fluxes are given by , where is the vector whose entries are given by independent, identically-distributed Gaussian random variables, and represents a random function when expanded in terms of the chosen basis of the function space . To characterize how stochastic features of the system are captured by the stochastic numerical methods we investigate the covariance structure of fluctuations using Monte-Carlo sampling with estimator . We take only samples with for typically covering about of our samples to exclude initial transients in the stochastic dynamics. When considering the stochastic systems we assume throughout ergodicity.

We remark that in our methods if, instead of using a nodal basis, we were to use a basis consisting of orthogonal polynomials on each element, then the resulting mass matrix would be the identity matrix, and the random fluxes described above would agree with those obtained using the fluctuation-dissipation principle. This would hold except in the case of weakly-enforced Dirichlet conditions, in which case the boundary penalization term must also be taken into account. In order to ensure that the units are consistent between our comparisons, we multiply the random fluxes by a geometric factor of , which represents the spatial resolution of the DG method.

We compare the results of our fluctuation-dissipation based SDGM and what would be obtained for random fluxes in the cases of periodic, Neumann, and Dirichlet boundary conditions. These present a strong test of how methods capture inherent statistical structure in these stochastic systems.

### 5.1 Periodic Boundary Conditions

We consider the case of periodic boundary conditions for the domain having a simple geometry given by . We enforce periodic boundary conditions on the edges using symmetry of the nodes. Relative to many other discretization methods, an advantage of the DG method is the ability to use general, unstructured meshes. To illustrate this property of the method, we use the mesh shown in Figure 1. We choose the finite element function space to consist of piecewise bilinear polynomials in each element, thus resulting in a method whose spatial convergence rate is . The function space is restricted to the space of mean-zero functions on to ensure non-singularity of the Laplacian operator.

We choose the time step to be given by , and approximate the covariance of the solution using Monte Carlo samples. We let denote the covariance of the solution obtained using forcing terms determined according to the fluctuation-dissipation principle, and let denote the covariance of the solution obtained using a random flux forcing function. We compare these two empirical covariance matrices with the target covariance, which is given by the inverse mass matrix, . The relative error for the fluctuation-dissipation case was 3.5%. We compare the resulting matrices in Figure 2. We notice that the matrix obtained using the fluctuation-dissipation principle displays excellent agreement with the target covariance matrix, with close to no correlations in between elements. On the other hand, the covariance matrix obtained using random fluxes does not accurately reproduce the target covariance, and displays significant spurious long-range correlations. For ease of comparison, we note that , and thus we expect and to well-approximate the identity matrix as well. We show both these matrices in Figure 4. We note that , resulting from the use of random fluxes, does not display good agreement with the identity matrix, illustrating the spurious correlations that arise from this strategy.

### 5.2 Neumann Boundary Conditions

Now, we extend the numerical tests beyond the case of periodic boundary conditions. We consider a curved geometry that is given by a warped annulus. We use a high-order polynomial basis. The mesh used for this problem is shown in Figure 5. This examples illustrates the use of Neumann boundary conditions, curved geometries via an isoparametric mapping, and the use of high-order polynomial spaces. As in the case of periodic boundary conditions, we restrict the finite element space to mean-zero functions, so that the Laplace operator is nonsingular. We note again that the minimal dissipation LDG method allows for the stabilization parameter to be identically zero on all edges in the mesh. Homogeneous Neumann conditions are enforced at the domain boundary.

We choose the time step to be and perform the same comparisons as in the previous case. We estimate the covariance matrix using samples. The relative error for the fluctuation-dissipation case was 4.8%. Comparisons of the matrices are shown in Figure 7. Additionally, in Figure 6, display a single row of the the matrices and , which represents the correlation between a specified nodal degree of freedom with all other nodal degrees of freedom. This figure illustrates the clean -correlation that arises from using the fluctuation-dissipation forcing terms. On the other hand, we see that using random flux forcing terms does not result in the desired -correlation.

### 5.3 Dirichlet Boundary Conditions

For a final test case, we consider , with homogeneous Dirichlet conditions imposed on . We discretize the geometry using a regular Cartesian grid, and use biquadratic polynomials. First we consider allowing fluctuations on the domain boundary, which we achieve by weakly enforcing the Dirichlet conditions using the penalty parameter, which we choose to be strictly positive on all edges lying on , and equal to zero on all interior edges. The target covariance for this case is, as in the preceding cases, . This case requires treatment of the stabilization term , as described in Section 4.2.2. Using Monte Carlo samples, we find the relative error between and in the matrix norm is 3.5%. A comparison of the empirical covariance matrices for this case is shown in Figure 8

We additionally consider the case where no fluctuations are permitted on the domain boundary. The target covariance in this case is given by , which is obtained from by setting whenever either or is the index of a boundary node. To ensure symmetry of the resulting covariance matrix , we enforce the homogeneous Dirichlet conditions strongly rather than through a penalty term. We compare the resulting covariance matrices in Figure 9.

## 6 Conclusions

We have introduced a general approach referred to as Fluctuation Dissipation Discretizations (FDD) for approximating SPDEs using fluctuation-dissipation balance. Our FDD framework provides for a given numerical method’s spatial and temporal discretizations a prescribed way to discretize the stochastic driving terms in order to take into account the dissipative properties of the discrete numerical operators and how this affects propagation of fluctuations in the system. Using these ideas, we have developed general Stochastic Discontinuous Galerkin Methods (SDGM) for approximating solutions of parabolic stochastic partial differential equations. We developed practical methods with linear-time computational complexity under -refinement for generating the needed random variates with the prescribed covariance structure for the DG discretizations on general geometries. We have also introduced methods for handling periodic, Neumann, or Dirichlet boundary conditions. We expect many of the ideas behind our FDD framework to provide helpful guidelines in the further development of SDGMs and other numerical discretizations for approximating stochastic partial differential equations.

## 7 Acknowledgments

We acknowledge support to P.J.A. and N.T. from research grant DOE ASCR CM4 DE-SC0009254. N.T. was also supported by the NSF-MSPRF program, and P.J.A. was also supported by NSF Grant DMS - 1616353. W.P. was supported by the Department of Defense through the National Defense Science and Engineering Graduate Fellowship Program and by the Natural Sciences and Engineering Research Council of Canada.

## References

- [1] E. J. Allen, S. J. Novosel, and Z. Zhang. Finite element and difference approximation of some linear stochastic partial differential equations. Stochastics and Stochastic Reports, 64(1-2):117–142, may 1998.
- [2] Douglas N. Arnold. An interior penalty finite element method with discontinuous elements. SIAM Journal on Numerical Analysis, 19(4):742–760, aug 1982.
- [3] Douglas N Arnold, Pavel B Bochev, Richard B Lehoucq, Roy A Nicolaides, and Mikhail Shashkov. Compatible Spatial Discretizations, volume 142. Springer Science & Business Media, 2007.
- [4] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, and L. Donatella Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39(5):1749–1779, 2002.
- [5] Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. Finite element exterior calculus: from hodge theory to numerical stability. Bulletin of the American Mathematical Society, 47(2):281–354, jan 2010.
- [6] Paul J. Atzberger. Velocity correlations of a thermally fluctuating brownian particle: A novel model of the hydrodynamic coupling. Physics Letters A, 351(4-5):225–230, mar 2006.
- [7] Paul J. Atzberger. Spatially adaptive stochastic numerical methods for intrinsic fluctuations in reaction-diffusion systems. Journal of Computational Physics, 229(9):3474–3501, may 2010.
- [8] Paul J. Atzberger. Stochastic eulerian lagrangian methods for fluid-structure interactions with thermal fluctuations. Journal of Computational Physics, 230(8):2821–2837, apr 2011.
- [9] Paul J. Atzberger. Incorporating shear into stochastic eulerian–lagrangian methods for rheological studies of complex fluids and soft materials. Physica D: Nonlinear Phenomena, 265:57–70, dec 2013.
- [10] Paul J. Atzberger, Peter R. Kramer, and Charles S. Peskin. A stochastic immersed boundary method for fluid-structure dynamics at microscopic length scales. Journal of Computational Physics, 224(2):1255–1292, jun 2007.
- [11] Florencio Balboa, John B Bell, Rafael Delgado-Buscalioni, Aleksandar Donev, Thomas G Fai, Boyce E Griffith, and Charles S Peskin. Staggered schemes for fluctuating hydrodynamics. Multiscale Modeling & Simulation, 10(4):1369–1408, 2012.
- [12] Florencio Balboa Usabiaga, John B. Bell, Rafael Delgado-Buscalioni, Aleksandar Donev, Thomas G. Fai, Boyce E. Griffith, and Charles S. Peskin. Staggered schemes for fluctuating hydrodynamics. Multiscale Modeling & Simulation, 10(4):1369–1408, jan 2012.
- [13] F. Bassi and S. Rebay. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier-Stokes equations. Journal of Computational Physics, 131(2):267–279, mar 1997.
- [14] Lorenzo Bertini, Stella Brassesco, Paolo Buttà, and Errico Presutti. Stochastic phase field equations: existence and uniqueness. Annales Henri Poincaré, 3(1):87–98, mar 2002.
- [15] F. Brezzi, G. Manzini, D. Marini, P. Pietra, and A. Russo. Discontinuous Galerkin approximations for elliptic problems. Numerical Methods for Partial Differential Equations, 16(4):365–378, 2000.
- [16] Herbert B. Callen and Theodore A. Welton. Irreversibility and generalized noise. Physical Review, 83(1):34–40, jul 1951.
- [17] Igor Cialenco. Statistical inference for SPDEs: an overview. Statistical Inference for Stochastic Processes, feb 2018.
- [18] Bernardo Cockburn and Bo Dong. An analysis of the minimal dissipation local discontinuous Galerkin method for convection-diffusion problems. Journal of Scientific Computing, 32(2):233–262, mar 2007.
- [19] Bernardo Cockburn, Suchung Hou, and Chi-Wang Shu. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case. Mathematics of Computation, 54(190):545–545, 1990.
- [20] Bernardo Cockburn, San-Yih Lin, and Chi-Wang Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems. Journal of Computational Physics, 84(1):90–113, 1989.
- [21] Bernardo Cockburn and Chi-Wang Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Mathematics of Computation, 52(186):411–411, 1989.
- [22] Bernardo Cockburn and Chi-Wang Shu. The Runge-Kutta local projection -discontinuous-Galerkin finite element method for scalar conservation laws. ESAIM: Mathematical Modelling and Numerical Analysis, 25(3):337–361, 1991.
- [23] Bernardo Cockburn and Chi-Wang Shu. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35(6):2440–2463, dec 1998.
- [24] Bernardo Cockburn and Chi-Wang Shu. The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. Journal of Computational Physics, 141(2):199–224, apr 1998.
- [25] Bernardo Cockburn and Chi-Wang Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing, 16(3):173–261, Sep 2001.
- [26] Giuseppe Da Prato and Arnaud Debussche. Stochastic cahn-hilliard equation. Nonlinear Analysis: Theory, Methods & Applications, 26(2):241–263, jan 1996.
- [27] Steven Delong, Yifei Sun, Boyce E. Griffith, Eric Vanden-Eijnden, and Aleksandar Donev. Multiscale temporal integrators for fluctuating hydrodynamics. Physical Review E, 90(6), dec 2014.
- [28] M. Desbrun, A.N. Hirani, and J.E. Marsden. Discrete exterior calculus for variational problems in computer vision and graphics. In 42nd IEEE International Conference on Decision and Control (IEEE Cat. No.03CH37475). IEEE.
- [29] Aleksandar Donev, Eric Vanden-Eijnden, Alejandro Garcia, and John Bell. On the accuracy of finite-volume schemes for fluctuating hydrodynamics. Communications in Applied Mathematics and Computational Science, 5(2):149–197, jun 2010.
- [30] Qiang Du and Tianyu Zhang. Numerical approximation of some linear stochastic partial differential equations driven by special additive noises. SIAM Journal on Numerical Analysis, 40(4):1421–1445, jan 2002.
- [31] Ronald Forrest Fox and George E. Uhlenbeck. Contributions to nonequilibrium thermodynamics. II. fluctuation theory for the Boltzmann equation. Physics of Fluids, 13(12):2881, 1970.
- [32] Andrew Gillette, Michael Holst, and Yunrong Zhu. Finite element exterior calculus for evolution problems. Journal of Computational Mathematics, 35(2):187–212, mar 2017.
- [33] Martin Hairer. An introduction to stochastic PDEs. arXiv, 2009.
- [34] Martin Hairer and Andrew M. Stuart. Weak error estimates for trajectories of spdes for spectral galerkin discretization. arXiv, 2016.
- [35] Jan S. Hesthaven and Tim Warburton. Nodal Discontinuous Galerkin Methods. Springer New York, 2008.
- [36] J. Hyman, J. Morel, M. Shashkov, and S. Steinberg. Mimetic finite difference methods for diffusion equations. Computational Geosciences, 6(3/4):333–352, 2002.
- [37] Panos Stinis Jing Li. Mori-zwanzig reduced models for uncertainty quantification. arXiv, 2018.
- [38] J. B. Johnson. Thermal agitation of electricity in conductors. Physical Review, 32(1):97–109, jul 1928.
- [39] Changho Kim, Andy Nonaka, John B. Bell, Alejandro L. Garcia, and Aleksandar Donev. Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach. The Journal of Chemical Physics, 146(12):124110, mar 2017.
- [40] Peter E. Kloeden and Eckhard Platen. Numerical Solution of Stochastic Differential Equations. Springer Berlin Heidelberg, 1992.
- [41] Timo Koski and Wilfried Loges. Asymptotic statistical inference for a stochastic heat flow problem. Statistics & Probability Letters, 3(4):185–189, jul 1985.
- [42] Peter R. Kramer, Charles S. Peskin, and Paul J. Atzberger. On the foundations of the stochastic immersed boundary method. Computer Methods in Applied Mechanics and Engineering, 197(25-28):2232–2249, apr 2008.
- [43] L. D. Landau and E. M. Lifshitz. Fluid Mechanics. Butterworth Heineman, Oxford, UK, 3 edition, 1986.
- [44] E.H. Lieb and M. Loss. Analysis. American Mathematical Society, 2001.
- [45] Ernesto A. B. F. Lima, Regina C. Almeida, and J. Tinsley Oden. Analysis and numerical solution of stochastic phase-field models of tumor growth. Numerical Methods for Partial Differential Equations, 31(2):552–574, oct 2014.
- [46] Guang Lin, Xiaoliang Wan, Chau-Hsing Su, and George Em Karniadakis. Stochastic computational fluid mechanics. Computing in Science & Engineering, 9(2):21–29, March 2007.
- [47] Wei Liu and Michael Röckner. Stochastic Partial Differential Equations: An Introduction. Springer International Publishing, 2015.
- [48] S. V. Lototsky. Statistical inference for stochastic parabolic equations: a spectral approach. Publicacions Matemàtiques, 53:3–45, jan 2009.
- [49] Hazime Mori. Transport, collective motion, and Brownian motion. Progress of Theoretical Physics, 33(3):423–455, mar 1965.
- [50] H. Nyquist. Thermal agitation of electric charge in conductors. Physical Review, 32(1):110–113, jul 1928.
- [51] Bernt Øksendal. Stochastic Differential Equations: An Introduction with Applications. Springer Berlin Heidelberg, 2003.
- [52] Lars Onsager. Reciprocal relations in irreversible processes. i. Physical Review, 37(4):405–426, feb 1931.
- [53] Pat Plunkett, Jonathan Hu, Christopher Siefert, and Paul J. Atzberger. Spatially adaptive stochastic methods for fluid–structure interactions subject to thermal fluctuations in domains with complex geometries. Journal of Computational Physics, 277:121–137, nov 2014.
- [54] Giuseppe Da Prato and Jerzy Zabczyk. Stochastic Equations in Infinite Dimensions. Cambridge University Press, 2009.
- [55] W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. Los Alamos Report LA-UR-73-479, 1973.
- [56] Linda E. Reichl. A Modern Course in Statistical Physics. Wiley-VCH Verlag GmbH & Co. KGaA, apr 2016.
- [57] J.M. Sancho, J. García-Ojalvo, and H. Guo. Non-equilibrium ginzburg-landau model driven by colored noise. Physica D: Nonlinear Phenomena, 113(2-4):331–337, mar 1998.
- [58] Gil Tabak and Paul J. Atzberger. Stochastic reductions for inertial fluid-structure interactions subject to thermal fluctuations. SIAM Journal on Applied Mathematics, 75(4):1884–1914, jan 2015.
- [59] B. Uma, T. N. Swaminathan, R. Radhakrishnan, D. M. Eckmann, and P. S. Ayyaswamy. Nanoparticle brownian motion and hydrodynamic interactions in the presence of flow fields. Physics of Fluids, 23(7), jul 2011.
- [60] Y. Wang, H. Lei, and P. J. Atzberger. Fluctuating hydrodynamic methods for fluid-structure interactions in confined channel geometries. Applied Mathematics and Mechanics, 39(1):125–152, dec 2017.
- [61] Y. Wang, J. K. Sigurdsson, and P. J. Atzberger. Fluctuating hydrodynamics methods for dynamic coarse-grained implicit-solvent simulations in LAMMPS. SIAM Journal on Scientific Computing, 38(5):S62–S77, jan 2016.
- [62] Zhongqiang Zhang and George Em Karniadakis. Numerical Methods for Stochastic Partial Differential Equations with White Noise. Springer, 2017.
- [63] Robert Zwanzig. Memory effects in irreversible thermodynamics. Physical Review, 124(4):983–992, nov 1961.