A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction Equation in High Contrast Domains

Abstract

Implementation of finite element method (FEM) needs special cares, particularly for essential boundary conditions that have an important effect on symmetry and number of unknowns in the linear systems. Moreover, avoiding numerical integration and using (off-line) calculated element integrals decrease the computational cost significantly. In this chapter we briefly present theoretical topics of FEM. Instead we focus on what is important (and how) to carefully implement FEM for equations that can be the core of a numerical simulator for a diffusion–advection-reaction problem. We consider general 2D and 3D domains, high contrast and heterogeneous diffusion coefficients and generalize the method to nonlinear parabolic equations. Although we use Matlab codes to simplify the explanation of the proposed method, we have implemented it in C++ to reveal the efficiency and examples are presented to admit it.

1. Introduction

In a domain $\mathrm{\Omega }$ in ${\mathbb{R}}^n,n=\mathrm{1,2,3},$ consider a partial differential equation of the form of

$$\nabla \cdot \left(-\kappa \nabla u\right)+\gamma \cdot \nabla u+\mu u=f\text{in}\mathrm{\Omega },$$E1

where the unknown function $u$, reaction coefficient $\mu$ and the given function $f$ are scalar functions $\left(\mathrm{\Omega }\to \mathbb{R}\right)$. Diffusion coefficient $\kappa$ and the (convergence free) advection coefficient $\gamma$ are defined in $\mathrm{\Omega }$ and give (normally piecewise constant) values in ${\mathbb{R}}^{n\times n}$ and ${\mathbb{R}}^n$, respectively. Eq. (1) is called a diffusion–advection-reaction equation and has immense applications in science and engineering. Transport of heat, momentum and energy, solid mechanics, CO₂ sequestration, computational fluid dynamics and biophysics are a few number of research fields that as a part of the solution strategy need to solve (1) numerically.

After several steps, numerical solution of Eq. (1) is obtained by solving a linear system

$$\mathcal{A}x=b$$E2

where $\mathcal{A}$ and $b$ are called stiffness matrix and right hand side vector, respectively. In practice each term of diffusion, advection or reaction is accumulated separately in the stiffness matrix.

Several challenges make finding the numerical solution of (1) difficult. 3 of the most important are 1) complex geometry of $\mathrm{\Omega }$, 2) heterogeneity and high contrast in coefficient $\kappa$ that produce very ill-conditioned linear systems and 3) large scale domains. A very common numerical method to solve (1) is FEM that uses a mesh representing the complex geometries accurately and overcomes the first issue. We can employ powerful linear solvers such as multigrid or multiscale methods to resolve the effect of heterogeneity and high contrast in $\kappa$ in the linear system. Finally a careful implementation in parallel machines reduces computational time of large scale simulations considerably. However, lots of effort and research are still needed to propose a method obtaining a reliable solution in a reasonable time for real problems. Other issues such as uncertainties in the data or nonlinear coupled system of equations should be addressed, as well.

Considering a positive definite $\kappa$ we explain (and implement) how to obtain a finite element solution of (1) through the following steps.

1. In Matlab we generate arbitrary 2D and simple 3D domains. The main domain is divided into a collection of subdomains, called elements where each element includes some nodes. Numbered nodes and elements generate a mesh. Actually by a mesh we mean 2 data structures, Cells (integer valued) that stores identifier or index of nodes in each element and Nodes (real valued) that stores coordinates of each node. Many software such as Matlab, GID and Gmesh generate reliable meshes. Corresponding to each element $e$, a positive definite matrix (a constant or a $2\times 2$ matrix in 2D or a $3\times 3$ matrix in 3D), named ${\kappa }_e$, will be set to form $\kappa$. In each element $e$, we can also set ${\mu }_e$ and ${\gamma }_e$ to form $\mu$ and $\gamma$ in Eq. (1), respectively. See Section 2 for details.

2. A brief introduction to FEM is presented that covers weak formulation (Section 3), shape functions and reference elements (Section 4).

3. Evaluating of element integrals and preparing table of calculated integrals are discussed in Section 4.1. Then with help of affine mappings (Section 4.2) we use a linear combination of element integrals to accumulate into the stiffness matrix which is explained in Section 5. Note that different types of elements (for example triangles and rectangles in 2D) might exist in the mesh and we can have an unstructured mesh, generally. Hence element integrals should be considered for any type of elements exist in the mesh.

4. Efficient implementation of boundary conditions, particularly Dirichlet boundary condition is presented in Section 5.2. We show that how a correct but careless implementation of Dirichlet boundary conditions decreases both accuracy and efficiency.

5. After assembly of the linear system (2), we use Matlab functions to solve the linear system and plot the result. Linear solvers are the core of a numerical simulator and their efficiency has a direct impact on the overall simulation. We refer to [1] for further discussions and references.

6. Generalization to more complex equations such as Eq. (25) now becomes straightforward and solving strategy is given in Section 6.

We finish this chapter by introducing the necessary topics to have an independent and efficient simulator in Section 7. Most of the topics presented here are fully discussed in [1, 2, 3] which are of great value for further reading.

2. Mesh generation

Implementation of FEM is began with mesh generation which is dividing the main domain into several subdomains or elements. Each element is finite (finite elements) and has a regular shape such as triangle or cuboid. Each element has a specific type, such as linear or bilinear, which is defined by shape functions and explained in Section 4. We prefer to use Lagrange triangles and parallelograms in 2D and tetrahedron and cuboid in 3D space as elements. The reason of such selection is to avoid numerical integration to evaluate definite integrals (of Lagrange elements) arise in FEM.

The main output of the mesh generation is two sets, nodes and elements that form the mesh. Nodes are not necessarily the vertices of the elements, for example in linear Crouzeix-Raviart element nodes are placed at edge midpoints or in quadratic Lagrange (triangle) elements midpoints and vertices together form the nodes.

In presence of complex geometries mesh generation can be a challenge for software mesh generators. For not very complicated geometry, Matlab, with a good quality, can generate triangular and tetrahedron mesh for 2D and 3D domains, respectively. We give an example, as shown in Part 1 of Figure 1.

Figure 1.

Creating a geometry and generating its mesh shown in Figure 2.

First with createpde(1) we create a raw model considering one partial differential equation. Then we import or create geometry for this model. For our purposes 3 stacked cylinders with radius 2 and heights 1, 3 and 2 would be fine as shown in Figure 2. pdegplot shows how Matlab has specified regions (cylinders) and faces by numbering them. Then generateMesh generates a mesh with linear (or quadratic by default) tetrahedron elements and can be viewed by pdeplot3D. Smaller value for Hmax gives a finer mesh. Setting $N_n$ and $N_e$ for number of nodes and elements, respectively, Nodes that stores 3 Cartesian coordinates of each node is a matrix of size $3\times N_n$. Since all elements are linear tetrahedron, Cells that stores index of nodes of each element, would be a $4\times N_e$ matrix. For example the first element of the mesh is formed by 4 nodes, stored in Cells(:,1) and Nodes(:,Cells(:,1)) returns 3D coordinates of that 4 nodes. So we can simply traverse over nodes and elements by their index. Moreover, we can extract element indices of a region with findElements or node indices of a region or a face by findNodes, where are necessary to set boundary conditions. For example we find node indices of bottom and top faces of the domain in Lines 21–22, since we will set boundary conditions on them.

Figure 2.

A 3D domain consists of 3 stacked cylinders (left) and the generated mesh (middle) and the solution of Eq. (21) in right.

3. Weak form

In a bounded domain $\mathrm{\Omega }$ with Lipschitz boundary $\mathrm{\Gamma }$, Green’s formula says for two regular functions $u$ and $v$ we have

$${\int }_{\mathrm{\Omega }}\nabla \cdot \left(-\kappa \nabla u\right)v={\int }_{\mathrm{\Omega }}\kappa \nabla u\cdot \nabla v-{\oint }_{\mathrm{\Gamma }}\nu \cdot \kappa \nabla uv$$E3

where $\nu$ is the (outward) normal vector. Regularity means that $u$ and $v$ have piecewise continuous (at least) first order partial derivatives to make the above integrals meaningful. So multiplying both side of (1) by $v$ and applying Green’s formula we obtain

$${\int }_{\mathrm{\Omega }}\kappa \nabla u\cdot \nabla v-{\oint }_{\mathrm{\Gamma }}\nu \cdot \kappa \nabla uv+{\int }_{\mathrm{\Omega }}\gamma \cdot \nabla uv+{\int }_{\mathrm{\Omega }}\mu uv={\int }_{\mathrm{\Omega }}fv.$$E4

Approximation of functions in FEM is done by means of basis functions. Assume ${\left\{{\phi }_j\right\}}_j$, $j=1,\dots ,N_n$ is a set of basis functions, introduced in Section 4, that we can write

$$u=\sum _{j=1}^{N_n}{u}_j{\phi }_j,$$E5

and our goal is to find unknown coefficients ${\left\{u_j\right\}}_j$. Substituting (5) in (4) and setting $v={\phi }_i$, $i=1,\cdots ,N_n,$ which is called standard Galerkin method we obtain

$$\sum _{j=1}^{N_n}{u}_j{\int }_{\mathrm{\Omega }}\kappa \nabla {\phi }_j\cdot \nabla {\phi }_i+\sum _{j=1}^{N_n}{u}_j{\int }_{\mathrm{\Omega }}\gamma \cdot \nabla {\phi }_j{\phi }_i+\sum _{j=1}^{N_n}{u}_j{\int }_{\mathrm{\Omega }}{\mathit{\mu \phi }}_j{\phi }_i+{\oint }_{\mathrm{\Gamma }}\nu \cdot \left(-\kappa \nabla u\right){\phi }_i={\int }_{\mathrm{\Omega }}f{\phi }_i$$E6

Eq. (6) consists of evaluation of (from left to right) diffusion, advection and reaction terms, boundary integral and right hand side that we calculate them in elements and then sum them up $\left({\int }_{\mathrm{\Omega }}={\sum }_e{\int }_e\right)$ to assemble the linear system (2). $u$ or its gradient in boundary integral are known, hence we do not write its expansion form.

4. Shape functions

In practice definite integrals of basis functions and their derivatives in (6) are evaluated in an element, hence restriction of basis functions over elements, called shape functions, contribute in computations. With help of affine mapping we map a typical element into a fixed element, called reference element and evaluate integrals in it. So only shape functions in reference element determine values of integrals in (6).

We use Lagrange elements which means corresponding to each node of an element a shape function is defined such that it is a polynomial in the element, takes value 1 at that specific node and 0 at other nodes of the element. Therefore a basis function in (5) becomes a continues function with value 1 at a specific node and 0 at neighboring nodes with trivial (zero) extension in the rest of the domain, as plotted in Figure 3. Note that the derivatives of Lagrange basis functions are not continues on the boundary of the elements.

Figure 3.

Left) a shape function of reference bilinear quadrilateral. Right) 2 basis functions where their restrictions in a triangle is linear and in a quadrilateral is bilinear.

Some examples of reference elements and their shape functions are presented in Figure 4. We see that in linear triangle element where nodes are the vertices of the triangle, the first shape function $S_1\left(r,s\right)=1-r-s$ corresponding to the first node $n_1=\left(0,0\right)$ has value 1 at $n_1$ and 0 at $n_2=\left(1,0\right)$ and $n_3=\left(0,1\right)$. It is called a linear element since shape functions are combination of first order polynomials $\left\{1,r,s\right\}$. Starting from $n_1$ nodes are numbered counterclockwise and we consider it for all Lagrange elements.

Figure 4.

Some of reference Lagrange elements and their shape functions.

For Lagrange elements it is easy to find shape functions. For example to find $S_1$ for bilinear quadrilateral element where nodes are vertices of the square, it suffices to consider the line passing $n_2$ and $n_3$ which is $1-r$ multiplied by the line passing $n_3$ and $n_4$ which is $1-s$. So $S_1$ has to be of the form $\alpha \left(1-r\right)\left(1-s\right)$ which gives 0 for all nodes except $n_1=\left(-1,-1\right)$. Substituting $n_1$ in equation of $S_1$, we can find $\alpha$ such that it gives 1 for $n_1$. This element is called bilinear since its shape functions are linear combination of $\left\{1,r,s,\mathit{rs}\right\}$.

For quadratic triangle element where vertices and midpoints form nodes, for shape function $S_1$ corresponding to $n_1$, we need the multiplication of two lines passing $n_2$ and $n_6$ which is $0.5-r-s$ and $n_3,n_4$ and $n_5$ which is $1-r-s$. So $S_1$ has to be of the form $\alpha \left(0.5-r-s\right)\left(1-r-s\right)$. Substituting $n_1$ in equation of $S_1$, we can find $\alpha$ such that it gives 1 for $n_1$. This element is called quadratic since its shape functions are linear combination of $\left\{1,r,s,\mathit{rs},r^2,s^2\right\}$. Note that we first numbered nodes at vertices and then at edge midpoints.

Exercise 1 Find shape functions of linear tetrahedron as presented in Figure 4.

Exercise 2 Quadratic tetrahedron element has nodes at midpoint of edges of linear tetrahedron element, so it has 6 other nodes (numbered 5 to 10) in addition to vertices (numbered 1 to 4 similar to linear tetrahedron). Find the 10 shape functions. For example for node at origin we have $S_1=2\left(1-r-s-t\right)\left(0.5-r-s-t\right)$ or for node at $\left(\mathrm{0,0,1}\right)$ we have $S_4=2t\left(t-0.5\right)$.

4.1 Element integrals

Definite integrals of shape functions and their derivatives in reference elements, called element integrals, play an important role in assembly of stiffness matrix. Denoting the reference element by $\widehat{e}$ which has $N_{\widehat{e}}$ nodes and partial derivatives of shape functions by ${\partial }_{r}S$ which means $\frac{\partial S}{\partial r}$ we define the matrix $D_{\mathit{rr}}$ (for diffusion term) such that its $\mathit{ij}$th entry is

$$D_{\mathit{rr}}^{\mathit{ij}}={\int }_{\widehat{e}}{\partial }_r{S}_i{\partial }_r{S}_j,i,j=1,\cdots ,N_{\widehat{e}}.$$E7

Similarly we define $D_{\mathit{rs}}$ and $D_{\mathit{ss}}$. In 3D space we also have to define $D_{\mathit{rt}},D_{\mathit{st}}$ and $D_{\mathit{tt}}$. In Figure 5 we show how to evaluate $D_{\mathit{rs}}$ for linear tetrahedron element. For example $D_{\mathit{rs}}$ and $D_{\mathit{st}}$ for linear tetrahedron are

Figure 5.

Element integrals ${\mathrm{D}}_{\mathit{rs}}$ for linear tetrahedron.

$$D_{\mathit{rs}}=\frac{1}{6}\left(\begin{array}{cccc}1& 0& -1& 0\\ -1& 0& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),D_{\mathit{st}}=\frac{1}{6}\left(\begin{array}{cccc}1& 0& 0& -1\\ 0& 0& 0& 0\\ -1& 0& 0& 1\\ 0& 0& 0& 0\end{array}\right).$$

Exercise 3 For quadratic tetrahedron element find $D_{\mathit{rs}}$ as

$$D_{\mathit{rs}}=\frac{1}{30}\left(\begin{array}{cccccccccc}3& 0& 1& 0& -1& -4& -1& 1& 0& 1\\ 1& 0& -1& 0& -3& 0& 1& 3& 0& -1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -4& 0& 0& 0& 4& 4& 0& -4& 0& 0\\ -1& 0& -3& 0& 4& 4& 4& -4& 0& -4\\ -1& 0& 1& 0& 4& 0& 8& -4& 0& -8\\ 1& 0& 3& 0& -4& -4& -4& 4& 0& 4\\ 1& 0& -1& 0& -4& 0& -8& 4& 0& 8\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$$

Other element integrals such as $D_r$ (for advection term) and $D$ (for reaction or right hand side terms) can be defined where their $\mathit{ij}$th entries are

$$D_r^{\mathit{ij}}={\int }_{\widehat{e}}{\partial }_r{S}_i{S}_j,D^{\mathit{ij}}={\int }_{\widehat{e}}{S}_i{S}_j,i,j=1,\cdots ,N_{\widehat{e}}.$$E8

Note that $D_{\mathit{rr}}$ and $D$ are symmetric matrices while $D_r$ is not. Moreover, since $D_{\mathit{rs}}=D_{\mathit{sr}}$ we calculate only one of them. All element integrals are evaluated once and saved for future use. In Figure 1 we load element integrals of linear tetrahedron that we have calculated and saved in LinearTetraElementIntegral.mat.

4.2 Affine mapping on reference elements

So far we introduced shape functions and calculated element integrals, all in reference element. Now we explain how to map an arbitrary element in physical space into the reference element and evaluate integrals in (6) only in terms of element integrals and without any numerical integration. Clearly change of variable is necessary to evaluate integrals when we use a mapping which is explained in next subsection.

Consider an arbitrary triangle, say $e$, in 2D space with vertices $v_i=\left(x_i,y_i\right)$, $i=\mathrm{1,2,3},$ and the reference linear triangle, $\widehat{e}$. Set

$$T=\left(\begin{array}{cc}{x}_2-x_1& x_3-x_1\\ y_2-y_1& y_3-y_1\end{array}\right),\eta =\left(\begin{array}{c}x\\ y\end{array}\right),\xi =\left(\begin{array}{c}r\\ s\end{array}\right),{\eta }_o=\left(\begin{array}{c}{x}_1\\ y_1\end{array}\right)$$E9

and define the affine (composition of a linear mapping and a translation) mapping $\widehat{e}\to e$ with the rule

$$\xi \to \eta =T\xi +{\eta }_o.$$E10

We see that the affine mapping (10) maps the nodes of reference linear triangle, ${\left\{n_i\right\}}_{i=\mathrm{1,2,3}}$, to the vertices of the triangle,

$$v_i=T{n}_i+v_1,i=\mathrm{1,2,3}.$$E11

Exercise 4 Show that the affine mapping in (9) and (10) also maps the reference quadratic triangle to an arbitrary triangle with 6 nodes (3 at vertices and 3 at edges midpoint).

Exercise 5 Show that the affine mapping,

$$\xi =\left(\begin{array}{c}r\\ s\end{array}\right)\to \eta =\left(\begin{array}{c}x\\ y\end{array}\right)=\frac{1}{2}\left(\begin{array}{cc}{x}_2-x_1& x_4-x_1\\ y_2-y_1& y_4-y_1\end{array}\right)\left(\begin{array}{c}r\\ s\end{array}\right)+\frac{1}{2}\left(\begin{array}{c}{x}_1+x_3\\ y_1+y_3\end{array}\right)$$E12

maps the reference bilinear quadrilateral to an arbitrary parallelogram. Remember in a parallelogram where $v_1$ and $v_3$ are opposite vertices, we have $x_1+x_3=x_2+x_4$ and $y_1+y_3=y_2+y_4$.

Exercise 6 Show that the affine mapping,

$$\xi =\left(\begin{array}{c}r\\ s\\ t\end{array}\right)\to \eta =\left(\begin{array}{c}x\\ y\\ z\end{array}\right)=\left(\begin{array}{ccc}{x}_2-x_1& x_3-x_1& x_4-x_1\\ y_2-y_1& y_3-y_1& y_4-y_1\\ z_2-z_1& z_3-z_1& z_4-z_1\end{array}\right)\left(\begin{array}{c}r\\ s\\ t\end{array}\right)+\left(\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right)$$E13

maps the reference linear and quadratic tetrahedron to an arbitrary tetrahedron.

4.3 Change of variables in integrals

Consider a scalar function (a typical basis function) $\phi =\phi \left(\eta \right)$, where $\eta =\left(x,y\right)$ in 2D space. Gradient of $\phi$ is ${\nabla }_{\eta }\phi ={\left({\partial }_x\phi ,{\partial }_y\phi \right)}^t$, where $t$ denotes the transpose of the vector. The chain rule and applying the affine mapping (10) on $\phi$ give

$${\nabla }_{\eta }\phi =T^{-t}{\nabla }_{\xi }\widehat{\phi }$$E14

where $\widehat{\phi }\left(\xi \right)=\phi \left(\eta \right)$ and $T^{-t}$ is the transpose of the inverse of $T$. Note that if $\phi$ is the restriction of a basis function in element $e$, then $\widehat{\phi }$ would be the corresponding shape function in $\widehat{e}$. Recalling $\mathit{d\eta }=\mid T\mid \mathit{d\xi }$, we can write

$${\int }_e{\kappa }_e{\nabla }_{\eta }{\phi }_i\cdot {\nabla }_{\eta }{\phi }_j\mathit{d\eta }={\int }_e{\left({\kappa }_e{\nabla }_{\eta }{\phi }_i\right)}^t{\nabla }_{\eta }{\phi }_j\mathit{d\eta }={\int }_{\widehat{e}}{\left({\nabla }_{\xi }{S}_i\right)}^t{T}^{-1}{\kappa }_e^t{T}^{-t}{\nabla }_{\xi }{S}_j\mid T\mid d\xi .$$E15

Setting matrix $M_e=\mid T\mid T^{-1}{\kappa }_e^t{T}^{-t}$ and dropping subscript $e$ as well as $\eta$ and $\xi$ from $\nabla$, then we have

$$\begin{array}{l}{\int }_e\kappa \nabla {\phi }_i\cdot \nabla {\phi }_j\mathit{d\eta }=M_{11}{\int }_{\widehat{e}}{\partial }_r{S}_i{\partial }_r{S}_j\mathit{d\xi }\\ +\left(M_{12}+M_{21}\right){\int }_{\widehat{e}}{\partial }_r{S}_i{\partial }_s{S}_j\mathit{d\xi }+M_{22}{\int }_{\widehat{e}}{\partial }_s{S}_i{\partial }_s{S}_j\mathit{d\xi }\\ =M_{11}{D}_{\mathit{rr}}+\left(M_{12}+M_{21}\right)D_{\mathit{rs}}+M_{22}{D}_{\mathit{ss}}\end{array}$$E16

Exercise 7 Show that in 3D space,

$$\begin{array}{l}{\int }_e\kappa \nabla {\phi }_i\cdot \nabla {\phi }_j\mathit{d\eta }=M_{11}{D}_{\mathit{rr}}+\left(M_{12}+M_{21}\right)D_{\mathit{rs}}+\left(M_{13}+M_{31}\right)D_{\mathit{rt}}\\ +M_{22}{D}_{\mathit{ss}}+\left(M_{23}+M_{32}\right)D_{\mathit{st}}\\ +M_{33}{D}_{\mathit{tt}},\end{array}$$E17

where $M=\mid T\mid T^{-1}{\kappa }^t{T}^{-t}$.

Exercise 8 Referring to notation (8) show that

$$1.{\int }_e{\phi }_i{\phi }_j\mathit{d\eta }=\mid T\mid D$$E18

2. for a vector $\gamma$, defined for each element $e$ in 3D space

$${\int }_e\gamma \cdot \nabla {\phi }_i{\phi }_j\mathit{d\eta }=M_1{D}_r+M_2{D}_s+M_3{D}_t$$E19

where vector $M=\mid T\mid T^{-1}\gamma .$

5. Assembly of the linear system

With help of element integrals, affine mapping and Eqs. (17)–(19) we can accumulate each term of (6) into the stiffness matrix of the linear system (2).

5.1 Diffusion term

Equation

$${\int }_{\mathrm{\Omega }}\kappa \nabla {\phi }_j\cdot \nabla {\phi }_i=\sum _e{\int }_e{\kappa }_e\nabla {\phi }_j\cdot \nabla {\phi }_i$$E20

suggests how to accumulate ${\int }_{\mathrm{\Omega }}\kappa \nabla u\cdot \nabla v$ into the stiffness matrix:

1. traverse elements and in each element map ${\int }_e{\kappa }_e\nabla {\phi }_i\cdot \nabla {\phi }_j$ into the reference element

2. with help of element integrals and Eq. (17) evaluate the desired term

3. map the local matrix into the stiffness matrix.

Exercise 9 Follow the above steps to accumulate advection and reaction terms into the stiffness matrix.

In Figure 6 we implemented accumulation of diffusion and reaction terms for linear or quadratic tetrahedron. Note that other elements only have a different affine mapping $T$ and dimension, if we had elements in 2D space. We assume that $\kappa$ is a diagonal matrix for each element, so a matrix storage of size $3\times N_c$ can store $\kappa$ of all elements. We also assume that reaction coefficient $\mu$ is constant, otherwise a vector of size $N_c$ can store $\mu$ for all elements and should be an input argument. The evaluated integrals are saved such that the indices of row and column and value of the integral are set in IM, JM and DDvec (for diffusion term), respectively. So Line 33 of Figure 1 is to assembly the stiffness matrix of the equation

Figure 6.

Assembly of diffusion–reaction terms for (linear or quadratic) tetrahedron.

$${\int }_{\mathrm{\Omega }}\kappa \nabla u\cdot \nabla v+0.1{\int }_{\mathrm{\Omega }}uv=0$$E21

where $\kappa$ is set in Lines 25–26. We also set boundary conditions at bottom and top faces of the boundary (with values 1 and 10, respectively) and explain how to implement it in next subsection to obtain the solution of the problem, plotted in Figure 2.

5.2 Boundary integral

Probably correct implementation of boundary integrals is the most important part of the FEM, since boundary conditions mainly determine the situation of physical problem. Neumann, Robin and Dirichlet are 3 types of boundary conditions that might be considered on different parts of the boundary. Although boundary conditions are set on the boundary, they only affect on boundary nodes.

Neumann boundary condition.$\nu \cdot \left(-\kappa \nabla u\right)=g$ is a Neumann boundary condition, so boundary integral in (6) becomes $g\oint {\phi }_i$. If $g$ is zero, then nothing has to be done for the boundary integrals. In fact incorporating only diffusion–advection-reaction terms and right hand side function means that we have considered a pure Neumann problem and this is exactly what we do in practice. After that we add requested boundary conditions to the linear system of (2). For example if $g$ is non-zero on part of the boundary, say ${\mathrm{\Gamma }}_N$, then $-g\oint {\phi }_i$ is known for indices that lie on ${\mathrm{\Gamma }}_N$ and should be added to the $i$th entry of $b$ in (2).

Dirichlet boundary condition.$u=g$ is a Dirichlet boundary condition, which means value of $u$ at some nodes is known. So if we decompose node indices into 2 sets, say $U$ and $K$, then we can write

$$u=\sum _{j\in U}{u}_j{\phi }_j+\sum _{j\in K}{u}_j{\phi }_j,$$E22

where $U$ and $K$ include indices of unknown and known values of $u$, respectively. If with permutation vector $\left[U,K\right]$ we reorder rows and columns of the linear system in (2), then we can write

$$\mathcal{A}u=\left(\begin{array}{cc}{\mathcal{A}}_{11}& {\mathcal{A}}_{12}\\ {\mathcal{A}}_{21}& {\mathcal{A}}_{22}\end{array}\right)\left(\begin{array}{c}{u}_U\\ u_K\end{array}\right)=\left(\begin{array}{c}{b}_U\\ b_K\end{array}\right)$$E23

The first line gives a smaller linear system

$${\mathcal{A}}_{11}{u}_U=b_U-{\mathcal{A}}_{12}{u}_K$$E24

which preserves the symmetry or positive definiteness of the original problem. In Lines 35–42 of Figure 1 we showed how to implement it in Matlab. Setting values 1 and 10 for bottom and top faces, we expect a smooth solution starting from 1 at the bottom and reaching to 10 at the top as shown in Figure 2.

The second approach to set Dirichlet condition that does not use permutation of the linear system is setting ${\mathcal{A}}_{21}$ and ${\mathcal{A}}_{22}$ in (23) to zero and identity matrices, receptively and solving the modified linear system. Although this approach gives the correct solution theoretically, it is very error-prone numerically as shown in a test case in Figure 7 and explained as follows.

Figure 7.

Left) Result of correct implementation of Dirchlet BC by using Eq. (24). Right) nonphysical (wrong) result due to using the modified linear system (the second approach to implement Dirchlet BC).

A domain with cuboid elements is refined in the middle and along the $x$-axis to simulate a fracture in the domain. We solve $\nabla \cdot \left(-\kappa \nabla u\right)=0$ where $\kappa$ in fracture is 100 times larger than rest of the domain. Setting ${10}^6$ and $5\times {10}^6$ as Dirichlet boundary condition on left and right faces (along $x$-axis) of the domain, respectively, a correct solution should be started from ${10}^6$ and monotonically reached to $5\times {10}^6$. The linear system is symmetric positive definite and preconditioned conjugate gradient method is employed to solve it. The correct solution shown in left image of Figure 7 is obtained by (24). However, the second approach gives the nonphysical (wrong) solution along the fracture as shown in right. Moreover, the correct solution is obtained 3 times faster than the wrong one, mainly because the second approach violates the symmetry (but preserves the positive definiteness) of the linear system, due to setting $A_{21}=0$.

On the other hand reduction in size of the linear system by using (24), has great advantage in applications that large number of nodes have Dirichlet (more precisely essential) boundary condition. For example, in pore scale modeling and to approximate the permeability of a sandstone model, we had to set Dirichlet boundary condition for $75\%$ of the nodes. Since the model included almost 10 million elements, solving a linear system with 2.5 million unknowns gave a significant speed up in computational time, in addition to significant improvement in accuracy.

Exercise 10 Robin boundary condition is a combination of Neumann and Dirichlet conditions which states $u-\mathit{\beta \nu }\cdot \kappa \nabla u=g$ on ${\mathrm{\Gamma }}_R$. $\beta$ is a non-zero function (normally a constant number) on ${\mathrm{\Gamma }}_R$. Considering $-\nu \cdot \kappa \nabla u=\left(g-u\right)/\beta$, explain how to implement it.

5.3 Right hand side term

Function $f$ in (1) is the source or sink term and hence normally is defined on a very small part of the domain. It can be defined as a constant number over an element or can be approximated by its nodal values in the form $f=\sum f_j{\phi }_j$ and therefore right hand side term becomes ${\sum }_{j=1}^{N_n}{f}_j\int {\phi }_j{\phi }_i$ with possibly too many $f_j=0$. We prefer nodal value approximation of $f$, since in applications that need to solve parabolic equation, $f$ has a value in all elements; see $q^{\left[k\right]}$ in Eqs. (27) and (28).

6. Generalization to parabolic equation

In applications such as simulation of compressible flows we need to solve the time dependent equation

$$\frac{\partial \rho \left(u\right)}{\partial t}+\nabla \cdot \left(-\rho \left(u\right)\kappa \nabla u\right)=f\left(u\right),\text{in}\mathrm{\Omega }$$E25

where $\rho$ and $f$ are nonlinear functions of $u$ and boundary condition and initial value are provided. A fully implicit time discretization of (25) gives

$$\frac{{\rho }^{n+1}-{\rho }^n}{\mathrm{\Delta }t}+\nabla \cdot \left(-{\rho }^{n+1}\kappa \nabla u^{n+1}\right)=f^{n+1}$$

which simplifies to

$$G\left(u\right)=\frac{1}{\mathrm{\Delta }t}\rho +\nabla \cdot \left(-\mathit{\rho \kappa }\nabla u\right)-f\left(u\right)-\frac{1}{\mathrm{\Delta }t}{\rho }^n=0.$$E26

In Eq. (26), the unknown $u$-dependent variables $\rho$ and $f$ should be evaluated at the current time $t^{n+1}$, while ${\rho }^n$ is known from the previous time step.

Newton–Raphson linearization is the most common method to solve the nonlinear Eq. (26) which generates a sequence $\left\{u^{\left[k\right]}\right\}$, $k=\mathrm{0,1,2},\dots$ by solving the following equation

$$\nabla \cdot \left(-{\rho }^{\left[k\right]}\kappa \nabla u\right)+{\mu }^{\left[k\right]}u=q^{\left[k\right]},$$E27

where

$${\mu }^{\left[k\right]}=\frac{1}{\mathrm{\Delta }t}{\left.\frac{\mathit{d\rho }}{\mathit{du}}\right|}_{u^{\left[k\right]}}-{\left.\frac{\mathit{df}}{\mathit{du}}\right|}_{u^{\left[k\right]}},q^{\left[k\right]}={\mu }^{\left[k\right]}{u}^{\left[k\right]}+f^{\left[k\right]}-\frac{1}{\mathrm{\Delta }t}{\rho }^n.$$E28

Having bounded derivatives and a good starting point $u^{\left[0\right]}$ (which is normally chosen $u^n$), the resulted sequence $u^{\left[k\right]}$ in Eq. (27) converges quadratically to $u^{n+1}$, the solution of Eq. (26). Note that (26) is a diffusion–reaction equation which is completely explained how to solve and (core) codes were provided, as well.

7. Conclusion and future works

In this introductory chapter we presented the framework of FEM in brief (but effective) and implemented a numerical solver for diffusion–advection-reaction equation. Accumulation of different terms and setting boundary conditions correctly as well as evaluating of definite integrals without numerical integration were explained and their codes were also given. Finally we showed that nonlinear parabolic equations can be solved by combining of Newton method and diffusion–reaction equation where the latter was explained in this chapter.

However, to have a reliable and efficient simulator several topics and challenges should be addressed. Assuming correct mathematical model, mesh generation and assigning (measured or suggested) values to the coefficients of the problem are very important to make close the numerical model and its solution to the real problem [4, 5]. In a real problem, 2D and 3D elements appear together. For example 2D elements are employed to model faults in a 3D reservoir. Moreover, the generated mesh is normally unstructured and different types of elements are included in the mesh. Hence to assembly the linear system, traversing elements and nodes should have been implemented very efficiently.

Solving linear systems is normally the most time consuming part of the simulation and efficient implementation of linear solvers particularly in parallel machines, are of great interest [6]. Multigrid and multiscale methods, particularly their algebraic form, have attracted interest to solve large scale linear systems since they have shown good scalability in addition to resolving low frequency parts of the error in solving linear systems [7, 8]. Standard Galerkin method that we used in this chapter is not a conservative scheme hence modifications or other methods such as discontinues Galerkin method would be necessary to solve problems in computational fluid dynamics [9, 10]. Proposing and implementing of advanced numerical algorithms to linearize and solve a system of nonlinear coupled initial-boundary value problem in a large scale domain become necessary. At last comparison with real data (if available) and quantifying uncertainties might force us to revisit all steps of the simulation again.

Acknowledgments

Hani Akbari would like to thank Prof. Lars K. Nielsen, scientific director at Novo Nordisk Foundation Center for Bio-Sustainability. Hani Akbari is grateful to DTU-Biosustain since this work was completed when he was a postdoc at DTU-Biosustain.