Algebraic Hybridization and Static Condensation with Application to Scalable H(div) Preconditioning
Veselin A. Dobrev Tzanio V. Kolev Chak S. Lee Vladimir Z. Tomov Panayot S. Vassilevski
math.NA 65F10, 65N20, 65N30









Algebraic Hybridization and Static Condensation with Application to
Scalable H(÷) Preconditioning

  This work was performed under the auspices of the U.S. Department of Energy by
  Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344, and
  sponsored by the U.S. Department of Energy, Office of Science, Office of
  Advanced Scientific Computing Research, Applied Mathematics
  program. LLNL-JRNL-732140.

    
V. Dobrev
Center for Applied Scientific Computing,
  Lawrence Livermore National Laboratory,
  P.O. Box 808, L-561,
  Livermore, CA 94551, U.S.A (dobrev1@llnl.gov,
  tzanio@llnl.gov, cslee@llnl.gov, tomov2@llnl.gov, panayot@llnl.gov).
T. Kolev
[2]
C. S. Lee
[2]
V. Tomov
[2]
P. S. Vassilevski
[2]

    December 30, 2023
=========================================================================================================================================================================================================================================================================================================================================================================================================================================




We propose an unified algebraic approach for static condensation and
hybridization, two popular techniques in finite element discretizations.  The
algebraic approach is supported by the construction of scalable solvers for
problems involving H(÷)-spaces discretized by conforming (Raviart-Thomas)
elements of arbitrary order. We illustrate through numerical experiments the
relative performance of the two (in some sense dual) techniques in comparison
with a state-of-the-art parallel solver, ADS <cit.>, available in
<cit.> and <cit.>.  The superior performance of the hybridization
technique is clearly demonstrated with increased benefit for higher order
elements.



static condensation, hybridization, algebraic multigrid, mixed finite elements,
ADS, H(÷) solvers, radiation-diffusion transport.



65F10, 65N20, 65N30


myheadings plain
DOBREV, KOLEV, LEE, TOMOV, and VASSILEVSKI
Algebraic hybridization and static condensation




§ INTRODUCTION


This paper proposes a unified algebraic approach of constructing reduced
problems by two popular techniques in the finite element method: static
condensation and hybridization. It also discusses approaches for constructing
scalable preconditioners for these reduced problems. The present work is
motivated by the solution of problems arising in variety of finite element
simulations and of particular interest is the case when the problems involve the
function space H(÷) (L_2-vector functions that have their weak divergence
in L_2).  Such problems naturally arise in the mixed finite element method for
Darcy equation <cit.>, Brinkman equation <cit.>, as
well as in the FOSLS (first order system least-squares) discretization approach,
<cit.>.  Also, certain formulations of radiation-diffusion transport
<cit.>, lead to problems involving the space H(÷). This is in
fact one of the practical examples from a more realistic simulation that we
consider in our numerical tests.

Highly scalable solvers for H(÷) problems have been designed in the past.
The most successful ones, of ADS-type <cit.>, were based on the
regular decomposition theory developed by Hiptmair and Xu <cit.>.
Although of optimal complexity, the current state-of-the-art ADS algorithms are
quite involved and require also going through solvers for H(), such as
AMS, <cit.>. Both solvers, ADS and AMS, require some
additional user input (e.g., discrete gradient), which makes them not completely
algebraic.  The results of the present paper can be viewed as a further
step towards the design of fully algebraic solvers for problems involving
H(÷).

To design faster, yet still scalable, alternative to ADS, we employ a
traditional finite element technique commonly used in mixed finite elements
referred to as hybridization. This is perhaps the first technique proposed
to solve the saddle-point problems arising in the mixed methods, since it leads
to symmetric and positive definite (SPD) reduced problems. Its early form
appeared in a paper by Fraeijis de Veubekein in 1965 <cit.> as an
efficient solution strategy for elasticity problems. In <cit.>,
Arnold and Brezzi studied the method from a theoretical perspective, and
obtained error estimates for all the unknowns involved; see also
<cit.>. The technique was further developed by Cockburn and Gopalakrishnan
to allow approximations by polynomials of different degree in different
subregions; see for example <cit.> and <cit.> for applications to
mixed finite element approximations of second order scalar elliptic problems and
Stokes flow.

Another very popular technique in finite elements is static condensation.
It is a reduction technique intended to reduce the size of a linear system of
finite element equations by eliminating at element or subdomain level internal
degrees of freedom (dofs). It was also first introduced in the structural
analysis literature in 1965, cf., <cit.>, and has
been widely used since then. Static condensation neglects the dynamic effect in
the reduction, hence the name. On element level, one can eliminate locally
degrees of freedom that are not shared by neighboring elements (a typical case
for high order elements). Thus, the reduced, Schur complement, problem contains
only shared (interface) dofs. If the elimination is performed on a subdomain
level, the resulting reduced system coincides with the one in the balancing
domain decomposition by constraints (BDDC) algorithm <cit.>. Recently, BDDC
preconditioners with deluxe scaling and adaptive selection of primal constraints
for H(÷) problems have been developed, see for instance <cit.>. The approach in the present paper, on the other
hand, solves the reduced system by algebraic multigrid (AMG).

For H(÷)-conforming elements (i.e., elements that have only shared dofs
through common faces), both static condensation and hybridization lead to
reduced problems of the same size (equal to the number of shared interface
dofs). If one applies these two approaches not on element level, but for
subdomains (union of elements), substantial reduction in size can be
achieved. Thus, these two approaches are of great practical interest and the
goal of the present paper is to study them in a common framework, emphasize
their similarities and differences, and most importantly, design new, or modify
existing solution techniques that are efficient and scalable, achieving
substantial savings compared to the state-of-the-art solvers for the original
problems. For problems obtained by static condensation, there has been success
in modifying the state-of-the-art solvers (AMS and ADS) to be directly
applicable to the reduced problem, cf., <cit.>. In this paper, we focus on the design of scalable solvers for
problems involving H(÷) in a general algebraic setting, extending the
preliminary results in <cit.>.

The remainder of the paper is organized as follows. In Section
<ref>, we describe the algebraic hybridization approach and
show how it can serve as an alternative to traditional finite element assembly.
Next, we deduce its relation to static condensation and the mixed hybridized
finite element method in Section <ref>.
Preconditioning strategies for the reduced systems obtained from hybridization
and static condensation are discussed in Section <ref>.
Then we describe the implementation details of the methods and present some
numerical examples in Section <ref> and Section
<ref> respectively. Finally, we conclude the paper in
Section <ref>.







§ TWO PERSPECTIVES OF FINITE ELEMENT ASSEMBLY



In this section we describe the traditional process of finite element assembly,
and introduce algebraic hybridization as an alternative process to obtain a
global linear system describing the same finite element discretization.




 §.§ Finite element assembly components



To fix notation, consider a bilinear form a(·,·) and a linear form
f(.) over a Hilbert space ℋ defining the Galerkin weak variational
formulation: Find u ∈ℋ, such that for all v ∈ℋ, we have

    a(u,v) = f(v).


In the finite element method, we introduce a triangulation =_h of the
given computational domain and decompose the forms into sums of their local
element versions, i.e a(·,·) = ∑_τ∈
a_τ(·,·) and f(.) = ∑_τ∈ f_τ(.).  We
further replace ℋ with a finite dimension space induced by _h and on each
element τ∈ we compute a local element stiffness matrix A_τ
and a load vector _τ by using the element basis functions in the
local bilinear and linear forms.  We are assuming that A_τ (which
correspond to a_τ(·,·)) are SPD. The latter holds for example if
the bilinear form a(u,v) contains a mass term, ∫β uv dx, with a
positive coefficient β=β(x) > 0.




 §.§ The traditional finite element assembly procedure



A central concept in the finite element method is that of assembly – the
accumulation of local stiffness matrices A_τ and load vectors _τ
into a global matrix A and a global vector  which form the (global)
linear system

    A  = ,

for the unknown vector of finite element degrees of freedom.  The n × n
matrix A is SPD and generally sparse, i.e., the number of its nonzero entries,
nnz(A), is Ø(n).

Let A = blockdiag(A_τ) be the n̂×n̂
block-diagonal matrix with all element stiffness matrices on its main diagonal,
and  be similarly the n̂× 1 vector containing all
element load vectors.  The finite element method produces a global-to-local
mapping, P, which maps global dofs into local ones for each element.  That is,
for a n ×n̂ matrix P, the traditional assembly process can be
written as the triple matrix product that relates A and A, and a
matrix-vector product that builds the global right-hand side,

    A = P^T A P ,
      and   = P^T  .

In the simplest case the action = P simply means copying the
values of the global dof vector  to the corresponding local values,
independently in each element. More complicated P can account for
non-conforming mesh refinement as we describe in some detail in a following
section.  Either way, (<ref>) is a formalized expression for the
traditional process of finite element assembly.




 §.§ Algebraic hybridization



It is classical in the finite element literature to introduce (<ref>) as the
global linear system that corresponds to the finite element discretization
specified by A,  and P.  As we describe next,
there is an alternative global linear system corresponding to the same solution,
which in certain cases is advantageous from the perspective of linear solvers.
To introduce this alternative, let C be another matrix such that

    (C) = (P),

i.e. a vector z satisfies C z = 0, if and only if there is a vector y,
such that z = P y.  In particular, we have the matrix equality C P = 0.

We assume that P is full column rank and that C is full row rank. Similarly
to P, the matrix C could also be provided algebraically by the finite
element space, e.g. the columns of C^T can be constructed by orthogonally
completing the basis given by the columns of P. (The matrix C is not unique
and there are other ways to construct it as illustrated later in the paper.)

The idea of algebraic hybridization is that instead of adding the different
element contributions of A, , we use C to enforce
equality constraints of the decoupled vector  through a new
Lagrange multiplier variable .  This process is known as hybridization in the context of mixed finite element method <cit.>, but
in this section we consider it on a purely algebraic level, independent of the
particular continuous problem or finite element discretization approach.  In
fact, the C-based reformulation can be described in more general terms than
classical assembly, which we do in the following main linear algebra result.


Let A, P and C be matrices of size n̂×n̂, n̂× n and m ×n̂ (m=n̂-n) respectively, such that A
is SPD, P is full column rank, C is full row rank, and (C) = (P).
Also, let  be a column vector of size n̂. Then the problem

    P^TA P  = P^T ,

is equivalent to the following, hybridized, saddle-point system

    [   A C^T;   C   0 ][ ;  ] =
    [  ; 0 ].

Specifically, the saddle-point system (<ref>) is
uniquely solvable and the solution  of
(<ref>) and  (the solution of the original
system), are related as

    = P .

The solution  is computable via standard (block-)Gaussian elimination:
we first compute the Schur complement system for the Lagrange multiplier

    H ≡ C A^-1 C^T  = C A^-1,

then, by back substitution, we have

    = A^-1(-C^T).



The saddle-point system is invertible due to our assumptions; namely,
A is SPD (hence invertible), and C^T is full column
rank, hence the (negative)
Schur complement H = C A^-1 C^T is also
SPD and invertible.
From A + C^T  = and
CP =0, we obtain

    P^T A = P^T  =: .

Now, since C  = 0 (the second equation of
(<ref>)), and from the assumption (C) = (P),
it follows that there is an  such that = P.
Then from = P^T A = P^T A
P, we see that this  is unique; namely,  is the unique
solution of the original problem (<ref>).


We can summarize the results of the theorem and the previous section by stating
that P^TA P and C A^-1 C^T are equally valid global
matrices for the same problem.  In the case of finite elements, the first one
corresponds to assembly, while the second one is the hybridized system for
Lagrange multipliers.

Note that the condition that P has full column rank is equivalent to the
existence of a left inverse of P: there is an n ×n̂ matrix R
such that R P = I, e.g. R=(P^T P)^-1 P^T.  We can use this restriction
matrix R to compute the solution  of (<ref>)
from the solution  of (<ref>) as = R. Furthermore, if instead of  we have access only to
the assembled right-hand side vector = P^T, we can use
= R^T in (<ref>) instead of
 since P^T = P^T =.

Note also that the first set of equations in (<ref>),

    A =  - C^T.

can be interpreted as a local version of the problem (in a weak form), where the
term C^T plays the role of a dual vector for Neumann-type boundary
conditions.
Indeed, if we integrate by parts locally, we get

    A =  - F

where F is the Neumann data for  tested against
the test functions in each element. Since the exact solution satisfies
(<ref>), we have P^T F = 0, i.e.
F = C^T for some  by (P^T)=(C^T).
One can then think of hybridization as prescribing Neumann data for local
problems from one set of  values and then imposing the fact that the
local solutions (fluxes) should match on their shared dofs to derive equations
for ; for more details, see <cit.>.

The result of Theorem <ref> can be extended to
more general settings that relax the conditions that A and P are
full column rank. This is illustrated in the next theorem.


Assume that A is a square matrix (without assuming that it is SPD or
even invertible) and that (C)=(P) (without assuming that P has full
column rank and C has full row rank). Then equations
(<ref>) and (<ref>) are equivalent in
the following sense


  * If  solves (<ref>) then there exists 
such that (=P,) solves (<ref>).

  * If (,) solves (<ref>) then there exists
 that solves (<ref>) such that
P =.

If, in addition, we assume that A is invertible then equations
(<ref>) and (<ref>) are also
equivalent in the following sense


  * If (,) solves (<ref>) then 
solves (<ref>).

  * If  solves (<ref>) then
(= A^-1(- C^T),) solves
(<ref>).



Assume that  solves (<ref>) and set
=P. We need to show that there exists  such that
(,) solves (<ref>). The second equation
in (<ref>) follows from C = C P  = 0. The
defining condition for  is

    A+ C^T  =  ,
         or    
    C^T  = - A=:  .

By (<ref>), we have that
P^T = 0, i.e. ∈(P^T). Since we have

    (P^T) = (P)^⊥⊆(C)^⊥ = (C^T)  ,

the existence of the required Lagrange multiplier  is proven.

Assume now that (,) solves (<ref>).
On one hand, the second equation in (<ref>) is
C =0, i.e. ∈(C) and since
(C)⊆(P), there exists  such that P=.
On the other hand, the first equation in (<ref>), implies

    P^T (A+ C^T  ) = P^T  ,
          or    
    P^T A= P^T  ,

showing that  solves (<ref>).

The proof of the second part of the theorem is straightforward.


Note that the uniqueness of  is not guaranteed in the settings of Theorem
<ref>, even if we assume that P has
full column rank and A is invertible. Similarly, the uniqueness of
 does not follow from the additional assumption that C has full row
rank and A is invertible.




 §.§ A main example of matrices P and C


Below we discuss an application of the main theorem to a slightly more general
case of finite element discretization.  This is the main example we have in mind
in the application to H(÷) bilinear form considered in the present paper.


In a finite element setting, we may think of the following P originating from
a generally non-conforming finite element space.

Assume, that we have elements (or more generally subdomains that are unions of
finite elements).  We also have interfaces between these substructures (original
elements or subdomains). There are degrees of freedom interior, “i”, to the
substructures. We also have dofs on the interfaces. In this setting we assume
that on the interfaces we have dofs that are master, “m”, and slave, “s”.

Introduce the matrices P and C with the described blocks of i, m, and
s-dofs,

    P =
    [  I_i    0;    0  I_m;    0 W_sm ]    and    
    C = [     0 -W_sm   I_s ].

The 3-by-2 block-form of P above, corresponds to the interior and master dofs
(for its columns) and interior, master and slave dofs (for its rows). The block
W represents the mapping that determines the slave dofs in terms of the master
ones. Finally, the matrix A, is the block-diagonal matrix, with
blocks assembled for the individual substructures (or simply the element
matrices).

As before, we assume that A is invertible, which is the case if we
have a mass term added to a semidefinite term, like the H(÷) bilinear form
we consider in this paper.



Consider the case when the solution space is H(÷)-conforming, and introduce
a multiplier finite element space, ℳ⊂ L^2(ℐ) on the
interface ℐ between elements (or substructures), see Figure
<ref> for illustration.

If we denote the broken solution space by 𝒱, then the
matrix C can be derived from the bilinear form:

    c(,λ) =
      ∫_ℐ·λ ,    ∈𝒱, λ∈ℳ .

Above, · stands for the jump of the normal component
of  across the interface ℐ.

Similarly to Example <ref>, for the same P, the matrix C
now has the M_s-weighted form

    C = [         0 -M_s W_sm       M_s ] ,

where M_s is a mass matrix between the multiplier space ℳ and the
normal trace space of 𝒱.





§ RELATION TO STATIC CONDENSATION AND HYBRIDIZED MIXED FINITE ELEMENTS



In this section we show that the hybridized matrix H, (<ref>), is closely related to the matrices obtained with two other
classical finite element approaches: that of algebraic static condensation
and that of traditional hybridization used in the mixed finite element method
for second order elliptic equations.




 §.§ Static condensation


In the setting of Example <ref>, the matrix A can
be partitioned into a two-by-two block from with its first block corresponding
to the interior dofs (with respect to the substructures). The second block
corresponds to the master and slave dofs combined, i.e., we have

    A =  [
    [ A_ii A_ib; A_bi A_bb ] ].

We notice that A_ii is block-diagonal, since the interior dofs in
one substructure do not couple with interior dofs in any other substructure. The
Schur complement

    S = A_bb -
        A_biA^-1_iiA_ib,

is SPD, hence invertible. It acts on vectors corresponding to the interface dofs
“b” (combined slave and master ones).
Now, consider the Schur complement H = C A^-1 C^T
from (<ref>) for the hybridized system
(<ref>).
Since C = [0,  -W,  I], and

    A^-1 =  [
    [    *    *;      S^-1 ] ],

we easily see that

    H = C A^-1C^T = [-W, I] S^-1 [-W, I]^T,

while the static condensation matrix, i.e., the reduced matrix for the
master dofs is

    S =  [
    [ I; W ] ]^TS [
    [ I; W ] ].

In the case of finite element matrices, when the coupling across elements (or
more generally, across substructures) occurs only via interface dofs (no vertex
and edge dofs), one can decouple each interface dof into exactly two copies, and
let W = -I.  Examples of such elements are the Raviart–Thomas
(H(÷)-conforming elements); and also, the H^1-non-conforming
Crouzeix-Raviart elements.

In this setting A is block-diagonal, and also S is
block-diagonal (substructure-by-substructure).  Thus, to form the hybridized
Schur complement H, we need to assemble the inverses of the local substructure
Schur complements, whereas to form the static condensation matrix S, we need
to assemble the local substructure Schur complements.

In other words, we may view the hybridization approach as a sort of a “dual”
technique to the static condensation approach.  This has important consequences
for the construction of solvers for the respective problems. For example, we can
expect that solvers for S will have the same nature as solvers for the
original matrix A, whereas solvers for H, will need preconditioners that are
effective in the dual to the original trace space.  This is the case for the
H(÷) problem, as observed in <cit.>.






 §.§ Hybridized mixed finite element method


We next consider traditional hybridization for mixed finite elements,
cf. <cit.>, applied in the settings of the H(÷) problem:
find ∈ H(÷;Ω) such that

    a_÷(,) = (α ÷, ÷) + (β, ) = (, )
            ∀ ∈ H(÷;Ω),

where the coefficients α and β are positive and in general
heterogeneous.  In this case,  has the following special form

    = ^T ^-1 + .

By introducing an additional unknown q=÷, a mixed formulation of
(<ref>) can be written as: find (,q)∈ H(÷;Ω)×
L^2(Ω) such that

    (β,)+(q,÷)     = (,)    ∀ ∈ H(÷;Ω),
    
    (÷,p) - (α^-1q,p)     = 0        ∀ p∈ L^2(Ω).

Notice that (<ref>) coincides with the mixed formulation of the
scalar second order elliptic problem

    -÷(β^-1∇ q) + α^-1q = g

by letting =β^-1∇ q. Hence, (<ref>) can be
discretized by the hybridized mixed finite element methods designed for
(<ref>), c.f. <cit.>, resulting in a discrete problem of the form

    [      ^T C^T;       -   0;   C   0   0 ][  ; q;   ]
    =
    [  ; 0; 0 ].

The corresponding reduced problem for the Lagrange multipliers reads

    H=
    [ C 0 ][    ^T;     -;    ]^-1[ C^T;   0 ] =
    [ C 0 ][    ^T;     -;    ]^-1[  ; 0 ],

which is identical to (<ref>) because of
(<ref>) and the fact that the (1,1) entry in the above
inverse is the inverse of the Schur complement + ^T ^-1.  In
other words, H = H, so the algebraic hybridization approach is
identical when applied to the mixed or second order version of the problem.




§ SCALABLE ALGEBRAIC H(÷) PRECONDITIONING
 

In this section we propose several scalable algebraic multigrid approaches for
finite element discretizations of the H(÷) form a_÷(·,·)
based on the methods discussed in the previous sections.  Due to their
scalability on challenging practical problems, we base our preconditioners on
the parallel algebraic multigrid approaches available in the hypre
library, <cit.>.  Specifically we take advantage of the BoomerAMG
<cit.>, AMS <cit.> and ADS <cit.> methods in
hypre, targeting H^1, H() and H(÷) discretizations
respectively.

We first consider the static condensation Schur complement S, (<ref>), and discuss approaches for its preconditioning.  Since
S acts in general on traces of functions from a given space (H(÷),
H(), or H^1), a suitable preconditioner for S will be some
approximation to the Schur complement of an (optimal) preconditioner for the
original matrix A.  Of course, for efficiency, we would like to build a
preconditioner directly for the (sparse) matrix S and not for the original
matrix A.  Such preconditioners are discussed in <cit.>
where it is shown, both theoretically and practically, that BoomerAMG and AMS
work well on Schur complements of H^1 and H() discretizations
respectively. Similar analysis for the H(÷) case was recently performed in
<cit.>.  In all three cases, the summary of the stable decomposition
theory is that if BoomerAMG/AMS/ADS works for a particular problem, it will also
work, purely algebraically, for a Schur complement of the problem.  Thus, the
H(÷) Schur complement S considered in this paper can be efficiently
preconditioned simply with a direct application of ADS.

We next consider preconditioning approaches for the hybridized matrix H.
Based on the connection with mixed methods from the previous section we argue
that one suitable preconditioner for the hybridized Schur complement H in the
case of the H(÷) problem is an H^1-based AMG built from the matrix H.
The rest of this section is devoted by providing arguments for this choice, by
exploring the “near-nullspace” of H, i.e. the set of its low-frequency
eigenmodes, that is critical in AMG theory.

Let us consider a norm · for the space of Lagrange multiplier,

    ^2 =
      ∑_τ∈𝒯_h(h_τ^2_∂τ+
      1/h_τ-m_τ()^2_∂τ),

where 𝒯_h is a finite element mesh and m_τ() =
1/|∂τ|∫_∂τ ds.  It is proven in
<cit.> that H is spectrally equivalent to ·,
so we can deduce that H has a small near-nullspace containing only the constant
functions, as opposed to the large near-nullspace of the original system A.
Therefore, the hybridized Schur complement H is “H^1-like”, which
motivates the use of H^1-based AMG for preconditioning.

In fact, by exploiting H, one can see that H has the explicit form:

    H = H = C(^-1-^-1^T(+^-1^T)^-1^-1)C^T.

In the above expression,  is an L^2 mass matrix while ^-1^T
resembles a weighted discrete Laplacian. Hence, ^-1^T ≫ and
so the near-nullspace of H is characterized by the nullspace of the matrix
C C^T, where

    = ^-1-^-1^T(^-1^T)^-1^-1.

Note that for any ∈(^T), there exists  such that =^T, so

    = ^-1-^-1^T(^-1^T)^-1^-1^T =
      ^-1-^-1^T = 0.

Thus, we have

    (^T)⊆().

On the other hand, let ∈(). Then

    ^-1-^-1^T(^-1^T)^-1^-1 = 0.

Consequently, =^T(^-1^T)^-1^-1. So for any
∈(), we have

    <,>     = <^T(^-1^T)^-1^-1,>
        = <(^-1^T)^-1^-1,> = 0

Hence,

    ()⊆()^⊥ = (^T).

The above inclusion relation and (<ref>) imply that

    () = (^T).

Therefore, in order to know what (C C^T) is, we should find out

    (C^T) ∩(^T).

Note that the constraint matrix C enforces the continuity of the master and
slave dofs on the interface (cf. Example <ref>), so
(C^T) is orthogonal to the subspace spanned by all the interior dofs. In
other words,  is in (C^T) only if every entry corresponding to an
interior dof is zero. So we are looking for the vectors in (^T) such
that entries corresponding to interior dofs are zero. Recall that  comes
from the bilinear form (÷,p). In the case of Raviart-Thomas elements,
constants are the only p such that (÷,p)=0 for all the 
associated with interior dofs (bubbles). Hence, we deduce that

    (C C^T) = {   |  C^T∈span{^T} },

where  is the coefficient vector of the constant 1 function in the
corresponding discrete space. This clearly shows that
((C C^T))≤ 1.




§ IMPLEMENTATION DETAILS


In this section we provide some details for the practical application of the
presented approaches as implemented in the finite element library MFEM
<cit.>.

The implementation is contained in class Hybridization which can be used
either directly, or more conveniently through the BilinearForm class.  The
prolongation operator P is not directly constructed, instead, its action can
be performed using the methods

in class FiniteElementSpace.  The constraints matrix C is constructed as
illustrated in Example <ref>, based on a specified
multiplier finite element space and a given bilinear form.

During a construction stage, the element matrices are reordered and stored in
the 2 × 2 block form (<ref>), see Figure <ref> for
illustration. Recall that the i- and b-block unknowns correspond to the
element interior and boundary (or interface) degrees of freedom.  Any rows and
columns associated with essential degrees of freedom (e-dofs) are
eliminated. Algebraically, the interior dofs can be defined as those
corresponding to zero columns in the C matrix. In other words, C has the 1
× 2 block form

    C = [   0 C_b ].

As a consequence, the matrix H of the hybridized system
(<ref>) can be computed as

    H = C A^-1 C^T = C_b S_b^-1 C_b^T  ,

where S_b = A_bb - A_biA_ii^-1A_ib is the
Schur complement of A which can be formed independently for each
element in the mesh. In MFEM, this is performed using element-wise block LU
factorizations of the form

    [ A_ii A_ib; A_bi A_bb ]
    =
    [  L_i    0; L_bi  I_b ][  U_i U_ib;    0  S_b ]

to compute the Schur complements S_b which are then themselves
LU-factorized, allowing simple multiplication of any vector or matrix with
S_b^-1.

Note that for Robin-type boundary conditions, one needs to assemble
contributions to the global system matrix A coming from integration over
boundary faces. In MFEM, such face matrices are assembled into the element
matrix of the element adjacent to the corresponding face. The rest of the
hybridization procedure remains the same.







§ NUMERICAL TESTS


To validate the parallel scalability of the proposed solvers, in this section we
present several weak scaling tests, where the number of degrees of freedom per
processor is held to be about the same while the number of processors is
increased. All the experiments except the ones in
Section <ref> are performed on the cluster Sierra at
Lawrence Livermore National Laboratory (LLNL). Sierra has a total of 1944 nodes
(Intel Xeon EP X5660 @ 2.80 GHz), each of which has 12 cores and 24 GB of
memory. The numerical results are generated using the libraries MFEM <cit.>
and hypre <cit.> developed at LLNL.




 §.§ Soft hard materials


Consider solving (<ref>) on the unit cube, i.e. Ω = [0,1]^3,
with the source function  being [1, 1, 1]^T. Let Ω_i =
[1/4,1/2]^3∪[1/2,3/4]^3. We consider α≡1 in
Ω, whereas

    β={[     1 Ω\Ω_i;  10^p   Ω_i ].,  p ∈{ -8, -4, 0, 4, 8}.


This example is intended to test the robustness of the solvers against
coefficient jumps. We discretize (<ref>) by Raviart-Thomas elements
(RT_k) on hexahedral meshes.  In order to better contrast the effect of higher
order discretizations, the initial meshes in the weak scaling tests are chosen
so that the problem size is the same for discretizations of different orders,
see Table <ref>.

The PCG stopping criterion is a reduction in the l_2 residual by a factor of
10^-12.  The time to solution and number of PCG iterations (in parentheses)
for the proposed approaches, as well as ADS, are reported in
Tables <ref>–<ref>. Here, the time to
solution for the hybridization approach includes the formation time of the
hybridization Schur complement H, the setup time of AMG (BoomerAMG
<cit.>), PCG solving time, and the time for recovering the original
solution by back substitution.  Similarly, the time to solution for the static
condensation approach includes the formation time of S
(cf. Section <ref>), ADS setup time, PCG solving
time, and back substitution. Lastly, the time to solution for “ADS” is simply
ADS setup time and PCG solving time.



Our results show that all the solvers have good weak scaling and robustness with
respect to coefficient jumps. In all the cases, while both the hybridization and
static condensation solvers give shorter solving time than ADS, the
hybridization approach is the fastest. Moreover, as the finite element order
goes higher, discrepancy in the solution time between the solvers becomes more
significant even though the problem size is the same. In particular, for the
RT_3 discretization, hybridization is in general 4 times faster than static
condensation, and 20 times faster than ADS.













 §.§ The crooked pipe problem


This example deals with large coefficient jumps and highly stretched elements,
see Figure <ref>. The coefficients α and β are both
piecewise-constant functions such that

    (α, β)={[    (1.641,0.2)               ; (0.00188,2000)                ].


The setup of the experiments is the same as in Section <ref>, and
the weak scaling results are reported in Table <ref>. Note that
we were running out of memory when setting up ADS for RT_3 (indicated by
“∗"), which shows that the hybridization and static condensation
approaches are more memory-friendly for high order discretizations.  The results
in Table <ref> again suggest that hybridization with AMG is the
most efficient solver among the three for this particular problem; it is also
observed that the iteration count in the hybridization approach is more steady
than the other two.






 §.§ Application to radiation diffusion




In this section we compare the proposed approaches when they are used for
solving the radiation diffusion equations in flux form, see <cit.>.
This model is often used in astrophysical and inertial confinement fusion (ICF)
simulations. We consider a problem with two materials (k=1,2 below), indicated
by blue and red colors in Figure <ref>, and there are no
mixed zones.  The governing equations are

    ρ_k de_k/dt   = c σ_k (E-B(T(e_k))), 
    dE/dt+∇· F     = -c ∑_k σ_k(E-B(T(e_k))), 
    1/3∇ E    = -∑_k σ_k/c F,

subject to the boundary conditions specified in Figure <ref>.
Here e_k, ρ_k, σ_k are the internal energy, density and opacity,
respectively, of material k, E is radiation energy, F is radiation flux,
B is the Planck's function for black-body radiation, and c is the speed of
light.  The temperature boundary condition heats up one of the ends of the pipe,
causing an increase of the radiation field E and its diffusion throughout the
domain.

The discrete approximations of E and e_k are in the finite elements space
Q_r ⊂ L^2, consisting of piecewise polynomials of degree r, while the
flux F is approximated in the Raviart-Thomas space RT_r.  For stability
reasons the time stepping of e_k and E is implicit, which results in a
nonlinear problem in each time step, as B(T(e_k)) is nonlinear.  Newton's
method is used to solve the resulting system. As E and e_k are
discontinuous, they can be eliminated locally during each Newton iteration,
which is inexpensive and inherently parallel. Upon eliminating E and e_k, an
H(÷) linear system for F is obtained, which we use to test the
performance of the two proposed approaches, as well as AMS (2D) and ADS (3D).









Because each method ultimately solves a different linear system, this example is
used to obtain better perception of the execution time needed for each method to
obtain similar solution accuracy with respect to the original system.  Starting
with identical initial conditions, each computation performs exactly one linear
solve, and then the F residual is calculated from the definition of the
nonlinear system. Execution times are presented for computations reaching
similar residuals. For this particular example, this is achieved by setting the
relative tolerances to 1e-6, 1e-5 and 1e-12 for AMS/ADS, static
condensation+AMS/ADS, and hybridization+AMG, respectively.

As observed in Tables <ref> and <ref>, even though the
hybridized system must be solved to a much lower relative tolerance, the
hybridization approach still has the shortest solving time, followed by the
static condensation and then ADS.  Except for RT_0 in the 2D case, the
hybridization approach is at least 10 times faster than AMS/ADS, which is a
substantial improvement.






 §.§ The SPE10 benchmark


Our last numerical example is the SPE10 benchmark coefficient <cit.> from
the reservoir simulation community. This coefficient is well-known for its
multi-scale heterogeneity and anisotropy. The initial mesh is a
60×220×85 hexahedral grid, with the size of each cell being
20×10×2. β is taken to be the SPE10 coefficient (see
Figure <ref>), and α≡1. The time spent on different
components of the solution process for the hybridization approach and ADS are
shown in Tables <ref> and <ref> respectively.  This
example also demonstrates the superiority in terms of performance of the
hybridization solver.








§ CONCLUSION
 

In this paper, two popular techniques in finite elements, hybridization and
static condensation, were studied and compared. Both of these two techniques
lead to a reduced system of the same size. The formation of the reduced system
involves inversion of independent local matrices, which can be done in
parallel. We also discuss suitable preconditioning strategies using various
algebraic multigrid methods for the respective reduced systems. In particular,
we demonstrate that if BoomerAMG/AMS/ADS work well for the original problem,
then they also work well for the reduced systems obtained from static
condensation. On the other hand, AMG for H^1-equivalent problems is actually a
well-suited preconditioner for the reduced system resulting from hybridization
of H(÷) problems. We presented several numerical examples comparing the
performance of the proposed approaches, as well as solving the original problem
directly by ADS/AMS. Our results show good weak scaling of all of the proposed
approaches. In general, for H(÷) problems, the approach using hybridization
and AMG is more efficient than the approach using static condensation and
ADS/AMS, while both of these approaches are faster than a direct application of
ADS/AMS to the original problem. In fact, we observed substantial savings in the
solve time if the hybridization approach is used, especially for higher order
discretizations. Implementation of the two approaches are freely available in
MFEM <cit.>.


siam    
                
                



       

