Linear system solvers

Linear system solvers and interface to third-party solvers. More...

Topics
	HYPRE linear system interface
	Massively parallel HYPRE library interface routines.
	LIS linear system interface
	Parallel LIS library interface routines.
	MUMPS linear solver library
	Interface for the MUMPS linear solver library.

Classes
type	type_solver::t_linear_solver
	Parent type for all solvers. More...
type	type_solver::t_direct_solver
	Parent type for all direct solvers. More...
type	type_solver::t_iterative_solver
	Parent type for all iterative solvers. More...

Functions
subroutine	mod_linear_system_solver_cell_interface::linear_system_solver_cell_interface (matrix, sol, rhs, solver, boundary_condition, stencil, ls_map, equation_name, coupled_equations, sol2, boundary_condition2, is_print_residual)
	Interface routine for cell solvers.
subroutine	mod_linear_system_solver_face_interface::linear_system_solver_face_interface (matrix, sol, rhs, solver, boundary_condition, face_stencil, face_ls_map, cell_ls_map)
	Fortran interface routine for face solvers (HYPRE SSTRUCT interface, Lis, Notus and MUMPS direct solvers).
subroutine	mod_scaled_residual::compute_scaled_residual_cell (matrix, rhs, solution, ls_map, stencil, scaled_residual)
	Compute the scaled residual of a cell-based linear system.
subroutine	mod_scaled_residual::compute_scaled_residual_face (matrix, rhs, solution, ls_map, stencil, scaled_residual)
	Compute the scaled residual of a face-based linear system.

Detailed Description

This directory provides a consistent interface to linear system solvers, which can be either built-in or third-party backend. A derived type system selects the solver to use, carry solver-specific parameters and carry result information. The routines linear_system_solver_cell_interface and linear_system_solver_face_interface solves a linear system by dispatching to the appropriate solver.

Currently, the built-in backends are:

• Notus BiCGStab solver.

and the third-party backends are:

• MUMPS (MUltifrontal Massively Parallel sparse direct Solver),
• HYPRE (Scalable Linear Solvers and Multigrid Methods),
• LIS (Library of Iterative Solvers for linear systems).

The table below presents available solvers for each interface

Backend	linear_system_solver_cell_interface	linear_system_solver_face_interface
Hypre_Struct	Yes	No
Hypre_SStruct	Yes	Yes
Mumps	Yes	Yes
Notus BiCGStab	Yes	Yes
Lis	Yes	Yes

Common interface for solvers

Notus uses inheritance of derived types to determine which solver to use. The two fundamental types are the type_solver::direct_solver type and the type_solver::iterative_solver type, which both are derived from the type_solver::base_solver type.

base_solver
 │
 ├── direct_solver
 │
 └── iterative_solver

Every type_solver::direct_solver has a residual attribute, which receives the solver's result residual.

Every type_solver::iterative_solver has max_iter and tolerance parameters, which holds the maximum number of iterations and the residual tolerance desired. The derived type also has a final_iter and a residual argument, which receives the actual number of iteration and the residual after the last iteration.

Notus' Jacobi preconditioner

Notus has its own preconditioner that can be applied prior to any solver backend. This is a left Jacobi preconditioner that simply divides each row by the diagonal coefficient.

Every solver has a is_jacobi_initial_preconditioner attribute to activate or deactivate the Notus' Jacobi preconditioner.

Function Documentation

compute_scaled_residual_cell()

subroutine mod_scaled_residual::compute_scaled_residual_cell	(	double precision, dimension(:), intent(in)	matrix,
		double precision, dimension(:), intent(in)	rhs,
		double precision, dimension(:,:,:), intent(in)	solution,
		type(t_ls_map), intent(in)	ls_map,
		type(t_cell_stencil), intent(in)	stencil,
		double precision, intent(out)	scaled_residual )

The scaled residual is equal to

\[ \frac{|Ax - b|_\infty}{|A|_\infty|x|_\infty} \]

where the infinity norm of the matrix \( A \) is equal to

\[ \max\left{\sum_{j=1}^n |A_{ij}| \mid| i \in [1,n] \right} \]

Parameters

[in]	matrix	matrix of the linear system.
[in]	rhs	right-hand side of the linear system (serialized).
[in]	solution	solution of the linear system in (i,j,k) form.
[in]	ls_map	linear system mapping structure.
[in]	stencil	stencil of the discretization.
[out]	scaled_redidual	scaled residual.

compute_scaled_residual_face()

subroutine mod_scaled_residual::compute_scaled_residual_face	(	double precision, dimension(:), intent(in)	matrix,
		double precision, dimension(:), intent(in)	rhs,
		type(t_face_field), intent(in)	solution,
		type(t_face_ls_map), intent(in)	ls_map,
		type(t_face_stencil), intent(in)	stencil,
		double precision, intent(out)	scaled_residual )

The scaled residual is equal to

\[ \frac{|Ax - b|_\infty}{|A|_\infty|x|_\infty} \]

where the infinity norm of the matrix \( A \) is equal to

\[ \max\left{\sum_{j=1}^n |A_{ij}| \mid| i \in [1,n] \right} \]

Parameters

[in]	matrix	matrix of the linear system.
[in]	rhs	right-hand side of the linear system (serialized).
[in]	solution	solution of the linear system in (i,j,k) form.
[in]	ls_map	linear system mapping structure.
[in]	stencil	stencil of the discretization.
[out]	scaled_redidual	scaled residual.

linear_system_solver_cell_interface()

subroutine mod_linear_system_solver_cell_interface::linear_system_solver_cell_interface	(	double precision, dimension(:), intent(inout)	matrix,
		double precision, dimension(:,:,:), intent(inout)	sol,
		double precision, dimension(:), intent(inout)	rhs,
		class(t_linear_solver), intent(inout)	solver,
		type(t_boundary_condition), intent(in)	boundary_condition,
		type(t_cell_stencil), intent(in)	stencil,
		type(t_ls_map), intent(in)	ls_map,
		character(len=*), intent(in)	equation_name,
		logical, intent(in), optional	coupled_equations,
		double precision, dimension(:,:,:), intent(inout), optional	sol2,
		type(t_boundary_condition), intent(in), optional	boundary_condition2,
		logical, intent(in), optional	is_print_residual )

This routine accepts the following classes of solvers:

• t_hypre_solver_Struct_interface
• t_hypre_solver_SStruct_interface
• lis_solver
• t_mumps_solver
• notus_solver_mg

The following steps are done:

• sol (array of dimension 3) is put in one dimension array solution
• Point Jacobi preconditioner is applied (left preconditioning)
• solvers are called
• solution is put in sol
• boundary condition are applied on boundaries ghost cells
• exchange beetween processors is done

Parameters

[in,out]	matrix	matrix of the linear system (vectorized)
[in,out]	sol	solution of the linear system (not vectorized)
[in,out]	rhs	right-hand-side of the linear system (vectorized)
[in,out]	solver	solver structure associated to the equation
[in]	boundary_condition	boundary condition
[in]	stencil	stencil structure of the equation
[in]	ls_map	mapping associated to the cell grid
[in]	equation_name	name of the equation (log message)
[in]	coupled_equations	optional logical to solve a 2-fields coupled linear system
[in,out]	sol2	optional solution of the second unknown
[in]	boundary_condition2	(optional) boundary condition associated to the second unknown
[in]	is_print_residual	(optional) print the residual

Note

The sol argument does not need to be vectorized at upper level, so it should be provided as a rank-3 array, and the vectorization occurs in the routine. On the opposite, rhs array is vectorized at a upper level.

linear_system_solver_face_interface()

subroutine mod_linear_system_solver_face_interface::linear_system_solver_face_interface	(	double precision, dimension(:), intent(inout)	matrix,
		type(t_face_field), intent(inout)	sol,
		double precision, dimension(:), intent(inout)	rhs,
		class(t_linear_solver), intent(inout)	solver,
		type(t_boundary_condition_face), intent(in)	boundary_condition,
		type(t_face_stencil), intent(in)	face_stencil,
		type(t_face_ls_map), intent(in)	face_ls_map,
		type(t_ls_map), intent(in)	cell_ls_map )

The following steps are performed:

• The rank 3 solution 'sol' (array of dimension 3) is copied into a rank 1 array.
• Jacobi preconditioner is applied (left preconditioning).
• Linear solvers are called.
• The rank 1 solution is copied into the 'sol' rank 3 array.
• Boundary conditions are applied on ghost boundary faces.
• The solution is exchanged between processes.

Parameters

[in,out]	matrix	matrix of the linear system (vectorized)
[in,out]	sol	solution of the linear system (not vectorized)
[in,out]	rhs	right-hand side of the linear system (vectorized)
[in,out]	solver	linear solver
[in]	boundary_condition	boundary conditions on faces
[in]	face_stencil	stencil structure
[in]	face_ls_map	mapping associated to the face grid for HYPRE solver
[in]	cell_ls_map	mapping associated to cell grid

Note: solution do not need to be vectorized at a upper level. This is the reason why the corresponding argument is a rank 3 array, that is vectorized inside this routine. On the opposite, rhs array is vectorized at a upper level.