Cell linear system

Linear system routines (construction, modification, finalization) for a scalar equation defined on cells. More...

Namespaces
module	type_ls_map
	Type declaration associated to the numerical grid mapping of a cell linear system.

Functions
subroutine	mod_add_diagonal_cell_matrix::add_diagonal_cell_matrix (matrix, stencil_size, diagonal_term, equation_ls_map)
	Add a term to the diagonal of a matrix for a cell-discretized equation.
subroutine	mod_add_line_cell_matrix::add_line_cell_matrix (matrix, matrix_stencil, equation_stencil, equation_ls_map, ns, i, j, k)
	Add a local stencil to a cell matrix, according to different stencil types.
subroutine	mod_add_rhs_cell::add_rhs_cell (rhs, grid_rhs, equation_ls_map)
	Add a term to the right-hand side for a cell-discretized equation.
subroutine	mod_finalize_cell_linear_system::finalize_cell_linear_system (matrix, rhs, stencil_size, equation_ls_map)
	Finalize a cell linear system.
double precision function	mod_impose_cell_value::get_nearest_cell_field_value (field, position)
	Get the field nearest cell's value at the given position.
subroutine	mod_impose_cell_value::impose_cell_value (matrix, rhs, position, value, stencil_size, equation_ls_map)
	Impose the linear system to yield the given value at the given position.
subroutine	mod_initialize_cell_linear_system::initialize_cell_linear_system (matrix, rhs, ls_map)
	Allocate and nullify matrix and right-hand-side of a cell linear system.
subroutine	mod_time_advance_explicitly::time_advance_explicitly (matrix, rhs, array, equation_ls_map, stencil_size, boundary_condition)
	Determine \( \phi^{n+1} \) for a fully explicit cell transport implementation.

Detailed Description

Linear system routines (construction, modification, finalization) for a scalar equation defined on cells.

Principle

The set of linear equations is translated into a linear system that is represented by a matrix \(A\) and a vector \(b\), usually referred as right-hand-side (rhs) such that:

\[ A \hat{\phi} = b. \]

where \(\hat{\phi}\) is a vectorized/serialized representation of the unknowns \(\phi_{i,j,k}\) on the grid. Linear solvers (Linear system solvers) are then used to solve

\[ \min_{\hat{\phi}}\,|A\hat{\phi}-b|. \]

Matrix \(A\) is mostly sparse, based on the stencil (Stencil type) used. For example, the classic second order star 2D stencil represented as:

 ━━━━━━╋━━━━━━━╋━━━━━━━
       ┃       ┃
       ┃  φ_T  ┃
       ┃       ┃
 ━━━━━━╋━━━━━━━╋━━━━━━━
┃      ┃       ┃       ┃
┃  φ_L ┃  φ_C  ┃  φ_R  ┃
┃      ┃       ┃       ┃
 ━━━━━━╋━━━━━━━╋━━━━━━━
       ┃       ┃
       ┃  φ_B  ┃
       ┃       ┃
 ━━━━━━╋━━━━━━━╋━━━━━━━

where \(C,L,R,B,T\) stand for, respectively, center, left, right, bottom, bottom and top. Each stencil associated coefficient is stored in a compact matrix form, such as all \(M_{i,j}\) lines:

\[ \left(\begin{array}{ccccc} & & \ldots\\ M_{i,j}^{C} & M_{i,j}^{L} & M_{i,j}^{R} & M_{i,j}^{B} & M_{i,j}^{T}\\ & & \ldots \end{array}\right). \]

The correspondance between grid, matrix and vector numbering is done via the t_ls_map (type_ls_map) and the t_cell_stencil types (type_stencil).

Strategy

In order to solve the cell scalar equation (Cell scalar), each term is discretized at each cell. To do so, we use a local stencil represented by a \(N^3\) matrix, where \(N\) is the stencil size ( \(N=3\) for the previous example). Once the cell neighbors' coefficients have been computed, they can be transferred into the matrix. The corresponding pseudo-algorithm reads:

M <- 0 # the linear system matrix to compute
For all terms in the equation do:
    For all cells i,j do:
        1 - Compute m_ij, the stencil matrix storing all term's coefficients
        2 - Transfer m_ij into M
Apply boundary conditions to M and rhs
Solve the linear system for phi

Code

The module provides the t_ls_map type (type_ls_map) that is passed to discretization routines for handling the matrix and rhs elements.

The module also provides the following routines to help building the linear system:

• initialize_cell_linear_system has to be called to allocate the matrix and the rhs.
• finalize_cell_linear_system has to be called to finalize the build-up, prior to calling the solver.
• add_diagonal_cell_matrix adds array elements to the diagonal terms of the matrix.
• add_rhs_cell adds array elements to the right-hand side.
• add_line_cell_matrix adds a stencil elements to a line of the matrix.
• impose_cell_value sets the linear system to impose a specific value at a specific position.
• time_advance_explicitly given the matrix and rhs, explicitly return the solution.

Routine discretize_cell_transport_equation can be used to assemble all terms of the equation (see solver_energy.f90 for instance)

Function Documentation

add_diagonal_cell_matrix()

subroutine mod_add_diagonal_cell_matrix::add_diagonal_cell_matrix	(	double precision, dimension(:), intent(inout)	matrix,
		integer, intent(in)	stencil_size,
		double precision, dimension(:,:,:), intent(in)	diagonal_term,
		type(t_ls_map), intent(in)	equation_ls_map )

Add a term to the diagonal of a matrix for a cell-discretized equation.

Parameters

[in,out]	matrix	the linear system's matrix
[in]	stencil_size	the stencil size
[in]	diagonal_term	the diagonal term to add
[in]	equation_ls_map	type_ls_map::t_ls_map

This routine adds a diagonal term such that:

\[ m_{ai} \leftarrow m_{ia} + d(i,j,k) \textrm{, with } ia \textrm{ the matrix index of the i,j,k cell.} \]

add_line_cell_matrix()

subroutine mod_add_line_cell_matrix::add_line_cell_matrix	(	double precision, dimension(:), intent(inout)	matrix,
		double precision, dimension(-ns:ns,-ns:ns,-ns:ns), intent(in)	matrix_stencil,
		type(t_cell_stencil), intent(in)	equation_stencil,
		type(t_ls_map), intent(in)	equation_ls_map,
		integer, intent(in)	ns,
		integer, intent(in)	i,
		integer, intent(in)	j,
		integer, intent(in)	k )

Add a local stencil to a cell matrix, according to different stencil types.

Parameters

[in,out]	matrix	the linear system's matrix
[in]	matrix_stencil	the local stencil [-ns:ns]^3
[in]	equation_stencil	See type_stencil::t_cell_stencil
[in]	equation_ls_map	See type_ls_map::t_ls_map
[in]	i,j,k	the grid position of the cell
[in]	ns	the local stencil size

Based on the equation's stencil type, add the [ns x ns]^3 local stencil to the linear system's matrix elements.

Warning

no check is done if some elements are outside the equation's stencil limits.

add_rhs_cell()

subroutine mod_add_rhs_cell::add_rhs_cell	(	double precision, dimension(:), intent(inout)	rhs,
		double precision, dimension(:,:,:), intent(in)	grid_rhs,
		type(t_ls_map), intent(in)	equation_ls_map )

Add a term to the right-hand side for a cell-discretized equation.

Parameters

[in,out]	rhs	the linear system's rhs
[in]	grid_rhs	the grid based rhs
[in]	equation_ls_map	type_ls_map::t_ls_map

Simply adds the cell grid array to the rhs.

finalize_cell_linear_system()

subroutine mod_finalize_cell_linear_system::finalize_cell_linear_system	(	double precision, dimension(:), intent(inout)	matrix,
		double precision, dimension(:), intent(inout)	rhs,
		integer, intent(in)	stencil_size,
		type(t_ls_map), intent(in)	equation_ls_map )

Finalize a cell linear system.

Parameters

[in,out]	matrix	the linear system's matrix
[in,out]	rhs	the linear system's rhs
[in]	stencil_size	the matrix stencil size
[in]	equation_ls_map	See type_ls_map::t_ls_map

The routine checks the if a diagonal element of the matrix is null (typically a corner of the domain). If so, impose a null value for the associated unknown by setting the diagonal to 1 and the rhs to 0.

Todo

improve unused nodes values

get_nearest_cell_field_value()

double precision function mod_impose_cell_value::get_nearest_cell_field_value	(	double precision, dimension(:,:,:), intent(in)	field,
		double precision, dimension(3), intent(in)	position )

Get the field nearest cell's value at the given position.

Parameters

[in]	field	the scalar field
[in]	position	the given position

Returns

the value of the cell nearest to position. If not in the current processor, return -hude(1d0)

impose_cell_value()

subroutine mod_impose_cell_value::impose_cell_value	(	double precision, dimension(:), intent(inout)	matrix,
		double precision, dimension(:), intent(inout)	rhs,
		double precision, dimension(3), intent(in)	position,
		double precision, intent(in)	value,
		integer, intent(in)	stencil_size,
		type(t_ls_map), intent(in)	equation_ls_map )

Impose the linear system to yield the given value at the given position.

Parameters

[in,out]	matrix	the linear system's matrix
[in,out]	rhs	the linear system's rhs
[in]	position	the 3D point where to impose the value
[in]	value	the value to impose
[in]	stencil_size	the matrix stencil size
[in]	equation_ls_map	See type_ls_map::t_ls_map

This routine modifies the linear system such as the solution will yield the given value at the given position. Pretty brutal, this routine has been implemented to fix ill-posed linear system such as the Laplace equation with Neumann boundary conditions everywhere.

Note

the value is imposed at the nearest cell to the given position.

See also

get_nearest_cell_indices

initialize_cell_linear_system()

subroutine mod_initialize_cell_linear_system::initialize_cell_linear_system	(	double precision, dimension(:), intent(out), allocatable	matrix,
		double precision, dimension(:), intent(out), allocatable	rhs,
		type(t_ls_map), intent(in)	ls_map )

Allocate and nullify matrix and right-hand-side of a cell linear system.

Parameters

[out]	matrix	the linear system's matrix
[out]	rhs	the linear system's rhs
[in]	ls_map	See type_ls_map::t_ls_map

time_advance_explicitly()

subroutine mod_time_advance_explicitly::time_advance_explicitly	(	double precision, dimension(:), intent(in)	matrix,
		double precision, dimension(:), intent(in)	rhs,
		double precision, dimension(nx,ny,nz), intent(out)	array,
		type(t_ls_map), intent(in)	equation_ls_map,
		integer, intent(in)	stencil_size,
		type(t_boundary_condition), intent(in)	boundary_condition )

Determine \( \phi^{n+1} \) for a fully explicit cell transport implementation.

Parameters

[in]	matrix	the linear system's matrix
[in]	rhs	the linear system's rhs
[out]	array	the resulting array
[in]	equation_ls_map	See type_ls_map::t_ls_map
[in]	stencil_size	the stencil size
[in]	boundary_condition	See type_boundary_condition::t_boundary_condition

The time_advance_explicitly routine solves \(\alpha_{i,j,k} \phi^{n+1} = RHS_{i,j,k} \) and applies the appropriate boundary conditions.