FVM and OpenFoam

In the previous session, we have used openFoam to resolve hydrothermal convection in oceanic crust. In this lecture, we will have a deeper look into how openFoam solves the governing equations and how the numerical schemes work.

OpenFoam uses the Finite Volume Method (FVM), which is a popular numerical technique used in computational fluid dynamics (CFD) and other areas of computational physics to solve partial differential equations (PDEs) that describe physical phenomena. It’s particularly well-suited for problems involving fluid flow, heat transfer, and associated physical processes. Here’s an overview of how FVM works:

Basic Concept

Domain Discretization The computational domain (the area or volume where the physical problem occurs) is divided into small, discrete control volumes (cells). These control volumes cover the entire domain without overlapping.

Integral Form of PDEs FVM operates on the integral form of the governing PDEs, as opposed to the differential form used in methods like the Finite Difference Method (FDM). The equations are integrated over each control volume.

Flux Calculations The key aspect of FVM is the calculation of fluxes of the conserved quantities like mass, momentum, energy across the faces of each control volume. These fluxes are used to determine how these quantities change within each control volume over time.

Boundary Conditions The method also incorporates boundary conditions, which define how the field behaves at the boundaries of the computational domain.

The general workflow in FVM is the following: Represent the variables like velocity, pressure, temperature at discrete locations in each control volume - either at the cell centers or at the cell faces. The integral form of the governing PDEs is discretized for each control volume. This involves expressing the rate of change of a quantity in a control volume in terms of the fluxes through its faces. The fluxes through the control volume faces are approximated using values of the variables at the centers or faces of the cells. Various schemes like upwind and central difference can be used for this approximation. The discretized equations form a system of algebraic equations. These equations are solved iteratively to find the field values (like pressure, velocity, etc.) throughout the domain. For steady-state problems, iterations continue until the solution converges. For transient problems, the solution is marched forward in time, updating the variables at each time step. One key advantage of the FVM is that it ensures conservation of physical quantities locally (in each control volume) and globally (across the entire domain); it can further handle complex geometries and irregular meshes. In summary, FVM is a versatile and powerful tool for solving complex physical problems numerically.

Openfoam uses what is called a cell-centered finite volume method. This means that the variables are stored at the center of the cells. The fluxes are calculated at the cell faces. The fluxes are then used to update the variables at the cell centers.

Example problem: heat diffusion

Let’s look at an example problem that solves for temperature diffusion, using the laplacianFoam solver.

Geometry and governing equations

The governing equation is the heat diffusion equation, which is given by:

(48)\[\rho c_p\frac{\partial T}{\partial t} - \nabla \cdot k \nabla T = 0\]

(49)\[\frac{\partial T}{\partial t} - \nabla \cdot D \nabla T =0\]

\(k\) is the thermal conductivity and \(D\) is the thermal diffusivity \(\frac{k}{\rho c_p}\) . We here assume that \(D\) is a constant value. See also Theory.

Below is a simple 2D setup for a heat diffusion test case. We will use it to learn about the details of how openFOAM’s FV method works.

Fig. 39 Model geometry, boundary conditions and initial condition.

What we need to solve ?

Temperature field \(T\) of the model region at time \(t\), therefore \(T\) is the only unknown variable of our problem. It is a scalar field at the center of each cell.

Domain discretization

Mesh structure and topology

We can use a simple blockMesh to discretize our domain. You can download a prepared case here: Regular box case.

Regular mesh

Fig. 40 2D regular mesh structure and topology information of OpenFOAM polyMesh (See 1. Read and plot mesh information in the notebook).

Unstructured mesh

Fig. 41 2D unstructured mesh topology information (Regular box case).

The mesh comprises cell centers, cell faces, and cell vertices. The cell centers are the points where the temperature field is stored. The cell faces are the boundaries of the cells. The cell vertices are the points that define the cell faces. Notice how the internal faces (the ones that are not boundary conditions) are numbered continuously starting from 0.

In addition, a connectivity is needed to describe the topology of the mesh. This connectivity describes which cells are connected to each other through the faces. In openFoam a highly efficent data structure is used to store the connectivity information. For each face, the owner cell and the neighbour cell are stored. The normal vector of face always points from the owner cell to the neighbour cell. The numbering of the vertices that make up the face is again done using the right-hand rule and the face normal always points from the cell with the lower index to the cell with the larger index. Boundary faces have only an owner cell.

You can find details on the mesh description in the OpenFoam manual. Have a look at it and check the structure of the files in constant/polymesh.

Tip

Check your understanding of the openFoam mesh structure and topology by creating the blockMesh mesh for the case we have downloaded and visualizing it using paraview. To create the labels, use the GenerateIds filter in paraview.

Spatial discretization: The diffusion term

The equation discretization step is performed over each element of the computational domain to yield an algebraic relation that connects the value of a variable in an element to the values of the variable in the neighboring elements [Moukalled et al., 2016].

Goal

Transform partial differential equation (49) into a set of algebraic equations: \(\mathbf{A} T = \mathbf{b}\).

1. Integrate equation (49) over element \(C\) (the green cell in Fig. 40) that enables recovering its integral balance form as,

(50)\[\iiint\limits_{V_C} \frac{\partial T}{\partial t} dV - \iiint\limits_{V_C} \nabla \cdot D\nabla T dV = 0\]

Note that the transient term \(\iiint\limits_{V_C} \frac{\partial T}{\partial t} dV\) will be processed in the later section. Let’s set it to zero for the moment and assume steady-state.

2. Transform the volume integral of the heat flux into a surface integral by applying the divergence theorem,

(51)\[- \iiint\limits_{V_C} \nabla \cdot D\nabla T dV = -\oint\limits_{\partial V_C} (D\nabla T)\cdot d\vec{S} = 0\]

The equation (51) is actually a heat balance over cell \(C\). It is basically the integral form of the original partial differential equation and involves no approximation.

The integrant is the diffusive heat flux,

(52)\[\vec{J}^{T,D} \equiv -D\nabla T\]

3. Transform the surface integral in equation (51) as a summation over the control volume faces (still no approximation),

(53)\[\sum\limits_{f\sim faces(V_C)} \iint\limits_{f}\vec{J}^{T,D}_f \cdot d\vec{S}_f =0\]

Here \(f\) denotes the boundary face of cell \(V_C\).

Now, the key problem is how to calculate flux integral on a face.

4. Evaluating flux integral on the faces

It’s worth noting that there are only two kinds of cells in the discretized domain, one is internal cell, the other is boundary cell. All faces of a internal cell are internal faces. The remain cells are boundary cells which contain at least one boundary face.

Internal cell

For a specific internal cell \(C\), e.g. cell 12, shown in Fig. 40, the equation (53) can be expanded as,

(54)\[\color{red}{\iint_{f_{23}}\vec{J}^{T,D}_{f_{23}} \cdot \vec{S}_{f_{23}}} + \iint_{f_{24}}\vec{J}^{T,D}_{f_{24}} \cdot \vec{S}_{f_{24}} + \iint_{f_{21}}\vec{J}^{T,D}_{f_{21}} \cdot \vec{S}_{f_{21}} + \iint_{f_{5}}\vec{J}^{T,D}_{f_{5}} \cdot \vec{S}_{f_{5}} = 0\]

Note that the surface normals \(\vec{S}\) always point out of the cell; for example, \(\vec{S}_{f_{23}}\) would have a positive and \(\vec{S}_{f_{21}}\) a negative sign.

4.1 Considering face \(f_{23}\), calculate the first term on the right hand side of the equation (54) (we have to introduce the first approximation at this step). Using a Gaussian quadrature the integral at the face \(f_{23}\), for example, becomes,

(55)\[\color{red}{\iint_{f_{23}}\vec{J}^{T,D}_{f_{23}} \cdot \vec{S}_{f_{23}}} = \iint_{f_{23}} (\vec{J}^{T,D}_{f_{23}} \cdot \vec{n}_{f_{23}}) dS_{f_{23}} \approx \color{orange}{\sum\limits_{ip\sim ip(f_{23})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega_{ip} S_{f_{23}}}\]

where \(S_f\) is the area of face \(f_{23}\), \(ip\) refers to a integration point and \(ip(f_{23})\) the number of integration points along surface \(f_{23}\), \(\omega_{ip}\) is the integral weights.

Fig. 42 Surface integration of fluxes using (a) one integration point, (b) two integration points, and (c) three integration points [Moukalled et al., 2016]. The total flux \(FluxT_f\) is computed as linear combinations of the two cell values plus a source term.

Important

Only the one integration point scheme, i.e. Fig. 42 (a), is implemented in OpenFOAM (at least before version 8). Therefore the first keyword of Laplacian scheme in system/fvScheme dictionary file of a case has only one option, it is Gauss.

4.2 Choose integral scheme or integral points

To simply explain the calculation process/logic and to also be consistent with OpenFOAM, we here adopt an one integration point scheme with weight \(\omega = 1\), thus equation (55) becomes,

(56)\[\color{orange}{\sum\limits_{ip\sim ip(f_{23})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega S_f } = -\left(D \frac{\partial T}{\partial x} \vec{i} + D\frac{\partial T}{\partial y} \vec{j} \right)_{f_{23}} \cdot \Delta y_{f_{23}} \vec{i} = -\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{23}} \Delta y_{f_{23}}}\]

where

• \(S_{f_{23}} = \Delta y_{f_{23}}\) is the area of face \(f_{23}\) by assuming \(\Delta z = 1\).

• \(\vec{S}_{f_{23}} = S_{f_{23}} \vec{n}_{f_{23}}\) is the surface vector of face \(f_{23}\).

• \(\vec{n}_{f_{23}} = \vec{i}\) is the normal vector of the face \(f_{23}\) directed out of the cell \(C\) (see Fig. 40).

4.3 Calculate gradient of \(T\) at face centroid, here introduce the second approximation, e.g. assuming linear variation of T and then the gradient term in the equation (56) can be approximated as,

\[-\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{23}} \Delta y_{f_{23}}} = -D\frac{T_{13} - T_{C}}{x_{13} - x_C}\Delta y_{f_{23}} = -\frac{D \Delta y_{f_{23}}} {\delta x_{f_{23}}} (T_{13} - T_{C})\]

where \(\delta x_{f_{23}}\) represents the distance between center of cells who share face \(f_{23}\), they are cell 12 and cell 13.

Now let’s look at the western face \(f_{21}\) . Using a single integration point, we can again write:

(57)\[\color{orange}{\sum\limits_{ip\sim ip(f_{21})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega S_f } = -\left(D \frac{\partial T}{\partial x} \vec{i} + D\frac{\partial T}{\partial y} \vec{j} \right)_{f_{21}} \cdot \Delta y_{f_{21}} \cdot-\vec{i} = \color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{21}} \Delta y_{f_{21}}}\]

Notice how the surface normal now has the opposite sign!

We continue to spell out the flux across face \(f_{21}\) as:

\[\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{21}} \Delta y_{f_{21}}} = D\frac{T_{C} - T_{11}}{x_{11} - x_C}\Delta y_{f_{21}} = -\frac{D \Delta y_{f_{21}}} {\delta x_{f_{21}}} (T_{11} - T_{C})\]

Do the same thing for the remain faces (\(f_5, f_{24}\)).

4.3 Construct coefficients of a specific cell \(C_{12}\)

Now we get all the coefficients, thus the equation (54) can be discretized and expressed as a matrix form,

(58)\[a_{C} T_C + a_{F_{23}} T_{F_{23}} + a_{F_{24}} T_{F_{24}} + a_{F_{21}} T_{F_{21}} + a_{F_{5}} T_{F_{5}} = 0\]

with

\[\begin{split}\begin{eqnarray} a_{23} & = & -\frac{D \Delta y_{f_{23}}} {\delta x_{f_{23}}} \\ a_{24} & = & -\frac{D \Delta y_{f_{24}}} {\delta x_{f_{24}}} \\ a_{21} & = & -\frac{D \Delta y_{f_{21}}} {\delta x_{f_{21}}} \\ a_{5} & = & -\frac{D \Delta y_{f_{5}}} {\delta x_{f_{5}}} \\ a_{C} & = & -(a_{23} + a_{24} + a_{21} + a_{5}) \\ \end{eqnarray}\end{split}\]

(59)\[\mathbf{A}_{N\times N}[12,:] \mathbf{T}_{N\times 1} = 0\]

Boundary cell

For a specific boundary cell, e.g. cell 19 shown in Fig. 40, the equation (53) can be rewritten as,

(60)\[\iint_{f_{18}}\vec{J}^{T,D}_{f_{18}} \cdot \vec{S}_{f_{18}} + \iint_{f_{35}}\vec{J}^{T,D}_{f_{35}} \cdot \vec{S}_{f_{35}} +\iint_{f_{37}}\vec{J}^{T,D}_{f_{37}} \cdot \vec{S}_{f_{37}} + \color{red}{\iint_{f_{51}}\vec{J}^{T,D}_{f_{51}} \cdot \vec{S}_{f_{51}}} = 0\]

The first three terms are normal fluxes across internal faces, so there is nothing special and we can just calculate them following steps for internal faces as before. The special thing is the boundary face \(f_{52}\), we can call it \(f_b\) for a general case. Now the problem is how to evaluate flux on the boundary face \(f_b\).

The fluxes on the interior faces are discretized as before, while the boundary flux is discretized with the aim of constructing a linearization with respect to the cell field \(T_C\), e.g. \(T_{C_{19}}\) of cell \(C_{19}\) shown in Fig. 40, thus

(61)\[\color{red}{\iint_{f_{b}}\vec{J}^{T,D}_{f_{b}} \cdot d\vec{S}_{f_{b}}} = a_{F_b} T_C + c_{F_b}\]

All right, now out goal is to determine the coefficients of \(a_{F_b}\) and \(c_{F_b}\) according to the boundary conditions. Generally there are two basic kinds of boundary conditions, they are Dirichlet boundary condition (also called the first type or fixed value boundary condition) and Von Neumann boundary condition (all called the second type or fixed flux boundary condition), respectively.

4.1 Fixed value boundary condition

Fixed value means the value of field \(T\) is given on the boundary face is give, i.e. \(\color{purple}{T_{f_b}}\) is given. There fore we can calculate flux on boundary face, \(f_{51}\) for example,

(62)\[\color{red}{\iint_{f_{51}}\vec{J}^{T,D}_{f_{51}} \cdot \vec{S}_{f_{51}}} = \iint_{f_{51}} (\vec{J}^{T,D}_{f_{51}} \cdot \vec{n}_{f_{51}}) dS_{f_{51}} \approx \color{orange}{\sum\limits_{ip\sim ip(f_{51})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega_{ip} S_{f_{51}}}\]

Here we still use one integration point scheme to explain the calculation process,

(63)\[\color{orange}{\sum\limits_{ip\sim ip(f_{51})} (\vec{J}^{T,D} \cdot \vec{n}_{ip})\omega S_f } = -\left(D \frac{\partial T}{\partial x} \vec{i} + D\frac{\partial T}{\partial y} \vec{j} \right)_{f_{51}} \cdot \Delta y_{f_{51}} \vec{i} = -\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{51}} \Delta y_{f_{51}}}\]

Now let’s calculate the gradient on the boundary face,

(64)\[-\color{blue}{\left( D\frac{\partial T}{\partial x} \right)_{f_{51}} \Delta y_{f_{51}}} = -D\frac{\color{purple}{T_{f_{51}}} - T_{C}}{x_{f_{51}} - x_C}\Delta y_{f_{51}} = -\frac{D \Delta y_{f_{51}}}{\delta x_{f_{51}}}(\color{purple}{T_{f_{51}}} - T_{C}) = a_{f_{51}}(-T_{C}) + c_{f_{51}}\]

Substituting equation (64) back to equation (63), equation (62) and equation (61), we can get Laplacian discretization coefficients \(a_{f_{51}}\) and \(c_{f_{51}}\),

(65)\[\begin{split}\begin{eqnarray} a_{F_{51}} &=& -D\frac{\Delta y_{f_{51}}}{\delta x_{f_{51}}}\\ c_{F_{51}} &=& a_{F_{51}}\color{purple}{T_{f_{51}}} \end{eqnarray}\end{split}\]

Notice how \(\color{purple}{T_{f_{51}}}\) is known (a boundary condition), so that the term \(-a_{F_{51}}\color{purple}{T_{f_{51}}}\) will end up on the right-hand side and not in the coefficient matrix.

4.2 Fixed flux boundary condition

Fixed flux means normal gradient to the boundary face of field \(T\) is given, i.e. \(\color{purple}{\vec{J}^{T,D}_{f_{b}} \cdot \vec{n}_{f_{b}} = q_{f_b}}\) is given. There fore we can calculate flux on boundary face, \(f_{51}\) for example,

(66)\[\color{red}{\iint_{f_{51}}\vec{J}^{T,D}_{f_{51}} \cdot \vec{S}_{f_{51}}} = \iint_{f_{51}} (\vec{J}^{T,D}_{f_{51}} \cdot \vec{n}_{f_{51}}) dS_{f_{51}} = \iint_{f_{51}} \color{purple}{q_{f_{51}}} dS_{f_{51}} = \color{purple}{q_{f_{51}}} S_{f_{51}} = \color{purple}{q_{f_{51}}} \Delta y_{f_{51}}\]

Substituting equation (66) back to equation (61), we can get Laplacian discretization coefficients \(a_{f_{51}}\) and \(c_{f_{51}}\),

(67)\[\begin{split}\begin{eqnarray} a_{F_{51}} &=& 0 \\ c_{F_{51}} &=& \color{purple}{q_{f_{51}}} \Delta y_{f_{51}} \end{eqnarray}\end{split}\]

4.3 Construct coefficients of the boundary cell \(C_{19}\)

Do the same thing as 4.1 or 4.2 (depends on the type of boundary condition) for all boundary faces of a boundary cell \(C\) and calculate the coefficients of \(a_{F_b}\) and \(c_{F_b}\), then the equation (60) can be expressed as,

(68)\[\begin{split} \begin{eqnarray} b_{C_{19}} & = & a_{F_{18}} (T_{F_{18}} - T_C) + a_{F_{35}} (T_{F_{35}} - T_C) + a_{F_{37}} (T_{F_{37}} - T_C) + a_{F_{51}} (-T_C) +c_{F_{51}} \\ & = & \sum\limits_{f\in[18, 35, 37]} a_{F_f} T_{F_f} + \left(-\sum\limits_{f\in[18, 35, 37, 51]} a_{F_f} \right) T_C + \sum\limits_{f\in[51]} c_{F_f} \\ & = & \sum\limits_{f\sim nb(\text{internal face})} a_{F_f} T_{F_f} + \left(-\sum\limits_{f\sim nb(\text{all face})} a_{F_f} \right) T_C + \sum\limits_{f\sim nb(\text{boundary face})} c_{F_f} \end{eqnarray}\end{split}\]

(69)\[\mathbf{A}_{N\times N}[19,:] \mathbf{T}_{N\times 1} = B_{19}\]

where \(B_{19} = b_{C_{19}} - \sum\limits_{f\sim nb(\text{boundary face})} c_{F_f}\).

5. Construct general matrix form of the equation (53)

(70)\[\begin{split}\begin{eqnarray} -\sum\limits_{f\sim nb(C)} (D\nabla T)_{F_f}\cdot \vec{S}_{F_f} & = & \sum\limits_{f\sim \text{internal faces}(C)} a_{F_f} T_{C} + \left(-\sum\limits_{f\sim nb(C)} a_{F_f} T_{F_f} \right) + \sum\limits_{f\sim \text{boundary faces}(C)} c_{F_{F_f}} \\ & = & a_CT_C - \sum\limits_{f\sim \text{internal faces}(C)} a_{F_f} T_{F_f} + \sum\limits_{f\sim \text{boundary faces}(C)} c_{F_{f}} \end{eqnarray}\end{split}\]

Matrix form of the equation (53) can be written as,

(71)\[\mathbf{A}_{N \times N} \mathbf{T}_{N\times1} = \mathbf{B}_{N\times1}\]

here \(N\) is the number of cell, \(\mathbf{A}\) is a sparse matrix of coefficients, \(\mathbf{T}\) is cell-centered temperature field. The matrix is visualized as Fig. 43. It should be noted that the coefficients matrix of Laplacian term is a symmetric matrix and the boundary conditions only affect diagonal entry of the coefficients matrix. Further, only the fixed value boundary condition affects the matrix and the fixed flux boundary condition affects the RHS(Right hand of side) \(B\).

Regular mesh

Fig. 43 Visualization of \(\mathbf{A}_{N \times N} \mathbf{T}_{N\times1} = \mathbf{b}_{N\times1}\). The coefficients of the selected cell \(C\) (Fig. 40) are marked by green rectangle.

Unstructured mesh

Fig. 44 Visualization of \(\mathbf{A}_{N \times N} \mathbf{T}_{N\times1} = \mathbf{b}_{N\times1}\). The coefficients of the selected cell \(C\) (Fig. 40) are marked by green rectangle.

Temporal Discretization: The Transient Term

After spatial discretization, the equation (50) can be expressed as,

(72)\[\frac{\partial T}{\partial t} V_C - L(T_C^t) =0\]

where \(V_C\) is the volume of the discretization cell and \(L(T_C^t)\) is the spatial discretization operator expressed at some reference time \(t\), which can be written as algebraic form (see also equation (70)),

(73)\[L(T_C^t) = a_C T_C^t - \sum\limits_{F\sim NB(C)} a_F T_F^t + \sum\limits_{F\sim NB(C)} c_F\]

where \(a_C\) is the diagonal coefficients of the matrix, \(a_F\) is the off-diagonal coefficients, and the \(c_F\) is the source coefficients as right hand side of matrix of system.

Tip

For specific cell \(C\):

• \(a_F\) has only contributed from internal faces.

• \(a_C\) is the negative summation of \(a_F\), contributed from all faces, but equation of the \(a_F\) for a boundary faces (equation (65) and equation (67)) is a little bit different from internal face.

• \(c_F\) only comes from boundary faces of cell \(C\) if it has boundary face, see also equation (65) and equation (67). \(c_F\) will contribute to RHS (\(\mathbf{B}\)) of the algebraic equation (71).

Let’s do the temporal discretization for the first term of equation (72),

\[\frac{T^{t+\Delta t/2} - T^{t - \Delta t/2}}{\Delta t} V_C + L(T_C^t) = 0\]

To derive the full discretized equation, an interpolation profile expressing the face values at (\(t-\Delta t/2\)) and (\(t+\Delta t\)) in terms of the element values at (\(t\)), (\(t-\Delta t\)), etc., is needed.Independent of the profile used, the flux will be linearized based on old and new values as,

(74)\[b_C = FluxC T_C + FluxC^o T_C^o + FluxV\]

where the superscript \(^{o}\) refers to old values. With the format of the linearization, the coefficient \(FluxC\) will be added to diagonal element and the coefficient \(FluxC^o T_C^o + FluxV\) will be added to the source or RHS (\(B\)).

First order implicit Euler scheme

(75)\[\frac{T^t - T^{t-\Delta t}}{\Delta t} V_C - L(T^{t}) = 0\]

Therefore the coefficients will be

(76)\[FluxC = \frac{V_C}{\Delta t}, ~ FluxC^o = - \frac{V_C}{\Delta t}, ~ FluxV =0\]

First order explicit Euler scheme

(77)\[\frac{T^{t+\Delta t} - T^t}{\Delta t} V_C - L(T^{t}) = 0\]

Note that now the new time is at \(t+\Delta t\) and that the spatial operator (\(L\)) of equation (77) has to be evaluated at time \(t\). Therefore, in order to find the value of \(T\) at time \(t+\Delta t\), we don’t need to solve a set of linear algebraic equations,

(78)\[T^{t+\Delta t} = \frac{\Delta t}{V_C} L(T^{t}) + T^t\]

Building the matrix: Mesh connectivity and internal data storage

Remember the mesh structure discussed above. The mesh connectivity is stored in the polyMesh folder, which containes the neighbour and owner files. The normal vector of the face always points from the owner cell to the neighbour cell. The numbering of the vertices that make up the face is again done using the right-hand rule and the face normal always points from the cell with the lower index to the cell with the larger index. Boundary faces have only an owner cell.

Is is straightforward to find the indices for setting up the flux balance for each cell, as we did in equation (58)?

No, because the owner/neighbor data structure does not provide the indices of all the cells that are connected to a given cell. Rather the owner/neighbor data structure is “face-centered” in that for each face, we know which cells are on each side. Hence, we can easily loop over all faces and compute the fluxes but we cannot easily loop over cells and compute the fluxes to/from the neighboring cells.

Tip

Check out how openfoam makes use of this face-centered logic to compute the average gradient of a field in the gaussGrad function in src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C (see below). The forAll loop is on the faces and the owner and neighbour functions are used to get the indices of the cells that are connected to a given face.

Listing 26 src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "gaussGrad.H"
#include "extrapolatedCalculatedFvPatchField.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

template<class Type>
Foam::tmp
<
   Foam::GeometricField
   <
      typename Foam::outerProduct<Foam::vector, Type>::type,
      Foam::fvPatchField,
      Foam::volMesh
   >
>
Foam::fv::gaussGrad<Type>::gradf
(
   const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,
   const word& name
)
{
   typedef typename outerProduct<vector, Type>::type GradType;

   const fvMesh& mesh = ssf.mesh();

   tmp<GeometricField<GradType, fvPatchField, volMesh>> tgGrad
   (
      new GeometricField<GradType, fvPatchField, volMesh>
      (
            IOobject
            (
               name,
               ssf.instance(),
               mesh,
               IOobject::NO_READ,
               IOobject::NO_WRITE
            ),
            mesh,
            dimensioned<GradType>
            (
               "0",
               ssf.dimensions()/dimLength,
               Zero
            ),
            extrapolatedCalculatedFvPatchField<GradType>::typeName
      )
   );
   GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad.ref();

   const labelUList& owner = mesh.owner();
   const labelUList& neighbour = mesh.neighbour();
   const vectorField& Sf = mesh.Sf();

   Field<GradType>& igGrad = gGrad;
   const Field<Type>& issf = ssf;

   forAll(owner, facei)
   {
      GradType Sfssf = Sf[facei]*issf[facei];

      igGrad[owner[facei]] += Sfssf;
      igGrad[neighbour[facei]] -= Sfssf;
   }

   forAll(mesh.boundary(), patchi)
   {
      const labelUList& pFaceCells =
            mesh.boundary()[patchi].faceCells();

      const vectorField& pSf = mesh.Sf().boundaryField()[patchi];

      const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];

      forAll(mesh.boundary()[patchi], facei)
      {
            igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];
      }
   }

   igGrad /= mesh.V();

   gGrad.correctBoundaryConditions();

   return tgGrad;
}


template<class Type>
Foam::tmp
<
   Foam::GeometricField
   <
      typename Foam::outerProduct<Foam::vector, Type>::type,
      Foam::fvPatchField,
      Foam::volMesh
   >
>
Foam::fv::gaussGrad<Type>::calcGrad
(
   const GeometricField<Type, fvPatchField, volMesh>& vsf,
   const word& name
) const
{
   typedef typename outerProduct<vector, Type>::type GradType;

   tmp<GeometricField<GradType, fvPatchField, volMesh>> tgGrad
   (
      gradf(tinterpScheme_().interpolate(vsf), name)
   );
   GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad.ref();

   correctBoundaryConditions(vsf, gGrad);

   return tgGrad;
}


template<class Type>
void Foam::fv::gaussGrad<Type>::correctBoundaryConditions
(
   const GeometricField<Type, fvPatchField, volMesh>& vsf,
   GeometricField
   <
      typename outerProduct<vector, Type>::type, fvPatchField, volMesh
   >& gGrad
)
{
   typename GeometricField
   <
      typename outerProduct<vector, Type>::type, fvPatchField, volMesh
   >::Boundary& gGradbf = gGrad.boundaryFieldRef();

   forAll(vsf.boundaryField(), patchi)
   {
      if (!vsf.boundaryField()[patchi].coupled())
      {
            const vectorField n
            (
               vsf.mesh().Sf().boundaryField()[patchi]
            / vsf.mesh().magSf().boundaryField()[patchi]
            );

            gGradbf[patchi] += n *
            (
               vsf.boundaryField()[patchi].snGrad()
            - (n & gGradbf[patchi])
            );
      }
   }
}


// ************************************************************************* //

OpenFoam uses a matrix format that is tightly integrated with the “face-centered” and the owner/neighbor structure. It is called the LDU matrix format. This format is a way of storing sparse matrices. A sparse matrix is one in which most elements are zero, which is a common case in CFD (check the example above). The LDU format is specifically tailored to store and manipulate these types of matrices efficiently.

Components of LDU:

• L (Lower triangle): This part of the matrix stores the coefficients that are below the main diagonal.

• D (Diagonal): This represents the main diagonal of the matrix. In many algorithms, the diagonal elements play a crucial role and are treated separately for reasons of numerical stability and efficiency.

• U (Upper triangle): This stores the coefficients above the main diagonal.

The LDU format is memory-efficient because it only stores non-zero elements. This is particularly beneficial in CFD where matrices can be very large, but only a small fraction of the elements are non-zero. OpenFOAM uses this structure to perform matrix operations like multiplication, addition, and especially the solving of linear systems, which is crucial in the iterative methods used in CFD simulations.

The LDU format uses scalar fields representing the diagonal, upper, and lower coefficients. The diagonal coefficients are indexed using the cell index. The upper and lower coefficients use face indices. In order to get the indices to match, mapping arrays between the face and cell indices are provided. The lowerAddr(), and upperAddr() functions of the lduAdressing class provide the owner and neighbor cell indices for each face.

A nice summary of the LDU format can be found in the OpenFoam Wiki.

Summary

In summary, the LDU format is a way of storing sparse matrices. It is specifically tailored to store and manipulate these types of matrices efficiently. OpenFOAM uses this structure to perform matrix operations like multiplication, addition, and especially the solving of linear systems, which is crucial in the iterative methods used in CFD simulations. Most of this happens under the hood but if you ever want to perform an operation on all rows of an lduMatrix, you need to know abou lduAdressing.