RZ example with Robin conditions

Temporarily this means disabling partial assembly for the RT mass term.

The PA functionality is implemented in PR #3539.

.gitignore

@@ -307,6 +307,12 @@ miniapps/solvers/ParaView
 miniapps/solvers/mesh.*
 miniapps/solvers/sol.*
 
+miniapps/hdiv-linear-solver/darcy
+miniapps/hdiv-linear-solver/grad_div
+miniapps/hdiv-linear-solver/residual
+miniapps/hdiv-linear-solver/rz
+miniapps/hdiv-linear-solver/ParaView
+
 miniapps/parelag/MultilevelHcurlHdivSolver
 miniapps/parelag/*.mesh
 

fem/bilininteg_mass_pa.cpp

@@ -56,8 +56,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    dim = mesh->Dimension();
    ne = fes.GetMesh()->GetNE();
    nq = ir->GetNPoints();
-   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::COORDINATES |
-                                    GeometricFactors::JACOBIANS, mt);
+   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS, mt);
    maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
@@ -74,7 +73,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
       const bool const_c = coeff.Size() == 1;
       const bool by_val = map_type == FiniteElement::VALUE;
       const auto W = Reshape(ir->GetWeights().Read(), Q1D,Q1D);
-      const auto J = Reshape(geom->J.Read(), Q1D,Q1D,2,2,NE);
+      const auto J = Reshape(geom->detJ.Read(), Q1D,Q1D,NE);
       const auto C = const_c ? Reshape(coeff.Read(), 1,1,1) :
                      Reshape(coeff.Read(), Q1D,Q1D,NE);
       auto v = Reshape(pa_data.Write(), Q1D,Q1D, NE);
@@ -84,11 +83,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
          {
             MFEM_FOREACH_THREAD(qy,y,Q1D)
             {
-               const double J11 = J(qx,qy,0,0,e);
-               const double J12 = J(qx,qy,1,0,e);
-               const double J21 = J(qx,qy,0,1,e);
-               const double J22 = J(qx,qy,1,1,e);
-               const double detJ = (J11*J22)-(J21*J12);
+               const double detJ = J(qx,qy,e);
                const double coeff = const_c ? C(0,0,0) : C(qx,qy,e);
                v(qx,qy,e) =  W(qx,qy) * coeff * (by_val ? detJ : 1.0/detJ);
             }
@@ -102,7 +97,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
       const bool const_c = coeff.Size() == 1;
       const bool by_val = map_type == FiniteElement::VALUE;
       const auto W = Reshape(ir->GetWeights().Read(), Q1D,Q1D,Q1D);
-      const auto J = Reshape(geom->J.Read(), Q1D,Q1D,Q1D,3,3,NE);
+      const auto J = Reshape(geom->detJ.Read(), Q1D,Q1D,Q1D,NE);
       const auto C = const_c ? Reshape(coeff.Read(), 1,1,1,1) :
                      Reshape(coeff.Read(), Q1D,Q1D,Q1D,NE);
       auto v = Reshape(pa_data.Write(), Q1D,Q1D,Q1D,NE);
@@ -114,18 +109,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
             {
                MFEM_FOREACH_THREAD(qz,z,Q1D)
                {
-                  const double J11 = J(qx,qy,qz,0,0,e);
-                  const double J21 = J(qx,qy,qz,1,0,e);
-                  const double J31 = J(qx,qy,qz,2,0,e);
-                  const double J12 = J(qx,qy,qz,0,1,e);
-                  const double J22 = J(qx,qy,qz,1,1,e);
-                  const double J32 = J(qx,qy,qz,2,1,e);
-                  const double J13 = J(qx,qy,qz,0,2,e);
-                  const double J23 = J(qx,qy,qz,1,2,e);
-                  const double J33 = J(qx,qy,qz,2,2,e);
-                  const double detJ = J11 * (J22 * J33 - J32 * J23) -
-                  /* */               J21 * (J12 * J33 - J32 * J13) +
-                  /* */               J31 * (J12 * J23 - J22 * J13);
+                  const double detJ = J(qx,qy,qz,e);
                   const double coeff = const_c ? C(0,0,0,0) : C(qx,qy,qz,e);
                   v(qx,qy,qz,e) = W(qx,qy,qz) * coeff * (by_val ? detJ : 1.0/detJ);
                }

fem/dgmassinv.hpp

@@ -87,6 +87,8 @@ public:
    ///
    /// If @ref iterative_mode is @a true, @a u is used as an initial guess.
    void Mult(const Vector &b, Vector &u) const;
+   /// Same as Mult() since the mass matrix is symmetric.
+   void MultTranspose(const Vector &b, Vector &u) const { Mult(b, u); }
    /// Not implemented. Aborts.
    void SetOperator(const Operator &op);
    /// Set the relative tolerance.

fem/restriction.hpp

@@ -53,17 +53,6 @@ protected:
    Array<int> indices;
    Array<int> gather_map;
 
-   friend class BatchedLORAssembly;
-   friend class BatchedLOR_ADS;
-   friend class BatchedLOR_AMS;
-
-   /// @name Low-level access to the underlying element-dof mappings
-   ///@{
-   const Array<int> &GatherMap() const { return gather_map; }
-   const Array<int> &Indices() const { return indices; }
-   const Array<int> &Offsets() const { return offsets; }
-   ///@}
-
 public:
    ElementRestriction(const FiniteElementSpace&, ElementDofOrdering);
    void Mult(const Vector &x, Vector &y) const override;
@@ -102,6 +91,13 @@ public:
    /// Performs either MultTranspose or AddMultTranspose depending on the
    /// boolean template parameter @a ADD.
    template <bool ADD> void TAddMultTranspose(const Vector &x, Vector &y) const;
+
+   /// @name Low-level access to the underlying element-dof mappings
+   ///@{
+   const Array<int> &GatherMap() const { return gather_map; }
+   const Array<int> &Indices() const { return indices; }
+   const Array<int> &Offsets() const { return offsets; }
+   ///@}
 };
 
 /// Operator that converts L2 FiniteElementSpace L-vectors to E-vectors.

miniapps/hdiv-linear-solver/change_basis.cpp

@@ -0,0 +1,359 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "change_basis.hpp"
+#include "fem/qinterp/dispatch.hpp"
+#include "general/forall.hpp"
+#include "linalg/dtensor.hpp"
+
+namespace mfem
+{
+
+/// @brief Compute the inverse of the matrix A and store the result in Ainv.
+///
+/// The input A is an array of size n*n, interpreted as a matrix with column
+/// major ordering.
+void ComputeInverse(const Array<double> &A, Array<double> &Ainv)
+{
+   Array<double> A2 = A;
+   const int n2 = A.Size();
+   const int n = sqrt(n2);
+   Array<int> ipiv(n);
+   LUFactors lu(A2.GetData(), ipiv.GetData());
+   lu.Factor(n);
+   Ainv.SetSize(n2);
+   lu.GetInverseMatrix(n, Ainv.GetData());
+}
+
+void SubcellIntegrals(int n, const Poly_1D::Basis &basis, Array<double> &B)
+{
+   const IntegrationRule &ir = IntRules.Get(Geometry::SEGMENT, n);
+   const double *gll_pts = poly1d.GetPoints(n, BasisType::GaussLobatto);
+   Vector u(n);
+   B.SetSize(n*n);
+   B = 0.0;
+
+   for (int i = 0; i < n; ++i)
+   {
+      const double h = gll_pts[i+1] - gll_pts[i];
+      // Loop over subcell quadrature points
+      for (int iq = 0; iq < ir.Size(); ++iq)
+      {
+         const IntegrationPoint &ip = ir[iq];
+         const double x = gll_pts[i] + h*ip.x;
+         const double w = h*ip.weight;
+         basis.Eval(x, u);
+         for (int j = 0; j < n; ++j)
+         {
+            B[i + j*n] += w*u[j];
+         }
+      }
+   }
+}
+
+void Transpose(const Array<double> &B, Array<double> &Bt)
+{
+   const int n = sqrt(B.Size());
+   Bt.SetSize(n*n);
+   for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) { Bt[i+j*n] = B[j+i*n]; }
+}
+
+ChangeOfBasis_L2::ChangeOfBasis_L2(FiniteElementSpace &fes)
+   : Operator(fes.GetTrueVSize()),
+     ne(fes.GetNE())
+{
+   auto *fec1 = dynamic_cast<const L2_FECollection*>(fes.FEColl());
+   MFEM_VERIFY(fec1, "Must be L2 finite element space");
+
+   const int btype = fec1->GetBasisType();
+
+   // If the basis types are the same, don't need to perform change of basis.
+   no_op = (btype == BasisType::IntegratedGLL);
+   if (no_op) { return; }
+
+   // Convert from the given basis to the "integrated GLL basis".
+   // The degrees of freedom are integrals over subcells.
+   const FiniteElement *fe = fes.GetFE(0);
+   auto *tbe = dynamic_cast<const TensorBasisElement*>(fe);
+   MFEM_VERIFY(tbe != nullptr, "Must be a tensor element.");
+   const Poly_1D::Basis &basis = tbe->GetBasis1D();
+
+   const int p = fes.GetMaxElementOrder();
+   const int pp1 = p + 1;
+
+   Array<double> B_inv;
+   SubcellIntegrals(pp1, basis, B_inv);
+
+   ComputeInverse(B_inv, B_1d);
+   Transpose(B_1d, Bt_1d);
+
+   // Set up the DofToQuad object, used in TensorValues
+   dof2quad.FE = fe;
+   dof2quad.mode = DofToQuad::TENSOR;
+   dof2quad.ndof = pp1;
+   dof2quad.nqpt = pp1;
+}
+
+void ChangeOfBasis_L2::Mult(const Vector &x, Vector &y) const
+{
+   if (no_op) { y = x; return; }
+   using namespace internal::quadrature_interpolator;
+   dof2quad.B.MakeRef(B_1d);
+   TensorValues<QVectorLayout::byVDIM>(ne, 1, dof2quad, x, y);
+}
+
+void ChangeOfBasis_L2::MultTranspose(const Vector &x, Vector &y) const
+{
+   if (no_op) { y = x; return; }
+   using namespace internal::quadrature_interpolator;
+   dof2quad.B.MakeRef(Bt_1d);
+   TensorValues<QVectorLayout::byVDIM>(ne, 1, dof2quad, x, y);
+}
+
+ChangeOfBasis_RT::ChangeOfBasis_RT(FiniteElementSpace &fes)
+   : Operator(fes.GetTrueVSize()),
+     fes(fes),
+     dim(fes.GetMesh()->Dimension()),
+     ne(fes.GetNE()),
+     p(fes.GetMaxElementOrder())
+{
+   auto op = fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+   elem_restr = dynamic_cast<const ElementRestriction*>(op);
+   MFEM_VERIFY(elem_restr != NULL, "Missing element restriciton.");
+
+   const auto *rt_fec = dynamic_cast<const RT_FECollection*>(fes.FEColl());
+   MFEM_VERIFY(rt_fec, "Must be RT finite element space.");
+
+   const int cb_type = rt_fec->GetClosedBasisType();
+   const int ob_type = rt_fec->GetOpenBasisType();
+
+   no_op = (cb_type == BasisType::GaussLobatto &&
+            ob_type == BasisType::IntegratedGLL);
+   if (no_op) { return; }
+
+   const int pp1 = p + 1;
+
+   Poly_1D::Basis &cbasis = poly1d.GetBasis(p, cb_type);
+   Poly_1D::Basis &obasis = poly1d.GetBasis(p-1, ob_type);
+
+   const double *cpts2 = poly1d.GetPoints(p, BasisType::GaussLobatto);
+
+   Bci_1d.SetSize(pp1*pp1);
+   Vector b(pp1);
+   for (int i = 0; i < pp1; ++i)
+   {
+      cbasis.Eval(cpts2[i], b);
+      for (int j = 0; j < pp1; ++j)
+      {
+         Bci_1d[i + j*pp1] = b[j];
+      }
+   }
+   SubcellIntegrals(p, obasis, Boi_1d);
+
+   ComputeInverse(Boi_1d, Bo_1d);
+   Transpose(Bo_1d, Bot_1d);
+   ComputeInverse(Bci_1d, Bc_1d);
+   Transpose(Bc_1d, Bct_1d);
+}
+
+const double *ChangeOfBasis_RT::GetOpenMap(Mode mode) const
+{
+   switch (mode)
+   {
+      case NORMAL: return Bo_1d.Read();
+      case TRANSPOSE: return Bot_1d.Read();
+      case INVERSE: return Boi_1d.Read();
+   }
+   return nullptr;
+}
+
+const double *ChangeOfBasis_RT::GetClosedMap(Mode mode) const
+{
+   switch (mode)
+   {
+      case NORMAL: return Bc_1d.Read();
+      case TRANSPOSE: return Bct_1d.Read();
+      case INVERSE: return Bci_1d.Read();
+   }
+   return nullptr;
+}
+
+void ChangeOfBasis_RT::MultRT_2D(const Vector &x, Vector &y, Mode mode) const
+{
+   const int DIM = dim;
+   const int NE = ne;
+   const int D1D = p + 1;
+   const int ND = (p+1)*p;
+   const double *BC = GetClosedMap(mode);
+   const double *BO = GetOpenMap(mode);
+   const auto X = Reshape(x.Read(), DIM*ND, ne);
+   auto Y = Reshape(y.Write(), DIM*ND, ne);
+
+   MFEM_FORALL(e, NE,
+   {
+      for (int c = 0; c < DIM; ++c)
+      {
+         const int nx = (c == 0) ? D1D : D1D-1;
+         const int ny = (c == 1) ? D1D : D1D-1;
+         const double *Bx = (c == 0) ? BC : BO;
+         const double *By = (c == 1) ? BC : BO;
+
+         for (int i = 0; i < ND; ++i)
+         {
+            Y(i + c*ND, e) = 0.0;
+         }
+         for (int iy = 0; iy < ny; ++ iy)
+         {
+            double xx[MAX_D1D];
+            for (int ix = 0; ix < nx; ++ix) { xx[ix] = 0.0; }
+            for (int jx = 0; jx < nx; ++jx)
+            {
+               const double val = X(jx + iy*nx + c*nx*ny, e);
+               for (int ix = 0; ix < nx; ++ix)
+               {
+                  xx[ix] += val*Bx[ix + jx*nx];
+               }
+            }
+            for (int jy = 0; jy < ny; ++jy)
+            {
+               const double b = By[jy + iy*ny];
+               for (int ix = 0; ix < nx; ++ix)
+               {
+                  Y(ix + jy*nx + c*nx*ny, e) += xx[ix]*b;
+               }
+            }
+         }
+      }
+   });
+}
+
+void ChangeOfBasis_RT::MultRT_3D(const Vector &x, Vector &y, Mode mode) const
+{
+   const int DIM = dim;
+   const int NE = ne;
+   const int D1D = p + 1;
+   const int ND = (p+1)*p*p;
+   const double *BC = GetClosedMap(mode);
+   const double *BO = GetOpenMap(mode);
+   const auto X = Reshape(x.Read(), DIM*ND, ne);
+   auto Y = Reshape(y.Write(), DIM*ND, ne);
+
+   MFEM_FORALL(e, NE,
+   {
+      for (int c = 0; c < DIM; ++c)
+      {
+         const int nx = (c == 0) ? D1D : D1D-1;
+         const int ny = (c == 1) ? D1D : D1D-1;
+         const int nz = (c == 2) ? D1D : D1D-1;
+         const double *Bx = (c == 0) ? BC : BO;
+         const double *By = (c == 1) ? BC : BO;
+         const double *Bz = (c == 2) ? BC : BO;
+
+         for (int i = 0; i < ND; ++i)
+         {
+            Y(i + c*ND, e) = 0.0;
+         }
+         for (int iz = 0; iz < nz; ++ iz)
+         {
+            double xy[MAX_D1D][MAX_D1D];
+            for (int iy = 0; iy < ny; ++iy)
+            {
+               for (int ix = 0; ix < nx; ++ix)
+               {
+                  xy[iy][ix] = 0.0;
+               }
+            }
+            for (int iy = 0; iy < ny; ++iy)
+            {
+               double xx[MAX_D1D];
+               for (int ix = 0; ix < nx; ++ix) { xx[ix] = 0.0; }
+               for (int ix = 0; ix < nx; ++ix)
+               {
+                  const double val = X(ix + iy*nx + iz*nx*ny + c*ND, e);
+                  for (int jx = 0; jx < nx; ++jx)
+                  {
+                     xx[jx] += val*Bx[jx + ix*nx];
+                  }
+               }
+               for (int jy = 0; jy < ny; ++jy)
+               {
+                  const double b = By[jy + iy*ny];
+                  for (int jx = 0; jx < nx; ++jx)
+                  {
+                     xy[jy][jx] += xx[jx] * b;
+                  }
+               }
+            }
+            for (int jz = 0; jz < nz; ++jz)
+            {
+               const double b = Bz[jz + iz*nz];
+               for (int jy = 0; jy < ny; ++jy)
+               {
+                  for (int jx = 0; jx < nx; ++jx)
+                  {
+                     Y(jx + jy*nx + jz*nx*ny + c*ND, e) += xy[jy][jx] * b;
+                  }
+               }
+            }
+         }
+      }
+   });
+}
+
+void ChangeOfBasis_RT::Mult(const Vector &x, Vector &y, Mode mode) const
+{
+   if (no_op) { y = x; return; }
+
+   const Operator *P = fes.GetProlongationMatrix();
+
+   if (IsIdentityProlongation(P))
+   {
+      x_l.MakeRef(const_cast<Vector&>(x), 0, fes.GetVSize());
+      y_l.MakeRef(y, 0, fes.GetVSize());
+   }
+   else
+   {
+      x_l.SetSize(fes.GetVSize());
+      y_l.SetSize(fes.GetVSize());
+      P->Mult(x, x_l);
+   }
+
+   x_e.SetSize(elem_restr->Height());
+   y_e.SetSize(elem_restr->Height());
+
+   elem_restr->Mult(x_l, x_e);
+
+   if (dim == 2) { MultRT_2D(x_e, y_e, mode); }
+   else { MultRT_3D(x_e, y_e, mode); }
+
+   elem_restr->MultLeftInverse(y_e, y_l);
+
+   const Operator *R = fes.GetRestrictionOperator();
+   if (R) { R->Mult(y_l, y); }
+   else { MFEM_VERIFY(P == NULL, "Invalid state."); }
+}
+
+void ChangeOfBasis_RT::Mult(const Vector &x, Vector &y) const
+{
+   Mult(x, y, NORMAL);
+}
+
+void ChangeOfBasis_RT::MultTranspose(const Vector &x, Vector &y) const
+{
+   Mult(x, y, TRANSPOSE);
+}
+
+void ChangeOfBasis_RT::MultInverse(const Vector &x, Vector &y) const
+{
+   Mult(x, y, INVERSE);
+}
+
+} // namespace mfem

miniapps/hdiv-linear-solver/change_basis.hpp

@@ -0,0 +1,87 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef CHANGE_BASIS_HPP
+#define CHANGE_BASIS_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+/// @brief Change of basis operator between L2 spaces.
+///
+/// This represents the change-of-basis operator from the given L2 space to a
+/// space using the IntegratedGLL basis.
+class ChangeOfBasis_L2 : public Operator
+{
+private:
+   const int ne; ///< Number of elements in the mesh.
+   mutable DofToQuad dof2quad; ///< 1D basis transformation.
+   Array<double> B_1d; ///< 1D basis transformation matrix.
+   Array<double> Bt_1d; ///< 1D basis transformation matrix traspose.
+   bool no_op; ///< If the basis types are the same, the operation is a no-op.
+public:
+   ChangeOfBasis_L2(FiniteElementSpace &fes);
+   void Mult(const Vector &x, Vector &y) const override;
+   void MultTranspose(const Vector &x, Vector &y) const override;
+};
+
+/// Change of basis operator between RT spaces.
+///
+/// This represents the change-of-basis operator from the given RT space to a
+/// space using Gauss-Lobatto as the "open" basis and IntegratedGLL as the
+/// "closed" basis.
+class ChangeOfBasis_RT : public Operator
+{
+public:
+   // Should be private, nvcc limitation...
+   enum Mode
+   {
+      NORMAL,
+      TRANSPOSE,
+      INVERSE
+   };
+private:
+   FiniteElementSpace &fes; ///< The finite element space.
+   const int dim; ///< Dimension of the mesh.
+   const int ne; ///< Number of elements.
+   const int p; ///< Polynomial degree.
+   const ElementRestriction *elem_restr; ///< Element restriction operator.
+   Array<double> Bc_1d; ///< 1D closed basis transformation matrix.
+   Array<double> Bci_1d; ///< 1D closed basis transformation matrix inverse.
+   Array<double> Bct_1d; ///< 1D closed basis transformation matrix transpose.
+   Array<double> Bo_1d; ///< 1D open basis transformation matrix.
+   Array<double> Boi_1d; ///< 1D open basis transformation matrix inverse.
+   Array<double> Bot_1d; ///< 1D open basis transformation matrix transpose.
+
+   mutable Vector x_l, y_l; ///< L-vector layout
+   mutable Vector x_e, y_e; ///< E-vector layout
+
+   bool no_op; ///< If the spaces are the same, the operation is a no-op.
+
+   void Mult(const Vector &x, Vector &y, Mode mode) const;
+   const double *GetOpenMap(Mode mode) const;
+   const double *GetClosedMap(Mode mode) const;
+public:
+   ChangeOfBasis_RT(FiniteElementSpace &fes);
+   void Mult(const Vector &x, Vector &y) const override;
+   void MultTranspose(const Vector &x, Vector &y) const override;
+   void MultInverse(const Vector &x, Vector &y) const;
+   // The following should be considered private, public because of compiler
+   // limitations
+   void MultRT_2D(const Vector &x, Vector &y, Mode mode) const;
+   void MultRT_3D(const Vector &x, Vector &y, Mode mode) const;
+};
+
+} // namespace mfem
+
+#endif

miniapps/hdiv-linear-solver/darcy.cpp

@@ -0,0 +1,212 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+//                     ---------------------------------
+//                     Poisson/Darcy Mixed Method Solver
+//                     ---------------------------------
+//
+// Solves a Poisson problem -Delta p = f using a mixed finite element
+// formulation). The right-hand side of the Poisson problem is the same as that
+// used in the LOR Solvers miniapp (see miniapps/solvers). Dirichlet boundary
+// conditions are enforced on all domain boundaries.
+//
+// Optionally, the equation alpha*p - Delta p = f can be solved by setting the
+// alpha parameter to a nonzero value.
+
+// This can be written in the form of a Darcy problem
+//
+//                           -u - grad(p) = 0
+//                       alpha*p + div(u) = f
+//
+// where natural boundary conditions are enforced on the flux u, and the
+// Dirichlet condition on p is enforced by modifying the right-hand side.
+//
+// The resulting saddle-point system is solved using MINRES with a matrix-free
+// block-diagonal preconditioner.
+//
+// See also example 5 and its parallel version.
+//
+// Sample runs:
+//
+//    darcy
+//    mpirun -np 4 darcy -m ../../data/fichera-q2.mesh
+
+#include "mfem.hpp"
+#include <iostream>
+#include <memory>
+
+#include "discrete_divergence.hpp"
+#include "hdiv_linear_solver.hpp"
+
+#include "../solvers/lor_mms.hpp"
+
+using namespace std;
+using namespace mfem;
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
+
+class RobinCoefficient : public Coefficient
+{
+   double Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      double xdata[3];
+      Vector xvec(xdata, 3);
+      T.Transform(ip, xvec);
+      const int dim = xvec.Size();
+
+      Vector n(dim);
+      CalcOrtho(T.Jacobian(), n);
+      n /= n.Norml2();
+
+      const double p_val = u(xvec);
+      const double x = pi*xvec[0];
+      const double y = pi*xvec[1];
+
+      if (dim == 2)
+      {
+         const double u_val = -pi*(n[0]*cos(x)*sin(y) + n[1]*sin(x)*cos(y));
+         return p_val - u_val;
+      }
+      else
+      {
+         MFEM_ABORT("Not implemented");
+      }
+      return 0.0;
+   }
+};
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *mesh_file = "../../data/star.mesh";
+   const char *device_config = "cpu";
+   int ser_ref = 1;
+   int par_ref = 1;
+   int order = 3;
+   double alpha = 0.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+   args.AddOption(&par_ref, "-rp", "--parallel-refine",
+                  "Number of times to refine the mesh in parallel.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&alpha, "-a", "--alpha", "Value of alpha coefficient.");
+   args.ParseCheck();
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
+   const int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+
+   const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
+   const int mt = FiniteElement::VALUE;
+   RT_FECollection fec_rt(order-1, dim, b1, b2);
+   L2_FECollection fec_l2(order-1, dim, b2, mt);
+   ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
+   ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
+
+   HYPRE_BigInt ndofs_rt = fes_rt.GlobalTrueVSize();
+   HYPRE_BigInt ndofs_l2 = fes_l2.GlobalTrueVSize();
+
+   if (Mpi::Root())
+   {
+      cout << "\nRT DOFs: " << ndofs_rt << "\nL2 DOFs: " << ndofs_l2 << endl;
+   }
+
+   Array<int> ess_rt_dofs; // empty
+
+   // f is the RHS, u is the exact solution
+   FunctionCoefficient f_coeff(f(alpha)), u_coeff(u);
+   // Coefficient to enforce Robin boundary condition
+   RobinCoefficient bc_coeff;
+
+   // Assemble the right-hand side for the scalar (L2) unknown.
+   ParLinearForm b_l2(&fes_l2);
+   b_l2.AddDomainIntegrator(new DomainLFIntegrator(f_coeff));
+   b_l2.UseFastAssembly(true);
+   b_l2.Assemble();
+
+   // Enforce Dirichlet boundary conditions on the scalar unknown by adding
+   // the boundary term to the flux equation.
+   ParLinearForm b_rt(&fes_rt);
+   b_rt.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(bc_coeff));
+   b_rt.UseFastAssembly(true);
+   b_rt.Assemble();
+
+   if (Mpi::Root()) { cout << "\nSaddle point solver... " << flush; }
+   tic_toc.Clear(); tic_toc.Start();
+
+   // Set up the block system of the form
+   //
+   // [  W     D ][ u ] = [  f  ]
+   // [ D^T   -M ][ q ] = [ g_D ]
+   //
+   // where W is the L2 mass matrix, D is the discrete divergence, and M is
+   // the RT mass matrix.
+   //
+   // If the coefficient alpha is set to zero, the system takes the form
+   //
+   // [  0     D ][ u ] = [  f  ]
+   // [ D^T   -M ][ q ] = [ g_D ]
+   //
+   // u is the scalar unknown, and q is the flux. f is the right-hand side from
+   // the Poisson problem, and g_D is the contribution to the right-hand side
+   // from the Dirichlet boundary condition.
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient alpha_coeff(alpha);
+   const auto solver_mode = HdivSaddlePointSolver::Mode::DARCY;
+   HdivSaddlePointSolver saddle_point_solver(
+      mesh, fes_rt, fes_l2, alpha_coeff, one, one, ess_rt_dofs, solver_mode);
+
+   const Array<int> &offsets = saddle_point_solver.GetOffsets();
+   BlockVector X_block(offsets), B_block(offsets);
+
+   b_l2.ParallelAssemble(B_block.GetBlock(0));
+   b_rt.ParallelAssemble(B_block.GetBlock(1));
+   B_block.SyncFromBlocks();
+
+   X_block = 0.0;
+   saddle_point_solver.Mult(B_block, X_block);
+   X_block.SyncToBlocks();
+
+   if (Mpi::Root())
+   {
+      cout << "Done.\nIterations: "
+           << saddle_point_solver.GetNumIterations()
+           << "\nElapsed: " << tic_toc.RealTime() << endl;
+   }
+
+   ParGridFunction x(&fes_l2);
+   x.SetFromTrueDofs(X_block.GetBlock(0));
+   const double error = x.ComputeL2Error(u_coeff);
+   if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+
+   return 0;
+}
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
+{
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
+   for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
+   return mesh;
+}

miniapps/hdiv-linear-solver/discrete_divergence.cpp

@@ -0,0 +1,276 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "general/forall.hpp"
+#include "discrete_divergence.hpp"
+
+namespace mfem
+{
+
+/// @brief Eliminates columns in the given HypreParMatrix.
+///
+/// This is similar to HypreParMatrix::EliminateBC, except that only the columns
+/// are eliminated.
+void EliminateColumns(HypreParMatrix &D, const Array<int> &ess_dofs)
+{
+
+   hypre_ParCSRMatrix *A_hypre = D;
+   D.HypreReadWrite();
+
+   hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A_hypre);
+   hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A_hypre);
+
+   HYPRE_Int diag_ncols = hypre_CSRMatrixNumCols(diag);
+   HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols(offd);
+
+   const int n_ess_dofs = ess_dofs.Size();
+
+   // Start communication to figure out which columns need to be eliminated in
+   // the off-diagonal block
+   hypre_ParCSRCommHandle *comm_handle;
+   HYPRE_Int *int_buf_data, *eliminate_col_diag, *eliminate_col_offd;
+   {
+      eliminate_col_diag = mfem_hypre_CTAlloc_host(HYPRE_Int, diag_ncols);
+      eliminate_col_offd = mfem_hypre_CTAlloc_host(HYPRE_Int, offd_ncols);
+
+      // Make sure A has a communication package
+      hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
+      if (!comm_pkg)
+      {
+         hypre_MatvecCommPkgCreate(A_hypre);
+         comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
+      }
+
+      // Which of the local columns are to be eliminated?
+      for (int i = 0; i < diag_ncols; i++)
+      {
+         eliminate_col_diag[i] = 0;
+      }
+
+      ess_dofs.HostRead();
+      for (int i = 0; i < n_ess_dofs; i++)
+      {
+         eliminate_col_diag[ess_dofs[i]] = 1;
+      }
+
+      // Use a matvec communication pattern to find (in eliminate_col_offd)
+      // which of the local offd columns are to be eliminated
+      HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
+      HYPRE_Int int_buf_sz = hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends);
+      int_buf_data = mfem_hypre_CTAlloc_host(HYPRE_Int, int_buf_sz);
+      HYPRE_Int *send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
+      for (int i = 0; i < int_buf_sz; ++i)
+      {
+         const int k = send_map_elmts[i];
+         int_buf_data[i] = eliminate_col_diag[k];
+      }
+      comm_handle = hypre_ParCSRCommHandleCreate(
+                       11, comm_pkg, int_buf_data, eliminate_col_offd);
+   }
+
+   // Eliminate columns in the diagonal block
+   {
+      Memory<HYPRE_Int> col_mem(eliminate_col_diag, diag_ncols, false);
+      const auto cols = col_mem.Read(GetHypreMemoryClass(), diag_ncols);
+      const int nrows_diag = hypre_CSRMatrixNumRows(diag);
+      const auto I = diag->i;
+      const auto J = diag->j;
+      auto data = diag->data;
+      MFEM_HYPRE_FORALL(i, nrows_diag,
+      {
+         for (int jj=I[i]; jj<I[i+1]; ++jj)
+         {
+            const int j = J[jj];
+            data[jj] *= 1 - cols[j];
+         }
+      });
+      col_mem.Delete();
+   }
+
+   // Wait for MPI communication to finish
+   hypre_ParCSRCommHandleDestroy(comm_handle);
+   mfem_hypre_TFree_host(int_buf_data);
+   mfem_hypre_TFree_host(eliminate_col_diag);
+
+   // Eliminate columns in the off-diagonal block
+   {
+      Memory<HYPRE_Int> col_mem(eliminate_col_offd, offd_ncols, false);
+      const auto cols = col_mem.Read(GetHypreMemoryClass(), offd_ncols);
+      const int nrows_offd = hypre_CSRMatrixNumRows(offd);
+      const auto I = offd->i;
+      const auto J = offd->j;
+      auto data = offd->data;
+      MFEM_HYPRE_FORALL(i, nrows_offd,
+      {
+         for (int jj=I[i]; jj<I[i+1]; ++jj)
+         {
+            const int j = J[jj];
+            data[jj] *= 1 - cols[j];
+         }
+      });
+      col_mem.Delete();
+   }
+
+   mfem_hypre_TFree_host(eliminate_col_offd);
+}
+
+void FormElementToFace2D(int order, Array<int> &element2face)
+{
+   const int o = order;
+   const int op1 = order + 1;
+
+   for (int iy = 0; iy < o; ++iy)
+   {
+      for (int ix = 0; ix < o; ++ix)
+      {
+         const int ivol = ix + iy*o;
+         element2face[0 + 4*ivol] = -1 - (ix + iy*op1); // left, x = 0
+         element2face[1 + 4*ivol] = ix+1 + iy*op1; // right, x = 1
+         element2face[2 + 4*ivol] = -1 - (ix + iy*o + o*op1); // bottom, y = 0
+         element2face[3 + 4*ivol] = ix + (iy+1)*o + o*op1; // top, y = 1
+      }
+   }
+}
+
+void FormElementToFace3D(int order, Array<int> &element2face)
+{
+   const int o = order;
+   const int op1 = order + 1;
+
+   const int n = o*o*op1; // number of faces per dimension
+
+   for (int iz = 0; iz < o; ++iz)
+   {
+      for (int iy = 0; iy < o; ++iy)
+      {
+         for (int ix = 0; ix < o; ++ix)
+         {
+            const int ivol = ix + iy*o + iz*o*o;
+            element2face[0 + 6*ivol] = -1 - (ix + iy*op1 + iz*o*op1); // x = 0
+            element2face[1 + 6*ivol] = ix+1 + iy*op1 + iz*o*op1; // x = 1
+            element2face[2 + 6*ivol] = -1 - (ix + iy*o + iz*o*op1 + n); // y = 0
+            element2face[3 + 6*ivol] = ix + (iy+1)*o + iz*o*op1 + n; // y = 1
+            element2face[4 + 6*ivol] = -1 - (ix + iy*o + iz*o*o + 2*n); // z = 0
+            element2face[5 + 6*ivol] = ix + iy*o + (iz+1)*o*o + 2*n; // z = 1
+         }
+      }
+   }
+}
+
+HypreParMatrix *FormDiscreteDivergenceMatrix(ParFiniteElementSpace &fes_rt,
+                                             ParFiniteElementSpace &fes_l2,
+                                             const Array<int> &ess_dofs)
+{
+   const Mesh &mesh = *fes_rt.GetMesh();
+   const int dim = mesh.Dimension();
+   const int order = fes_rt.GetMaxElementOrder();
+
+   const int n_rt = fes_rt.GetNDofs();
+   const int n_l2 = fes_l2.GetNDofs();
+
+   SparseMatrix D_local;
+   D_local.OverrideSize(n_l2, n_rt);
+
+   D_local.GetMemoryI().New(n_l2 + 1);
+   // Each row always has 2*dim nonzeros (one for each face of the element)
+   const int nnz = n_l2*2*dim;
+   auto I = D_local.WriteI();
+   MFEM_FORALL(i, n_l2+1, I[i] = 2*dim*i; );
+
+   const int nel_ho = mesh.GetNE();
+   const int nface_per_el = dim*pow(order, dim-1)*(order+1);
+   const int nvol_per_el = pow(order, dim);
+
+   // element2face is a mapping of size (2*dim, nvol_per_el) such that with a
+   // macro element, subelement i (in lexicographic ordering) has faces (also
+   // in lexicographic order) given by the entries (j, i).
+   Array<int> element2face;
+   element2face.SetSize(2*dim*nvol_per_el);
+
+   if (dim == 2) { FormElementToFace2D(order, element2face); }
+   else if (dim == 3) { FormElementToFace3D(order, element2face); }
+   else { MFEM_ABORT("Unsupported dimension.") }
+
+   const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
+   const auto *R_rt = dynamic_cast<const ElementRestriction*>(
+                         fes_rt.GetElementRestriction(ordering));
+   const auto gather_rt = Reshape(R_rt->GatherMap().Read(), nface_per_el, nel_ho);
+   const auto e2f = Reshape(element2face.Read(), 2*dim, nvol_per_el);
+
+   // Fill J and data
+   D_local.GetMemoryJ().New(nnz);
+   D_local.GetMemoryData().New(nnz);
+
+   auto J = D_local.WriteJ();
+   auto V = D_local.WriteData();
+
+   // Loop over L2 DOFs
+   MFEM_FORALL(i, n_l2,
+   {
+      const int i_loc = i%nvol_per_el;
+      const int i_el = i/nvol_per_el;
+
+      for (int k = 0; k < 2*dim; ++k)
+      {
+         const int sjv_loc = e2f(k, i_loc);
+         const int jv_loc = (sjv_loc >= 0) ? sjv_loc : -1 - sjv_loc;
+         const int sgn1 = (sjv_loc >= 0) ? 1 : -1;
+         const int sj = gather_rt(jv_loc, i_el);
+         const int j = (sj >= 0) ? sj : -1 - sj;
+         const int sgn2 = (sj >= 0) ? 1 : -1;
+
+         J[k + 2*dim*i] = j;
+         V[k + 2*dim*i] = sgn1*sgn2;
+      }
+   });
+
+   // Create a block diagonal parallel matrix
+   OperatorHandle D_diag(Operator::Hypre_ParCSR);
+   D_diag.MakeRectangularBlockDiag(fes_rt.GetComm(),
+                                   fes_l2.GlobalVSize(),
+                                   fes_rt.GlobalVSize(),
+                                   fes_l2.GetDofOffsets(),
+                                   fes_rt.GetDofOffsets(),
+                                   &D_local);
+
+   HypreParMatrix *D;
+   // Assemble the parallel gradient matrix, must be deleted by the caller
+   if (IsIdentityProlongation(fes_rt.GetProlongationMatrix()))
+   {
+      D = D_diag.As<HypreParMatrix>();
+      D_diag.SetOperatorOwner(false);
+      HypreStealOwnership(*D, D_local);
+   }
+   else
+   {
+      OperatorHandle Rt(Transpose(*fes_l2.GetRestrictionMatrix()));
+      OperatorHandle Rt_diag(Operator::Hypre_ParCSR);
+      Rt_diag.MakeRectangularBlockDiag(fes_l2.GetComm(),
+                                       fes_l2.GlobalVSize(),
+                                       fes_l2.GlobalTrueVSize(),
+                                       fes_l2.GetDofOffsets(),
+                                       fes_l2.GetTrueDofOffsets(),
+                                       Rt.As<SparseMatrix>());
+      D = RAP(Rt_diag.As<HypreParMatrix>(),
+              D_diag.As<HypreParMatrix>(),
+              fes_rt.Dof_TrueDof_Matrix());
+   }
+   D->CopyRowStarts();
+   D->CopyColStarts();
+
+   // Eliminate the boundary conditions
+   EliminateColumns(*D, ess_dofs);
+
+   return D;
+}
+
+} // namespace mfem

miniapps/hdiv-linear-solver/discrete_divergence.hpp

@@ -0,0 +1,32 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_DISCRETE_DIVERGENCE_HPP
+#define MFEM_DISCRETE_DIVERGENCE_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+/// @brief Eliminates columns in the given HypreParMatrix.
+///
+/// This is similar to HypreParMatrix::EliminateBC, except that only the columns
+/// are eliminated.
+void EliminateColumns(HypreParMatrix &D, const Array<int> &ess_dofs);
+
+HypreParMatrix *FormDiscreteDivergenceMatrix(ParFiniteElementSpace &fes_rt,
+                                             ParFiniteElementSpace &fes_l2,
+                                             const Array<int> &ess_dofs);
+
+} // namespace mfem
+
+#endif

miniapps/hdiv-linear-solver/grad_div.cpp

@@ -0,0 +1,280 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+//                     ---------------------------------
+//                     H(div) saddle-point system solver
+//                     ---------------------------------
+//
+// Solves the grad-div problem u - grad(div(u)) = f using a variety of solver
+// techniques. This miniapp supports solving this problem using a variety of
+// matrix-free and matrix-based preconditioning methods, inclding:
+//
+// * Matrix-free block-diagonal preconditioning for the saddle-point system.
+// * ADS-AMG preconditioning.
+// * Low-order-refined ADS-AMG preconditioning (matrix-free).
+// * Hybridization with AMG preconditioning.
+//
+// The problem setup is the same as in the LOR solvers miniapps (in the
+// miniapps/solvers directory). Dirichlet conditions are enforced on the normal
+// component of u.
+//
+// Sample runs:
+//
+//    grad_div -sp -ams -lor -hb
+//    mpirun -np 4 grad_div -sp -ams -lor -hb -m ../../data/fichera-q2.mesh -rp 0
+
+#include "mfem.hpp"
+#include <iostream>
+#include <memory>
+#include "hdiv_linear_solver.hpp"
+#include "../solvers/lor_mms.hpp"
+
+using namespace std;
+using namespace mfem;
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
+void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X);
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *mesh_file = "../../data/star.mesh";
+   const char *device_config = "cpu";
+   int ser_ref = 1;
+   int par_ref = 1;
+   int order = 3;
+   bool use_saddle_point = false;
+   bool use_ams = false;
+   bool use_lor_ams = false;
+   bool use_hybridization = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+   args.AddOption(&par_ref, "-rp", "--parallel-refine",
+                  "Number of times to refine the mesh in parallel.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&use_saddle_point,
+                  "-sp", "--saddle-point", "-no-sp", "--no-saddle-point",
+                  "Enable or disable saddle-point solver.");
+   args.AddOption(&use_ams, "-ams", "--ams", "-no-ams", "--no-ams",
+                  "Enable or disable AMS solver.");
+   args.AddOption(&use_lor_ams, "-lor", "--lor-ams", "-no-lor", "--no-lor-ams",
+                  "Enable or disable LOR-AMS solver.");
+   args.AddOption(&use_hybridization,
+                  "-hb", "--hybridization", "-no-hb", "--no-hybridization",
+                  "Enable or disable hybridization solver.");
+   args.ParseCheck();
+
+   if (!use_saddle_point && !use_ams && !use_lor_ams && !use_hybridization)
+   {
+      if (Mpi::Root()) { cout << "No solver enabled. Exiting.\n"; }
+      return 0;
+   }
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
+   const int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+
+   const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
+   RT_FECollection fec_rt(order-1, dim, b1, b2);
+   ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
+
+   Array<int> ess_rt_dofs;
+   fes_rt.GetBoundaryTrueDofs(ess_rt_dofs);
+
+   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(true)), u_vec_coeff(dim, u_vec);
+
+   ParLinearForm b(&fes_rt);
+   b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff));
+   b.UseFastAssembly(true);
+   b.Assemble();
+
+   ConstantCoefficient alpha_coeff(1.0);
+   ConstantCoefficient beta_coeff(1.0);
+
+   ParGridFunction x(&fes_rt);
+   x.ProjectCoefficient(u_vec_coeff);
+
+   cout.precision(4);
+   cout << scientific;
+
+   if (use_saddle_point)
+   {
+      if (Mpi::Root()) { cout << "\nSaddle point solver... " << flush; }
+      tic_toc.Clear(); tic_toc.Start();
+
+      const int mt = FiniteElement::INTEGRAL;
+      L2_FECollection fec_l2(order-1, dim, b2, mt);
+      ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
+
+      HdivSaddlePointSolver saddle_point_solver(
+         mesh, fes_rt, fes_l2, alpha_coeff, beta_coeff, ess_rt_dofs,
+         HdivSaddlePointSolver::Mode::GRAD_DIV);
+
+      const Array<int> &offsets = saddle_point_solver.GetOffsets();
+
+      BlockVector X_block(offsets), B_block(offsets);
+      B_block.GetBlock(0) = 0.0;
+      b.ParallelAssemble(B_block.GetBlock(1));
+      B_block.GetBlock(1) *= -1.0;
+      B_block.SyncFromBlocks();
+
+      x.ParallelProject(X_block.GetBlock(1));
+      saddle_point_solver.SetBC(X_block.GetBlock(1));
+
+      X_block = 0.0;
+      saddle_point_solver.Mult(B_block, X_block);
+
+      if (Mpi::Root())
+      {
+         cout << "Done.\nIterations: "
+              << saddle_point_solver.GetNumIterations()
+              << "\nElapsed: " << tic_toc.RealTime() << endl;
+      }
+
+      X_block.SyncToBlocks();
+      x.SetFromTrueDofs(X_block.GetBlock(1));
+      const double error = x.ComputeL2Error(u_vec_coeff);
+      if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+   }
+
+   if (use_ams)
+   {
+      if (Mpi::Root()) { cout << "\nAMS solver... " << flush; }
+      tic_toc.Clear(); tic_toc.Start();
+
+      ParBilinearForm a(&fes_rt);
+      a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
+      a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
+      a.Assemble();
+
+      OperatorHandle A;
+      Vector B, X;
+      b.Assemble();
+      x.ProjectCoefficient(u_vec_coeff);
+      a.FormLinearSystem(ess_rt_dofs, x, b, A, X, B);
+      HypreParMatrix &Ah = *A.As<HypreParMatrix>();
+
+      std::unique_ptr<Solver> prec;
+      if (dim == 2) { prec.reset(new HypreAMS(Ah, &fes_rt)); }
+      else  { prec.reset(new HypreADS(Ah, &fes_rt)); }
+
+      SolveCG(Ah, *prec, B, X);
+      x.SetFromTrueDofs(X);
+      const double error = x.ComputeL2Error(u_vec_coeff);
+      if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+   }
+
+   if (use_lor_ams)
+   {
+      const int b2_lor = BasisType::IntegratedGLL;
+      RT_FECollection fec_rt_lor(order-1, dim, b1, b2_lor);
+      ParFiniteElementSpace fes_rt_lor(&mesh, &fec_rt_lor);
+
+      ParLinearForm b_lor(&fes_rt_lor);
+      b_lor.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff));
+      b_lor.UseFastAssembly(true);
+      b_lor.Assemble();
+
+      if (Mpi::Root()) { cout << "\nLOR-AMS solver... " << flush; }
+      tic_toc.Clear(); tic_toc.Start();
+
+      ParBilinearForm a(&fes_rt_lor);
+      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
+      a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
+      a.Assemble();
+
+      ParGridFunction x_lor(&fes_rt_lor);
+      x_lor.ProjectCoefficient(u_vec_coeff);
+
+      OperatorHandle A;
+      Vector B, X;
+      a.FormLinearSystem(ess_rt_dofs, x_lor, b_lor, A, X, B);
+
+      std::unique_ptr<Solver> prec;
+      if (dim == 2) { prec.reset(new LORSolver<HypreAMS>(a, ess_rt_dofs)); }
+      else { prec.reset(new LORSolver<HypreADS>(a, ess_rt_dofs)); }
+
+      SolveCG(*A, *prec, B, X);
+      a.RecoverFEMSolution(X, b_lor, x_lor);
+      const double error = x_lor.ComputeL2Error(u_vec_coeff);
+      if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+   }
+
+   if (use_hybridization)
+   {
+      if (Mpi::Root()) { cout << "\nHybridization solver... " << flush; }
+      tic_toc.Clear(); tic_toc.Start();
+
+      DG_Interface_FECollection fec_hb(order-1, dim);
+      ParFiniteElementSpace fes_hb(&mesh, &fec_hb);
+
+      ParBilinearForm a(&fes_rt);
+      a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
+      a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
+      a.EnableHybridization(&fes_hb, new NormalTraceJumpIntegrator, ess_rt_dofs);
+      a.Assemble();
+
+      OperatorHandle A;
+      Vector B, X;
+      b.Assemble();
+      x.ProjectCoefficient(u_vec_coeff);
+      a.FormLinearSystem(ess_rt_dofs, x, b, A, X, B);
+
+      HypreBoomerAMG amg_hb(*A.As<HypreParMatrix>());
+      amg_hb.SetPrintLevel(0);
+
+      SolveCG(*A, amg_hb, B, X);
+      a.RecoverFEMSolution(X, b, x);
+      const double error = x.ComputeL2Error(u_vec_coeff);
+      if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+   }
+
+   return 0;
+}
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
+{
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
+   for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
+   return mesh;
+}
+
+void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X)
+{
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetAbsTol(0.0);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(500);
+   cg.SetPrintLevel(0);
+   cg.SetOperator(A);
+   cg.SetPreconditioner(P);
+   X = 0.0;
+   cg.Mult(B, X);
+   if (Mpi::Root())
+   {
+      cout << "Done.\nIterations: " << cg.GetNumIterations()
+            << "\nElapsed: " << tic_toc.RealTime() << endl;
+   }
+};

miniapps/hdiv-linear-solver/hdiv_linear_solver.cpp

@@ -0,0 +1,387 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "general/forall.hpp"
+#include "hdiv_linear_solver.hpp"
+#include "discrete_divergence.hpp"
+
+namespace mfem
+{
+
+/// Replace x[i] with 1.0/x[i] for all i.
+void Reciprocal(Vector &x)
+{
+   const int n = x.Size();
+   double *d_x = x.ReadWrite();
+   MFEM_FORALL(i, n, d_x[i] = 1.0/d_x[i]; );
+}
+
+/// Return a new HypreParMatrix with given diagonal entries
+HypreParMatrix *MakeDiagonalMatrix(Vector &diag,
+                                   const ParFiniteElementSpace &fes)
+{
+   const int n = diag.Size();
+
+   SparseMatrix diag_spmat;
+   diag_spmat.OverrideSize(n, n);
+   diag_spmat.GetMemoryI().New(n+1, Device::GetDeviceMemoryType());
+   diag_spmat.GetMemoryJ().New(n, Device::GetDeviceMemoryType());
+   diag_spmat.GetMemoryData().New(n, Device::GetDeviceMemoryType());
+
+   {
+      int *I = diag_spmat.WriteI();
+      int *J = diag_spmat.WriteJ();
+      double *A = diag_spmat.WriteData();
+      const double *d_diag = diag.Read();
+      MFEM_FORALL(i, n+1, I[i] = i;);
+      MFEM_FORALL(i, n,
+      {
+         J[i] = i;
+         A[i] = d_diag[i];
+      });
+   }
+
+   HYPRE_BigInt global_size = fes.GlobalTrueVSize();
+   HYPRE_BigInt *row_starts = fes.GetTrueDofOffsets();
+   HypreParMatrix D(MPI_COMM_WORLD, global_size, row_starts, &diag_spmat);
+   return new HypreParMatrix(D); // make a deep copy
+}
+
+const IntegrationRule &GetMassIntRule(FiniteElementSpace &fes_l2)
+{
+   Mesh *mesh = fes_l2.GetMesh();
+   const FiniteElement *fe = fes_l2.GetFE(0);
+   return MassIntegrator::GetRule(*fe, *fe, *mesh->GetElementTransformation(0));
+}
+
+HdivSaddlePointSolver::HdivSaddlePointSolver(
+   ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
+   Coefficient &L_coeff_, Coefficient &R_coeff_, Coefficient &B_coeff_,
+   const Array<int> &ess_rt_dofs_, Mode mode_)
+   : minres(mesh.GetComm()),
+     order(fes_rt_.GetMaxElementOrder()),
+     fec_l2(order - 1, mesh.Dimension(), b2, mt),
+     fes_l2(&mesh, &fec_l2),
+     fec_rt(order - 1, mesh.Dimension(), b1, b2),
+     fes_rt(&mesh, &fec_rt),
+     ess_rt_dofs(ess_rt_dofs_),
+     basis_l2(fes_l2_),
+     basis_rt(fes_rt_),
+     convert_map_type(fes_l2_.GetFE(0)->GetMapType() == FiniteElement::VALUE),
+     mass_l2(&fes_l2),
+     mass_rt(&fes_rt),
+     L_coeff(L_coeff_),
+     R_coeff(R_coeff_),
+     B_coeff(B_coeff_),
+     mode(mode_),
+     qs(mesh, GetMassIntRule(fes_l2)),
+     W_coeff_qf(qs),
+     W_mix_coeff_qf(qs),
+     W_coeff(W_coeff_qf),
+     W_mix_coeff(W_mix_coeff_qf)
+{
+   // If the user gives zero L coefficient, switch mode to DARCY_ZERO
+   auto *L_const_coeff = dynamic_cast<ConstantCoefficient*>(&L_coeff);
+   zero_l2_block = (L_const_coeff && L_const_coeff->constant == 0.0);
+
+   if (mode == Mode::GRAD_DIV)
+   {
+      MFEM_VERIFY(!zero_l2_block,
+                  "Mode::GRAD_DIV incompatible with zero coefficient.");
+   }
+
+   mass_l2.AddDomainIntegrator(new MassIntegrator(W_coeff));
+   mass_l2.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+   mass_rt.AddDomainIntegrator(new VectorFEMassIntegrator(&R_coeff));
+   mass_rt.AddBoundaryIntegrator(new MassIntegrator(B_coeff));
+   // mass_rt.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+   D.reset(FormDiscreteDivergenceMatrix(fes_rt, fes_l2, ess_rt_dofs));
+   Dt.reset(D->Transpose());
+
+   // Versions without BCs needed for elimination
+   D_e.reset(FormDiscreteDivergenceMatrix(fes_rt, fes_l2, empty));
+   mass_rt.FormSystemMatrix(empty, R_e);
+
+   offsets.SetSize(3);
+   offsets[0] = 0;
+   offsets[1] = fes_l2.GetTrueVSize();
+   offsets[2] = offsets[1] + fes_rt.GetTrueVSize();
+
+   minres.SetAbsTol(0.0);
+   minres.SetRelTol(1e-12);
+   minres.SetMaxIter(500);
+   minres.SetPrintLevel(IterativeSolver::PrintLevel().None());
+   minres.iterative_mode = false;
+
+   R_diag.SetSize(fes_rt.GetTrueVSize());
+   L_diag.SetSize(fes_l2.GetTrueVSize());
+
+   S_inv.SetPrintLevel(0);
+
+   if (mode == Mode::DARCY && !zero_l2_block)
+   {
+      ParBilinearForm mass_l2_unweighted(&fes_l2);
+      QuadratureFunction det_J_qf(qs);
+      QuadratureFunctionCoefficient det_J_coeff(det_J_qf);
+      if (convert_map_type)
+      {
+         const auto flags = GeometricFactors::DETERMINANTS;
+         auto *geom = fes_l2.GetMesh()->GetGeometricFactors(qs.GetIntRule(0), flags);
+         det_J_qf = geom->detJ;
+         mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator(det_J_coeff));
+      }
+      else
+      {
+         mass_l2_unweighted.AddDomainIntegrator(new MassIntegrator);
+      }
+      mass_l2_unweighted.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      mass_l2_unweighted.Assemble();
+      const int n_l2 = fes_l2.GetTrueVSize();
+      L_diag_unweighted.SetSize(n_l2);
+      mass_l2_unweighted.AssembleDiagonal(L_diag_unweighted);
+   }
+
+   Setup();
+}
+
+HdivSaddlePointSolver::HdivSaddlePointSolver(
+   ParMesh &mesh_, ParFiniteElementSpace &fes_rt_,
+   ParFiniteElementSpace &fes_l2_, Coefficient &L_coeff_, Coefficient &R_coeff_,
+   const Array<int> &ess_rt_dofs_, Mode mode_)
+   : HdivSaddlePointSolver(mesh_, fes_rt_, fes_l2_, L_coeff_, R_coeff_, zero,
+                           ess_rt_dofs_, mode_)
+{ }
+
+HdivSaddlePointSolver::HdivSaddlePointSolver(
+   ParMesh &mesh, ParFiniteElementSpace &fes_rt_, ParFiniteElementSpace &fes_l2_,
+   Coefficient &R_coeff_, const Array<int> &ess_rt_dofs_)
+   : HdivSaddlePointSolver(mesh, fes_rt_, fes_l2_, zero, R_coeff_, zero,
+                           ess_rt_dofs_, Mode::DARCY)
+{ }
+
+void HdivSaddlePointSolver::Setup()
+{
+   const auto flags = GeometricFactors::DETERMINANTS;
+   auto *geom = fes_l2.GetMesh()->GetGeometricFactors(qs.GetIntRule(0), flags);
+
+   if (!zero_l2_block) { L_coeff.Project(W_coeff_qf); }
+   // In "grad-div mode", the transformation matrix is scaled by the coefficient
+   // of the mass and divergence matrices.
+   // In "Darcy mode", the transformation matrix is unweighted.
+   if (mode == Mode::GRAD_DIV) { W_mix_coeff_qf = W_coeff_qf; }
+   else { W_mix_coeff_qf = 1.0; }
+
+   // The transformation matrix has to be "mixed" value and integral map type,
+   // which means that the coefficient has to be scaled like the Jacobian
+   // determinant.
+   if (convert_map_type)
+   {
+      const int n = W_mix_coeff_qf.Size();
+      const double *d_detJ = geom->detJ.Read();
+      double *d_w_mix = W_mix_coeff_qf.ReadWrite();
+      double *d_w = W_coeff_qf.ReadWrite();
+      const bool zero_l2 = zero_l2_block;
+      MFEM_FORALL(i, n,
+      {
+         const double detJ = d_detJ[i];
+         if (!zero_l2) { d_w[i] *= detJ*detJ; }
+         d_w_mix[i] *= detJ;
+      });
+   }
+
+   L_inv.reset(new DGMassInverse(fes_l2, W_mix_coeff));
+
+   if (zero_l2_block)
+   {
+      A_11.reset();
+   }
+   else
+   {
+      mass_l2.Assemble();
+      mass_l2.AssembleDiagonal(L_diag);
+      mass_l2.FormSystemMatrix(empty, L);
+
+      A_11.reset(new RAPOperator(*L_inv, *L, *L_inv));
+
+      if (mode == GRAD_DIV)
+      {
+         L_diag_unweighted.SetSize(L_diag.Size());
+
+         BilinearForm mass_l2_mix(&fes_l2);
+         mass_l2_mix.AddDomainIntegrator(new MassIntegrator(W_mix_coeff));
+         mass_l2_mix.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+         mass_l2_mix.Assemble();
+         mass_l2_mix.AssembleDiagonal(L_diag_unweighted);
+      }
+
+      const double *d_L_diag_unweighted = L_diag_unweighted.Read();
+      double *d_L_diag = L_diag.ReadWrite();
+      MFEM_FORALL(i, L_diag.Size(),
+      {
+         const double d = d_L_diag_unweighted[i];
+         d_L_diag[i] /= d*d;
+      });
+   }
+
+   // Reassmble the RT mass operator with the new coefficient
+   mass_rt.Update();
+   mass_rt.Assemble();
+   mass_rt.FormSystemMatrix(ess_rt_dofs, R);
+
+   // Form the updated approximate Schur complement
+   mass_rt.AssembleDiagonal(R_diag);
+
+   // Update the mass RT diagonal for essential DOFs
+   {
+      const int *d_I = ess_rt_dofs.Read();
+      double *d_R_diag = R_diag.ReadWrite();
+      MFEM_FORALL(i, ess_rt_dofs.Size(), d_R_diag[d_I[i]] = 1.0;);
+   }
+
+   // Form the approximate Schur complement
+   {
+      Reciprocal(R_diag);
+      std::unique_ptr<HypreParMatrix> R_diag_inv(MakeDiagonalMatrix(R_diag, fes_rt));
+      if (zero_l2_block)
+      {
+         S.reset(RAP(R_diag_inv.get(), Dt.get()));
+      }
+      else
+      {
+         std::unique_ptr<HypreParMatrix> D_Minv_Dt(RAP(R_diag_inv.get(), Dt.get()));
+         std::unique_ptr<HypreParMatrix> L_diag_inv(MakeDiagonalMatrix(L_diag, fes_l2));
+         S.reset(ParAdd(D_Minv_Dt.get(), L_diag_inv.get()));
+      }
+   }
+
+   // Reassemble the preconditioners
+   R_inv.reset(new OperatorJacobiSmoother(mass_rt, ess_rt_dofs));
+   S_inv.SetOperator(*S);
+
+   // Set up the block operators
+   A_block.reset(new BlockOperator(offsets));
+   // Omit the (1,1)-block when the L coefficient is identically zero.
+   if (A_11) { A_block->SetBlock(0, 0, A_11.get()); }
+   A_block->SetBlock(0, 1, D.get());
+   A_block->SetBlock(1, 0, Dt.get());
+   A_block->SetBlock(1, 1, R.Ptr(), -1.0);
+
+   D_prec.reset(new BlockDiagonalPreconditioner(offsets));
+   D_prec->SetDiagonalBlock(0, &S_inv);
+   D_prec->SetDiagonalBlock(1, R_inv.get());
+
+   minres.SetPreconditioner(*D_prec);
+   minres.SetOperator(*A_block);
+}
+
+void HdivSaddlePointSolver::EliminateBC(Vector &b) const
+{
+   const int n_ess_dofs = ess_rt_dofs.Size();
+   if (fes_l2.GetParMesh()->ReduceInt(n_ess_dofs) == 0) { return; }
+
+   const int n_l2 = offsets[1];
+   const int n_rt = offsets[2]-offsets[1];
+   Vector bE(b, 0, n_l2);
+   Vector bF(b, n_l2, n_rt);
+
+   // SetBC must be called first
+   MFEM_VERIFY(x_bc.Size() == n_rt || n_ess_dofs == 0, "BCs not set");
+
+   // Create a vector z that has the BC values at essential DOFs, zero elsewhere
+   z.SetSize(n_rt);
+   z.UseDevice(true);
+   z = 0.0;
+   const int *d_I = ess_rt_dofs.Read();
+   const double *d_x_bc = x_bc.Read();
+   double *d_z = z.ReadWrite();
+   MFEM_FORALL(i, n_ess_dofs,
+   {
+      const int j = d_I[i];
+      d_z[j] = d_x_bc[j];
+   });
+
+   // Convert to the IntegratedGLL basis used internally
+   w.SetSize(n_rt);
+   basis_rt.MultInverse(z, w);
+
+   // Eliminate the BCs in the L2 RHS
+   D_e->Mult(-1.0, w, 1.0, bE);
+
+   // Eliminate the BCs in the RT RHS
+   // Flip the sign because the R block appears with multiplier -1
+   z.SetSize(n_rt);
+   R_e->Mult(w, z);
+   bF += z;
+
+   // Insert the RT BCs into the RHS at the essential DOFs.
+   const double *d_w = w.Read();
+   double *d_bF = bF.ReadWrite(); // Need read-write access to set subvector
+   MFEM_FORALL(i, n_ess_dofs,
+   {
+      const int j = d_I[i];
+      d_bF[j] = -d_w[j];
+   });
+
+   // Make sure the monolithic RHS is updated
+   bE.SyncAliasMemory(b);
+   bF.SyncAliasMemory(b);
+}
+
+void HdivSaddlePointSolver::Mult(const Vector &b, Vector &x) const
+{
+   w.SetSize(fes_l2.GetTrueVSize());
+   b_prime.SetSize(b.Size());
+   x_prime.SetSize(x.Size());
+
+   // Transform RHS to the IntegratedGLL basis
+   Vector bE_prime(b_prime, offsets[0], offsets[1]-offsets[0]);
+   Vector bF_prime(b_prime, offsets[1], offsets[2]-offsets[1]);
+
+   const Vector bE(const_cast<Vector&>(b), offsets[0], offsets[1]-offsets[0]);
+   const Vector bF(const_cast<Vector&>(b), offsets[1], offsets[2]-offsets[1]);
+
+   z.SetSize(bE.Size());
+   basis_l2.MultTranspose(bE, z);
+   basis_rt.MultTranspose(bF, bF_prime);
+   // Transform by the inverse of the L2 mass matrix
+   L_inv->Mult(z, bE_prime);
+
+   // Update the monolithic transformed RHS
+   bE_prime.SyncAliasMemory(b_prime);
+   bF_prime.SyncAliasMemory(b_prime);
+
+   // Eliminate the RT essential BCs
+   EliminateBC(b_prime);
+
+   // Solve the transformed system
+   minres.Mult(b_prime, x_prime);
+
+   // Transform the solution back to the user's basis
+   Vector xE_prime(x_prime, offsets[0], offsets[1]-offsets[0]);
+   Vector xF_prime(x_prime, offsets[1], offsets[2]-offsets[1]);
+
+   Vector xE(x, offsets[0], offsets[1]-offsets[0]);
+   Vector xF(x, offsets[1], offsets[2]-offsets[1]);
+
+   z.SetSize(bE.Size()); // Size of z may have changed in EliminateBC
+   L_inv->Mult(xE_prime, z);
+
+   basis_l2.Mult(z, xE);
+   basis_rt.Mult(xF_prime, xF);
+
+   // Update the monolithic solution vector
+   xE.SyncAliasMemory(x);
+   xF.SyncAliasMemory(x);
+}
+
+} // namespace mfem

miniapps/hdiv-linear-solver/hdiv_linear_solver.hpp

@@ -0,0 +1,190 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef HDIV_LINEAR_SOLVER_HPP
+#define HDIV_LINEAR_SOLVER_HPP
+
+#include "mfem.hpp"
+#include "change_basis.hpp"
+#include <memory>
+
+namespace mfem
+{
+
+/// @brief Solve the H(div) saddle-point system using MINRES with matrix-free
+/// block-diagonal preconditioning.
+///
+/// See HdivSaddlePointSolver::HdivSaddlePointSolver for the problem
+/// description.
+class HdivSaddlePointSolver : public Solver
+{
+public:
+   /// Which type of saddle-point problem is being solved?
+   enum Mode
+   {
+      GRAD_DIV, ///< Grad-div problem.
+      DARCY     ///< Darcy/mixed Poisson problem.
+   };
+private:
+   MINRESSolver minres;
+
+   static constexpr int b1 = BasisType::GaussLobatto;
+   static constexpr int b2 = BasisType::IntegratedGLL;
+   static constexpr int mt = FiniteElement::INTEGRAL;
+
+   const int order;
+
+   // L2 and RT spaces, using the interpolation-histopolation bases
+   L2_FECollection fec_l2;
+   ParFiniteElementSpace fes_l2;
+
+   RT_FECollection fec_rt;
+   ParFiniteElementSpace fes_rt;
+
+   const Array<int> &ess_rt_dofs; ///< Essential BCs (in the RT space only).
+
+   // Change of basis operators
+   ChangeOfBasis_L2 basis_l2;
+   ChangeOfBasis_RT basis_rt;
+
+   /// Whether conversion from map type VALUE to INTEGRAL is required.
+   const bool convert_map_type;
+
+   ParBilinearForm mass_l2, mass_rt;
+
+   // Components needed for the block operator
+   OperatorHandle L, R, R_e; ///< Mass matrices.
+   std::unique_ptr<HypreParMatrix> D, Dt, D_e; ///< Divergence matrices.
+   std::shared_ptr<DGMassInverse> L_inv; ///< Inverse of the DG mass matrix.
+   std::shared_ptr<Operator> A_11; ///< (1,1)-block of the matrix
+
+   /// Diagonals of the mass matrices
+   Vector L_diag, R_diag, L_diag_unweighted;
+
+   // Components needed for the preconditioner
+
+   /// Jacobi preconditioner for the RT mass matrix.
+   std::unique_ptr<OperatorJacobiSmoother> R_inv;
+   std::unique_ptr<HypreParMatrix> S; ///< Approximate Schur complement.
+   HypreBoomerAMG S_inv; ///< AMG preconditioner for #S.
+
+   Array<int> offsets, empty;
+   /// The 2x2 block operator.
+   std::unique_ptr<BlockOperator> A_block;
+   /// The block-diagonal preconditioner.
+   std::unique_ptr<BlockDiagonalPreconditioner> D_prec;
+
+   Coefficient &L_coeff, &R_coeff, &B_coeff;
+
+   const Mode mode;
+   bool zero_l2_block = false;
+   QuadratureSpace qs;
+   QuadratureFunction W_coeff_qf, W_mix_coeff_qf;
+   QuadratureFunctionCoefficient W_coeff, W_mix_coeff;
+
+   ConstantCoefficient zero = ConstantCoefficient(0.0);
+
+   // Work vectors
+   mutable Vector b_prime, x_prime, x_bc, w, z;
+public:
+   /// @brief Creates a solver for the H(div) saddle-point system.
+   ///
+   /// The associated matrix is given by
+   ///
+   ///     [  L    B ]
+   ///     [ B^T  -R ]
+   ///
+   /// where L is the L2 mass matrix, R is the RT mass matrix, and B is the
+   /// divergence form (VectorFEDivergenceIntegrator).
+   ///
+   /// Essential boundary conditions in the RT space are given by @a
+   /// ess_rt_dofs_. (Rows and columns are eliminated from R and columns are
+   /// eliminated from B).
+   ///
+   /// The L block has coefficient @a L_coeff_ and the R block has coefficient
+   /// @a R_coeff_.
+   ///
+   /// The parameter @a mode_ determines whether the block system corresponds to
+   /// a grad-div problem or a Darcy problem. Specifically, if @a mode_ is
+   /// Mode::GRAD_DIV, then the B and B^T blocks are also scaled by @a L_coeff_,
+   /// and if @a mode_ is Mode::DARCY, then the B and B^T blocks are unweighted.
+   ///
+   /// Mode::GRAD_DIV corresponds to the grad-div problem
+   ///
+   ///     alpha u - grad ( beta div ( u )) = f,
+   ///
+   /// where alpha is @a R_coeff_ and beta is @a L_coeff_.
+   ///
+   /// Mode::DARCY corresponds to the Darcy-type problem
+   ///
+   ///     alpha p - div ( beta grad ( p )) = f,
+   ///
+   /// where alpha is @a L_coeff and beta is @a R_coeff_. In this case, the
+   /// coefficient alpha is allowed to be zero (see also @link
+   /// HdivSaddlePointSolver(ParMesh&, ParFiniteElementSpace&,
+   /// ParFiniteElementSpace&, Coefficient&, const Array<int>&) the zero-block
+   /// HdivSaddlePointSolver constructor@endlink).
+   HdivSaddlePointSolver(ParMesh &mesh_,
+                         ParFiniteElementSpace &fes_rt_,
+                         ParFiniteElementSpace &fes_l2_,
+                         Coefficient &L_coeff_,
+                         Coefficient &R_coeff_,
+                         Coefficient &B_coeff_,
+                         const Array<int> &ess_rt_dofs_,
+                         Mode mode_);
+
+   /// Same as the main constructor, but with B_coeff set to zero.
+   HdivSaddlePointSolver(ParMesh &mesh_,
+                         ParFiniteElementSpace &fes_rt_,
+                         ParFiniteElementSpace &fes_l2_,
+                         Coefficient &L_coeff_,
+                         Coefficient &R_coeff_,
+                         const Array<int> &ess_rt_dofs_,
+                         Mode mode_);
+
+   /// @brief Creates a linear solver for the case when the L2 diagonal block is
+   /// zero (for Darcy problems).
+   ///
+   /// Equivalent to passing ConstantCoefficient(0.0) as @a L_coeff_ and
+   /// Mode::DARCY as @a mode_ to the @link HdivSaddlePointSolver(ParMesh&,
+   /// ParFiniteElementSpace&, ParFiniteElementSpace&, Coefficient &,
+   /// Coefficient&, const Array<int>&, Mode) the primary constructor@endlink.
+   HdivSaddlePointSolver(ParMesh &mesh_,
+                         ParFiniteElementSpace &fes_rt_,
+                         ParFiniteElementSpace &fes_l2_,
+                         Coefficient &R_coeff_,
+                         const Array<int> &ess_rt_dofs_);
+
+   /// @brief Build the linear operator and solver. Must be called when the
+   /// coefficients change.
+   void Setup();
+   /// Sets the Dirichlet boundary conditions at the RT essential DOFs.
+   void SetBC(const Vector &x_rt) { x_bc = x_rt; }
+   /// @brief Solve the linear system for L2 (scalar) and RT (flux) unknowns.
+   ///
+   /// If the problem has essential boundary conditions (i.e. if @a ess_rt_dofs
+   /// is not empty), then SetBC() must be called before Mult().
+   void Mult(const Vector &b, Vector &x) const override;
+   /// No-op.
+   void SetOperator(const Operator &op) override { }
+   /// Get the number of MINRES iterations.
+   int GetNumIterations() const { return minres.GetNumIterations(); }
+   /// Eliminates the BCs (called internally, not public interface).
+   void EliminateBC(Vector &) const;
+   /// Return the offsets of the block system.
+   const Array<int> &GetOffsets() const { return offsets; }
+   /// Returns the internal MINRES solver.
+   MINRESSolver &GetMINRES() { return minres; }
+};
+
+} // namespace mfem
+
+#endif

miniapps/hdiv-linear-solver/makefile

@@ -0,0 +1,59 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/hdiv-linear-solver/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+HDIV_HEADERS = hdiv_linear_solver.hpp discrete_divergence.hpp change_basis.hpp \
+					../solvers/lor_mms.hpp
+HDIV_SRC = hdiv_linear_solver.cpp change_basis.cpp discrete_divergence.cpp
+HDIV_OBJ = $(HDIV_SRC:.cpp=.o)
+
+MINIAPPS = grad_div darcy residual rz
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+.PRECIOUS: %.o
+
+all: $(MINIAPPS)
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+$(MINIAPPS):%: %.o $(HDIV_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $(HDIV_OBJ) -o $@ $< $(MFEM_LIBS)
+
+%.o: $(SRC)%.cpp $(HDIV_HEADERS) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf mesh.* sol.* ParaView

miniapps/hdiv-linear-solver/residual.cpp

@@ -0,0 +1,205 @@
+#include "mfem.hpp"
+#include <iostream>
+
+#include "hdiv_linear_solver.hpp"
+#include "discrete_divergence.hpp"
+
+using namespace std;
+using namespace mfem;
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
+
+double f(const Vector &xvec);
+double g(const Vector &xvec);
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *mesh_file = "../../data/star.mesh";
+   const char *device_config = "cpu";
+   int ser_ref = 1;
+   int par_ref = 1;
+   int order = 3;
+   bool mt_value = true;
+   bool darcy = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+   args.AddOption(&par_ref, "-rp", "--parallel-refine",
+                  "Number of times to refine the mesh in parallel.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.AddOption(&mt_value, "-val", "--value", "-int", "--integral",
+                  "Map type integral or value.");
+   args.AddOption(&darcy, "-da", "--darcy", "-g", "--grad-div",
+                  "Grad-div or Darcy problem");
+   args.ParseCheck();
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
+   const int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+
+   const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
+   const int mt = mt_value ? FiniteElement::VALUE : FiniteElement::INTEGRAL;
+   RT_FECollection fec_rt(order-1, dim, b1, b2);
+   L2_FECollection fec_l2(order-1, dim, b2, mt);
+   ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
+   ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
+
+   HYPRE_BigInt ndofs_rt = fes_rt.GlobalTrueVSize();
+   HYPRE_BigInt ndofs_l2 = fes_l2.GlobalTrueVSize();
+
+   if (Mpi::Root())
+   {
+      cout << "\nRT DOFs: " << ndofs_rt << "\nL2 DOFs: " << ndofs_l2 << endl;
+   }
+
+   Array<int> ess_rt_dofs;
+
+   FunctionCoefficient a_coeff(f);
+   FunctionCoefficient b_coeff(g);
+   ConstantCoefficient one(1.0);
+
+   Coefficient &div_coeff = darcy ? (Coefficient&)one : (Coefficient&)a_coeff;
+
+   // Solve the system with the saddle-point solver
+   const auto solver_mode = darcy ? HdivSaddlePointSolver::Mode::DARCY
+                            : HdivSaddlePointSolver::Mode::GRAD_DIV;
+   HdivSaddlePointSolver saddle_point_solver(
+      mesh, fes_rt, fes_l2, a_coeff, b_coeff, ess_rt_dofs, solver_mode);
+
+   const Array<int> &offsets = saddle_point_solver.GetOffsets();
+   BlockVector X_block(offsets), B_block(offsets);
+
+   saddle_point_solver.GetMINRES().SetAbsTol(1e-18);
+   saddle_point_solver.GetMINRES().SetRelTol(1e-20);
+   saddle_point_solver.GetMINRES().SetPrintLevel(
+      IterativeSolver::PrintLevel().FirstAndLast());
+   X_block = 0.0;
+   B_block.Randomize(1);
+   B_block.GetBlock(0) = 0.0;
+   if (Mpi::Root()) { std::cout << "Saddle point solver... " << std::endl; }
+   saddle_point_solver.Mult(B_block, X_block);
+
+   // Form the matrix-based system
+   ParBilinearForm w(&fes_l2);
+   w.AddDomainIntegrator(new MassIntegrator(a_coeff));
+   w.Assemble();
+   w.Finalize();
+   std::unique_ptr<HypreParMatrix> W(w.ParallelAssemble());
+
+   ParMixedBilinearForm b(&fes_rt, &fes_l2);
+   b.AddDomainIntegrator(new VectorFEDivergenceIntegrator(div_coeff));
+   b.Assemble();
+   b.Finalize();
+   std::unique_ptr<HypreParMatrix> B(b.ParallelAssemble());
+   std::unique_ptr<HypreParMatrix> Bt(B->Transpose());
+
+   ParBilinearForm m(&fes_rt);
+   m.AddDomainIntegrator(new VectorFEMassIntegrator(b_coeff));
+   m.Assemble();
+   m.Finalize();
+   std::unique_ptr<HypreParMatrix> M(m.ParallelAssemble());
+
+   BlockOperator A(offsets);
+   A.SetBlock(0, 0, W.get());
+   A.SetBlock(0, 1, B.get());
+   A.SetBlock(1, 0, Bt.get());
+   A.SetBlock(1, 1, M.get(), -1.0);
+
+   // Compute the residual
+   BlockVector Y_block(offsets);
+   A.Mult(X_block, Y_block);
+   Y_block -= B_block;
+
+   auto nrm2 = [](const Vector &x)
+   {
+      return sqrt(InnerProduct(MPI_COMM_WORLD, x, x));
+   };
+
+   const double resnorm1 = nrm2(Y_block)/nrm2(B_block);
+   if (Mpi::Root()) { std::cout << "Linear residual norm: " << resnorm1 << "\n\n"; }
+
+   // Solve the system with a matrix-based solver (see ex5p)
+   HypreParVector Md(MPI_COMM_WORLD, M->GetGlobalNumRows(),
+                     M->GetRowStarts());
+   M->GetDiag(Md);
+   std::unique_ptr<HypreParMatrix> MinvBt(B->Transpose());
+   MinvBt->InvScaleRows(Md);
+   std::unique_ptr<HypreParMatrix> S(ParMult(B.get(), MinvBt.get()));
+
+   HypreDiagScale M_inv(*M);
+   HypreBoomerAMG S_inv(*S);
+   S_inv.SetPrintLevel(0);
+
+   BlockDiagonalPreconditioner D(offsets);
+   D.SetDiagonalBlock(0, &S_inv);
+   D.SetDiagonalBlock(1, &M_inv);
+
+   X_block = 0.0;
+   MINRESSolver minres(MPI_COMM_WORLD);
+   minres.SetAbsTol(1e-18);
+   minres.SetRelTol(1e-20);
+   minres.SetMaxIter(500);
+   minres.SetOperator(A);
+   minres.SetPreconditioner(D);
+   minres.SetPrintLevel(IterativeSolver::PrintLevel().FirstAndLast());
+   if (Mpi::Root()) { std::cout << "Matrix-based solver... " << std::endl; }
+   minres.Mult(B_block, X_block);
+
+   A.Mult(X_block, Y_block);
+   Y_block -= B_block;
+
+   const double resnorm2 = nrm2(Y_block)/nrm2(B_block);
+   if (Mpi::Root()) { std::cout << "Linear residual norm: " << resnorm2 << "\n\n"; }
+
+   return 0;
+}
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
+{
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
+   for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
+   return mesh;
+}
+
+double f(const Vector &xvec)
+{
+   const int dim = xvec.Size();
+   const double x = xvec[0], y = xvec[1];
+   if (dim == 2)
+   {
+      return 2*(2.0 + sin(x)*sin(y));
+   }
+   else // dim == 3
+   {
+      const double z = xvec[2];
+      return 3*(2.0 + sin(x)*sin(y)*sin(z));
+   }
+}
+
+double g(const Vector &xvec)
+{
+   const int dim = xvec.Size();
+   const double x = xvec[0], y = xvec[1];
+   if (dim == 2)
+   {
+      return 2*(2.0 + cos(x)*cos(y));
+   }
+   else // dim == 3
+   {
+      const double z = xvec[2];
+      return 3*(2.0 + cos(x)*cos(y)*cos(z));
+   }
+}

miniapps/hdiv-linear-solver/rz.cpp

@@ -0,0 +1,207 @@
+#include "mfem.hpp"
+#include <iostream>
+#include <memory>
+
+#include "discrete_divergence.hpp"
+#include "hdiv_linear_solver.hpp"
+
+#include "../solvers/lor_mms.hpp"
+
+using namespace std;
+using namespace mfem;
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
+
+double one_over_r(const Vector &xvec)
+{
+   // xvec = [z, r]
+   const double r = xvec[1];
+   return r == 0.0 ? 0.0 : 1.0/r;
+}
+
+double f_rz(const Vector &xvec)
+{
+   const double z = xvec[0];
+   const double r = xvec[1];
+
+   // alpha is the coefficient in the equation -Delta(u) + alpha*u = f
+   const double alpha = 1.0;
+   const double f = -cos(z)*(4*sin(r) + 5*r*cos(r) - (2 + alpha)*r*r*sin(r));
+
+   // scale integral by r because of coordinate transformation
+   return r*f;
+}
+
+double u_rz(const Vector &xvec)
+{
+   const double z = xvec[0];
+   const double r = xvec[1];
+
+   return r*r*sin(r)*cos(z);
+}
+
+class RobinCoefficient : public Coefficient
+{
+   double Eval(ElementTransformation &T, const IntegrationPoint &ip) override
+   {
+      double xdata[3];
+      Vector xvec(xdata, 3);
+      T.Transform(ip, xvec);
+      const int dim = xvec.Size();
+
+      Vector n(dim);
+      CalcOrtho(T.Jacobian(), n);
+      n /= n.Norml2();
+
+      const double p_val = u_rz(xvec);
+      const double z = xvec[0];
+      const double r = xvec[1];
+
+      if (dim == 2)
+      {
+         const double dpdz = -r*r*sin(r)*sin(z);
+         const double dpdr = r*cos(z)*(r*cos(r) + 2*sin(r));
+         const double u_val = n[0]*dpdz + n[1]*dpdr;
+         return p_val + u_val;
+      }
+      else
+      {
+         MFEM_ABORT("Not implemented");
+      }
+      return 0.0;
+   }
+};
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const char *mesh_file = "rz.mesh";
+   const char *device_config = "cpu";
+   int ser_ref = 1;
+   int par_ref = 1;
+   int order = 3;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ser_ref, "-rs", "--serial-refine",
+                  "Number of times to refine the mesh in serial.");
+   args.AddOption(&par_ref, "-rp", "--parallel-refine",
+                  "Number of times to refine the mesh in parallel.");
+   args.AddOption(&order, "-o", "--order", "Polynomial degree.");
+   args.ParseCheck();
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
+   const int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+
+   const int b1 = BasisType::GaussLobatto, b2 = BasisType::GaussLegendre;
+   const int mt = FiniteElement::VALUE;
+   RT_FECollection fec_rt(order-1, dim, b1, b2);
+   L2_FECollection fec_l2(order-1, dim, b2, mt);
+   ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
+   ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
+
+   HYPRE_BigInt ndofs_rt = fes_rt.GlobalTrueVSize();
+   HYPRE_BigInt ndofs_l2 = fes_l2.GlobalTrueVSize();
+
+   if (Mpi::Root())
+   {
+      cout << "\nRT DOFs: " << ndofs_rt << "\nL2 DOFs: " << ndofs_l2 << endl;
+   }
+
+   Array<int> ess_rt_dofs; // empty
+
+   // f is the RHS, u is the exact solution
+   FunctionCoefficient f_coeff(f_rz), u_coeff(u_rz);
+
+   // Assemble the right-hand side for the scalar (L2) unknown.
+   ParLinearForm b_l2(&fes_l2);
+   // f_coeff has to include the r scaling for the coordinate transformation
+   b_l2.AddDomainIntegrator(new DomainLFIntegrator(f_coeff));
+   b_l2.UseFastAssembly(true);
+   b_l2.Assemble();
+
+   // Coefficient to enforce Robin boundary condition
+   RobinCoefficient bc_coeff;
+
+   // Enforce Robin boundary conditions by adding the boundary term to the flux
+   // equation.
+   ParLinearForm b_rt(&fes_rt);
+   b_rt.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(bc_coeff));
+   b_rt.UseFastAssembly(true);
+   b_rt.Assemble();
+
+   if (Mpi::Root()) { cout << "\nSaddle point solver... " << flush; }
+   tic_toc.Clear(); tic_toc.Start();
+
+   // Have to scale the RT mass matrix by (1/r)
+   FunctionCoefficient r_recip_coeff(one_over_r);
+   // Have to scale the L2 mass matrix by r
+   FunctionCoefficient r_coeff([](const Vector &xvec) { return xvec[1]; });
+
+   const auto solver_mode = HdivSaddlePointSolver::Mode::DARCY;
+   HdivSaddlePointSolver saddle_point_solver(
+      mesh, fes_rt, fes_l2, r_coeff, r_recip_coeff, r_recip_coeff, ess_rt_dofs, solver_mode);
+
+   const Array<int> &offsets = saddle_point_solver.GetOffsets();
+   BlockVector X_block(offsets), B_block(offsets);
+
+   b_l2.ParallelAssemble(B_block.GetBlock(0));
+   b_rt.ParallelAssemble(B_block.GetBlock(1));
+   B_block.SyncFromBlocks();
+
+   X_block = 0.0;
+   saddle_point_solver.Mult(B_block, X_block);
+   X_block.SyncToBlocks();
+
+   if (Mpi::Root())
+   {
+      cout << "Done.\nIterations: "
+           << saddle_point_solver.GetNumIterations()
+           << "\nElapsed: " << tic_toc.RealTime() << endl;
+   }
+
+   ParGridFunction x(&fes_l2);
+   x.SetFromTrueDofs(X_block.GetBlock(0));
+   ParGridFunction flux(&fes_rt);
+   flux.SetFromTrueDofs(X_block.GetBlock(1));
+
+   const double error = x.ComputeL2Error(u_coeff);
+   if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
+
+   ParGridFunction u_ex(&fes_l2), er(&fes_l2);
+   u_ex.ProjectCoefficient(u_coeff);
+   er = x;
+   er -= u_ex;
+
+   ParaViewDataCollection pv("RZ", &mesh);
+   pv.SetPrefixPath("ParaView");
+   pv.SetHighOrderOutput(true);
+   pv.SetLevelsOfDetail(order + 1);
+   pv.RegisterField("u", &x);
+   pv.RegisterField("flux", &flux);
+   pv.RegisterField("exact", &u_ex);
+   pv.RegisterField("error", &er);
+   pv.SetCycle(0);
+   pv.SetTime(0);
+   pv.Save();
+
+   return 0;
+}
+
+ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
+{
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
+   for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
+   return mesh;
+}

miniapps/hdiv-linear-solver/rz.mesh

@@ -0,0 +1,35 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+1
+1 3 0 1 2 3
+
+boundary
+4
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 0
+
+vertices
+4
+2
+0 0
+12 0
+12 7.5
+0 7.5

miniapps/solvers/lor_mms.hpp

@@ -12,8 +12,6 @@
 #ifndef MFEM_LOR_MMS_HPP
 #define MFEM_LOR_MMS_HPP
 
-extern bool grad_div_problem;
-
 namespace mfem
 {
 
@@ -23,34 +21,36 @@ static constexpr double pi = M_PI, pi2 = M_PI*M_PI;
 // defined below.
 double u(const Vector &xvec)
 {
-   int dim = xvec.Size();
-   double x = pi*xvec[0], y = pi*xvec[1];
+   const int dim = xvec.Size();
+   const double x = pi*xvec[0], y = pi*xvec[1];
    if (dim == 2) { return sin(x)*sin(y); }
-   else { double z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
+   else { const double z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
 }
 
-double f(const Vector &xvec)
+std::function<double(const Vector &)> f(double mass_coeff)
 {
-   int dim = xvec.Size();
-   double x = pi*xvec[0], y = pi*xvec[1];
-
-   if (dim == 2)
+   return [mass_coeff](const Vector &xvec)
    {
-      return sin(x)*sin(y) + 2*pi2*sin(x)*sin(y);
-   }
-   else // dim == 3
-   {
-      double z = pi*xvec[2];
-      return sin(x)*sin(y)*sin(z) + 3*pi2*sin(x)*sin(y)*sin(z);
-   }
+      const int dim = xvec.Size();
+      const double x = pi*xvec[0], y = pi*xvec[1];
+      if (dim == 2)
+      {
+         return mass_coeff*sin(x)*sin(y) + 2*pi2*sin(x)*sin(y);
+      }
+      else // dim == 3
+      {
+         const double z = pi*xvec[2];
+         return mass_coeff*sin(x)*sin(y)*sin(z) + 3*pi2*sin(x)*sin(y)*sin(z);
+      }
+   };
 }
 
 // Exact solution for definite Maxwell and grad-div problems with RHS
 // corresponding to f_vec below.
 void u_vec(const Vector &xvec, Vector &u)
 {
-   int dim = xvec.Size();
-   double x = pi*xvec[0], y = pi*xvec[1];
+   const int dim = xvec.Size();
+   const double x = pi*xvec[0], y = pi*xvec[1];
    if (dim == 2)
    {
       u[0] = cos(x)*sin(y);
@@ -58,47 +58,50 @@ void u_vec(const Vector &xvec, Vector &u)
    }
    else // dim == 3
    {
-      double z = pi*xvec[2];
+      const double z = pi*xvec[2];
       u[0] = cos(x)*sin(y)*sin(z);
       u[1] = sin(x)*cos(y)*sin(z);
       u[2] = sin(x)*sin(y)*cos(z);
    }
 }
 
-void f_vec(const Vector &xvec, Vector &f)
+std::function<void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
 {
-   int dim = xvec.Size();
-   double x = pi*xvec[0], y = pi*xvec[1];
-   if (grad_div_problem)
+   return [grad_div_problem](const Vector &xvec, Vector &f)
    {
-      if (dim == 2)
+      const int dim = xvec.Size();
+      const double x = pi*xvec[0], y = pi*xvec[1];
+      if (grad_div_problem)
       {
-         f[0] = (1 + 2*pi2)*cos(x)*sin(y);
-         f[1] = (1 + 2*pi2)*cos(y)*sin(x);
+         if (dim == 2)
+         {
+            f[0] = (1 + 2*pi2)*cos(x)*sin(y);
+            f[1] = (1 + 2*pi2)*cos(y)*sin(x);
+         }
+         else // dim == 3
+         {
+            const double z = pi*xvec[2];
+            f[0] = (1 + 3*pi2)*cos(x)*sin(y)*sin(z);
+            f[1] = (1 + 3*pi2)*cos(y)*sin(x)*sin(z);
+            f[2] = (1 + 3*pi2)*cos(z)*sin(x)*sin(y);
+         }
       }
-      else // dim == 3
+      else
       {
-         double z = pi*xvec[2];
-         f[0] = (1 + 3*pi2)*cos(x)*sin(y)*sin(z);
-         f[1] = (1 + 3*pi2)*cos(y)*sin(x)*sin(z);
-         f[2] = (1 + 3*pi2)*cos(z)*sin(x)*sin(y);
+         if (dim == 2)
+         {
+            f[0] = cos(x)*sin(y);
+            f[1] = sin(x)*cos(y);
+         }
+         else // dim == 3
+         {
+            const double z = pi*xvec[2];
+            f[0] = cos(x)*sin(y)*sin(z);
+            f[1] = sin(x)*cos(y)*sin(z);
+            f[2] = sin(x)*sin(y)*cos(z);
+         }
       }
-   }
-   else
-   {
-      if (dim == 2)
-      {
-         f[0] = cos(x)*sin(y);
-         f[1] = sin(x)*cos(y);
-      }
-      else // dim == 3
-      {
-         double z = pi*xvec[2];
-         f[0] = cos(x)*sin(y)*sin(z);
-         f[1] = sin(x)*cos(y)*sin(z);
-         f[2] = sin(x)*sin(y)*cos(z);
-      }
-   }
+   };
 }
 
 } // namespace mfem

miniapps/solvers/lor_solvers.cpp

@@ -74,8 +74,6 @@
 using namespace std;
 using namespace mfem;
 
-bool grad_div_problem = false;
-
 int main(int argc, char *argv[])
 {
    const char *mesh_file = "../../data/star.mesh";
@@ -109,7 +107,6 @@ int main(int argc, char *argv[])
    else if (string(fe) == "l") { L2 = true; }
    else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
 
-   if (RT) { grad_div_problem = true; }
    double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh mesh(mesh_file, 1, 1);
@@ -117,8 +114,8 @@ int main(int argc, char *argv[])
    MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
    for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
 
-   FunctionCoefficient f_coeff(f), u_coeff(u);
-   VectorFunctionCoefficient f_vec_coeff(dim, f_vec), u_vec_coeff(dim, u_vec);
+   FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
+   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
 
    int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
    unique_ptr<FiniteElementCollection> fec;

miniapps/solvers/plor_solvers.cpp

@@ -72,8 +72,6 @@
 using namespace std;
 using namespace mfem;
 
-bool grad_div_problem = false;
-
 int main(int argc, char *argv[])
 {
    Mpi::Init();
@@ -112,7 +110,6 @@ int main(int argc, char *argv[])
    else if (string(fe) == "l") { L2 = true; }
    else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
 
-   if (RT) { grad_div_problem = true; }
    double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh serial_mesh(mesh_file, 1, 1);
@@ -126,8 +123,8 @@ int main(int argc, char *argv[])
    if (mesh.ncmesh && (RT || ND))
    { MFEM_ABORT("LOR AMS and ADS solvers are not supported with AMR meshes."); }
 
-   FunctionCoefficient f_coeff(f), u_coeff(u);
-   VectorFunctionCoefficient f_vec_coeff(dim, f_vec), u_vec_coeff(dim, u_vec);
+   FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
+   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
 
    int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
    unique_ptr<FiniteElementCollection> fec;