Open access peer-reviewed chapter

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

By Hani Akbari

Submitted: November 27th 2020Reviewed: May 6th 2021Published: June 3rd 2021

DOI: 10.5772/intechopen.98291

Downloaded: 62


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.


  • finite element method
  • diffusion–advection-reaction equation
  • boundary condition
  • implementation

1. Introduction

In a domain Ωin Rn,n=1,2,3,consider a partial differential equation of the form of


where the unknown function u, reaction coefficient μand the given function fare scalar functions ΩR. Diffusion coefficient κand the (convergence free) advection coefficient γare defined in Ωand give (normally piecewise constant) values in Rn×nand Rn, 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, CO2 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


where Aand bare 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 Ω, 2) heterogeneity and high contrast in coefficient κ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 κ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 κ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 andNodes(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×2matrix in 2D or a 3×3matrix in 3D), named κe, will be set to form κ. In each element e, we can also set μeand γeto form μand γ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 inFigure 2.

First withcreatepde(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.pdegplotshows how Matlab has specified regions (cylinders) and faces by numbering them. ThengenerateMeshgenerates a mesh with linear (or quadratic by default) tetrahedron elements and can be viewed bypdeplot3D. Smaller value forHmaxgives a finer mesh. Setting Nnand Nefor number of nodes and elements, respectively,Nodesthat stores 3 Cartesian coordinates of each node is a matrix of size 3×Nn. Since all elements are linear tetrahedron,Cellsthat stores index of nodes of each element, would be a 4×Nematrix. For example the first element of the mesh is formed by 4 nodes, stored inCells(:,1)andNodes(:,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 withfindElementsor node indices of a region or a face byfindNodes, 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 ofEq. (21)in right.


3. Weak form

In a bounded domain Ωwith Lipschitz boundary Γ, Green’s formula says for two regular functions uand vwe have


where νis the (outward) normal vector. Regularity means that uand vhave piecewise continuous (at least) first order partial derivatives to make the above integrals meaningful. So multiplying both side of (1) by vand applying Green’s formula we obtain


Approximation of functions in FEM is done by means of basis functions. Assume φjj, j=1,,Nnis a set of basis functions, introduced in Section 4, that we can write


and our goal is to find unknown coefficients ujj. Substituting (5) in (4) and setting v=φi, i=1,,Nn,which is called standard Galerkin method we obtain


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 Ω=eeto assemble the linear system (2). uor 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 S1rs=1rscorresponding to the first node n1=00has value 1 at n1and 0 at n2=10and n3=01. It is called a linear element since shape functions are combination of first order polynomials 1rs. Starting from n1nodes 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 S1for bilinear quadrilateral element where nodes are vertices of the square, it suffices to consider the line passing n2and n3which is 1rmultiplied by the line passing n3and n4which is 1s. So S1has to be of the form α1r1swhich gives 0 for all nodes except n1=11. Substituting n1in equation of S1, we can find αsuch that it gives 1 for n1. This element is called bilinear since its shape functions are linear combination of 1rsrs.

For quadratic triangle element where vertices and midpoints form nodes, for shape function S1corresponding to n1, we need the multiplication of two lines passing n2and n6which is 0.5rsand n3,n4and n5which is 1rs. So S1has to be of the form α0.5rs1rs. Substituting n1in equation of S1, we can find αsuch that it gives 1 for n1. This element is called quadratic since its shape functions are linear combination of 1rsrsr2s2. Note that we first numbered nodes at vertices and then at edge midpoints.

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

Exercise 2Quadratic 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 S1=21rst0.5rstor for node at 0,0,1we have S4=2tt0.5.

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 êwhich has Nênodes and partial derivatives of shape functions by rSwhich means Srwe define the matrix Drr(for diffusion term) such that its ijth entry is


Similarly we define Drsand Dss. In 3D space we also have to define Drt,Dstand Dtt. In Figure 5 we show how to evaluate Drsfor linear tetrahedron element. For example Drsand Dstfor linear tetrahedron are

Figure 5.

Element integralsDrsfor linear tetrahedron.


Exercise 3For quadratic tetrahedron element find Drsas


Other element integrals such as Dr(for advection term) and D(for reaction or right hand side terms) can be defined where their ijth entries are


Note that Drrand Dare symmetric matrices while Dris not. Moreover, since Drs=Dsrwe 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 inLinearTetraElementIntegral.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 vi=xiyi, i=1,2,3,and the reference linear triangle, ê. Set


and define the affine (composition of a linear mapping and a translation) mapping êewith the rule


We see that the affine mapping (10) maps the nodes of reference linear triangle, nii=1,2,3, to the vertices of the triangle,


Exercise 4Show 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 5Show that the affine mapping,


maps the reference bilinear quadrilateral to an arbitrary parallelogram. Remember in a parallelogram where v1and v3are opposite vertices, we have x1+x3=x2+x4and y1+y3=y2+y4.

Exercise 6Show that the affine mapping,


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) φ=φη, where η=xyin 2D space. Gradient of φis ηφ=xφyφt, where tdenotes the transpose of the vector. The chain rule and applying the affine mapping (10) on φgive


where φ̂ξ=φηand Ttis the transpose of the inverse of T. Note that if φis the restriction of a basis function in element e, then φ̂would be the corresponding shape function in ê. Recalling =T, we can write


Setting matrix Me=TT1κetTtand dropping subscript eas well as ηand ξfrom , then we have


Exercise 7Show that in 3D space,


where M=TT1κtTt.

Exercise 8Referring to notation (8) show that


2. for a vector γ, defined for each element ein 3D space


where vector M=TT1γ.


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



suggests how to accumulate Ωκuvinto the stiffness matrix:

  1. traverse elements and in each element map eκeφiφjinto 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 9Follow 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 Tand dimension, if we had elements in 2D space. We assume that κis a diagonal matrix for each element, so a matrix storage of size 3×Nccan store κof all elements. We also assume that reaction coefficient μis constant, otherwise a vector of size Nccan store μ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 inIM,JMandDDvec(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.


where κ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.νκu=gis a Neumann boundary condition, so boundary integral in (6) becomes gφi. If gis 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 gis non-zero on part of the boundary, say ΓN, then gφiis known for indices that lie on ΓNand should be added to the ith entry of bin (2).

Dirichlet boundary condition.u=gis a Dirichlet boundary condition, which means value of uat some nodes is known. So if we decompose node indices into 2 sets, say Uand K, then we can write


where Uand Kinclude indices of unknown and known values of u, respectively. If with permutation vector UKwe reorder rows and columns of the linear system in (2), then we can write


The first line gives a smaller linear system


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 A21and A22in (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 usingEq. (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 κu=0where κin fracture is 100 times larger than rest of the domain. Setting 106and 5×106as Dirichlet boundary condition on left and right faces (along x-axis) of the domain, respectively, a correct solution should be started from 106and monotonically reached to 5×106. 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 A21=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 10Robin boundary condition is a combination of Neumann and Dirichlet conditions which states uβνκu=gon ΓR. βis a non-zero function (normally a constant number) on ΓR. Considering νκu=gu/β, explain how to implement it.

5.3 Right hand side term

Function fin (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=fjφjand therefore right hand side term becomes j=1Nnfjφjφiwith possibly too many fj=0. We prefer nodal value approximation of f, since in applications that need to solve parabolic equation, fhas a value in all elements; see qkin 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


where ρand fare nonlinear functions of uand boundary condition and initial value are provided. A fully implicit time discretization of (25) gives


which simplifies to


In Eq. (26), the unknown u-dependent variables ρand fshould be evaluated at the current time tn+1, while ρnis known from the previous time step.

Newton–Raphson linearization is the most common method to solve the nonlinear Eq. (26) which generates a sequence uk, k=0,1,2,by solving the following equation




Having bounded derivatives and a good starting point u0(which is normally chosen un), the resulted sequence ukin Eq. (27) converges quadratically to un+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.



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.

© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Hani Akbari (June 3rd 2021). A Numerical Simulator Based on Finite Element Method for Diffusion-Advection-Reaction Equation in High Contrast Domains, Recent Advances in Numerical Simulations, Francisco Bulnes and Jan Peter Hessling, IntechOpen, DOI: 10.5772/intechopen.98291. Available from:

chapter statistics

62total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

A Modified Spectral Relaxation Method for Some Emden-Fowler Equations

By Gerald Tendayi Marewo

Related Book

First chapter

Introductory Chapter: Advanced Communication and Nano-Processing of Quantum Signals

By Francisco Bulnes

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us