Latest revision as of 14:43, 6 July 2021

Abstract

This paper proposes an efficient and robust algorithm for solving a physical orthotropy problem. The algorithm is based on choosing the most efficient restriction operator and on an incomplete LU decomposition suited for each orthotropy direction. Local Fourier Analysis (LFA) is carried out in order to increase the efficiency of the multigrid method. Pure diffusion with orthotropy aligned to the coordinate axis x is the model considered. Equations are discretized by Finite Difference Method with uniform grid and second-order numerical scheme. Problems are solved with geometric multigrid method, correction scheme, V-cycle and standard coarsening ratio. The asymptotic convergence factor is calculated for different multigrid components, such as restriction operators, prolongation operators and solvers. Based on the optimum components obtained by LFA, we carried out experiments to analyze the complexity and computational cost of the algorithm proposed. The main conclusion is that the methodology proposed is efficient for the resolution of problems with strong orthotropy.

Keywords: Physical orthotropy, multigrid components, diffusion, local Fourier analysis

1. Introduction

Computational Fluid Dynamics (CFD) is a branch of Computational Science that studies numerical methods used for simulating fluid flow problems. It is known that these methods usually have a high computational cost. In general, this happens as the problems that have to be solved require the resolution of algebraic equation systems whose coefficient matrices are large and sparse [1].

Linear systems are obtained by discretizing the mathematical model, which consists in approximating, through algebraic equations, each term of the mathematical model for each grid node or point. This process leads to an algebraic equation system of the form

where ${\displaystyle A\in {R^{NxN}}}$ is the coefficient matrix, ${\textstyle T\in {R^{N}}}$ is the variable of interest and ${\textstyle b\in {R^{N}}}$ is the independent vector.

In CFD, the methods that are traditionally employed in this process are: Finite Difference Method (FDM) [2], Finite Volume Method (FVM) [3] and Finite Element Method (FEM) [4].

The algebraic equation system given by Eq. (1) can be solved using direct or iterative methods. In this work, iterative methods were employed due to the aforementioned characteristics of linear systems in CFD (sparse and large coefficients matrices), for which iterative methods are recommended. The multigrid method is one of the most effective methods to accelerate the convergence of iterative methods when solving linear or nonlinear systems, isotropic or anisotropic problems, among others [5-7].

Anisotropic problems are fairly common in Engineering and appear in many phenomena, such as when a material has different heat conduction behaviors in different directions. In this case, the coefficients of the differential equations are distinct among themselves and generate what is called physical anisotropy. Anisotropies can also appear from discretization of grids with different spacing in each direction, for instance, boundary layer problems. This is called geometric anisotropy [5,6]. The efficiency of the multigrid method has not yet been fully achieved for problems with strong anisotropy, either physical or geometric [8].

Physical anisotropic convection problems have been investigated by Rabi and de Lemos [9], who discretized two-dimensional pure diffusion, pure advection and advection-diffusion equations by applying the finite volume method. The multigrid method was employed using correction scheme and V- and W-cycles. The authors presented a study on the different speed ranges, number of grids, number of smoothing steps at each grid level and different solvers. They concluded that there was a significant reduction in the computational effort required for increasing the values of the components of the advection velocity.

Wienands and Joppich [10] presented an in-depth study on the Local Fourier Analysis (LFA) and its application on several problems, including anisotropic problems. The authors calculated the convergence factor of the multigrid method for an anisotropic diffusion equation with different solvers and restriction and prolongation operators.

Johannsen [11] solved an anisotropic diffusion problem using finite volume method for discretizing the equations and 9-point incomplete LU (ILU) decomposition for solving the systems of linear equations. The author employed LFA to demonstrate the superior smoothing properties of ILU.

Oliveira et al. [8] evaluated geometric anisotropy for different grids and aspect ratios. They also assessed some components of the multigrid method, such as: solvers, type of restriction, number of levels and number of inner iterations, among several coarsening algorithms. They concluded that Partial Weighting (PW) had a good performance.

Vinogradova and Krukier [12] solved a three-dimensional advection-diffusion problem with intermediate anisotropy using FDM for discretizing equations and ILU for the resolution of the linear systems. They concluded that the proposed methodology is efficient, however, the coefficients of the mixed derivatives present limitations, which is a disadvantage and has no physical significance.

Vassoler-Rutz et al. [13] analyzed the effect of physical anisotropy on the multigrid method for two anisotropic diffusion problems. They used FAS scheme, V-cycle as well as Modified Strongly Implicit (MSI) and Gauss-Seidel (GS) solvers. They concluded that, for strong anisotropies, the complexity order of the multigrid method is not suitable. 

Several works on the implementation of the multigrid method found in the literature demonstrate that the choice of the multigrid components is crucial for the convergence or not of the method. Trottenberg et al. [6] state that this choice is difficult and thus small changes can considerably improve convergence. In this sense, LFA can help this choice as it allows to predict the performance of the multigrid method, since it provides estimates of the convergence rates based on the variation of the multigrid method components.

Pinto et al. [14] solved an anisotropic diffusion problem using ILU in triangular grids. They used LFA to highlight the good smoothing properties of the solver and the asymptotic convergence of the multigrid.

Franco et al. [15] performed LFA in transient problems and obtained the critical value of the parameter that represents the level of space-time anisotropy for 1D and 2D Fourier equations.

Oliveria et al. [16] solved an 2D anisotropic diffusion equation. The equation was discretized by the Finite Difference Method (FDM) and Central Differencing Scheme (CDS). Correction Scheme (CS). An xy-zebra-GS smoother was proposed, which proved to be efficient and robust for the different anisotropy coefficients. They concluded that, the convergence factors calculated empirically and by LFA are in agreement.

A particular case of anisotropy is denominated orthotropy, which happens when the anisotropy occurs in orthogonal directions. In this work, an efficient and robust method for solving physical orthotropy problems using LFA is proposed. A two-dimensional diffusion mathematical model is considered, in which physical orthotropy appears in the coefficients and it will be denominated diffusion orthotropy. Equations were discretized using FDM with second-order central difference scheme.

The asymptotic convergence factor (${\textstyle \rho _{loc}}$) of the multigrid method was calculated by assessing the ILU solver in different directions [7] as well as several restriction and prolongation operators. The results obtained via LFA were used to assess the influence of the diffusion orthotropy on the computational cost and took into account the CPU time and number of operations in each V-cycle and at the restriction step.

The remainder of this work is organized as follows: section 2 presents the mathematical and numerical models; section 3 discusses considerations regarding the LFA used; section 4 presents the results of the convergence analysis and complexity analysis; and section 5 presents the conclusion.

2. Mathematical and Numerical Models

For the problem presented below, the calculus domain used is given by ${\textstyle {0}\leq {x}\leq {1}}$, ${\textstyle {0}\leq {y}\leq {1}}$ and the discretization of the equations is done using a uniform grid with a number of points given by ${\textstyle N=N_{x}.N_{y}}$ , where ${\textstyle N_{x}}$ and ${\textstyle N_{y}}$ are the number of points in the directions ${\displaystyle x}$ and ${\displaystyle y}$, respectively, including the boundaries.

2.1. Mathematical Model and Discretization

A model that exemplifies physical and geometric anisotropy for a two-dimensional diffusion equation is given by Trottenberg et al. (2001) as shown below,

where ${\displaystyle g=\cos(\alpha )}$, ${\displaystyle w=sen\,(\alpha )}$, ${\displaystyle 0\leq \alpha \leq {\frac {\pi }{2}}}$ and ${\displaystyle 0<\varepsilon <<1}$ or ${\displaystyle \varepsilon >>1}$.

For ${\displaystyle \alpha ={\frac {\pi }{2}}}$, it is considered that the expression given by Eq. (2) is aligned with the axis of the coordinate y, so it becomes

From here, the diffusion orthotropic problem will be assessed by means of the two-dimensional diffusion equation given by Eq. (4) [5,6]

where T is the temperature, ${\textstyle T_{xx}}$ is the second derivative of T as a function of ${\textstyle x}$, ${\textstyle T_{yy}}$ is the second derivative as a function of ${\textstyle y}$ and ${\textstyle \varepsilon >0}$.

The source term S and the analytical solution T are given by

Eq. (4) was discretized using FDM with CDS, resulting in

where ${\displaystyle T}$ is the unknown of the system.

Figure 1(b) depicts the notation of the grid points in Figure 1(a). The points P (central), W (west), E (east), N (north) and S (south) in Figure 1(b) correspond to the points ${\textstyle (i,j),(i-1,j),(i+1,j),(i,j+1),(i,j-1)}$ in Figure 1(a), respectively.

Figure 1. Points of a uniform two-dimensional grid

The classic 5-point finite difference is not convergent for cases of general anisotropy, such as the example given by the full tensor

Therefore, this methodology cannot be generalized to any type of anisotropic problem. In this paper, we analyze a specific case of anisotropy given by the tensor

which represents a case of orthotropy.

The discretization of Eq. (4) results in Eq. (6), and for the inner points, considering ${\textstyle h_{x}^{}={\frac {1}{N_{x}^{}-1}}}$ and ${\textstyle h_{y}^{}={\frac {1}{N_{y}^{}-1}}}$

For the boundaries north, south, east and west,${\textstyle {{a}_{P}}=1}$, ${\displaystyle {{a}_{W}}={{a}_{E}}={{a}_{N}}={{a}_{S}}=0}$.

2.2. Multigrid Method and Computational Details

The multigrid method accelerates the convergence rate of iterative methods. It consists in employing a group of grids with different refinement levels. At each refinement level of the grid, the more oscillatory errors are smoothed, and only low frequency errors remain. After passing to another grid, the remaining low frequency errors become more oscillatory. The efficiency of this process, called smoothing, depends on the choice of a suitable solver.

For employing multigrid, besides a solver with good smoothing properties, grid transfer operators are required (restriction and prolongation).

In [8] and [17], restriction through partial weighting in the directions x and y, henceforth denoted by ${\displaystyle {\text{P}}{{\text{W}}_{\text{x}}}}$ and ${\displaystyle {\text{P}}{{\text{W}}_{\text{y}}}}$, respectively, applied in problems involving geometric anisotropy was proposed. These operators are given in stencil notation as

In this work, such restriction was combined with 7-point incomplete LU decomposition (henceforth denoted by ILU). According to [18], this decomposition has a better convergence factor than 5-point incomplete LU decomposition for orthotropic problems.

In Eq. (4), discretized using FDM, the stencil for the 5-point Laplacian operator is given by

The ${\displaystyle {\mbox{ILU}}}$ decomposition of the same operator will be given by

note that in this case, ${\textstyle b=f=0}$.

Thus, the ${\displaystyle {\mbox{ILU}}}$ decomposition is represented by:

where ${\displaystyle {{L}_{h}}}$ is the stencil of a lower triangular matrix, ${\displaystyle {{U}_{h}}}$ is the stencil of an upper triangular matrix and ${\displaystyle {{R}_{h}}}$ is the residual matrix.

The iterative process for solving Eq. (4) can be:

Depending on the ordination of the grid points, different directions can be obtained for the ${\displaystyle {\mbox{ILU}}}$ decomposition. In lexicographical order, ${\textstyle {\mbox{ILU}}_{\mbox{EN}}}$ [18], is given by

Another example of ordination for ${\displaystyle {\mbox{ILU}}}$, ${\textstyle {\mbox{ILU}}_{\mbox{NE}}}$ [18], is given by

The linear system given by Eq. (1) was solved using geometric multigrid method [5,6] with correction scheme (CS), V-cycle and zero initial guess.

The coarsening ratio is given by ${\displaystyle r=2}$ (standard coarsening) [18]. The grid transfer operators employed were: Injection Restriction (INJ), Half Weighting (HW), Full Weighting (FW) [6], Partial Weighting in ${\textstyle x}$ (${\textstyle {\mbox{PW}}_{\mbox{x}}}$ ), Partial Weighting in ${\textstyle y}$ (${\textstyle {\mbox{PW}}_{\mbox{y}}}$) as well as prolongation by bilinear interpolation and 7-point interpolation. The systems of equations obtained by means of discretization were resolved using 7-point ILU solvers in different directions (${\textstyle {\mbox{ILU}}_{\mbox{EN}}}$, ${\textstyle {\mbox{ILU}}_{\mbox{NE}}}$, among others).

The stop criterion used to interrupt the iterative process is based on the nondimensionlized residual norm. The residual of the system of algebraic equation is defined by

where ${\displaystyle T^{m}}$ is the solution of the unknown in the iteration ${\displaystyle m}$.

Considering ${\displaystyle {{L}^{m}}={{\left\|r_{}^{m}\right\|}_{1}}}$ and ${\displaystyle {{L}^{0}}={{\left\|R_{}^{0}\right\|}_{1}}}$, if ${\displaystyle {\frac {{L}^{m}}{{L}^{0}}}\leq tol}$ the iterative process is interrupted if ${\displaystyle tol={{10}^{-10}}}$.

Double precision arithmetic was used for the simulations. The numerical codes were implemented in the Fortran 2003 language, using the Intel 9.1 Visual Fortran application. All numerical results were obtained in a computer with Intel Core i7 2.66 GHz processor, 16 GB RAM and Windows 8 operating system, 64-bit version.

3. Local Fourier analysis

LFA allows to predict the performance of the multigrid method as it supplies estimates of the convergence rate of its components. Therefore, LFA becomes a powerful tool in quantitative analysis and in the research of efficient multigrid methods.

In order to perform the LFA, general discrete linear operators with constant coefficients are considered, which are defined in an infinite grid ${\displaystyle G_{h}}$ , where the influence of the boundaries can be dismissed [6].

Consider the grid functions of the form of

where ${\displaystyle x}$ varies in the infinite grid ${\displaystyle G_{h}}$ and ${\displaystyle \theta }$ is a parameter that characterize the frequency of the function concerning to grid ${\displaystyle G_{h}}$ .

For ${\displaystyle -\pi \leq \theta <\pi }$ , every function of the grid ${\displaystyle {{\phi }_{h}}\left(\theta ,x\right)}$ are eigenfunctions of a discrete operator that can be written as a stencil.

Thus,

where

and ${\displaystyle {{s}_{k}}}$ is the stencil notation of the operator, were ${\displaystyle k\in \{(-1,-1),\,(-1,0),(-1,1),...,(1,-1),(1,0),(1,1)\}}$ for  the 5-point stencil.

In order to smooth as well as to analyze the two grids, it is necessary to distinguish between components of low and high frequency of ${\displaystyle G_{h}}$ and ${\displaystyle G_{2h}}$ .

It is known that [6] only

are visible in ${\displaystyle G_{h}}$ .

For each ${\textstyle {\bar {\theta }}\in \left[-{\frac {\pi }{2}},{\frac {\pi }{2}}\right)\times \left[-{\frac {\pi }{2}},{\frac {\pi }{2}}\right)}$, three other frequency components ${\displaystyle {{\phi }_{h}}\left(\theta ,x\right)}$ with ${\textstyle \theta \in [-\pi ,\pi )\times [-\pi ,\pi )}$ coincide in ${\displaystyle G_{2h}}$ with ${\displaystyle {{\phi }_{h}}\left({\bar {\theta }},x\right)}$ and are not visible in ${\displaystyle G_{2h}}$ . Therefore, the low- and high-frequency components are defined as follows:

${\displaystyle \phi }$ is a low-frequency component${\displaystyle \Leftrightarrow \theta \in {{T}^{low}}=\left[-{\frac {\pi }{2}},{\frac {\pi }{2}}\right)\times \left[-{\frac {\pi }{2}},{\frac {\pi }{2}}\right)}$

${\displaystyle \phi }$ is a high-frequency component${\displaystyle \Leftrightarrow \theta \in {{T}^{high}}=\left[-\pi ,\pi \right)\times \left[-\pi ,\pi \right)\backslash \left[-{\frac {\pi }{2}},{\frac {\pi }{2}}\right)\times \left[-{\frac {\pi }{2}},{\frac {\pi }{2}}\right)}$ (Figure 2).

Figure 2. Low- (inner white area) and high-frequency (hatched area) areas

Considering the frequencies

the correction operator of the coarse grid is given by ${\displaystyle K_{h}^{2h}}$ ans is represented by a 4x4 matrix ${\displaystyle {\hat {K}}_{h}^{2h}}$ , as follows

where ${\displaystyle {\hat {I}}_{h}^{}}$ is represented by a 4x4 identity matrix.

${\displaystyle {{\hat {L}}_{h}}(\theta )}$ is the 4x4 matrix:

where ${\displaystyle {{\tilde {L}}_{h}}}$ are eigenvalues, evaluated by ${\displaystyle {{\tilde {L}}_{h}}(\theta )=\sum \limits _{\kappa }{{{s}_{\kappa }}{{e}^{i{{\theta }^{\alpha }}\kappa }}}}$, and ${\displaystyle {{s}_{\kappa }}\in R}$ are stencil coefficients [7].

The discrete Laplace operator given by Eq. (9) is represented by:

where ${\displaystyle {{\tilde {L}}_{1}}=(2+2\varepsilon )-2(\varepsilon \cos({{\theta }_{1}})+\cos({{\theta }_{2}}))}$, ${\displaystyle {{\tilde {L}}_{2}}=(2+2\varepsilon )-2(\varepsilon \cos({{\bar {\theta }}_{1}})+\cos({{\bar {\theta }}_{2}}))}$, ${\displaystyle {{\tilde {L}}_{3}}=(2+2\varepsilon )-2(\varepsilon \cos({{\bar {\theta }}_{1}})+\cos({{\theta }_{2}}))}$ and ${\displaystyle {{\tilde {L}}_{4}}=(2+2\varepsilon )-2(\varepsilon \cos({{\theta }_{1}})+\cos({{\bar {\theta }}_{2}}))}$.

The restriction operator ${\displaystyle {\hat {I}}_{h}^{2h}(\theta )}$ is a 1x4 matrix, and is given by:

For the ${\displaystyle {\mbox{INJ}}}$ restricion operator, ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{h}^{2h}}}}\,\left({{\theta }^{\alpha }}\right)=1}$, for the ${\displaystyle {\mbox{HW}}}$ restriction operator, ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{h}^{2h}}}}\,\left({{\theta }^{\alpha }}\right)={\frac {1}{4}}(2+\cos {{\bar {\theta }}_{1}}+\cos {{\bar {\theta }}_{2}})}$, for the ${\displaystyle {\mbox{FW}}}$ restriction operator, ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{h}^{2h}}}}\,\left({{\theta }^{\alpha }}\right)={\frac {1}{4}}(1+\cos {{\bar {\theta }}_{1}})(1+\cos {{\bar {\theta }}_{2}})}$, for the ${\textstyle {\mbox{PW}}_{\mbox{x}}}$ (see Eq.(8)), ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{h}^{2h}}}}\,\left({{\theta }^{\alpha }}\right)={\frac {1}{2}}(1+\cos {{\bar {\theta }}_{1}})}$ and for the ${\textstyle {\mbox{PW}}_{\mbox{y}}}$ (see Eq.(8)), ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{h}^{2h}}}}\,\left({{\theta }^{\alpha }}\right)={\frac {1}{2}}(1+\cos {{\bar {\theta }}_{2}})}$.

The prolongation operator ${\displaystyle {\hat {I}}_{2h}^{h}(\theta )}$ is a 4x1 matrix, and is given by:

For the bilinear prolongation operator, ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{2h}^{h}}}}\,\left({{\theta }^{\alpha }}\right)=(1+\cos {{\bar {\theta }}_{1}})(1+\cos {{\bar {\theta }}_{2}})}$ and for the 7-point prolongation operator, ${\displaystyle {\overset {}{\mathop {{\tilde {I}}_{2h}^{h}}}}\,\left({{\theta }^{\alpha }}\right)=\left(1+\cos {{\bar {\theta }}_{1}}+\cos {{\bar {\theta }}_{2}}+\cos({{\bar {\theta }}_{1}}-{{\bar {\theta }}_{2}})\right)}$ [10].

The operator of the grid coarse ${\displaystyle {{\left({{\hat {L}}_{2h}}(2\theta )\right)}^{-1}}}$ is a 1x1 matrix, and for the discrete Laplace operators, ${\displaystyle {{\tilde {L}}_{2h}}}$ is represented by ${\displaystyle {{\tilde {L}}_{2h}}\,(2\theta )=\sum \limits _{k}{{{s}_{k,2h}}{{e}^{i2\,\theta \,k}}}}$or ${\displaystyle {{\tilde {L}}_{2h}}={\frac {(2+2\varepsilon )-2(\varepsilon \cos(2{{\theta }_{1}})+\cos(2{{\theta }_{2}}))}{2{{h}^{2}}}}}$.

A representation for the operator of two grids ${\textstyle M_{h}^{2h}}$ can be obtained by a matrix ${\displaystyle {\hat {M}}_{h}^{2h}(\theta )}$ of the form

where ${\displaystyle {\hat {K}}_{h}^{2h}(\theta )}$ is given by Eq. (21) and ${\displaystyle {{\hat {S}}_{h}}(\theta )}$ is a 4x4 matrix and represents the smoothing operator ${\displaystyle {{S}_{h}}(\theta )}$ given by:

In order to perform the LFA using the ${\textstyle {\mbox{ILU}}}$ solver, the smoothing operator ${\displaystyle {{S}_{h}}}$, according to [6], is given by

with

where ${\displaystyle {{\lambda }_{h}}(\theta )=1}$, for ${\textstyle {\mbox{ILU}}_{\mbox{EN}}}$, ${\displaystyle \lambda _{h}^{R}(\theta )\,={{p}_{1}}{{e}^{i(2{{\theta }_{1}}-{{\theta }_{2}})}}+{{p}_{3}}+{{p}_{2}}{{e}^{i(-2{{\theta }_{1}}+{{\theta }_{2}})}}}$ and for ${\textstyle {\mbox{ILU}}_{\mbox{NE}}}$, ${\displaystyle \lambda _{h}^{R}(\theta )\,={{p}_{1}}{{e}^{i(-{{\theta }_{1}}+2{{\theta }_{2}})}}+{{p}_{3}}+{{p}_{2}}{{e}^{i({{\theta }_{1}}-2{{\theta }_{2}})}}.}$.

The asymptotic convergence factor ${\displaystyle \rho \left(M_{h}^{2h}\right)}$ can be calculated by

where ${\displaystyle \Lambda =\left\{\theta \in {{T}^{low}}:\,{{\tilde {L}}_{h}}(\theta )=0\,\,{\text{or}}\,\,{{\tilde {L}}_{2h}}(\theta )=0\right\}}$ and ${\displaystyle \rho \left({\hat {M}}_{h}^{2h}(\theta )\right)}$ is the spectral radius of the 4x4 matrix ${\displaystyle {\hat {M}}_{h}^{2h}(\theta )}$.

In this work, LFA was used to determine the asymptotic convergence factor of the multigrid method ${\displaystyle {\bigl (}\rho \left(M_{h}^{2h}\right)={{\rho }_{loc}}{\bigr )}}$ combining ${\displaystyle {\mbox{ILU}}}$ solvers in several directions (such as ${\textstyle {\mbox{ILU}}_{\mbox{EN}}}$ and ${\textstyle {\mbox{ILU}}_{\mbox{NE}}}$), ${\displaystyle {\mbox{FW}}}$, HW, INJ, ${\textstyle {\mbox{PW}}_{\mbox{x}}}$ and ${\textstyle {\mbox{PW}}_{\mbox{y}}}$ restriction operators and bilinear and 7-point prolongation operators.

4. Numerical results

An orthotropic diffusion equation was solved using 7-point ILU solver in different directions. Several restriction operators and two prolongation operators were employed. We proposed an algorithm that presents the lowest asymptotic convergence factor values and the lowest computational cost for the multigrid method.

Equation (4) was assessed for ${\displaystyle \varepsilon ={{10}^{\kappa }}}$ and ${\displaystyle \varepsilon ={{10}^{-\kappa }}}$, with ${\displaystyle \kappa \in \mathrm {K} =\{0,\,\,1,\,\,2,\,\,3,\,\,4,\,\,5,\,\,6,\,\,7\}}$. When ${\displaystyle \varepsilon ={{10}^{\kappa }}}$ or ${\displaystyle \varepsilon ={{10}^{-\kappa }}}$ in this work, there is symmetric orthotropy. For instance, ${\displaystyle \varepsilon ={{10}^{2}}}$ is an orthotropy symmetric to ${\displaystyle \varepsilon ={{10}^{-2}}}$.