Boundary condition treatment

Namespaces
module	enum_cell_bc_scheme
	Enumeration modules associated to cell scalar boundary condition scheme.
module	enum_face_bc_scheme
	Enumeration modules associated to face-centered vector field boundary condition scheme.

Functions
subroutine	mod_add_cell_bc::add_cell_bc (matrix, rhs, boundary_condition, equation_stencil, equation_ls_map)
	Set boundary condition to a cell-discretized linear system (with ghost cells) More...
subroutine	mod_add_cell_bc_noghost::add_cell_bc_noghost (matrix, rhs, boundary_condition, equation_stencil, equation_ls_map)
	Set boundary condition to a cell-discretized linear system (without ghost cells) More...

Detailed Description

The boundary conditions can be Dirichlet, Neumann and Robin type. Since it is possible to express the first two boundary conditions as special Robin one, we only discretize the Robin BC (coefficient are set in the UI part of the code). It writes:

\begin{align} g = \alpha S_b + \beta \frac {\partial S_b} {\partial n} \end{align}

Imposition of the BC can be done following two methods:

1) the ghost nodes (outside the physical domain) are used to discretized the boundary condition. Lines corresponding to these nodes are added to the linear system.

2) the ghost nodes are not used, the linear system is thus smaller, but matrix coefficients of the inner nodes are modified according to the boundary conditions.

Function Documentation

add_cell_bc()

subroutine mod_add_cell_bc::add_cell_bc	(	double precision, dimension(:), intent(inout)	matrix,
		double precision, dimension(:), intent(inout)	rhs,
		type(t_boundary_condition), intent(in)	boundary_condition,
		type(t_cell_stencil), intent(in)	equation_stencil,
		type(t_ls_map), intent(in)	equation_ls_map
	)

The Robin boundary condition is discretized:

\begin{align} g = \alpha \frac {\partial S_b} {\partial n} + \beta S_b \end{align}

The first row of ghost cells outside the physical domain are used to impose boundary conditions. A second-order discretization scheme is used as follows (right boundary):

--x--x--|··x··  i: first interior cell node,  e: exterior cell node,
  j  i  b  e    j: second interior cell node, b: boundary face node.

• Dirichlet "part": the ghost node value \( S_e \), located at the exterior cell node, is computed using the unknown value from the interior cell node and the given boundary value at the boundary face node as:

  \begin{align} S_b = \tfrac{1}{2}(S_e + S_i) \end{align}

• **Dirichlet (only) **: the given boundary flux at the boundary face is expressed from the unknown values at the interior and exterior cell nodes as:

  \begin{align} S_b = 0.375S_e + 0.750S_i - 0.125S_j \alpha = 0 \beta = 1 \end{align}

• Neumann "part": the given boundary flux at the boundary face is expressed from the unknown values at the interior and exterior cell nodes as:

  \begin{align} \alpha \frac {\partial S_b} {\partial n} = \alpha \frac {(S_e - S_i)} {\Delta x} \end{align}

The second row of ghost cells outside the physical domain may be used by a second order scheme, whose stencil is equal to 9 points (13 points in 3D). For this purpose, a second-order extrapolation is done (to be improved)

add_cell_bc_noghost()

subroutine mod_add_cell_bc_noghost::add_cell_bc_noghost	(	double precision, dimension(:), intent(inout)	matrix,
		double precision, dimension(:), intent(inout)	rhs,
		type(t_boundary_condition), intent(in)	boundary_condition,
		type(t_cell_stencil), intent(in)	equation_stencil,
		type(t_ls_map), intent(in)	equation_ls_map
	)

The Robin boundary condition is discretized:

\begin{align} g = \alpha \frac {\partial S_b} {\partial n} + \beta S_b \end{align}

The boundary condition is discretized with a centered scheme using exterior and interior nodes.

-----|·····  i: interior cell node,  e: exterior cell node,
  i  b  e    b: boundary face node.

The value of the unknwon on the exterior node can be expressed

\begin{align} S_e = A + B S_i \end{align}

with

\begin{align} A = \frac {g} {\alpha/\Delta + \beta/2} \end{align}

and

\begin{align} B = \frac {\alpha/\Delta - \beta/2} {\alpha/\Delta + \beta/2} \end{align}

Then, on the line of the matrix corresponding to the interior node, the matrix coefficient of the element corresponding to the interior node is set to zero, and matrix elements of the interior node are changed thanks to the relation find above. (as well as the RHS)

For instance, in 1D, on the left boundary:

\begin{align} a_0 S_i + a_{left} S_{i-1} + a_{right} S_{i+1} = Rhs_i \end{align}

is changed by:

\begin{align} (a_0 + Ba_{left}) S_i + a_{right} S_{i+1} = Rhs_i - a_{left} A \end{align}

Special care must be given is the stencil size is equal to 2. A second order extrapolation is done:

\begin{align} S_{i2}=2S_i - Se \end{align}

To be improved.