ADS

The Auxiliary-space Divergence Solver (ADS) is a parallel unstructured solver similar to AMS, but targeting \(H(div)\) instead of \(H(curl)\) problems. Its usage and options are very similar to those of AMS, and in general the relationship between ADS and AMS is analogous to that between AMS and AMG.

Specifically ADS was designed for the scalable solution of linear systems arising in the finite element discretization of the variational problem

(1)\[\mbox{Find } {\mathbf u} \in {\mathbf W}_h \>:\qquad (\alpha\, \nabla \cdot {\mathbf u}, \nabla \cdot {\mathbf v}) + (\beta\, {\mathbf u}, {\mathbf v}) = ({\mathbf f}, {\mathbf v})\,, \qquad \mbox{for all } {\mathbf v} \in {\mathbf W}_h \,,\]

where \({\mathbf W}_h\) is the lowest order Raviart-Thomas (face) finite element space, and \(\alpha>0\) and \(\beta>0\) are scalar, or SPD matrix variable coefficients. It is based on the auxiliary space methods for \(H(div)\) problems proposed in [HiXu2006].

Overview

Let \({\mathbf A}\) and \({\mathbf b}\) be the stiffness matrix and the load vector corresponding to (1). Then the resulting linear system of interest reads,

(2)\[{\mathbf A}\, {\mathbf x} = {\mathbf b} \,.\]

The coefficients \(\alpha\) and \(\beta\) are naturally associated with the “stiffness” and “mass” terms of \({\mathbf A}\). Besides \({\mathbf A}\) and \({\mathbf b}\), ADS requires the following additional user input:

  1. The discrete curl matrix \(C\) representing the faces of the mesh in terms of its edges. \(C\) has as many rows as the number of faces in the mesh, with each row having nonzero entries \(+1\) and \(-1\) in the columns corresponding to the edges composing the face. The sign is determined based on the orientation of the edges relative to the face. We require that \(C\) includes all (interior and boundary) faces and edges.

  2. The discrete gradient matrix \(G\) representing the edges of the mesh in terms of its vertices. \(G\) has as many rows as the number of edges in the mesh, with each row having two nonzero entries: \(+1\) and \(-1\) in the columns corresponding to the vertices composing the edge. The sign is determined based on the orientation of the edge. We require that \(G\) includes all (interior and boundary) edges and vertices.

  3. Vectors \(x\), \(y\), and \(z\) representing the coordinates of the vertices of the mesh.

Internally, ADS proceeds with the construction of the following additional objects:

  • \(A_C\) – the curl-curl matrix \(C^{\,T} {\mathbf A} C\).

  • \({\mathbf \Pi}\) – the matrix representation of the interpolation operator from vector linear to face finite elements.

  • \({\mathbf A}_{{\mathbf \Pi}}\) – the vector nodal matrix \({\mathbf \Pi}^{\,T} {\mathbf A} {\mathbf \Pi}\).

  • \(B_C\) and \({\mathbf B}_{{\mathbf \Pi}}\) – efficient (AMS/AMG) solvers for \(A_C\) and \({\mathbf A}_{{\mathbf \Pi}}\).

The solution procedure then is a three-level method using smoothing in the original face space and subspace corrections based on \(B_C\) and \({\mathbf B}_{{\mathbf \Pi}}\). We can employ a number of options here utilizing various combinations of the smoother and solvers in additive or multiplicative fashion.

Sample Usage

ADS can be used either as a solver or as a preconditioner. Below we list the sequence of hypre calls needed to create and use it as a solver. We start with the allocation of the HYPRE_Solver object:

HYPRE_Solver solver;
HYPRE_ADSCreate(&solver);

Next, we set a number of solver parameters. Some of them are optional, while others are necessary in order to perform the solver setup.

The user is required to provide the discrete curl and gradient matrices \(C\) and \(G\). ADS expects a matrix defined on the whole mesh with no boundary faces, edges or nodes excluded. It is essential to not impose any boundary conditions on \(C\) or \(G\). Regardless of which hypre conceptual interface was used to construct the matrices, one can always obtain a ParCSR version of them. This is the expected format in the following functions.

HYPRE_ADSSetDiscreteCurl(solver, C);
HYPRE_ADSSetDiscreteGradient(solver, G);

Next, ADS requires the coordinates of the vertices in the mesh as three hypre parallel vectors. The corresponding function call reads:

HYPRE_ADSSetCoordinateVectors(solver, x, y, z);

The remaining solver parameters are optional. For example, the user can choose a different cycle type by calling

HYPRE_ADSSetCycleType(solver, cycle_type); /* default value: 1 */

The available cycle types in ADS are:

  • cycle_type=1: multiplicative solver \((01210)\)

  • cycle_type=2: additive solver \((0+1+2)\)

  • cycle_type=3: multiplicative solver \((02120)\)

  • cycle_type=4: additive solver \((010+2)\)

  • cycle_type=5: multiplicative solver \((0102010)\)

  • cycle_type=6: additive solver \((1+020)\)

  • cycle_type=7: multiplicative solver \((0201020)\)

  • cycle_type=8: additive solver \((0(1+2)0)\)

  • cycle_type=11: multiplicative solver \((013454310)\)

  • cycle_type=12: additive solver \((0+1+3+4+5)\)

  • cycle_type=13: multiplicative solver \((034515430)\)

  • cycle_type=14: additive solver \((01(3+4+5)10)\)

Here we use the following convention for the three subspace correction methods: \(0\) refers to smoothing, \(1\) stands for AMS based on \(B_C\), and \(2\) refers to a call to BoomerAMG for \({\mathbf B}_{{\mathbf \Pi}}\). The values \(3\), \(4\) and \(5\) refer to the scalar subspaces corresponding to the \(x\), \(y\) and \(z\) components of \(\mathbf \Pi\).

The abbreviation \(xyyz\) for \(x,y,z \in \{0,1,2,3,4,5\}\) refers to a multiplicative subspace correction based on solvers \(x\), \(y\), \(y\), and \(z\) (in that order). The abbreviation \(x+y+z\) stands for an additive subspace correction method based on \(x\), \(y\) and \(z\) solvers. The additive cycles are meant to be used only when ADS is called as a preconditioner. In our experience the choices cycle_type=1,5,8,11,13 often produced fastest solution times, while cycle_type=7 resulted in smallest number of iterations.

Additional solver parameters, such as the maximum number of iterations, the convergence tolerance and the output level, can be set with

HYPRE_ADSSetMaxIter(solver, maxit);     /* default value: 20 */
HYPRE_ADSSetTol(solver, tol);           /* default value: 1e-6 */
HYPRE_ADSSetPrintLevel(solver, print);  /* default value: 1 */

More advanced parameters, affecting the smoothing and the internal AMS and AMG solvers, can be set with the following three functions:

HYPRE_ADSSetSmoothingOptions(solver, 2, 1, 1.0, 1.0);
HYPRE_ADSSetAMSOptions(solver, 11, 10, 1, 3, 0.25, 0, 0);
HYPRE_ADSSetAMGOptions(solver, 10, 1, 3, 0.25, 0, 0);

We note that the AMS cycle type, which is the second parameter of HYPRE_ADSSetAMSOptions should be greater than 10, unless the high-order interface of HYPRE_ADSSetInterpolations described in the next subsection is being used.

After the above calls, the solver is ready to be constructed. The user has to provide the stiffness matrix \({\mathbf A}\) (in ParCSR format) and the hypre parallel vectors \({\mathbf b}\) and \({\mathbf x}\). (The vectors are actually not used in the current ADS setup.) The setup call reads,

HYPRE_ADSSetup(solver, A, b, x);

It is important to note the order of the calling sequence. For example, do not call HYPRE_ADSSetup before calling each of the functions HYPRE_ADSSetDiscreteCurl, HYPRE_ADSSetDiscreteGradient and HYPRE_ADSSetCoordinateVectors.

Once the setup has completed, we can solve the linear system by calling

HYPRE_ADSSolve(solver, A, b, x);

Finally, the solver can be destroyed with

HYPRE_ADSDestroy(&solver);

More details can be found in the files ads.h and ads.c located in the parcsr_ls directory.

High-order Discretizations

Similarly to AMS, ADS also provides support for (arbitrary) high-order \(H(div)\) discretizations. Since the robustness of ADS depends on the performance of AMS and BoomerAMG on the associated (high-order) auxiliary subspace problems, we note that the convergence may not be optimal for large polynomial degrees \(k \geq 1\).

In the high-order ADS interface, the user does not need to provide the coordinates of the vertices, but instead should construct and pass the Raviart-Thomas and Nedelec interpolation matrices \({\mathbf \Pi}_{RT}\) and \({\mathbf \Pi}_{ND}\) which map (high-order) vector nodal finite elements into the (high-order) Raviart-Thomas and Nedelec space. In other words, these are the (parallel) matrix representation of the interpolation mappings from \(\mathrm{P}_k^3 / \mathrm{Q}_k^3\) into \(\mathrm{RT}_{k-1}\) and \(\mathrm{ND}_k\), see [HiXu2006], [KoVa2009]. We require these matrices as inputs, since in the high-order case their entries very much depend on the particular choice of the basis functions in the finite element spaces, as well as on the geometry of the mesh elements. The columns of the \({\mathbf \Pi}\) matrices should use a node-based numbering, where the \(x\)/\(y\)/\(z\) components of the first node (vertex or high-order degree of freedom) should be listed first, followed by the \(x\)/\(y\)/\(z\) components of the second node and so on (see the documentation of HYPRE_BoomerAMGSetDofFunc). Furthermore, each interpolation matrix can be split into \(x\), \(y\) and \(z\) components by defining \({\mathbf \Pi}^x \varphi = {\mathbf \Pi} (\varphi,0,0)\), and similarly for \({\mathbf \Pi}^y\) and \({\mathbf \Pi}^z\).

The discrete gradient and curl matrices \(G\) and \(C\) should correspond to the mappings \(\varphi \in \mathrm{P}_k^3 / \mathrm{Q}_k^3 \mapsto \nabla \varphi \in \mathrm{ND}_k\) and \({\mathbf u} \in \mathrm{ND}_k \mapsto \nabla \times {\mathbf u} \in \mathrm{RT}_{k-1}\), so even though their values are still independent of the mesh coordinates, they will not be \(\pm 1\), but will be determined by the particular form of the high-order basis functions and degrees of freedom.

With these matrices, the high-order setup procedure is simply

HYPRE_ADSSetDiscreteCurl(solver, C);
HYPRE_ADSSetDiscreteGradient(solver, G);
HYPRE_ADSSetInterpolations(solver, RT_Pi, NULL, NULL, NULL,
                                   ND_Pi, NULL, NULL, NULL);

We remark that the above interface calls can also be used in the lowest-order case (or even other types of discretizations), but we recommend calling the previously described HYPRE_ADSSetCoordinateVectors instead, since this allows ADS to handle the construction and use of the interpolations internally.

Specifying the monolithic \({\mathbf \Pi}_{RT}\) limits the ADS cycle type options to those less than 10. Alternatively one can separately specify the \(x\), \(y\) and \(z\) components of \({\mathbf \Pi}_{RT}\).

HYPRE_ADSSetInterpolations(solver, NULL, RT_Pix, RT_Piy, RT_Piz,
                                   ND_Pi, NULL, NULL, NULL);

which enables the use of ADS cycle types with index greater than 10. The same holds for \({\mathbf \Pi}_{ND}\) and its components, e.g. to enable the subspace AMS cycle type greater then 10 we need to call

HYPRE_ADSSetInterpolations(solver, NULL, RT_Pix, RT_Piy, RT_Piz,
                                   NULL, ND_Pix, ND_Piy, ND_Piz);

Finally, both \({\mathbf \Pi}\) and their components can be passed to the solver:

HYPRE_ADSSetInterpolations(solver, RT_Pi, RT_Pix, RT_Piy, RT_Piz
                                   ND_Pi, ND_Pix, ND_Piy, ND_Piz);

which will duplicate some memory, but allows for experimentation with all available ADS and AMS cycle types.