A Parallelizable Space-Time Approach in the Uzawa Solver with Multigrid for Poroelasticity

Poroelasticity refers to the dynamic interrelation between fluid flow and solid deformation within porous media, characterized by solid structures containing fluid-filled pores. External loads exerted on the porous medium influence the entire system. Pressure variations within the fluid-filled pores subsequently influence fluid movement, while the solid material undergoes elastic deformation in response [1].

Karl von Terzaghi’s (1883-1963) work on the one-dimensional theory of soil consolidation, rooted in empirical observations, and Maurice A. Biot’s (1905-1985) approach, which incorporates the compressibility of the solid component within porous media, jointly established the framework for understanding the coupling between porous elastic solid deformation and the linear viscous fluid flow. Since then, it has been established that the principles of continuous poroelasticity and Darcy’s law are derived from the microscopic Navier equations, which dictate the deformation of the elastic solid constituting the porous medium’s solid component, and the Navier-Stokes equations, which govern flow at the pore scale [2].

The study of behavior in poroelasticity is fundamental for developing accurate models that describe porous materials across various contexts. Poroelasticity facilitates modeling applications such as subsurface exploration for oil [3, 4], groundwater flow using brine in porous media [5, 6], assessing displacements of geotechnical structures [7, 8], examining mechanical properties of bones and biological tissues [9-11], analyzing noise generated by pile driving in the wind power industry through seabed modeling [12], agricultural engineering [13] and numerous other applications.

The mathematical model of the poroelasticity problem encompasses two key relationships: the displacement and pressure of the fluid within the porous medium, and the structural displacement of the porous matrix, as elucidated by Biot’s theory. This results in a system of equations comprising the momentum balance equation for the porous system and the mass conservation equation for the fluid within the pores.

One approach to solving the resultant system of differential equations is through numerical methods, necessitating the discretization of the differential equations. This discretization can be achieved using various techniques, including the finite difference method [14], and finite volume method [15], among others.

Saddle point problems frequently arise in the numerical analysis of solutions to partial differential equations (PDEs) with constraints, such as diffusion-advection problems, Stokes problems, and Darcy problems. Effectively addressing saddle point problems necessitates the utilization of smoothers capable of efficiently managing the coupling between displacement and pressure [16].

Several studies have been conducted in the domain of poroelasticity problem resolution. Oosterlee and Gaspar [17] introduced the Vanka smoother, which addresses variables in a mutually dependent manner concerning both the Stokes equations and the incompressible poroelasticity problem, focusing on the 2D spatial context. Their findings indicated superior performance, particularly on collocated meshes.

Borregales et al. [18] proposed a methodology aimed at resolving the 2D and 3D poroelasticity problem, leveraging an updated version of the Fixed-Stress method. This approach integrates the temporal variable as an additional dimension to effectively parallelize the Biot problem. Through a series of numerical experiments and convergence analyses, they provided empirical evidence of the method’s robustness.

The Parallel Full Approximation Scheme in Space and Time (PFASST) methodology operates by interleaving Spectral Deferred Corrections (SDC) iterations with varied resolutions, similar to the Parareal algorithm. It incorporates a Full Approximation Scheme (FAS) correction to enhance the accuracy of coarser SDC iterations. Parallelization along the time direction facilitates the computation of solutions at subsequent time steps based on solutions obtained at preceding time points, employing spectral corrections [19].

Franco and Pinto [20] employed colored ordering in conjunction with the Vanka smoother within the Finite Difference Method framework to address the poroelasticity problem. Their approach also integrated the Multigrid method and Waveform Relaxation temporal sweep technique.

Adler et al. [21] introduced a finite element discretization approach to tackle the poroelasticity problem. They developed block preconditioners coupled with Multigrid and Fixed-Stress techniques, demonstrating a robust methodology capable of accommodating fluctuations in physical parameters, mesh sizes, and time steps inherent to the poroelasticity problem.

The Uzawa solver stands as an efficient iterative technique for addressing saddle point systems. Due to its simplicity and low memory requirements, it has garnered extensive utilization [22]. Subsequently, recent investigations into the classical iteration of this method or its integration with other methodologies, aimed at enhancing convergence, are outlined.

Miao and Zhang [22] introduced the Uzawa Single-Step Iteration (Uzawa-SSI) method as a solution for both non-Hermitian nonsingular and singular saddle point problems. Their study involved a comparative analysis with other Uzawa variants, highlighting the simplicity of convergence conditions in Uzawa-HSS (Uzawa Hermitian and Skew-Hermitian Splitting). They provided numerical evidence showcasing the efficiency of the Uzawa-SSI method. However, it is worth noting that this method entails the evaluation of a parameter $\lambda$, which can pose challenges in practice.

Axelsson and Karátson [23] proposed a Krylov-Uzawa method tailored for the Stokes problem, aimed at accelerating convergence. Their analysis established that the Krylov-Uzawa iteration may demonstrate superlinear convergence.

Keram and Huang [24] devised a Uzawa-type algorithm leveraging the mixed finite element method for resolving the steady-state natural convection model. Comparing it with the conventional Uzawa algorithm, they noted a reduction in CPU time with the proposed approach.

Wu and Gao [25] introduced a diagonally preconditioned Uzawa Splitting method tailored for saddle point problems, offering a sufficient condition to guarantee convergence. Their findings illustrated that the proposed method outperforms certain existing methodologies in the literature. Furthermore, they observed that the preconditioner can enhance the efficiency of associated Krylov subspace methods.

Uzawa methods coupled with the Multigrid technique offer the benefit of accelerated convergence. This notion is supported by several references:

Luo et al. [26] implemented the Uzawa smoother in solving the 2D poroelasticity equations utilizing finite volume discretization, in conjunction with the Multigrid method. Their proposed smoother combines Symmetric Gauss-Seidel iterations for displacements with Richardson iterations for the Schur complement within the pressure field. The relaxation parameters were determined through Local Fourier Analysis. Numerical simulations corroborated the efficacy and reliability of the proposed approach.

Badea [27] showcased the convergence of Uzawa and Arrow-Hurwicz algorithms integrated with Multigrid for resolving general saddle point problems, characterized by constrained optimizations.

In a related study, Spies et al. [28] conducted a comparative analysis of three relaxation scheme options – Braess-Sarazin, Vanka, and Schur-Uzawa – utilizing a block triangular preconditioner and Multigrid method to address saddle point problems associated with the Stokes equations and the Oseen equation.

In contrast to previous works, this article proposes a solution to the poroelasticity problem through a novel parallelizable spatio-temporal scanning approach, combining the Uzawa method with Multigrid techniques.

The structure of the paper unfolds as follows: Section 2 delineates the mathematical and numerical models underpinning the problem; Section 3 delineates the Multigrid method; Section 4 introduces the novel time sweep method; Section 5 elucidates the numerical outcomes of the simulations; and finally, Section 6 encapsulates the conclusions drawn from the study.

The Eqs. (9)-(11) give rise to large-scale sparse linear systems, for which direct methods are not recommended [14]. Classical iterative techniques for solving such linear systems, such as Gauss-Seidel methods, excel at eliminating only the oscillatory components of errors associated with dependent variables (high-frequency Fourier modes). However, they typically struggle to efficiently eliminate smooth error components (low-frequency Fourier modes). As a result, these methods are often referred to as smoothers, given their tendency to smooth out errors. Consequently, at the outset of the iterative process, a high convergence rate is observed, stemming from the rapid elimination of oscillatory components [35].

The Multigrid method’s effectiveness lies in its ability to solve systems of equations, both linear and nonlinear, by employing multiple grid levels. This approach accelerates the convergence of the iterative method used for smoothing by going through the entire spectrum of Fourier error modes, thereby transforming smooth error modes into more oscillatory ones.

The Multigrid cycle combines iterations on fine and coarse grids, starting at the finest grid level (referred to as h), smoothing the system of equations, and then restricting to the immediately coarser grid (denoted as 2h). This process iterates until reaching the coarsest grid level, followed by prolongation and return to the immediately finer grid. The method involves alternating restriction and interpolation operations between fine and coarse grids.

As outlined by Wesseling [36], the sequence in which the various grids are traversed determines the Multigrid cycle. The most prevalent cycles include the V, W, and F types. The Wcycle is adopted in this article for its robustness, as it is the preferred choice for problems involving poroelasticity [17, 37]. Figure 3 illustrates a W -cycle for four levels of grids, i.e., $n g l=4$, where $n g l$ is the number of grid levels.

3.png

Figure 3. Representation of the W-cycle for 4 levels of grid

In the context of transferring information among different grids, the choice between residue and/or solution transfer depends on the nature of the problem. For linear problems, the correction scheme (CS) is deemed more suitable as it exclusively transfers the residue. Conversely, for nonlinear problems, the full approximation scheme (FAS) is preferred as it transfers both residue and solution [35, 36]. Given that the problem under consideration in this work is fundamentally linear, only the CS scheme is employed. Algorithm 1 outlines a W-cycle with CS.

In this algorithm, $A$ represents the coefficient matrix, while $f$ denotes the vector of independent terms obtained by discretizing the Eqs. (5)-(7).

Here, $\Lambda$ represents a generic variable, and $\bar{\Lambda}$ is its approximation. $I_h^H$ and $I_H^h$ represent, respectively, the space restriction and prolongation operators (utilized at each time step), while $v_1$ and $v_2$ indicate the numbers of pre- and postsmoothing steps, respectively.

Within the Multigrid method, various smoothers can be employed, including classical solvers such as Gauss-Seidel with different formulations like Zebra, Red Black, and incomplete LU decomposition [32]. Zebra-type smoothers involve partitioning the domain into planes, sequentially solving each odd plane before addressing each even plane independently [19].

The Uzawa method is a smoother that involves decomposing the coefficient matrix obtained in the discretization (Eqs. (9-11)) into two matrices, as described in the system of equations Eq. (12) [38, 26]:

$\left[\begin{array}{cc}A & B^t \\ B & -C\end{array}\right]=\left[\begin{array}{cc}M_A & 0 \\ B & -\omega^{-1} I\end{array}\right]-\left[\begin{array}{cc}M_A-A & -B^t \\ 0 & C-\omega^{-1} I\end{array}\right]$   (12)

where, $\omega$ is a positive term called the relaxation parameter and $M_A$ represents a preconditioner for the positive definite matrix $A$ that makes the process less expensive for each iteration.

The matrix $B$ represents the negative discrete divergence matrix, $B^t$ stands for the discrete gradient matrix, and $I$ denotes the identity matrix. $C$ is a diagonal matrix primarily comprising negligible values due to its composition of hydraulic permeability, which tends to be very small. Additionally, it integrates the time discretization parameter $\tau$ as a multiplication factor, which can be arbitrarily small.

As a solver, Uzawa utilizes the Symmetric Gauss-Seidel method, which can be described as an incomplete LU decomposition [26].

The algorithm for the standard Uzawa method involves updating the variables $\boldsymbol{u}=(u, v)$ (displacement) and $p$ (pressure) vectors based on their values at the previous iteration $\boldsymbol{u}^0$ and $p^0$ as described below:

$\boldsymbol{u}=A^{-1}\left(U-B^t \boldsymbol{u}^0\right)$   (13)

$p=p^0+\omega(B \boldsymbol{u}-\mathcal{P})$    （14）

Therefore, the decoupled Uzawa method is equivalent to one iteration of the Richardson method, which can be employed to determine an optimal $\omega$ in terms of the maximum and minimum eigenvalues of the Schur complement matrix $S=B A^{-1} B^t$ [39, 40].

Luo et al. [26] conducted a systematic analysis of the relaxation parameter $\omega$ using Local Fourier Analysis, convergence, and acceleration techniques to enhance efficiency and accuracy in solving complex poroelastic problems.

The algorithm for the solution using the Uzawa smoother initially smooths the $\boldsymbol{u}$ vector for all control volumes, using $p^0$ and $\boldsymbol{u}^0$. The systems of this iterative process make use of the Symmetric Gauss-Seidel method as a solver. Subsequently, with this data, the pressure residual is determined and then $p$ is updated.

Above is Algorithm 2 for solving the poroelasticity problem using the Uzawa method with Multigrid:

Unlike the classical Uzawa solver, here we propose Multigrid in the inner step (and in this case, limited to two cycles). In this algorithm, $I t_{\text {int}}$ represents the number of inner smoothings, $i t_{u v}$ denotes the number of smoothings of the displacement variables $u$ and $v$ that are performed within the W-cycle, $\operatorname{res}(p)$ corresponds to the residue of the pressure, and $\omega$ depicts the relaxation parameter that is calculated on the finest grid in the Multigrid cycle, defined for $\tau=1$ by [26]:

$\omega=\frac{2 h^2(\lambda+2 \mu)}{10 K(\lambda+2 \mu)}$    (15)

In this section, we present some numerical experiments aiming to demonstrate the robustness, efficiency, and good performance of the proposed Space-time Multigrid method outlined in the previous section. For this paper, we relied on Wang [1], which utilizes the values of the constants $E$ (Young's modulus) and $K$ (hydraulic conductivity) with realistic values, chosen as follows:

Case 1: $K=1$ and $E=1$ for verification purposes;

Case 2: $K=10^{-9}$ and $E=10^6$ representing a material equivalent to clay;

Case 3: $K=10^{-6}$ and $E=10^4$ representing a material equivalent to silt;

Case 4: $K=10^{-3}$ and $E=10^6$ representing a material equivalent to sand.

For comparison purposes, we used the same relaxation parameter $\omega$, time step $\Delta \mathrm{t}$ and mesh size $h$ parameters as in the study [26].

5.1 Computational details

For the proposed problem, we utilized the W-cycle Multigrid method, employing a full weighting restriction operator and a piecewise constant interpolation operator [15]. The choice of the W-cycle is consistent with reference [37], which asserts that the W-cycle is more robust than the V-cycle for poroelasticity problems. We set $v_1=v_2=v=1$ for both pre- and post-smoothing. For relaxation in the plane, the Uzawa solver is applied with $i t_{M C}=2$ for the variables $u$ and $v$, where $i t_{M C}$ denotes the number of Multigrid cycles. The parallelizable Symmetrical Gauss-Seidel Red-Black method is used. The number of grid levels ($n g l$) for Multigrid was limited to $n g l=13$, i.e., $N_{x y}=N_x \times N_y=2^{13} \times 2^{13}=$ $8192 \times 8192$ volumes in the spatial mesh, and in Singlegrid with $N_{x y}=N_x \times N_y=2^{10} \times 2^{10}=1024 \times 1024$ volumes in the spatial mesh. These limits were set due to memory and computational time usage, respectively. Additionally, the final time was defined as $t_f=10 \mathrm{~s}$. The stopping criterion was based on the reduction of the dimensionless residual from the initial estimate, set to $10^{-8}$, as follows:

$r e s=\left\|r e s_u\right\|_{\infty}+\left\|r e s_v\right\|_{\infty}+\left\|r e s_p\right\|_{\infty}$,     (16)

where $r e s_u$, $res_v$ and $r e s_p$ represent, respectively, the residue obtained in the variables $u, v$, and $p$ [20, 46].

The computational code was implemented using MATLAB R2023b - Academic. The tests were conducted on a computer equipped with a 13th Gen Intel (R) Core (TM) i5-13600K processor running at 3500 Mhz, featuring 14 cores and 20 logical processors, with 64 GB of RAM, and a 64-bit Microsoft Windows 11 Home operating system.

The verification of the code used in this article was carried out by comparing the solutions with those presented in Luo et al. [26] and Franco et al. [30].

The results were obtained using the implicit Euler method, which is unconditionally stable [14], and tested in conjunction with Multigrid.

5.2 Singlegrid and Multigrid comparison

In this work, to compare the performances of the Singlegrid and Multigrid methodologies, the average convergence factor, the order of complexity, and speed-up are employed.

5.2.1 Convergence factor

The convergence factor, which approximates the asymptotic convergence factor, is calculated by [35]:

$\rho=\frac{\left\|\operatorname{res}_{M C}\right\|_{\infty}}{\left\|\operatorname{res}_{M C-1}\right\|_{\infty}}$   (17)

where, $\operatorname{res}_{M C}$ represents the value of the residual in the $MC$-th Multigrid cycle or iteration of the Singlegrid. Thus, the average convergence factor is given by:

$\rho_m=\sqrt[M C]{\frac{\left\|r e s_{M C}\right\|_{\infty}}{\left\|r e s_0\right\|_{\infty}}}$    (18)

where, $r e s_0$ denotes the value of the residual in the initial estimate.

In Figure 4, the comparison of the average convergence factor $\rho_m$ with mesh refinement for the Multigrid method (MG) and Singlegrid (SG) is presented, considering the constants $K$ and $E$ for the 4 cases studied. Here, $n g l$ denotes the number of grid levels in the Multigrid method and represents $2^{{ngl}}$ control volumes in each spatial direction.

4.png

Figure 4. Comparison of $\rho_m$ between Singlegrid and Multigrid for the 4 cases studied

The results obtained for the 4 cases are similar; that is, $\rho_m$ for Singlegrid converges to the constant 1, while for Multigrid, it consistently remains close to 0 [35, 36].

Accurate solutions can be obtained by combining the Multigrid method with Repeated Richardson Extrapolation (RRE) [47].

5.png

Figure 5. Comparison of $\rho_m$ for Multigrid defined for the 4 cases studied

Figure 5 depicts the average convergence factor $\left(\rho_m\right)$ related to $2^{\text {ngl}}$ control volumes for the Multigrid method across the various cases studied.

It is observed that $\rho_m$ for Cases 2, 3, and 4 (realistic values of $K$ and $E$ ) for Multigrid is below 0.1 for the most refined meshes, demonstrating high efficiency. The convergence factor $\rho$ for the first iteration of the Uzawa solver is very efficient, obtaining values tending to 0. However, subsequent iterations do not exhibit the same efficiency, and the average convergence factor $\rho_m$ tends to 0.31 for Case 1, and tends to zero in the other cases studied.

The results of $\rho_m$ in this work are better compared to Luo et al. [26] and Franco et al. [30].

Table 3 presents the number of W-cycles required until the stopping criterion is reached with mesh refinement.

Table 3. Number of W-cycles

From Table 3, it becomes evident that in every case examined, the number of iterations stabilizes with mesh refinement. This is a highly desirable characteristic.

Luo et al. [37] presented the number of cycles W(1,1) ranging from 30 to 45 for different values of hydraulic conductivity, which is three times higher than the number of cycles reported in this paper.

5.2.2 Complexity

Considering the computational time $t_{c p u}$ of the simulations, it is possible to perform a geometric fit to assess the complexity of the algorithm utilizing the equation [14]:

$t_{c p u}=c . \mathbf{N}^\gamma$  （19）

where, $c$ represents a coefficient depending on the method and solver adopted, $\gamma$ denotes the complexity order of the solver associated with the employed method (graphically, it is the slope of the fitting curve), and $\boldsymbol{N}=N_x \times N_y \times N_t$ is the total number of unknowns of the problem. Theoretically, values of $\gamma$ close to 1 and $c$ tending to 0 represent better performances of the employed algorithm. Particularly in Multigrid, they mean that CPU time increases linearly with the mesh size. Table 4 presents the results obtained for all cases studied.

Table 4. Parameters $\gamma$ and c of Eq. (19)

The complexity presented agrees with the findings in the literature [35, 36].

5.2.3 Computational time and speed-up

In this section, our objective is to assess the performance of both the Singlegrid and Multigrid methods, employing computational time $\left(t_{c p u}\right)$ and speed-up $(S)$ as pivotal metrics.

In Figure 6, the computational time $t_{c p u}$, in seconds, is depicted with mesh refinement using both the Singlegrid and Multigrid methods for Cases 1 and 4. It is evident that the Singlegrid time is notably higher compared to Multigrid one and escalates further with refinement. This trend renders it impractical to execute simulations for problems using Singlegrid with more than $N_x \times N_y=2^{10} \times 2^{10}$ control volumes in the mesh.

6.png

Figure 6. CPU Time for Singlegrid (SG) and Multigrid (MG) with mesh refinement

7.png

Figure 7. Speed-up for Case 4

Finally, we evaluate the speed-up (S), which is a metric adopted to measure the increase in the speed of the Multigrid algorithm compared to Singlegrid, given by [48]:

$S=\frac{t_{c p u}(S G)}{t_{\text {cpu }}(M G)}$  （20）

Note that from Figure 7, for example, for Case $4\left(K=10^{-3}\right.$ and $E=10^6$), with $n g l=10$ we have $t_{c p u}(S G)=$ 389514.11 s, approximately 108 hours, $t_{c p u}(M G)=79.79 \mathrm{~s}$, and $S=4881.74$. This implies a substantial advantage of the Multigrid method, as it is approximately 4881 times faster than the Singlegrid, as evident from Figure 7 for Case 4.

Additionally, the notable increasing behavior of the curve indicates the rise of $S$ with the increase of the number of grid levels $n g l$, a highly desirable property. Analogous results were achieved for the other cases studied.

The process used in this paper is the same as presented by Luo et al. [26], applied to the Time Stepping scan, and yielded better results when applied to the new space-time sweep.

In the study conducted by Luo et al. [37], a comparison is made between the Uzawa and Vanka methods, noting that the Uzawa method saves 50% of operations compared to the Vanka method, which justifies the preference for Uzawa. Furthermore, Uzawa is 30% more computationally efficient.

We extend our appreciation to PPGMNE for granting us the opportunity to further our understanding of the subject matter. Additionally, we would like to acknowledge the encouragement and support received from the Pontifical Catholic University of Parana and the Federal University of Parana.