A finite element method for stochastic diffusion equations using fluctuating hydrodynamics

Abstract

We present a finite element approach for diffusion problems with thermal fluctuations based on a fluctuating hydrodynamics model. The governing equations are stochastic partial differential equations with a fluctuating forcing term. We propose a discrete formulation of the fluctuating forcing term that has the correct covariance matrix up to a standard discretization error. Furthermore, we derive a linear mapping to transform the finite element solution into an equivalent discrete solution that is free of the artificial correlations introduced by the spatial discretization. The method is validated by applying it to two diffusion problems: a second-order diffusion equation and a fourth-order diffusion equation. The theoretical (continuum) solution to the first case presents spatially decorrelated fluctuations, while the second case presents fluctuations correlated over a finite length. In both cases, the numerical solution presents a structure factor that approximates well the continuum one.

1 Introduction

With microfluidics technologies becoming increasingly widespread and multiple applications already reaching the market level volpatti2014commercialization , the great challenges of technology now shift towards the nanoscale, with potential for multiple breakthroughs and ground-breaking innovations. In analogy with microfluidics, the rapidly growing field of research that studies the flow and transport phenomena in fluidic environments of nanoscopic dimensions has been termed “nanofluidics” bocquet2020nanofluidics . The blooming of nanofluidics was enabled in the last decade by the improvements and development of new experimental and computational techniques and has shown multiple phenomena that have no counterparts at larger scales faucher2019critical ; laine2020nanotribology .

Fluid flow, solute transport and chemical reactions look very different at the nanoscale than they do at the macroscopic level. At the small scales, the thermal fluctuations of solvent and solute molecules are relevant and cannot be ignored. Even in situations where one can average over many molecules, quantities like mass, temperature and momentum fluctuate around their mean value de2006hydrodynamic . These fluctuations are very well understood around equilibrium de2006hydrodynamic but can lead to unexpected and highly nontrivial phenomena when they occur in systems driven out of equilibrium donev2011diffusive ; peraud2017fluctuation ; vutukuri2020active .

As most of the applications in energy harvesting, membrane technology and biomedicine involve systems driven out from equilibrium, understanding the effects of thermal fluctuations in these instances is of utmost importance. Thermal fluctuations are also particularly large inside eukaryotic cells, where some molecules are present in very small numbers and need to be recruited in specific locations to achieve their biological functions brown2013linking . Deterministic models fail to predict many of the important consequences of thermal fluctuations, which calls for the development of efficient numerical tools that can include them within an efficient framework.

To do so, one possibility is to modify the deterministic partial differential equations (PDEs) that are typically used to model phenomena at larger scales to include stochasticity. One such approach to incorporate the effect of the thermal fluctuations in the transport equations was proposed by Landau and Lifschitz landau2013statistical , and is based on adding a stochastic forcing term that satisfies the fluctuation-dissipation balance. The transport equations that result from adding a fluctuating forcing term are stochastic partial differential equations (SPDEs). This approach to include thermal fluctuations into deterministic PDEs has been called fluctuating hydrodynamics de2006hydrodynamic and it has been used to study stochastic drift-diffusion equations (e.g. donev2011diffusive ; leonard2013stochastic ; chaudhri2014modeling ; peraud2016low ; donev2019fluctuating ), stochastic reaction-diffusion equations (e.g. atzberger2010spatially ; bhattacharjee2015fluctuating ; kim2017stochastic ), the Brownian motion of colloidal particles (e.g. usabiaga2013inertial ; keaveny2014fluctuating ; delmotte2015simulating ; de2016finite ; sprinkle2019brownian ; westwood2022generalised ), fluctuations of viscoelastic fluids hutter2018fluctuating ; hutter2020fluctuating , stochastic thin-film equations (e.g. tsekov1993effect ; sprittles2023rogue ; liu2023thermal ), diffusion within membranes and liquid-liquid interfaces (e.g. wang2013dynamic ; rower2022surface ), and bubble nucleation and growth (e.g. gallo2018thermally ; gallo2023nanoscale ). The advantage of a fluctuating hydrodynamics approach compared to discrete particle methods is apparent when one needs to study mesoscopic systems that involve a vast number of particles and long timescales compared to the characteristic diffusion time of individual particles.

In this work, we propose a finite element framework to solve diffusion equations using fluctuating hydrodynamics. The governing equations are SPDEs of the kind

$$\frac{\partial{u}}{\partial t}-\nabla\cdot\left(D\nabla F(u)\right)=\nabla\cdot\left(\sqrt{2Du}\boldsymbol{\zeta}\right)\,.$$ (1)

where $u$ is the number density of a diffusing species, $D$ is its diffusion coefficient and $\boldsymbol{\zeta}$ represents white noise in space and time, with stochastic properties given by

$$\left<\zeta_{i}(\mathbf{x},t)\right>=0\,$$ (2)

and

$$\left<\zeta_{i}(\mathbf{x},t)\zeta_{j}(\mathbf{x}^{\prime},t^{\prime})\right>=\delta_{ij}\delta(\mathbf{x}-\mathbf{x}^{\prime})\delta(t-t^{\prime})\,.$$ (3)

In the case of $F(u)=u$, Equation (1) can be formally derived from the motion of independent point-sized Brownian walkers using Ito calculus dean1996langevin ; donev2014dynamic , and it is nothing more than a rewriting of the equations of motion of a collection of independent Brownian walkers in terms of their density, $u$. Stochastic equations of the form of Eq. (1) also arise in the context of Dynamic Density Functional Theory donev2014dynamic ; archer2004dynamical . While an ensemble-average of the stochastic equations describing the trajectories of the Brownian walkers would lead to a deterministic PDE describing the time evolution of the particle distribution, an SPDE like Equation (1) describes the evolution of a coarse-grained density distribution instead. Detailed discussions on the derivation of deterministic and stochastic equations in this context and their different interpretations can be found in archer2004dynamical ; donev2014dynamic ; chavanis2015generalized ; chavanis2019generalized .

In Equation (1) we have abused notation to highlight its formal similarity to deterministic PDEs, but, in fact, neither the solution $u$ nor the white-noise source term $\boldsymbol{\zeta}$ can be interpreted as differentiable point-wise functions. Still, as mentioned above, numerical solutions to equations based on fluctuating hydrodynamics such as Equation (1) are widely used to capture the physical effects from the effect of thermal fluctuations. Furthermore, recent advances in the theoretical analysis front seem to confirm that the truncation of the fluctuations for small wavelengths, including the natural truncation resulting from standard spatial discretization techniques, leads to well-posed fluctuating hydrodynamics equations cornalba2023dean .

When solving such equations, one is interested in the expected value of the solution, but also in its second-order statistical moments. For instance, the so-called structure factor, related to the spatial Fourier transform of the autocorrelation function of the solution, can be measured experimentally de2006hydrodynamic . However, both the choice of spatial discretization and temporal integration schemes influence the fluctuation-dissipation balance and can introduce artificial spatial correlations in the solution, thus leading to non-physical second-order statistical moments. delong2013temporal ; donev2010accuracy

For this reason, the effect of time integrators on the solution of fluctuating-hydrodynamics equations has been investigated in depth delong2013temporal ; donev2010accuracy . As for the spatial discretization, finite volume methods (e.g. donev2010accuracy ; balboa2012staggered ; russo2021finite ; kim2017stochastic ) can be designed so as not to introduce additional artificial correlations, since each degree of freedom typically corresponds to a single finite volume, for which the solution can be integrated independently of the rest of finite volumes. This is not the case, however, for finite element methods in general, which yield solutions with artificial correlations that depend on the choice of shape functions de2015finite . This makes finite element solutions difficult to interpret physically and nonlinear terms in the equations difficult to handle mathematically. To overcome this limitation for general finite element methods, it has been suggested that specialized finite element discretization schemes should be used, which are designed based on physical considerations instead of just numerical analysis de2015finite . The rationale behind this strategy is based on the idea that discretization and physical coarse-graining are deeply connected.

In the present work, we start from a different perspective and propose an approach that yields a finite element solution to Equation (1) in which both the expected value of the solution and its second-order statistical moments are well approximated. In particular, we propose to first obtain a numerical solution using any suitable spatial discretization that is selected based on numerical analysis only and then apply a linear transformation that removes from the numerical solution the artificial correlations introduced by the chosen spatial discretization. With the proposed approach, the choice of spatial discretization techniques is not restricted anymore to those that produce any particular spatial correlations, and we can use the finite element method and its powerful mathematical framework to deal with complex geometries and boundary conditions without limitations. The resulting solution has a physical interpretation that is detached from the numerical discretization, just as numerical solutions to deterministic PDEs do.

We work in a finite element framework with a standard Galerkin discretization. Nevertheless, our conclusions are generalizable to other numerical approaches that involve a spatial discretization. We apply our method to a second-order diffusion equation and to a fourth-order diffusion equation. The two test cases are complementary, since one presents only short-ranged fluctuations, while the other has a solution that is spatially correlated over a finite length. Both cases were already used in Ref. de2015finite to test a Petrov-Galerkin finite element implementation.

The outline of this paper is as follows: In Section 2 we describe the finite element discretization for a second-order fluctuating diffusion equation that has a solution with short-ranged correlations. We also propose an implementation for the fluctuating forcing term that has the correct covariance up to a normal discretization error. In Section 3, we present finite element results for a one-dimensional boundary value problem of the second-order equation, and we discuss the presence of artificial spatial correlations in the solution due to the effects of the spatial discretization and temporal integration. In Section 4, we propose a linear mapping to remove from the numerical solution the artificial correlations introduced by the spatial discretization, and we apply it to the results of the same second-order diffusion equation. In Section 5, we apply the proposed method to a one-dimensional fourth-order diffusion equation, to test its performance in the case of a solution that is spatially correlated over a finite length. Finally, we summarize our conclusions in Section 6.

2 Numerical approximation

In this section, we describe the numerical schemes used to solve Equation (1) numerically. We focus here on the case $F(u)=u$, for which numerical results will be presented in Sections 3 and 4. Nevertheless, most of the discretization details presented in this section are also relevant for the case $F(u)=u+\left(\ell_{0}/2\pi\right)^{2}\nabla^{2}u\,,$ which will be discussed in Section 5.

2.1 Weak form and finite element discretization

We consider the second-order equation

$$\frac{\partial{u}}{\partial t}-D\nabla^{2}u=\nabla\cdot\left(\sqrt{2Du}\boldsymbol{\zeta}\right)\,$$ (4)

with constant diffusion coefficient $D$. This corresponds to Dean’s equation dean1996langevin , which describes the fluctuations of the density in a system of particles undergoing independent Brownian motions. The concentration $u$ represents therefore a coarse-grained number of particles per unit volume.

A weak formulation of Equation (4) can be obtained by multiplying it by the conjugate of a test function $v$ and integrating over the domain $\Omega$. After integration by parts, one obtains

$$\int_{\Omega}{v^{*}\frac{\partial{u}}{\partial t}}\,d^{3}\mathbf{x}+\int_{\Omega}{D\nabla v^{*}\cdot\nabla u}\,d^{3}\mathbf{x}=-\int_{\Omega}{\nabla v^{*}\cdot\left(\sqrt{2Du}\boldsymbol{\zeta}\right)}\,d^{3}\mathbf{x}+\\ +\int_{\partial\Omega}{v^{*}\left(\sqrt{2Du}\boldsymbol{\zeta}\right)\cdot\mathbf{n}}\,d^{2}\mathbf{x}+\int_{\partial\Omega}{Dv^{*}\nabla u\cdot\mathbf{n}}\,d^{2}\mathbf{x}\,,$$ (5)

where $\partial\Omega$ is the boundary of the domain and $\mathbf{n}$ its normal unitary vector.

The domain $\Omega$ is discretized by dividing it into non-overlapping elements $\Omega_{e}$, such that $\Omega=\cup_{e}\Omega_{e}$ with $\Omega_{e}\cap\Omega_{e^{\prime}}=\varnothing$ for $e\neq e^{\prime}$. The solution space, as well as the test function space, are discretized using a given basis of shape functions. Here we use a standard Galerkin discretization with Lagrange polynomials as shape functions $\phi_{i}=\phi_{i}(\mathbf{x})$ for both the test function $v$ and the solution $u$, such that

$$u(\mathbf{x},t)\approx\sum_{i=1}^{N_{\text{dof}}}{\widetilde{u}_{i}(t)\phi_{i}(\mathbf{x})}\,,$$ (6)

and

$$v(\mathbf{x},t)\approx\sum_{i=1}^{N_{\text{dof}}}{\widetilde{v}_{i}(t)\phi_{i}(\mathbf{x})}\,.$$ (7)

where $N_{\text{dof}}$ is the total number of degrees of freedom.

The goal of the finite element method is to find the values $\widetilde{u}_{i}(t)$ that satisfy the discretized weak equation (5) for any values $\widetilde{v}_{i}(t)$. Using Equations (5), (6) and (7), we can express this problem as a linear system of equations in matrix form

$$\mathbf{M}\frac{\partial\widetilde{\mathbf{u}}}{\partial t}+\mathbf{D}\widetilde{\mathbf{u}}=\mathbf{f}+\mathbf{f}_{BC}\,,$$ (8)

where $\mathbf{M}$ is the so-called mass matrix

$$M_{ij}=\int_{\Omega}{\phi_{i}\phi_{j}}d^{3}\mathbf{x}\,$$ (9)

$\mathbf{D}$ is the diffusion matrix

$$D_{ij}=\int_{\Omega}{D\nabla\phi_{i}\cdot\nabla\phi_{j}}d^{3}\mathbf{x}\,,$$ (10)

$\widetilde{\mathbf{u}}$ is an array that contains the coefficients $\widetilde{u}_{i}$ for each degree of freedom, $\mathbf{f}_{BC}$ is an array that incorporates the effect of the boundary conditions in the system and $\mathbf{f}$ is the random excitation term defined such that

$$f_{i}(t)=-\int_{\Omega}{\nabla\phi_{i}\cdot\sqrt{2Du}\mathbf{\boldsymbol{\zeta}}(t)}d^{3}\mathbf{x}\,.$$ (11)

The covariance matrix of the components $f_{i}(t)$ can be expressed as

$$\left<f_{i}(t)f_{j}(t^{\prime})\right>={2D\delta(t-t^{\prime})}\int_{\Omega}{u(\mathbf{x})\nabla\phi_{i}(\mathbf{x})\cdot\nabla\phi_{j}(\mathbf{x})}d^{3}\mathbf{x}\,,$$ (12)

where $\left<.\right>$ indicates the expected value. The implementation of $\mathbf{f}$, which represents the effect of the thermal fluctuations, is discussed in Section 2.2, and the time integration schemes are presented in Section 2.3.

2.2 Finite element formulation for the thermal fluctuations

To find a discrete approximation of $\mathbf{f}$, we seek to postulate a formulation that satisfies Equation (12). Here we propose an approach inspired by previous work in the context of the Stokes equations with fluctuating hydrodynamics de2016finite . To this end, we make use of the same numerical integration schemes that are used in the finite element implementation, and that allow the integration of a given function $g(\mathbf{x})$ as

$$\int_{\Omega}g(\mathbf{x})d^{3}\mathbf{x}\approx\sum_{k=1}^{N_{g}}{w_{k}g(\mathbf{x}_{k})}\,,$$ (13)

where $N_{g}$ is the number of integration points, $g(\mathbf{x}_{k})$ is the value of the function evaluated at the integration point and $w_{k}$ is the weight corresponding to integration point $k$. If an isoparametric formulation is used, $w_{k}$ would include both the weight corresponding to the integration rule in the reference element as well as the Jacobian for the mapping from the reference element to the physical element. We assume that $w_{k}\geq 0$. With this, we can postulate the following formulation for array $\mathbf{f}$ before the time discretization, as defined by Equation (11):

$$f_{i}=-\sqrt{2D}\sum_{k=1}^{N_{g}}{\sqrt{w_{k}u(\mathbf{x}_{k},t)}\boldsymbol{\zeta}_{k}(t)\cdot\nabla\phi_{i}(\mathbf{x}_{k})}\,,$$ (14)

where the components of $\boldsymbol{\zeta}_{k}$ are stochastic Gaussian process with $\left<\zeta_{k,m}\right>=0$ and $\left<\zeta_{k,m}(t)\zeta_{l,m}(t^{\prime})\right>=\delta(t-t^{\prime})\delta_{kl}$. The formulation given by Equation (14) leads to the following covariance matrix

$$\begin{split}\left<f_{i}\right.&\left.(t)f_{j}(t^{\prime})\right>=\\ \,&={2D}\sum_{m=1}^{3}{\sum_{k=1}^{N_{g}}{\sum_{l=1}^{N_{g}}{\sqrt{w_{k}w_{l}u(\mathbf{x}_{k},t)u(\mathbf{x}_{l},t^{\prime})}\frac{\partial\phi_{i}(\mathbf{x}_{k})}{\partial x_{m}}\frac{\partial\phi_{j}(\mathbf{x}_{l})}{\partial x_{m}}\left<\zeta_{k,m}(t)\mathbf{\zeta}_{l,m}(t^{\prime})\right>}}}\\ \,&={2D\delta(t-t^{\prime})}\sum_{k=1}^{N_{g}}{{u(\mathbf{x}_{k},t)}\nabla\phi_{i}(\mathbf{x}_{k})\cdot\nabla\phi_{j}(\mathbf{x}_{k}){w_{k}}}\,.\end{split}$$ (15)

This covariance matrix is identical to that of Equation (12) up to an error introduced by the numerical integration scheme, thus proving that the postulated formulation for $\mathbf{f}$ in Equation (14) has the correct second-order statistical moments. To preserve the fluctuation-dissipation balance, the integration rule used to compute the diffusion matrix $\mathbf{D}$ and the forcing term $\mathbf{f}$ should be the same de2016finite . In the particular case that the value of the concentration $u$ is large enough to linearize the fluctuating forcing term, the co-variance of the stochastic noise term for degrees of freedom $i$ and $j$ becomes

$$\left<f_{i}(t)f_{j}(t^{\prime})\right>=2u_{0}\,\delta(t-t^{\prime}){D}_{ij}\,$$ (16)

where $D_{ij}$ is the $(i,j)$ component of the diffusion matrix $\mathbf{D}$ and $u_{0}$ represents the average concentration. In A, an alternative to Equation (14) to define the linearized $\mathbf{f}$ is described, based on a decomposition of matrix $\mathbf{D}$. This alternative implementation of the linearized forcing term has computational advantages and is equivalent to Equation (14) if the average concentration $u_{0}$ in the domain is large enough.

2.3 Time integration and fluctuation-dissipation balance

The effect of the temporal integrators for fluctuating-hydrodynamics equations has been explored extensively donev2010accuracy ; delong2013temporal . In this work, we use one-stage time integration schemes. For Equation (8), the result for each time iteration can be expressed as

$$\widetilde{\mathbf{u}}^{n+1}=\left[\mathbf{M}+(1-\alpha)\Delta t\mathbf{D}\right]^{-1}\left[\left(\mathbf{M}-\alpha\Delta t\mathbf{D}\right)\,\widetilde{\mathbf{u}}^{n}+\sqrt{\Delta t}\,\mathbf{f}^{n}+\Delta t\mathbf{f}_{BC}^{n}\right]\,,$$ (17)

where the superindex $(.)^{n}$ indicates that the term is evaluated at the time step $n$, and $\alpha$ is a constant that indicates whether the time integration is explicit ($\alpha=1$), implicit ($\alpha=0$) or semi-implicit ($\alpha=1/2$). Based on Equation (14) in time, the stochastic forcing term in Equation (17) can be expressed as

$$f_{i}^{n}=-\sqrt{2D}\sum_{k=1}^{N_{g}}{\sqrt{w_{k}u(\mathbf{x}_{k},t_{n})}\boldsymbol{z}^{n}_{k}\cdot\nabla\phi_{i}(\mathbf{x}_{k})}\,,$$ (18)

where the components of $\boldsymbol{z}_{k}$ are stochastic Gaussian process with $\left<z_{k,m}\right>=0$ and $\left<z_{k,m}(t)z_{l,m}(t^{\prime})\right>=\delta_{kl}$. We note that the solution, $u(\mathbf{x}_{k},t_{n})$, at the integration point, $\mathbf{x}_{k}$, has to be evaluated at time $t_{n}$ because the stochastic term has been derived using an Ito interpretation dean1996langevin ; donev2014dynamic ; kloeden1992stochastic .

If the fluctuating term can be linearized, the alternative expression provided by Equation (79) can be used. The elements of the covariance matrix of the discrete forcing term given by Equation (18) can be expressed as

$$\left<f_{i}^{n}f_{j}^{m}\right>={2D\delta_{nm}}\sum_{k=1}^{N_{g}}{{u(\mathbf{x}_{k},t)}\nabla\phi_{i}(\mathbf{x}_{k})\cdot\nabla\phi_{j}(\mathbf{x}_{k}){w_{k}}}\,.$$ (19)

In Equation (17), the value $\alpha=1/2$ corresponds to the implicit midpoint method (Crank-Nicholson scheme). This scheme is not only stable regardless of the value of the time step, but it also produces a solution with static spatial correlations that are independent of the time step as well donev2010accuracy ; delong2013temporal . Indeed, if we consider for now an infinite domain, so that we can neglect the effect of the boundary conditions, and we take the covariance of Equation (17), we obtain

$$\mathbf{D}\,\mathbf{C}\,\mathbf{M}^{T}+\mathbf{M}\,\mathbf{C}\,\mathbf{D}^{T}+\Delta t(1-2\alpha)\mathbf{D}\,\mathbf{C}\,\mathbf{D}^{T}=\mathbf{C}_{ff}\,,$$ (20)

where $\mathbf{C}_{ff}$ is the covariance matrix of the forcing term $\mathbf{f}$, given by Equation (19), and $\mathbf{C}$ is the covariance matrix of the discrete solution

$$\mathbf{C}=\left<\left(\mathbf{\widetilde{u}}-\left<\mathbf{\widetilde{u}}\right>\right)\left(\mathbf{\widetilde{u}}^{T}-\left<\mathbf{\widetilde{u}}^{T}\right>\right)\right>\,.$$ (21)

Equation (20) expresses the fluctuation-dissipation balance at the discrete level. The term that depends on the time step vanishes if $\alpha=1/2$. By contrast, any other value of $\alpha$, including $\alpha=1$ (fully explicit) and $\alpha=0$ (fully implicit), will not make this term vanish and, as a result, will lead to a solution with spatial correlations that depend on the ratio

$$\beta=\frac{D\Delta t}{\Delta x^{2}}\,,$$ (22)

where $\Delta t$ and $\Delta x$ are the time step and the reference cell size, respectively.

The choice of time integration scheme influences, therefore, the fluctuation-dissipation balance, and, therefore, the spatial correlations of the solution. In this work, we use the Crank-Nicholson scheme ($\alpha=1/2$) unless otherwise stated. With this scheme, Equation (20) becomes

$$\mathbf{D}\,\mathbf{C}\,\mathbf{M}^{T}+\mathbf{M}\,\mathbf{C}\,\mathbf{D}^{T}=\mathbf{C}_{ff}\,,$$ (23)

which shows that the fluctuation-dissipation balance is independent of the time step. In the linearized case, Equation (23) becomes

$$\mathbf{D}\,\mathbf{C}\,\mathbf{M}^{T}+\mathbf{M}\,\mathbf{C}\,\mathbf{D}^{T}=2u_{0}\mathbf{D}\,.$$ (24)

It is worth noting that in general, the fluctuation-dissipation balance is modified by the presence of boundary conditions, which we are not considering in this simplified analysis.

We see from Equation (24) that, even with a choice of a temporal integrator that does not introduce additional correlations in the solution, the covariance matrix of the solution will depend on the mass matrix $\mathbf{M}$. As a result, the discrete solution may display spatial correlations that are not physical but induced by the spatial discretization. Indeed, if the diffusion and mass matrices are symmetric, the covariance matrix

$$\mathbf{C}=u_{0}\mathbf{M}^{-1}$$ (25)

satisfies Equation (24). If $\mathbf{M}$ is diagonal, $\mathbf{C}$ will also be diagonal, and the solution for each degree of freedom will be uncorrelated with the solution for the rest of degrees of freedom; conversely, if $\mathbf{M}$ is not diagonal, as is in general the case with finite-element methods because of the way the shape functions are defined in the elements, the solution for the different degrees of freedom will be correlated. It is worth noting that these artificial correlations do not represent a numerical error, but are instead a natural consequence of the choice of shape functions. Still, even if these artificial correlations do not represent a numerical error, they make the solution difficult to interpret and to compare to experimental measurements, and they can make non-linear terms of fluctuating-hydrodynamics equations difficult to treat. In Section 3, we show finite element results for a (bounded) 1D problem with periodic boundary conditions that present discretization-induced spatial correlations in the solution. Furthermore, in Section 4, we propose a change of basis to transform the discrete solution into an equivalent discrete solution that is free of these artificial correlations.

2.4 Coarse-graining and discretization

Stochastic diffusion equations such as the ones that we are considering in this work can be seen as coarse-grained descriptions of the underlying system of particles. In fact, the macroscopic diffusion equations can be obtained from a description of the microscopic system through a coarse-graining procedure de2011coarse ; espanol2015coupling . In the case of diffusion equations such as Equation (1), the concentration $u$ is a statistical representation of the distribution of the particles moving in the domain. In particular, $u$ satisfies

$$\int_{\Omega}{\gamma_{i}(\mathbf{x})u(\mathbf{x},t)}d^{3}\mathbf{x}=\sum_{p=1}^{N_{p}}{\gamma_{i}(\mathbf{x}_{p}(t))}\,,$$ (26)

where $N_{p}$ is the number of particles of the system, $\mathbf{x}_{p}$ their positions and $\gamma_{i}$ is a suitable coarse-graining function associated with the discrete coarse-graining volume $i$.

Using the approximation given by Equation (6), we obtain

$$\int_{\Omega}{\gamma_{i}(\mathbf{x})\sum_{j}\phi_{j}(\mathbf{x})\widetilde{u}_{j}(t)}d^{3}\mathbf{x}=\sum_{p=1}^{N_{p}}{\gamma_{i}(\mathbf{x}_{p}(t))}=\widetilde{n}_{i}\,.$$ (27)

In this work, we want to stress that, although the spatial correlations of the coarse-grained variable $\widetilde{\mathbf{n}}$ depend on the functions $\gamma_{i}$, the spatial correlations of the discrete concentration $\widetilde{\mathbf{u}}$ depend on the choice of shape functions $\phi_{i}$ used to approximate $u$ as expressed by Equation (6), and not on the coarse-graining functions $\gamma_{i}$. We will discuss this further in Section 4.1, where we propose an approach to find a discrete approximation of the continuum concentration $u$ that is free of discretization-induced spatial correlations.

3 Finite element results for a 1D diffusion problem

In this section, we present numerical results for a one-dimensional version of Equation (4), using the discretization schemes described in Section 2, and illustrate the effects of the spatiotemporal discretization on the spatial correlations of the discrete solution.

3.1 One-dimensional boundary value problem

We consider a one-dimensional diffusion problem with thermal fluctuations, governed by equation

$$\frac{\partial u({x},t)}{\partial t}-D\frac{\partial^{2}u({x},t)}{\partial x^{2}}=\frac{\partial}{\partial x}\left(\sqrt{2Du}\mathbf{{\zeta}}({x},t)\right)\,,$$ (28)

where $D$ is the diffusion coefficient, which we assume to be constant, and $\zeta$ is a Gaussian white noise with

$$\left<\zeta(x,t)\right>=0\,$$ (29)

and

$$\left<\zeta(x,t)\zeta(x^{\prime},t^{\prime})\right>=\delta(x-x^{\prime})\delta(t-t^{\prime})\,.$$ (30)

The results that will be presented here correspond to dimensionless variables, with a domain length $L=1$ and reference time $T=L^{2}/D=1$. Periodic boundary conditions

$$u(x=L,t)=u(x=0,t)\,,$$ (31)

are prescribed at the ends of the computational domain, and the initial conditions are defined as

$$u(x,t=0)=u_{0}\,.$$ (32)

The periodic boundary conditions guarantee the conservation of mass (or, equivalently, of the number of particles) in the domain, so that the average concentration in the domain is always $u_{0}$.

If we linearize the equation around a large value of the average concentration $u_{0}$, the equal-time autocorrelation function of the fluctuations at points $x$ and $x^{\prime}$ for an infinite domain can be expressed as de2015finite

$$\left<u^{\prime}(x,t)u^{\prime}(x^{\prime},t)\right>=u_{0}\delta(x-x^{\prime})\,,$$ (33)

where $u^{\prime}=u-\left<u\right>=u-u_{0}$, and the so-called static structure factor, which is the Fourier transform of the equal-time correlation function de2006hydrodynamic , becomes

$$S(k)=\left<\mathcal{U}(k,t)\mathcal{U}(-k,t)\right>=u_{0}$$ (34)

where $k$ is the wavenumber, and where the Fourier transform $\mathcal{U}(k)$ of $u(x)$ is defined as

$$\mathcal{U}(k)=\int_{-\infty}^{\infty}{u(x)\exp{(-ikx)}}dx\,.$$ (35)

3.2 Discrete model, simulation setup and discrete structure factor

We apply the numerical approximation described in Section 2 to Equation (28). The one-dimensional domain of length $L$ is split into $N_{e}$ identical elements. Lagrange polynomials are used as shape functions, both linear ($p1$-elements) and quadratic ($p2$-elements). When $p1$-elements are used, the total number of nodes is $N_{n}=N_{e}+1$; when $p2$-elements are used, the number of nodes is $N_{n}=2N_{e}+1$. Since the nodes are always equidistant, the distance between two nodes regardless of the type of element is $\Delta x=L/(N_{n}-1)$. The number of degrees of freedom is $N_{\text{dof}}\equiv N_{n}$.

We use $\widetilde{u}^{0}_{i}=u_{0}$ for every $i=1,..,N_{\text{dof}}$ as the initial condition, where the value of the reference concentration $u_{0}$ is selected to be high enough for the behavior of the problem to be approximately linear. We then simulate an initial number of time steps $L^{2}/(D\Delta t)$, so that the solution reaches an equilibrium state before starting to collect results, and, finally, we run the simulations for a number $N_{t}$ of time steps, for which data are collected.

The periodic boundary conditions are implemented by adding Lagrange multipliers to the discrete finite element system (8), to force the solution at the extremes of the domain to be identical. These boundary conditions ensure that the total number of particles in the domain remains constant (up to a certain numerical error). This is verified in the computations by integrating the results for the concentration over the domain. For the computations presented in this work, the standard deviation of the total computed number of particles in the domain varies with the discretization parameters, but it is typically smaller or of the order of $\sim$110-6$\,u_{0}L$.

As a consequence of the mass conservation in the domain, the structure factor for wavenumber $k=0$ is zero, unlike the structure factor of the case with an infinite domain given by Equation (34). For the rest of discrete wavenumbers $k_{m}=2\pi(m-1)/L$, with $m>1$, and provided that the discretization does not introduce artificial correlations in the solution, the structure factor of the solution should approximate that of Equation (34), $S=u_{0}$. This is equivalent to saying that, although the solution for the continuum equation in an infinite domain has a Gaussian probability distribution, the theoretical discrete solution in our bounded one-dimensional domain follows a multinomial distribution corresponding to $N_{p}=u_{0}L$ particles in $N_{n}-1$ bins. It is worth noting that, as the discretization is refined (i.e. as $N_{n}$ increases), the multinomial distribution becomes closer to a Gaussian distribution.

However, the discretization introduces artificial correlations in the solution, so the structure factor of the discrete FE results will not necessarily converge to the theoretical one. An analysis of the one-dimensional discretized equation with uniform linear ($p1$) Lagrange elements shows instead that the numerical results with $p1$-elements will converge to a solution with a discrete structure factor that can be expressed as de2015finite

$$S_{\text{th, FE-p1}}(k_{m})=\frac{9u_{0}}{\left[2+\cos{(k_{m}\Delta x)}\right]^{2}}\sum_{j\in\mathbb{Z}}{\left(\frac{\sin{\left(\frac{k_{m}\Delta x}{2}-\pi j\right)}}{\frac{k\Delta x}{2}-\pi j}\right)^{4}}\,,$$ (36)

where $\mathbb{Z}$ represents the set of integer numbers. The structure factor given by Equation(36) only converges to the continuum structure factor given by Equation (34) in the limit $k_{m}\Delta x\rightarrow 0$.

To evaluate our numerical approach, we compare the theoretical structure factor with the structure factor of the numerical solution, computed as

$$S(k_{m})=\left<{\widetilde{\mathcal{U}}}(k_{m}){{\widetilde{\mathcal{U}}}^{*}}(k_{m})\right>\,,$$ (37)

where $\widetilde{\mathcal{U}}(k_{m})$ is the discrete Fourier transform of the normalized discrete concentration, defined as

$${\widetilde{\mathcal{U}}}(k_{m})=\frac{1}{\sqrt{L}}\sum_{n=1}^{N_{n}-1}{\widetilde{u}}_{n}\Delta\widetilde{L}_{n}\exp\left(-\frac{i2\pi(m-1)(n-1)}{N_{n}-1}\right)\,,$$ (38)

${\widetilde{\mathcal{U}}}^{*}(k_{m})$ is its complex conjugate and $\Delta\widetilde{L}_{n}$ is the equivalent volume (equivalent length in 1D) associated to node $n$,

$$\Delta\widetilde{L}_{n}=\int_{\Omega}{\phi_{n}(x)}dx\,.$$ (39)

3.3 Results with a semi-implicit time integrator

Figure 1 shows the discrete structure factor based on the finite element results using linear Lagrange elements and a semi-implicit integration scheme. As explained in Section 2.3, a semi-implicit integration scheme has the advantage of making the spatial correlations independent from the time step. As shown in the figure, the discrete structure factor converges towards the theoretical structure factor of the discrete weak equation, given by Equation (36). The error

$$e_{\text{FE}}=\frac{1}{N_{k}}\sum_{m=1}^{N_{k}}{\frac{\left|S(k_{m})-S_{\text{FE-p1}}(k_{m})\right|}{\left|S_{\text{FE-p1}}(k_{m})\right|}}\,,$$ (40)

where $k_{m}=2\pi m/L$ (for $m=1,2,..,(N_{n}-1)/2$) are the $N_{k}=(N_{n}-1)/2$ discrete wavenumbers, is reduced with increasing number of elements $N_{e}$ following a power law ${e_{FE}}\propto N_{e}^{-1}$.

Figure 1: Left: Structure factor of the discrete solution; theoretical structure factor for the discrete equation (solid line) and computed with $N_{n}=51$ (dashed) and $N_{n}=101$ (dotted). Right: error ${e_{FE}}$ in the computed discrete structure factor vs the number of elements ($\circ$). All computations have $p1$-elements, $\Delta t=$110-4$$, $N_{t}=$1106$$, $u_{0}=10000$. It can be seen that the solution obtained from the FE computation converges to the theoretical solution of the discretized equation as the number of elements increases.

However, as explained in Section 3.2, this theoretical discrete structure factor towards which the finite element-solution converges is different from the continuum one, which is a constant ($S/u_{0}=1$). This is due to the fact that the spatial discretization introduces artificial correlations in the discrete solution; as a result, the structure factor only converges to that of the continuum in the low wavenumber limit. In Figure 2, the same results are represented, but as a function of $kL$ instead of $k\Delta x$. We see that, as the discretization is refined, the structure factor is close to the continuum one for a larger interval of wavenumbers. Nevertheless, artificial correlations are always present in the solution, and the resulting variance cannot be interpreted independently of the discretization.

Figure 2: Structure factor of the discrete solution as a function of $kL$; computed with $N_{n}=51$ (dashed), $N_{n}=101$ (dotted) and $N_{n}=201$ (dash-dotted). All computations with $p1$-elements, $\Delta t=$110-4$$, $N_{t}=$1106$$, $u_{0}=10000$.

It is worth noting that the obtained results agree with the finite element solution obtained by de la Torre and co-workers de2015finite , who use a Petrov-Galerkin discretization and also obtain a discrete structure factor that converges to that given by Equation (36). The reason for this is that the artificial correlations only depend on the shape functions used to approximate the solution $u$ in Equation (6), which are linear Lagrange polynomials in both works.

Figure 3 shows results obtained with both linear and quadratic elements. It can be seen that increasing the order contributes to approximating better the structure factor towards the continuum one. Although the convergence is still limited to the low wavenumber limit, the results obtained with quadratic elements present a larger interval of wavenumbers for which the error in the structure factor is below a given threshold.

Figure 3: Structure factor of the discrete solution computed with $p1$-elements (solid line) and $p2$-elements (dash-dot). The computations have been run with $N_{n}=101$, $\Delta t=$110-4$$, $N_{t}=$1106$$, $u_{0}=10000$.

Figure 4 shows the results obtained with quadratic elements and different mesh sizes. Like the results obtained with linear elements, the results computed with quadratic elements do not converge towards the continuum as the mesh is refined, indicating the presence of artificial correlations induced by the discretization in the solution.

Figure 4: Structure factor of the discrete solution obtained with $p2$-elements, computed with $N_{n}=51$ (dashed) and $N_{n}=101$ (dotted). The computations have been run with $\Delta t=$110-4$$, $N_{t}=$1106$$, $u_{0}=10000$.

The results presented in this section have been focused on the static structure factor because it serves the purpose of showing the presence of artificial spatial correlations (or eventual lack thereof) in the solution, which is the focus of this work. Nevertheless, in B, complementary results are presented showing the dynamic structure factor, which contains information on the time evolution of the solution as well.

3.4 Results with fully implicit and fully explicit time integrators

The results shown so far for the semi-implicit time integrator present artificial correlations that are related to the spatial discretization. However, as explained in Section 2.3, if instead of the semi-implicit integration scheme, one uses a fully explicit or a fully implicit scheme, additional spatial correlations dependent on the dimensionless time step $\beta=D\Delta t/\Delta x^{2}$ appear in the solution donev2010accuracy . Figures 5 and 6 show the results obtained for linear and quadratic elements, using fully explicit and fully implicit integration schemes. The results confirm that indeed the obtained structure factor is influenced by additional correlations that depend on the dimensionless time step $\beta$. It is worth noting that the influence of the time step on the stability is an independent issue altogether: the implicit method is unconditionally stable, while the explicit becomes unstable if the time step exceeds a certain threshold.

The results obtained with the fully explicit and fully implicit integration schemes are consistent with what other authors have shown before donev2010accuracy ; delong2013temporal , and illustrate the advantage of the semi-implicit integration scheme, as it does not require decreasing the time step to yield the correct spatial structure. In the rest of this work, only results obtained with the semi-implicit scheme ($\alpha=1/2$) are presented.

Figure 5: Structure factor of the discrete solution with $p1$-elements, computed with $N_{n}=51$ (dashed) and $N_{n}=101$ (dotted). Left: fully implicit time integration ($\alpha=0$). Right: fully explicit time integration ($\alpha=1$). All computations with $\Delta t=$110-5$$, $N_{t}=$1106$$, $u_{0}=10000$. As a reference, numerical results obtained with $\alpha=1/2$ are also displayed (solid line), which are equivalent to the case with $\beta=0$.

Figure 6: Structure factor of the discrete solution with $p2$-elements, computed with $N_{n}=51$ (dashed) and $N_{n}=101$ (dotted). Left: fully implicit time integration ($\alpha=0$). Right: fully explicit time integration ($\alpha=1$). All computations with $\Delta t=$110-5$$, $N_{t}=$1106$$, $u_{0}=10000$. As a reference, numerical results obtained with $\alpha=1/2$ are also displayed (solid line), which are equivalent to the case with $\beta=0$.

4 A linear mapping to remove artificial correlations

As explained in Sections 2.3 and 2.4, the spatial discretization can introduce artificial correlations in the numerical solution. These artificial correlations do not represent a numerical error; they are instead the natural consequence of the choice of shape functions for the spatial discretization, and they are thus part of a consistent solution in the context of that basis. In spite of this, these discretization-related correlations make the discrete solution difficult to interprete, and can also make non-linear terms in fluctuating-hydrodynamics equations difficult to handle. In this section, we propose a linear transformation between the numerical solution and an equivalent discrete solution free of these artificial correlations. The goal is to perform a change of basis to obtain an equivalent discrete solution $\widetilde{\widetilde{\mathbf{u}}}$ that approximates the solution as

$$u(\mathbf{x},t)\approx\sum_{i=1}^{N_{\text{dof}}}{\widetilde{\widetilde{u}}_{i}(t)\psi_{i}(\mathbf{x})}\,,$$ (41)

where $\psi_{i}$ are basis functions defined so as not to introduce artificial spatial correlations in the finite element solution, that is, so as to satisfy

$$\int_{\Omega}\psi_{i}(\mathbf{x})\psi_{j}(\mathbf{x})d^{3}\mathbf{x}=0\,,\,\,\,\,\text{if }i\neq j\,.$$ (42)

4.1 Spatial decorrelation matrix

In this section, we show that a linear relationship exists between the discrete solution $\widetilde{\mathbf{u}}$ and an equivalent discrete solution $\widetilde{\widetilde{\mathbf{u}}}$ that is free of the artificial correlations introduced by the spatial discretization. We start by re-writing here Equation (6)

$$u(\mathbf{x},t)\approx\sum_{i=1}^{N_{\text{dof}}}{\widetilde{u}_{i}(t)\phi_{i}(\mathbf{x})}\,,$$

where $\phi_{i}(\mathbf{x})$ are the basis functions used to approximate the solution, and Equation (41), which provides an alternative approximation using shape functions $\psi_{i}$,

$$u(\mathbf{x},t)\approx\sum_{i=1}^{N_{\text{dof}}}{\widetilde{\widetilde{u}}_{i}(t)\psi_{i}(\mathbf{x})}\,.$$

In order to satisfy Equation (42) the shape functions $\psi_{i}(\mathbf{x})$ are defined such that matrix $\mathbf{M}_{\psi\psi}$

$$({M}_{\psi\psi})_{ij}=\int_{\Omega}{\psi_{i}\psi_{j}}d^{3}\mathbf{x}\,$$ (43)

is diagonal, with

$$({M}_{\psi\psi})_{ii}=\int_{\Omega}{\psi_{i}\psi_{i}}d^{3}\mathbf{x}=\Delta\widetilde{\widetilde{V}}_{i}\,,$$ (44)

where $\Delta\widetilde{\widetilde{V}}_{i}$ is an equivalent discrete volume associated with each degree of freedom $i$,

$$\Delta\widetilde{\widetilde{V}}_{i}=\int_{\Omega}{\psi_{i}}d^{3}\mathbf{x}\,.$$ (45)

Note that we have not defined the shape functions $\psi_{i}$ beyond saying that they lead to a diagonal mass matrix, so $\Delta\widetilde{\widetilde{{V}_{i}}}$ is undefined so far. We will discuss later in this section which definition of $\Delta\widetilde{\widetilde{{V}_{i}}}$ is most advantageous (see Equation (64)).

Since the discrete solutions ${\widetilde{\mathbf{u}}}$ and $\widetilde{\widetilde{\mathbf{u}}}$ both approximate the same continuum function, we require them to satisfy

$$\sum_{j=1}^{N_{\text{dof}}}{{\widetilde{u}}_{j}(t){{\phi}}_{j}(\mathbf{x})}=\sum_{j=1}^{N_{\text{dof}}}{\widetilde{\widetilde{u}}_{j}(t)\psi_{j}(\mathbf{x})}\,.$$ (46)

Multiplying by $\phi_{i}$ and integrating over the domain leads to

$$\int_{\Omega}{{{\phi}}_{i}\sum_{j=1}^{N_{\text{dof}}}{{\phi}}_{j}}\widetilde{u}_{j}d^{3}\mathbf{x}=\int_{\Omega}{{{\phi}}_{i}\sum_{j=1}^{N_{\text{dof}}}\psi_{j}}\widetilde{\widetilde{u}}_{j}d^{3}\mathbf{x}\,,$$ (47)

which is equivalent to the following equation in matrix form

$$\mathbf{A}_{\phi}^{T}\widetilde{\mathbf{u}}=\mathbf{A}_{\psi}\,\widetilde{\widetilde{\mathbf{u}}}\,,$$ (48)

where $\mathbf{A}_{\psi}$ is a diagonal matrix satisfying

$$\mathbf{M}_{\psi\psi}=\mathbf{A}_{\psi}\,\mathbf{A}_{\psi}\,,$$ (49)

and $\mathbf{A}_{\phi}^{T}$ is the transpose of $\mathbf{A}_{\phi}$, satisfying

$$\mathbf{M}_{\phi\phi}=\mathbf{A}_{\phi}\,\mathbf{A}_{\phi}^{T}\,,$$ (50)

with matrix $\mathbf{M}_{\phi\phi}$ defined as

$$({M}_{\phi\phi})_{ij}=\int_{\Omega}{\phi_{i}\phi_{j}}d^{3}\mathbf{x}\,.$$ (51)

It is worth pointing out that, since in Section 2 we have used a standard Galerkin discretization with the same shape functions for test function and solution, $\mathbf{M}_{\phi\phi}$ and the mass matrix $\mathbf{M}$ of the finite element system are the same in our numerical setup. However, this is not necessarily the case for a general discretization, in which different shape functions may be used to approximate the test function. Therefore, $\mathbf{M}_{\phi\phi}$ is not necessarily equal to the mass matrix of the finite element system leading to the solution $\widetilde{\mathbf{u}}$.

Matrix $\mathbf{A}_{\phi}$, resulting from the decomposition given by Equation (50), can be expressed as

$$\mathbf{A}_{\phi}=\mathbf{U}_{\phi}\sqrt{\mathbf{\Sigma}_{\phi}}\mathbf{U}_{\phi}^{T}\mathbf{U}_{*}^{T}\,,$$ (52)

where $\mathbf{U}_{*}$ and $\mathbf{U}_{\phi}$ are unitary matrices and $\mathbf{\Sigma}_{\phi}$ is a diagonal matrix satisfying

$$\mathbf{M}_{\phi\phi}=\mathbf{U}_{\phi}\mathbf{\Sigma}_{\phi}\mathbf{U}_{\phi}^{T}\,.$$ (53)

Equation (50) holds for any unitary matrix $\mathbf{U}_{*}$; the decomposition (50) is therefore not unique. However, not every unitary matrix $\mathbf{U}_{*}$ will in addition guarantee that the mass is conserved. Mass conservation requires

$$\sum_{j=1}^{N_{\text{dof}}}\int_{\Omega}{{{\phi}}_{j}}\widetilde{u}_{j}d^{3}\mathbf{x}=\sum_{j=1}^{N_{\text{dof}}}\int_{\Omega}{\psi_{j}}\widetilde{\widetilde{u}}_{j}d^{3}\mathbf{x}\,,$$ (54)

which in matrix form can be written as

$$\Delta{\widetilde{\mathbf{V}}}^{T}\,\widetilde{\mathbf{u}}=\Delta\widetilde{\widetilde{\mathbf{V}}}^{T}\,\widetilde{\widetilde{\mathbf{u}}}\,,$$ (55)

where $\Delta\widetilde{\widetilde{\mathbf{V}}}$ and $\Delta{\widetilde{\mathbf{V}}}$ are vectors with the corresponding equivalent volume associated with each degree of freedom, with $\Delta\widetilde{\widetilde{\mathbf{V}}}$ defined by Equation (45) and $\Delta{\widetilde{\mathbf{V}}}$ defined by

$$\Delta\widetilde{V}_{i}=\int_{\Omega}{{{\phi}}_{i}}d^{3}\mathbf{x}\,.$$ (56)

If we now use Equations (48) and (52) inside Equation (55), we obtain

$$\Delta{\widetilde{\mathbf{V}}}^{T}\,\widetilde{\mathbf{u}}=\Delta\widetilde{\widetilde{\mathbf{V}}}^{T}\,\mathbf{A}_{\psi}^{-1}\,\mathbf{U}_{*}\,\mathbf{U}_{\phi}\,\sqrt{\mathbf{\Sigma}_{\phi}}\,\mathbf{U}_{\phi}^{T}\,\widetilde{\mathbf{u}}\,.$$ (57)

Given that this identity must be satisfied for any vector $\widetilde{\mathbf{u}}$, we obtain

$$\mathbf{U}_{*}\,\widetilde{\mathbf{y}}=\widetilde{\widetilde{\mathbf{y}}}\,,$$ (58)

where ${\widetilde{\mathbf{y}}}$ and $\widetilde{\widetilde{\mathbf{y}}}$ are vectors defined as

$${\widetilde{\mathbf{y}}}=\mathbf{U}_{\phi}\,\sqrt{\mathbf{\Sigma}^{-1}_{\phi}}\,\mathbf{U}_{\phi}^{T}\,\Delta\widetilde{\mathbf{V}}\,,$$ (59)

and

$$\widetilde{\widetilde{\mathbf{y}}}=\mathbf{A}_{\psi}^{-1}\,\Delta\widetilde{\widetilde{\mathbf{V}}}\,.$$ (60)

Therefore, matrix $\mathbf{U}_{*}$ can be defined as the rotation matrix that transforms vector ${\widetilde{\mathbf{y}}}$ into $\widetilde{\widetilde{\mathbf{y}}}$.

With these definitions, we can rewrite Equation (48) as

$$\widetilde{\widetilde{\mathbf{u}}}=\mathbf{A}_{\psi}^{-1}\,\mathbf{A}_{\phi}^{T}\,\widetilde{\mathbf{u}}=\mathbf{A}_{\psi}^{-1}\,\mathbf{U}_{*}\,\mathbf{U}_{\phi}\,\sqrt{\mathbf{\Sigma}_{\phi}}\,\mathbf{U}_{\phi}^{T}\,\widetilde{\mathbf{u}}\,,$$ (61)

which provides a mass-preserving linear relationship between the discrete solution $\widetilde{\mathbf{u}}$ and a discrete solution $\widetilde{\widetilde{\mathbf{u}}}$ that does not have artificial spatial correlations introduced by the spatial discretization:

$$\widetilde{\widetilde{\mathbf{u}}}=\mathbf{Q}\,\widetilde{\mathbf{u}}\,,$$ (62)

where matrix $\mathbf{Q}$ is given by

$$\mathbf{Q}=\mathbf{A}_{\psi}^{-1}\,\mathbf{A}_{\phi}^{T}=\mathbf{A}_{\psi}^{-1}\,\mathbf{U}_{*}\,\mathbf{U}_{\phi}\,\sqrt{\mathbf{\Sigma}_{\phi}}\,\mathbf{U}_{\phi}^{T}\,$$ (63)

and can be interpreted as a decorrelation matrix that removes the artificial spatial correlations from the discrete solution. This decorrelation matrix can be applied as a postprocessing step to remove the artificial spatial correlations from the numerical solution.

As mentioned above, $\Delta\widetilde{\widetilde{\mathbf{V}}}$ is undefined so far. A convenient definition is given by

$$\Delta\widetilde{\widetilde{\mathbf{V}}}_{i}\equiv\widetilde{\mathbf{y}}_{i}^{2}\,.$$ (64)

This definition ensures that $\widetilde{\mathbf{y}}$ and $\widetilde{\widetilde{\mathbf{y}}}$ are the same, and, therefore, $\mathbf{U}_{*}$ becomes the identity matrix. As a result, $\mathbf{A}_{\phi}$ becomes symmetric, and the decorrelation matrix becomes

$$\mathbf{Q}=\mathbf{A}_{\psi}^{-1}\,\mathbf{A}_{\phi}^{T}=\mathbf{A}_{\psi}^{-1}\,\mathbf{U}_{\phi}\,\sqrt{\mathbf{\Sigma}_{\phi}}\,\mathbf{U}_{\phi}^{T}\,.$$ (65)

An advantage of this definition is that the decorrelation matrix can be approximated by a sparse matrix with a higher sparsity degree than what is obtained when the rotation matrix $\mathbf{U}_{*}$ is a dense matrix. In Section 4.2, we discuss more details about the implementation and computational aspects of the mapping given by Equation (65).

So far, we have not considered any effect of the boundary conditions on the decorrelation matrix, and, in general, the analysis presented in this section is valid for any boundary condition. However, in some cases, it may be desirable to use a decorrelation matrix based on a matrix $\mathbf{M}_{\phi\phi}$ and vector $\Delta\widetilde{\mathbf{V}}$ that have been modified to incorporate the effect of the boundary conditions. For instance, since here we are solving problems with periodic boundary conditions, in order to maintain the periodicity of the discrete solution, only $N_{\text{dof}}-1$ elements of $\widetilde{\mathbf{u}}$ are mapped, thereby excluding the solution of one of the extremes of the domain. Matrix $\mathbf{M}_{\phi\phi}$ and vector $\Delta\widetilde{\mathbf{V}}$ become, respectively, an $(N_{\text{dof}}-1)\times(N_{\text{dof}}-1)$ matrix and a vector of length $(N_{\text{dof}}-1)$, and the elements in each eliminated row or column are added in suitable positions of the row or column corresponding to the other extreme of the domain, so as to reflect the periodicity of the domain in the arrays.

It is worth noting that the linear mapping given by Equation (62) depends only on the shape functions used to approximate the solution. This means that the choice of finite-element discretization does not need to be based on whether it produces the correct structure factor, because, once a discrete solution is obtained, an equivalent solution with a correct structure factor can be found through the proposed linear transformation. Moreover, this makes explicit that the coarse-graining and the spatial discretization used to solve the diffusion equation numerically are two separate things, and that the spatial discretization should therefore be designed on the basis of numerical analysis exclusively.

In Section 4.3, we present results for the one-dimensional boundary value problem discussed in Section 3, which are obtained by applying the linear mapping presented in this section to remove the artificial correlations introduced by the spatial discretization.

4.2 Implementation and computational aspects

The decorrelation matrix proposed in Section 4.1 is in general a dense matrix. Although this is not an issue for the 1D problems solved in this work, the mapping can lead to a high computational cost for models with a large number of degrees of freedom. To preserve the computational advantages of the finite element method, the matrix should be sparse.

Figure 7: Structure factor of the discrete solution after mapping with the decorrelation matrix; computed with a dense decorrelation matrix (dotted) and with a sparse decorrelated matrix with threshold $\epsilon=$110-5$$ ($\times$). Left: $p1$-elements; right: $p2$-elements. All computations with semi-implicit time integration, $N_{\text{dof}}=101$, $\Delta t=$110-4$$, $N_{t}=$1106$$, $u_{0}=10000$.

If we use the definition of the decorrelation matrix $\mathbf{Q}$ given by Equation (65), with the equivalent volumes $\Delta\widetilde{\widetilde{V}}_{i}$ defined by Equation (64), the matrix can be approximated by a sparse matrix, by setting to zero the elements of the matrix that have an absolute value below a threshold $\epsilon$. This is a reasonable approximation in our case, because the matrix $\mathbf{Q}$ that we obtain already has many elements with small value, corresponding to degrees of freedom belonging to two different elements that are not contiguous and that are relatively far from each other. It is worth pointing out that the sparse matrix that we obtain in this way is still a full-rank matrix.

Figure 7 shows the structure factor obtained with such truncated matrices with a threshold $\epsilon=$110-5$$. The truncated matrix gives results that are very close to those of the original dense matrix (with a maximum added relative error in the computed structure factor around $510-5$, which is significantly smaller than the typical numerical error in the solution).

Figure 8: Number $N_{\text{nz}}$ of non-zero elements as a function of the number of degrees of freedom $N_{\text{dof}}$ for mass matrix $\mathbf{M}$ with $p1$-elements ($\triangle$) and $p2$-elements ($\blacktriangle$), truncated decorrelation matrix $\mathbf{Q}$ ($\epsilon=$110-5$$) with $p1$-elements ($\circ$) and $p2$-elements ($\bullet$).

Figure 8 shows the number of non-zero elements in the decorrelation matrix with both $p1$- and $p2$-elements, computed as well with threshold $\epsilon=$110-5$$. Although the number of non-zero elements in the decorrelation matrices is higher than in the respective mass matrices, they both increase linearly with the number of degrees-of-freedom $N_{\text{dof}}$. This indicates that the computational complexity of the problem is not increased by the mapping, thus preserving the efficiency of the computational method.

An example of implementation of a finite-element method including a decorrelation step for the problems discussed in previous sections consists of the following steps:

1. Discretize the domain: define the elements and node connectivity, select the shape functions

2. Assemble system matrices $\mathbf{M}$ and $\mathbf{D}$

3. Implement the components of the boundary conditions (e.g. add Lagrange multipliers for periodic boundary conditions) and source terms (e.g. compute matrix $\mathbf{A_{D}}$ for the fluctuating source term) that are not time-dependent

4. Compute the decorrelation matrix $\mathbf{Q}$:

   1. Compute $\mathbf{U}_{\phi}$ and $\mathbf{\Sigma}_{\phi}$ through a singular value decomposition of $\mathbf{M}_{\phi\phi}$ (note that, since $\mathbf{M}_{\phi\phi}$ is symmetric and positive definite, the eigenvalues and singular values are identical)

   2. Compute $\Delta\widetilde{\widetilde{\mathbf{V}}}$ using Equation (64)

   3. Compute $\mathbf{A}_{\psi}$ as a diagonal matrix with diagonal elements given by $\sqrt{\Delta\widetilde{\widetilde{\mathbf{V}}}_{i}}$

   4. Compute the decorrelation matrix as $\mathbf{Q}=\mathbf{A}_{\psi}^{-1}\,\mathbf{U}_{\phi}\,\sqrt{\mathbf{\Sigma}_{\phi}}\,\mathbf{U}_{\phi}^{T}$

   5. Make $\mathbf{Q}$ sparse by setting elements with absolute value below a given tolerance to zero.

5. Initialize the solution $\widetilde{\mathbf{u}}^{0}$

6. Start the time-dependent loop, for $n=1$ to $n=N_{t}$ do:

   1. Compute the fluctuating source term

   2. Solve the linear system to obtain the finite-element solution array $\widetilde{\mathbf{u}}^{n}$ as a function of $\widetilde{\mathbf{u}}^{n-1}$ (see Section 2.3)

   3. Compute the mapped solution as $\widetilde{\widetilde{\mathbf{u}}}^{n}=\mathbf{Q}\,\widetilde{\mathbf{u}}^{n}$ and store it

4.3 Results using the spatial decorrelation matrix

Figure 9: Structure factor of the discrete solution after mapping with the decorrelation matrix computed with $N_{n}=51$ (dashed) and $N_{n}=101$ (dotted), and theoretical structure factor for the continuum equation (solid line). Top left: $\Delta t=$110-3$$, $N_{t}=$1106$$; top right: $\Delta t=$110-4$$, $N_{t}=$1106$$; bottom left: $\Delta t=$110-5$$, $N_{t}=$1107$$; bottom right: $\Delta t=$110-4$$, $N_{t}=$1107$$. All computations with $p1$-elements, semi-implicit time integration, $u_{0}=10000$.

Figures (9) and (10) show the results obtained applying the linear mapping (61) given by Equation (62) to the finite element results presented in Section 3. Only the results obtained with the semi-implicit integration scheme have been used. The discrete structure factor for the decorrelated results is computed as

$$S(k_{m})=\left<{\widetilde{\widetilde{\mathcal{U}}}}(k_{m}){\widetilde{\widetilde{\mathcal{U}}}^{*}}(k_{m})\right>\,,$$ (66)

with

$${\widetilde{\widetilde{\mathcal{U}}}}(k_{m})=\frac{1}{\sqrt{L}}\sum_{n=1}^{N_{n}-1}\widetilde{\widetilde{u}}_{n}\Delta\widetilde{\widetilde{L}}_{n}\exp\left(-\frac{i2\pi(m-1)(n-1)}{N_{n}-1}\right)\,,$$ (67)

where $\Delta\widetilde{\widetilde{L}}_{n}$ is the equivalent volume (equivalent length in 1D) related to shape function $\psi$ for node $n$, which can be computed using Equation (64).

The results presented in Figures (9) have been computed for two different meshes (number of nodes $N_{n}=51$ and $N_{n}=101$) with linear shape functions, and for different values of the time step $\Delta t$ and of the number of time steps $N_{t}$ for which solution data are collected. It can be seen that the linear mapping succeeds in yielding a structure factor that closely approximates that of the continuum ($S/u_{0}=1$). Once the artificial correlations are removed, and since the solution is completely decorrelated in space, refining the mesh does not improve the solution significantly, and, for the largest time step $\Delta t=$110-3$$ (see top left plot), refining the mesh actually increases the error, since the dimensionless time step $\beta=\Delta tD/\Delta x^{2}$ becomes larger. For small-enough values of the time step, increasing the number of time steps $N_{t}$ and the total time of the simulations $t_{\text{tot}}=N_{t}\Delta t$ is what contributes most to decreasing the error.

Figure (10) shows results obtained with both linear shape functions ($p1$-elements) and quadratic shape functions ($p2$-elements). Just like refining the mesh does not decrease the numerical error in this case, since the solution is white noise, increasing the order does not influence the accuracy of the results significantly.

Figure 10: Structure factor of the discrete solution after mapping with the decorrelation matrix computed with $N_{n}=51$ (dashed) and $N_{n}=101$ (dotted), and theoretical structure factor for the continuum equation (solid line). Left: $p1$-elements; right: $p2$-elements. All computations with semi-implicit time integration, $\Delta t=$110-4$$, $N_{t}=$1106$$, $u_{0}=10000$.

In B, complementary results for the computed dynamic structure structure of the mapped solution are shown, which illustrate that the method not only succeeds in removing artificial spatial correlations from the solution, but is also able to capture the time evolution of the solution.

5 FE solution of a fourth-order 1D diffusion problem

Following de2015finite , we consider next a stochastic diffusion equation of the form given by Equation (1) with

$$F(u)=u+\left(\frac{\ell_{0}}{2\pi}\right)^{2}\nabla^{2}u\,,$$ (68)

which leads to the fourth-order equation

$$\frac{\partial{u}}{\partial t}-D\nabla^{2}\left(u+\left(\frac{\ell_{0}}{2\pi}\right)^{2}\nabla^{2}u\right)=\nabla\cdot\left(\sqrt{2Du}\boldsymbol{\zeta}\right)\,.$$ (69)

Unlike Equation (4), the solution of which has only spatially decorrelated fluctuations, Equation (69) has a solution that exhibits spatial correlations with a finite correlation length defined by $\ell_{0}$. This allows us to investigate the effect of the mapping for different ratios of cell size to correlation length ${\Delta x}/{\ell_{0}}$.

5.1 One-dimensional boundary value problem

We consider the one-dimensional version of Equation (69)

$$\frac{\partial u({x},t)}{\partial t}-D\frac{\partial^{2}}{\partial x^{2}}\left(u({x},t)+\left(\frac{\ell_{0}}{2\pi}\right)^{2}\frac{\partial^{2}u({x},t)}{\partial x^{2}}\right)=\frac{\partial}{\partial x}\left(\sqrt{2Du_{0}}{\zeta}({x},t)\right)\,,$$ (70)

where $u_{0}$ and $\ell_{0}$ are constants. The rest of the variables are defined in the same way as for the second-order problem discussed in Section 3.

Instead of solving Eq. (70) directly, we replace the fourth-order PDE with the coupled set of second-order equations

$$\frac{\partial u}{\partial t}-D\frac{\partial^{2}}{\partial x^{2}}\left(u+\left(\frac{\ell_{0}}{2\pi}\right)^{2}w\right) =\frac{\partial}{\partial x}\left(\sqrt{2Du_{0}}\mathbf{\boldsymbol{\zeta}}\right)$$ (71a)

$$w-\frac{\partial^{2}u}{\partial x^{2}} =0\,.$$ (71b)

This allows us to use the same kind of discretization described in Section 2. Periodic boundary conditions are prescribed at the boundaries of the domain,

$$u(x=L,t)=u(x=0,t)\,,$$ (72)

$$w(x=L,t)=w(x=0,t)\,,$$ (73)

and the initial conditions are defined as

$$u(x,t=0)=u_{0}\,,$$ (74)

$$w(x,t=0)=0\,.$$ (75)

The theoretical static structure factor of the solution to this problem is given by

$$S(k)=u_{0}\frac{1}{1+\frac{k^{2}}{k_{0}^{2}}}\,,$$ (76)

with the cutoff wavenumber defined as

$$k_{0}=\frac{2\pi}{\ell_{0}}\,.$$ (77)

Unlike the case presented in Sec. 3, the covariance of the fluctuating solution to this equation for points $x_{1}$ and $x_{2}$ is not proportional to the delta function, but to an exponentially decaying function $\propto u_{0}k_{0}e^{-k_{0}(x_{1}-x_{2})}$ with characteristic length $\ell_{0}$ de2015finite .

This problem has a limit case when $\Delta x/\ell_{0}\gg 1$. In this case, the correlation length of the fluctuations is much smaller than the discretization length; this corresponds to the spatially decorrelated case described in Section 3, which needs the mapping proposed in Section 4 to remove the artificial correlations introduced by the discretization. The other limit case occurs when $\Delta x/\ell_{0}\rightarrow 0$. In this case, the discretization captures the correlation length of the fluctuations completely, the variance of the fluctuations not resolved by the mesh tends to zero, and, therefore, the direct FE solution without mapping converges to the continuum solution. In the following section, we discuss the results for different ratios $\Delta x/\ell_{0}$.

Like in Section 3, all the variables are made dimensionless variables, with the reference length $L=1$ and the reference time $T=L^{2}/D=1$.

5.2 Discretization and results

The simulation setup for this problem is based on the same discretization schemes and numerical procedures described in Sections 2 and 3.2 for the second-order boundary value problem. Figures 11 and 12 show results for the computed structure factor based on the direct FE results and on the results after the linear mapping to remove artificial correlations, for three different ratios $\Delta x/\ell_{0}=L/(N_{n}-1)/\ell_{0}$ and for both $p1$-elements (Figure 11) and $p2$-elements (Figure 12). The numerical results are close to the theoretical curve given by Eq. (76) for all discretizations when the linear mapping provided by Equation (62) is applied. The direct FE results without the mapping step tend to converge to the theoretical curve too as the mesh is refined for both $p1$- and $p2$-elements, with the results obtained with quadratic elements approximating better the theoretical curve for the structure factor than those obtained with linear elements. However, if no mapping is applied, a finer discretization is needed, as the characteristic length of the mesh needs to be significantly smaller than the typical correlation length to avoid the presence of significant artificial correlations. Therefore, applying the decorrelation matrix seems to be advantageous even in the case where there are spatial correlations over a finite length $\ell_{0}>0$, and its effect is larger for increasing $\Delta x/\ell_{0}$.

Figure 11: Structure factor of the FE results before (left) and after applying the decorrelation matrix (right): theoretical structure factor (solid line), $\Delta x/\ell_{0}=0.25$ (dash-dot), $\Delta x/\ell_{0}=0.5$ (dotted) and $\Delta x/\ell_{0}=1$ (dashed). The results have been computed with $p1$-elements, $\Delta t=$110-4$$, $N_{t}=$5106$$, $u_{0}=10000$.

Figure 12: Structure factor of the FE results before (left) and after applying the decorrelation matrix (right): theoretical structure factor (solid line), $\Delta x/\ell_{0}=0.25$ (dash-dot), $\Delta x/\ell_{0}=0.5$ (dotted) and $\Delta x/\ell_{0}=1$ (dashed). The results have been computed with $p2$-elements, $\Delta t=$110-4$$, $N_{t}=$5106$$, $u_{0}=10000$.

6 Discussion and conclusions

In this work, we have presented a numerical approach based on the finite element method to solve stochastic diffusion equations with thermal fluctuations. In particular, we have proposed a general formulation for the discrete fluctuating forcing term that is not restricted to any particular kind of element or shape function, and we have also derived a mass-conserving linear mapping to remove from the discrete solution the artificial correlations introduced by the spatial discretization.

The proposed mapping provides a solution that has a straightforward physical interpretation regardless of the spatial discretization, so it enables the use of finite element schemes without restrictions. Our analysis makes explicit the fact that the spatial discretization and the coarse-graining are two separate things, and that the artificial correlations due to the spatial discretization depend exclusively on the basis functions used to approximate the finite element solution. The mapping can be applied to any spatial discretization, including those defined on unstructured meshes. While the performance of the mapping has been demonstrated for linear equations, the approach is valid for nonlinear equations as well. In fact, the existence of an equivalent mapped solution that is free of artificial correlations can make treating non-linear terms in fluctuating-hydrodynamics equations easier.

In the case of problems for which the fluctuations are spatially correlated over a finite length, if the spatial discretization is refined enough with respect to the correlation length, the finite element method will eventually provide a solution with physically-meaningful spatial correlations. Nevertheless, the mapping succeeds in removing the artificial correlations related to underresolved fluctuations, thus providing physically meaningful results for coarser meshes.

Appendix A Implementation of the linearized stochastic forcing term

Based on Equation (16), an implementation of the stochastic forcing term may be proposed through a decomposition of the diffusion matrix, which for our standard Galerkin discretization is symmetric, and can therefore always be decomposed into

$$\mathbf{D}=\mathbf{A}_{D}\,\mathbf{A}_{D}^{T}\,,$$ (78)

where matrix $\mathbf{A}_{D}^{T}$ is the transpose of matrix $\mathbf{A}_{D}$. After integrating in time (see also Section 2.3), we can compute the array $\mathbf{f}^{n}$ at every time step $n$ as

$${f}^{n}_{i}=\sqrt{2u_{0}}\,a^{n}_{i}\,,$$ (79)

with array $\mathbf{a}^{n}=\left\{a^{n}_{i}\right\}_{i=1}^{N_{\text{dof}}}$ defined as

$$\mathbf{a}^{n}=\mathbf{A}_{D}\,\mathbf{z}^{n}\,,$$ (80)

where $\mathbf{z}(t)$ is an array of random numbers of length $N_{\text{dof}}$ following a Gaussian distribution of expected value $0$ and variance $1$. Matrix $\mathbf{A}_{D}$ can be computed once before the time-stepping starts, and then $\mathbf{f}^{n}$ can be computed at each time step based on it. The decomposition given by Equation (78) is not unique, but a convenient definition is

$$\mathbf{A}_{D}=\mathbf{U}_{D}\sqrt{\boldsymbol{\Sigma}_{D}}\mathbf{U}_{D}^{T}\,,$$ (81)

where $\mathbf{U}_{D}$ is a unitary matrix and ${\boldsymbol{\Sigma}_{D}}$ is a diagonal matrix, satisfying

$$\mathbf{D}=\mathbf{U}_{D}{\boldsymbol{\Sigma}_{D}}\mathbf{U}_{D}^{T}\,.$$ (82)

The definition of $\mathbf{A}_{D}$ provided by Equation (81) makes $\mathbf{A}_{D}$ symmetric. Moreover, while in general a matrix $\mathbf{A}_{D}$ satisfying Equation (78) will be a dense matrix, the definition provided by Equation (81) makes it possible to approximate it by a sparse matrix, which is important to preserve the computational efficiency of the finite element method.

Appendix B Dynamic structure factor

In this appendix, we present complementary results to those presented in Sections 3 and 4.3. In particular, we present results for the dynamic structure factor, which, unlike the static structure factor discussed in previous sections, contains information about the time evolution of the solution. The dynamic structure factor is defined as the Fourier transform of the time-dependent correlation function. For the problem discussed in Section 3, the analytical expression for the dynamic structure factor is given by

$$S_{\text{dyn}}(k,\tau)=\left<\mathcal{U}(k,t)\mathcal{U}(-k,t-\tau)\right>=u_{0}e^{-Dk^{2}\tau}\,.$$ (83)

Figure 13 shows the dynamic structure factor corresponding to the finite element solution (before mapping) for two values of the wavenumber $k$ as a function of the time lag $\tau$. As also explained in Section 3, for a given value of $k$, the results converge to their theoretical value as the mesh is refined.

Figure 13: Dynamic structure factor of the finite-element solution (before applying the decorrelation matrix) as a function of the time lag for $k=20\,(2\pi/L)$ (left) and $k=30\,(2\pi/L)$ (right): computed with $N_{n}=51$ (dashed), $N_{n}=101$ (dotted) and $N_{n}=201$ (dash-dot), and theoretical dynamic structure factor (solid line). All computations with $p1$-elements, $\Delta t=$110-5$$, $N_{t}=$1106$$, $u_{0}=10000$, $\alpha=1/2$.

Figure 14 displays the dynamic structure factor corresponding to the mapped finite element solution for the same two values of $k$. A key difference with respect to the finite element results before mapping (Figure 13) is that the mapped results present a good agreement with the theoretical value for $\tau$ values close to zero regardless of the size of the mesh. For larger values of $\tau$, the dynamic structure factor of the mapped results converges as well towards the theoretical one as the mesh is refined.

Figure 14: Dynamic structure factor of the mapped solution (after applying the decorrelation matrix) as a function of the time lag for $k=20\,(2\pi/L)$ (left) and $k=30\,(2\pi/L)$ (right): computed with $N_{n}=51$ (dashed), $N_{n}=101$ (dotted) and $N_{n}=201$ (dash-dot), and theoretical dynamic structure factor (solid line). All computations with $p1$-elements, $\Delta t=$110-5$$, $N_{t}=$1106$$, $u_{0}=10000$, $\alpha=1/2$.

Figures 15 and 16 show the effect of the time step $\Delta t$ and of the total simulated time $N_{t}\Delta t$ on the computed structure factor based on the mapped solution, and illustrate the convergence of the numerical solution towards the theoretical one as the time step decreases and the total simulated time increases. For values of $\tau$ below a threshold, decreasing the time step leads to a better agreement with the theoretical solution, as seen in Figure 15. For values of $\tau$ above a certain threshold, the value of the dynamic structure factor is small, and the residual error due to the finite time $N_{t}\Delta t$ that is simulated dominates the solution. As shown in Figure 16, this error decreases as the total simulated time $N_{t}\Delta t$ increases.

Figure 15: Influence of the time step on the dynamic structure factor of the mapped solution for $k=20\,(2\pi/L)$: computed with $\Delta t=$110-5$$ (dashed) and $\Delta t=$510-5$$ (dotted), and theoretical dynamic structure factor (solid line). All computations with $p1$-elements, $N_{n}=201$, $N_{t}\Delta t=10L^{2}/D$, $u_{0}=10000$, $\alpha=1/2$.

Figure 16: Influence of the total simulation time on the dynamic structure factor of the mapped solution for $k=20\,(2\pi/L)$: computed with $N_{t}\Delta t=10L^{2}/D$ (dashed) and $N_{t}\Delta t=100L^{2}/D$ (dotted), and theoretical dynamic structure factor (solid line). All computations with $p1$-elements, $N_{n}=201$, $\Delta t=$110-5$$, $u_{0}=10000$, $\alpha=1/2$.

These results show that the numerical approach based on applying the decorrelation matrix to the finite element results, as proposed in Section 4, can capture the time evolution of the solution.