Poroelasticity problem: numerical difficulties and efficient multigrid solution

This work contains some of the more relevant results obtained by the author regarding the numerical solution of the Biot’s consolidation problem. The emphasis here is on the stable discretization and the highly efficient solution of the resulting algebraic system of equations, which is of saddle point type. On the one hand, a stabilized linear finite element scheme providing oscillation-free solutions for this model is proposed and theoretically analyzed. On the other hand, a monolithic multigrid method is considered for the solution of the resulting system of equations after discretization by using the stabilized scheme. Since this system is of saddle point type, special smoothers of “Vanka”-type have to be considered. This multigrid method is designed with the help of an special local Fourier analysis that takes into account the specific characteristics of the considered block-relaxations. Results from this analysis are presented and compared with those experimentally computed.


Introduction to poroelasticity problem
Poroelasticity theory describes the interaction between the deformation of an elastic porous material and the fluid flow inside of it. This coupling was already taken into account in the early one-dimensional work of Terzaghi, considered the father of soil mechanics, which was based on laboratory experimentation, see [38]. However, Maurice Biot, known as the founder of the theory of poroelasticity, was who established the three-dimensional mathematical formulation in the forties, see for example [2]. Nowadays, the analysis and numerical simulation B Carmen Rodrigo carmenr@unizar.es 1 Departamento de Matemática Aplicada, Universidad de Zaragoza, Zaragoza, Spain of Biot's models has become a trend topic due to their wide range of applications, and therefore, the study of this type of problems is of great interest to scientists and engineers. Biot's models are used today in a great variety of fields, ranging from geomechanics and petrol engineering, where these models have been applied ever since its inception, to biomechanics or even food processing more recently. Some examples of applications in geosciences include petroleum production, solid waste disposal, carbon sequestration, soil consolidation, glaciers dynamics, subsidence, liquefaction and hydraulic fracturing, for instance. In biomechanics, the poroelastic theory can be used to describe tumor-induced stresses in the brain (see [34]), which can cause deformation of the surrounding tissue, and bone deformation under a mechanical load (see [37]), for example.
We consider the quasi-static Biot model for soil consolidation. We assume the porous medium to be linearly elastic, homogeneous, isotropic and saturated by an incompressible Newtonian fluid. According to Biot's theory [2], the consolidation process must satisfy the following system of equations: Constitutive equation: σ = λtr( )I + 2μ , in , Compatibility condition: Darcy's law: Continuity equation: where λ and μ are the Lamé coefficients, α is the Biot-Willis constant which we will assume equal to one, κ is the permeability of the porous medium, η the viscosity of the fluid, I is the identity tensor, u is the displacement vector, p is the pore pressure, σ and are the effective stress (the part of the total stress that is not carried by the fluid, which is the stress applied to the grains of the porous medium) and strain tensors for the porous medium and q is the percolation velocity of the fluid relative to the soil. The right hand term g is the density of applied body forces and the source term f represents a forced fluid extraction or injection process. Here, we consider a bounded open subset ⊂ R n , n ≤ 3 with regular boundary . This mathematical model can also be written in terms of the displacements of the solid matrix and the pressure of the fluid. We assume g(x, t) = 0 for simplicity in the presentation, and the system can be written as follows, To complete the formulation of the problem we must add appropriate boundary and initial conditions. For instance, p = 0, σ n = t, on 1 , where n is the unit outward normal to the boundary and 1 ∪ 2 = , with 1 and 2 disjoint subsets of with non null measure. For the initial time, t = 0, the following incompressibility condition is fulfilled Let L 2 ( ) be the Hilbert space of square integrable scalar valued functions defined on , and let H 1 ( ) denote the subspace of L 2 ( ) of functions with first (distributional) derivatives in L 2 ( ). To present the variational formulation of the problem, we introduce the following functional spaces Q = q ∈ H 1 ( ), q = 0 on 1 and U = u ∈ (H 1 ( )) n , u = 0 on 2 . Then, considering the bilinear forms a(u, v) = 2μ n i, j=1 the variational formulation for problem (6)- (7) with the boundary and initial conditions (8)- (9) consists of the following: For each t ∈ (0, T ], find (u(t), p(t)) ∈ U × Q such that with the initial condition (∇ · u(0), q) = 0, ∀ q ∈ L 2 ( ), (12) and where

Stable discretization of the model
The numerical solution of poroelastic problems is usually approached using finite element methods, see for example the monograph of Lewis and Schrefler in [20] and the papers in [18,19,21,23]. It is well-known that the standard finite element solution of the poroelasticity equations can present strong nonphysical oscillations in the fluid pressure, see for instance [1,7,8,15,29]. For example, this is the case when linear finite elements are used to approximate both displacement and pressure unknowns. As in Stokes problems, approximation spaces for the vector and scalar fields, satisfying LBB condition (see [6]) can be used. This approach has been theoretically investigated by Murad et al. [25][26][27]. However, the LBB property is not a sufficient condition to ensure a free-oscillatory behavior of the solution. These nonphysical oscillations may be removed by using a very fine grid in space, which is not practical. Therefore, in order to avoid the nonphysical oscillations of the discrete solution, for example, one can add certain stabilization terms to the Galerkin formulation. We have applied this strategy in [1] to provide a stable scheme by using linear finite element approximations for both unknowns. More concretely, a time-dependent artificial term was added to the flow equation. The stabilization parameter, which depends on the elastic properties of the solid and on the size of the triangulation, was given a priori, and we showed its optimality in the one-dimensional case. This scheme provides solutions without oscillations independently of the chosen discretization parameters. In [1], any theoretical result on this stabilized scheme was given, but recently, we have presented a convergence proof in [30]. Also in this same paper, by using the same technique, we also present a stabilization of the Stokes-stable MINI element scheme for the poroelasticity problem.

One-dimensional problem
To illustrate the cause of the pressure oscillatory behaviour, we consider the following onedimensional test problem which corresponds to a column of height H of a porous medium saturated by an incompressible fluid, bounded by impermeable and rigid lateral walls and bottom, and supporting a load σ 0 on the top which is free to drain. This problem can be written as with boundary and initial conditions where E is the Young's modulus and K = κ/η is the hydraulic conductivity. It can be easily seen that problem (13) is decoupled, giving rise to the following heat-type equation for the pressure ∂ ∂t This means that therefore the pressure solution is monotone, and we should use a monotone scheme in order to keep this property. In order to discretize problem (13), we consider a non-uniform partition of spatial domain In this way, the domain is given by the disjoint union of elements We assume that the physical parameters E(x) and K (x) are constant on each element T i , and they are denoted here by E i and K i respectively. We consider an uniform time discretization by using a backward Euler method. Regarding the space discretization, we consider linear finite elements for both displacement and pressure. In this case, the following linear system of equations has to be solved on each time step where m ≥ 1, and τ is the time discretization parameter. It is clear that the pressure at time level m must satisfy the following equation where C l = −G T l A −1 l G l is a tridiagonal matrix such that for an interior node x i it is given by Notice that the scheme associated with the above equation should be an appropriate discretization for problem (14). Depending on the relation between the space and time discretization parameters, the off-diagonal elements of matrix C l + τ A p could be positive and therefore the cause of possible non-physical oscillations in the approximation of the pressure. To avoid these instabilities, the following restriction holds, For example, in the case of an uniform-grid of size h and constant values of the parameters E and K in the whole domain, such restriction becomes h 2 < 4E K τ . To confirm these unstable behaviour, we solve system (13) in the computational domain (0, 1) by using linear finite elements considering K E τ = 10 −6 . In this case, it is necessary a mesh of at least 500 nodes to fulfill the restriction. In Fig. 1a, b we show the corresponding approximation of the pressure at the first time step, for two different values of h, that is, h = 1/32 and h = 1/500. Besides, we have plotted the analytical solution of the problem (see [1]). We can observe that strong non-physical oscillations appear for this type of finite element approximations, when the space discretization parameter is not small enough. It is clear that this is due to a lack of monotonicity of the scheme. Other authors have related these oscillations to the locking effect and/or the fact that the pair of finite element does not satisfy an inf-sup condition. However, since our test is an one-dimensional problem, elastic locking can not appear, and therefore, in general, this can not be the only cause of this oscillatory behaviour. Regarding the LBB condition, we are going to show now that this property is not enough to avoid the presence of the spurious oscillations. We consider the Taylor-Hood finite element method, which approximates the displacement by continuous piecewise quadratic functions and the pressure by continuous piecewise linear functions. It is well-known that this pair of finite elements satisfies the LBB condition. Following similar computations as for the P1-P1 case, and considering again the backward Euler scheme, we obtain the following lineal system of equations on each time step ⎡ where A l , G l correspond again to the linear basis functions whereas A b , G b are associated with the bubble basis functions. In this case, the pressure at time level m satisfies the equation ), (20) where C l is as in (17) and We can observe that the off-diagonal entries of matrix C b have the right sign, but again depending on the values of the involved parameters, the whole matrix C l + C b + τ A p can still have positive off-diagonal terms. To avoid this, on each element the restriction max 0≤i≤n−1 must be fulfilled. Summarizing, we remark that the use of quadratic finite elements for displacement contributes in a positive way to the reduction of the oscillations, but it is still not enough. To illustrate this behaviour, we consider again system (19) on an uniform grid of size h and constant coefficients E and K . In this particular case, the restriction (21) is simplified to h 2 < 6E K τ , and when E K τ = 10 −6 it is deduced that 409 nodes are needed to ensure a non-oscillatory behaviour. In Fig. 2 we show the corresponding approximation of the pressure at the first time step, for two different values of h, that is, h = 1/32 and h = 1/409. Notice again that in the first case non-physical oscillations appear, and then we can conclude that the LBB property is not enough to ensure the monotonicity of the scheme.
To avoid the restrictions (18) for P1-P1 and (21) for P2-P1, which can yield to the necessity of using a very fine grid, we have proposed stable schemes providing oscillationfree solutions independently of the chosen parameters. Since schemes (16) and (20) should be suitable discretizations of the heat-type Eq. (14), the idea is to add artificial terms in order to recover the standard monotone linear finite element discretization of such equation. With this purpose, we define the following tridiagonal matrix where ε = 1/4 for the linear finite element pair and ε = 1/6 for the Taylor-Hood method. Then, it is clear that the perturbation of scheme (16) or the perturbation of (20) results in the monotone standard discretization by linear finite element method with masslumping. Notice that this perturbation corresponds to add to the variational formulation of the second equation of system (13) the term Finally, in Fig. 3 we show the approximation for the pressure obtained by using the proposed stabilization scheme for both the linear finite element pair and the Taylor Hood method, and considering h = 1/32. Notice that for this value of the spatial discretization parameter, not satisfying the restriction, we obtained oscillatory solutions. However, with the stabilized schemes (23) and (24) we obtain oscillation-free solutions. Next, we illustrate the appearance of non-physical oscillations in the pressure field when a low-permeability is assumed in a region of the domain. We consider a porous material on which a low-permeable layer (K = 10 −8 ) is placed between two layers with unit permeability (K = 1), as shown in Fig. 4, see [15]. The boundary of the squared domain is split in two disjoint subsets 1 and 2 on which we assume the following boundary conditions: on the top, which is free to drain, a uniform load is applied, that is, whereas at the sides and bottom that are rigid the boundary is considered to be impermeable, that is, Zero initial conditions are considered for both variables, and the time step is chosen as τ = 1. This test can be reduced to a one-dimensional problem. Therefore, in the following simulations we will show the numerical solutions corresponding to one vertical line in the domain as displayed in Fig. 4. First we approximate the proposed model problem by using linear finite elements for displacements and pressure. If no stabilization term is added to the discrete formulation, the approximation for the pressure field that is obtained by using 32 elements in the grid is shown in Fig. 5a. We can observe that strong spurious oscillations appear in the part corresponding to the low-permeable layer. However, if the proposed stabilized scheme is used for the simulation with the same number of nodes, the oscillations are completely eliminated and the method gives rise to the real physical solution for the pressure, as we can see in Fig. 5b.
Next, the MINI element method is considered to approximate our model problem. We consider the same number of elements in the mesh. Similarly to the previous case, when no stabilization parameter is included in the formulation, the oscillatory behaviour of the pressure approximation appears, as shown in Fig. 6a. Notice that the oscillations are much weaker than in the case of P1-P1 elements, but still are not eliminated by using this pair of finite elements. Again, if the proposed stabilization is considered for the solution of the problem, an oscillation-free approximation for the pressure field is obtained as displayed in

Multi-dimensional problem
We now consider the finite element approximation of problem (10)-(12) by using linear finite elements for both unknowns or the well-known MINI element. Let T h be a regular triangularization of satisfying the usual admissibility assumption. In the two-dimensional case, we can define finite element approximations for U and Q as U h and Q h , respectively. We can write the semidiscrete Galerkin approximation to the problem (10)- (12) as such that In a standard way in the finite element framework, we can consider the usual basis functions {ψ i } n u i=1 for the displacements and those for the pressure {φ i } n p i=1 corresponding either to P1-P1 or MINI element schemes. In this way, we find that the semi-discrete formulation (28)- (29) can be expressed as a system of linear ordinary differential equations Here U represents the unknown vector (u 1 (t), u 2 (t), . . . , u n u (t)) and P the unknown vector ( p 1 (t), p 2 (t), . . . , p n p (t)), and U t and P t the corresponding vectors with the time derivatives of the unknowns. The matrix A is the elasticity matrix, −A p is the diffusive matrix, B is the gradient matrix and B is its transpose, that is minus the divergence matrix. H and F are the right hand side vectors with components To discretize in time, we consider the backward Euler scheme: which yields the following linear equations system on each time step where τ is the time discretization parameter and BU 0 = 0.
In the multi-dimensional case, we cannot proceed as in 1D since the Schur complement becomes a full matrix and therefore it is very difficult to obtain any result about the monotonicity of the discrete schemes. However, we can extend the idea proposed for the onedimensional case to obtain stabilized discretizations for (6)- (8). The final scheme results to be the corresponding discretization of the problem where a perturbation term − ∂ ∂t ∇ ·(β∇ p) has been added to the flow equation, with parameter The corresponding semidiscrete Galerkin approximation is as follows, where β = β η κ . Considering the backward Euler scheme to discretize in time, we get the discrete variational problem: which in matrix-form reads as In order to assess the performance of the proposed method, we solve numerically a poroelastic problem on a cylindrical shell of deformable porous material with a uniform load on the inner boundary. The outer boundary is constrained by a rigid body that offers no resistance to the passage of the fluid, so the boundary conditions for the displacement and for the pressure are u = 0 and p = 0 respectively. In the inner surface of the cylinder, there is a fixed pressure, p = 1, which produces a uniform load σ n = (cos θ, sin θ). A scheme of the geometry and boundary conditions are given in Fig. 7.
The material properties of the porous medium are given in Table 1 where λ and μ are related to the Young's modulus E and the Poisson's ratio ν by . Figure 8a shows how standard finite elements P1-P1 lead to spurious oscillations in the pressure approximation at time T = 10 −6 using one time step. These nonphysical oscillations in the pressure are eliminated completely adding the stabilization term h 2 a p ( ∂ ∂t p h , q h )/4(λ+ 2μ) to the flow equation as can be seen in Fig. 8b.

Stability and convergence of the discrete schemes
In this section we study the stability of the proposed stabilized schemes. By stability, here we mean bounds on the inverse of the discrete operator (for a fixed time step). This is analyzed by proving an inf-sup condition for the different discretizations that we consider, namely stabilized P1-P1 and MINI element schemes. From now on, we consider that operator , and C is bounded, selfadjoint and positive definite. Since C may take different form for different discretizations, we do not specify its definition here. Notice that A is bounded, selfadjoint and positive definite, and B is also bounded. Moreover, we define the following norm

Theorem 1 We consider the discretization of poroelasticity problem given by an operator of the form A C , as previously explained. Then, A C is an isomorphism if and only if for any
Proof First we assume that (40) holds and we show that A C is an isomorphism. For this purpose, we introduce the bilinear form From the continuity of operators A and B, we can easily show that operator A C is bounded in ||| · |||. Moreover, from the inf-sup condition (40), for any Without loss of generality, we may assume that w h A = p h , and then we have, with w h defined as above and some θ > 0 that will be determined later. Then, we can write, Applying the inequality ab ≥ − 1 2θ Also, by using (41), it is satisfied that, In this way, On the other hand, the triangle inequality implies that with γ 1 depending only on γ B . Hence, and we can conclude that A C is an isomorphism. Next, we assume the invertibility of A C and we show that condition (40) is fulfilled. For Then, the invertibility of A C implies that We can obtain estimates for Bv h,q Q h and Cq h Q h in order to deduce (40 Finally, the inf-sup condition (40) easily follows by combining the last two estimates.
A more general version of the previous result can be found in [30] for a general operator where V and Q are Hilbert spaces and V and Q their dual spaces, operator A : V → V is bounded, selfadjoint and positive definite, operator B : V → Q is also bounded and C : Q → Q is a bounded operator, sefadjoint and positive (semi)definite.
Notice that if the inf-sup condition in (40) is satisfied with C = 0, then it is also satisfied for any C bounded, selfadjoint and positive (semi)definite. This implies that a stable finite element pair for Stokes equations is also stable for poroelasticity. From this, we can conclude that MINI element is stable for poroelasticity, since it is well-known that this finite element scheme satisfies an inf-sup condition for Stokes system, see [3] for example. The discrete space of functions corresponding to MINI element scheme satisfies that is a direct sum of the space of piece-wise linear continuous vector valued functions and the space of bubble functions.
Then, we can write the discrete problem (31)-(32) in the following block form: By eliminating the equation corresponding to bubble functions, we obtain the following operator, Therefore, changing accordingly the right hand sides, the systems of equations corresponding to (42) and (43) are equivalent and solvability for MINI element implies solvability for P1-P1 + S b . We then have the following theorem from which arises that S b improves the "monotonicity" of the P1-P1 scheme.
Theorem 2 S b is spectrally equivalent to h 2 L, where h is the grid-size and L is the discretization of Laplace operator by piece-wise linear continuous finite elements, that is, S b h 2 L with constants only depending on the shape regularity of the mesh.
This theorem can be proven by using the results in the Appendix of [30], which hold for one, two and three spatial dimensions. Moreover, this statement implies that if the stiffness matrix for linear continuous elements is an M-matrix, then S b is an M-matrix as well, and therefore adding the term S b improves the properties of monotonicity of the resulting discrete problem. Finally, the following result, together with the previous theorem, justifies the addition of the proposed stabilization terms to both the MINI element and the P1-P1 discretizations.

Theorem 3
Suppose that A C , as in Theorem 1, is an isomorphism, and that D is spectrally equivalent to C, namely α 0 q h C ≤ q h D ≤ α 1 q h C for some positive constants α 0 and α 1 . Then A D is an isomorphism.
From this result, we can state that any operator C, spectrally equivalent to τ A p + S b will result in a stable discretization of the Biot's model. The perturbations spectrally equivalent to S b are of the form where C T , T ∈ T h are constants independent of the mesh size h or τ , and therefore we can conclude that an inf-sup condition is satisfied for stabilized P1-P1. Once the stability of the two proposed stabilized schemes has been deduced, we can obtain estimates for the analysis of the fully discretized time dependent Biot's model. A detailed error analysis of the problem can be found in [30], and next we summarize the main results presented in that paper. In the next, we will denote by c a generic constant independent of time step, mesh size and other parameters.
We can write the fully discretized scheme at time t n , n = 1, 2, . . . as The initial data u 0 h and p 0 h are given by the following stabilized Stokes equation, Although we do not consider this case here, the initial data can also be defined satisfying divu 0 h = 0, p 0 h = 0. In this case, similar convergence results can be obtained (see [30]). We define the following norm on the finite element spaces: (u, p) τ,h := u 2 a + τ p 2 a p + εh 2 ∇ p 2 1/2 , and the following elliptic projections u h and p h for t > 0, such that In this way, we can write the discretization error as follows, and study separately estimates for each of the terms. For the error of the elliptic projections (by using P1-P1 or mini-element), we have, for any t, ρ u a ≤ ch(|u| 2 + |p| 1 ), Besides this results, we also have the corresponding estimates for ∂ t ρ u and ∂ t ρ p . We can prove the following auxiliary estimates of the error between the elliptic projection and the numerical solution: . Furthermore, we also have the following estimate in the L 2 -norm, The next step is to find estimates for the terms corresponding to w j u and w j p :

Lemma 5 Let u(t) and p(t) be the solution to the model Biot's problem, w
Finally, assuming extra regularity of u(t) and p(t), as usual for convergence analysis of the finite element method, we have the following theorem for both MINI element and P1-P1 stabilized element, which provides estimates for the final error.

Theorem 6 Let u(t) and p(t) be the solution of the Biot's system. Then we have the following estimate for displacement u(t),
For pore pressure p(t), we have the estimate and we also have the following error estimate in L 2 -norm, Previous error estimates consists of two parts: the error for t > 0, which gives optimal convergence order; and the error in the approximation of the initial data, which can be split in two terms, one depending on the errors due to the elliptic projection and the other one corresponding to the errors due to the choice of initial conditions. Due to the choice of our initial data, satisfying a stabilized Stokes equation, the initial errors strongly depend on the regularity of the initial data, and in particular the assumptions on the regularity of the initial pore pressure are crucial. By assuming p(0) ∈ H 1 0 ( ), the standard error estimates for the elliptic projection and stabilized Stokes equation show that the initial data errors are appropriately bounded, and therefore optimal order of convergence is achieved for the discrete problem. In this way, the overall convergence rate of the stabilized MINI element is optimal. However, if we assume that p(0) is merely in L 2 ( ), then we cannot expect that the errors in the initial data are of optimal order, and, therefore, the overall convergence rate of the stabilized MINI element is not optimal as well.
of the models. Biot's models lead to computationally complex problems for which traditional simulations become too expensive. Therefore, fast numerical algorithms have to be designed for their solution. Here, we deal with the design of a monolithic multigrid solver for the highly efficient solution of the resulting algebraic systems of equations. Special smoothers suitable for saddle point problems, as the Vanka type relaxations, are considered. By using this approach, an efficient multigrid solver is obtained for the solution of the quasi-static Biot's model for soil consolidation.
Since their development in the 70's, multigrid methods (see [14,36,39], for example) have been proved to be among the most efficient numerical algorithms for solving the large sparse systems of equations arising from the discretization of partial differential equations, achieving asymptotically optimal complexity at least for elliptic problems. These methods are iterative solvers mainly based on the acceleration of the convergence of classical iterative methods by using solutions obtained on coarser meshes as corrections. Two ideas are involved in the development of multigrid methods: the first one is the fact that some classical iterative methods have a strong smoothing effect on the components of the error corresponding to the high frequencies (high oscillating error components), and the second one is that a smooth error can be well represented on a coarser grid. These observations suggest the following structure of a two-grid cycle: 1. Perform ν 1 iterations of an iterative relaxation method S h on the fine grid (pre-smoothing) 2. Compute the defect of the current fine grid approximation 3. Restrict the defect to the coarse grid by using a restriction operator R h,2h 4. Solve the coarse grid defect equation 5. Interpolate the correction to the fine grid using a prolongation operator P 2h,h 6. Add the interpolated correction to the current fine grid approximation 7. Perform ν 2 iterations of an iterative relaxation method S h on the fine grid (postsmoothing) Following this algorithm, the two-grid error transformation operator is given by where I h denotes the identity and the subscript "2h" indicates that the coarse grid is obtained by doubling the mesh size in each space direction, which is called "standard coarsening". Instead of inverting L 2h , the coarse-grid equation can be solved by recursive application of this procedure, yielding a multigrid method. From the previous algorithm, it is clear that many details are open for discussion and decision, since all the components have to be properly chosen.
Regarding the coarse-grid correction part of the algorithm, standard inter-grid transfer operators, dictated by the geometry of the triangular grid, are considered here. This results in 7-point restriction operator and bilinear interpolation of neighboring coarse-grid unknowns. Finally, on the coarse grids, we apply direct coarse-grid discretization. However, the efficiency and robustness of a multigrid algorithm is usually strongly influenced by the smoother. Moreover, for the problem we are dealing with, an additional difficulty appears, since it results in a system of saddle point type aspect. An overview of multigrid methods for discretizations on rectangular grids of this type of problems is presented in [28], where coupled or Vankarelaxation and decoupled distributive relaxation methods appear as the most suitable for this kind of problems. Due to the fact that for some systems of equations it is a challenge to design an efficient distributive relaxation scheme, Vanka smoothers seems to be the best option. These smoothers consist of decomposing the mesh into small subdomains and treating them separately. Therefore, one relaxation step consists of a loop over all subdomains, solving for each one the system arising from the corresponding equations. Next, we give a more detailed description of the iterative method. We consider a linear system of equations A h u h = f h , which, in our case, arises from the discretization of a PDE problem. Vector u h is composed of unknowns corresponding to m different variables. More concretely, N i unknowns of each variable i are considered. Let B be the subset of unknowns involved in an arbitrary block, that is, are the global indexes of the n i unknowns corresponding to variable i. In order to obtain the matrix A B h of the system to solve associated with block B, we introduce the matrix V B representing the projection operator from the vector of all unknowns to the vector of the unknowns involved in the block, as the following block-diagonal matrix Here, each block V i B is a (n i × N i )−matrix, whose jth-row is the k i ( j)th-row of the identity matrix of order N i . In this way, matrix A B h can be defined as Therefore, this type of smoother results in a multiplicative Schwarz method with iteration matrix where N B is the number of blocks or small systems to be solved in a relaxation step of the iterative method. Very often in practice, instead of solving the local problems exactly, one can replace A B h with an approximation A B h , obtaining the so-called multiplicative Schwarz method with inexact local solver, with iteration matrix given by Therefore, many variants of box-type smoothers can be considered. They can differ in the choice of the subdomains which are solved simultaneously, and in the way in which the local systems to be solved are built. Notice that the different subdomains can also be visited in different orderings, for instance, they can also be treated with some patterning scheme, yielding to a multicolored version of these relaxation schemes. All this makes wider the variety of this type of relaxations. This class of smoothers was introduced by Vanka in [41] to solve the finite difference discretization on rectangular grids of the Navier-Stokes equations. Since then, much literature can be found about the application of this type of smoothers, mainly in the field of Computational Fluid Dynamics [16,17,40], but also in the context of Computational Solid Mechanics, for example in [44]. Here, we consider a Vanka smoother suitable for the stabilized P1-P1 discretization of poroelasticity problem on triangular grids, in which all unknowns appearing in the discrete divergence operator in the pressure equation are simultaneously updated. This way of building the blocks is very common in Vanka-relaxations used for Stokes and Navier-Stokes problems. This approach implies that for the two-dimensional problem twelve unknowns corresponding to displacements and one pressure unknown, see Fig. 9, are relaxed simultaneously, making necessary to solve a 13 × 13−system for each grid point. Then, we Fig. 9 Unknowns simultaneously updated in point-wise box Gauss-Seidel smoother, and example of the overlapping between two arbitrary blocks. Circles denote velocity unknowns whereas the square refers to pressure degrees of freedom iterate over all grid points in lexicographic order, and for each of them the corresponding box is solved. In this full variant, all the unknowns in the system are considered coupled. Therefore, the need of solving such systems makes these smoothers expensive. A simplified variant can be considered by coupling each displacement unknown in the system only with itself and with the corresponding pressure unknown. In this way, the solution of the resulting systems becomes much cheaper, making the use of this diagonal version very appealing in practice. However, it is observed that the diagonal version of the point-box smoother can be less robust with respect to some applications. In our case, the multigrid based on the diagonal Vanka smoother provides a very efficient solver for the poroelasticity equations.
On the other hand, it is well-known that the performance of multigrid methods strongly depends on the interplay between the smoothing and the coarse-grid correction part of the algorithm. This two principles can be combined by two basic approaches to multigrid solvers. In algebraic multigrid no information is used concerning the grid on which the governing PDE is discretized, and therefore it is more suitable when unstructured grids are considered. However, an alternative to this strategy is to consider a hierarchy of semi-structured grids, in which a geometric multigrid that takes advantage of the geometry of the grid, can be efficiently implemented. An initial totally unstructured grid is then considered to represent the geometry of the domain, and a regular refinement process is applied on each element of the input grid, resulting in a hierarchy of globally unstructured grids composed of structured patches, where the geometric multigrid can be implemented based on stencil-based operations. We have proposed this strategy for example in [33] for other problems like elasticity equations. This approach is also advantageous in the sense that the implementation uses stencil-based operations, since the local nature of the multigrid operators and the semi-structured character of the grid allow it, as well as it allows the use of low-cost memory storage of the discrete operator based on stencil formulation. For poroelasticity equations, we presented the efficient implementation of a multigrid finite element method on semi-structured grids in [11]. In this paper, the implementation on semi-structured grids of the linear finite element scheme is explained, focusing on the stencil-based implementation on each regular patch of the semistructured grid. Moreover, an efficient procedure to construct these stencils by means of a reference hexagon is also described. Furthermore, due to the block-structure of the considered grids, each triangle of the coarsest triangulation can be treated as a different block, and then different components can be chosen on these patches. In particular, in [9], we have proposed the use of different smoothers for triangles which have different geometries, but this choice can be done in a very efficient way with the help of the local Fourier analysis.

Local Fourier analysis for Vanka smoother based multigrid
Local Fourier analysis (LFA, or local mode analysis) is a commonly used approach for analyzing the convergence properties of geometric multigrid methods. In this analysis an infinite regular grid is considered and boundary conditions are not taken into account. LFA was introduced by Brandt in [4] and afterward extended in [5]. A good introduction can be found in the paper by Stüben and Trottenberg [36] and in the books by Wesseling [42], Trottenberg et al. [39], and Wienands and Joppich [43]. LFA was generalized to triangular grids in [10], for discretizations based on linear finite element methods. Afterwards, we extended this generalization to systems of partial differential equations [12,13].
Local Fourier analysis assumes that all the operations involved in the multigrid algorithm are local processes, that is, the operations performed on each unknown only depend on the information on nearby neighbors, what allows to neglect the effect of boundary conditions. By imposing some assumptions on the discrete operator: linear differential operators with constant coefficients, a basis of complex exponential eigenfunctions of the operator, called Fourier components, can be obtained. Summarizing, this analysis mainly simplifies the computation of the spectral radius of the k-grid iteration matrix, by considering the matrix representation of this operator with respect to such basis of eigenfunctions, which results in a block-diagonal matrix.
To perform the LFA on the considered regular triangular grids, we establish a nonorthogonal unitary basis of R 2 , {e 1 , e 2 }, in the directions of two of the edges determining the triangular grid. We also consider, for the frequency space, its reciprocal basis {e 1 , e 2 }, i.e., (e i , e j ) = δ i j , 1 ≤ i, j ≤ 2, where (·, ·) is the usual inner product in R 2 and δ i j is the Kronecker's delta, see [12]. Considering that the coordinates of a point in the basis, {e 1 , e 2 } are x = (x 1 , x 2 ), we can extend our computational grid to a regular infinite grid, where h = (h 1 , h 2 ) is a grid spacing in the coordinate system {e 1 , e 2 }. Considering the scalar Fourier modes, ϕ h (θ , x) = e iθ 1 x 1 e iθ 2 x 2 , their vector counterparts are ϕ h (θ , x) := ϕ h (θ , x) · 1, where 1 = (1, . . . , 1) t ∈ R q , (with q equal to the number of variables) x ∈ G h , and θ ∈ h = (−π/ h 1 , π/h 1 ] × (−π/ h 2 , π/h 2 ]. They give rise to the Fourier space, From the expression of the back Fourier transformation, it follows that each vector discrete function u h (x) can be written as a formal linear combination of the Fourier modes, which are linearly independent discrete functions. By considering vector operators L h defined on the infinite grid and satisfying the assumptions of the LFA, it is fulfilled that the Fourier modes ϕ h (θ , x) are formal eigenfunctions of L h . More precisely, the following relation reads where L h (θ ) is the so-called symbol of operator L h , which denotes the representation of the discrete operator on the Fourier space. Most classical iterative methods can be expressed by means of a splitting of the operator L h of the form L h = L + h + L − h , where L − h relates the part of the operator corresponding to the unknowns which have been relaxed before the current approximation, and L + h to those that are going to be updated in the current or in the following steps. In this way, the iteration matrix of the smoothing operator is given by S h = −(L + h ) −1 L − h , and it is easy to prove that the Fourier modes are also its eigenvectors, satisfying a relation like in (51), with its Fourier symbol denoted as S h (θ ). However, the overlapping block smoothers considered here do not come from such a decomposition of the discrete operator. The distinction with respect to other smoothers is that they update some variable more than once, due to the overlapping of the local subdomains which are simultaneously solved, and this fact has to be taken into account in the analysis because it causes that some intermediate errors appear apart from the initial and final errors. For this reason, Vanka smoothers require a special strategy to perform the local Fourier analysis. To our knowledge, there are only few papers dealing with local Fourier analysis for overlapping smoothers, all of them for discretizations on rectangular grids. This analysis was performed in [35] for the staggered finite-difference discretizations of the Stokes equations, and in [24] for a mixed finite element discretization of the Laplace equation. In [32] an LFA for overlapping block smoothers on triangular grids is presented. This tool was applied to linear finite element discretizations for poroelasticity problems. Later, in [22], the analysis for such overlapping block smoothers is performed on rectangular grids for finite element discretizations of the grad-div, curl-curl and Stokes equations. In [31], we present and extend this analysis to general discretizations on triangular grids, including some special techniques for the case of edge-based discretizations. Two model problems are chosen to show this analysis, but we keep in mind that it can be carried over to a variety of other problems and other overlapping smoothers. The considered problems are the discretization by stabilized linear finite elements of the Stokes problem, and the low-order Nédélec's edge elements for the curl-curl equation. Since the analysis of the Vanka smoother considered here for the poroelasticity problem in analogous to the one presented in [32] and similar to that explained in [31] for Stokes equations, we refer to the reader to these works in order to avoid to include here very technical and involved calculations. Instead, we present some results obtained with this analysis for both the full and the diagonal versions of the Vanka smoother. From now on, in all the results, for simplicity of notation we will denote k = τ κ η .
We start presenting some results for equilateral triangular grids. In Table 2, the two-grid convergence factors ρ 2g are displayed together with the experimentally measured W-cycle convergence factors ρ h , using nine refinement levels, which have been obtained with a random initial guess and a right-hand side zero to avoid round-off errors. For a value of k = 10 −8 , the results obtained for different number of pre-and post-smoothing steps, ν 1 , ν 2 , are shown for both full and diagonal point-wise box-smoothers. In this table, we can observe the good correspondence between the experimentally measured factors and the predicted ones, what indicates that the proposed analysis gives very accurate predictions of the convergence factors of the box-relaxation based multigrid method. Also it is seen that the behavior of full and diagonal box smoothers, in this case, is identical, giving rise to almost the same values.
Next, we analyze the influence of parameter k on the behavior of box-type smoothers. To this end, two pre-and one post-smoothing steps are considered, and the two-grid convergence factors ρ 2g , predicted by LFA, together with those factors experimentally computed, ρ h , are displayed in Table 3, for different values of k, varying from 10 −4 to 10 −13 , by using both full and diagonal point-wise box-smoothers. The experimentally obtained and the predicted results match perfectly for any value of parameter k and, although the obtained results are not independent of k, they are very satisfactory for any value of such parameter.
Due to the fact that the obtained results for full and diagonal point-wise box-smoothers are very similar, the latest is preferred for the computations since this approach is significantly cheaper than the full variant, and results in a very good performance.  Fig. 10 Spectral radius ρ(M 2h h ) for the diagonal point-wise box-smoother with two pre-and one postsmoothing steps for different triangles in function of two of their angles, for poroelasticity problem with parameter k = 10 −10 Despite the good behavior of the diagonal-Vanka based multigrid for equilateral triangulations, for meshes characterized by a small angle, the obtained convergence factors deteriorate, as we can see in Fig. 10. It shows the spectral radius ρ 2g = ρ(M 2h h ) for diagonal point-wise box-smoother with two pre-and one post-smoothing steps, for a value of k equal to 10 −10 , for a wide range of values of angles α and β, which are geometric parameters determining the shape of the triangular grid.   Fig. 11 Coarsest triangulation of the computational domain and obtained grid in the hierarchy after four refinement levels In order to improve the results obtained for a very anisotropic grid, it is possible to consider the use of a line-wise Vanka smoother. This relaxation simultaneously updates all the unknowns involved in the blocks corresponding to the pressure points in a whole line. It is also possible to perform a local Fourier analysis for this line smoother, and some results are presented here in order to show how these smoothers can improve the obtained convergence factors. To this end, we consider an isosceles triangular domain with smallest angle of 10 • . In Table 4, the predicted two-grid convergence factors, ρ 2g by considering two pre-and one post-smoothing steps, are shown for both line-wise and diagonal point-wise box-smoothers and for different values of parameter k. It is observed how the behavior of the diagonal point-wise box-smoother deteriorates whereas line-wise box-relaxation gives a significant improvement.
Finally, in order to demonstrate the good behavior of the proposed box-relaxation based multigrid, we present some results for the previously considered poroelastic problem on a cylindrical shell. In order to apply the proposed box-relaxation based geometric multigrid, a hierarchy of grids must be defined as in previous experiments. Then, for the unstructured coarsest grid, we have chosen the triangulation shown in left Fig. 11, which is composed of 18 quite regular triangles. From this triangulation, a regular refinement process is applied giving rise to a hierarchy of meshes. In particular, due to the curvature of the boundary, the refinement is performed leading to a good fitting of the grid to the real domain boundary as more refinement is made. For instance, the obtained grid after four refinement levels is  Fig. 11, as an example, where it can be seen how well the refined grid approximates the real boundary. Then, this problem is solved by applying the proposed multigrid method, considering the diagonal point-wise box-smoother, and by using an F-cycle with two pre-and one postsmoothing steps. In Fig. 12, the obtained convergence of the multigrid algorithm for this poroelastic problem for different numbers of refinement levels is displayed. Satisfactory results are obtained, since for such complicated problem the residual is reduced to a value of 10 −6 after about 16 iterations of the multigrid algorithm.