On the Accuracy of Finite-Volume Schemes for Fluctuating Hydrodynamics

# On the Accuracy of Finite-Volume Schemes for Fluctuating Hydrodynamics

Aleksandar Donev Lawrence Livermore National Laboratory, P.O.Box 808, Livermore, CA 94551-9900 Center for Computational Science and Engineering, Lawrence Berkeley National Laboratory, Berkeley, CA, 94720    Eric Vanden-Eijnden Courant Institute of Mathematical Sciences, New York University, New York, NY 10012    Alejandro L. Garcia Department of Physics, San Jose State University, San Jose, California, 95192    John B. Bell Center for Computational Science and Engineering, Lawrence Berkeley National Laboratory, Berkeley, CA, 94720
###### Abstract

This paper describes the development and analysis of finite-volume methods for the Landau-Lifshitz Navier-Stokes (LLNS) equations and related stochastic partial differential equations in fluid dynamics. The LLNS equations incorporate thermal fluctuations into macroscopic hydrodynamics by the addition of white-noise fluxes whose magnitudes are set by a fluctuation-dissipation relation. Originally derived for equilibrium fluctuations, the LLNS equations have also been shown to be accurate for non-equilibrium systems. Previous studies of numerical methods for the LLNS equations focused primarily on measuring variances and correlations computed at equilibrium and for selected non-equilibrium flows. In this paper, we introduce a more systematic approach based on studying discrete equilibrium structure factors for a broad class of explicit linear finite-volume schemes. This new approach provides a better characterization of the accuracy of a spatio-temporal discretization as a function of wavenumber and frequency, allowing us to distinguish between behavior at long wavelengths, where accuracy is a prime concern, and short wavelengths, where stability concerns are of greater importance. We use this analysis to develop a specialized third-order Runge Kutta scheme that minimizes the temporal integration error in the discrete structure factor at long wavelengths for the one-dimensional linearized LLNS equations. Together with a novel method for discretizing the stochastic stress tensor in dimension larger than one, our improved temporal integrator yields a scheme for the three-dimensional equations that satisfies a discrete fluctuation-dissipation balance for small time steps and is also sufficiently accurate even for time steps close to the stability limit.

## I Introduction

Recently the fluid dynamics community has considered increasingly complex physical, chemical, and biological phenomena at the microscopic scale, including systems for which significant interactions occur across multiple scales. At a molecular scale, fluids are not deterministic; the state of the fluid is constantly changing and stochastic, even at thermodynamic equilibrium. As simulations of fluids push toward the microscale, these random thermal fluctuations play an increasingly important role in describing the state of the fluid, especially when investigating systems where the microscopic fluctuations drive a macroscopic phenomenon such as the evolution of instabilities, or where the thermal fluctuations drive the motion of suspended microscopic objects in complex fluids. Some examples in which spontaneous fluctuations can significantly affect the dynamics include the breakup of droplets in jets Moseler:00 (); Eggers:02 (); Kang:07 (), Brownian molecular motors Astumian:02 (); Oster:02 (); Broeck:04 (); Meurs:04 (), Rayleigh-Bernard convection (both single species Wu:95 () and mixtures Quentin:95 (), Kolmogorov flows Bena:00 (); Bena:99 (); Mansour:99 (), Rayleigh-Taylor mixing Kadau:04 (); FluidMixing_DSMC (), combustion and explosive detonation Nowakowski:03 (); Lemarchand:04 (), and reaction fronts Moro:04 ().

Numerical schemes based on a particle representation of a fluid (e.g., molecular dynamics, direct simulation Monte Carlo DSMCReview_Garcia ()) inherently include spontaneous fluctuations due to the irregular dynamics of the particles. However, by far the most common numerical schemes in computational fluid dynamics are based on solving partial differential equations. To incorporate thermal fluctuations into macroscopic hydrodynamics, Landau and Lifshitz introduced an extended form of the compressible Navier-Stokes equations obtained by adding white-noise stochastic flux terms to the standard deterministic equations. While they were originally developed for equilibrium fluctuations, specifically the Rayleigh and Brillouin spectral lines in light scattering, the validity of the Landau-Lifshitz Navier-Stokes (LLNS) equations for non-equilibrium systems has been assessed LLNS_Espanol () and verified in molecular simulations Garcia:91 (); Mansour:87 (); Mareschal:92 (). The LLNS system is one of the more complex examples in a broad family of PDEs with stochastic fluxes. Many members of this family arise from the LLNS equations in a variety of approximations (e.g., stochastic heat equation) while others are stochastic variants of well-known PDEs, such as the stochastic Burger’s equation Bell:06 (), which can be derived from the continuum limit of an asymmetric excluded random walk.

Several numerical approaches for fluctuating hydrodynamics have been proposed. The earliest work by Garcia et al. Garcia:87 () developed a simple scheme for the stochastic heat equation and the linearized one-dimensional LLNS equations. Ladd et al. have included stress fluctuations in (isothermal) Lattice-Boltzmann methods for some time Ladd:93 (), and recently a better theoretical foundation has been established LB_ThermalFluctuations (); LB_StatMech (). Moseler and Landman Moseler:00 () included the stochastic stress tensor of the LLNS equations in the lubrication equations and obtain good agreement with their molecular dynamics simulation in modeling the breakup of nano-jets. Sharma and Patankar Sharma:04 () developed a fluid-structure coupling between a fluctuating incompressible solver and suspended Brownian particles. Coveney, De Fabritiis, Delgado-Buscalioni and co-workers have also used the isothermal LLNS equations in a hybrid scheme, coupling a continuum fluctuating solver to a Molecular Dynamics simulation of a liquid FluctuatingHydro_Coveney (); FluctuatingHydroMD_Coveney (); FluctuatingHydroHybrid_MD (). Atzberger and collaborators AtzbergerETAL:2007 () have developed a version of the immersed boundary method that includes fluctuations in a pseudo-spectral method for the incompressible Navier-Stokes equations. Voulgarakis and Chu StagerredFluctHydro () developed a staggered scheme for the isothermal LLNS equations as part of a multiscale method for biological applications, and a similar staggered scheme was also described in Ref. Delgado:08 ().

Recently, Bell et al. Bell:07 () introduced a centered scheme for the LLNS equations based on interpolation schemes designed to preserve fluctuations combined with a third-order Runge-Kutta (RK3) temporal integrator. In that work, the principal diagnostic used for evaluation of the numerical method was the accuracy of the local (cell) variance and spatial (cell-to-cell) correlation structure for equilibrium and selected non-equilibrium scenarios (e.g., constant temperature gradient). The metric established by those types of tests is, in some sense, simultaneously too crude and too demanding. It is too crude in the sense that it provides only limited information from detailed simulations that cannot be directly linked to specific properties of the scheme. On the other hand, such criteria are too demanding in the sense that they place requirements on the discretization integrated over all wavelengths, requiring that the method perform well at high wavenumbers where a deterministic PDE solver performs poorly. Furthermore, although Bell et al. Bell:07 () demonstrate that RK3 is an effective algorithm, compared with other explicit schemes for the compressible Navier-Stokes equations, the general development of schemes for the LLNS equations has been mostly trial-and-error.

Here, our goal is to establish a more rational basis for the analysis and development of explicit finite volume scheme for SPDEs with a stochastic flux. The approach is based on analysis of the structure factor (equilibrium fluctuation spectrum) of the discrete system. The structure factor is, in essence, the stationary spatio-temporal correlations of hydrodynamic fluctuations as a function of spatial wavenumber and temporal frequency; the static structure factor is the integral over frequency (i.e., the spatial spectrum). By analyzing the structure factor for a numerical scheme, we are able to develop notions of accuracy for a given discretization at long wavelengths. Furthermore, in many cases the theoretical analysis for the structure factor is tractable (with the aid of symbolic manipulators) allowing us to determine optimal coefficients for a given numerical scheme. We perform this optimization as a two-step procedure. First, a spatial discretization is developed that satisfies a discrete form of the fluctuation-dissipation balance condition. Then, a stable temporal integrator is proposed and the covariances of the random numbers are chosen so as to maximize the order of temporal accuracy of the small-wavenumber static structure factor. We focus primarily on explicit schemes for solving the LLNS equations because even at the scales where thermal fluctuations are important, the limitation on time step imposed by stability is primarily due to the hyperbolic terms. That is, when the cell size is comparable to the length scale for molecular transport (e.g., mean free path in a dilute gas) the time step for these compressible hydrodynamic equations is limited by the acoustic CFL condition. At even smaller length scales the viscous terms further limit the time step yet the validity of a continuum representation for the fluid starts to break down at those atomic scales.

The paper is divided into roughly two parts: The first half (sections II-IV) defines notation, develops the formalism, and derives the expressions for analyzing a general class of linear stochastic PDEs from the LLNS family of equations. The main result in the first half, how to evaluate the structure factor for a numerical scheme, appears in section III.2. The second half applies this analysis to systems of increasing complexity, starting with the stochastic heat equation (section V.1), followed by the LLNS system in one dimension (section VI) and three dimensions (section VII). The paper closes with a summary and concluding remarks.

## Ii Landau-Lifshitz Navier-Stokes Equations

We consider the accuracy of explicit finite-volume methods for solving the Landau-Lifshitz Navier-Stokes (LLNS) system of stochastic partial differential equations (SPDEs) in dimensions, given in conservative form by

 ∂tU=−∇⋅[F(U)−Z(U,r,t)], (1)

where is a vector of conserved variables that are a function of the spatial position and time . The conserved variables are the densities of mass , momentum , and energy , expressed in terms of the primitive variables, mass density , velocity and temperature ; here is the internal energy density. The deterministic flux is taken from the traditional compressible Navier-Stokes-Fourier equations and can be split into hyperbolic and diffusive fluxes:

 F(U)=FH(U)+FD(U),

where

 FH=⎡⎢⎣ρvρvvT+PI(e+P)v⎤⎥⎦ and FD=−⎡⎢⎣0σσ⋅v+ξ⎤⎥⎦,

is the pressure, the viscous stress tensor is for (we have assumed zero bulk viscosity) and for , and the heat flux is . We denote the adjoint (conjugate transpose) of a matrix or linear operator with . As postulated by Landau-Lifshitz Landau:Fluid (); LLNS_Espanol (), the stochastic flux

 Z=⎡⎢⎣0ΣΣ⋅v+Ξ⎤⎥⎦

is composed of the stochastic stress tensor and stochastic heat flux vector , assumed to be mutually uncorrelated random Gaussian fields with a covariance

 ⟨Σ(r,t)Σ⋆(r′,t′)⟩= CΣδ(t−t′)δ(r−r′), where C(Σ)ij,kl=2¯ηkB¯¯¯¯T(δikδjl+δilδjk−2dfδijδkl) ⟨Ξ(r,t)Ξ⋆(r′,t′)⟩= CΞδ(t−t′)δ(r−r′), where C(Ξ)i,j=2¯μkB¯¯¯¯T2δij (2)

In the LLNS system, the hyperbolic or advective fluxes are responsible for transporting the conserved quantities at the speed of sound or fluid velocity, without dissipation. On the other hand, the diffusive or dissipative fluxes are the ones responsible for damping the thermal fluctuations generated by the stochastic or fluctuating fluxes. At equilibrium a steady state is reached in which a fluctuation-dissipation balance condition is satisfied.

In the original formulation, Landau and Lifshitz only considered adding stochastic fluxes to the linearized Navier-Stokes equations, which leads to a well-defined system of SPDEs whose equilibrium solutions are random Gaussian fields. Derivations of the equations of fluctuating hydrodynamics through careful asymptotic expansions of the underlying microscopic (particle) dynamics give equations for the Gaussian fluctuations around the solution to the usual deterministic Navier-Stokes equations LebowitzHydroReview (), in the spirit of the Central Limit Theorem. Therefore, numerical solutions should, in principle, consist of two steps: First solving the nonlinear deterministic equations for the mean solution, and then solving the linearized equations for the fluctuations around the mean. If the fluctuations are small perturbations, it makes sense numerically to try to combine these two steps into one and simply consider non-linear equations with added thermal fluctuations. There is also hope that this might capture effects not captured in the two-system approach, such as the effect of fluctuations on the very long-time dynamics of the mean (e.g., shock drift Bell:07 ()) or hydrodynamic instabilities Wu:95 (); Moseler:00 (); FluidMixing_DSMC ().

The linearized equations of fluctuating hydrodynamics can be given a well-defined interpretation with the use of generalized functions or distributions DaPratoBook (). However, the non-linear fluctuating hydrodynamic equations (1) must be treated with some care since they have not been derived from first-principles LLNS_Espanol () and are in fact mathematically ill-defined due to the high irregularity of white-noise fluctuating stresses GardinerBook (). More specifically, because the solution of these equations is itself a distribution the interpretation of the nonlinear terms requires giving a precise meaning to products of distributions, which cannot be defined in general and requires introducing some sort of regularization. Although written formally as an SPDE, the LLNS equations are usually interpreted in a finite volume context, where the issues of regularity, at first sight, disappear. However, in finite volume form the level of fluctuations becomes increasingly large as the volume shrinks and the non-linear terms diverge leading to an “ultraviolet catastrophe” of the kind familiar in other fields of physics. Furthermore, because the noise terms are Gaussian, it is possible for rare events to push the system to states that are not thermodynamically valid such as negative or . For that reason, we will focus on the linearized LLNS equations, which can be given a well-defined interpretation. Since the fluctuations are expected to be a small perturbation of the deterministic solution, the nonlinear equations should behave similarly to the linearized equations anyway, at least near equilibrium for sufficiently large cells.

To simplify the exposition we assume the fluid to be a mono-atomic ideal gas; the generalization of the results for an arbitrary fluid is tedious but straightforward. For an ideal gas the equation of state may be written as , where is the isothermal speed of sound. The internal energy density is , where is the heat capacity at constant volume, which may be written as where is the number of degrees of freedom of the molecules (for monoatomic gases there are translational degrees of freedom), and is the heat capacity at constant pressure. For analytical calculations it is convenient to convert the LLNS system from conserved variables to primitive variables, since the primitive variables are uncorrelated at equilibrium and the equations (1) simplify considerably,

 Dtρ= ∇ ⋅(ρv) ρ(Dtv)= −∇P+∇⋅(σ+Σ) ρcp(DtT)= DtP+∇⋅(ξ+Ξ)+(σ+Σ):∇v, (3)

where denotes the familiar advective derivative. Note that in the fully non-linear numerical implementation, however, we continue to use the conserved variables to ensure that the physical conservation laws are strictly obeyed.

Linearizing (3) around a reference uniform equilibrium state , , , and dropping the deltas for notational simplicity,

 U=⎡⎢⎣δρδvδT⎤⎥⎦→⎡⎢⎣ρvT⎤⎥⎦,

we obtain the linearized LLNS system for the equilibrium thermal fluctuations,

 ∂tU=−∇⋅[FU−Z]=−∇⋅[FHU+FD∇U−Z], (4)

where

 FHU=⎡⎢ ⎢⎣ρ0v+ρv0(c20ρ−10ρ+c20T−10T)I+v0vTc20c−1vv+Tv0⎤⎥ ⎥⎦ % and FD∇U=⎡⎢ ⎢⎣c0ρ−10η0¯¯¯¯¯∇vρ−10c−1vμ0∇T⎤⎥ ⎥⎦,

and denotes a symmetrized traceless gradient, . Here is a random Gaussian field with a covariance

 ⟨Z(r,t)Z⋆(r′,t′)⟩=CZδ(t−t′)δ(r−r′),

where the covariance matrix is block diagonal,

 CZ=⎡⎢ ⎢⎣0000ρ−20CΣ000ρ−20c−2vCΞ⎤⎥ ⎥⎦,

and and are given in Eq. (2). Equation (4) is a system of linear SPDEs with additive noise that can be analyzed within a general framework, as we develop next. We note that the stochastic “forcing” in (4) is essentially a divergence of white noise, modeling conservative intrinsic (thermal) fluctuations LebowitzHydroReview (), rather than the more common external fluctuations modeled through white noise forcing StochHeatEq_Weak ().

The next two sections develop the tools for analyzing explicit finite volume schemes for linearized SPDEs, such as the LLNS system, specifically how to predict the equilibrium spectrum of the fluctuations (i.e., structure factor) from the spatial and temporal discretization used by the numerical algorithm. These analysis tools are demonstrated for simple examples in Section V.1 and applied to the LLNS system in Sections VI and VII.

## Iii Explicit Methods for Linear Stochastic Partial Differential Equations

In this section, we develop an approach for analyzing the behavior of explicit discretizations for a broad class of SPDEs, motivated by the linearized form of the LLNS equations. In particular, we consider a general linear SPDE for the stochastic field of the form

 dU(t)=LU(t)dt+KdB(t), (5)

with periodic boundary conditions on the torus , where (the generator) and (the filter) are time-independent linear operators, and is a cylindrical Wiener process (Brownian sheet), and the initial condition at is . As common in the physics literature, we will abuse notation and often write

 ∂tU=LU+KW,

where is spatio-temporal white noise, i.e., a random Gaussian field with zero mean and covariance

 ⟨W(r,t)W⋆(r′,t′)⟩=δ(t−t′)δ(r−r′). (6)

The so-called mild solution DaPratoBook () of (5) is a generalized process

 U(t)=etLU0+∫t0e(t−s)LKdB(s), (7)

where the integral denotes a stochastic convolution. If the operator is dissipative, that is, for all , then at long times the solution to (5) is a Gaussian process with mean zero and covariance

 CU(t)=⟨U(t′)U⋆(t′+t)⟩=∫0−∞e−sLKK⋆e(t−s)L⋆ds,t≥0. (8)

This means that (5) has a unique invariant measure (equilibrium or stationary distribution) that is Gaussian with mean zero and covariance given in Eq. (8).

In general, the field is only a generalized function of the spatial coordinate and cannot be evaluated pointwise. For the cases we will consider here, specifically, translationally-invariant problems where and are differential operators, this difficulty can be avoided by transforming (5) to Fourier space via the Fourier series transform

 U(r,t)= ∑k∈ˆVeik⋅rˆU(k,t) (9) ˆU(k,t)=1V ∫r∈Ve−ik⋅rU(r,t)dr, (10)

where is the volume of the system, and each wavevector is expressed in terms of the integer wave index , giving the set of discrete wavevectors

 ˆV={k=2πκ/H|κ∈Zd}.

In Fourier space, the SPDE (5) becomes an infinite system of uncoupled stochastic ordinary differential equations (SODEs),

 dˆU(t)=ˆLˆU(t)dt+ˆKdˆB(t), (11)

one SODE for each . The invariant distribution of (11) is a zero-mean Gaussian random process, characterized fully by the covariance obtained from the spatial Fourier transform of (8),

 S(k,t)=V⟨ˆU(k,t′)ˆU⋆(k,t′+t)⟩=12π∫∞−∞eiωtS(k,ω)dω, (12)

where the dynamic structure factor (space-time spectrum) is

 S(k,ω)=V⟨ˆU(k,ω)ˆU⋆(k,ω)⟩=(ˆL−iω)−1(ˆKˆK⋆)(ˆL⋆+iω)−1, (13)

which follows directly from the space-time Fourier transform of the SPDE (5). By integrating the dynamic spectrum over all frequencies , one gets the static structure factor

 S(k)=S(k,t=0)=12π∫∞−∞S(k,ω)dω, (14)

which is the spatial spectrum of an equilibrium snapshot of the fluctuating field and is the Fourier equivalent of . Note that the static structure factor of spatial white noise (a snapshot of ) is unity independent of the wavevector, .

### iii.1 Discretization

For the types of equations we will consider in this paper, the invariant measure is spatially white, specifically, is diagonal and independent of . The associated fluctuating field cannot be evaluated pointwise, therefore, it is more natural to use finite-volume cell averages, denoted here by . In the deterministic setting, for uniform periodic grids there is no important difference between finite-volume and finite-difference methods. Our general approach can likely be extended also to analysis of stochastic finite-element discretizations, however, such methods have yet to be developed for the LLNS equations and here we focus on finite-volume methods. For notational simplicity, we will discuss problems in one spatial dimension (), with (mostly) obvious generalizations to higher dimensions.

Space is discretized into identical cells of length , and the value stored in cell is the average of the corresponding variable over the cell

 Uj(t)=1Δx∫jΔx(j−1)ΔxU(x,t)dx. (15)

Time is discretized with a time step , approximating cell averages of pointwise in time with ,

 Unj≈Uj(nΔt),

where enumerates the time steps. The white noise cannot be evaluated pointwise in either space or time and is discretized using a spatio-temporal average

 ¯¯¯¯¯¯¯Wnj(t)=1ΔxΔt∫(n+1)ΔtnΔt∫jΔx(j−1)ΔxW(x,t)dxdt, (16)

which is a normal random variable with zero mean and variance , independent between different cells and time steps. Note that for certain types of equations the dynamic structure factor may be white in frequency as well. In this case, a pointwise-in-time discretization is not appropriate and one can instead use a spatio-temporal average as done for white noise in (16).

We will study the accuracy of explicit linear finite-volume schemes for solving the SPDE (5). Rather generally, such methods are specified by a linear recursion of the form

 Un+1=(I+LΔt)Un+√ΔtΔxKWn, (17)

where and are consistent stencil discretizations of the continuum differential operators and (note that and may involve powers of in general). Here

 Wn=(ΔxΔt)12¯¯¯¯¯¯¯Wn (18)

is a vector of standard normal variables with mean zero and variance one.

Without the random forcing, the deterministic equation and the associated discretization can be studied using classical tools and notions of stability, consistency, and convergence. Under the assumption that the discrete generator is dissipative, the initial condition will be damped and the equilibrium solution will simply be a constant. The addition of the random forcing, however, leads to a non-trivial invariant measure (equilibrium distribution) of determined by an interplay between the (discretized) fluctuations and dissipation. Because of the dissipative nature of the generator, any memory of the initial condition will eventually disappear and the long time dynamics is guaranteed to follow an ergodic trajectory that samples the unique invariant measure. In order to characterize the accuracy of the stochastic integrator, we will analyze how well the discrete invariant measure (equilibrium distribution) reproduces the invariant measure of the continuum SPDE (this is a form of weak convergence). Note that due to ergodicity, ensemble averages can either be computed by averaging the power spectrum of the fields over multiple samples or averaging over time (after sufficiently many initial equilibration steps). In the theory we will consider the limit and then average over different realizations of the noise to obtain the discrete structure factors. In numerical calculations, we perform temporal averaging.

Regardless of the details of the iteration (17), will always be a Gaussian random vector generated anew at each step using a random number generator. The discretized field is therefore a linear combination of Gaussian variates and it is therefore a Gaussian vector-valued stochastic process. In particular, the invariant measure (equilibrium distribution) of is fully characterized by the covariance

 C(U)j,j′,n=limNs→∞⟨UNsj(UNs+nj′)⋆⟩, (19)

which we would like to compare to the covariance of the continuum Gaussian field given by (8). This comparison is best done in the Fourier domain by using the spatial discrete Fourier transform, defined for a spatially-discrete field [for example, or ] via

 Uj= ∑k∈ˆVdˆUkeijΔk (20) ˆUk= 1VNc−1∑j=0Uj+1e−ijΔkΔx, (21)

where we have denoted the discrete dimensionless wavenumber , and the wave index is now limited to the first values,

 ˆVd={k=2πκ/H|0≤κ

Since the fields are real-valued, there is a redundancy in the Fourier coefficients because of the Hermitian symmetry between and (essentially, the second half of the wave indices correspond to negative ), and thus we will only consider , giving a (Nyquist) cutoff wavenumber .

What we would like to compare is the Fourier coefficients of the numerical approximation, , with the Fourier coefficients of the continuum solution, . The invariant measure of has zero mean and is characterized by the covariance obtained from the spatial Fourier transform of (19),

 Sk,n=VlimNs→∞⟨ˆUNsk(ˆUNs+nk)⋆⟩. (22)

From the definition of the discrete Fourier transform it follows that for small , i.e., smooth Fourier basis functions on the scale of the discrete grid, converges to the Fourier coefficient of the continuum field. Therefore, is the discrete equivalent (numerical approximation) to the continuum structure factor . We define a discrete approximation to be weakly consistent if

 limΔx,Δt→0Sk,n=⌊t/Δt⌋=S(k,t),

for any chosen and . This means that, given a sufficiently fine discretization, the numerical scheme can accurately reproduce the structure factor for a desired wave index and time lag. An alternative view is that a convergent scheme reproduces the slow (compared to ) and large-scale (compared to ) fluctuations, that is, it accurately reproduces the dynamic structure factor for small and . Our goal here is to quantify this for several numerical methods for solving stochastic conservation laws and optimize the numerical schemes by tuning parameters to obtain the best possible approximation to for small and .

Much of our analysis will be focused on the discrete static structure factor

Note that for a spatially-white field , the finite-volume averages are independent Gaussian variates with mean zero and variance , and the discrete Fourier coefficients are independent Gaussian variates with mean zero and variance . As a measure of the accuracy of numerical schemes for solving Eq. (5), we will compare the discrete static structure factors with the continuum prediction , for all of the discrete wavenumbers (i.e., pointwise in Fourier space). It is expected that any numerical scheme will produce some artifacts at the largest wavenumbers because of the strong corrections due to the discretization; however, small wavenumbers ought to have much smaller errors because they evolve over time scales and length scales much larger than the discretization step sizes. Specifically, we propose to look at the series expansions

 Sk−S(k)=O(Δtp1kp2)

and optimize the numerical schemes by maximizing the powers and . Next we describe the general formalism used to obtain explicit expressions for the discrete structure factors for a general explicit method, and then illustrate the formalism on some simple examples, before attacking the more complex equations of fluctuating hydrodynamics.

### iii.2 Analysis of Linear Explicit Methods

Regardless of the details of a particular scheme and the particular linear SPDE being solved, at the end of the timestep a typical explicit scheme makes a linear combination of the values in the neighboring cells and random variates to produce an updated value,

 Un+1j=Unj+Δj=wD∑Δj=−wDΦΔjUnj+Δj+Δj=wS∑Δj=−wSΨΔjWnj+Δj, (23)

where and are the deterministic and stochastic stencil widths. The particular forms of the matrices of coefficients and depend on the scheme, and will involve powers of and . Here we assume that for each the random increment is an independent vector of normal variates with covariance constant for all of the cells and thus wavenumbers, where is the total number of random numbers utilized per cell per stage. Computer algebra systems can be used to obtain explicit formulas for the matrices in (23); we have made extensive use of Maple for the calculations presented in this paper.

Assuming a translation invariant scheme, the iteration (23) can easily be converted from real space to an iteration in Fourier space,

 ˆUn+1k=ˆUnk+Δj=wD∑Δj=−wDΦΔjˆUnkexp(iΔjΔk)+Δj=wS∑Δj=−wSΨΔjˆWnkexp(iΔjΔk), (24)

where different wavenumbers are not coupled to each other. In general, any linear explicit method can be represented in Fourier space as a recursion of the form

 ˆUn+1k=MkˆUnk+NkˆWnk, (25)

where the explicit form of the matrices and depend on the particular scheme and typically contain various powers of , , and , and . By iterating this recurrence relation, we can easily obtain (assuming )

 ˆUn+1k=n∑l=0(Mk)lNkˆWn−lk,

from which we can calculate

 Snk=V⟨(ˆUnk)(ˆUnk)⋆⟩=n−1∑l=0(Mk)l(ΔxNkCWN⋆k)(M⋆k)l=n−1∑l=0(Mk)l˜C(M⋆k)l.

In order to calculate this sum explicitly, we will use the following identity

to obtain a linear system for the entries of the matrix . If the deterministic method is stable, which means that all eigenvalues of the matrix are below unity for all wavenumbers, then in the limit the first term on the right hand side will vanish, to give

 MkSkM⋆k−Sk=−ΔxNkCWN⋆k. (26)

This is a linear system of equations for the equilibrium static structure factor produced by a given scheme, where the number of unknowns is equal to the square of the number of variables (field components). By simply deleting the subscripts one obtains a more general but much larger linear system MultigridMC_Goodman () for the real space equilibrium covariance of a snapshot of the discrete field ,

 MCUM⋆−CU=−ΔxNC(Nc)WN⋆,

where is the covariance matrix of the random increments. Note that this relation continues to hold even for schemes that are not translation invariant such as generalizations to non-periodic boundary conditions; however, the number of unknowns is now the square of the total number of degrees of freedom so that explicit solutions will in general not be possible. Based on standard wisdom for deterministic schemes, it is expected that schemes that perform well under periodic boundary conditions will also perform well in the presence of boundaries when the discretization is suitably modified only near the boundaries.

A similar approach to the one illustrated above for the static structure factor can be used to evaluate the discrete dynamic structure factor

 Sk,ω=limNs→∞V(NsΔt)⟨ˆUNsk,ω(ˆUNsk,ω)⋆⟩

from the time-discrete Fourier transform

 ˆUNsk,ω=1NsNs∑l=0exp(−ilΔω)ˆUlk,

where , and the frequency is less than the Nyquist cutoff, . The calculation yields

 Sk,ω=[I−exp(−iΔω)Mk]−1(ΔxΔtNkCWN⋆k)[I−exp(iΔω)M⋆k]−1. (27)

Equation (27) can be seen as discretized forms of the continuum version (13) in the limits , (the corresponding correlations in the time-domain are given in Ref. MultigridMC_Goodman ()).

Equations (26) and (27) are the main result of this section and we have used it to obtain explicit expressions for and for several equations and schemes. In the next sections we will illustrate the above formalism for several simple examples of stochastic conservation laws.

#### iii.2.1 Discrete Fluctuation-Dissipation Balance

Let us first consider the static structure factors for very small time steps. In the limit , temporal terms of order two or more can be ignored so that all time-integration methods behave like an explicit first-order Euler iteration as in (17),

 ˆUn+1k=(I+ΔtˆL(0)k)ˆUnk+√ΔtΔxˆK(0)kˆWk, (28)

where can be thought of as the spatial discretization of the generator , and is the spatial discretization of the filtering operator . Comparing to (25) we can directly identify and and substitute these into Eq. (26). Keeping only terms of order on both sides we obtain the condition

 ˆL(0)kS(0)k+S(0)k(ˆL(0)k)⋆=−ˆK(0)kCW(ˆK(0)k)⋆, (29)

where (see also a related real-space derivation using Ito’s calculus in Ref. AMR_ReactionDiffusion_Atzberger (), as well as Section VIII in Ref. MultigridMC_Goodman ()). It can be shown that if is definite, Eq. (29) has a unique solution. Assuming that is as given in Eq. (18), i.e., that , and that the spatial discretizations of the generator and filter operators satisfy a discrete fluctuation-dissipation balance

 ˆL(0)k+(ˆL(0)k)⋆=−ˆK(0)k(ˆK(0)k)⋆, (30)

we see that is the solution to equation (29), that is, at equilibrium the discrete fields are spatially-white. The discrete fluctuation-dissipation balance condition can also be written in real space,

 L(0)+(L(0))⋆=K(0)(K(0))⋆. (31)

The condition (31) is the discrete equivalent of the continuum fluctuation-dissipation balance condition FluctuationDissipation_Kubo ()

 L+L⋆=−KK⋆, (32)

which ensures that , i.e., that the invariant measure of the SPDE is spatially-white. We observe that adding a skew adjoint component to does not alter the fluctuation-dissipation balance above, as is the case with non-dissipative (advective) terms. Numerous equations LebowitzHydroReview () modeling conservative thermal systems satisfy condition (32), including the linearized LLNS equations (with some additional prefactors). In essence, the fluctuations injected at all scales by the spatially white forcing are filtered by and then dissipated by at just equal rates.

## Iv Linear Stochastic Conservation Laws

The remainder of this paper is devoted to the study of the accuracy of finite-volume methods for solving linear stochastic PDEs in conservation form,

 ∂tU=−∇⋅[(AU−C∇U)−EW], (33)

where , and are constants, and is Gaussian spatio-temporal white noise. The white noise forcing and its divergence here need to be interpreted in the (weak) sense of distributions since they lack the regularity required for the classical definitions. The linearization of the LLNS equations (1) leads to a system of the form (33), as do a number of other classical PDEs LebowitzHydroReview (), such as the stochastic advection-diffusion equation

 ∂tT=−a⋅∇T+μ∇2T+√2μ∇⋅W, (34)

where is a scalar stochastic field, is the advective velocity, , is the diffusion coefficient, and . The simplest case is the stochastic heat equation, obtained by taking .

A key feature the type of system considered here is that the noise is intrinsic to the system and appears in the flux as opposed to commonly treated systems that include an external stochastic forcing term, such as the form of a stochastic heat equation considered in Ref. StochHeatEq_Weak (). Since white noise is more regular than the spatial derivative of white noise, external noise leads to more regular equilibrium fields (e.g., continuous functions in one dimension). Intrinsic noise, on the other hand, leads to very irregular equilibrium fields. Notationally, it is convenient to write (33) as,

 ∂tU=−D(AU−CGU−EW), (35)

defining the divergence and gradient operators, . In the types of equations that appear in hydrodynamics, such as the LLNS equations, the operator is skew-adjoint, (hyperbolic or advective flux), (dissipative or diffusive flux), and , i.e., . Therefore, the generator and filter satisfy the fluctuation-dissipation balance condition (32) and the equilibrium distribution is spatially-white. Note that even though advection makes some of the eigenvalues of complex, the generator is dissipative and (33) has a unique invariant measure because the real part of all of the eigenvalues of is negative except for the unique zero eigenvalue.

It is important to point out that discretizations of the continuum operators do not necessarily satisfy the discrete fluctuation-dissipation condition (31). One way to ensure the condition is satisfied is to discretize the diffusive components of the generator and the filter using a discrete divergence and discrete gradient so that the discrete fluctuation-dissipation balance condition holds. If, however, the discretization of the advective component of the generator is not skew-adjoint, this can perturb the balance (30). Notably, various upwinding methods lead to discretizations that are not skew-adjoint. The correction to the structure factor due to a non-zero can easily be obtained from Eq. (29), and in one dimension the result is simply

 ΔS(0)k=−ΔL(A)kL(D)k+ΔL(A)k. (36)

We will use centered differences for the advective generator in this work, which ensures a skew-adjoint , and our focus will therefore be on satisfying the discrete fluctuation-dissipation balance for the diffusive terms.

### iv.1 Finite-Volume Numerical Schemes

We consider here rather general finite-volume methods for solving the linear SPDE (33) in one dimension,

 ∂tU=−∂∂x[F(U)−Z]=−∂∂x[(A−C∂∂x)U−EW] (37)

with periodic boundaries, where we have denoted the modified (potentially correlated) white noise flux with . As for classical finite-volume methods for the deterministic case, we start from the PDE and integrate the left and right hand sides over a given cell over a given time step , and use integration by parts to obtain the formally exact

 (38)

where the deterministic discrete fluxes and dimensionless discrete stochastic fluxes are calculated on the boundaries of the cells (points in one dimension, edges in two dimensions, and faces in three dimensions), indexed here with half-integers. These fluxes represent the total rate of transport through the interface between two cells over a given finite time interval , and (38) is nothing more than a restatement of conservation. The classical interpretation of pointwise evaluation of the fluxes is not appropriate because white noise forcing lacks the regularity of classical smooth forcing and cannot be represented in a finite basis. Instead, just as we projected the fluctuating fields using finite-volume averaging, we ought to project the fluxes to a finite representation as well through spatio-temporal averaging, as done in Eq. (16). For the purposes of our analysis, one can simply think of the discrete fluxes as an approximation that has the same spectral properties as the corresponding continuum Gaussian fields over the wavevectors and frequencies represented by the finite discretization.

The goal of numerical methods is to approximate the fluxes as best as possible. In general, within each time step of a scheme there may be stages or substeps; for example, in the classic MacCormack method there is a predictor and a corrector stage (), and in the three-stage Runge-Kutta method of Williams et al. Bell:07 () there are three stages (). Each stage is of the conservative form (38),

 Un+sNstj=s−1∑s′=0α(s)s′Un+s′Nstj−ΔtΔx(F(s)j+12−F(s)j−12)+Δt1/2Δx3/2(Z(s)j+12−Z(s)j−12), (39)

where the ’s are some coefficients, , and each of the stage fluxes are partial approximations of the continuum flux. For the stochastic integrators we discuss here, the deterministic fluxes are calculated the same way as they would be in the corresponding deterministic scheme. In general, the stochastic fluxes can be expressed in terms of independent unit normal variates that are sampled using a random number generator. The stochastic fluxes in each stage may be the same, may be completely independent, or they may have non-trivial correlations between stages.

Note that it is possible to avoid non-integer indices by re-indexing the fluxes in Eq. (38) and writing it in a form consistent with (23),

 (40)

However, when considering the order of accuracy of the stencils and also fluctuation-dissipation balance in higher dimensions, it will become important to keep in mind that the fluxes are evaluated on the faces (edges or half-grid points) of the grid, and therefore we will keep the half-integer indices. Note that for face-centered values, such as fluxes, it is best to add a phase factor in the definition of the Fourier transform, even though such pure phase shifts will not affect the correlation functions and structure factors.

Before we analyze schemes for the complex LLNS equations, we present an illustrative explicit calculation for the one-dimensional stochastic heat equation.

## V Example: Stochastic Heat Equation

We now illustrate the general formalism presented in Section IV for the simple case of an Euler and predictor-corrector scheme for solving the stochastic heat equation in one dimension,

 υt=μυxx+√2μWx, (41)

where is a scalar field and is the mass or heat diffusion coefficient. The solution in the Fourier domain is trivial, giving

 S(k,ω)=2μk2ω2+μ2k4, and S(k)=1. (42)

### v.1 Static Structure Factor

We first study a simple second-order spatial discretization of the dissipative fluxes

 Fj+12=μΔx(uj+1−uj),

combined with an Euler integration in time, to give a simple numerical method for solving the SPDE (41),

 un+1j=unj+μΔtΔx2(unj−1−2unj+unj+1)+√2μΔt1/2Δx3/2(Wnj+12−Wnj−12), (43)

where and the ’s are independent unit normal random numbers with zero mean generated anew at every time step (here ). From (43), we can extract the recursion coefficients appearing in (25),

 Mk=1+β(e−iΔk−2+eiΔk)=1+2β(cosΔk−1),
 Nk=√2μΔt1/2Δx3/2(eiΔk/2−e−iΔk/2),

where

 β=μΔtΔx2

denotes a dimensionless diffusive time step (ratio of the time step to the diffusive CFL limit). Together with , (26) becomes a scalar equation for the discrete structure factor,

 (MkM⋆k−1)Sk=−ΔxNkN⋆k,

with dimensionless solution

 (44)

The time-dependent result can also easily be derived from (27),

 Snk=(1−e−t/τ)Sk, where t=nΔt

where is the familiar relaxation time for wavenumber , showing that the smallest wavenumbers take a long time to reach the equilibrium distribution.

Equation (44) is a vivid illustration of the typical result for schemes for stochastic transport equations based on finite difference stencils, also shown in Fig. 1. Firstly, we see that for small we have that , showing that the smallest wavenumbers are correctly handled by the discretization for any time step. Also, this shows that the error in the structure factor is of order , i.e., of order , as expected for the Euler scheme, whose weak order of convergence is one for SODEs. Finally, it shows that the error grows quadratically with (from symmetry arguments, only even powers will appear). By looking at the largest wavenumber, , we see that , from which we instantly see the CFL stability condition , which guarantees that the structure factor is finite and positive for all . Furthermore, we see that for , the structure factor is approximately unity for all wavenumbers. That is, a sufficiently small step will indeed reproduce the proper equilibrium distribution.

By contrast, a two-stage predictor-corrector scheme for the diffusion equation,

 ~unj= unj+μΔtΔx2(unj−1−2unj+unj+1)+√2μΔt1/2Δx3/2(Wnj+12−Wnj−12) (predictor) un+1j= 12[unj+~unj+μΔtΔx2(~unj−1−2~unj+~unj+1)+√2μΔt1/2Δx3/2(Wnj+12−Wnj−12)] (corrector), (45)

achieves much higher accuracy, namely, a structure factor that deviates from unity by a higher order in both and ,

 PC-1RNG: Sk≈1−β2Δk4/4,

as illustrated in Fig. 1. We can also use different stochastic fluxes in the predictor and the corrector schemes (i.e., use random numbers per cell per stage), with an added pre-factor of to compensate for the variance reduction of the averaging between the two stages,

 ~unj= unj+μΔtΔx2(unj−1−2unj+unj+1)+2√μΔt1/2Δx3/2(W(n,P)j+12−W(n,P)j−12) (% predictor) un+1j= 12[unj+~unj+μΔtΔx2(~unj−1−2~unj+~unj+1)+2√μΔt1/2Δx3/2(W(n,C)j+12−W(n,C)j−12)] (corrector). (46)

For the scheme (46) the analysis reveals an even greater spatio-temporal accuracy of the static structure factors, namely, third order temporal accuracy

 PC-2RNG: Sk≈1+β3Δk6/8.

This illustrates the importance of the handling of the stochastic fluxes in multi-stage algorithms, as we will come back to shortly. Note that the analysis we presented here for explicit methods can easily be extended to implicit and semi-implicit schemes as well, as illustrated in Appendix .1 for the Crank-Nicolson method for the stochastic heat equation.

Previous studies Bell:07 (); FluctuatingHydro_Coveney () have measured the accuracy of numerical schemes through the variance of the fields in real space, which, by Parseval’s theorem, is related to the integral of the structure factor over all wavenumbers. For the Euler scheme (43) for the stochastic heat equation this can be calculated analytically,

showing first-order temporal accuracy (in the weak sense). For the predictor-corrector scheme (45), on the other hand, . It is important to note, however, that using the variance as a measure of accuracy of stochastic real-space integrators is both too rough and also too stringent of a test. It does not give insights into how well the equipartition is satisfied for the different modes, and, at the same time, it requires that the structure factor be good even for the highest wavenumbers, which is unreasonable to ask from a finite-stencil scheme.

For pseudo-spectral methods, as studied for the incompressible fluctuating Navier-Stokes equation in Ref. StochasticImmersedBoundary (); StochasticImmersedBoundary_Theory (), one can modify the spectrum of the stochastic forcing so as to balance the numerical stencil artifacts, and one can also use an (exact) exponential temporal integrator in Fourier space to avoid the artifacts of time stepping. However, for finite-volume schemes, a more reasonable approach is to keep the stochastic fluxes uncorrelated between disjoint cells (which is actually physical), and instead of looking at the variance, focus on the accuracy of the static structure factor for small wavenumbers. Specifically, basic schemes will typically have , while multi-step schemes will typically achieve or higher temporal order, or even .

### v.2 Dynamic Structure Factor

It is also constructive to study the full dynamic structure factor for a given numerical scheme, especially for small wavenumbers and low frequencies. This is significantly more involved in terms of analytical calculations and the results are analytically more complicated, especially for multi-stage methods and more complex equations. For the Euler scheme (43) the solution to Eq. (27) is

 Sk,ω=2χ1χ−12μk22Δt−2(1−cosΔω)+χ21χ−12μ2k4,

where and . This shows that the dynamic structure factor does not converge to the correct answer for all wavenumbers even in the limit , namely

 limβ→0Sk,ω=2χ1μk2ω