mass matrix, finite elements, additive Schwarz, block preconditioning, Schur complement, algebraic multigrid

1 Introduction↩︎

Let \(V_h\) be a conforming finite-element space with basis \(\{\phi_i\}\). The mass matrix \(M_{ij} = (\phi_j,\phi_i)_\Omega\) is the coefficient representation of the \(L^2(\Omega)\) Riesz map on \(V_h\). Inverting or approximately inverting mass matrices arises often in finite-element simulation of partial differential equations (PDEs), both in explicit time integration of dynamic PDEs and in block preconditioning for mixed systems of PDEs. For conforming finite-element spaces, the mass inverse is typically dense, and although often considered one of the least interesting aspects of finite-element simulations, it can have a significant effect on total computational cost and wallclock time. For explicit integration, applying the mass inverse can be a dominant computational cost, and the only globally coupled component across all spatial cells and MPI cores. In block preconditioning, the local/sparse approximation of a mass inverse within an approximate Schur complement can determine the quality of the Schur complement approximation. This in turn largely determines the convergence of preconditioned fixed-point and Krylov iterations; see, e.g., [1], [2].

One line of research avoids the need for mass inverse approximations or preconditioners by directly using lumped (diagonal) mass matrices in the discretization. Analyzing such an approach becomes a question of quadrature and functional approximation [3], [4] rather than algebraic and spectral properties of the assembled matrix. With sufficient assumptions on mesh and finite-element space, mass lumping can work well, though it is not always viable, with higher-order discretizations [5], [6], hyperbolic equations [7], Sobolev spaces other than \(H^1\) [8], [9], and non-tensor-product or curvilinear meshes taking extra care. Recent work has also considered deferred correction approaches to avoid mass solves [10].

Relative to solvers for differential operators, little attention has been given to fast, general-purpose solvers for mass matrices; see related commentary in [11]. In fact, diagonal (Jacobi) approximations are commonplace in practice, including in many state-of-the-art codes. There are likely several reasons for this. Firstly, diagonal approximations are simple and cheap, and in many cases they are simply “good enough,” with many conforming elements enjoying a reasonable spectral equivalence between the mass and its diagonal, independent of mesh resolution. Moreover, a diagonal approximation is natural in the block-preconditioning setting, where maintaining a sparse approximate Schur complement is typically important; in fact, a diagonal inverse in this setting can perform better than its standalone mass-preconditioning behavior would suggest.⁵ Lastly, mass matrices typically do not have a very nuanced structure, and this likely inhibits more effective approximations, like multilevel solvers.

Although mass preconditioning literature is relatively small, existing work gives several important points of comparison. Wathen [12] showed that continuous Galerkin mass matrices are spectrally equivalent to their diagonal, independent of mesh resolution, with spectral bounds built from unassembled element matrices. This relation to element matrices leads to independence of mesh resolution, but not finite-element order [13]. Recent work has focused on low-order-refined algorithms that precondition high-order conforming finite-element discretizations with appropriately defined low-order discretizations [14], [15]. Such an approach has demonstrated strong results, but requires a more involved implementation, especially if one wanted to use it in the context of block preconditioning. Specialized \(p\)-robust mass preconditioners for scalar \(H^1\)-conforming spaces have been developed that rely on carefully constructed high-order bases [11], [16], [17], and fast factorization techniques have also been developed for high-order Bernstein local element matrices [18]. A special interpolation–histopolation basis is also used in [19], wherein high-order Raviart–Thomas mass matrices are then spectrally equivalent to their diagonal, independent of polynomial degree, and this approximation is used to construct Schur complement approximations in solving \(H(\operatorname{div})\) saddle-point problems.

In this paper we show that there is significant room for improvement over diagonal mass inverse approximations, particularly for higher-order discretizations and for approximation of Sobolev spaces other than \(H^1\). In contrast to the specialized and basis-function oriented approaches mentioned above, here we propose a largely algebraic, and implementationally basis-function agnostic, framework for sparse approximation of mass inverses that applies to any conforming finite-element space, and is naturally suited for both efficient explicit time integration and approximate Schur complements. Our approach is based on a sparse approximate inverse formed by applying exact element-local inverse mass matrices on the underlying broken or unassembled finite-element space and averaging the result back to the conforming space. The resulting BRoken-space Additive Schwarz (BRAS) algorithm yields a mass inverse approximation in the sparsity pattern of the original mass matrix, and provides a significantly better approximation than the matrix diagonal or scalings thereof. Its construction only requires the same element-level information needed to assemble the mass matrix, namely the element mass matrices and the local-to-global DOF map.

The remainder of this paper is structured as follows. 2 introduces the BRAS preconditioner and provides supporting theory for it. 3 studies BRAS for standalone mass systems \(M\boldsymbol{u}=\boldsymbol{f}\), reporting condition numbers, solver iteration counts, and solve times. Despite the BRAS inverse being more expensive to apply than a diagonal inverse, in all cases tested the reduced iteration counts provided by BRAS offset its higher apply cost, yielding 1.5–6\(\times\) speedups in preconditioned conjugate gradient (CG) solve time across \(H^1\), \(H(\operatorname{curl})\), and \(H(\operatorname{div})\) finite elements on two- and three-dimensional simplicial meshes, and finite element orders \(p\in[1,4]\). Aspects of using approximate mass inverses in block preconditioning are considered in 4, with accompanying numerical results given in 5. For mixed Poisson and biharmonic problems, block preconditioning based on Schur complements built from BRAS inverses yield 1.5–3\(\times\) speedups over those using standard diagonal mass inverse approximations. Concluding remarks are provided in 6.

2 BRAS: Broken-space additive Schwarz mass inverse↩︎

2.1 Conforming mass matrix inverse↩︎

Let \(\Omega\subset\mathbb{R}^d\) be a bounded domain and let \(\mathcal{T}_h\) be a finite-element mesh. Let \(V_h\) be a conforming finite-element space, scalar- or vector-valued, with global conforming basis \(\{\phi_i\}_{i=1}^n\). The corresponding mass matrix is the \(L^2(\Omega)\) Gram matrix \[\label{eq:mass95matrix95def} M_{ij} = (\phi_j,\phi_i)_{L^2(\Omega)} = \int_\Omega \phi_j(x)\cdot\phi_i(x)\,\,\mathrm{d}x, \qquad i,j=1,\ldots,n.\tag{1}\] For scalar-valued spaces, the dot product denotes ordinary multiplication; for vector-valued spaces, it denotes the Euclidean dot product of the vector-valued basis functions. The matrix \(M\) represents the discrete \(L^2\) Riesz map. If \(u_h=\sum_j (\boldsymbol{u})_j\phi_j\), then the coefficient vector \(M\boldsymbol{u}\) represents the functional \(v_h \mapsto (u_h,v_h)_\Omega.\) Thus \(M^{-1}\) maps a coefficient representation of a discrete \(L^2\) functional back to the coefficient vector of the corresponding finite-element function. The inverse is global even though the bilinear form itself is elementwise local.

Given a symmetric positive definite (SPD) inverse approximation \(\widetilde{M}^{-1}\approx M^{-1}\), the implied approximate mass matrix is \({\widetilde{M}}\), and the natural measure of approximation is the generalized eigenproblem \(M\boldsymbol{v} = \lambda {\widetilde{M}}\boldsymbol{v}.\) If the generalized eigenvalues \(\{\lambda\}\) are contained in a small interval near unity, this implies rapid convergence of fixed-point or preconditioned CG applied to \(M\boldsymbol{u}=\boldsymbol{f}\), and also controls the conditioning of Schur complement approximations involving \(\widetilde{M}^{-1}\), as discussed in [sec:block_preconditioning] [sec:block_preconditioning_numerics]. The baseline approximation for a mass preconditioner and approximate inverse is the diagonal inverse \(M_{\rm diag}^{-1} = \mathop{\mathrm{diag}}(M)^{-1}\) or slight modifications or weightings thereof. This is the cheapest SPD approximation of \(M^{-1}\), but ignores all off-diagonal coupling, which is an increasing fraction of the connections in \(M\) for increasing \(p\).

2.2 Broken-space mass inverse↩︎

Let \(V_h\) be a conforming finite-element space. Define the broken companion space \[V_h^{\rm br}= \bigoplus_{K\in\mathcal{T}_h} V(K),\] where \(V(K)\) is the local finite-element space on the element \(K\). Let \(M_{\rm br}\) be the mass matrix on \(V_h^{\rm br}\). Since \(V_h^{\rm br}\) is broken over the mesh, \(M_{\rm br}\) is block diagonal: \[M_{\rm br}= \operatorname{diag}_{K\in\mathcal{T}_h} M_K,\] where \(M_K\) is the local element mass matrix. The conforming mass matrix is then assembled from broken pieces. Let \(n=\dim V_h\) and \(n_{\rm br}=\dim V_h^{\rm br}\). Each broken DOF is associated with a unique conforming DOF; denote this map by \(\pi:\{1,\ldots,n_{\rm br}\}\to\{1,\ldots,n\}\). Let \[P:V_h\to V_h^{\rm br}\] be the conforming-to-broken injection defined in coefficient form by \[(P\boldsymbol{u})_j = \boldsymbol{u}_{\pi(j)}, \qquad j=1,\ldots,n_{\rm br}.\] Equivalently, \(P_{ji}=1\) if \(i=\pi(j)\), and \(P_{ji}=0\) otherwise. Thus \(P\) represents a conforming function in the broken basis by copying each conforming DOF to all of its element-local broken copies. An illustrative example is shown in 1. Then the conforming mass matrix is assembled as \[\label{eq:M95PMbrP} M = P^{\top}M_{\rm br}P = \sum_{K\in\mathcal{T}_h} P_K^{\top}M_KP_K,\tag{2}\] where \(P_K\) denotes the local interpolation to the element \(K\).

The BRAS inverse approximation we present in this paper consists of inverting local element matrices in the broken space, and mapping them back to the global conforming space with a careful weighting by multiplicity: \[\label{eq:mhatinvbras95def} \boxed{ M_{\rm bras}^{-1} = (P^\top P)^{-1}P^{\top}M_{\rm br}^{-1}P (P^\top P)^{-1}. }\tag{3}\] Here \(D := P^\top P \in \mathbb{R}^{n \times n}\) is a diagonal matrix, whose \(i\)th entry counts the number of broken copies of the \(i\)th conforming DOF, \[D=P^{\top}P=\mathop{\mathrm{diag}}(d_1,\ldots,d_n), \qquad d_i=\#\{j:\pi(j)=i\}.\] \[Note that the inverse \(M_{\rm br}^{-1} = \operatorname{diag}_{K\in\mathcal{T}_h} M_K^{-1}\) is well defined, since element-local mass matrices are SPD. Thus BRAS is assembled analogously to the conforming mass matrix itself, replacing each local mass block by its inverse and applying multiplicity scaling. As such, in principle any finite-element library would be able to explicitly assemble BRAS using the same kernels as it uses to assemble conforming mass matrices. Note that for full column-rank P, it follows immediately from \eqref{eq:mhatinvbras95def} that M_{\rm bras}^{-1} is SPD. The BRAS inverse can also be written as a symmetrically diagonally scaled elementwise sum, \label{eq:mass95and95bras95element95sums} \begin{align} M_{\rm bras}^{-1} &= D^{-1} \left( \sum_{K\in\mathcal{T}_h} P_K^{\top}M_K^{-1}P_K \right) D^{-1}. \end{align}\tag{4}\] Here we see that \(M_{\rm bras}^{-1}\) has the same nonzero pattern as \(M\). Therefore one application of the assembled BRAS matrix has the same local dense matrix-vector work as one application of the assembled mass matrix.

2.3 Theory↩︎

We define two SPD (or symmetric negative definite) matrices \(A,B\in\mathbb{R}^{n\times n}\) in the matrix pencil \((A,B)\) to be spectrally equivalent if there exist constants \(0<\alpha \leq \beta\) such that for all \(\boldsymbol{u}\in\mathbb{R}^n\backslash \{\boldsymbol{0} \}\) \[\label{eq:spec-equiv} \alpha \leq \frac{\langle A\boldsymbol{u}, \boldsymbol{u}\rangle }{\langle B\boldsymbol{u}, \boldsymbol{u}\rangle } \leq \beta.\tag{5}\] Let \(\beta\) and \(\alpha\) be the maximum and minimum constants such that 5 holds, respectively, corresponding to the maximum and minimum generalized eigenvalues of the matrix pencil \((A,B)\). Then we define the spectral condition number (distinct from the more common \(\ell^2\)-condition number) of the matrix pencil \((A,B)\) as \[\label{eq:spectral-kappa} \widehat{\kappa}(A,B) \mathrel{\vcenter{:}}= \frac{\beta}{\alpha} = \frac{\lambda_{\max}(B^{-1}A)}{\lambda_{\min}(B^{-1}A)},\tag{6}\] where \(\widehat{\kappa}(A,B) = \widehat{\kappa}(B,A)\). In this section we prove a sharp bound on \(\alpha\) for the BRAS-preconditioned mass operator; as discussed in 3, this is useful in practice because certain fixed-point iterations require reliable spectral intervals. Moreover, we characterize \(\beta\) through a broken-space projection problem; numerical results in 3 suggest that \(M_{\rm bras}^{-1}\) may give spectral condition numbers bounded independently of, or weakly dependent on, \(p\) in certain cases, although we do not prove such a bound in this work.

Consider the conforming mass solve \(M\boldsymbol{u}=\boldsymbol{f}\) for \(\boldsymbol{u} \in \mathbb{R}^{n}\). By 2 , \(M=P^{\top}M_{\rm br}P\), so this problem can also be posed in terms of broken coefficients \(\boldsymbol{w} \in \mathbb{R}^{n_{\rm br}}\) satisfying \(P^\top M_{\rm br}\boldsymbol{w} = \boldsymbol{f}\) subject to the constraint that \(\boldsymbol{w} \in \operatorname{range}(P)\). Thus, let \({\cal C} \in \mathbb{R}^{m \times n_{\rm br}}\) be a full-row-rank constraint matrix on the broken coefficient space such that \(\operatorname{ker}({\cal C})=\operatorname{range}(P)\), wherein \({\cal C}\boldsymbol{w}=\boldsymbol{0}\) exactly enforces that the broken copies are consistent with a single conforming coefficient vector. Here the number of constraints is \(m = n_{\rm br} - n\). Define the broken right-hand side as \(\boldsymbol{f}_{\rm br} = P D^{-1}\boldsymbol{f}\), where \(D=P^{\top}P\), and note that this lift has the correct assembled action because \(P^{\top}P D^{-1}\boldsymbol{f}=\boldsymbol{f}\). The conforming solve \(M\boldsymbol{u}=\boldsymbol{f}\) can therefore be written as the constrained broken-space saddle-point problem \[\label{eq:broken95kkt} \begin{bmatrix} M_{\rm br}& {\cal C}^{\top} \\ {\cal C} & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{w} \\ \boldsymbol{\mu} \end{bmatrix} = \begin{bmatrix} \boldsymbol{f}_{\rm br} \\ \boldsymbol{0} \end{bmatrix}.\tag{7}\] Indeed, the constraint gives \(\boldsymbol{w}=P\boldsymbol{u}\) for some conforming vector \(\boldsymbol{u}\). Multiplying the first block row of 7 by \(P^{\top}\) eliminates the multiplier term, since \({\cal C}P=0\), and gives \(P^{\top}M_{\rm br}P\boldsymbol{u}=\boldsymbol{f}\), which is exactly the conforming system. The following Lemma uses this formulation to directly relate \(M_{\rm bras}^{-1}\) to \(M^{-1}\).

The lower bound in ?? can be sharp; numerically, we find that equality holds in many cases, and even when it does not we often have \(\lambda_{\min} \approx 1\). Let \(\Delta_{\cal C}:=-Z_{\cal C}S_{\cal C}^{-1}Z_{\cal C}^{\top}\), so that \(M_{\rm bras}^{-1}=M^{-1}+\Delta_{\cal C}\) and \(\Delta_{\cal C}\succeq0\). Then \(\lambda_{\min}(M_{\rm bras}^{-1}M)=1\) precisely when the correction \(\Delta_{\cal C}\) has a non-trivial nullspace, and there are at least two mechanisms by which this can occur, with the first being purely dimensional. Since \(-S_{\cal C}^{-1}\) is SPD, \(\operatorname{ker}(\Delta_{\cal C})=\operatorname{ker}(Z_{\cal C}^{\top})\). Since \(Z_{\cal C}^{\top}\in\mathbb{R}^{m\times n}\), rank-nullity implies that \(\operatorname{ker}Z_{\cal C}^{\top}\) is non-trivial whenever \(m<n\). Recalling that \(m=n_{\rm br}-n\), this condition is equivalent to \(n_{\rm br}<2n\). To make this more concrete, if the \(i\)th conforming DOF has \(d_i\geq1\) broken copies, then \[n_{\rm br}=\sum_i d_i, \quad n=\sum_i 1 \qquad \Longrightarrow \qquad m-n = n_{\rm br}-2n = \sum_i(d_i-2).\] Hence, degrees of freedom (DOFs) with multiplicity \(d_i = 1\) contribute negatively towards \(m - n\), DOFs with multiplicity \(d_i = 2\) contribute neutrally towards \(m - n\), and DOFs with \(d_i > 2\) contribute positively towards \(m - n\). Natural candidates for \(m<n\) are therefore spaces whose DOFs are concentrated on entities with low sharing multiplicity: cell-interior and boundary-facet DOFs have \(d_i=1\), while interior-facet DOFs (edges in 2D, faces in 3D) typically have \(d_i=2\). This makes \(H(\operatorname{div})\) a natural candidate, where DOFs are typically on facets and cell interiors. By contrast, low-order \(H^1\) spaces in two and three dimensions often place DOFs on shared vertices, for which \(d_i>2\) typically, so the dimensional argument may fail.

A second mechanism for \(\lambda_{\min}(M_{\rm bras}^{-1} M)=1\) occurs when \(M_{\rm br}P\boldsymbol{1} \in\operatorname{range}(P)\). Let \(Q:=P D^{-1}P^{\top}\) be the \(\ell^2\)-orthogonal projector onto \(\operatorname{range}(P)\) (see also 1 below). If \(M_{\rm br}P\boldsymbol{1}\in\operatorname{range}(P)\), then \(QM_{\rm br}P\boldsymbol{1}=M_{\rm br}P\boldsymbol{1}\), and hence \[M_{\rm bras}^{-1}M\boldsymbol{1} = D^{-1}P^{\top}M_{\rm br}^{-1}QM_{\rm br}P\boldsymbol{1} = D^{-1}P^{\top}M_{\rm br}^{-1}M_{\rm br}P\boldsymbol{1} = D^{-1}P^{\top}P\boldsymbol{1} = \boldsymbol{1} .\] A sufficient condition for \(M_{\rm br}P\boldsymbol{1}\in\operatorname{range}(P)\) is that \(M_K\boldsymbol{1}_K=c\boldsymbol{1}_K\) with \(c\) independent of \(K\), since then \(M_{\rm br}P\boldsymbol{1}=cP\boldsymbol{1}\). This arises, for example, for scalar \(P_1\) nodal elements on uniform affine simplicial meshes, for which \(c=|K|/(d+1)\) is element-independent.

The viewpoint described above surrounding 7 is also related to domain decomposition methods such as finite-element tearing and interconnecting (FETI) and balancing domain decomposition by constraints (BDDC), which introduce duplicated interface variables and constraints between nonoverlapping subdomains; see [20]. These methods are also closely connected to the algebraic hybridization approach of Dobrev et al. [21].

2.4 Piola mappings and metric-dependent BRAS bounds↩︎

Let \(F_K:\widehat K\to K\) denote a sufficiently smooth, orientation-preserving map from a reference cell \(\widehat K\) to a physical mesh cell \(K\), with Jacobian matrix \(J_K\) and corresponding determinant \(\det (J_K)\). We write \(\boldsymbol{x} = F_K(\widehat{\boldsymbol{x}})\), and suppress the dependence of \(J_K\) on \(\widehat{\boldsymbol{x}}\). Element mass matrices are computed by mapping physical-cell \(L^2(K)\) inner products back to reference-cell integrals. Reference-cell integrands naturally inherit a factor of \(\det (J_K)\) due to this change of variables, \(\,\mathrm{d}\boldsymbol{x} = \det(J_K) \,\mathrm{d}\boldsymbol{\widehat{x}}\). Furthermore, the reference-cell integrand depends on the finite-element space, because the reference-to-physical map must preserve the trace quantity used to impose conformity; see [22]. For \(H(\mathop{\mathrm{curl}})\) elements, continuity of tangential traces is enforced, leading to the covariant Piola map, \(\boldsymbol{\Phi}( \boldsymbol{x} ) = J_K^{-\top} \widehat{\boldsymbol{\Phi}}( \widehat{\boldsymbol{x}} )\), with \(\boldsymbol{\Phi}(\boldsymbol{x})\) a vector field on the physical cell and \(\boldsymbol{\widehat{\Phi}}(\boldsymbol{\widehat{x}})\) its analog on the reference cell. For \(H(\relax)\) elements, continuity of normal fluxes is enforced, leading to the contravariant Piola map, \(\boldsymbol{\Phi}( \boldsymbol{x} ) = \tfrac{1}{\det (J_K)} J_K \widehat{\boldsymbol{\Phi}}( \widehat{\boldsymbol{x}} ).\)

Applying the above change of variables, we find that corresponding physical mass inner products on a single element pull back as follows: \[\tag{8} \begin{align} (u,v)_{L^2(K)} &= \int_{\widehat K} \widehat u\, \det (J_K)\, \widehat v\,\,\mathrm{d}\widehat x, && u,v\in H^1(K), \tag{9} \\ (\boldsymbol{u},\boldsymbol{v})_{L^2(K)} &= \int_{\widehat K} \widehat{\boldsymbol{u}}^{\top} \left(\det (J_K)J_K^{-1}J_K^{-\top}\right) \widehat{\boldsymbol{v}}\,\,\mathrm{d}\widehat{\boldsymbol{x}}, && \boldsymbol{u},\boldsymbol{v}\in H(\mathop{\mathrm{curl}};K), \tag{10} \\ (\boldsymbol{u},\boldsymbol{v})_{L^2(K)} &= \int_{\widehat K} \widehat{\boldsymbol{u}}^{\top} \left(\frac{J_K^\top J_K}{\det (J_K)}\right) \widehat{\boldsymbol{v}}\,\,\mathrm{d}\widehat{\boldsymbol{x}}, && \boldsymbol{u},\boldsymbol{v}\in H(\relax;K). \tag{11} \end{align}\] Thus scalar \(H^1\) mass matrices depend on the element map only through the scalar Jacobian factor \(\det (J_K)\), while \(H(\mathop{\mathrm{curl}})\) and \(H(\relax)\) mass matrices depend on tensor-valued Piola metric weights.

The following results compare these metric effects algebraically. We use \(G_K\) to denote the scalar or tensor coefficient that appears in the pulled-back reference-cell mass integrand in 8 . Let \(M^I\) denote the mass matrix assembled on the reference cells with the identity coefficient, and let \(M^G\) denote the mass matrix assembled on the same reference finite element space with coefficients \(G_K\). Equivalently, \(M^G\) is the physical mass matrix on the mapped mesh. The matrices \(M_{\rm br}^I\) and \(M_{\rm br}^G\) denote the corresponding unassembled mass matrices on the associated broken space.

In particular, for the pullbacks in 8 , the metric contrast \(\beta/\alpha\) may be bounded by \[\frac{\beta}{\alpha} \leq \begin{cases} \displaystyle \frac{\sup_{K\in\mathcal{T}_h} \det (J_K)}{\inf_{K\in\mathcal{T}_h} \det (J_K)}, & H^1, \\[1.25ex] \displaystyle \frac{\sup_{K\in\mathcal{T}_h}\lambda_{\max}(\det (J_K)J_K^{-1}J_K^{-\top})}{\inf_{K\in\mathcal{T}_h}\lambda_{\min}(\det (J_K)J_K^{-1}J_K^{-\top})}, & H(\mathop{\mathrm{curl}}), \\[2ex] \displaystyle \frac{\sup_{K\in\mathcal{T}_h}\lambda_{\max}(J_K^\top J_K/\det (J_K))}{\inf_{K\in\mathcal{T}_h}\lambda_{\min}(J_K^\top J_K/\det (J_K))}, & H(\relax). \end{cases}\]

3 Numerical results: Mass matrices↩︎

This section studies solve performance for the mass system \(M\boldsymbol{u}=\boldsymbol{f}\) using preconditioned CG and damped fixed-point iterations. Although consistently slower than CG as a standalone solver, fixed-point iterations can still be useful to avoid the parallel communication costs associated with CG dot products or inside of larger coupled block preconditioners. In addition, our interest is not only in solving an isolated mass system, but also in using BRAS as an approximate inverse inside larger preconditioners; see [sec:block_preconditioning] [sec:block_preconditioning_numerics]. In that setting, the standalone condition number, and the fixed-point behavior it controls, are often a cleaner measure of preconditioner quality than CG iterations.

Let \(\widetilde{M}^{-1}\) denote the inverse approximation used as a preconditioner. Given a damping parameter \(\omega>0\), the fixed-point iteration is \[\boldsymbol{u}_{k+1} = \boldsymbol{u}_k + \omega \widetilde{M}^{-1} \left(\boldsymbol{f} - M\boldsymbol{u}_k\right), \qquad \widetilde{M}^{-1} \in \bigl\{ M_{\rm diag}^{-1}, \, M_{\rm bras}^{-1} \bigr\}.\] The fixed-point iteration is convergent for any \(\omega\in(0,2/\lambda_{\max}(\widetilde{M}^{-1}M))\). The choice minimizing the spectral radius of the error propagation matrix \(I-\omega\widetilde{M}^{-1}M\) is \[\omega_{\rm opt}(\widetilde{M}^{-1}) = \frac{2}{\lambda_{\min}(\widetilde{M}^{-1}M)+\lambda_{\max}(\widetilde{M}^{-1}M)} .\] Thus any convergent choice requires an estimate of \(\lambda_{\max}(\widetilde{M}^{-1}M)\), while the optimal choice also requires an estimate of \(\lambda_{\min}(\widetilde{M}^{-1}M)\). In general, maximum eigenvalues are much cheaper to estimate than minimum eigenvalues: power and Arnoldi-type iterations are usually sufficient for the former, while the latter often require inverse iterations and therefore solves with \(M\) shifted by a diagonal matrix. This motivates the following pragmatic weight choices used in the experiments here: \[\omega_{\rm diag} = \frac{1}{\lambda_{\max}(M_{\rm diag}^{-1} M)}, \qquad \omega_{\rm bras} = \frac{2}{1+\lambda_{\max}(M_{\rm bras}^{-1} M)} .\] For BRAS this choice often yields the optimal weight or not far from it, using the lower spectral bound \(\lambda_{\min}(M_{\rm bras}^{-1} M) \geq 1\) from 1 that we note is generally fairly tight. For the diagonal method we use only the upper spectral estimate, since we are unaware of any analogous lower spectral estimate. The largest eigenvalue in both weights is estimated numerically from the corresponding preconditioned operator using the power iteration.

All finite elements are implemented using the Firedrake library [23], and we use Firedrake element-family names in the tables and text. For \(H^1\), “Lagrange” and “Bernstein” denote \(H^1\)-conforming piecewise-\(P_p\) spaces represented using nodal and Bernstein bases, respectively; the nodal basis uses Firedrake’s spectral variant, based on Gauss–Lobatto–Legendre (GLL) points. For \(H(\operatorname{curl})\), “\(\mathrm{N1curl}\)” and “\(\mathrm{N2curl}\)” denote the first- and second-kind Nédélec elements, respectively, and for \(H(\operatorname{div})\) “\(\mathrm{RT}\)” and “\(\mathrm{BDM}\)” denote the Raviart–Thomas and Brezzi–Douglas–Marini elements, respectively. All meshes use simplices (triangles in 2D and tetrahedra in 3D). For solve tests, the number of mesh elements is chosen so that the conforming problem has approximately \(10^6\) DOFs, while spectral condition number calculations use a smaller problem with approximately \(5\times 10^3\) DOFs to enable direct calculation. Solves are done on a random right-hand side \(\boldsymbol{f}\). Both CG and fixed-point iterations are initialized with \(\boldsymbol{u}_0=\boldsymbol{0}\) and are terminated when the residual \(\boldsymbol{r}_k=\boldsymbol{f}-M\boldsymbol{u}_k\) satisfies \(\|\boldsymbol{r}_k\|_2/\|\boldsymbol{r}_0\|_2 \le 10^{-10}\), or after \(5000\) iterations. Solve times do not include construction of the preconditioner or estimation of \(\lambda_{\max}\).

1.8pt

1.8pt

1.8pt

[tab:hdiv_mass_solve,tab:hcurl_mass_solve,tab:h1_mass_solve] show spectral condition numbers, solve times, and iteration counts of \(M_{\rm diag}\) and \(M_{\rm bras}\) for two sets of basis functions associated with \(H^1, H(\operatorname{curl})\) and \(H(\operatorname{div})\), respectively, and finite-element orders \(p\in[1,4]\). We report both solve time and iteration counts due to the different cost of applying \(M_{\rm diag}^{-1}\) vs. \(M_{\rm bras}^{-1}\), with the latter costing the same as application of \(M\). To compare, we consider speedup \(S\) as the ratio of the \(M_{\rm diag}\)-based solve execution time to the \(M_{\rm bras}\)-based solve execution time, and performance improvement \(I\) as the percentage by which the solve is faster, where \(I = 100(S-1)\). Throughout this paper, all reported solve times are mean values over three solves. The tables show that BRAS reduces the spectral condition number across the tested spaces, with the largest relative improvements occurring in high-order spaces. Consistent with this, across all scenarios tested, preconditioning with \(M_{\rm bras}^{-1}\) yields meaningful speedups over \(M_{\rm diag}^{-1}\). For CG solves, almost all cases see at least a 40% speedup, with certain instances such as 4th-order Bernstein bases seeing a \(4\times\) speedup (or \(\approx 300\%\) improvement in performance). Fixed-point improvements are even more dramatic, with \(M_{\rm bras}^{-1}\) yielding at least a \(3\times\) speedup (or \(\approx 200\%\) improvement in performance) in almost all cases, and reaching as high as 20–40\(\times\) faster in certain cases of 4th-order Bernstein and BDM elements. In several cases the BRAS fixed-point solve is actually faster than diagonal-preconditioned CG. We also note that in four of the six \(p=4\) tests in 3D, the fixed-point iteration did not converge in 5000 iterations with \(M_{\rm diag}^{-1}\), while \(M_{\rm bras}^{-1}\) yielded convergence on the order of 100 iterations.

2 shows the speedup provided by BRAS by plotting the solve-time ratio \(t_{\rm diag}/t_{\rm bras}\) for the CG and fixed-point runs in each table. In all cases considered, BRAS yields faster CG and fixed-point solve times than diagonal preconditioning, starting at 30–50% improvements for \(p=1\), and generally increasing with \(p\), significantly in some cases. We again point out that BRAS improvements in the fixed-point setting are much stronger than those for CG.

Wathen and Rees [24] observed that, for diagonally preconditioned low-order \(H^1\) mass matrices, Chebyshev iteration gives convergence very similar to CG when accurate spectral bounds are available [24]. We observe the same behavior in the broader setting considered here, including for both diagonal and BRAS preconditioning. In particular, across essentially all of our tests, the observed CG convergence rates are remarkably close to the standard worst-case estimate \[\frac{\|\boldsymbol{e}_k\|_{A_{{\widetilde{M}}}}}{\|\boldsymbol{e}_0\|_{A_{{\widetilde{M}}}}} \leq 2 \left( \frac{\sqrt{\widehat{\kappa}}-1}{\sqrt{\widehat{\kappa}}+1} \right)^k,\] with \(A_{{\widetilde{M}}}={\widetilde{M}}^{-1/2}M{\widetilde{M}}^{-1/2}\). This is precisely the Chebyshev upper bound obtained from the spectral interval of the preconditioned operator. This makes Chebyshev iteration an attractive alternative to CG in practice, because it is linear, and it avoids the global dot products required by CG, which can be advantageous on parallel machines. To our knowledge, this observation has not previously been emphasized for high-order \(H^1\), \(H(\mathop{\mathrm{curl}})\), and \(H(\relax)\) mass matrices.

Next we recall the theory of 2.4, which considered mass conditioning under mesh deformation. To provide numerical support for 2 therein, we consider preconditioning an \(H(\relax)\) mass on a mesh obtained by rotating and stretching the unit-square simplicial mesh. The Jacobian determinant in this example satisfies \(\det (J_K)=1\) on every cell, so that it isolates the effect of the tensor-valued metric weight appearing in the Piola pullback for the \(H(\relax)\) mass matrix. Given a stretch factor \(r\geq 1\), we consider the coordinate map \[\begin{align} \label{eq:rot95stretch95map} F_r(\boldsymbol{x}) = \boldsymbol{c} + R_{\pi/4} \begin{pmatrix} r & 0 \\ 0 & r^{-1} \end{pmatrix} R_{\pi/4}^{\top} (\boldsymbol{x}-\boldsymbol{c}), \qquad \boldsymbol{c}=(1/2,1/2), \end{align}\tag{12}\] where \(R_{\pi/4}\) denotes rotation by \(\pi/4\).

The Jacobian for this affine map is \(J_K=DF_r\), with \(\det (J_K)=1\), and the singular values of \(J_K\) are \(r\) and \(r^{-1}\). Therefore the \(H(\relax)\) Piola metric \(G_K=J_K^\top J_K/\det (J_K)\) has eigenvalues \(r^2\) and \(r^{-2}\), so the metric contrast appearing in 2 satisfies \[\frac{\beta}{\alpha} \leq \frac{\sup_{K\in\mathcal{T}_h}\lambda_{\max}(J_K^\top J_K/\det (J_K))}{\inf_{K\in\mathcal{T}_h}\lambda_{\min}(J_K^\top J_K/\det (J_K))} = \frac{\lambda_{\max}(J_K^\top J_K)}{\lambda_{\min}(J_K^\top J_K)} = r^4 .\] Thus, the values \(r=1,3,5,7\) measured in the experiment below correspond to metric contrasts of \(1,81,625,\) and \(2401\), respectively.

The general estimate in 2 is stated for \(\lambda_{\max}\), not directly for the spectral condition number. However, in this particular example the minimum eigenvalue satisfies \(\lambda_{\min}(M_{\rm bras}^{-1}M)=1\) on both the reference and stretched meshes, for the reasons discussed after 1. Hence the spectral condition number bound stated in 2 also applies in this case. The resulting spectral condition numbers are shown in 3 for polynomial degrees \(p=1,3,5,7\). The BRAS spectral condition numbers grow essentially like \(r^4\), consistent with the metric-contrast factor \(\beta / \alpha \leq r^4\) above; thus, BRAS does not remove the Piola metric dependence. However, at each fixed value of \(r\), the BRAS spectral condition numbers are substantially smaller than their diagonal counterparts, and the BRAS curves remain nearly flat as the polynomial degree is increased. By contrast, the diagonal spectral condition numbers deteriorate substantially with \(p\), roughly doubling with each increment of \(p\) shown.

4 Block preconditioning↩︎

We now show why the quality of a mass inverse approximation matters when the inverse appears inside a Schur-complement preconditioner. Consider the nonsingular symmetric saddle-point matrix \[\label{eq:saddle95matrix95general9522} \mathcal{A} = \begin{bmatrix} M & G^{\top} \\ G & -C \end{bmatrix},\tag{13}\] where \(C\) is SPSD and \(M\) is an SPD mass matrix. This type of matrix is commonplace throughout scientific computing [25]. Eliminating the first block gives the positive Schur complement \[\label{eq:positive95schur95general} -S = C + GM^{-1}G^{\top} \succ 0 .\tag{14}\] The Schur complement is the central object in block preconditioning for 13 . For standard block diagonal and block triangular preconditioners, the convergence of the resulting fixed-point or Krylov iteration is governed by the conditioning and spectral distribution of the preconditioned Schur complement, \({\widetilde{S}}^{-1}S\), for Schur complement approximation \({\widetilde{S}}\approx S\) [1], [2], [25].

In the setting considered here, the Schur complement approximation is obtained by replacing the exact mass inverse in 14 by a sparse SPD approximation \(\widetilde{M}^{-1}\): \[\label{eq:approx95positive95schur95general} -{\widetilde{S}} \mathrel{\vcenter{:}}= C + G\widetilde{M}^{-1}G^{\top} \succ 0 .\tag{15}\] Below we show that spectral bounds for the preconditioned mass operator \(\widetilde{M}^{-1}M\) directly imply spectral bounds for \({\widetilde{S}}^{-1}S\). Thus, a stronger approximation to \(M^{-1}\) may give a better conditioned preconditioned Schur complement, which in turn improves the convergence of the outer iteration on the saddle-point system. Note that the same question of mass inverse approximation arises in nonsymmetric mixed formulations and saddle-point problems as well; here we assume symmetry to enable the spectral analysis in the following section.

4.1 Spectral mass bounds imply spectral Schur complement bounds↩︎

For approximate mass inverse \(\widetilde{M}^{-1}\), let \(\alpha_C(\widetilde{M}^{-1}) = \lambda_{\min}( S_C^{-1} {\widetilde{S}}_C )\) and \(\beta_C(\widetilde{M}^{-1}) = \lambda_{\max}( S_C^{-1} {\widetilde{S}}_C )\) denote the constants of spectral equivalence as in 5 for \(({\widetilde{S}}_C, S_C)\), so that \(\widehat{\kappa}({\widetilde{S}}_C, S_C) = \widehat{\kappa}(S_C,{\widetilde{S}}_C) = \frac{\beta_C(\widetilde{M}^{-1})}{\alpha_C(\widetilde{M}^{-1})}\). Further define related constants \[\label{eq:schur95a0-c0} \begin{align} \alpha_0(\widetilde{M}^{-1}) := \min_{\boldsymbol{q} \neq \boldsymbol{0}} \frac{\boldsymbol{q}^{\top}G\widetilde{M}^{-1}G^{\top}\boldsymbol{q}}{\boldsymbol{q}^{\top}GM^{-1} G^{\top}\boldsymbol{q}}, \qquad \beta_0(\widetilde{M}^{-1}) := \max_{\boldsymbol{q}\neq \boldsymbol{0}} \frac{\boldsymbol{q}^{\top}G\widetilde{M}^{-1}G^{\top}\boldsymbol{q}}{\boldsymbol{q}^{\top}GM^{-1} G^{\top}\boldsymbol{q}}, \end{align}\tag{16}\] which correspond to \(\alpha_C(\widetilde{M}^{-1})\) and \(\beta_C(\widetilde{M}^{-1})\) when \(C=0\), as well as a spectral equivalence of \((\widetilde{M}, M)\) restricted to the subspace defined by \(\operatorname{range}(G^\top)\). In the following results we show that \(\alpha_C\) and \(\beta_C\) are controlled by \(\alpha_0\) and \(\beta_0\), and that those in turn are controlled by the spectrum of the preconditioned mass matrix.

Note in the case \(C=0\), the bound \(\widehat{\kappa}(S_0,{\widetilde{S}}_{0,\mathrm{bras}}) \le \widehat{\kappa}(M,M_{\rm bras})\) can be strict: the Schur-complement Rayleigh quotients \(\alpha_0, \beta_0\) only sample the preconditioned mass \(M_{\rm bras}^{-1}M\) on \(\operatorname{range}(G^{\top})\). Thus, if either of the extremal modes of \(M_{\rm bras}^{-1}M\) lie outside this range we get the strict inequality \(\widehat{\kappa}(S_0,{\widetilde{S}}_{0,\mathrm{bras}}) < \widehat{\kappa}(M,M_{\rm bras})\). In fact, experiments in 5.1 show that \(\operatorname{range}(G^{\top})\) may very effectively filter out the extremal modes of \(M_{\rm bras}^{-1}M\), where in some cases we observe \(\widehat{\kappa}(S_0,{\widetilde{S}}_{0,\mathrm{bras}}) \ll \widehat{\kappa}(M,M_{\rm bras})\).

5 Numerical results: Block preconditioning↩︎

This section tests Schur-complement-based block preconditioners for saddle-point and nonsymmetric mixed systems, wherein we use either a diagonal or BRAS mass inverse approximation. The approximate Schur complements are inverted only approximately using algebraic multigrid (AMG) inside an outer Krylov method. Note that diagonal and BRAS inverses generally lead to different approximate Schur complement matrices, with the BRAS-based matrices being relatively denser. In the experiments below we report results for small polynomial degrees \(p\), where the AMG methods considered here give reliable Schur-complement preconditioners for both the diagonal and BRAS choices. Developing \(p\)-robust solvers for these approximate Schur complements is a separate issue. We note that if such \(p\)-robust solvers were available, the theory of 4.1, together with the mass-only results of 3, suggests that the gains from BRAS over diagonal preconditioning would increase with \(p\). We first consider a mixed Poisson problem in 5.1 and then a mixed biharmonic problem in 5.2.

5.1 Mixed Poisson in \(H(\operatorname{div}) \times L^2\)↩︎

Consider the Poisson problem
\(-\relax\mathop{\mathrm{grad}}u = f\) on \(\Omega=(0,1)^d\), \(d=2,3\), with homogeneous Dirichlet data for \(u\) on \(\partial\Omega\). We write this equation in first-order form as \[\label{eq:poisson95first95order} \boldsymbol{\sigma} - \mathop{\mathrm{grad}}u = 0 \quad \mathrm{in}\; \Omega, \qquad \operatorname{div}\boldsymbol{\sigma} = -f \quad \mathrm{in}\; \Omega, \qquad u = 0 \quad \mathrm{on}\; \partial \Omega.\tag{17}\] The mixed weak problem is: find \((\boldsymbol{\sigma},u)\in H(\operatorname{div};\Omega)\times L^2(\Omega)\) such that \[\label{eq:mixed95poisson95weak} \begin{align} (\boldsymbol{\sigma},\boldsymbol{\tau})_\Omega + (u,\operatorname{div}\boldsymbol{\tau})_\Omega &= 0 &&\forall \boldsymbol{\tau}\in H(\operatorname{div};\Omega), \\ (\operatorname{div}\boldsymbol{\sigma},v)_\Omega &= -(f,v)_\Omega &&\forall v\in L^2(\Omega). \end{align}\tag{18}\] We use the \(\mathrm{BDM}_p/\mathrm{DG}_{p-1}\) mixed finite element pair, denoted by \(\Sigma_h\subset H(\operatorname{div};\Omega)\) and \(Q_h\subset L^2(\Omega)\), respectively; see [22]. The discrete saddle-point problem is \[\label{eq:mixed95poisson95saddle} \begin{bmatrix} M & B^{\top} \\ B & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\sigma} \\ \boldsymbol{u} \end{bmatrix} = \begin{bmatrix} \boldsymbol{0} \\ - \boldsymbol{f} \end{bmatrix},\tag{19}\] where \(M_{ij} = (\boldsymbol{\phi}_j,\boldsymbol{\phi}_i)_\Omega\), \(B_{ij} = (\operatorname{div}\boldsymbol{\phi}_j,\chi_i)_\Omega\), and \((\boldsymbol{f})_i = ( f, \chi_i )_{\Omega}\). Eliminating \(\boldsymbol{\sigma}\) gives the positive \(u\)-block Schur complement \(-S=BM^{-1}B^{\top} \succ 0\), which is generically dense since \(M^{-1}\) is generically dense.

We consider simplicial meshes, consisting of triangles in two dimensions and tetrahedra in three. The two discrete spaces used here have different sizes, depending on \(p\) and the spatial dimension; ignoring boundary effects, the ratio between the number of \(\Sigma_h\)-DOFs and \(Q_h\)-DOFs is approximately \[\begin{align} d=2:\; \frac{\dim\Sigma_h}{\dim Q_h} &\approx \begin{cases} 3, & p=1,\\ 2.5, & p=2, \end{cases} & \qquad d=3:\; \frac{\dim\Sigma_h}{\dim Q_h} &\approx \begin{cases} 6, & p=1,\\ 4.5, & p=2. \end{cases} \end{align}\] Thus, the \(\boldsymbol{\sigma}\)-block represents a substantial fraction of the block system 19 , especially at low order, which is an important observation for the preconditioning strategy we describe below.

We solve 19 with MINRES [26] using a block-diagonal preconditioner \(\mathcal{P}_{X_{\boldsymbol{\sigma}},X_u}^{-1}\) that approximates the leading mass block and Schur complement as follows:⁶ \[\label{eq:mixed95poisson95split95minres95prec} \mathcal{P}_{X_{\boldsymbol{\sigma}},X_u}^{-1} := \begin{bmatrix} X_{\boldsymbol{\sigma}} & 0 \\ 0 & \mathcal{V}_{X_u} \end{bmatrix} \approx \begin{bmatrix} M^{-1} & 0 \\ 0 & -S^{-1} \end{bmatrix} =: \mathcal{P}_{\rm ideal}^{-1}.\tag{20}\] Here, \(X_{\boldsymbol{\sigma}},X_u\in\{M_{\rm diag}^{-1},M_{\rm bras}^{-1}\}\), and \(\mathcal{V}_{X_u}\) denotes one AMG V-cycle for the approximate Schur complement \(-{\widetilde{S}}_{X_u}=BX_uB^{\top}\). Note that there is an additional indirect effect on convergence related to the efficacy of AMG on the approximate Schur complement formed with \(M_{\rm diag}^{-1}\) vs. \(M_{\rm bras}^{-1}\), but we do not consider this comparison in detail.

4 presents spectral condition numbers of the associated preconditioned operators on some fixed, small meshes. The quantity \(\widehat{\kappa}( M, X_{\boldsymbol{\sigma}}^{-1})\) measures the mass approximation that appears in the \(\boldsymbol{\sigma}\)-block, while \(\widehat{\kappa}(S, {\widetilde{S}}_{X_u})\) measures the Schur complement approximation assuming exact inversion of \({\widetilde{S}}_{X_u}\), rather than one AMG V-cycle. The spectral condition number \(\widehat{\kappa}( S, {\widetilde{S}}_{X_u})\) is covered by the \(C=0\) special case of 3; rewriting that result to the current notation, we have: \[\label{eq:schur95condition95number95zero95C95bound95copy} \widehat{\kappa}( S,{\widetilde{S}}_{X_u}) = \frac{ \displaystyle \max_{\boldsymbol{q}\neq \boldsymbol{0}} \frac{\boldsymbol{q}^{\top}B X_u B^{\top}\boldsymbol{q}}{\boldsymbol{q}^{\top}B M^{-1} B^{\top}\boldsymbol{q}} }{ \displaystyle \min_{\boldsymbol{q} \neq \boldsymbol{0}} \frac{\boldsymbol{q}^{\top}B X_u B^{\top}\boldsymbol{q}}{\boldsymbol{q}^{\top}B M^{-1} B^{\top}\boldsymbol{q}} } \leq \frac{ \displaystyle \max_{\boldsymbol{g}\neq \boldsymbol{0}} \frac{\boldsymbol{g}^{\top} X_u\boldsymbol{g}}{\boldsymbol{g}^{\top} M^{-1}\boldsymbol{g}} }{ \displaystyle \min_{\boldsymbol{g} \neq \boldsymbol{0}} \frac{\boldsymbol{g}^{\top} X_u \boldsymbol{g}}{\boldsymbol{g}^{\top} M^{-1} \boldsymbol{g}} } = \widehat{\kappa}( M, X_u^{-1}).\tag{21}\] That is, the spectral condition number of the preconditioned Schur complement is bounded by that of the preconditioned mass, and more sharply depends only on the action of the preconditioned mass \(X_u M\) on \(\operatorname{range}(B^{\top})\). This restricted nature is clearly visible in 4, where spectral condition numbers of preconditioned Schur complements are typically much smaller than the corresponding preconditioned mass spectral condition numbers \(\widehat{\kappa}( S,{\widetilde{S}}_{X_u}) \ll \widehat{\kappa}( M, X_u^{-1})\), particularly in three dimensions and at higher order. While BRAS gives a uniformly better approximation to \(M^{-1}\) than diagonal scaling in the full mass comparison, i.e., \(\widehat{\kappa}( M, M_{\rm bras}) \ll \widehat{\kappa}( M, M_{\rm diag})\), its improvement in the Schur complement approximation is more nuanced. In the precision displayed, for \(p=1\) the diagonal/BRAS values are \(2.9/2.9\) in two dimensions and \(3.0/3.0\) in three dimensions. Thus, at \(p=1\), the BRAS approximation improves the \(\boldsymbol{\sigma}\)-block but does not improve the Schur complement approximation in the \(u\) block, despite the fact that the BRAS-based Schur complement is more expensive to build and to use in AMG. This motivates the “hybrid” preconditioning choice of \(X_{\boldsymbol{\sigma}}=M_{\rm bras}^{-1}\), \(X_u=M_{\rm diag}^{-1}\). Note that for \(p\ge2\) BRAS uniformly improves the condition number of the preconditioned Schur-complement.

3pt

Next we consider the numerical solution of 19 with MINRES and the preconditioner \(\mathcal{P}_{X_{\boldsymbol{\sigma}},X_u}^{-1}\) from 20 . We use the MINRES implementation provided by SciPy [27], with tolerance \(\texttt{rtol=1e-10} \| \boldsymbol{f} \|\), and zero initial iterate, and the source function \(f\) chosen randomlu. We consider three preconditioning variations, given by “diag”, corresponding to \(X_{\boldsymbol{\sigma}}=X_u=M_{\rm diag}^{-1}\), “BRAS”, corresponding to \(X_{\boldsymbol{\sigma}}=X_u=M_{\rm bras}^{-1}\), and “hybrid” corresponding to \(X_{\boldsymbol{\sigma}}=M_{\rm bras}^{-1}\), \(X_u=M_{\rm diag}^{-1}\). The AMG hierarchy for \({\widetilde{S}}_{X_u}\) is constructed using classical Ruge–Stüben AMG from the PyAMG library [28]. We use classical strength of connection, RS coarsening with a second pass, classical interpolation, and symmetric Gauss–Seidel pre- and post-smoothing. The strength-of-connection threshold is selected from the fixed candidate set \(\{0.15,0.25\}\). For each fixed choice of dimension, \(p\), and Schur-complement approximation \(X_u\), we choose the strength threshold to be the one that yields the fastest convergence of standalone AMG on \({\widetilde{S}}_{X_u}\) across a variety of representative meshes.

4 compares the three preconditioning approaches: diag, hybrid, and BRAS. The left column reports \(t_{\rm diag}/t_{\rm hybrid}\), the middle column reports \(t_{\rm diag}/t_{\rm bras}\), and the right column reports the corresponding MINRES iteration counts. Thus values larger than one in the first two columns indicate faster solves than diag. For \(p=1\), the hybrid preconditioner gives the strongest performance, reducing MINRES iteration counts and solve times relative to diag. BRAS reduces MINRES iterations over diag, but not solve time, which is consistent with 4, because there BRAS improves the \(\boldsymbol{\sigma}\)-block mass approximation, but gives essentially the same conditioning for the Schur complement as diagonal. For \(p=2\), the hybrid preconditioner still improves substantially over diag, but BRAS gives the lowest iteration counts and the best solve times on the largest cases, with 1.3–2\(\times\) speedup over the diagonal reference. This is consistent with the spectral condition number results from 4, where at \(p=2\) BRAS improves both the \(\boldsymbol{\sigma}\)-block mass approximation and the \(u\)-block Schur-complement approximation.

5.2 Mixed biharmonic in \(H^1 \times H^1\)↩︎

We next consider a two-dimensional, scalar mixed formulation arising from an anisotropic fourth-order tensor-divergence operator, which can be seen as a scalar version of the vector-valued electron viscosity term within Ohm’s law in extended magnetohydrodynamics [29]. Discretizing the PDE problem in time with a diagonally implicit Runge-Kutta method results in a semi-discrete system of the form \[\begin{align} \label{eq:biharmonic95semi} \chi + \delta t \, \nabla^\perp \cdot \boldsymbol{E} = f, \qquad \boldsymbol{E} = -\nabla \cdot \big(\Pi_e(\boldsymbol{j}_e)\big), \qquad \boldsymbol{j}_e = \nabla^\perp \chi, \end{align}\tag{22}\] with constant \({\delta t > 0}\) the time-step size (possibly up to some proportionality constant). Here, \(\chi\) is an unknown scalar field (corresponding to a magnetic field component), \(\boldsymbol{E}\) is the electric field, \(\boldsymbol{j}_e\) is the current density, and \(f\) is a known function. The electric field is equipped with perfect conducting wall boundary conditions \(\boldsymbol{E} \cdot \boldsymbol{n}^\perp = 0\), and, for simplicity, we assume a zero traction condition \({{\Pi}_e(\boldsymbol{j}_e) \cdot \boldsymbol{n} = \boldsymbol{0}}\), for boundary normal vector \(\boldsymbol{n}\). Finally, we set the anisotropic electron stress tensor to \[\Pi_e(\boldsymbol{j}_e) = -A \big(A \, : \, \nabla \boldsymbol{j}_e \big), \qquad A := \big(\boldsymbol{b} \otimes \boldsymbol{b} - \tfrac{1}{3}I\big)\] with \(\boldsymbol{b}\) a prescribed, unit-length anisotropy direction field.⁷ One can show that 22 simplifies to the fourth-order equation \[\chi+\delta t \, D^\ast_A D_A \chi=f, \qquad \mathrm{where} \: D_A u = A \, :\, \nabla(\nabla^\perp u).\] In our tests, we set \(\delta t = h\), for spatial mesh size \(h\), in view of balancing spatial and temporal discretization errors; explicit time integration would require a CFL condition of the form \(\delta t \lesssim h^4\), so the linear systems considered here are numerically stiff.

We now pose this problem in mixed form by introducing the auxiliary variable \[\zeta = \sqrt{\delta t} A \, : \, \nabla \boldsymbol{j}_e = \sqrt{\delta t} A \, : \, \nabla (\nabla^\perp \chi),\] and the associated second-order coupling \[c(v,u) = \sqrt{\delta t} \int_\Omega \nabla^\perp v\cdot \nabla\cdot\bigl(u\, A \bigr) \,\,\mathrm{d}x .\] Given \(V_h \subset H^{1}(\Omega)\), the mixed weak problem is to find \((\chi_h,\zeta_h)\in V_h\times V_h\) such that \[(\eta_h,\chi_h)_\Omega-c(\eta_h,\zeta_h) = (\eta_h,f_h)_\Omega, \qquad (\gamma_h,\zeta_h)_\Omega+c(\chi_h,\gamma_h) = 0\] for all \((\eta_h,\gamma_h)\in V_h\times V_h\). Here \(V_h \subset H^1(\Omega)\) is an unrestricted finite-element space, and the aforementioned homogeneous boundary conditions are incorporated naturally in \(c\) through integration by parts. With \(M_{ij}=(\phi_j,\phi_i)_\Omega\) and \(C_{ij}=c(\phi_i,\phi_j)\), this gives \[\begin{align} \label{eq:biharmonic952x2} \begin{bmatrix} M & C^{\top} \\ -C & M \end{bmatrix} \begin{bmatrix} \boldsymbol{\zeta} \\ \boldsymbol{\chi} \end{bmatrix} = \begin{bmatrix} \boldsymbol{0} \\ \boldsymbol{f} \end{bmatrix}. \end{align}\tag{23}\] We solve this system with FGMRES [30] using a lower triangular preconditioner, \[\begin{align} \mathcal{P}_X = \begin{bmatrix} M & 0 \\ -C & {\widetilde{S}}_X \end{bmatrix}, \qquad {\widetilde{S}}_X := M+ C X C^{\top} \approx M+ C M^{-1}C^{\top} =: S. \end{align}\] For the mass inverse approximations, we consider \(X \in \{M_{\rm diag}^{-1}, \, M_{\rm bras}^{-1} \}\), and we choose \(\boldsymbol{f}\) to be a random vector. In terms of FGMRES, we initialize with the zero vector, iterate until \(\|\boldsymbol{r}_k\| \leq 10^{-10} \| \boldsymbol{r}_0 \|\), and we restart every 30 iterations.

We apply \({\cal P}_X^{-1}\) inexactly as follows. The mass inverse \(M^{-1}\) in the \((1,1)\) block is approximated by CG preconditioned by \(X\), iterated until \(\texttt{rtol=1e-6}\) or the number of iterations exceeds 40.⁸ To approximate \({\widetilde{S}}_X^{-1}\) we use CG preconditioned by one AMG V-cycle, and iterate until \(\texttt{rtol=1e-3}\) or the number of iterations exceeds 20. Specifically, we use the least-squares algebraic-multigrid domain-decomposition (LS-AMG-DD) solver [31], which we have shown to be robust on a different fourth-order problem discretized directly in \(H^2\) [32]. The LS-AMG-DD parameters used are as follows: 1 pass of aggregation for coarsening, pre- and post smoothing given by overlapping multiplicative Schwarz and its adjoint, respectively, the local spectral cutoff is \(\tau_{\rm cut} = 1.5 \max_j {\rm mult}_{\omega}(j)\), the maximum number of levels in the hierarchy is capped at five, and the coarsening is halted if a coarse-level matrix density exceeds 25%.

We remark that to apply LS-AMG-DD to a matrix \({\widetilde{S}}\) requires a Gram factor \(G_{{\widetilde{S}}}\), such that \({\widetilde{S}} = G_{{\widetilde{S}}}^\top G_{{\widetilde{S}}}\) is a Gram matrix. Given Gram factors \(G_M\) and \(G_X\) for the mass \(M\) and its approximate inverse \(X \approx M^{-1}\), such that \(M = G_M^\top G_M\) and \(X = G_X^\top G_X\), we write the approximate Schur complement as a Gram matrix as \[\begin{align} {\widetilde{S}}_X = G_{{\widetilde{S}}}^{\top}G_{{\widetilde{S}}}, \qquad G_{{\widetilde{S}}} = \begin{bmatrix} G_M \\ G_X C^{\top} \end{bmatrix}. \end{align}\] Recalling that the conforming mass matrix can be written as \(M = P^\top M_{\rm br} P\), we compute its Gram factor as \(G_M = M_{\rm br}^{1/2} P\), using the approach outlined in [32]. The Gram factors for the approximate mass inverses are computed as follows: For \(X=M_{\rm diag}^{-1}\), we take \(G_X=M_{\rm diag}^{-1/2}\); for \(X=M_{\rm bras}^{-1}\), we take \(G_X = M_{\rm br}^{-1/2}P D^{-1}\), according to 3 .

Numerical results are shown in 5 for Lagrange basis functions of degrees \(p = 1,2,3\) using GLL nodes. The domain \(\Omega = (0, 1)^2\) is meshed with quadrilateral elements rather than triangular elements as in all other two-dimensional tests in the paper. The direction vector is \(\boldsymbol{b} = ( \cos \theta, \sin \theta )\) with \(\theta = \pi/6 = 30^\circ\). The first column shows the solve time of the diagonally-preconditioned approach relative to the BRAS-preconditioned approach; in all cases, BRAS yields a faster solve time than diagonal by a factor of \(\sim 2\)–\(3 \times\). Considering the middle column, the total number of FGMRES iterations is significantly decreased when using the BRAS strategy relative to the diagonal strategy, and the iteration counts are more tightly clustered as a function of \(p\) for the BRAS approach. Finally, considering the last column, the total number of LS-AMG-DD iterations is significantly smaller for the BRAS approach, and also has much improved scaling with respect to the number of DOFs in the problem.

6 Conclusions↩︎

We introduce BRAS, a sparse algebraic approximate inverse for conforming finite-element mass matrices. BRAS is obtained by applying exact broken element mass inverses and averaging back to the conforming space. The preconditioner is local, SPD, basis independent at the algebraic level, and has the same element-adjacency sparsity as the conforming mass matrix. We characterize the preconditioned mass spectrum using a broken-space projector and show how mass-inverse approximation bounds pass to approximate Schur complements. The numerical results show that BRAS substantially reduces mass-matrix spectral condition numbers and Krylov iteration counts relative to diagonal scaling across \(H^1\), \(H(\operatorname{curl})\), and \(H(\operatorname{div})\) spaces. For both explicit mass matrix solves and mixed block preconditioning problems, BRAS formulations lead to consistent solve-time improvements over diagonal preconditioning in essentially all of the tested cases, typically on the order of 1.5–3\(\times\).

Acknowledgments↩︎

The authors used OpenAI’s ChatGPT during the preparation of this manuscript, including for exploratory mathematical discussion, coding, and drafting of the exposition. All mathematical claims, proofs, computations, citations, and final wording were reviewed, verified, and edited by the authors, who take full responsibility for the manuscript. Discussions with Steven Walton on aspects of mass preconditioning are gratefully acknowledged.

References↩︎