version 0.6.0

Topics

 Schemes
 Scheme submodules.
 

Namespaces

module  mod_finite_differences_computer
 Computer for evaluating Finite Difference Schemes.
 
module  mod_finite_differences_3d_stencils
 Stencil filling for 3D schemes.
 
module  mod_finite_differences_init
 Initialization of declared finite differences schemes.
 
module  type_fd_scheme
 The Finite Difference Scheme type definition.
 
module  type_fd_weno_scheme
 The Finite Difference Weno Scheme type definition.
 
module  type_discrete_stencil
 The Finite Difference Scheme type definition.
 

Detailed Description

This module permits to easily build and use finite difference schemes on the Notus grid or any discretized function. They can be used for implicit or explicit use. Various schemes are written for:

Note
All the schemes are written for non uniform grid spacing ; no uniform version is given as the computational cost is almost identical while being numerically the same.

Introduction

The Finite Difference Scheme concept

A finite difference scheme (noted fd_scheme) is defined as the list of coefficients (named weights and noted \(w_i\)) that corresponds to a discrete approximation of a certain derivative. It is based on the mathematical definition of the \(M^{th}\) order approximation of the \(n^{th}\) derivative of a function \(\phi\):

\[ \frac{d^n \phi}{d \xi^n} ( \xi_0 ) \simeq \sum_i w_i \phi( \xi_i ) + O(\Delta \xi^M) = \sum_i w_i \phi_i + O(\Delta \xi^M) \]

where \( \xi \) is any direction (being \(x\), \(y\) or \(z\)), \( \phi_i = \phi(\xi_i) \) are the discrete values of the function \(\phi\). By convention, we note \(\xi_0\) the position to which the finite difference scheme is related.

Code organization

Each scheme is defined in a separate module/file with an explicit name, eg. the scheme first O2 for the first derivative ( \(d \phi / d \xi\)) with second order precision (O2). By convention, the associated functions are named fd_scheme_derivative_order_type, where:

Usage

Schemes can be directly used for implicit formulation or can be applied for explicit formulation (see below).

The Finite Difference Scheme type

A scheme is defined as an array of real values. Each scheme is defined on a specific stencil (a mask on the grid, see the paragraph below) defined as effective elements surrounding a grid element. To each stencil element is associated a value called the weight. These values represent the factor (noted \(w_i\)) associated to each grid element, as presented in the main equation above. As we use a vector representation for weights, null values may appear inside the stencil.

The Fortran type is defined in type_fd_scheme.

Scheme's stencil

The scheme stencil (ie. boundaries) is representend by the min and max indices outside of which the weights are all null. Thus, weights are only stored inside this stencil. The indices are relative to the element where the finite difference will be computed. The stencil's indices are stored in fd_scheme%index_start and fd_scheme%index_end. They are also specified in the documentation of each scheme.

Example
The stencil of the fd_scheme_first_o2_centered scheme is \( [-2 ; 0] \). It implies that the scheme is defined on the grid elements \( ( \phi_{i-2} ; \phi_{i-1} ; \phi_{i} ) \).

Schemes

All the instances of the schemes are (and has to be) created with the procedure fd_initialize()

Every schemes are built (returned) by the fd_scheme* functions.

Example
fd_scheme_first_o2_centered returns the following array: \( \frac{(-1 ; 0 ; 1)}{2h} \) where \(h\) is the constant grid spacing given as input argument
Example
t_fd_scheme_first_o2_backward returns the following array: \( \frac{(1 ; -4 ; 3)}{2h} \) where \(h\) is the constant grid spacing given as input argument

The developer can define her/his own schemes. Please refer to the Schemes documentation for more information.

Implicit usage

In an implicit formulation, one only needs the coefficients associated to each scheme. Thus, the developer will need to declare a variable of type fd_scheme and get the wanted scheme through the build functions detailed above. The coefficients are then accessible in the fd_scheme%weight array. Remember that the weights are indexed relatively to the incriminated discretization point.

Example
type(t_fd_scheme) :: my_scheme
double precision, dimension(2) :: steps = [ 0.1d0 0.12d0 ]
my_scheme = t_fd_scheme_first_o2_backward( steps )
will store in my_scheme%weight(-2), my_scheme%weight(-1) and my_scheme%weight(0) the numerical values of the coefficients that would correspond to \(\phi_{-2}\), \(\phi_{-1}\) and \(\phi_0\), given that the steps between those points are \(0.1\) and \(0.2\).

Explicit usage

An explicit formulation requires that the developer evaluates the finite difference through given values (usualy the values taken from the grid). To do so, there are two ways:

  1. the easy way: simply use the Finite Differences Computer module by supplying the wanted scheme :
double precision, dimension(2) :: steps = [ 0.1d0 0.12d0 ]
double precision, dimension(3) :: values = [ 3.1d0 4.1d0 5.9d0 ]
result = fd_compute( t_fd_scheme_first_o2_backward, steps, values )
  1. the optimized way: if the weights you are using are constant (on a uniform grid for example), you can save computational time by storing the finite difference scheme:
type(t_fd_scheme) :: my_scheme
double precision, dimension(2) :: steps = [ 0.1d0 0.12d0 ]
double precision, dimension(3) :: values = [ 3.1d0 4.1d0 5.9d0 ]
my_scheme = t_fd_scheme_first_o2_backward( steps )
result = fd_compute( my_scheme, values )
result = my_scheme%apply(values)

Detailed example

In this section we furnish a complete Fortran example for numericaly solving the 1D heat equation.

The 1D head equation

Consider the following heat equation:

\[ \frac{\partial \phi}{\partial t} - \mu \frac{\partial^2 \phi}{\partial x^2} = 0 \]

For the formulations below, we assume that the temporal derivative is discretized with a second order backward scheme:

\[ \frac{\partial \phi}{\partial t} \simeq \alpha \phi^{n+1}_i + \beta \phi^{n}_i + \gamma \phi^{n-1}_i \]

with possible variable time step and where superscripts denotes temporal indices and subscripts denotes spatial indices.

Note
one has to note that, in this case, the reference temporal index is \(n+1\) (see the numerical codes for implicit and explicit below).

The following code is used as a preamble for both implicit and explicit formulation, defining the known values:

use mod_finite_differences
type(t_fd_scheme) :: temporal_scheme
type(t_fd_scheme) :: spatial_scheme
integer, parameter :: M = 10 ! Number of spatial points
double precision, dimension(M) :: phi_n, phi_np1, phi_nm1
double precision, dimension(M-1) :: dx ! spatial steps
double precision, dimension(M,M) :: A ! The implicit matrix
double precision, dimension(M) :: b ! The explicit vector (right-hand term)
double precision :: dt_n ! time step from t^n to t^{n+1}
double precision :: dt_nm1 ! time step from t^{n-1} to t^{n}
double precision, parameter :: mu = ... ! Some value
integer :: i, t
phi = ... ! Fill phi
dx = ... ! fill dx

where the wanted result will be stored in phi_np1.

Implicit formulation

The spatial discretization is treated implicitely with a second order centered scheme (assuming that the step is uniform, for simplicity) as:

\[ \frac{\partial^2 \phi}{\partial x^2} \simeq = \frac{\phi^{n+1}_{i-1} - 2 \phi^{n+1}_{i} + \phi^{n+1}_{i+1}}{\Delta x} \]

The complete numerical scheme reads:

\[ \alpha \phi^{n+1}_i - \mu \frac{\phi^{n+1}_{i-1} - 2 \phi^{n+1}_{i} + \phi^{n+1}_{i+1}}{\Delta x} = - \beta \phi^n_i - \gamma \phi^{n-1}_i \]

with the implicit part on the left-hand side of the equation and the explicit part on the right-hand side. The corresponding Fortran code would be:

dt_n = ...
dt_nm1 = ...
do t=2,maxiteration
a = 0.d0 ! Reset the matrix
b = 0.d0 ! Reset the vector
temporal_scheme = t_fd_scheme_first_o2_backward( [ dt_n, dt_nm1 ] )
do i=2,9 ! For all indices but boundaries
spatial_scheme = fd_scheme_second_o2_centered( dx(i-1:i+1) )
a(i,i) = temporal_scheme%weight(0)
a(i,i-1:i+1) = a(i,i-1:i+1) - mu * spatial_scheme%weight(-1:+1)
b(i) = -temporal_scheme%weight(-1)*phi_n(i) - temporal_scheme%weight(-2)*phi_nm1(i)
end do
phi_nm1 = phi_n
phi_n = phi_np1
dt_nm1 = dt_n
dt_n = ...
end do

where the particular treatment needed for the first time step and the spatial boundaries haven't been detailed.

Explicit formulation

The explicit formulation is quite similar to the implicit one presented above. The spatial discretization is treated explicitely with a second order centered scheme (assuming that the step is uniform, for simplicity) as:

\[ \frac{\partial^2 \phi}{\partial x^2} \simeq = \frac{\phi^{n}_{i-1} - 2 \phi^{n}_{i} + \phi^{n}_{i+1}}{\Delta x} \]

The complete numerical scheme reads:

\[ \phi^{n+1}_i = \alpha^{-1} \left( - \beta \phi^n_i - \gamma \phi^{n-1}_i + \mu \frac{\phi^{n}_{i-1} - 2 \phi^{n}_{i} + \phi^{n}_{i+1}}{\Delta x} \right) \]

The corresponding Fortran code would be:

dt_n = ...
dt_nm1 = ...
do t=2,maxiteration
temporal_scheme = t_fd_scheme_first_o2_backward( [ dt_n, dt_nm1 ] )
do i=2,9 ! For all indices but boundaries
phi_np1(i) = -temporal_scheme%weight(-1)*phi_n(i) - temporal_scheme%weight(-2)*phi_nm1(i)
phi_np1(i) = phi_np1(i) + mu * fd_compute( fd_scheme_second_o2_centered, dx(i-1:i+1), phi_n(i-1:i+1) )
phi_np1(i) = phi_np1(i)/temporal_scheme%weight(0)
end do
phi_nm1 = phi_n
phi_n = phi_np1
dt_nm1 = dt_n
dt_n = ...
end do

where the particular treatment needed for the first time step and the spatial boundaries haven't been detailed.