## 1 Introduction

Models describing multiphysics and multiscale processes are ubiquitous in numerical simulations. Classical applications include those arising in aeronautics, meterology, biology, material and environmental sciences which are modeled by Navier-Stokes [30, 35], reaction-diffusion [32, 34] or ADR equations [2, 31]. The individual physics or scale components have very different properties that are reflected in their discretization; for example, for ADR systems, the discrete advection has a relatively slow (or ‘nonstiff’) dynamics while the diffusion and chemistry are fast (or ‘stiff’) evolving processes [12, 33]. The discretization in time of slow processes is advantageous with an explicit method since the nonstiff terms are (often) nonlinear and computationally more expensive to evaluate whereas implicit methods are more appropriate for stiff processes because of their favorable stability properties, allowing one to select longer time steps [40]. IMEX integrators have been proposed as an attractive alternative (compared with fully explicit or fully implicit time integration methods) where one combines the implicit (explicit) integration for the fast (slow) scale [7]. IMEX linear multistep methods are investigated in [11]

, while IMEX Runge-Kutta schemes have been classified in

[1]. The literature on higher order IMEX methods are limited since these are difficult to construct due to large number of order constraints, increasing stability restrictions with increasing order of accuracy [3] as well as the malaise of the factitious ‘computational mode’ that all multi-time level discretization schemes, with three or more time levels suffer [30].Similarly, the spatial discretization schemes must resolve all scales present in the flow for accurate numerical solution [17]. Additionally, the numerical propagation properties of every individually resolved scale (i. e., the numerical group velocity and the numerical phase speed) must be the same as the physical propagation properties of the respective scale, so-called the DRP property [25]. The Fourier spectral methods are known for their ability to provide numerical accuracy and dispersion error-free results with all scales resolved and aliasing problems avoided [9]. But the application of such methods are largely limited to simple geometries and boundary conditions with uniformly spaced nodes [4], with some isolated cases attempted in complex domians [8, 19]. Alternative to spectral method are the higher order explicit discretization methods [38]

and methods based on Padé approximation

[15], synonymous with the compact schemes. All compact schemes for evaluating theth derivative of an unknown vector

require solving the set of linear algebraic equations given by or(1) |

Compact schemes have excellent resolution property [26], even when their formal order of accuracy is lower than the traditional explicit discretization methods [28]. Although compact schemes have fewer points in the stencil, their implicit nature gives a global approximation of the derivatives [30]. Other relevant numerical properties of compact schemes have been analyzed earlier, including the DRP which is essential in solving time-dependent problems accurately [27], the de-aliasing or the unphysical pile-up in the energy spectrum for high wavenumbers [29] and higher spectral resolution on coarser meshes [21]. However, in this article we examine the dispersion properties of the combined, spatiotemporal IMEX-Compact schemes using the error propagation equation, described next.

Traditional stability analysis attributed to von Neumann for linear problems assume error to follow the same dynamics as given by the governing difference equation for the signal [5]. Recent studies have identified mechanisms for solution breakdown of linear equation to be related to the numerical properties of discrete computing, i. e., the dissipation, dispersion and phase errors, and that the governing error propagation equation is different from the governing dynamics of the signal [25, 22]. The numerical schemes for non-periodic problems discussed in §2, require a similar analysis and for this purpose we have selected the 1D linear ADR equation,

(2) |

where are constant advection speed, diffusivity and growth rate, respectively. The technique for analyzing discrete computational methods using Fourier-Laplace transforms is invoked in this article [25], where the unknown is expressed in the wave number () plane by, , such that an amplification factor can be introduced as . For direct numerical simulations (DNS) one must have a neutrally stable method (). If we define the numerical solution of equation (2) by and the numerical error as

, then the correct error propagation, partial differential equation (PDE) is given by (refer §

6 for the detailed derivation),(3) |

where denotes the Fourier amplitude of the initial condition, is the amplification factor of the numerical discretization method and are the temporal and the spatial step-size, respectively. One notes that equation (3) clubs error based on generic stability (through the numerical amplification factor) and the DRP properties of the method of discretization (e. g. the scaled group velocity and the absolute error in the phase speed), unlike the one obtained using the modified equation approach [10], where the truncation error is accounted by collating and representing the discretized terms in the difference equation by their equivalent differential forms.

The aim of this article is to investigate the DRP properties of four, high accuracy (two-time level) spatiotemporal discretization using the state of the art error propagation equation. The two-time stepping methods are selected since they are devoid of the hazards of the contrived computational modes [30] and the complication related to the linear resonance stability [36]. §2 details the spectral analysis of these four numerical methods including the Explicit-OUCS3-CD scheme (§2.1), the Implicit-OUCS3-Lele scheme (§2.2), the IMEX-OUCS3-Lele scheme (§2.3) and the IMEX-NCCD scheme (§2.4). §3 examines the DRP properties of these four discretization techniques using the error propagation equation (3). The spectral features of these four numerical methods (and in particular the IMEX-NCCD scheme) is benchmarked via the solution of the nonlinear, parabolic-elliptic PKS chemotaxis model in §4. §5 concludes with a brief discussion of the implication of these results.

## 2 Spectral analysis for linear 1D ADR equation

The basis of the spectral analysis for the numerical methods formulated for equation (2), is to transform the information from the physical to the spectral plane using its Fourier-Laplace transform. The spectral resolution of the discretization methods discussed below (i. e., the absolute value of the ratio of the numerical to the physical amplification factor, , the scaled group velocity, , and the absolute value of the error in the phase speed, , refer §6) is investigated for all wavenumbers within the range where the upper bound of this range is determined by the Nyquist criterion [20]. is the grid spacing. The non-dimensional parameters introduced in the analysis are the Courant-Friedrichs-Lewy (CFL) number, , the Peclet number, , and the Damköhler number, . For illustration purpose, and is fixed, the domain is discretized using equidistant points and the plots in §2.1-2.4 are shown for an interior node, .

### 2.1 Explicit-OUCS3-CD

In the first numerical method, the temporal discretization of equation (2) is achieved via an explicit two-stage Runge Kutta/Heun’s method (also known as the explicit trapezoidal rule) [2]. The advection term is spatially discretized using the OUCS3 [26] and the diffusion component is discretized by the central CD scheme.

Assuming a domain with equidistant points with spacing where a function

is defined, the second order accurate OUCS3 estimates the first spatial derivative of the solution at the

th node, [30]. The inner stencil of the scheme is given by(4) |

where , , , with the coefficients and . For an improved numerical stability, the one-sided boundary stencil () and the near boundary stencil () are proposed as follows,

(5a) | |||

(5b) | |||

(5c) | |||

(5d) |

where and . The point in all of the compact schemes discussed in this and the next sections is of little concern since the governing PDF is not discretized at this node, only the boundary condition is specified there. In OUCS3, the point shows instability across all wavenumbers and therefore this calculated derivative is eventually replaced by the CD scheme locally. The solution of equation (2) at the th node () at each stage of the temporal discretization as a function of the solution at previous time-step, , is outlined as follows,

(6a) | |||

(6b) |

where is the matrix of the OUCS3 scheme given by equations (1), (4), (5). The numerical amplification factor, , is given by,

(7) |

where

(8) |

In figure 0(a), the amplification ratio, , are plotted as contours in the indicated range of and . One notes that in the limit of vanishingly small values of and , the stability of this numerical scheme is reaction-rate dependent (i. e., in the limit, and , we have , refer equations (7, 8, 31). However, for intermediate values of the CFL number, , this numerical method is stable (or ). In figure 0(b), the scaled group velocity contours show significant dispersion effects, for almost all the selected values of and (i. e., ), that would invalidate the long time integration results even when the neutral stability (i. e., ) is ensured for those wavenumbers and CFL numbers. In fact the numerical solution travels in the wrong direction (or ) within the range or within leading to numerical instabiltites in the form of waves (detailed analysis in §3) [39]. Although the error in the phase speed disappear for very small values of (as highlighted in figure 0(c)), we conclude through figures 0(b), 0(c) that the combined effects of the dispersion errors cannot be simply eliminated or reduced by grid refinement, as suggested in [10].

### 2.2 Implicit-OUCS3-Lele

In the second case, equation (2) is numerically time-integrated using the implicit mid-point rule [36]. The advection and the diffusion terms in equation (2) are spatially discretized using the OUCS3 scheme, equations (4, 5) and the Lele’s scheme, respectively. The interior, boundary and the near boundary stencils of the Lele’s scheme for the second derivative [15] are given by,

(9a) | |||

(9b) | |||

(9c) | |||

(9d) | |||

(9e) |

Lele’s compact scheme [15] has superior resolution (e. g. when compared with the repeated operation of evaluating the first derivative twice using OUCS3) and has lower dispersion error for all the nodes [30]. The solution of equation (2) at the th node, , is expressed as a function of the known at the previous time-step, , as follows,

(10) |

where and are the matrices of the OUCS3 scheme (equation (4, 5)) and the Lele’s scheme (equation (9)), respectively. The numerical amplification factor is specified as follows,

(11) |

The contours of the amplification ratio are outlined in figure 1(a) in the selected range of and . This discretization scheme has a comparatively larger stability region than the explicit scheme discussed in §2.1 (i. e., , for all and ), a feature attributed due to the favorable temporal stability property of implicit methods for stiff equations [36]. The scaled group velocity contours in figure 1(b) indicate that the region of spurious propagation (via waves) is limited within the range (for all ) or within (and ). Further, the dispersion errors via the scaled group velocity is inconsequential in the limit of small wavenumbers (i. e., V in the limiting case of ). Figure 1(c) suggests that the error in the phase speed vanishes in limit of negligibly small wavenumbers (i. e., when ), an observation parallel to the numerical method discussed in §2.1 (refer figure 0(c)). Overall, the stability and the propagation characteristics indicate that this discretization scheme is marginally better than the first case in determining the numerical solution of the 1-D linear advection-diffusion-reaction equation (2).

### 2.3 IMEX-OUCS3-Lele

In the third case, the temporal discretization of equation (2) is accomplished via a combination of the explicit two-stage Runge Kutta/Heun’s method and the implicit mid-point rule [36]; the (nonstiff) advection term is treated explicitly while the (stiff) diffusion and reaction terms are handled implicitly. Spatial discretization is achieved through the OUCS3 (Lele) for first (second) derivatives, i. e. equation (4, 5) (equation (9)). Using this combination of spatiotemporal discretization, the unknown solution of equation (2) at the th node, , at each stage of the temporal discretization as a function of the known variable at the previous time-step, , is described as follows,

(12a) | |||

(12b) |

where and are the matrices of the OUCS3 scheme and the Lele’s scheme, respectively. The numerical amplification factor is prescribed by,

(13) |

where

(14) |

The amplification ratio contours in this case indicate a region of stability, and (figure 2(a)), while the spurious propagation features (i. e., , figure 2(b)) is restricted inside the domain, or within the region, . As expected, the phase speed error in this case is inconsequential for small wavenumbers ( in the limit , figure 2(c)). In summary, this numerical method is comparable to the second case (§2.2) but with slightly reduced CFL range for temporal stability.

### 2.4 Imex-Nccd

In this numerical method, the IMEX temporal discretization (refer §2.3) combined with the NCCD scheme for the spatial discretization of the advection and the diffusion terms is utilized to analyze equation (2) [27, 29]. The NCCD scheme is used to simultaneously evaluate the first and the second spatial derivatives () from the following discrete equations for the boundary and the interior nodes in terms of the known values, ,

(15a) | ||||

(15b) | ||||

(15c) | ||||

(15d) | ||||

(15e) | ||||

(15f) |

To avoid the instability and the attenuation near the inflow and the outflow, the first derivative of the unknowns near boundary () are eventually replaced with the locally explicit stencil (5c, 5d), while the second derivative of the unknowns near boundary () are replaced with [30],

(16a) | |||

(16b) |

Writing the NCCD interior and boundary stencils given by equations (5, 15, 16) as a system of linear algebraic equations,

(17a) | |||

(17b) |

and on solving equations (17a, 17b) simultaneously we arrive at,

(18a) | |||

(18b) |

where,

(19a) | |||

(19b) |

The numerical amplification factor is evaluated using equations (13, 14) with replaced with the matrices of the NCCD scheme given by equations (19a,19b), respectively.

The superior numerical resolution feature of this method is evident from the amplification ratio contours in figure 3(a), this scheme is neutrally stable (i. e., ) within the region, - a property which is absolutely essential in Direct Numerical Simulation (DNS) of multiscale models [4]. The spurious propagation characteristics (or the region of negative group velocity, figure 3(b)) is limited to a relatively smaller range, , while the region where the dispersion effects due to phase speed error (figure 3(c)) is diminished in comparison with the other three numerical methods (§2.1-2.3). To recapitulate, the IMEX-NCCD scheme shows excellent DRP properties and it is undoubtedly the preferred numerical method to solve the 1D linear advection-diffusion-reaction equation (2), a conclusion which is further elucidated via the analysis through the error propagation equation, described next.

## 3 Error dynamics

In this section, we identify the reasons for the dispersion error in the computational solution of equation (2), solved via the numerical methods discussed in §2. The analyses in this section is based on the error propagation equation (34). An extreme form of the dispersion error has been identified as the waves [39] which are small amplitude, parasitic waves propagating in the direction opposite to the physical direction of propagation [23], leading to numerical instabilities as well as unphysical bypass transition [21]. If is negative the numerical waves propagate upstream despite the physical requirement of downstream movement (since the physical advection speed, in equation (2)). In other words, is the necessary condition for the generation of waves. Further elaboration on the role of waves in the solution of equation (2) is considered through the propagation of a wave-packet which is given by the following initial condition,

(20) |

where controls the width of the packet and x, and is the center of the wave-packet, the central wavenumber and half domain length, respectively. The in silico analyses of equation (2) together with the initial condition (20) and Dirichlet boundary condition (, equation (20)) is executed with the fixed parameters . The results reported in §3.1-3.3 are presented at (non-dimensional) simulation time , using a time-step of numerical integration, and (with the exception of the results described in §3.2) the domain is discretized uniformly with points. Two trials cases are considered where the width of the initial wave-packet is (a) compact, , and (b) outstretched, . The initial conditions (equation (20)) along with their Fourier spectra are displayed in figures 4(a) and 4(b), respectively.

### 3.1 Role of initial conditions

The solution of equations (2, 20) using the Explicit-OUCS3-CD scheme (§2.1) for () is shown figure 4(c) (figure 4(d)). The reason behind the negligible presence of waves for the case of compact initial wave-packet (figure 4(c)) versus a substaintial presence of the same for the outstretched initial wave-packet case (i. e., the boxed region, figure 4(d)

) can be explained via the respective Fourier transform of the initial condition (figure

4(b)). From figure 0(b), we have noted the region of negative scaled group velocity for this numerical scheme to be (and ). The Fourier spectra reveals that the amplitude of the outstretched initial wave-packet (dashed curve, figure 4(b)) is significantly large within this range of wavenumbers, , and it is these components of initial conditions which are responsible for waves. The examples in figure 4(c), 4(d) are depicted for fixed value of the parameters, , and for these values the Explicit-OUCS3-CD scheme predicts the scaled group velocity, . Thus, the numerical solution propagates at a speed slower than the exact solution while, for certain values of wavenumbers within the range (and ), the waves travel at speed faster than the exact solution (refer figure 0(b) for the contour values of within the wavenumber range ) leading to the oscillations as shown in figure 4(d).### 3.2 Role of grid resolution

Spurious generation of flow structures in 1D advection equation, due to inadequate grid resolution leading to the creation and propagation of waves has been communicated in [23]. A similar analysis the 1D ADR equation is outlined in this section. Figure 4(c), 4(e) presents the solution of equations (2, 20) using the Explicit scheme (§2.1) for with the domain discretized uniformly using and points, respectively, such that the central wavenumber of the initial wave-packet is fixed (equation (20)). Although the central wavenumber of the initial wave-packet is identical, the band-width of the waves resolved in both cases are different. A coarser grid (i. e., ) is unable to capture the initial conditions accurately (solid curve, figure 4(a)) resulting in a larger range of excited wavenumbers (as compared with the finer grid), in particular those wavenumbers that correspond to negative scaled group velocity thereby generating large amplitude waves. The in silico studies reported in [23] emphasize that merely introducing finer mesh is not adequate for accurate computing and one must be punctilious towards the basic DRP properties of the spatiotemporal discretization technique, as discussed in the next section.

### 3.3 Role of different numerical methods

Figure 4(c), 4(f) highlights the numerical solution of the 1D ADR equation (2) for an initial wave-packet of width (equation (20)) using the Explicit scheme (§2.1) and the IMEX-NCCD scheme (§2.4), respectively. At initial excitation (and ), the Explicit scheme estimates the propagation characteristics of the signal poorly (i. e., ) while the IMEX-NCCD scheme predicts this value relatively accurately (). Another source of error predicted by the error propagation equation (34) is the absolute error in the phase speed, which is observed to be ) for the Explicit-OUCS3-CD (IMEX-NCCD) scheme. A substantial phase error of the Explicit scheme is responsible for the dispersion effects resulting in an asymmetry in the computed solution of equation (2) (figure 4(c)) albeit the exact solution is symmetric (figure 4(a)). The amplification ratio of the initial excitation is estimated to be for the Explicit-OUCS3-CD / IMEX-NCCD scheme, respectively. We note that the solution obtained by IMEX-NCCD scheme has negligible waves (figure 4(f)). The presence of waves implies that the energy of the waves (or physical waves) is siphoned off to the waves ensuing a diminished amplitude of the numerical solution (e. g., compare the amplitude of the computed solution of the Explicit scheme (figure 4(c)) versus the IMEX-NCCD scheme (figure 4(f))). The fact that the amplification ratio of the IMEX-NCCD scheme remains very close to one indicates the near absence of waves.

We summarize our discussion by indicating two sources of error which are particularly perplexing in the DNS of ADR equations, the existence of waves for those set of numerical parameters for which the spatiotemporal discretization is neutrally stable. In such a situation the waves do not attenuate and have to be eliminated by deploying an explicit filter [24]. Another aspect of dispersion error is related to the Gibbs’ phenomena which occurs as a consequence of sharp discontinuity in the solution and which causes fictitious oscillations [25, 22], a problem which can be remedied using the high accuracy DRP methods. Both these sources of error are discernible in nonlinear ADR equation whose solution (for a particular case) is explored next.

## 4 IMEX method for chemotaxis model

In this section we demonstrate the suitability of the IMEX-NCCD method developed in §2.4 for solving the PKS chemotaxis model [18]. Chemotaxis refers to the mechanism by which cellular motion occurs in response to an external, chemical stimulus [13]. The PKS model is governed by a system of ADR equations, which in the classical 2D case reads as,

(21) |

where denotes the cell density and the chemoattractant concentration, respectively. is the chemotactic sensitivity parameter. and represents the 2D gradient and the Laplacian operators, respectively. A common property of all chemotaxis systems is their ability to accurately capture concentrations that mathematically result in a rapid growth in small neighborhoods of the solution. For example, the solution of equation (21) may ‘blow up’ in finite time, provided that the initial mass of the cells, , is above a critical threshold which is for radially symmetric cases [14]. Steep gradients and spikes may give rise to nonphysical oscillations if the numerical scheme is not guaranteed to satisfy the discrete maximum principle (DMP), resulting in negative cell densities [37]. To capture the singular, spiky behavior of the solution of (21), a high accuracy, positivity-preserving, Finite Volume (FV) scheme (satisfying DMP) is deployed by rewriting equation (21) as follows,

(22) |

where and are the chemotactic velocities. The subscripts in equation (22) denote the partial derivative with the respective variables. We consider the model (22) in a square domain, , and introduce a Cartesian mesh consisting of cells which are assumed to be of uniform size (), such that for all and or all . On this mesh, a general FV scheme for the PKS system (22) will have the following form,

(23a) | ||||

(23b) |

The values and are the approximate point variable values at the cell center, is the discrete Laplacian operator, and are the numerical fluxes in the x- and the y- direction, respectively, and these are evaluated at the cell edges. The derivatives, and are approximated by appropriate compact / explicit finite difference (FD) scheme at the cell centers and then linearly extrapolated at the cell edges (refer §2 for the details of the spatial discretization utilized). The point values and are computed as follows,

(25a) | ||||

(25b) |

where the positivity-preserving generalized minmod limiter, given by

(26) |

ensures the positivity of the reconstructed point values. The parameter in equation (25) controls the amount of numerical viscosity present in the scheme, larger values of correspond to less dissipative but usually more oscillatory reconstructions. An explicit FV-FD numerical scheme for solving the 2-D Keller-Segel model is outlined in [6].

### 4.1 Numerical results

An initial exploration to test the efficacy of the IMEX-NCCD method was conducted using an example taken from [6]. We consider the initial-boundary value problem for the PKS system (21) with and (refer (25)) on a square domain subject to a radially symmetric gaussian initial condition:

(27) |

As it was highlighted in [6], the solution with the above mentioned initial data is expected to develop a like singularity at the center of the computational domain in a very short time. Figure 6 compares the cell densities computed via the FV-Explicit-OUCS3-CD method (figures 5(a), 5(c)) versus the FV-IMEX-NCCD scheme (figures 5(b), 5(d)) on a uniform mesh with , at two different computational times (Figures 5(a), 5(b)) and (Figures 5(c), 5(d)), respectively. In both of these cases the zero Neumann boundary conditions are used, which is implemented using the standard ghost point technique when the CD discretization is used. Further, the time steps and were chosen for the respective methods to ensure the positivity of the computed densities. As outlined in figure 6 although both methods capture the ‘needle-like’ structure, the blowup phenomena is better resolved by the IMEX-NCCD method since the solution obtained via this method is oscillation-free (compare the figure insets of figure 5(a) versus figure 5(b) or that of figure 5(c) versus figure 5(d)). This is because, centered approximations (e. g. the CD discretization) for hyperbolic PDEs lead to well-known difficulties with nonphysical oscillations (or Gibbs’ phenomena) arising near discontinuities or steep gradients [16]. A detailed study relating the origin of these oscillations with the time-step of the numerical methods will be reported later.

## 5 Summary and conclusion

This introductory investigation compares and reports the dispersion error of four, two-time level, DRP schemes using the novel error propagation equation for 1D linear ADR equation (details given in §6). §2 presented a comprehensive spectral analysis of these four numerical methods, highlighting regions of temporal stability, positive scaled group velocity and minimal phase errors in the wavenumber () - CFL number () plane. The role of the initial conditions, inadequate grid resolution or the choice of the numerical methods in the generation of the dispersion errors via waves is outlined in §3. The superior spectral resolution of the IMEX-NCCD method is then benchmarked via the solution of the nonlinear, parabolic-elliptic PKS chemotaxis model in §4. The blowup phenomena in the numerical solution of the PKS model was accurately captured by the IMEX-NCCD scheme, devoid of any amplitude fluctuations. In a nutshell, the numerical solution obtained by the IMEX-NCCD method has negligible dispersion error due to waves and is ideally suited to arrest the oscillations via the Gibbs’ phenomena.

## 6 Appendix A: Derivation of the error propagation equation for 1D linear ADR equation

The analysis of the space-time discretization of the 1D linear ADR equation (2), with initial condition given by, for , is used as a model for PDE replicating multiscale processes. Using the spectral representation , equation (2) is transformed as

(28) |

where is the wavenumber. Using the transformed initial condition, , one obtains the solution of equation (28) as

(29) |

If we represent by its Fourier-Laplace transform as,

Comments

There are no comments yet.