Implicit diffusion schemes of a scalar equation defined on cells. More...

Functions
subroutine	mod_harmonic_diffusion_coef_interpolation::harmonic_diffusion_coef_interpolation (face_diffusion_coef, cell_diffusion_coef)
	Harmonic interpolatation of the diffusion coefficient to face nodes.

subroutine	mod_add_cell_diffusion_term_centered_o2::add_cell_diffusion_term_centered_o2 (matrix, flux_coef, equation_stencil, equation_ls_map)
	Add the diffusive term of an equation defined on cells to a matrix.

subroutine	mod_apply_cell_diffusion_scheme::apply_cell_diffusion_o2_centered_scheme (face_diffusion_coefficient, matrix_line, ns, i, j, k, extra_inside_gradient_coef2)
	Discretize the diffusive term of an equation defined on cells.

subroutine	mod_interpolate_diffusion_coef::interpolate_diffusion_coef (face_diffusion_coef, cell_diffusion_coef)
	Interpolate diffusion coefficient to face nodes.

subroutine	mod_set_diffusive_flux_coef::set_diffusive_flux_coef (flux_coef)
	Convert face diffusion coefficient to face diffusive flux coefficient in place.

subroutine	mod_thermal_resistance::add_face_resistance (face_conductivity, face_resistance)
	Adds face diffusive resistance to face conductivity.

Detailed Description

Discretization of the diffusion term

In Notus, the diffusion terms \( \nabla \cdot ( D \nabla \phi) \) is discretized using a finite volume approach.

       ┃ φ_0+1 ┃
       ┃       ┃
 ━━━━━━╋━━━∧━━━╋━━━━━━  ┬
       ┃       ┃        │
 φ_-10 > φ_00  > φ_+10  │ δy_0
       ┃       ┃        │
 ━━━━━━╋━━━∧━━━╋━━━━━━  ┴
       ┃       ┃
       ┃ φ_0-1 ┃

       ├───────┤
          δx_0

Here, \( D \) is the diffusion coefficient, and \( \phi \) is the implicit variable. The diffusion term at the cell center is approximated as the average of the diffusion term over the cell \( (0,0) \), and the Green-Ostrogradsky theorem allow us to write the integral as

\begin{align} \nabla \cdot ( D \nabla \phi)_{0,0} = \frac{1}{\delta x_0 \delta y_0} \iint_{\text{cell}} \nabla \cdot ( D \nabla \phi) &= \frac{1}{\delta x_0 \delta y_0} \left( \int_{\text{left}} D \nabla \phi \cdot \mathbf{n} + \int_{\text{right}} D \nabla \phi \cdot \mathbf{n} + \int_{\text{bottom}} D \nabla \phi \cdot \mathbf{n} + \int_{\text{top}} D \nabla \phi \cdot \mathbf{n} \right) \,, \end{align}

and using face values to approximates the fluxes gives

\begin{align} \nabla \cdot ( D \nabla \phi)_{0,0} = \frac{1}{\delta x_0} (D \nabla \phi)_{+\frac12,0} -\frac{1}{\delta x_0} (D \nabla \phi)_{-\frac12,0} +\frac{1}{\delta y_0} (D \nabla \phi)_{0,+\frac12} -\frac{1}{\delta y_0} (D \nabla \phi)_{0,-\frac12} \,, \end{align}

where each term \( (D \nabla \phi)_{+\frac12,0} \) represent a face diffusive flux density, representable by a face_field. Gradients are discretized using a seconth-order centered method

\begin{align} \nabla \cdot ( D \nabla \phi)_{0,0} =& \frac{1}{\delta x_0} \frac{D_{+\frac12,0}}{\delta x_{+\frac12}} \phi_{+1,0} -\frac{1}{\delta x_0} \left(\frac{D_{+\frac12,0}}{\delta x_{+\frac12}} + \frac{D_{-\frac12,0}}{\delta x_{-\frac12}}\right) \phi_{0,0} +\frac{1}{\delta x_0} \frac{D_{-\frac12,0}}{\delta x_{-\frac12}} \phi_{-1,0} \\ +& \frac{1}{\delta y_0} \frac{D_{0,+\frac12}}{\delta y_{+\frac12}} \phi_{0,+1} -\frac{1}{\delta y_0} \left(\frac{D_{0,+\frac12}}{\delta y_{+\frac12}} + \frac{D_{0,-\frac12}}{\delta y_{-\frac12}}\right) \phi_{0,0} +\frac{1}{\delta y_0} \frac{D_{0,-\frac12}}{\delta y_{-\frac12}} \phi_{0,-1} \,. \end{align}

where each term \( \frac{D_{+\frac12,0}}{\delta x_{+\frac12}} = F_{+\frac12,0} \) represent a face diffusive flux coefficient (flux_coef). The routine add_cell_diffusion_term_centered_o2 adds a diffusion term to its matrix argument for any given face diffusive flux coefficient. The routine discretize_cell_diffusion_term adds a diffusion term ot its matrix argument for any given diffusion coefficient.

Computation of face diffusive flux coefficients

The face diffusive flux coefficient \( F_{+\frac12,0} \) is computed using steady solutions of the diffusion equation.

 ┏━━━┯━━━┳━━━┯━━━┓  ┬
 ┃   │   ┃   │   ┃  │
 ┃  φ_0  >  φ_1  ┃  │ δy_0
 ┃   │   ┃   │   ┃  │
 ┗━━━┷━━━┻━━━┷━━━┛  ┴

 ├───────┼───────┤
    δx_0   δx_1

The diffusive flux density is constant through the two cell centers

\begin{align} (D \nabla \phi)_\frac12 &=& F_\frac12 &(\phi_1 - \phi_0) & &=& F_0 &(\phi_\frac12 - \phi_0) & &=& F_1 &(\phi_1 - \phi_\frac12) \\ &=& \frac{D_\frac12}{\delta x_\frac12} &(\phi_1 - \phi_0) & &=& \frac{D_0}{\tfrac12\delta x_0} &(\phi_\frac12 - \phi_0) & &=& \frac{D_1}{\tfrac12\delta x_1} &(\phi_1 - \phi_\frac12)\,. \end{align}

Hence the face diffusion coefficient \( D_\frac12 \) is obtained by harmonic interpolation of the cell diffusion coefficients

\begin{align} D_\frac12 = \dfrac{ \tfrac12\delta x_0 + \tfrac12\delta x_1 }{ \dfrac{\tfrac12\delta x_0}{D_0} + \dfrac{\tfrac12\delta x_1}{D_1} }\,, \end{align}

and the face diffusive flux coefficient by

\begin{align} F_\frac12 = \dfrac{ D_\frac12 }{ \delta x_\frac12 } \,. \end{align}

The routine interpolate_diffusion_coef interpolates the diffusion coefficient from cell to face, while the routine set_diffusive_flux_coef converts face diffusion coefficient to face diffusive flux coefficient.

Function Documentation

◆ add_cell_diffusion_term_centered_o2()

subroutine mod_add_cell_diffusion_term_centered_o2::add_cell_diffusion_term_centered_o2	(	double precision, dimension(:), intent(inout)	matrix,
		type(t_face_field), intent(in)	flux_coef,
		type(t_cell_stencil), intent(in)	equation_stencil,
		type(t_ls_map), intent(in)	equation_ls_map )

The term writes:

\begin{align} \nabla \cdot ( D \nabla \phi)_{0,0} =& \frac{1}{\delta x_0} F_{+\frac12} \phi_{+\frac12,0} -\frac{1}{\delta x_0} \left(F_{+\frac12} + F_{-\frac12}\right) \phi_{0,0} +\frac{1}{\delta x_0} F_{-\frac12} \phi_{-\frac12,0} \\ +& \frac{1}{\delta y_0} G_{+\frac12} \phi_{0,+\frac12} -\frac{1}{\delta y_0} \left(G_{+\frac12} + G_{-\frac12}\right) \phi_{0,0} +\frac{1}{\delta y_0} G_{-\frac12} \phi_{0,-\frac12} \,. \end{align}

where \( \mathbf{F} = (F, G) \) (flux_coef) is the face diffusive flux coefficients, and \( \delta x_0 \) and \( \delta y_0 \) is the size of the cell \( (0,0) \).

Note

As described in Cell diffusion term, a central schema of order 2 is used (regular mesh).
No treatment of boundary condition is done here.

Todo: check order with irregular grids

◆ apply_cell_diffusion_o2_centered_scheme()

subroutine mod_apply_cell_diffusion_scheme::apply_cell_diffusion_o2_centered_scheme	(	type(t_face_field), intent(in)	face_diffusion_coefficient,
		double precision, dimension(-ns:ns,-ns:ns,-ns:ns), intent(inout)	matrix_line,
		integer, intent(in)	ns,
		integer, intent(in)	i,
		integer, intent(in)	j,
		integer, intent(in)	k,
		double precision, dimension(:,:,:), intent(in), optional	extra_inside_gradient_coef2 )

The term writes:

\begin{align} \nabla \cdot ( D \nabla \phi)_{0,0} =& \frac{1}{\delta x_0} F_{+\frac12} \phi_{+\frac12,0} -\frac{1}{\delta x_0} \left(F_{+\frac12} + F_{-\frac12}\right) \phi_{0,0} +\frac{1}{\delta x_0} F_{-\frac12} \phi_{-\frac12,0} \\ +& \frac{1}{\delta y_0} G_{+\frac12} \phi_{0,+\frac12} -\frac{1}{\delta y_0} \left(G_{+\frac12} + G_{-\frac12}\right) \phi_{0,0} +\frac{1}{\delta y_0} G_{-\frac12} \phi_{0,-\frac12} \,. \end{align}

where \( \mathbf{F} = (F, G) \) (flux_coef) is the face diffusive flux coefficients, and \( \delta x_0 \) and \( \delta y_0 \) is the size of the cell \( (0,0) \).

Note

As described in Cell diffusion term, a central schema of order 2 is used (regular mesh).
No treatment of boundary condition is done here.

Todo: check order with irregular grids

◆ harmonic_diffusion_coef_interpolation()

subroutine mod_harmonic_diffusion_coef_interpolation::harmonic_diffusion_coef_interpolation	(	type(t_face_field), intent(inout)	face_diffusion_coef,
		double precision, dimension(nx,ny,nz), intent(in)	cell_diffusion_coef )

The face diffusion coefficients is obtained by harmonic interpolation, as explained in Cell diffusion term.

\begin{align} D_\frac12 = \dfrac{ \tfrac12\delta x_0 + \tfrac12\delta x_1 }{ \dfrac{\tfrac12\delta x_0}{D_0} + \dfrac{\tfrac12\delta x_1}{D_1} }\,. \end{align}

◆ interpolate_diffusion_coef()

subroutine mod_interpolate_diffusion_coef::interpolate_diffusion_coef	(	type(t_face_field), intent(inout)	face_diffusion_coef,
		double precision, dimension(nx,ny,nz), intent(in)	cell_diffusion_coef )

The face diffusion coefficients is obtained by harmonic interpolation, as explained in Cell diffusion term.

\begin{align} D_\frac12 = \dfrac{ \tfrac12\delta x_0 + \tfrac12\delta x_1 }{ \dfrac{\tfrac12\delta x_0}{D_0} + \dfrac{\tfrac12\delta x_1}{D_1} }\,. \end{align}

◆ set_diffusive_flux_coef()

subroutine mod_set_diffusive_flux_coef::set_diffusive_flux_coef ( type(t_face_field), intent(inout) flux_coef )

The face diffusive flux coefficient is obtained by

\begin{align} F_\frac12 &= \frac{D_{\frac12,0,0}}{\delta x_\frac12}\,, & G_\frac12 &= \frac{D_{0,\frac12,0}}{\delta y_\frac12}\,, & H_\frac12 &= \frac{D_{0,0,\frac12}}{\delta z_\frac12}\,. \end{align}

Functions