How can parameter-free preconditioning handle nearly incompressible poroelasticity and elasticity?

Why the locking regime makes the linear system nearly singular

The paper studies the Biot model of linear poroelasticity, where solid displacement and fluid pressure are coupled, discretized with the lowest-order weak Galerkin method in space and the implicit Euler scheme in time. Each time step produces a saddle point linear system that must be solved iteratively.

The difficulty is the locking regime: as the Lamé parameter λ grows, the solid becomes nearly incompressible and the divergence of the displacement tends to zero. The divergence term then dominates the first equation, and with pure Dirichlet boundary conditions the resulting algebraic system becomes nearly singular. The paper also notes the related Poisson phenomenon, in which schemes that fail to preserve the correct asymptotic behavior show deteriorating convergence as the mesh is refined.

This matters because a nearly singular system is hard for iterative solvers regardless of how well the discretization avoids locking. The paper therefore separates two questions: obtaining a locking-free discretization, and solving the resulting linear systems efficiently when λ is large.

The regularization comes from an identity, not an added model

The authors introduce a numerical pressure variable and rewrite the original two-field system as a three-field formulation. They then add a term of the form −ρ w wᵀ z_h = 0 to the third block equation, where w is defined in the paper and ρ is a regularization constant.

The key point is that the condition wᵀ z_h = 0 is an inherent identity of the original formulation. Adding a term that vanishes on the exact solution therefore regularizes the algebraic system without perturbing the solution being computed. The paper states that the regularized system preserves the original solution and that the resulting Schur complement has eigenvalues uniformly bounded above and below by positive constants, staying away from zero.

The regularization parameter is chosen as ρ = C λ min(M°_p), with C a constant independent of the mesh size h and the physical parameters. In the numerical experiments the authors use ρ = 0.1 λ min(M°_p), and they also test ρ = 0.01 λ min(M°_p) and ρ = λ min(M°_p). Because the constant is tied to the problem rather than tuned per mesh, the approach is described as parameter-free preconditioning.

Block Schur complement preconditioners become effective again

For the regularized three-field system, the paper applies standard inexact block Schur complement preconditioning. Two variants are constructed: a block diagonal preconditioner P_d and a block triangular preconditioner P_t, each built from approximations of the leading block, a middle block, and a diagonal pressure block.

The analysis shows that MINRES applied to the preconditioned system with P_d, and GMRES applied with P_t, converge essentially independently of the mesh size h, the locking parameter λ, and the combination κΔt, in terms of both convergence factors and asymptotic error constants. The paper states this for the regularized systems and notes that the corresponding proof for the three-field preconditioners is omitted for brevity.

Implementation requires inverting the diagonal blocks. The leading block and the middle block are handled with CG using an incomplete Cholesky decomposition with threshold dropping and a drop tolerance of 10⁻³; the term ρ w wᵀ is omitted when forming that decomposition. The pressure block is diagonal, so its inversion is trivial. The authors report that PCG typically converges in just a few iterations.

Iteration counts stay nearly flat as λ and the mesh change

The 2D and 3D experiments use c₀ = 1, μ = 1, κ ∈ {10⁻⁶, 10⁻³, 1}, λ ∈ {1, 10², 10⁴, 10⁶}, and Δt ∈ {10⁻⁶, 10⁻³, 10⁻²}, with a relative tolerance of 10⁻⁸, at most 1000 iterations, a zero initial guess, and GMRES restart 30.

In the 2D regularized results, GMRES converges in roughly 25 iterations across all tested cases, regardless of λ, Δt, κ, and h. MINRES counts are higher but also stable: for example, at κ = 1 and Δt = 10⁻⁶ the counts move from 44 to 48 as the mesh grows from 918 to 58,608 elements, and from 44 to 58 as λ rises from 1 to 10⁶. The 3D table shows the same pattern at larger counts, with GMRES around 24 to 45 iterations and MINRES around 44 to 88.

The paper notes that the iteration count changes slightly as λ increases from 1 to 10², then remains nearly constant for λ = 10², 10⁴, and 10⁶. The authors also report that the solving time scales as O(N^{3/2}) across mesh sizes, consistent with nearly constant outer iteration counts combined with O(h⁻¹) growth in inner PCG iterations, each costing O(N).

Regularization mainly stabilizes GMRES as λ grows

The paper compares the regularized system with the unregularized one in the 2D example. For λ = 1, the unregularized formulation needs nearly the same number of iterations as the regularized formulation. For larger λ, the unregularized formulation needs slightly more iterations, and the count increases with λ.

The clearest contrast is in GMRES: with regularization it converges in approximately 25 iterations across all tested cases, while without regularization the GMRES count increases with λ and nearly doubles as λ rises from 1 to 10⁶. For MINRES the improvement from regularization is less pronounced, but the counts remain more stable with regularization than without. The paper also reports that GMRES takes slightly less time to converge with regularization than without.

This comparison is the paper's most decision-relevant result: the regularization is not merely a theoretical device, but a change that keeps the outer solver's iteration count from growing with the locking parameter, particularly for the triangular-preconditioned GMRES path.

Jumps in permeability are tested, but the analysis assumes constant κ

A second experiment varies the permeability with a jump: κ = κ₁ on 0.25 ≤ x ≤ 0.75 and κ = 1 elsewhere, with κ₁ ∈ {10⁻¹⁰, 10⁻⁶, 10⁻²}, Δt = 10⁻⁶, λ ∈ {1, 10², 10⁴, 10⁶}, and α ∈ {1, 10⁻², 10⁻⁴}. The paper notes that jumps in κ may occur within elements because triangular meshes are used.

The convergence analysis for MINRES and GMRES applies only to the case of constant κ. The block preconditioners remain applicable algorithmically for variable and possibly discontinuous κ, and the experiment examines whether they still perform well. The iteration counts remain nearly constant as κ₁ varies with other parameters fixed, and stable under mesh refinement. For fixed κ₁, variations in λ and α produce only modest changes: on a mesh with 58,608 elements, MINRES needs 52 iterations when α = 10⁻² and λ = 10⁶, versus 26 iterations when α = 10⁻⁴ and λ = 1.

The authors also examine how inner PCG accuracy affects the outer solver. MINRES iteration counts remain stable when the inner PCG relative tolerance is at most 10⁻⁶, which is looser than the 10⁻⁸ outer tolerance; for robustness they set the inner tolerance to 0.01 times the outer one. The paper identifies extension to higher-order weak Galerkin discretizations and spatially varying parameters as future research.

The priority is a solution-preserving fix, not a new discretization

The paper's central claim is narrower than a general solver recipe: for linear poroelasticity discretized with the lowest-order weak Galerkin method and implicit Euler, the near-singularity in the locking regime can be removed by a regularization drawn from an identity of the original system, and the exact solution is preserved. The payoff is that standard inexact block Schur complement preconditioners then deliver MINRES and GMRES convergence essentially independent of mesh size and the locking parameter.

The evidence that should guide a reader's judgment is the comparison between regularized and unregularized runs. The regularization buys the most for GMRES, holding it near 25 iterations in 2D across λ from 1 to 10⁶, while the unregularized count nearly doubles over that range. MINRES gains less, and at λ = 1 the two formulations are nearly indistinguishable.

The scope conditions are explicit: the convergence analysis covers constant κ, while the jump experiment is numerical; the tests use the lowest-order weak Galerkin method with implicit Euler and Dirichlet boundary conditions for both displacement and pressure; and the regularization constant is chosen as a multiple of λ min(M°_p), with several tested multiples giving comparable counts. Within that scope, the paper's priority is clear: fix the algebraic near-singularity without altering the solution, then let standard block preconditioners carry the solver.

Sources used in this answer