# AMS¶

AMS (the Auxiliary-space Maxwell Solver) is a parallel unstructured Maxwell solver for edge finite element discretizations of the variational problem

Here \({\mathbf V}_h\) is the lowest order Nedelec (edge) finite element space, and \(\alpha>0\) and \(\beta \ge 0\) are scalar, or SPD matrix coefficients. AMS was designed to be scalable on problems with variable coefficients, and allows for \(\beta\) to be zero in part or the whole domain. In either case the resulting problem is only semidefinite, and for solvability the right-hand side should be chosen to satisfy compatibility conditions.

AMS is based on the auxiliary space methods for definite Maxwell problems proposed in [HiXu2006]. For more details, see [KoVa2009].

## 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,

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

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.

The representations of the constant vector fields \((1,0,0)\), \((0,1,0)\), and \((0,0,1)\) in the \({\mathbf V}_h\) basis, given as three vectors: \(G_x\), \(G_y\), and \(G_z\). Note that since no boundary conditions are imposed on \(G\), the above vectors can be computed as \(G_x = G x\), \(G_y = G y\) and \(G_z = G z\), where \(x\), \(y\), and \(z\) are vectors representing the coordinates of the vertices of the mesh.

In addition to the above quantities, AMS can utilize the following (optional) information:

The Poisson matrices \(A_\alpha\) and \(A_\beta\), corresponding to assembling of the forms \((\alpha\, \nabla u, \nabla v)+(\beta\, \nabla u, \nabla v)\) and \((\beta\, \nabla u, \nabla v)\) using standard linear finite elements on the same mesh.

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

\(A_G\) – a matrix associated with the mass term which is either \(G^T {\mathbf A} G\) or the Poisson matrix \(A_\beta\) (if given).

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

\({\mathbf A}_{{\mathbf \Pi}}\) – a matrix associated with the stiffness term which is either \({\mathbf \Pi}^{\,T} {\mathbf A} {\mathbf \Pi}\) or a block-diagonal matrix with diagonal blocks \(A_\alpha\) (if given).

\(B_G\) and \({\mathbf B}_{{\mathbf \Pi}}\) – efficient (AMG) solvers for \(A_G\) and \({\mathbf A}_{{\mathbf \Pi}}\).

The solution procedure then is a three-level method using smoothing in the original edge space and subspace corrections based on \(B_G\) 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. If \(\beta\) is identically zero one can skip the subspace correction associated with \(G\), in which case the solver is a two-level method.

## Sample Usage¶

AMS 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. See example
code `ex15.c`

for a complete implementation. We start with the allocation of
the `HYPRE_Solver`

object:

```
HYPRE_Solver solver;
HYPRE_AMSCreate(&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.

AMS offers the option to set the space dimension. By default we consider the dimension to be \(3\). The only other option is \(2\), and it can be set with the function given below. We note that a 3D solver will still work for a 2D problem, but it will be slower and will require more memory than necessary.

```
HYPRE_AMSSetDimension(solver, dim);
```

The user is required to provide the discrete gradient matrix \(G\). AMS
expects a matrix defined on the whole mesh with no boundary edges/nodes
excluded. It is essential to **not** impose any boundary conditions on
\(G\). Regardless of which hypre conceptual interface was used to construct
\(G\), one can obtain a ParCSR version of it. This is the expected format in
the following function.

```
HYPRE_AMSSetDiscreteGradient(solver, G);
```

In addition to \(G\), we need one additional piece of information in order to construct the solver. The user has the option to choose either the coordinates of the vertices in the mesh or the representations of the constant vector fields in the edge element basis. In both cases three hypre parallel vectors should be provided. For 2D problems, the user can set the third vector to NULL. The corresponding function calls read:

```
HYPRE_AMSSetCoordinateVectors(solver,x,y,z);
```

or

```
HYPRE_AMSSetEdgeConstantVectors(solver, one_zero_zero, zero_one_zero, zero_zero_one);
```

The vectors `one_zero_zero`

, `zero_one_zero`

and `zero_zero_one`

above
correspond to the constant vector fields \((1,0,0)\), \((0,1,0)\) and
\((0,0,1)\).

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

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

The available cycle types in AMS 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 BoomerAMG based on \(B_G\), 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 AMS 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_AMSSetMaxIter(solver, maxit); /* default value: 20 */
HYPRE_AMSSetTol(solver, tol); /* default value: 1e-6 */
HYPRE_AMSSetPrintLevel(solver, print); /* default value: 1 */
```

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

```
HYPRE_AMSSetSmoothingOptions(solver, 2, 1, 1.0, 1.0);
HYPRE_AMSSetAlphaAMGOptions(solver, 10, 1, 3, 0.25, 0, 0);
HYPRE_AMSSetBetaAMGOptions(solver, 10, 1, 3, 0.25, 0, 0);
```

For (singular) problems where \(\beta = 0\) in the whole domain, different (in fact simpler) version of the AMS solver is offered. To allow for this simplification, use the following hypre call

```
HYPRE_AMSSetBetaPoissonMatrix(solver, NULL);
```

If \(\beta\) is zero only in parts of the domain, the problem is still singular, but the AMS solver will try to detect this and construct a non-singular preconditioner. Though this often works well in practice, AMS also provides a more robust version for solving such singular problems to very low convergence tolerances. This version takes advantage of additional information: the list of nodes which are interior to a zero-conductivity region provided by the function

```
HYPRE_AMSSetInteriorNodes(solver, HYPRE_ParVector interior_nodes);
```

A node is interior, if its entry in the `interior_nodes`

array is \(1.0\).
Based on this array, a restricted discrete gradient operator \(G_0\) is
constructed, and AMS is then defined based on the matrix \({\mathbf
A}+\delta G_0^TG_0\) which is non-singular, and a small \(\delta>0\)
perturbation of \({\mathbf A}\). When iterating with this preconditioner, it
is advantageous to project on the compatible subspace \(Ker(G_0^T)\). This
can be done periodically, or manually through the functions

```
HYPRE_AMSSetProjectionFrequency(solver, int projection_frequency);
HYPRE_AMSProjectOutGradients(solver, HYPRE_ParVector x);
```

Two additional matrices are constructed in the setup of the AMS method—one
corresponding to the coefficient \(\alpha\) and another corresponding to
\(\beta\). This may lead to prohibitively high memory requirements, and the
next two function calls may help to save some memory. For example, if the
Poisson matrix with coefficient \(\beta\) (denoted by `Abeta`

) is
available then one can avoid one matrix construction by calling

```
HYPRE_AMSSetBetaPoissonMatrix(solver, Abeta);
```

Similarly, if the Poisson matrix with coefficient \(\alpha\) is available
(denoted by `Aalpha`

) the second matrix construction can also be avoided by
calling

```
HYPRE_AMSSetAlphaPoissonMatrix(solver, Aalpha);
```

Note the following regarding these functions:

Both of them change their input. More specifically, the diagonal entries of the input matrix corresponding to eliminated degrees of freedom (due to essential boundary conditions) are penalized.

It is assumed that their essential boundary conditions of \({\mathbf A}\),

`Abeta`

and`Aalpha`

are on the same part of the boundary.`HYPRE_AMSSetAlphaPoissonMatrix`

forces the AMS method to use a simpler, but weaker (in terms of convergence) method. With this option, the multiplicative AMS cycle is not guaranteed to converge with the default parameters. The reason for this is the fact the solver is not variationally obtained from the original matrix (it utilizes the auxiliary Poisson–like matrices`Abeta`

and`Aalpha`

). Therefore, it is recommended in this case to use AMS as preconditioner only.

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 AMS setup.) The setup call reads,

```
HYPRE_AMSSetup(solver, A, b, x);
```

It is important to note the order of the calling sequence. For example, do
**not** call `HYPRE_AMSSetup`

before calling `HYPRE_AMSSetDiscreteGradient`

and one of the functions `HYPRE_AMSSetCoordinateVectors`

or
`HYPRE_AMSSetEdgeConstantVectors`

.

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

```
HYPRE_AMSSolve(solver, A, b, x);
```

Finally, the solver can be destroyed with

```
HYPRE_AMSDestroy(&solver);
```

More details can be found in the files `ams.h`

and `ams.c`

located in the
`parcsr_ls`

directory.

## High-order Discretizations¶

In addition to the interface for the lowest-order Nedelec elements described in the previous subsections, AMS also provides support for (arbitrary) high-order Nedelec element discretizations. Since the robustness of AMS depends on the performance of 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 AMS interface, the user does not need to provide the
coordinates of the vertices (or the representations of the constant vector
fields in the edge basis), but instead should construct and pass the Nedelec
interpolation matrix \({\mathbf \Pi}\) which maps (high-order) vector nodal
finite elements into the (high-order) Nedelec space. In other words,
\({\mathbf \Pi}\) is the (parallel) matrix representation of the
interpolation mapping from \(\mathrm{P}_k^3\)/\(\mathrm{Q}_k^3\) into
\(\mathrm{ND}_k\), see [HiXu2006], [KoVa2009]. We require this matrix as
an input, since in the high-order case its entries very much depend on the
particular choice of the basis functions in the edge and nodal spaces, as well
as on the geometry of the mesh elements. The columns of \({\mathbf \Pi}\)
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`

).

Similarly to the Nedelec interpolation, the discrete gradient matrix \(G\) should correspond to the mapping \(\varphi \in \mathrm{P}_k^3 / \mathrm{Q}_k^3 \mapsto \nabla \varphi \in \mathrm{ND}_k\), so even though its 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_AMSSetDimension(solver, dim);
HYPRE_AMSSetDiscreteGradient(solver, G);
HYPRE_AMSSetInterpolations(solver, 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 such as those based on the second
family of Nedelec elements), but we recommend calling the previously described
`HYPRE_AMSSetCoordinateVectors`

instead, since this allows AMS to handle the
construction and use of \({\mathbf \Pi}\) internally.

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

```
HYPRE_AMSSetInterpolations(solver, NULL, Pix, Piy, Piz);
```

which enables the use of AMS cycle types with index greater than 10. By definition, \({\mathbf \Pi}^x \varphi = {\mathbf \Pi} (\varphi,0,0)\), and similarly for \({\mathbf \Pi}^y\) and \({\mathbf \Pi}^z\). Each of these matrices has the same sparsity pattern as \(G\), but their entries depend on the coordinates of the mesh vertices.

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

```
HYPRE_AMSSetInterpolations(solver, Pi, Pix, Piy, Piz);
```

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

## Non-conforming AMR Grids¶

AMS could also be applied to problems with adaptive mesh refinement (AMR) posed on non-conforming quadrilateral/hexahedral meshes, see [GrKo2015] for more details.

On non-conforming grids (assuming also arbitrarily high-order elements), each
finite element space has two versions: a conforming one,
e.g. \(\mathrm{Q}_k^{c} / \mathrm{ND}_k^c\), where the *hanging* degrees of
freedom are constrained by the conforming (*real*) degrees of freedom, and a
non-conforming one, e.g. \(\mathrm{Q}_k^{nc} / \mathrm{ND}_k^{nc}\) where
the non-conforming degrees of freedom (hanging and real) are unconstrained.
These spaces are related with the conforming prolongation and the pure
restriction operators \(P\) and \(R\), as well as the conforming and
non-conforming version of the discrete gradient operator as follows:

Since the linear system is posed on \(\mathrm{ND}_k^c\), the user needs to
provide the conforming discrete gradient matrix \(G_c\) to AMS, using
`HYPRE_AMSSetDiscreteGradient`

. This matrix is defined by the requirement
that the above diagram commutes from \(\mathrm{Q}_k^{c}\) to
\(\mathrm{ND}_k^{nc}\), corresponding to the definition

i.e. the conforming gradient is computed by starting with a conforming nodal \(\mathrm{Q}_k\) function, interpolating it in the hanging nodes, computing the gradient locally and representing it in the Nedelec space on each element (the non-conforming discrete gradient \(G_{nc}\) in the above formula), and disregarding the values in the hanging \(\mathrm{ND}_k\) degrees of freedom.

Similar considerations imply that the conforming Nedelec interpolation matrix \({\mathbf \Pi}_{c}\) should be defined as

with \({\mathbf \Pi}_{nc}\) computed element-wise as in the previous subsection. Note that in the low-order case, \({\mathbf \Pi}_{c}\) can be computed internally in AMS based only \(G_c\) and the conforming coordinates of the vertices \(x_c\)/\(y_c\)/\(z_c\), see [GrKo2015].