Numerical Verification of PolyDG Algebraic Solvers for the Pseudo-Stress Stokes Problem

This work focuses on the development of efficient solvers for the pseudo-stress formulation of the unsteady Stokes problem, discretised by means of a discontinuous Galerkin method on polytopal grids (PolyDG). The introduction of the pseudo-stress variable is motivated by the growing interest in non-Newtonian flow models and coupled interface problems, where the stress field plays a fundamental role in the physical description. The space-time discretisation of the problem is obtained by combining the PolyDG approach in space with the implicit Euler method for time integration. The resulting linear system, characterised by a symmetric, positive, definite matrix, exhibits deteriorating convergence with standard solvers as the time step decreases. To address this issue, we investigate two tailored strategies: deflated Conjugate Gradient, which mitigates the effect of the most problematic eigenmodes, and collective Block-Jacobi, which exploits the block structure of the system matrix. Numerical experiments show that both approaches yield iteration counts effectively independent of $Δt$, ensuring robust performance with respect to the time step. Future work will focus on extending this robustness to the spatial discretisation parameter $h$ by integrating multigrid strategies with the time-robust solvers developed in this study.

💡 Research Summary

The paper addresses the development of robust linear solvers for the pseudo‑stress formulation of the unsteady Stokes equations discretised with a discontinuous Galerkin method on general polytopal meshes (PolyDG). By introducing the pseudo‑stress variable σ(u,p)=μ∇u−pI, the authors target applications where the stress tensor is a primary quantity, such as non‑Newtonian fluids and coupled interface problems. Spatial discretisation employs PolyDG with element‑wise polynomial degree p≥1 on arbitrary polygonal elements, while temporal integration uses the fully implicit Euler scheme (θ=1). At each time step the semi‑discrete system leads to a symmetric positive‑definite linear system

 A* σ^{n+1}_h = M σ^n_h + Δt f^{n+1}, with A* = M + Δt A,

where M is the mass matrix (positive semi‑definite) and A is the stiffness matrix (positive definite). Because M possesses a non‑trivial kernel, the condition number of A* scales like 1/Δt; consequently, standard Conjugate Gradient (CG) or CG preconditioned with a classical block‑Jacobi diagonal deteriorates dramatically as the time step becomes small.

Two tailored strategies are investigated:

1. Deflated CG – The kernel of M is explicitly characterised and used as a deflation subspace S. An A*-orthogonal projector π_{A*}(S) is built from a basis V of ker(M). The original system is split into a low‑dimensional component that can be solved directly (via VᵀAV) and a deflated system (I−π)A(I−π) x̂ = (I−π)ᵀf*. CG is then applied to the deflated system, which is symmetric positive semi‑definite. The deflation removes the eigen‑modes responsible for the 1/Δt growth, yielding iteration counts essentially independent of Δt. The inversion of VᵀA*V is performed with a direct solver; analytically it behaves like Δt²(B₁+B₃), i.e., a Laplacian‑type operator, suggesting that multigrid could replace the direct step in future work.

2. Collective Block‑Jacobi Preconditioner – Instead of treating each scalar component of σ separately (standard block‑Jacobi), the four components associated with a single element are grouped into a 4×4 block. This “collective” block diagonal matrix is block‑diagonally dominant and captures the strong intra‑element coupling. When used as a preconditioner for CG, it dramatically reduces the Δt‑dependence of the iteration count, achieving near‑constant iteration numbers across a wide range of time steps and mesh refinements.

Numerical experiments are carried out on the unit square with polynomial degree p=3 and a sequence of polytopal meshes containing 50, 100, 200, 400, and 800 elements (h≈0.2462 → 0.0637). Time steps ranging from 10⁻² down to 10⁻¹⁰ are tested. Table 1 shows that the condition number of A* grows proportionally to 1/Δt for small Δt, while the collective Block‑Jacobi preconditioned matrix exhibits a markedly lower condition number. Tables 2 and 3 report iteration counts for four solvers: standard CG, Deflated CG, CG with standard Block‑Jacobi, and CG with collective Block‑Jacobi. Standard CG requires thousands of iterations for Δt≤10⁻⁶ and refined meshes, confirming the ill‑conditioning. Deflated CG reduces the count to a few dozen (≈60–100) regardless of Δt or h. CG preconditioned with the collective Block‑Jacobi also achieves a stable count (≈55–65) across all test cases, whereas the standard Block‑Jacobi preconditioner still shows a noticeable Δt‑dependence.

The results demonstrate that both the deflation technique and the collective block‑Jacobi preconditioner effectively neutralise the time‑step induced deterioration of the linear system. Consequently, the solvers exhibit iteration counts that are essentially independent of Δt, fulfilling the primary goal of time‑step robustness. The authors conclude by outlining future work: integrating multigrid methods to efficiently invert the deflation correction (VᵀA*V) and extending the robustness to the spatial discretisation parameter h, thereby achieving a fully time‑ and space‑robust solver for pseudo‑stress PolyDG discretisations.