High Order Asymptotic Preserving and Classical Semi-implicit RK Schemes for the Euler-Poisson System in the Quasineutral Limit

The rest of this section is devoted to a brief summary of some analytical results pertaining to the quasineutral limit of the EP system (2.1)-(2.3), deemed relevant for the analysis to be carried out later, along the lines of [25]. As discussed in the introduction, the multi-scale nature of the EP system is elucidated through the parameter $\varepsilon$, when there is a disparity in the scale of the physical Debye length compared to the macroscopic length scale. The quasineutral limit model is obtained by letting $\varepsilon\to 0$ in the system (2.1)-(2.3). Further, when $\varepsilon\ll 1$, the flow regime can be considered to be close to the quasineutral limit.

Let us assume that the dependent variables $\rho,q$ and $\phi$ yield the following limits as $\varepsilon\to 0$:

Passing to the limit $\varepsilon\to 0$ in (2.1)-(2.3), cf. [25], we obtain the ensuing incompressible Euler system as the quasineutral limit of the Euler-Poisson system,:

Based on the above result we put together the following definition, which also concerns itself with the choice of the appropriate set of initial conditions for the EP system (2.1)-(2.3) and, reveals the multi-scale nature of the solution in accordance with the formal limits considered in (2.4).

In the incompressible Euler limit system (2.5)-(2.7), the hydrodynamic pressure is identified as the negative of the limiting electric potential $\phi_{(0)}$. It is a mixed elliptic-hyperbolic system for the unknowns $u_{(0)}$ and $\phi_{(0)}$. The electric potential $\phi_{(0)}$ plays the role a Lagrange multiplier for the incompressibility constraint (2.5). Taking divergence of the momentum equation (2.6), and using the divergence-free condition (2.5) on the velocity $u_{(0)}$, see also [15]; yields the following elliptic equation for $\phi_{(0)}$:

The equation (2.10) can be used to eliminate the divergence constraint (2.5). In fact, the divergence-free condition (2.5) needs to be satisfied only at time $t=0$. We summarise a few results from [25] as follows.

As is shown in the case of the limiting incompressible system (2.11)-(2.14), similarly, the scaled Euler-Poisson sytem (2.1)-(2.3) can be also shown to have a reformulated counterpart; see [25] for a formal derivation of the same.

To elucidate the parallel between the EP system and index-1 DAEs, we first linearise the EP system (2.1)-(2.3), following the lines of [26], about a stationary homogeneous state $\bar{\rho}=1,\ \bar{q}=0,\ \partial_{x}\bar{\phi}=0$. Let the perturbations from the considered stationary state be represented by an infinitesimal $\delta\ll 1$. An asymptotic expansion of the unknowns around the constant state with respect to $\delta$ is given by;

Considering only the first order terms in $\delta$, and dropping the primes for the sake of convenience, we obtain the following linearised EP system:

where $c_{s}=\sqrt{p^{\prime}(1)}$ is the speed of sound. We take the partial Fourier transform, in space of the equations (3.2)-(3.4) to obtain the following system of equations:

Here, $\hat{\rho},\hat{u}$ and $\hat{\phi}$ are the Fourier transformed variables corresponding to $\rho,u$, and $\phi$, respectively. Furthermore, to comprehend the analogy between the Fourier transformed linearised EP system (3.5)-(3.7), the Fourier space variable $\xi$ is treated as a fixed parameter.

Motivated by the conclusion of the above proposition for the linearised EP system (3.2)-(3.4), we rewrite the non-linear one-fluid EP system (2.1)-(2.3) as a set of DAEs via the following definition.

We underline the fact that the method of lines approach is used to write the system of PDEs (2.1)-(2.3) in a DAE form given by (3.13)-(3.14). Analogous treatment of PDEs using the method of lines philosophy to obtain a system of ODEs is a regular practice in the literature; see, e.g. [1, 3, 68] and the references therein for examples.

Having recast the EP system as a set of DAEs in (3.13)-(3.14), we now turn towards the numerical approximation. It has to be noted that even though implicit RK schemes were primarily designed for stiff systems of ODEs, these methods can be systematically extended to even broader classes of equations, such as DAEs. What follows is a brief exposure to the so-called indirect approach [40] to design implicit RK schemes for index-1 DAEs. Dedicated to this aim, we consider an arbitrary set DAEs of index 1 in the form

where $f$ and $g$ are sufficiently smooth functions, and $y$ and $z$ are vectors of appropriate dimensions with respect to the functions $f$ and $g$.

We now consider an $s$-stage implicit RK method defined by the triple $(A,c,\omega)$, where $A=(a_{k,l})$ are the coefficients, $c=(c_{k})$ are the intermediate times and $\omega=(\omega_{k})$ are the weights. The indirect approach [40] defines an RK scheme for (3.17)-(3.18) with intermediate stages $(Y^{k},Z^{k})$ in the following way:

where the last row of the coefficient matrix $A$ is taken to be equal to the vector $\omega$; in other words, the implicit RK scheme is stiffly accurate (SA); see [40]. Hence, the numerical solution $(y^{n+1},z^{n+1})$ is defined by $y^{n+1}=Y^{s}$ and

Note that the indirect approach enables us to ensure that the manifold

is preserved by all the intermediate stages and automatically by the update stage. For a detailed discussion on the derivation, validity and stability analysis of the direct approach of implicit RK schemes applied to index-1 DAEs, interested readers may refer to [40].

The indirect approach presented in [40], and briefed above, in its core, takes implicit RK schemes for ODEs and derives implicit schemes for DAEs. Yes indeed an implicit scheme does fair well mathematically in the sense that it yields the asymptotic properties, such as stability and consistency which leads to the AP property. But, as a drawback, because of the fully implicit nature of the method, it becomes computationally very expensive when tackling problems, especially with nonlinearities, such as the EP system (3.13)-(3.14). In order to reduce the overdose of implicitness offered by fully implicit RK schemes via the indirect approach platform (3.19)-(3.21), in the following, we consider two semi-implicit RK approaches. The ultimate aim is to develop RK time steppings for index-1 DAEs and in turn for the EP system (3.13)-(3.14) as well, with an optimal portion of implicitness.

Semi-implicit or IMEX schemes trade unconditional stability of fully implicit schemes in order to draw some gains in the computational efficiency front. IMEX-RK schemes are a staple for the CFD community; see, e.g. [1, 3, 7, 8, 9, 43, 66, 68, 78] for IMEX schemes applied to tackle problems plagued with stiffness. Introduction of an adequate dosage of implicitness through IMEX schemes provides an additional amount of stability which is usually required for methods resolving singular parameters, while maintaining a pardonable level of computational performance. In this work we consider two types of IMEX-RK schemes, namely

• •

  additively partitioned IMEX-RK schemes [2, 67];

• •

  semi-implicit IMEX-RK schemes [9].

The former one, the class of additively partitioned IMEX-RK schemes, is primarily designed for stiff ODEs in which the identifiable stiff and non-stiff parts are summed together. To this end, consider the following system of additively partitioned stiff ODEs:

where $y=y(t)\in\mathbb{R}^{n},\ f,g\colon\mathbb{R}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$, and $\varepsilon>0$ is a small parameter known as the stiffness parameter. Here, $f$ and $g$ are referred to as, respectively, the non-stiff and stiff parts.

There exist a large section of stiff problems amongst ODEs wherein identifying an additive partition of the right hand side into stiff and non-stiff subparts may not be possible. There are even systems where it is possible to find an additive splitting of the stiff terms, however, stiffness can be associated with a nonlinear term. Consequently, in this scenario the implementation of additive IMEX-RK schemes may lead to a need for solving highly nonlinear equations, which in-turn requires the usage of expensive, iterative Newton solvers.

In the light of the aforestated predicaments the second type of schemes that we consider here, i.e. the semi-implicit IMEX-RK schemes, have been recently presented in [9]. The authors have devised an efficient platform to derive IMEX-RK schemes for problems where a multiplicative stiff-nonstiff splitting can be recognised. Further, the stiffness can be associated to some part of the unknown vector, and rest of it might be non-stiff. In such scenarios, these schemes have an advantage that either a linear solver can be employed or a Newton iteration for a problem with reduced non-linearity can be engaged; see [10, 11, 33] for some related works.

Consider a general stiff-ODE in the form

where in $H$ the first argument represents the non-stiff part of $y$ and the second argument the stiff-part. Note that the right hand side is assumed to have a stiff dependence only through the last argument. In order to distinguish these, we use subscripts and write

It has to be noted that both $y_{E}$ and $y_{I}$ represent the same variable $y$ itself. The subscripts are used to only distinguish between the stiff and non-stiff components.The class of problems characterised by the system (3.26)-(3.27) is said to be in a generalised partitioned form. The additively partitioned stiff-problems considered above also fall under this general class. In what follows, we iterate a semi-implicit time discretisation technique to numerically approximate the solution of (3.26)-(3.27) as presented in [9]. The discretisation is done in such a way that the variable $y$ without a stiff dependence is treated explicitly, and the other part implicitly. The formal definition of the scheme is given below, where we consider only a DIRK variant; the interested reader is referred to [9] for other variants and a detailed exposure to the scheme.

Here, the vectors $Y_{I}$ and $Y_{E}$ correspond to the implicit and explicit parts of the intermediate solutions, and the corresponding update formulae are implicit and explicit, respectively.

As a consequence of the above remark, the numerical solution of the SI-IMEX-RK schemes under consideration in this paper boils down to the solution of the last implicit intermediate stage $s$, i.e.,

The order conditions of SI-IMEX schemes and additively partitioned IMEX-RK schemes are more or less the same [9]. Moreover the SI-IMEX schemes are also represented using a double Butcher tableaux, cf. Figure 1.

In order to design RK schemes for the EP system (2.1)-(2.3), henceforth we consider only its DAE formulation (3.13)-(3.14). As a first step, we identify an additive stiff-non-stiff splitting of the function $H$ defined in (3.13). Note that the equation (3.15) already gives a natural splitting of $H$ as the sum of the nonstiff hyperbolic flux $F$ and the stiff source term $S$, cf. the asymptotic analysis presented in Section 2.1. This splitting corresponds to the so-called classical splitting reported in the literature; see [25, 31] for further details. Making use of the indirect approach from Section 3 with the above choice of $F$ and $S$, we define an additively partitioned IMEX-RK-DAE scheme as follows.

The numerical solution is computed in the following chronology. For each $k=1,2,\dots,s$, first, the density $\rho^{k}$ is explicitly updated using the mass conservation update in the system (4.1). Using the updated value of $\rho^{k}$ thus obtained, the linear elliptic problem (4.2) is solved for $\phi^{k}$. The momentum $q^{k}$ can then be obtained using the computed values of $\rho^{k}$ and $\phi^{k}$ from (4.1), which is now an explicit update.

Assuming that the data $\rho^{n},q^{n}$, and $\phi^{n}$ at time $t^{n}$ for the above scheme (4.1)-(4.3) are well-prepared in the sense of Definition 2.1, for each $k=1,2,\ldots,s$, taking the limit $\varepsilon\to 0$ in the Poisson equation (4.2) we obtain

Using the above equation in the mass updates for each $k=1,2,\ldots,s$ we obtain

Therefore, in the limit $\varepsilon\to 0$, the additively partitioned IMEX-RK-DAE scheme (4.1)-(4.3) boils down to the following limiting scheme:

Note that the above limiting scheme doesn’t provide a way for consistently caring forward the divergence-free condition to the numerical solution at time $t^{n+1}$. Precisely speaking, the IMEX-RK-DAE schemes (4.1)-(4.3) only offer the divergence-free condition up-till the $(s-1)^{th}$ stage via (4.8), and fail to provide it for $u^{s}_{(0)}=u^{n+1}_{(0)}$. A closer look at the limiting scheme (4.6)-(4.8) also reveals that, it lacks a way to obtain all the unknowns. Being fettered by the divergence-free condition on the leading order velocity, the IMEX-RK-DAE schemes fail to be consistent with the incompressible limit system (2.5)-(2.7) when $\varepsilon\to 0$ and also donot provide valid recursion to update all the variable at time $t^{n+1}$. Therefore, the IMEX-RK-DAE scheme cannot be AP.

A similar analysis of a first order version of the additively partitioned IMEX-RK-DAE scheme (4.1)-(4.3) can be found in [25]. In [25], the first order scheme is termed as the classical scheme which can be obtained in the present setting by using the IMEX-RK coefficients given in Figure 2. Consequently, the IMEX-RK-DAE schemes can be conceived as natural higher order extensions of the classical scheme.

In [26], based on knowledge from plasma physics, it has been reported that the stability constraints for the classical scheme are extremely prohibitive in the quasineutral regime $\varepsilon\ll 1$. Therefore, combining this observation and the above presented proof of inconsistency, one can conclude that the IMEX-RK-DAE schemes derived from additively partitioned IMEX-RK schemes fail to qualify to being asymptotic preserving.

As a first step towards designing an SI-IMEX-RK scheme based on the indirect approach for the EP system (3.13)-(3.14), we specify the stiff and non-stiff terms. Consequently, it leads to the identification of terms that are to be treated implicitly or explicitly in $H$ and $G$. In the following, tracing the lines of SI-IMEX schemes presentation in Section 3, we subscript variables with $I$ and $E$, for those which are to be treated implicitly and explicitly, respectively. Based on the asymptotic analysis presented in Section 2, we recast the DAE form (3.13)-(3.14) of the EP system via the stiff and non-stiff conservative variables $U_{I}:=(\rho_{I},q_{I})^{T}$ and $U_{E}:=(\rho_{E},q_{E})^{T}$, and the stiff potential $\phi_{I}$ into,

where

Using the indirect approach, we define the SI-IMEX-RK-DAE scheme for the above DAE formulation of the EP system as follows.

Even though the above scheme seems to have more implicit operations to be performed than the one based on additively partitioned approach, we will show that it can be equivalently reformulated so that it has the need to carry out same number of implicit operations as the scheme given in Definition 4.1.

We note that the above proof also serves as a step by step algorithm to sequentially obtain the numerical solution of the reformulated SI-IMEX-RK-DAE (4.16)-(4.22) scheme. We see that, each stage of the SI-IMEX-RK-DAE scheme comprises of three sub-steps: first, an explicit evaluation, second, solution of an elliptic equation, third and last, another explicit evaluation to update the solution of an implicit step. In other words, the numerical solution that we compute in this paper is not the solution of the actual SI-IMEX-RK-DAE scheme (4.12)-(4.15); rather it is the solution of a reformulated Poisson equation (4.20) coupled with the mass (4.22) and momentum (4.21) updates of the original numerical scheme (4.12)-(4.15) which are now explicit in nature.

As observed for the additively partitioned IMEX-RK-DAE scheme here also we see that by making a suitable choice of the RK tableaux; see Figure 3, the first order AP scheme presented in [25, 26] can be obtained as:

That being the case we come to the understanding that the SI-IMEX-RK-DAE scheme (4.12)-(4.15) serves as a natural higher order extension to the first order scheme (4.33)-(4.35).

We refer the reader to [25, 26] and the references therein for a detailed study of the first order scheme; particularly, its computational cost compared to the so-called classical scheme; cf. [25, 26]. As stated there it is imperative to understand that the computational time expense typically associated with any amount of implicitness doesn’t affect our scheme adversely.