version 0.6.0
Discretization of PDE

Topics

 Explicit discretization
 Routines for discretizing explicit terms.
 
 Field level operators
 High-level functions for interpolating/extrapolating/differenciating whole fields.
 
 Immersed Boundaries
 Implement boundaries of arbitrary shapes in PDEs.
 
 Implicit discretization
 Routines that fill the matrix and the right-hand side of linear systems associated to the discretization of PDE.
 
 Lagrangian particles
 Lagrangian particles transported by the flow.
 
 Multiphase methods
 Interface advection equation resolution methods (Volume of Fluid, Moment of Fluid, Level Set, etc.)
 
 Node level schemes
 Node level (cell, face, edge, vertex, point) discrete schemes.
 
 Partial Differential Equations
 First level routines to discretize PDE term by term.
 
 Time discretization
 Routines relative to time discretization of PDEs.
 

Namespaces

module  enum_discretization_type
 Discretization types identifiers.
 

Detailed Description

This directory contains the lower level of Notus numerical methods:

General approach

Let's consider the general form of a Partial Differential Equation

\begin{align} F_1(\phi, \partial_x \phi, \ldots, t) + F_2(\phi, \partial_x \phi, \ldots, t) + \cdots = 0\,, \end{align}

where \( \phi \) is the unknown, and functions \( F_1, F_2, \ldots \) represents terms of the PDE (advection term, diffusion term, etc.) The discretization consists to transform the PDEs into a linear system

\begin{align} A \phi = b \end{align}

where \( A \) and \( b \) are the matrix and the right-hand side of the linear system, respectively. In the process, the derivatives of \( \phi \) are transformed in some combinations of values of \( \phi \). The filling of the linear system is done term-by-term

\begin{align} A &= A_1 + A_2 + \cdots\,, & b &= b_1 + b_2 + \cdots\,. \end{align}

Relatively to the temporal discretization (see the IMEX paragraph below), the various terms including spatial derivatives falls in one of these two categories:

The directories Explicit discretization and Implicit discretization contains routines add_<some-name>_term that fills the linear system.

When the PDE contains time derivatives, the unknown is generally \( \phi^{n+1} \) and the linear system usually contains \( \phi^{n} \) in its definition. The time derivative term is discretized differently, and the corresponding routines are in Time discretization.

Additionally, there are other directories that covers specific topics:

of the equations.

In Notus, most PDEs are discretized implicitly except for the phase advection equations and optionally the advection term of the others equations.

Mixted IMplicit-EXplicit schemes (IMEX)

Theory

Mixing implicit and explicit terms in the solving of a PDE is not trivial. There exists a family of methods known under the name of IMEX (for IMplicit-Explicit). Various schemes exist (eg. implicit diffusion coupled with explicit advection) for various order of precision. The choice of one variant or another is usualy based on the efficiency/accuracy ratio and stability matters that all highly depend on the considered problem. For example, stiff problems (ie. for us, when diffusion is higher than advection) might require a second order (ie. Crank_Nicolson) in time discretization of the diffusion term. On the other hand, stiff problems generally have higher restriction on the time step.

The general form of a PDE can be written:

\[ \frac{\partial \phi}{\partial t} = F(\phi) + \frac{1}{\epsilon} R(\phi) + f(t) \]

where \( \epsilon > 0\) (and where \( \epsilon \rightarrow 0 \) denotes stiff problems) and \(f(t)\) is independant of \(\phi\). We can rewrite the previous equation by seperating the operators depending on the explicit (using known values of \(\phi\)) and implicit (using unknown values of \(\phi\)) treatment as:

\[ \frac{\partial \phi}{\partial t} = I(\phi) + E(\phi) + f(t) \]

where \(I(\phi)\) (resp. \(E(\phi)\) stands for the implicit (resp. explicit) terms.

Consider the advection-diffusion equation as an example:

\[ \frac{\partial \phi}{\partial t} = \nu \Delta \phi - c\frac{\partial\phi}{\partial x} \]

It is quite common to use an explicit scheme for the advection term. These schemes are particularly efficient when advection is dominant and/or when high-order spatial schemes are necessary for computing \( \nabla\phi \). However, numerical stability issues can induce very low time steps (see CFL) which can lead to increasing computational time. On the other side, the diffusion term can be efficiently treated with an implicit scheme (here, explicit schemes suffer for even more restrictive time steps without particular gain in accuracy). Hence, it is natural to temporaly discretize the equation as follows (where we seek to compute \(\phi^{n+1}\)):

\[ \frac{\phi^{n+1} - \phi^{n}}{\delta t} = \nu \Delta \phi^{n+1} - c^n\frac{\partial\phi^n}{\partial x} \]

This scheme is clearly first order in time and is submitted to the CFL condition because of the explicit treatment of the advection term. It is called SBDF1 (for Semi-implicit Backward Dfiferentiation Formula). As stated above, we can imagine a wide variety of schemes depending on the combinasions un/known values of \(\phi\) at various times. For a detailed presentation of different schemes, we refer the reader to [Rosales2018].

Remember\n
There are two main families of methods: Multi-Steps (MS) methods and Runge-Kutta (RK) methods. The first considers previous (backward) time values (memory hungry) while the second ones consider forward in time discretization of the integrals (computational time hungry).

Implementation in Notus

In Notus, the user has the choice on the various discretization methods: fully implicit, 1st order IMEX or 2nd order MS-IMEX. A particularity resides in the fact that high-order Runge-Kutta integration methods can also be used for the explicit advection term (be aware that, for now, no implicit Runge-Kutta scheme has been implemented). The classic second order MS-IMEX method with an explicit scheme for the advection, namely the SBDF2 scheme, is written:

\[ \frac{\alpha \phi^{n+1} + \beta \phi^{n} + \gamma \phi^{n-1}}{\delta t} = \nu \Delta \phi^{n+1} - \widetilde{Ext}(c\frac{\partial\phi}{\partial x},t^{n+1}) \]

(see compute_time_coefficients for the definition of the associated coefficients of the temporal derivative). Here the term \(\widetilde{Ext}(c\frac{\partial\phi}{\partial x},t^{n+1})\) stands for the second order extrapolation (of the known computed values) of the advection term towards time \(t^{n+1}\); this is a particularity of the SBDF2 scheme. If one was using a (first order) Euler approach for the explicit term (ie. computing the approximation of the advection term at point values \(t^{n-1}\) and \(t^{n}\)), the classical extrapolation scheme would be:

\[ \widetilde{Ext}(c\frac{\partial\phi}{\partial x},t^{n+1}) = 2 c^n\frac{\partial\phi^n}{\partial x} - c^{n-1}\frac{\partial\phi^{n-1}}{\partial x} \]

However, as we want to be able to use high-order Runge-Kutta explicit integration of the advection term, the interpolation scheme is a bit different. Actualy, a high-order integration of the advection scheme will give us the mean value of the advection term in the temporal interval such that:

\[ \overline{E}(t^n,t^{n+1}) = \frac{1}{\delta t} \int_{t^n}^{t^{n+1}} \frac{\partial \phi}{\partial t} = \frac{\phi^*(t^{n+1}) - \phi(t^{n})}{\delta t} = \frac{1}{\delta t} \int_{t^n}^{t^{n+1}} -c\frac{\partial\phi}{\partial x} \]

where \(\phi^*(t^{n+1})\) is the numerical result of the Runge-Kutta integration of the advection scheme towards time \(t^{n+1}\). Hence, we are in fact dealing with mean values of the explicit term and not point values. Thus, the reconstruction of the extrapolated explicit term at time \(t^{n+1}\) is:

\[ \widetilde{Ext}(c\frac{\partial\phi}{\partial x},t^{n+1}) = \frac{3}{2} \overline{E}(t^{n-1},t^{n}) - \frac{1}{2} \overline{E}(t^n,t^{n+1}) \]

Obviously, the extrapolation is done correctly (with different weights) when the time step is variable.

Note
For the first iteration, we are explicitely reducing the time order discretization to 1 as the function at previous time step is not known. A Runge-Kutta approach would in that case be interesting for an overall second order accuracy.
Keep in mind
When adding an external source term to your equation, be aware that as to conserve the accuracy of the second order temporal scheme, one has to consider an appropriate integration of this term and how to correctly integrate it in the associated SBDF2 method.
Validation
The validation of the various advection schemes are provided in advection_schemes.

Discussions

While a Crank-Nicolson 2nd order scheme for the diffusion term might seems to be an interesting choice when using a first order discretization of the temporal derivative, it appears to give doubtly pertinent physical results when using large time-steps (citation needed).

Choosing the schemes in Notus UI

The temporal accuracy can be changed through the time_order_discretization keyword in the numerical_parameters block. Implicit/explicit schemes for the advection terms are also selectable for each equation through the advection_scheme keyword.

Note
We refer the reader to the associated NTS files documentation for a detailed list of .

References

[Rosales2018]: Unconditional Stability for Multistep ImEx Schemes – Theory. Rodolfo R. Rosales, Benjamin Seibold, David Shirokoff, Dong Zhou. February 19, 2018