Merge branch 'master' into diag-smoothers-dev

fix 'make test' issues.

.gitignore

@@ -392,6 +392,12 @@ miniapps/spde/ParaView
 
 miniapps/tribol/contact-patch-test
 
+miniapps/diag-smoothers/lpq-jacobi
+miniapps/diag-smoothers/abs-l1-jacobi
+miniapps/diag-smoothers/mg-lpq-jacobi
+miniapps/diag-smoothers/mg-abs-l1-jacobi
+miniapps/diag-smoothers/meshes
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests

CHANGELOG

@@ -57,6 +57,11 @@ New and updated examples and miniapps
 - Added a command line option to all miniapps (`-p` or `--send-port`) for
   specifying the GLVis server socket port (19916 by default).
 
+- Added miniapps to demonstrate the utilization of L(p,q)-Jacobi preconditioners
+  in full assembly and absolute-value L(1)-Jacobi preconditioners in partially
+  assembled operators. This includes Multigrid wrappers to demonstrate the
+  effectiveness of these Jacobi-type operators as smoothers.
+
 GPU computing
 -------------
 - Added support for GPU-accelerated batched linear algebra (using cuBLAS,
@@ -128,6 +133,9 @@ API changes
 -----------
 - API change: in class GridFunction, 'fec' was renamed to 'fec_owned'.
 
+- API addition: in class `Operator`, added virtual functions: `AbsMult`,
+  `AbsMultTranspose`, `AddAbsMult`, and `AddAbsMultTranspose`.
+
 - API change: `RiemannSolver` was renamed to `NumericalFlux` (the old name has
   been been depracated through typedef)
 

fem/bilinearform_ext.cpp

@@ -525,7 +525,8 @@ void PABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
    A.Reset(oper); // A will own oper
 }
 
-void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
+void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
+                                           bool useAbs) const
 {
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
 
@@ -557,11 +558,13 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
       {
          if (integrators[i]->Patchwise())
          {
+            MFEM_ASSERT(!useAbs, "AbsMult not implemented with NURBS!");
             integrators[i]->AddMultNURBSPA(x, y);
          }
          else
          {
-            integrators[i]->AddMultPA(x, y);
+            if (!useAbs) { integrators[i]->AddMultPA(x, y); }
+            else { integrators[i]->AddAbsMultPA(x, y); }
          }
       }
    }
@@ -570,14 +573,32 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
       if (iSz)
       {
          Array<Array<int>*> &elem_markers = *a->GetDBFI_Marker();
-         elem_restrict->Mult(x, localX);
+
+         auto H1elem_restrict = dynamic_cast<const ElementRestriction*>(elem_restrict);
+         if (H1elem_restrict && useAbs)
+         {
+            H1elem_restrict->AbsMult(x, localX);
+         }
+         else
+         {
+            elem_restrict->Mult(x,localX);
+         }
+
          localY = 0.0;
          for (int i = 0; i < iSz; ++i)
          {
             AddMultWithMarkers(*integrators[i], localX, elem_markers[i],
-                               elem_attributes, false, localY);
+                               elem_attributes, false, localY, useAbs);
+         }
+
+         if (H1elem_restrict && useAbs)
+         {
+            H1elem_restrict->AbsMultTranspose(localY, y);
+         }
+         else
+         {
+            elem_restrict->MultTranspose(localY,y);
          }
-         elem_restrict->MultTranspose(localY, y);
       }
       else
       {
@@ -589,6 +610,7 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    const int iFISz = intFaceIntegrators.Size();
    if (int_face_restrict_lex && iFISz>0)
    {
+      MFEM_ASSERT(!useAbs, "AbsMult not implemented in face integrators!");
       // When assembling interior face integrators for DG spaces, we need to
       // exchange the face-neighbor information. This happens inside member
       // functions of the 'int_face_restrict_lex'. To avoid repeated calls to
@@ -650,6 +672,7 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    const bool has_bdr_integs = (n_bdr_face_integs > 0 || n_bdr_integs > 0);
    if (bdr_face_restrict_lex && has_bdr_integs)
    {
+      MFEM_ASSERT(!useAbs, "AbsMult not implemented in bdr integrators!");
       Array<Array<int>*> &bdr_markers = *a->GetBBFI_Marker();
       Array<Array<int>*> &bdr_face_markers = *a->GetBFBFI_Marker();
       bdr_face_restrict_lex->Mult(x, bdr_face_X);
@@ -827,22 +850,39 @@ void PABilinearFormExtension::AddMultWithMarkers(
    const Array<int> *markers,
    const Array<int> &attributes,
    const bool transpose,
-   Vector &y) const
+   Vector &y,
+   bool useAbs) const
 {
    if (markers)
    {
       tmp_evec.SetSize(y.Size());
       tmp_evec = 0.0;
-      if (transpose) { integ.AddMultTransposePA(x, tmp_evec); }
-      else { integ.AddMultPA(x, tmp_evec); }
+      if (!useAbs)
+      {
+         if (transpose) { integ.AddMultTransposePA(x, tmp_evec); }
+         else { integ.AddMultPA(x, tmp_evec); }
+      }
+      else
+      {
+         if (transpose) { integ.AddAbsMultTransposePA(x, tmp_evec); }
+         else { integ.AddAbsMultPA(x, tmp_evec); }
+      }
       const int ne = attributes.Size();
       const int nd = x.Size() / ne;
       AddWithMarkers_(ne, nd, tmp_evec, *markers, attributes, y);
    }
    else
    {
-      if (transpose) { integ.AddMultTransposePA(x, y); }
-      else { integ.AddMultPA(x, y); }
+      if (!useAbs)
+      {
+         if (transpose) { integ.AddMultTransposePA(x, y); }
+         else { integ.AddMultPA(x, y); }
+      }
+      else
+      {
+         if (transpose) { integ.AddAbsMultTransposePA(x, y); }
+         else { integ.AddAbsMultPA(x, y); }
+      }
    }
 }
 
@@ -924,8 +964,13 @@ void EABilinearFormExtension::Assemble()
    }
 }
 
-void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
+void EABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
+                                           bool useTranspose, bool useAbs) const
 {
+   auto el_rest = dynamic_cast<const ElementRestriction*>(elem_restrict);
+   MFEM_ASSERT(useAbs?(el_rest!=nullptr):true,
+               "elem_restrict is not ElementRestriction*!");
+   // MFEM_ASSERT(DeviceCanUseCeed() && useAbs, "AbsMult not implemented with CEED!");
    // Apply the Element Restriction
    const bool useRestrict = !DeviceCanUseCeed() && elem_restrict;
    if (!useRestrict)
@@ -933,6 +978,11 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
       y.UseDevice(true); // typically this is a large vector, so store on device
       y = 0.0;
    }
+   else if (useAbs)
+   {
+      el_rest->AbsMult(x, localX);
+      localY = 0.0;
+   }
    else
    {
       elem_restrict->Mult(x, localX);
@@ -940,25 +990,55 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    }
    // Apply the Element Matrices
    {
+      Vector abs_ea_data(ea_data.Size());
+      if (useAbs)
+      {
+         abs_ea_data = ea_data;
+         abs_ea_data.PowerAbs(1.0);
+      }
       const int NDOFS = elemDofs;
       auto X = Reshape(useRestrict?localX.Read():x.Read(), NDOFS, ne);
       auto Y = Reshape(useRestrict?localY.ReadWrite():y.ReadWrite(), NDOFS, ne);
-      auto A = Reshape(ea_data.Read(), NDOFS, NDOFS, ne);
-      mfem::forall(ne*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+      auto A = Reshape(useAbs?abs_ea_data.Read():ea_data.Read(), NDOFS, NDOFS, ne);
+      if (!useTranspose)
       {
-         const int e = glob_j/NDOFS;
-         const int j = glob_j%NDOFS;
-         real_t res = 0.0;
-         for (int i = 0; i < NDOFS; i++)
+         mfem::forall(ne*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
          {
-            res += A(i, j, e)*X(i, e);
-         }
-         Y(j, e) += res;
-      });
+            const int e = glob_j/NDOFS;
+            const int j = glob_j%NDOFS;
+            real_t res = 0.0;
+            for (int i = 0; i < NDOFS; i++)
+            {
+               res += A(i, j, e)*X(i, e);
+            }
+            Y(j, e) += res;
+         });
+      }
+      else
+      {
+         mfem::forall(ne*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+         {
+            const int e = glob_j/NDOFS;
+            const int j = glob_j%NDOFS;
+            real_t res = 0.0;
+            for (int i = 0; i < NDOFS; i++)
+            {
+               res += A(j, i, e)*X(i, e);
+            }
+            Y(j, e) += res;
+         });
+      }
       // Apply the Element Restriction transposed
       if (useRestrict)
       {
-         elem_restrict->MultTranspose(localY, y);
+         if (useAbs)
+         {
+            el_rest->AbsMultTranspose(localY, y);
+         }
+         else
+         {
+            elem_restrict->MultTranspose(localY, y);
+         }
       }
    }
 
@@ -967,7 +1047,9 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    const int iFISz = intFaceIntegrators.Size();
    if (int_face_restrict_lex && iFISz>0)
    {
+      MFEM_ASSERT(!useAbs, "AbsMult not implemented with Face integrators");
       // Apply the Interior Face Restriction
+      // TODO: AbsMult if needed
       int_face_restrict_lex->Mult(x, int_face_X);
       if (int_face_X.Size()>0)
       {
@@ -978,7 +1060,65 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
          auto Y = Reshape(int_face_Y.ReadWrite(), NDOFS, 2, nf_int);
          if (!factorize_face_terms)
          {
-            auto A_int = Reshape(ea_data_int.Read(), NDOFS, NDOFS, 2, nf_int);
+            Vector abs_ea_data_int(ea_data_int.Size());
+            if (useAbs)
+            {
+               abs_ea_data_int = ea_data_int;
+               abs_ea_data_int.PowerAbs(1.0);
+            }
+            auto A_int = Reshape(useAbs?abs_ea_data_int.Read():ea_data_int.Read(), NDOFS,
+                                 NDOFS, 2, nf_int);
+            if (!useTranspose)
+            {
+               mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+               {
+                  const int f = glob_j/NDOFS;
+                  const int j = glob_j%NDOFS;
+                  real_t res = 0.0;
+                  for (int i = 0; i < NDOFS; i++)
+                  {
+                     res += A_int(i, j, 0, f)*X(i, 0, f);
+                  }
+                  Y(j, 0, f) += res;
+                  res = 0.0;
+                  for (int i = 0; i < NDOFS; i++)
+                  {
+                     res += A_int(i, j, 1, f)*X(i, 1, f);
+                  }
+                  Y(j, 1, f) += res;
+               });
+            }
+            else
+            {
+               mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+               {
+                  const int f = glob_j/NDOFS;
+                  const int j = glob_j%NDOFS;
+                  real_t res = 0.0;
+                  for (int i = 0; i < NDOFS; i++)
+                  {
+                     res += A_int(j, i, 0, f)*X(i, 0, f);
+                  }
+                  Y(j, 0, f) += res;
+                  res = 0.0;
+                  for (int i = 0; i < NDOFS; i++)
+                  {
+                     res += A_int(j, i, 1, f)*X(i, 1, f);
+                  }
+                  Y(j, 1, f) += res;
+               });
+            }
+         }
+         Vector abs_ea_data_ext(ea_data_ext.Size());
+         if (useAbs)
+         {
+            abs_ea_data_ext = ea_data_ext;
+            abs_ea_data_ext.PowerAbs(1.0);
+         }
+         auto A_ext = Reshape(useAbs?abs_ea_data_ext.Read():ea_data_ext.Read(), NDOFS,
+                              NDOFS, 2, nf_int);
+         if (!useTranspose)
+         {
             mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
             {
                const int f = glob_j/NDOFS;
@@ -986,36 +1126,39 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
                real_t res = 0.0;
                for (int i = 0; i < NDOFS; i++)
                {
-                  res += A_int(i, j, 0, f)*X(i, 0, f);
+                  res += A_ext(i, j, 0, f)*X(i, 0, f);
                }
-               Y(j, 0, f) += res;
+               Y(j, 1, f) += res;
                res = 0.0;
                for (int i = 0; i < NDOFS; i++)
                {
-                  res += A_int(i, j, 1, f)*X(i, 1, f);
+                  res += A_ext(i, j, 1, f)*X(i, 1, f);
                }
-               Y(j, 1, f) += res;
+               Y(j, 0, f) += res;
             });
          }
-         auto A_ext = Reshape(ea_data_ext.Read(), NDOFS, NDOFS, 2, nf_int);
-         mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+         else
          {
-            const int f = glob_j/NDOFS;
-            const int j = glob_j%NDOFS;
-            real_t res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
+            mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
             {
-               res += A_ext(i, j, 0, f)*X(i, 0, f);
-            }
-            Y(j, 1, f) += res;
-            res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
-            {
-               res += A_ext(i, j, 1, f)*X(i, 1, f);
-            }
-            Y(j, 0, f) += res;
-         });
+               const int f = glob_j/NDOFS;
+               const int j = glob_j%NDOFS;
+               real_t res = 0.0;
+               for (int i = 0; i < NDOFS; i++)
+               {
+                  res += A_ext(j, i, 1, f)*X(i, 0, f);
+               }
+               Y(j, 1, f) += res;
+               res = 0.0;
+               for (int i = 0; i < NDOFS; i++)
+               {
+                  res += A_ext(j, i, 0, f)*X(i, 1, f);
+               }
+               Y(j, 0, f) += res;
+            });
+         }
          // Apply the Interior Face Restriction transposed
+         // TODO: AbsMultTranspose if needed
          int_face_restrict_lex->AddMultTransposeInPlace(int_face_Y, y);
       }
    }
@@ -1025,7 +1168,9 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    const int bFISz = bdrFaceIntegrators.Size();
    if (!factorize_face_terms && bdr_face_restrict_lex && bFISz>0)
    {
+      MFEM_ASSERT(!useAbs, "AbsMult not implemented with Face integrators");
       // Apply the Boundary Face Restriction
+      // TODO: AbsMult if needed
       bdr_face_restrict_lex->Mult(x, bdr_face_X);
       if (bdr_face_X.Size()>0)
       {
@@ -1034,147 +1179,44 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
          const int NDOFS = faceDofs;
          auto X = Reshape(bdr_face_X.Read(), NDOFS, nf_bdr);
          auto Y = Reshape(bdr_face_Y.ReadWrite(), NDOFS, nf_bdr);
-         auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
-         mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+         Vector abs_ea_data_bdr(ea_data_bdr.Size());
+         if (useAbs)
          {
-            const int f = glob_j/NDOFS;
-            const int j = glob_j%NDOFS;
-            real_t res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
-            {
-               res += A(i, j, f)*X(i, f);
-            }
-            Y(j, f) += res;
-         });
-         // Apply the Boundary Face Restriction transposed
-         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
-      }
-   }
-}
-
-void EABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
-{
-   // Apply the Element Restriction
-   const bool useRestrict = !DeviceCanUseCeed() && elem_restrict;
-   if (!useRestrict)
-   {
-      y.UseDevice(true); // typically this is a large vector, so store on device
-      y = 0.0;
-   }
-   else
-   {
-      elem_restrict->Mult(x, localX);
-      localY = 0.0;
-   }
-   // Apply the Element Matrices transposed
-   {
-      const int NDOFS = elemDofs;
-      auto X = Reshape(useRestrict?localX.Read():x.Read(), NDOFS, ne);
-      auto Y = Reshape(useRestrict?localY.ReadWrite():y.ReadWrite(), NDOFS, ne);
-      auto A = Reshape(ea_data.Read(), NDOFS, NDOFS, ne);
-      mfem::forall(ne*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
-      {
-         const int e = glob_j/NDOFS;
-         const int j = glob_j%NDOFS;
-         real_t res = 0.0;
-         for (int i = 0; i < NDOFS; i++)
-         {
-            res += A(j, i, e)*X(i, e);
+            abs_ea_data_bdr = ea_data_bdr;
+            abs_ea_data_bdr.PowerAbs(1.0);
          }
-         Y(j, e) += res;
-      });
-      // Apply the Element Restriction transposed
-      if (useRestrict)
-      {
-         elem_restrict->MultTranspose(localY, y);
-      }
-   }
-
-   // Treatment of interior faces
-   Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
-   const int iFISz = intFaceIntegrators.Size();
-   if (int_face_restrict_lex && iFISz>0)
-   {
-      // Apply the Interior Face Restriction
-      int_face_restrict_lex->Mult(x, int_face_X);
-      if (int_face_X.Size()>0)
-      {
-         int_face_Y = 0.0;
-         // Apply the interior face matrices transposed
-         const int NDOFS = faceDofs;
-         auto X = Reshape(int_face_X.Read(), NDOFS, 2, nf_int);
-         auto Y = Reshape(int_face_Y.ReadWrite(), NDOFS, 2, nf_int);
-         if (!factorize_face_terms)
+         auto A = Reshape(useAbs?abs_ea_data_bdr.Read():ea_data_bdr.Read(), NDOFS, NDOFS,
+                          nf_bdr);
+         if (!useTranspose)
          {
-            auto A_int = Reshape(ea_data_int.Read(), NDOFS, NDOFS, 2, nf_int);
-            mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+            mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
             {
                const int f = glob_j/NDOFS;
                const int j = glob_j%NDOFS;
                real_t res = 0.0;
                for (int i = 0; i < NDOFS; i++)
                {
-                  res += A_int(j, i, 0, f)*X(i, 0, f);
+                  res += A(i, j, f)*X(i, f);
                }
-               Y(j, 0, f) += res;
-               res = 0.0;
-               for (int i = 0; i < NDOFS; i++)
-               {
-                  res += A_int(j, i, 1, f)*X(i, 1, f);
-               }
-               Y(j, 1, f) += res;
+               Y(j, f) += res;
             });
          }
-         auto A_ext = Reshape(ea_data_ext.Read(), NDOFS, NDOFS, 2, nf_int);
-         mfem::forall(nf_int*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+         else
          {
-            const int f = glob_j/NDOFS;
-            const int j = glob_j%NDOFS;
-            real_t res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
+            mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
             {
-               res += A_ext(j, i, 1, f)*X(i, 0, f);
-            }
-            Y(j, 1, f) += res;
-            res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
-            {
-               res += A_ext(j, i, 0, f)*X(i, 1, f);
-            }
-            Y(j, 0, f) += res;
-         });
-         // Apply the Interior Face Restriction transposed
-         int_face_restrict_lex->AddMultTransposeInPlace(int_face_Y, y);
-      }
-   }
-
-   // Treatment of boundary faces
-   Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
-   const int bFISz = bdrFaceIntegrators.Size();
-   if (!factorize_face_terms && bdr_face_restrict_lex && bFISz>0)
-   {
-      // Apply the Boundary Face Restriction
-      bdr_face_restrict_lex->Mult(x, bdr_face_X);
-      if (bdr_face_X.Size()>0)
-      {
-         bdr_face_Y = 0.0;
-         // Apply the boundary face matrices transposed
-         const int NDOFS = faceDofs;
-         auto X = Reshape(bdr_face_X.Read(), NDOFS, nf_bdr);
-         auto Y = Reshape(bdr_face_Y.ReadWrite(), NDOFS, nf_bdr);
-         auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
-         mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
-         {
-            const int f = glob_j/NDOFS;
-            const int j = glob_j%NDOFS;
-            real_t res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
-            {
-               res += A(j, i, f)*X(i, f);
-            }
-            Y(j, f) += res;
-         });
+               const int f = glob_j/NDOFS;
+               const int j = glob_j%NDOFS;
+               real_t res = 0.0;
+               for (int i = 0; i < NDOFS; i++)
+               {
+                  res += A(j, i, f)*X(i, f);
+               }
+               Y(j, f) += res;
+            });
+         }
          // Apply the Boundary Face Restriction transposed
+         // TODO: AbsMult if needed
          bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
       }
    }

fem/bilinearform_ext.hpp

@@ -54,6 +54,11 @@ public:
       MFEM_ABORT("AssembleDiagonal not implemented for this assembly level!");
    }
 
+   void AbsMult(const Vector &x, Vector &y) const override
+   {
+      MFEM_ABORT("AbsMult not implemented for this assembly level!");
+   }
+
    virtual void FormSystemMatrix(const Array<int> &ess_tdof_list,
                                  OperatorHandle &A) = 0;
    virtual void FormLinearSystem(const Array<int> &ess_tdof_list,
@@ -91,12 +96,15 @@ public:
                          Vector &x, Vector &b,
                          OperatorHandle &A, Vector &X, Vector &B,
                          int copy_interior = 0) override;
-   void Mult(const Vector &x, Vector &y) const override;
+   void Mult(const Vector &x, Vector &y) const override { MultInternal(x,y); }
+   void AbsMult(const Vector &x, Vector &y) const override
+   { MultInternal(x,y,true); }
    void MultTranspose(const Vector &x, Vector &y) const override;
    void Update() override;
 
 protected:
    void SetupRestrictionOperators(const L2FaceValues m);
+   void MultInternal(const Vector &x, Vector &y, bool useAbs = false) const;
 
    /// @brief Accumulate the action (or transpose) of the integrator on @a x
    /// into @a y, taking into account the (possibly null) @a markers array.
@@ -110,12 +118,14 @@ protected:
    /// @param attributes Array of element or boundary element attributes.
    /// @param transpose Compute the action or transpose of the integrator .
    /// @param y Output E-vector
+   /// @param useAbs Apply absolute-value operator
    void AddMultWithMarkers(const BilinearFormIntegrator &integ,
                            const Vector &x,
                            const Array<int> *markers,
                            const Array<int> &attributes,
                            const bool transpose,
-                           Vector &y) const;
+                           Vector &y,
+                           bool useAbs = false) const;
 
    /// @brief Performs the same function as AddMultWithMarkers, but takes as
    /// input and output face normal derivatives.
@@ -152,8 +162,17 @@ public:
    EABilinearFormExtension(BilinearForm *form);
 
    void Assemble() override;
-   void Mult(const Vector &x, Vector &y) const override;
-   void MultTranspose(const Vector &x, Vector &y) const override;
+   void Mult(const Vector &x, Vector &y) const override
+   { MultInternal(x,y,false); }
+   void AbsMult(const Vector &x, Vector &y) const override
+   { MultInternal(x,y, false, true); }
+   void MultTranspose(const Vector &x, Vector &y) const override
+   { MultInternal(x,y,true); }
+   void AbsMultTranspose(const Vector &x, Vector &y) const override
+   { MultInternal(x,y,true, true); }
+protected:
+   void MultInternal(const Vector &x, Vector &y, bool useTranspose,
+                     bool useAbs = false) const;
 };
 
 /// Data and methods for fully-assembled bilinear forms

fem/bilininteg.cpp

@@ -103,18 +103,37 @@ void BilinearFormIntegrator::AddMultPA(const Vector &, Vector &) const
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AddAbsMultPA(const Vector &, Vector &) const
+{
+   MFEM_ABORT("BilinearFormIntegrator:AddAbsMultPA:(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AddMultNURBSPA(const Vector &, Vector &) const
 {
    MFEM_ABORT("BilinearFormIntegrator::AddMultNURBSPA(...)\n"
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AddAbsMultNURBSPA(const Vector &, Vector &) const
+{
+   MFEM_ABORT("BilinearFormIntegrator::AddAbsMultNURBSPA(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AddMultTransposePA(const Vector &, Vector &) const
 {
    MFEM_ABORT("BilinearFormIntegrator::AddMultTransposePA(...)\n"
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AddAbsMultTransposePA(const Vector &,
+                                                   Vector &) const
+{
+   MFEM_ABORT("BilinearFormIntegrator::AddAbsMultTransposePA(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AssembleMF(const FiniteElementSpace &fes)
 {
    MFEM_ABORT("BilinearFormIntegrator::AssembleMF(...)\n"
@@ -127,12 +146,25 @@ void BilinearFormIntegrator::AddMultMF(const Vector &, Vector &) const
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AddAbsMultMF(const Vector &, Vector &) const
+{
+   MFEM_ABORT("BilinearFormIntegrator::AddAbsMultMF(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AddMultTransposeMF(const Vector &, Vector &) const
 {
    MFEM_ABORT("BilinearFormIntegrator::AddMultTransposeMF(...)\n"
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AddAbsMultTransposeMF(const Vector &,
+                                                   Vector &) const
+{
+   MFEM_ABORT("BilinearFormIntegrator::AddAbsMultTransposeMF(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AssembleDiagonalMF(Vector &)
 {
    MFEM_ABORT("BilinearFormIntegrator::AssembleDiagonalMF(...)\n"
@@ -400,6 +432,14 @@ void SumIntegrator::AddMultPA(const Vector& x, Vector& y) const
    }
 }
 
+void SumIntegrator::AddAbsMultPA(const Vector& x, Vector& y) const
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AddAbsMultPA(x, y);
+   }
+}
+
 void SumIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
 {
    for (int i = 0; i < integrators.Size(); i++)
@@ -408,6 +448,14 @@ void SumIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
    }
 }
 
+void SumIntegrator::AddAbsMultTransposePA(const Vector &x, Vector &y) const
+{
+   for (int i = 0; i < integrators.Size(); i++)
+   {
+      integrators[i]->AddAbsMultTransposePA(x, y);
+   }
+}
+
 void SumIntegrator::AssembleMF(const FiniteElementSpace &fes)
 {
    for (int i = 0; i < integrators.Size(); i++)

fem/bilininteg.hpp

@@ -78,9 +78,13 @@ public:
        called. */
    void AddMultPA(const Vector &x, Vector &y) const override;
 
+   virtual void AddAbsMultPA(const Vector &x, Vector &y) const;
+
    /// Method for partially assembled action on NURBS patches.
    virtual void AddMultNURBSPA(const Vector&x, Vector&y) const;
 
+   virtual void AddAbsMultNURBSPA(const Vector&x, Vector&y) const;
+
    /// Method for partially assembled transposed action.
    /** Perform the transpose action of integrator on the input @a x and add the
        result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
@@ -90,6 +94,8 @@ public:
        called. */
    virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
 
+   virtual void AddAbsMultTransposePA(const Vector &x, Vector &y) const;
+
    /// Method defining element assembly.
    /** The result of the element assembly is added to the @a emat Vector if
        @a add is true. Otherwise, if @a add is false, we set @a emat. */
@@ -113,6 +119,8 @@ public:
        called. */
    void AddMultMF(const Vector &x, Vector &y) const override;
 
+   virtual void AddAbsMultMF(const Vector &x, Vector &y) const;
+
    /** Perform the transpose action of integrator on the input @a x and add the
        result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
        represent the element-wise discontinuous version of the FE space.
@@ -121,6 +129,8 @@ public:
        called. */
    virtual void AddMultTransposeMF(const Vector &x, Vector &y) const;
 
+   virtual void AddAbsMultTransposeMF(const Vector &x, Vector &y) const;
+
    /// Assemble diagonal and add it to Vector @a diag.
    virtual void AssembleDiagonalMF(Vector &diag);
 
@@ -481,8 +491,12 @@ public:
 
    void AddMultTransposePA(const Vector &x, Vector &y) const override;
 
+   void AddAbsMultTransposePA(const Vector &x, Vector &y) const override;
+
    void AddMultPA(const Vector& x, Vector& y) const override;
 
+   void AddAbsMultPA(const Vector& x, Vector& y) const override;
+
    void AssembleMF(const FiniteElementSpace &fes) override;
 
    void AddMultMF(const Vector &x, Vector &y) const override;
@@ -2304,8 +2318,12 @@ public:
 
    void AddMultPA(const Vector&, Vector&) const override;
 
+   void AddAbsMultPA(const Vector&, Vector&) const override;
+
    void AddMultTransposePA(const Vector&, Vector&) const override;
 
+   void AddAbsMultTransposePA(const Vector&, Vector&) const override;
+
    void AddMultNURBSPA(const Vector&, Vector&) const override;
 
    void AddMultPatchPA(const int patch, const Vector &x, Vector &y) const;
@@ -2398,8 +2416,12 @@ public:
 
    void AddMultPA(const Vector&, Vector&) const override;
 
+   void AddAbsMultPA(const Vector&, Vector&) const override;
+
    void AddMultTransposePA(const Vector&, Vector&) const override;
 
+   void AddAbsMultTransposePA(const Vector&, Vector&) const override;
+
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe,
                                          const ElementTransformation &Trans);
@@ -2795,6 +2817,7 @@ public:
    using BilinearFormIntegrator::AssemblePA;
    void AssemblePA(const FiniteElementSpace &fes) override;
    void AddMultPA(const Vector &x, Vector &y) const override;
+   void AddAbsMultPA(const Vector &x, Vector &y) const override;
    void AssembleDiagonalPA(Vector& diag) override;
 
    const Coefficient *GetCoefficient() const { return Q; }
@@ -2912,6 +2935,7 @@ public:
    void AssemblePA(const FiniteElementSpace &trial_fes,
                    const FiniteElementSpace &test_fes) override;
    void AddMultPA(const Vector &x, Vector &y) const override;
+   void AddAbsMultPA(const Vector &x, Vector &y) const override;
    void AddMultTransposePA(const Vector &x, Vector &y) const override;
    void AssembleDiagonalPA(Vector& diag) override;
 
@@ -3170,8 +3194,12 @@ public:
 
    void AddMultPA(const Vector &x, Vector &y) const override;
 
+   void AddAbsMultPA(const Vector &x, Vector &y) const override;
+
    void AddMultTransposePA(const Vector &x, Vector &y) const override;
 
+   void AddAbsMultTransposePA(const Vector &x, Vector &y) const override;
+
    /** Compute the stress corresponding to the local displacement @a $u$ and
        interpolate it at the nodes of the given @a fluxelem. Only the symmetric
        part of the stress is stored, so that the size of @a flux is equal to

fem/integ/bilininteg_curlcurl_pa.cpp

@@ -202,4 +202,80 @@ void CurlCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
+void CurlCurlIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
+{
+   Vector abs_pa_data(pa_data);
+   Array<real_t> absBo(mapsO->B);
+   Array<real_t> absBc(mapsC->B);
+   Array<real_t> absBto(mapsO->Bt);
+   Array<real_t> absBtc(mapsC->Bt);
+   Array<real_t> absGc(mapsC->G);
+   Array<real_t> absGtc(mapsC->Gt);
+
+   auto abs_val = static_cast<real_t(*)(real_t)>(std::abs);
+
+   abs_pa_data.PowerAbs(1.0);
+   absBo.Apply(abs_val);
+   absBc.Apply(abs_val);
+   absBto.Apply(abs_val);
+   absBtc.Apply(abs_val);
+   absGc.Apply(abs_val);
+   absGtc.Apply(abs_val);
+
+   if (dim == 3)
+   {
+      if (Device::Allows(Backend::DEVICE_MASK))
+      {
+         const int ID = (dofs1D << 4) | quad1D;
+         switch (ID)
+         {
+            case 0x23:
+               return internal::SmemPACurlCurlApply3D<2,3>(
+                         dofs1D, quad1D,
+                         symmetric, ne,
+                         absBo, absBc, absBto, absBtc,
+                         absGc, absGtc, abs_pa_data, x, y, true);
+            case 0x34:
+               return internal::SmemPACurlCurlApply3D<3,4>(
+                         dofs1D, quad1D,
+                         symmetric, ne,
+                         absBo, absBc, absBto, absBtc,
+                         absGc, absGtc, abs_pa_data, x, y, true);
+            case 0x45:
+               return internal::SmemPACurlCurlApply3D<4,5>(
+                         dofs1D, quad1D,
+                         symmetric, ne,
+                         absBo, absBc, absBto, absBtc,
+                         absGc, absGtc, abs_pa_data, x, y, true);
+            case 0x56:
+               return internal::SmemPACurlCurlApply3D<5,6>(
+                         dofs1D, quad1D,
+                         symmetric, ne,
+                         absBo, absBc, absBto, absBtc,
+                         absGc, absGtc, abs_pa_data, x, y, true);
+            default:
+               return internal::SmemPACurlCurlApply3D<0,0>(
+                         dofs1D, quad1D, symmetric, ne,
+                         absBo, absBc, absBto, absBtc,
+                         absGc, absGtc, abs_pa_data, x, y, true);
+         }
+      }
+      else
+      {
+         internal::PACurlCurlApply3D<0,0>(dofs1D, quad1D, symmetric, ne,
+                                          absBo, absBc, absBto, absBtc, absGc, absGtc,
+                                          abs_pa_data, x, y, true);
+      }
+   }
+   else if (dim == 2)
+   {
+      internal::PACurlCurlApply2D(dofs1D, quad1D, ne, absBo, absBto,
+                                  absGc, absGtc, abs_pa_data, x, y, true);
+   }
+   else
+   {
+      MFEM_ABORT("Unsupported dimension!");
+   }
+}
+
 } // namespace mfem

fem/integ/bilininteg_diffusion_kernels.hpp

@@ -482,19 +482,6 @@ inline void SmemPADiffusionDiagonal3D(const int NE,
    });
 }
 
-void PADiffusionApply(const int dim,
-                      const int D1D,
-                      const int Q1D,
-                      const int NE,
-                      const bool symm,
-                      const Array<real_t> &B,
-                      const Array<real_t> &G,
-                      const Array<real_t> &Bt,
-                      const Array<real_t> &Gt,
-                      const Vector &D,
-                      const Vector &X,
-                      Vector &Y);
-
 #ifdef MFEM_USE_OCCA
 // OCCA PA Diffusion Apply 2D kernel
 void OccaPADiffusionApply2D(const int D1D,

fem/integ/bilininteg_diffusion_pa.cpp

@@ -164,6 +164,47 @@ void DiffusionIntegrator::AssemblePatchPA(const int patch,
    SetupPatchPA(patch, mesh);  // For full quadrature, unitWeights = false
 }
 
+void DiffusionIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
+{
+   if (DeviceCanUseCeed())
+   {
+      MFEM_ABORT("What to do if Ceed?");
+      ceedOp->AddMult(x, y);
+   }
+   else
+   {
+      Vector abs_pa_data(pa_data);
+      abs_pa_data.PowerAbs(1.0);
+      Array<real_t> absB(maps->B);
+      Array<real_t> absG(maps->G);
+      Array<real_t> absBt(maps->Bt);
+      Array<real_t> absGt(maps->Gt);
+      auto abs_val = static_cast<real_t(*)(real_t)>(std::abs);
+      absB.Apply(abs_val);
+      absG.Apply(abs_val);
+      absBt.Apply(abs_val);
+      absGt.Apply(abs_val);
+
+      ApplyPAKernels::Run(dim, dofs1D, quad1D, ne, symmetric,
+                          absB, absG, absBt, absGt,
+                          abs_pa_data, x, y, dofs1D, quad1D);
+   }
+}
+
+void DiffusionIntegrator::AddAbsMultTransposePA(const Vector &x,
+                                                Vector &y) const
+{
+   if (symmetric)
+   {
+      AddAbsMultPA(x, y);
+   }
+   else
+   {
+      MFEM_ABORT("DiffusionIntegrator::AddAbsMultTransposePA only implemented "
+                 "in the symmetric case.")
+   }
+}
+
 // This version uses full 1D quadrature rules, taking into account the
 // minimum interaction between basis functions and integration points.
 void DiffusionIntegrator::AddMultPatchPA(const int patch, const Vector &x,

fem/integ/bilininteg_elasticity_kernels.cpp

@@ -88,6 +88,26 @@ void ElasticityAddMultPA(const int dim, const int nDofs,
    }
 }
 
+void ElasticityAddAbsMultPA(const int dim, const int nDofs,
+                            const FiniteElementSpace &fespace, const CoefficientVector &lambda,
+                            const CoefficientVector &mu, const GeometricFactors &geom,
+                            const DofToQuad &maps, const Vector &x, QuadratureFunction &QVec, Vector &y)
+{
+   switch (dim)
+   {
+      case 2:
+         ElasticityAddMultPA_<2>(nDofs, fespace, lambda, mu, geom, maps, x,
+                                 QVec, y, true);
+         break;
+      case 3:
+         ElasticityAddMultPA_<3>(nDofs, fespace, lambda, mu, geom, maps, x,
+                                 QVec, y, true);
+         break;
+      default:
+         MFEM_ABORT("Only dimensions 2 and 3 supported.");
+   }
+}
+
 void ElasticityAssembleDiagonalPA(const int dim, const int nDofs,
                                   const CoefficientVector &lambda,
                                   const CoefficientVector &mu, const GeometricFactors &geom,

fem/integ/bilininteg_elasticity_kernels.hpp

@@ -67,6 +67,11 @@ void ElasticityAddMultPA(const int dim, const int nDofs,
                          const CoefficientVector &mu, const GeometricFactors &geom,
                          const DofToQuad &maps, const Vector &x, QuadratureFunction &QVec, Vector &y);
 
+void ElasticityAddAbsMultPA(const int dim, const int nDofs,
+                            const FiniteElementSpace &fespace, const CoefficientVector &lambda,
+                            const CoefficientVector &mu, const GeometricFactors &geom,
+                            const DofToQuad &maps, const Vector &x, QuadratureFunction &QVec, Vector &y);
+
 /// @brief Elasticity component kernel for AddMultPA.
 ///
 /// Performs y += Ax. Implemented for byNODES ordering only, and does not use
@@ -144,7 +149,7 @@ template<int dim, int i_block = -1, int j_block = -1>
 void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
                           const CoefficientVector &lambda, const CoefficientVector &mu,
                           const GeometricFactors &geom, const DofToQuad &maps, const Vector &x,
-                          QuadratureFunction &QVec, Vector &y)
+                          QuadratureFunction &QVec, Vector &y, bool useAbs = false)
 {
    static_assert((i_block < 0) == (j_block < 0),
                  "i_block and j_block must both be non-negative or strictly negative.");
@@ -163,7 +168,8 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
                                                  ir);
    E_To_Q_Map->SetOutputLayout(QVectorLayout::byNODES);
    // interpolate physical derivatives to quadrature points.
-   E_To_Q_Map->PhysDerivatives(x, QVec);
+   if (!useAbs) { E_To_Q_Map->PhysDerivatives(x, QVec); }
+   else { E_To_Q_Map->AbsPhysDerivatives(x, QVec); }
 
    const int numPoints = ir.GetNPoints();
    const int numEls = fespace.GetNE();
@@ -172,6 +178,7 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
    const auto J = Reshape(geom.J.Read(), numPoints, d, d, numEls);
    auto Q = Reshape(QVec.ReadWrite(), numPoints, d, qSize, numEls);
    const real_t *ipWeights = ir.GetWeights().Read();
+
    mfem::forall_2D(numEls, numPoints, 1, [=] MFEM_HOST_DEVICE (int e)
    {
       // for(int p = 0; p < numPoints, )
@@ -206,7 +213,8 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
             const int iIndex = isComponent ? 0 : i;
             div += gradx(iIndex,i);
          }
-         const real_t w = ipWeights[p] /det(invJ);
+         const real_t w = (!useAbs) ? ipWeights[p] /det(invJ) :
+                          std::abs(ipWeights[p] /det(invJ));
          for (int m = 0; m < d; m++)
          {
             for (int q = qLower; q < qUpper; q++)
@@ -220,8 +228,16 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
                {
                   for (int a = 0; a < d; a++)
                   {
-                     contraction += 2*((a == q)*invJ(m,j_block) + (j_block==q)*invJ(m,a))*(gradx(0,
-                                                                                                 a));
+                     if (!useAbs)
+                     {
+                        contraction += 2*((a == q)*invJ(m,j_block) +
+                                          (j_block==q)*invJ(m,a))*(gradx(0,a));
+                     }
+                     else
+                     {
+                        contraction += 2*((a == q)*std::abs(invJ(m,j_block)) +
+                                          (j_block==q)*std::abs(invJ(m,a)))*(gradx(0,a));
+                     }
                   }
                }
                else
@@ -230,15 +246,31 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
                   {
                      for (int b = 0; b < d; b++)
                      {
-                        contraction += ((a == q)*invJ(m,b) + (b==q)*invJ(m,a))
-                                       *(gradx(a,b) + gradx(b, a));
+                        if (!useAbs)
+                        {
+                           contraction += ((a == q)*invJ(m,b) + (b==q)*invJ(m,a))
+                                          *(gradx(a,b) + gradx(b, a));
+                        }
+                        else
+                        {
+                           contraction += ((a == q)*std::abs(invJ(m,b)) + (b==q)*std::abs(invJ(m,a)))
+                                          *(gradx(a,b) + gradx(b, a));
+                        }
                      }
                   }
                }
                // lambda*div(u)*div(v) + 2*mu*sym(grad(u))*sym(grad(v))
                // contraction = 4*sym(grad(u))sym(grad(v))
                const int qIndex = isComponent ? 0 : q;
-               Q(p,m,qIndex,e) = w*(lamDev(p, e)*invJ(m,q)*div + 0.5*muDev(p, e)*contraction);
+               if (!useAbs)
+               {
+                  Q(p,m,qIndex,e) = w*(lamDev(p, e)*invJ(m,q)*div + 0.5*muDev(p, e)*contraction);
+               }
+               else
+               {
+                  Q(p,m,qIndex,e) = w*(std::abs(lamDev(p, e)*invJ(m,q))*div +
+                                       0.5*std::abs(muDev(p, e))*contraction);
+               }
             }
          }
       }

fem/integ/bilininteg_elasticity_pa.cpp

@@ -75,6 +75,36 @@ void ElasticityIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
    AddMultPA(x, y); // Operator is symmetric
 }
 
+void ElasticityIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
+{
+   DofToQuad abs_maps;
+
+   abs_maps.FE = maps->FE;
+   abs_maps.IntRule = maps->IntRule;
+   abs_maps.mode = maps->mode;
+   abs_maps.ndof = maps->ndof;
+   abs_maps.nqpt = maps->nqpt;
+
+   abs_maps.B = maps->B;
+   abs_maps.Bt = maps->Bt;
+   abs_maps.G = maps->G;
+   abs_maps.Gt = maps->Gt;
+   auto abs_val = static_cast<real_t(*)(real_t)>(std::abs);
+   abs_maps.B.Apply(abs_val);
+   abs_maps.G.Apply(abs_val);
+   abs_maps.Bt.Apply(abs_val);
+   abs_maps.Gt.Apply(abs_val);
+
+   internal::ElasticityAddAbsMultPA(vdim, ndofs, *fespace, *lambda_quad, *mu_quad,
+                                    *geom, abs_maps, x, *q_vec, y);
+}
+
+void ElasticityIntegrator::AddAbsMultTransposePA(const Vector &x,
+                                                 Vector &y) const
+{
+   AddAbsMultPA(x, y); // Operator is symmetric
+}
+
 void ElasticityComponentIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
    fespace = &fes;

fem/integ/bilininteg_hcurl_kernels.cpp

@@ -662,7 +662,8 @@ void PACurlCurlApply2D(const int D1D,
                        const Array<real_t> &gct,
                        const Vector &pa_data,
                        const Vector &x,
-                       Vector &y)
+                       Vector &y,
+                       bool useAbs)
 {
 
    auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
@@ -717,7 +718,8 @@ void PACurlCurlApply2D(const int D1D,
 
             for (int qy = 0; qy < Q1D; ++qy)
             {
-               const real_t wy = (c == 0) ? -Gc(qy,dy) : Bo(qy,dy);
+               const int sign = useAbs ? 1 : -1;
+               const real_t wy = (c == 0) ? (sign*Gc(qy,dy)) : Bo(qy,dy);
                for (int qx = 0; qx < Q1D; ++qx)
                {
                   curl[qy][qx] += gradX[qx] * wy;
@@ -760,7 +762,8 @@ void PACurlCurlApply2D(const int D1D,
             }
             for (int dy = 0; dy < D1Dy; ++dy)
             {
-               const real_t wy = (c == 0) ? -Gct(dy,qy) : Bot(dy,qy);
+               const int sign = useAbs ? 1 : -1;
+               const real_t wy = (c == 0) ? sign*Gct(dy,qy) : Bot(dy,qy);
 
                for (int dx = 0; dx < D1Dx; ++dx)
                {

fem/integ/bilininteg_hcurl_kernels.hpp

@@ -828,7 +828,7 @@ inline void SmemPACurlCurlAssembleDiagonal3D(const int d1d,
    }); // end of element loop
 }
 
-// PA H(curl) curl-curl Apply 2D kernel
+// PA H(curl) curl-curl Apply 2D kernel (and AbsApply 2D kernel)
 void PACurlCurlApply2D(const int D1D,
                        const int Q1D,
                        const int NE,
@@ -838,9 +838,10 @@ void PACurlCurlApply2D(const int D1D,
                        const Array<real_t> &gct,
                        const Vector &pa_data,
                        const Vector &x,
-                       Vector &y);
+                       Vector &y,
+                       bool useAbs = false);
 
-// PA H(curl) curl-curl Apply 3D kernel
+// PA H(curl) curl-curl Apply 3D kernel (and AbsApply 3D kernel)
 template<int T_D1D = 0, int T_Q1D = 0>
 inline void PACurlCurlApply3D(const int d1d,
                               const int q1d,
@@ -854,7 +855,8 @@ inline void PACurlCurlApply3D(const int d1d,
                               const Array<real_t> &gct,
                               const Vector &pa_data,
                               const Vector &x,
-                              Vector &y)
+                              Vector &y,
+                              bool useAbs = false)
 {
    MFEM_VERIFY(T_D1D || d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
                "Error: d1d > HCURL_MAX_D1D");
@@ -970,7 +972,8 @@ inline void PACurlCurlApply3D(const int d1d,
                   {
                      // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
                      curl[qz][qy][qx][1] += gradXY[qy][qx][1] * wDz; // (u_0)_{x_2}
-                     curl[qz][qy][qx][2] -= gradXY[qy][qx][0] * wz;  // -(u_0)_{x_1}
+                     if (!useAbs) { curl[qz][qy][qx][2] -= gradXY[qy][qx][0] * wz; } // -(u_0)_{x_1}
+                     else { curl[qz][qy][qx][2] += gradXY[qy][qx][0] * wz; }    // +(u_0)_{x_1}
                   }
                }
             }
@@ -1038,7 +1041,8 @@ inline void PACurlCurlApply3D(const int d1d,
                   for (int qx = 0; qx < Q1D; ++qx)
                   {
                      // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
-                     curl[qz][qy][qx][0] -= gradXY[qy][qx][1] * wDz; // -(u_1)_{x_2}
+                     if (!useAbs) { curl[qz][qy][qx][0] -= gradXY[qy][qx][1] * wDz; } // -(u_1)_{x_2}
+                     else { curl[qz][qy][qx][0] += gradXY[qy][qx][1] * wDz; }      // +(u_1)_{x_2}
                      curl[qz][qy][qx][2] += gradXY[qy][qx][0] * wz;  // (u_1)_{x_0}
                   }
                }
@@ -1109,7 +1113,8 @@ inline void PACurlCurlApply3D(const int d1d,
                   {
                      // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
                      curl[qz][qy][qx][0] += gradYZ[qz][qy][1] * wx;  // (u_2)_{x_1}
-                     curl[qz][qy][qx][1] -= gradYZ[qz][qy][0] * wDx; // -(u_2)_{x_0}
+                     if (!useAbs) { curl[qz][qy][qx][1] -= gradYZ[qz][qy][0] * wDx; } // -(u_2)_{x_0}
+                     else { curl[qz][qy][qx][1] += gradYZ[qz][qy][0] * wDx; }      // +(u_2)_{x_0}
                   }
                }
             }
@@ -1210,8 +1215,16 @@ inline void PACurlCurlApply3D(const int d1d,
                   {
                      // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
                      // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
-                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
-                       e) += (gradXY21[dy][dx] * wDz) - (gradXY12[dy][dx] * wz);
+                     if (!useAbs)
+                     {
+                        Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                          e) += (gradXY21[dy][dx] * wDz) - (gradXY12[dy][dx] * wz);
+                     }
+                     else
+                     {
+                        Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                          e) += (gradXY21[dy][dx] * wDz) + (gradXY12[dy][dx] * wz);
+                     }
                   }
                }
             }
@@ -1280,8 +1293,16 @@ inline void PACurlCurlApply3D(const int d1d,
                   {
                      // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
                      // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
-                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
-                       e) += (-gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
+                     if (!useAbs)
+                     {
+                        Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                          e) += (-gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
+                     }
+                     else
+                     {
+                        Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                          e) += (gradXY20[dy][dx] * wDz) + (gradXY02[dy][dx] * wz);
+                     }
                   }
                }
             }
@@ -1353,8 +1374,16 @@ inline void PACurlCurlApply3D(const int d1d,
                   {
                      // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
                      // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
-                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
-                       e) += (gradYZ10[dz][dy] * wx) - (gradYZ01[dz][dy] * wDx);
+                     if (!useAbs)
+                     {
+                        Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                          e) += (gradYZ10[dz][dy] * wx) - (gradYZ01[dz][dy] * wDx);
+                     }
+                     else
+                     {
+                        Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc,
+                          e) += (gradYZ10[dz][dy] * wx) + (gradYZ01[dz][dy] * wDx);
+                     }
                   }
                }
             }
@@ -1363,7 +1392,7 @@ inline void PACurlCurlApply3D(const int d1d,
    }); // end of element loop
 }
 
-// Shared memory PA H(curl) curl-curl Apply 3D kernel
+// Shared memory PA H(curl) curl-curl Apply 3D kernel (and AbsApply 3D kernel)
 template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPACurlCurlApply3D(const int d1d,
                                   const int q1d,
@@ -1377,7 +1406,8 @@ inline void SmemPACurlCurlApply3D(const int d1d,
                                   const Array<real_t> &gct,
                                   const Vector &pa_data,
                                   const Vector &x,
-                                  Vector &y)
+                                  Vector &y,
+                                  bool useAbs = false)
 {
    MFEM_VERIFY(T_D1D || d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
                "Error: d1d > HCURL_MAX_D1D");
@@ -1531,7 +1561,14 @@ inline void SmemPACurlCurlApply3D(const int d1d,
                         }
 
                         curl[qy][qx][1] += v; // (u_0)_{x_2}
-                        curl[qy][qx][2] -= u;  // -(u_0)_{x_1}
+                        if (!useAbs)
+                        {
+                           curl[qy][qx][2] -= u;  // -(u_0)_{x_1}
+                        }
+                        else
+                        {
+                           curl[qy][qx][2] += u;  // +(u_0)_{x_1}
+                        }
                      }
                      else if (c == 1)  // y component
                      {
@@ -1558,7 +1595,14 @@ inline void SmemPACurlCurlApply3D(const int d1d,
                            }
                         }
 
-                        curl[qy][qx][0] -= v; // -(u_1)_{x_2}
+                        if (!useAbs)
+                        {
+                           curl[qy][qx][0] -= v; // -(u_1)_{x_2}
+                        }
+                        else
+                        {
+                           curl[qy][qx][0] += v; // +(u_1)_{x_2}
+                        }
                         curl[qy][qx][2] += u; // (u_1)_{x_0}
                      }
                      else // z component
@@ -1587,7 +1631,14 @@ inline void SmemPACurlCurlApply3D(const int d1d,
                         }
 
                         curl[qy][qx][0] += v; // (u_2)_{x_1}
-                        curl[qy][qx][1] -= u; // -(u_2)_{x_0}
+                        if (!useAbs)
+                        {
+                           curl[qy][qx][1] -= u; // -(u_2)_{x_0}
+                        }
+                        else
+                        {
+                           curl[qy][qx][1] += u; // +(u_2)_{x_0}
+                        }
                      }
                   } // qx
                } // qy
@@ -1644,16 +1695,37 @@ inline void SmemPACurlCurlApply3D(const int d1d,
                            // \hat{\nabla}\times\hat{u} is [0, (u_0)_{x_2}, -(u_0)_{x_1}]
                            // (u_0)_{x_2} * (op * curl)_1 - (u_0)_{x_1} * (op * curl)_2
                            const real_t wx = sBo[dx][qx];
-                           dxyz1 += (wx * c2 * wcy * wcDz) - (wx * c3 * wcDy * wcz);
+                           if (!useAbs)
+                           {
+                              dxyz1 += (wx * c2 * wcy * wcDz) - (wx * c3 * wcDy * wcz);
+                           }
+                           else
+                           {
+                              dxyz1 += (wx * c2 * wcy * wcDz) + (wx * c3 * wcDy * wcz);
+                           }
                         }
 
                         // \hat{\nabla}\times\hat{u} is [-(u_1)_{x_2}, 0, (u_1)_{x_0}]
                         // -(u_1)_{x_2} * (op * curl)_0 + (u_1)_{x_0} * (op * curl)_2
-                        dxyz2 += (-wy * c1 * wcx * wcDz) + (wy * c3 * wDx * wcz);
+                        if (!useAbs)
+                        {
+                           dxyz2 += (-wy * c1 * wcx * wcDz) + (wy * c3 * wDx * wcz);
+                        }
+                        else
+                        {
+                           dxyz2 += (wy * c1 * wcx * wcDz) + (wy * c3 * wDx * wcz);
+                        }
 
                         // \hat{\nabla}\times\hat{u} is [(u_2)_{x_1}, -(u_2)_{x_0}, 0]
                         // (u_2)_{x_1} * (op * curl)_0 - (u_2)_{x_0} * (op * curl)_1
-                        dxyz3 += (wcDy * wz * c1 * wcx) - (wcy * wz * c2 * wDx);
+                        if (!useAbs)
+                        {
+                           dxyz3 += (wcDy * wz * c1 * wcx) - (wcy * wz * c2 * wDx);
+                        }
+                        else
+                        {
+                           dxyz3 += (wcDy * wz * c1 * wcx) + (wcy * wz * c2 * wDx);
+                        }
                      } // qx
                   } // qy
                } // dx

fem/integ/bilininteg_mass_pa.cpp

@@ -232,10 +232,38 @@ void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
+void MassIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
+{
+   if (DeviceCanUseCeed())
+   {
+      MFEM_ABORT("AddAbsMultPA not implemented with CEED!");
+      ceedOp->AddMult(x, y);
+   }
+   else
+   {
+      Vector abs_pa_data(pa_data);
+      abs_pa_data.PowerAbs(1.0);
+      Array<real_t> absB(maps->B);
+      Array<real_t> absBt(maps->Bt);
+      auto abs_val = static_cast<real_t(*)(real_t)>(std::abs);
+      absB.Apply(abs_val);
+      absBt.Apply(abs_val);
+
+      ApplyPAKernels::Run(dim, dofs1D, quad1D, ne, absB, absBt, abs_pa_data,
+                          x, y, dofs1D, quad1D);
+   }
+}
+
 void MassIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
 {
    // Mass integrator is symmetric
    AddMultPA(x, y);
 }
 
+void MassIntegrator::AddAbsMultTransposePA(const Vector &x, Vector &y) const
+{
+   // Mass integrator is symmetric
+   AddAbsMultPA(x, y);
+}
+
 } // namespace mfem

fem/integ/bilininteg_vectorfemass_pa.cpp

@@ -313,6 +313,127 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
+void VectorFEMassIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
+{
+   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
+   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
+   const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
+   const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
+
+   Vector abs_pa_data(pa_data);
+   abs_pa_data.PowerAbs(1.0);
+
+   Array<real_t> absBo(mapsO->B);
+   Array<real_t> absBc(mapsC->B);
+   Array<real_t> absBto(mapsO->Bt);
+   Array<real_t> absBtc(mapsC->Bt);
+   Array<real_t> absBto_t(mapsOtest->Bt);
+   Array<real_t> absBtc_t(mapsCtest->Bt);
+   auto abs_val = static_cast<real_t(*)(real_t)>(std::abs);
+   absBo.Apply(abs_val);
+   absBc.Apply(abs_val);
+   absBto.Apply(abs_val);
+   absBtc.Apply(abs_val);
+   absBto_t.Apply(abs_val);
+   absBtc_t.Apply(abs_val);
+
+   if (dim == 3)
+   {
+      if (trial_curl && test_curl)
+      {
+         if (Device::Allows(Backend::DEVICE_MASK))
+         {
+            const int ID = (dofs1D << 4) | quad1D;
+            switch (ID)
+            {
+               case 0x23:
+                  return internal::SmemPAHcurlMassApply3D<2,3>(
+                            dofs1D, quad1D, ne, symmetric,
+                            absBo, absBc, absBto, absBtc,
+                            abs_pa_data, x, y);
+               case 0x34:
+                  return internal::SmemPAHcurlMassApply3D<3,4>(
+                            dofs1D, quad1D, ne, symmetric,
+                            absBo, absBc, absBto, absBtc,
+                            abs_pa_data, x, y);
+               case 0x45:
+                  return internal::SmemPAHcurlMassApply3D<4,5>(
+                            dofs1D, quad1D, ne, symmetric,
+                            absBo, absBc, absBto, absBtc,
+                            abs_pa_data, x, y);
+               case 0x56:
+                  return internal::SmemPAHcurlMassApply3D<5,6>(
+                            dofs1D, quad1D, ne, symmetric,
+                            absBo, absBc, absBto, absBtc,
+                            abs_pa_data, x, y);
+               default:
+                  return internal::SmemPAHcurlMassApply3D(
+                            dofs1D, quad1D, ne, symmetric,
+                            absBo, absBc, absBto, absBtc,
+                            abs_pa_data, x, y);
+            }
+         }
+         else
+         {
+            internal::PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric,
+                                         absBo, absBc, absBto, absBtc,
+                                         abs_pa_data, x, y);
+         }
+      }
+      else if (trial_div && test_div)
+      {
+         internal::PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric,
+                                   absBo, absBc, absBto, absBtc,
+                                   abs_pa_data, x, y);
+      }
+      else if (trial_curl && test_div)
+      {
+         const bool scalarCoeff = !(DQ || MQ);
+         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                                          true, false, absBo, absBc, absBto_t,
+                                          absBtc_t, abs_pa_data, x, y);
+      }
+      else if (trial_div && test_curl)
+      {
+         const bool scalarCoeff = !(DQ || MQ);
+         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                                          false, false, absBo, absBc, absBto_t,
+                                          absBtc_t, abs_pa_data, x, y);
+      }
+      else
+      {
+         MFEM_ABORT("Unknown kernel.");
+      }
+   }
+   else // 2D
+   {
+      if (trial_curl && test_curl)
+      {
+         internal::PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric,
+                                      absBo, absBc, absBto, absBtc,
+                                      abs_pa_data, x, y);
+      }
+      else if (trial_div && test_div)
+      {
+         internal::PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric,
+                                   absBo, absBc, absBto, absBtc,
+                                   abs_pa_data, x, y);
+      }
+      else if ((trial_curl && test_div) || (trial_div && test_curl))
+      {
+         const bool scalarCoeff = !(DQ || MQ);
+         internal::PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                                          trial_curl, false,
+                                          absBo, absBc, absBto_t, absBtc_t,
+                                          abs_pa_data, x, y);
+      }
+      else
+      {
+         MFEM_ABORT("Unknown kernel.");
+      }
+   }
+}
+
 void VectorFEMassIntegrator::AddMultTransposePA(const Vector &x,
                                                 Vector &y) const
 {

fem/pfespace.hpp

@@ -465,7 +465,11 @@ public:
 
    void Mult(const Vector &x, Vector &y) const override;
 
+   void AbsMult(const Vector &x, Vector &y) const override { Mult(x,y); }
+
    void MultTranspose(const Vector &x, Vector &y) const override;
+
+   void AbsMultTranspose(const Vector &x, Vector &y) const override { MultTranspose(x,y); }
 };
 
 /// Auxiliary device class used by ParFiniteElementSpace.
@@ -516,7 +520,11 @@ public:
 
    void Mult(const Vector &x, Vector &y) const override;
 
+   void AbsMult(const Vector &x, Vector &y) const override { Mult(x,y); }
+
    void MultTranspose(const Vector &x, Vector &y) const override;
+
+   void AbsMultTranspose(const Vector &x, Vector &y) const override { MultTranspose(x,y); }
 };
 
 }

fem/quadinterpolator.cpp

@@ -497,6 +497,25 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
                                   Vector &q_val,
                                   Vector &q_der,
                                   Vector &q_det) const
+{
+   MultInternal(e_vec, eval_flags, q_val, q_der, q_det);
+}
+
+void QuadratureInterpolator::AbsMult(const Vector &e_vec,
+                                     unsigned eval_flags,
+                                     Vector &q_val,
+                                     Vector &q_der,
+                                     Vector &q_det) const
+{
+   MultInternal(e_vec, eval_flags, q_val, q_der, q_det, true);
+}
+
+void QuadratureInterpolator::MultInternal(const Vector &e_vec,
+                                          unsigned eval_flags,
+                                          Vector &q_val,
+                                          Vector &q_der,
+                                          Vector &q_det,
+                                          bool ABS) const
 {
    using namespace internal::quadrature_interpolator;
 
@@ -528,12 +547,36 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
    MFEM_ASSERT(fespace->GetMesh()->GetNumGeometries(
                   fespace->GetMesh()->Dimension()) == 1,
                "mixed meshes are not supported");
+   MFEM_ASSERT(ABS?use_tensor_eval:true,
+               "AbsMult only implemented for tensor elements!");
+
+   // Create abs_maps, make B, Bt, G, Gt positive
+   DofToQuad abs_maps;
+   if (ABS)
+   {
+      abs_maps.FE = maps.FE;
+      abs_maps.IntRule = maps.IntRule;
+      abs_maps.mode = maps.mode;
+      abs_maps.ndof = maps.ndof;
+      abs_maps.nqpt = maps.nqpt;
+
+      abs_maps.B = maps.B;
+      abs_maps.Bt = maps.Bt;
+      abs_maps.G = maps.G;
+      abs_maps.Gt = maps.Gt;
+      auto abs_val = static_cast<real_t(*)(real_t)>(std::abs);
+      abs_maps.B.Apply(abs_val);
+      abs_maps.G.Apply(abs_val);
+      abs_maps.Bt.Apply(abs_val);
+      abs_maps.Gt.Apply(abs_val);
+   }
+   const DofToQuad &maps_ = ABS ? abs_maps : maps;
 
    if (use_tensor_eval)
    {
       if (eval_flags & VALUES)
       {
-         TensorEvalKernels::Run(dim, q_layout, vdim, nd, nq, ne, maps.B.Read(),
+         TensorEvalKernels::Run(dim, q_layout, vdim, nd, nq, ne, maps_.B.Read(),
                                 e_vec.Read(), q_val.Write(), vdim, nd, nq);
       }
       if (eval_flags & (DERIVATIVES | PHYSICAL_DERIVATIVES))
@@ -542,20 +585,20 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
          const real_t *J = phys ? geom->J.Read() : nullptr;
          const int s_dim = phys ? sdim : dim;
          GradKernels::Run(dim, q_layout, phys, vdim, nd, nq, ne,
-                          maps.B.Read(), maps.G.Read(), J, e_vec.Read(),
+                          maps_.B.Read(), maps_.G.Read(), J, e_vec.Read(),
                           q_der.Write(), s_dim, vdim, nd, nq);
       }
       if (eval_flags & DETERMINANTS)
       {
-         DetKernels::Run(dim, vdim, nd, nq, ne, maps.B.Read(),
-                         maps.G.Read(), e_vec.Read(), q_det.Write(), nd,
+         DetKernels::Run(dim, vdim, nd, nq, ne, maps_.B.Read(),
+                         maps_.G.Read(), e_vec.Read(), q_det.Write(), nd,
                          nq, &d_buffer);
       }
    }
    else // use_tensor_eval == false
    {
       EvalKernels::Run(dim, vdim, maps.ndof, maps.nqpt, ne,vdim,q_layout,
-                       geom, maps,e_vec, q_val,q_der,q_det,eval_flags);
+                       geom, maps_,e_vec, q_val,q_der,q_det,eval_flags);
    }
 }
 
@@ -592,6 +635,13 @@ void QuadratureInterpolator::PhysDerivatives(const Vector &e_vec,
    Mult(e_vec, PHYSICAL_DERIVATIVES, empty, q_der, empty);
 }
 
+void QuadratureInterpolator::AbsPhysDerivatives(const Vector &e_vec,
+                                                Vector &q_der) const
+{
+   Vector empty;
+   AbsMult(e_vec, PHYSICAL_DERIVATIVES, empty, q_der, empty);
+}
+
 void QuadratureInterpolator::Determinants(const Vector &e_vec,
                                           Vector &q_det) const
 {

fem/quadinterpolator.hpp

@@ -113,6 +113,13 @@ public:
    void Mult(const Vector &e_vec, unsigned eval_flags,
              Vector &q_val, Vector &q_der, Vector &q_det) const;
 
+   void AbsMult(const Vector &e_vec, unsigned eval_flags,
+                Vector &q_val, Vector &q_der, Vector &q_det) const;
+
+   void MultInternal(const Vector &e_vec, unsigned eval_flags,
+                     Vector &q_val, Vector &q_der, Vector &q_det,
+                     bool ABS = false) const;
+
    /// Interpolate the values of the E-vector @a e_vec at quadrature points.
    void Values(const Vector &e_vec, Vector &q_val) const;
 
@@ -124,6 +131,8 @@ public:
        @a e_vec at quadrature points. */
    void PhysDerivatives(const Vector &e_vec, Vector &q_der) const;
 
+   void AbsPhysDerivatives(const Vector &e_vec, Vector &q_der) const;
+
    /** @brief Compute the determinants of the derivatives (with respect to
        reference coordinates) of the E-vector @a e_vec at quadrature points. */
    void Determinants(const Vector &e_vec, Vector &q_det) const;

fem/restriction.cpp

@@ -678,6 +678,29 @@ void ConformingFaceRestriction::Mult(const Vector& x, Vector& y) const
    });
 }
 
+void ConformingFaceRestriction::MultUnsigned(const Vector& x, Vector& y) const
+{
+   if (nf==0) { return; }
+   // Assumes all elements have the same number of dofs
+   const int nface_dofs = face_dofs;
+   const int vd = vdim;
+   const bool t = byvdim;
+   auto d_indices = scatter_indices.Read();
+   auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
+   auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
+   mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
+   {
+      const int s_idx = d_indices[i];
+      const int idx = (s_idx >= 0) ? s_idx : -1 - s_idx;
+      const int dof = i % nface_dofs;
+      const int face = i / nface_dofs;
+      for (int c = 0; c < vd; ++c)
+      {
+         d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
+      }
+   });
+}
+
 static void ConformingFaceRestriction_AddMultTranspose(
    const int ndofs,
    const int face_dofs,

fem/restriction.hpp

@@ -65,9 +65,14 @@ public:
 
    /// Compute Mult without applying signs based on DOF orientations.
    void MultUnsigned(const Vector &x, Vector &y) const;
+
+   void AbsMult(const Vector &x, Vector &y) const override { MultUnsigned(x,y); };
+
    /// Compute MultTranspose without applying signs based on DOF orientations.
    void MultTransposeUnsigned(const Vector &x, Vector &y) const;
 
+   void AbsMultTranspose(const Vector &x, Vector &y) const override { MultTransposeUnsigned(x,y); };
+
    /// Compute MultTranspose by setting (rather than adding) element
    /// contributions; this is a left inverse of the Mult() operation
    void MultLeftInverse(const Vector &x, Vector &y) const;
@@ -176,6 +181,13 @@ public:
    */
    void Mult(const Vector &x, Vector &y) const override = 0;
 
+   virtual void MultUnsigned(const Vector &x, Vector &y) const
+   {
+      MFEM_ABORT("MultUnsigned not implemented yet...");
+   }
+
+   void AbsMult(const Vector &x, Vector &y) const override { MultUnsigned(x,y); }
+
    /** @brief Add the face degrees of freedom @a x to the element degrees of
        freedom @a y.
 
@@ -224,6 +236,12 @@ public:
       AddMultTranspose(x, y);
    }
 
+   void AbsMultTranspose(const Vector &x, Vector &y) const override
+   {
+      y = 0.0;
+      AddMultTransposeUnsigned(x,y);
+   }
+
    /** @brief For each face, sets @a y to the partial derivative of @a x with
               respect to the reference coordinate whose direction is
               perpendicular to the face on the reference element.
@@ -320,6 +338,8 @@ public:
                      ElementDofOrdering. */
    void Mult(const Vector &x, Vector &y) const override;
 
+   void MultUnsigned(const Vector &x, Vector &y) const override;
+
    using FaceRestriction::AddMultTransposeInPlace;
 
    /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
@@ -343,6 +363,19 @@ public:
    void AddMultTransposeUnsigned(const Vector &x, Vector &y,
                                  const real_t a = 1.0) const override;
 
+   using Operator::AddAbsMultTranspose;
+
+   void AddAbsMultTranspose(const Vector &x, Vector &y) const
+   {
+      AddMultTransposeUnsigned(x,y);
+   }
+
+   void AbsMultTranspose(const Vector &x, Vector &y) const override
+   {
+      y = 0.0;
+      AddMultTransposeUnsigned(x,y);
+   }
+
 private:
    /** @brief Compute the scatter indices: L-vector to E-vector, and the offsets
        for the gathering: E-vector to L-vector.

general/array.cpp

@@ -140,6 +140,15 @@ int Array<T>::IsSorted() const
    return 1;
 }
 
+template <class T>
+void Array<T>::Apply(std::function<T(T)> function)
+{
+   static_assert(std::is_arithmetic<T>::value, "Apply to arithmetric types!");
+   const bool use_device = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_device);
+   mfem::forall_switch(use_device, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = function(y[i]); });
+}
 
 template <class T>
 void Array2D<T>::Load(const char *filename, int fmt)
@@ -175,6 +184,7 @@ void Array2D<T>::Print(std::ostream &os, int width_)
    }
 }
 
+
 template class Array<char>;
 template class Array<int>;
 template class Array<long long>;

general/array.hpp

@@ -22,6 +22,7 @@
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
+#include <functional>
 #include <type_traits>
 #include <initializer_list>
 
@@ -104,6 +105,9 @@ public:
    /// Destructor
    inline ~Array() { data.Delete(); }
 
+   /// Apply function
+   void Apply(std::function<T(T)> function);
+
    /// Assignment operator: deep copy from 'src'.
    Array<T> &operator=(const Array<T> &src) { src.Copy(*this); return *this; }
 

linalg/densemat.cpp

@@ -149,6 +149,13 @@ void DenseMatrix::Mult(const Vector &x, Vector &y) const
    Mult(x.GetData(), y.GetData());
 }
 
+void DenseMatrix::PowAbsMult(const real_t p, const Vector &x, Vector &y) const
+{
+   MFEM_ASSERT(height == y.Size() && width == x.Size(),
+               "incompatible dimensions");
+   kernels::PowAbsMult(height, width, p, Data(), x.GetData(), y.GetData());
+}
+
 real_t DenseMatrix::operator *(const DenseMatrix &m) const
 {
    MFEM_ASSERT(Height() == m.Height() && Width() == m.Width(),

linalg/densemat.hpp

@@ -153,6 +153,9 @@ public:
    /// Matrix vector multiplication.
    void Mult(const Vector &x, Vector &y) const override;
 
+   /// Power absolute-value matrix vector multiplication.
+   virtual void PowAbsMult(const real_t p, const Vector &x, Vector &y) const;
+
    /// Multiply a vector with the transpose matrix.
    void MultTranspose(const real_t *x, real_t *y) const;
 

linalg/hypre.cpp

@@ -2037,6 +2037,58 @@ void HypreParMatrix::AbsMultTranspose(real_t a, const Vector &x,
    HypreRead();
 }
 
+void HypreParMatrix::PowAbsMult(real_t p, real_t a, const Vector &x,
+                                real_t b, Vector &y) const
+{
+   MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
+               << ", expected size = " << Width());
+   MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
+               << ", expected size = " << Height());
+   MFEM_ASSERT(p > 0.0, "Non-positive powers not implemented!");
+
+   auto x_data = x.HostRead();
+   auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
+
+   HostRead();
+   if (p == 1.0)
+   {
+      internal::hypre_ParCSRMatrixAbsMatvec(A, a, const_cast<real_t*>(x_data),
+                                            b, y_data);
+   }
+   else
+   {
+      internal::hypre_ParCSRMatrixPowAbsMatvec(A, p, a, const_cast<real_t*>(x_data),
+                                               b, y_data);
+   }
+   HypreRead();
+}
+
+void HypreParMatrix::PowAbsMultTranspose(real_t p, real_t a, const Vector &x,
+                                         real_t b, Vector &y) const
+{
+   MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
+               << ", expected size = " << Height());
+   MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
+               << ", expected size = " << Width());
+   MFEM_ASSERT(p > 0.0, "Non-positive powers not implemented!");
+
+   auto x_data = x.HostRead();
+   auto y_data = (b == 0.0) ? y.HostWrite() : y.HostReadWrite();
+
+   HostRead();
+   if (p == 1.0)
+   {
+      internal::hypre_ParCSRMatrixAbsMatvecT(A, a, const_cast<real_t*>(x_data),
+                                             b, y_data);
+   }
+   else
+   {
+      internal::hypre_ParCSRMatrixPowAbsMatvecT(A, p, a, const_cast<real_t*>(x_data),
+                                                b, y_data);
+   }
+   HypreRead();
+}
+
 HypreParMatrix* HypreParMatrix::LeftDiagMult(const SparseMatrix &D,
                                              HYPRE_BigInt* row_starts) const
 {

linalg/hypre.hpp

@@ -767,6 +767,7 @@ public:
 
    void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const override
    { Mult(a, x, 1.0, y); }
+
    void AddMultTranspose(const Vector &x, Vector &y,
                          const real_t a = 1.0) const override
    { MultTranspose(a, x, 1.0, y); }
@@ -778,10 +779,32 @@ public:
        of the matrix A. */
    void AbsMult(real_t a, const Vector &x, real_t b, Vector &y) const;
 
+   void AbsMult(const Vector &x, Vector &y) const override
+   { AbsMult(1.0, x, 0.0, y); }
+
+   void AddAbsMult(const Vector &x, Vector &y, const real_t a = 1.0) const override
+   { AbsMult(a, x, 1.0, y); }
+
    /** @brief Computes y = a * |At| * x + b * y, using entry-wise absolute
        values of the transpose of the matrix A. */
    void AbsMultTranspose(real_t a, const Vector &x, real_t b, Vector &y) const;
 
+   void AbsMultTranspose(const Vector &x, Vector &y) const override
+   { AbsMultTranspose(1.0, x, 0.0, y); }
+
+   void AddAbsMultTranspose(const Vector &x, Vector &y,
+                            const real_t a) const override
+   { AbsMultTranspose(a, x, 1.0, y); }
+
+   /** @brief Computes y = a * |A|**p * x + b * y, using entry-wise absolute values
+       of the matrix A. */
+   void PowAbsMult(real_t p, real_t a, const Vector &x, real_t b, Vector &y) const;
+
+   /** @brief Computes y = a * |At|**p * x + b * y, using entry-wise absolute
+       values of the transpose of the matrix A. */
+   void PowAbsMultTranspose(real_t p, real_t a, const Vector &x, real_t b,
+                            Vector &y) const;
+
    /** @brief The "Boolean" analog of y = alpha * A * x + beta * y, where
        elements in the sparsity pattern of the matrix are treated as "true". */
    void BooleanMult(int alpha, const int *x, int beta, int *y)

linalg/hypre_parcsr.cpp

@@ -1207,6 +1207,196 @@ void hypre_CSRMatrixAbsMatvecT(hypre_CSRMatrix *A,
    }
 }
 
+void hypre_CSRMatrixPowAbsMatvec(hypre_CSRMatrix *A,
+                                 HYPRE_Real p,
+                                 HYPRE_Real alpha,
+                                 HYPRE_Real *x,
+                                 HYPRE_Real beta,
+                                 HYPRE_Real *y)
+{
+   HYPRE_Real       *A_data   = hypre_CSRMatrixData(A);
+   HYPRE_Int        *A_i      = hypre_CSRMatrixI(A);
+   HYPRE_Int        *A_j      = hypre_CSRMatrixJ(A);
+   HYPRE_Int         num_rows = hypre_CSRMatrixNumRows(A);
+
+   HYPRE_Int        *A_rownnz = hypre_CSRMatrixRownnz(A);
+   HYPRE_Int         num_rownnz = hypre_CSRMatrixNumRownnz(A);
+
+   HYPRE_Real       *x_data = x;
+   HYPRE_Real       *y_data = y;
+
+   HYPRE_Real        temp, tempx;
+
+   HYPRE_Int         i, jj;
+
+   HYPRE_Int         m;
+
+   HYPRE_Real        xpar=0.7;
+
+   /*-----------------------------------------------------------------------
+    * Do (alpha == 0.0) computation - RDF: USE MACHINE EPS
+    *-----------------------------------------------------------------------*/
+
+   if (alpha == 0.0)
+   {
+      for (i = 0; i < num_rows; i++)
+      {
+         y_data[i] *= beta;
+      }
+      return;
+   }
+
+   /*-----------------------------------------------------------------------
+    * y = (beta/alpha)*y
+    *-----------------------------------------------------------------------*/
+
+   temp = beta / alpha;
+
+   if (temp != 1.0)
+   {
+      if (temp == 0.0)
+      {
+         for (i = 0; i < num_rows; i++)
+         {
+            y_data[i] = 0.0;
+         }
+      }
+      else
+      {
+         for (i = 0; i < num_rows; i++)
+         {
+            y_data[i] *= temp;
+         }
+      }
+   }
+
+   /*-----------------------------------------------------------------
+    * y += abs(A)*x
+    *-----------------------------------------------------------------*/
+
+   /* use rownnz pointer to do the abs(A)*x multiplication
+      when num_rownnz is smaller than num_rows */
+
+   if (num_rownnz < xpar*(num_rows))
+   {
+      for (i = 0; i < num_rownnz; i++)
+      {
+         m = A_rownnz[i];
+
+         tempx = 0;
+         for (jj = A_i[m]; jj < A_i[m+1]; jj++)
+         {
+            tempx += std::pow(std::abs(A_data[jj]), p)*x_data[A_j[jj]];
+         }
+         y_data[m] += tempx;
+      }
+   }
+   else
+   {
+      for (i = 0; i < num_rows; i++)
+      {
+         tempx = 0;
+         for (jj = A_i[i]; jj < A_i[i+1]; jj++)
+         {
+            tempx += std::pow(std::abs(A_data[jj]), p)*x_data[A_j[jj]];
+         }
+         y_data[i] += tempx;
+      }
+   }
+
+   /*-----------------------------------------------------------------
+    * y = alpha*y
+    *-----------------------------------------------------------------*/
+
+   if (alpha != 1.0)
+   {
+      for (i = 0; i < num_rows; i++)
+      {
+         y_data[i] *= alpha;
+      }
+   }
+
+}
+
+void hypre_CSRMatrixPowAbsMatvecT(hypre_CSRMatrix *A,
+                                  HYPRE_Real p,
+                                  HYPRE_Real alpha,
+                                  HYPRE_Real *x,
+                                  HYPRE_Real beta,
+                                  HYPRE_Real *y)
+{
+   HYPRE_Real       *A_data   = hypre_CSRMatrixData(A);
+   HYPRE_Int        *A_i      = hypre_CSRMatrixI(A);
+   HYPRE_Int        *A_j      = hypre_CSRMatrixJ(A);
+   HYPRE_Int         num_rows = hypre_CSRMatrixNumRows(A);
+   HYPRE_Int         num_cols = hypre_CSRMatrixNumCols(A);
+
+   HYPRE_Real       *x_data = x;
+   HYPRE_Real       *y_data = y;
+
+   HYPRE_Int         i, j, jj;
+
+   HYPRE_Real        temp;
+
+   if (alpha == 0.0)
+   {
+      for (i = 0; i < num_cols; i++)
+      {
+         y_data[i] *= beta;
+      }
+      return;
+   }
+
+   /*-----------------------------------------------------------------------
+    * y = (beta/alpha)*y
+    *-----------------------------------------------------------------------*/
+
+   temp = beta / alpha;
+
+   if (temp != 1.0)
+   {
+      if (temp == 0.0)
+      {
+         for (i = 0; i < num_cols; i++)
+         {
+            y_data[i] = 0.0;
+         }
+      }
+      else
+      {
+         for (i = 0; i < num_cols; i++)
+         {
+            y_data[i] *= temp;
+         }
+      }
+   }
+
+   /*-----------------------------------------------------------------
+    * y += abs(A)^T*x
+    *-----------------------------------------------------------------*/
+
+   for (i = 0; i < num_rows; i++)
+   {
+      for (jj = A_i[i]; jj < A_i[i+1]; jj++)
+      {
+         j = A_j[jj];
+         y_data[j] += std::pow(std::abs(A_data[jj]), p) * x_data[i];
+      }
+   }
+
+   /*-----------------------------------------------------------------
+    * y = alpha*y
+    *-----------------------------------------------------------------*/
+
+   if (alpha != 1.0)
+   {
+      for (i = 0; i < num_cols; i++)
+      {
+         y_data[i] *= alpha;
+      }
+   }
+}
+
 /* Based on hypre_CSRMatrixMatvec in hypre's csr_matvec.c */
 void hypre_CSRMatrixBooleanMatvec(hypre_CSRMatrix *A,
                                   HYPRE_Bool alpha,
@@ -1607,6 +1797,148 @@ void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
    mfem_hypre_TFree_host(y_tmp);
 }
 
+void hypre_ParCSRMatrixPowAbsMatvec(hypre_ParCSRMatrix *A,
+                                    HYPRE_Real p,
+                                    HYPRE_Real alpha,
+                                    HYPRE_Real *x,
+                                    HYPRE_Real beta,
+                                    HYPRE_Real *y)
+{
+   hypre_ParCSRCommHandle *comm_handle;
+   hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
+   hypre_CSRMatrix   *diag   = hypre_ParCSRMatrixDiag(A);
+   hypre_CSRMatrix   *offd   = hypre_ParCSRMatrixOffd(A);
+
+   HYPRE_Int          num_cols_offd = hypre_CSRMatrixNumCols(offd);
+   HYPRE_Int          num_sends, i, j, index;
+
+   HYPRE_Real        *x_tmp, *x_buf;
+
+   x_tmp = mfem_hypre_CTAlloc_host(HYPRE_Real, num_cols_offd);
+
+   /*---------------------------------------------------------------------
+    * If there exists no CommPkg for A, a CommPkg is generated using
+    * equally load balanced partitionings
+    *--------------------------------------------------------------------*/
+   if (!comm_pkg)
+   {
+      hypre_MatvecCommPkgCreate(A);
+      comm_pkg = hypre_ParCSRMatrixCommPkg(A);
+   }
+
+   num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
+   x_buf = mfem_hypre_CTAlloc_host(
+              HYPRE_Real, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
+
+   index = 0;
+   for (i = 0; i < num_sends; i++)
+   {
+      j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
+      for ( ; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
+      {
+         x_buf[index++] = x[hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j)];
+      }
+   }
+
+   comm_handle = hypre_ParCSRCommHandleCreate(1, comm_pkg, x_buf, x_tmp);
+
+   hypre_CSRMatrixPowAbsMatvec(diag, p, alpha, x, beta, y);
+
+   hypre_ParCSRCommHandleDestroy(comm_handle);
+
+   if (num_cols_offd)
+   {
+      hypre_CSRMatrixPowAbsMatvec(offd, p, alpha, x_tmp, 1.0, y);
+   }
+
+   mfem_hypre_TFree_host(x_buf);
+   mfem_hypre_TFree_host(x_tmp);
+}
+
+/* Based on hypre_ParCSRMatrixMatvecT in par_csr_matvec.c */
+void hypre_ParCSRMatrixPowAbsMatvecT(hypre_ParCSRMatrix *A,
+                                     HYPRE_Real p,
+                                     HYPRE_Real alpha,
+                                     HYPRE_Real *x,
+                                     HYPRE_Real beta,
+                                     HYPRE_Real *y)
+{
+   hypre_ParCSRCommHandle *comm_handle;
+   hypre_ParCSRCommPkg    *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
+   hypre_CSRMatrix        *diag = hypre_ParCSRMatrixDiag(A);
+   hypre_CSRMatrix        *offd = hypre_ParCSRMatrixOffd(A);
+   HYPRE_Real             *y_tmp;
+   HYPRE_Real             *y_buf;
+
+   HYPRE_Int               num_cols_offd = hypre_CSRMatrixNumCols(offd);
+
+   HYPRE_Int               i, j, jj, end, num_sends;
+
+   y_tmp = mfem_hypre_TAlloc_host(HYPRE_Real, num_cols_offd);
+
+   /*---------------------------------------------------------------------
+    * If there exists no CommPkg for A, a CommPkg is generated using
+    * equally load balanced partitionings
+    *--------------------------------------------------------------------*/
+   if (!comm_pkg)
+   {
+      hypre_MatvecCommPkgCreate(A);
+      comm_pkg = hypre_ParCSRMatrixCommPkg(A);
+   }
+
+   num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
+   y_buf = mfem_hypre_CTAlloc_host(
+              HYPRE_Real, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
+
+   if (num_cols_offd)
+   {
+      // Disable the use of offdT for now, until we implement
+      // hypre_CSRMatrixAbsMatvec on device.
+#if MFEM_HYPRE_VERSION >= 21100 && 0
+      if (A->offdT)
+      {
+         // offdT is optional. Used only if it's present.
+         hypre_CSRMatrixPowAbsMatvec(A->offdT, p, alpha, x, 0., y_tmp);
+      }
+      else
+#endif
+      {
+         hypre_CSRMatrixPowAbsMatvecT(offd, p, alpha, x, 0., y_tmp);
+      }
+   }
+
+   comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg, y_tmp, y_buf);
+
+   // Disable the use of diagT for now, until we implement
+   // hypre_CSRMatrixAbsMatvec on device.
+#if MFEM_HYPRE_VERSION >= 21100 && 0
+   if (A->diagT)
+   {
+      // diagT is optional. Used only if it's present.
+      hypre_CSRMatrixPowAbsMatvec(A->diagT, p, alpha, x, beta, y);
+   }
+   else
+#endif
+   {
+      hypre_CSRMatrixPowAbsMatvecT(diag, p, alpha, x, beta, y);
+   }
+
+   hypre_ParCSRCommHandleDestroy(comm_handle);
+
+   for (i = 0; i < num_sends; i++)
+   {
+      end = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1);
+      for (j = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i); j < end; j++)
+      {
+         jj = hypre_ParCSRCommPkgSendMapElmt(comm_pkg, j);
+         y[jj] += y_buf[j];
+      }
+   }
+
+   mfem_hypre_TFree_host(y_buf);
+   mfem_hypre_TFree_host(y_tmp);
+}
+
 /* Based on hypre_ParCSRMatrixMatvec in par_csr_matvec.c */
 void hypre_ParCSRMatrixBooleanMatvec(hypre_ParCSRMatrix *A,
                                      HYPRE_Bool alpha,

linalg/hypre_parcsr.hpp

@@ -169,6 +169,38 @@ void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
                                   HYPRE_Real beta,
                                   HYPRE_Real *y);
 
+/// Computes y = alpha * |A|**p * x + beta * y, using entry-wise absolute values of matrix A
+void hypre_CSRMatrixPowAbsMatvec(hypre_CSRMatrix *A,
+                                 HYPRE_Real p,
+                                 HYPRE_Real alpha,
+                                 HYPRE_Real *x,
+                                 HYPRE_Real beta,
+                                 HYPRE_Real *y);
+
+/// Computes y = alpha * |At|**p * x + beta * y, using entry-wise absolute values of the transpose of matrix A
+void hypre_CSRMatrixPowAbsMatvecT(hypre_CSRMatrix *A,
+                                  HYPRE_Real p,
+                                  HYPRE_Real alpha,
+                                  HYPRE_Real *x,
+                                  HYPRE_Real beta,
+                                  HYPRE_Real *y);
+
+/// Computes y = alpha * |A|**p * x + beta * y, using entry-wise absolute values of matrix A
+void hypre_ParCSRMatrixPowAbsMatvec(hypre_ParCSRMatrix *A,
+                                    HYPRE_Real p,
+                                    HYPRE_Real alpha,
+                                    HYPRE_Real *x,
+                                    HYPRE_Real beta,
+                                    HYPRE_Real *y);
+
+/// Computes y = alpha * |At|**p * x + beta * y, using entry-wise absolute values of the transpose of matrix A
+void hypre_ParCSRMatrixPowAbsMatvecT(hypre_ParCSRMatrix *A,
+                                     HYPRE_Real p,
+                                     HYPRE_Real alpha,
+                                     HYPRE_Real *x,
+                                     HYPRE_Real beta,
+                                     HYPRE_Real *y);
+
 /** The "Boolean" analog of y = alpha * A * x + beta * y, where elements in the
     sparsity pattern of the CSR matrix A are treated as "true". */
 void hypre_CSRMatrixBooleanMatvec(hypre_CSRMatrix *A,

linalg/kernels.hpp

@@ -188,6 +188,54 @@ void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
    }
 }
 
+/** @brief Power absolute-value Matrix vector multiplication: y = |A|^p x,
+    where the matrix A is of size @a height x @a width with given @a data,
+    while @a x and @a y specify the data of the input and output vectors. */
+template<typename TA, typename TX, typename TY>
+MFEM_HOST_DEVICE inline
+void PowAbsMult(const int height, const int width, const real_t p,
+                const TA *data, const TX *x, TY *y)
+{
+   if (width == 0)
+   {
+      for (int row = 0; row < height; row++)
+      {
+         y[row] = 0.0;
+      }
+      return;
+   }
+   const TA *d_col = data;
+   TX x_col = x[0];
+   for (int row = 0; row < height; row++)
+   {
+      if (p == 1.0)
+      {
+         y[row] = x_col*std::abs(d_col[row]);
+      }
+      else
+      {
+         y[row] = x_col*std::pow(std::abs(d_col[row]),p);
+      }
+   }
+   d_col += height;
+   for (int col = 1; col < width; col++)
+   {
+      x_col = x[col];
+      for (int row = 0; row < height; row++)
+      {
+         if (p == 1.0)
+         {
+            y[row] += x_col*std::abs(d_col[row]);
+         }
+         else
+         {
+            y[row] += x_col*std::pow(std::abs(d_col[row]),p);
+         }
+      }
+      d_col += height;
+   }
+}
+
 /** @brief Matrix transpose vector multiplication: y = At x, where the matrix A
     is of size @a height x @a width with given @a data, while @a x and @a y
     specify the data of the input and output vectors. */
@@ -217,6 +265,42 @@ void MultTranspose(const int height, const int width, const TA *data,
    }
 }
 
+/** @brief Power absolute-value matrix transpose vector multiplication: y = |At|^p x,
+    where the matrix A is of size @a height x @a width with given @a data, while @a x
+    and @a y specify the data of the input and output vectors. */
+template<typename TA, typename TX, typename TY>
+MFEM_HOST_DEVICE inline
+void PowAbsMultTranspose(const int height, const int width, const real_t p,
+                         const TA *data, const TX *x, TY *y)
+{
+   if (height == 0)
+   {
+      for (int row = 0; row < width; row++)
+      {
+         y[row] = 0.0;
+      }
+      return;
+   }
+   TY *y_off = y;
+   for (int i = 0; i < width; ++i)
+   {
+      TY val = 0.0;
+      for (int j = 0; j < height; ++j)
+      {
+         if (p == 1.0)
+         {
+            val += x[j] * std::abs(data[i * height + j]);
+         }
+         else
+         {
+            val += x[j] * std::pow(std::abs(data[i * height + j],p));
+         }
+      }
+      *y_off = val;
+      y_off++;
+   }
+}
+
 /// Symmetrize a square matrix with given @a size and @a data: A -> (A+A^T)/2.
 template<typename T>
 MFEM_HOST_DEVICE inline

linalg/operator.cpp

@@ -63,6 +63,21 @@ void Operator::AddMultTranspose(const Vector &x, Vector &y,
    y.Add(a, z);
 }
 
+void Operator::AddAbsMult(const Vector &x, Vector &y, const real_t a) const
+{
+   mfem::Vector z(y.Size());
+   AbsMult(x, z);
+   y.Add(a, z);
+}
+
+void Operator::AddAbsMultTranspose(const Vector &x, Vector &y,
+                                   const real_t a) const
+{
+   mfem::Vector z(y.Size());
+   AbsMultTranspose(x, z);
+   y.Add(a, z);
+}
+
 void Operator::ArrayMult(const Array<const Vector *> &X,
                          Array<Vector *> &Y) const
 {
@@ -645,18 +660,92 @@ void ConstrainedOperator::ConstrainedMult(const Vector &x, Vector &y,
    }
 }
 
+void ConstrainedOperator::ConstrainedAbsMult(const Vector &x, Vector &y,
+                                             const bool transpose) const
+{
+   const int csz = constraint_list.Size();
+   if (csz == 0)
+   {
+      if (transpose)
+      {
+         A->AbsMultTranspose(x, y);
+      }
+      else
+      {
+         A->AbsMult(x, y);
+      }
+      return;
+   }
+
+   z = x;
+
+   auto idx = constraint_list.Read();
+   // Use read+write access - we are modifying sub-vector of z
+   auto d_z = z.ReadWrite();
+   mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i) { d_z[idx[i]] = 0.0; });
+
+   if (transpose)
+   {
+      A->AbsMultTranspose(z, y);
+   }
+   else
+   {
+      A->AbsMult(z, y);
+   }
+
+   auto d_x = x.Read();
+   // Use read+write access - we are modifying sub-vector of y
+   auto d_y = y.ReadWrite();
+   switch (diag_policy)
+   {
+      case DIAG_ONE:
+         mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
+         {
+            const int id = idx[i];
+            d_y[id] = d_x[id];
+         });
+         break;
+      case DIAG_ZERO:
+         mfem::forall(csz, [=] MFEM_HOST_DEVICE (int i)
+         {
+            const int id = idx[i];
+            d_y[id] = 0.0;
+         });
+         break;
+      case DIAG_KEEP:
+         // Needs action of the operator diagonal on vector
+         mfem_error("ConstrainedOperator::AbsMult #1");
+         break;
+      default:
+         mfem_error("ConstrainedOperator::AbsMult #2");
+         break;
+   }
+}
+
 void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
 {
    constexpr bool transpose = false;
    ConstrainedMult(x, y, transpose);
 }
 
+void ConstrainedOperator::AbsMult(const Vector &x, Vector &y) const
+{
+   constexpr bool transpose = false;
+   ConstrainedAbsMult(x, y, transpose);
+}
+
 void ConstrainedOperator::MultTranspose(const Vector &x, Vector &y) const
 {
    constexpr bool transpose = true;
    ConstrainedMult(x, y, transpose);
 }
 
+void ConstrainedOperator::AbsMultTranspose(const Vector &x, Vector &y) const
+{
+   constexpr bool transpose = true;
+   ConstrainedAbsMult(x, y, transpose);
+}
+
 void ConstrainedOperator::AddMult(const Vector &x, Vector &y,
                                   const real_t a) const
 {
@@ -664,6 +753,13 @@ void ConstrainedOperator::AddMult(const Vector &x, Vector &y,
    y.Add(a, w);
 }
 
+void ConstrainedOperator::AddAbsMult(const Vector &x, Vector &y,
+                                     const real_t a) const
+{
+   AbsMult(x, w);
+   y.Add(a, w);
+}
+
 RectangularConstrainedOperator::RectangularConstrainedOperator(
    Operator *A,
    const Array<int> &trial_list,

linalg/operator.hpp

@@ -88,11 +88,23 @@ public:
    /// Operator application: `y=A(x)`.
    virtual void Mult(const Vector &x, Vector &y) const = 0;
 
+   /** @brief Action of the absolute-value operator: `y=|A|(x)`. The default
+       behavior in class Operator is to generate an error. If the Operator is a
+       composition of several operators, the composition unfold into a product
+       of absolute-value operators too. */
+   virtual void AbsMult(const Vector &x, Vector &y) const
+   { mfem_error("Operator::AbsMult() is not overridden!"); }
+
    /** @brief Action of the transpose operator: `y=A^t(x)`. The default behavior
        in class Operator is to generate an error. */
    virtual void MultTranspose(const Vector &x, Vector &y) const
    { mfem_error("Operator::MultTranspose() is not overridden!"); }
 
+   /** @brief Action of the transpose absolute-value operator: `y=|A|^t(x)`.
+      The default behavior in class Operator is to generate an error. */
+   virtual void AbsMultTranspose(const Vector &x, Vector &y) const
+   { mfem_error("Operator::AbsMultTranspose() is not overridden!"); }
+
    /// Operator application: `y+=A(x)` (default) or `y+=a*A(x)`.
    virtual void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
 
@@ -100,6 +112,13 @@ public:
    virtual void AddMultTranspose(const Vector &x, Vector &y,
                                  const real_t a = 1.0) const;
 
+   /// Operator application: `y+=|A|(x)` (default) or `y+=a*|A|(x)`.
+   virtual void AddAbsMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
+
+   /// Operator transpose application: `y+=|A|^t(x)` (default) or `y+=a*|A|^t(x)`.
+   virtual void AddAbsMultTranspose(const Vector &x, Vector &y,
+                                    const real_t a = 1.0) const;
+
    /// Operator application on a matrix: `Y=A(X)`.
    virtual void ArrayMult(const Array<const Vector *> &X,
                           Array<Vector *> &Y) const;
@@ -930,6 +949,10 @@ public:
    void Mult(const Vector & x, Vector & y) const override
    { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
 
+   /// Operator-wise absolute-value application.
+   void AbsMult(const Vector & x, Vector & y) const override
+   { P.AbsMult(x, Px); A.AbsMult(Px, APx); Rt.AbsMultTranspose(APx, y); }
+
    /// Approximate diagonal of the RAP Operator.
    /** Returns the diagonal of A, as returned by its AssembleDiagonal method,
        multiplied be P^T.
@@ -950,6 +973,10 @@ public:
    /// Application of the transpose.
    void MultTranspose(const Vector & x, Vector & y) const override
    { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
+
+   /// Operator-wise absolute-value application of the transpose
+   void AbsMultTranspose(const Vector & x, Vector & y) const override
+   { Rt.AbsMult(x, APx); A.AbsMultTranspose(APx, Px); P.AbsMultTranspose(Px, y); }
 };
 
 
@@ -1043,15 +1070,27 @@ public:
        the vectors, and "_i" -- the rest of the entries. */
    void Mult(const Vector &x, Vector &y) const override;
 
+   void AbsMult(const Vector &x, Vector &y) const override;
+
    void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const override;
 
+   void AddAbsMult(const Vector &x, Vector &y,
+                   const real_t a = 1.0) const override;
+
    void MultTranspose(const Vector &x, Vector &y) const override;
 
+   void AbsMultTranspose(const Vector &x, Vector &y) const override;
+
    /** @brief Implementation of Mult or MultTranspose.
     *  TODO - Generalize to allow constraining rows and columns differently.
    */
    void ConstrainedMult(const Vector &x, Vector &y, const bool transpose) const;
 
+   /** @brief Implementation of AbsMult or AbsMultTranspose.
+    *  TODO - Generalize to allow constraining rows and columns differently.
+   */
+   void ConstrainedAbsMult(const Vector &x, Vector &y, const bool transpose) const;
+
    /// Destructor: destroys the unconstrained Operator, if owned.
    ~ConstrainedOperator() override { if (own_A) { delete A; } }
 };

linalg/solvers.cpp

@@ -329,6 +329,128 @@ void OperatorJacobiSmoother::Mult(const Vector &x, Vector &y) const
    });
 }
 
+#ifdef MFEM_USE_MPI
+OperatorLpqJacobiSmoother::OperatorLpqJacobiSmoother(const real_t p_order,
+                                                     const real_t q_order,
+                                                     const real_t dmpng)
+   : damping(dmpng),
+     p_order(p_order),
+     q_order(q_order),
+     ess_tdof_list(nullptr),
+     oper(nullptr)
+{ }
+
+OperatorLpqJacobiSmoother::OperatorLpqJacobiSmoother(const HypreParMatrix &A,
+                                                     const Array<int> &ess_tdofs,
+                                                     const real_t p_order,
+                                                     const real_t q_order,
+                                                     const real_t dmpng)
+   : Solver(A.Height()),
+     dinv(height),
+     damping(dmpng),
+     p_order(p_order),
+     q_order(q_order),
+     ess_tdof_list(&ess_tdofs),
+     residual(height)
+{
+
+   oper = &A;
+   Setup();
+}
+
+void OperatorLpqJacobiSmoother::Mult(const Vector &x, Vector &y) const
+{
+   // For empty MPI ranks, height may be 0:
+   // MFEM_VERIFY(Height() > 0, "The diagonal hasn't been computed.");
+   MFEM_ASSERT(x.Size() == Width(), "invalid input vector");
+   MFEM_ASSERT(y.Size() == Height(), "invalid output vector");
+
+   if (iterative_mode)
+   {
+      MFEM_VERIFY(oper, "iterative_mode == true requires the forward operator");
+      oper->Mult(y, residual);  // r = A y
+      subtract(x, residual, residual); // r = x - A y
+   }
+   else
+   {
+      residual = x;
+      y.UseDevice(true);
+      y = 0.0;
+   }
+   auto DI = dinv.Read();
+   auto R = residual.Read();
+   auto Y = y.ReadWrite();
+   mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
+   {
+      Y[i] += DI[i] * R[i];
+   });
+}
+
+void OperatorLpqJacobiSmoother::Setup()
+{
+   auto hypre_mat = dynamic_cast<const HypreParMatrix *>(oper);
+   if (!hypre_mat) { MFEM_ABORT("Only implemented with HypreParMatrix!"); }
+
+   // Make vector
+   Vector right(hypre_mat->Height());
+   Vector temp(hypre_mat->Height());
+   Vector left(hypre_mat->Height());
+
+   // D^{-q} 1
+   right = 1.0;
+   if (q_order !=0.0)
+   {
+      hypre_mat->GetDiag(right);
+      right.PowerAbs(-q_order);
+   }
+
+   // |A|^p D^{-q} 1
+   temp = 0.0;
+   hypre_mat->PowAbsMult(p_order, 1.0, right, 0.0, temp);
+
+   // D^{1+q-p} |A|^p D^{-q} 1
+   left = temp;
+   if (1.0 + q_order - p_order != 0.0)
+   {
+      hypre_mat->GetDiag(left);
+      left.PowerAbs(1.0 + q_order - p_order);
+      left *= temp;
+   }
+
+   residual.UseDevice(true);
+   const real_t delta = damping;
+   auto D = left.Read();
+   auto DI = dinv.Write();
+
+   mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
+   {
+      if (D[i] == 0.0)
+      {
+         MFEM_ABORT_KERNEL("Zero diagonal entry in OperatorLpqJacobiSmoother");
+      }
+      DI[i] = delta / D[i];
+   });
+   if (ess_tdof_list && ess_tdof_list->Size() > 0)
+   {
+      auto I = ess_tdof_list->Read();
+      mfem::forall(ess_tdof_list->Size(), [=] MFEM_HOST_DEVICE (int i)
+      {
+         DI[I[i]] = delta;
+      });
+   }
+}
+
+real_t OperatorLpqJacobiSmoother::CheckSpectralBoundConstant()
+{
+   Vector ones(oper->Height());
+   Vector diag_comp(oper->Height());
+   ones = 1.0;
+   oper->AbsMult(ones, diag_comp);
+   diag_comp *= dinv;
+   return diag_comp.Max();
+}
+#endif // MFEM_USE_MPI
+
 OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
                                                      const Vector &d,
                                                      const Array<int>& ess_tdofs,
@@ -597,6 +719,7 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
       mfem::out << "   Iteration : " << setw(3) << right << 0 << "  ||Br|| = "
                 << nom << (print_options.first_and_last ? " ..." : "") << '\n';
    }
+   Monitor(0, nom, r, x);
 
    r0 = std::max(nom*rel_tol, abs_tol);
    if (nom <= r0)
@@ -645,17 +768,19 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
          done = true;
       }
 
-      if (++i > max_iter)
-      {
-         done = true;
-      }
-
       if (print_options.iterations || (done && print_options.first_and_last))
       {
          mfem::out << "   Iteration : " << setw(3) << right << (i-1)
                    << "  ||Br|| = " << setw(11) << left << nom
                    << "\tConv. rate: " << cf << '\n';
       }
+      Monitor(i, nom, r, x);
+
+      if (++i > max_iter)
+      {
+         done = true;
+      }
+
 
       if (done) { break; }
    }
@@ -673,6 +798,8 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
    }
 
    final_norm = nom;
+
+   Monitor(final_iter, final_norm, r, x, true);
 }
 
 void SLI(const Operator &A, const Vector &b, Vector &x,
@@ -750,7 +877,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       mfem::out << "   Iteration : " << setw(3) << 0 << "  (B r, r) = "
                 << nom << (print_options.first_and_last ? " ...\n" : "\n");
    }
-   Monitor(0, nom, r, x);
+   Monitor(0, sqrt(nom), r, x);
 
    if (nom < 0.0)
    {
@@ -830,7 +957,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
                    << betanom << std::endl;
       }
 
-      Monitor(i, betanom, r, x);
+      Monitor(i, sqrt(betanom), r, x);
 
       if (betanom <= r0)
       {

linalg/solvers.hpp

@@ -386,6 +386,43 @@ public:
    void Setup(const Vector &diag);
 };
 
+/// Generalized L(p,q)-Jacobi smoothing for bilinear form or a matrix.
+#ifdef MFEM_USE_MPI
+class OperatorLpqJacobiSmoother : public Solver
+{
+public:
+   OperatorLpqJacobiSmoother(const real_t p_order,
+                             const real_t q_order,
+                             const real_t dmpng = 1.0);
+
+   OperatorLpqJacobiSmoother(const HypreParMatrix &A,
+                             const Array<int> &ess_tdofs,
+                             const real_t p_order,
+                             const real_t q_order,
+                             const real_t dmpng =  1.0);
+
+   ~OperatorLpqJacobiSmoother() {}
+
+   void Mult(const Vector &x, Vector &y) const;
+
+   void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
+
+   void SetOperator(const Operator &op) { oper = &op; Setup(); };
+
+   real_t CheckSpectralBoundConstant();
+
+   void Setup();
+private:
+   Vector dinv;
+   const real_t damping;
+   const real_t p_order;
+   const real_t q_order;
+   const Array<int> *ess_tdof_list; // not owned; may be NULL
+   mutable Vector residual;
+   const Operator *oper; // not owned
+};
+#endif
+
 /// Chebyshev accelerated smoothing with given vector, no matrix necessary
 /** Potentially useful with tensorized operators, for example. This is just a
     very basic Chebyshev iteration, if you want tolerances, iteration control,

linalg/sparsemat.cpp

@@ -1224,6 +1224,107 @@ void SparseMatrix::AbsMultTranspose(const Vector &x, Vector &y) const
    }
 }
 
+void SparseMatrix::PowAbsMult(real_t p, const Vector &x, Vector &y) const
+{
+   MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
+               << ") must match matrix width (" << width << ")");
+   MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
+               << ") must match matrix height (" << height << ")");
+   MFEM_ASSERT(p > 0.0, "Non-positive powers not implemented!");
+
+   if (p == 1.0) { AbsMult(x,y); return; }
+
+   if (Finalized()) { y.UseDevice(true); }
+   y = 0.0;
+
+   if (!Finalized())
+   {
+      const real_t *xp = x.HostRead();
+      real_t *yp = y.HostReadWrite();
+
+      // The matrix is not finalized, but multiplication is still possible
+      for (int i = 0; i < height; i++)
+      {
+         RowNode *row = Rows[i];
+         real_t b = 0.0;
+         for ( ; row != NULL; row = row->Prev)
+         {
+            b += std::pow(std::abs(row->Value), p) * xp[row->Column];
+         }
+         *yp += b;
+         yp++;
+      }
+      return;
+   }
+
+   const int height = this->height;
+   const int nnz = J.Capacity();
+   auto d_I = Read(I, height+1);
+   auto d_J = Read(J, nnz);
+   auto d_A = Read(A, nnz);
+   auto d_x = x.Read();
+   auto d_y = y.ReadWrite();
+   mfem::forall(height, [=] MFEM_HOST_DEVICE (int i)
+   {
+      real_t d = 0.0;
+      const int end = d_I[i+1];
+      for (int j = d_I[i]; j < end; j++)
+      {
+         d += std::pow(std::abs(d_A[j]), p) * d_x[d_J[j]];
+      }
+      d_y[i] += d;
+   });
+}
+
+void SparseMatrix::PowAbsMultTranspose(real_t p, const Vector &x,
+                                       Vector &y) const
+{
+   MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
+               << ") must match matrix height (" << height << ")");
+   MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
+               << ") must match matrix width (" << width << ")");
+   MFEM_ASSERT(p > 0.0, "Non-positive powers not implemented!");
+
+   if (p == 1.0) { AbsMultTranspose(x,y); return; }
+
+   y = 0.0;
+
+   if (!Finalized())
+   {
+      real_t *yp = y.GetData();
+      // The matrix is not finalized, but multiplication is still possible
+      for (int i = 0; i < height; i++)
+      {
+         RowNode *row = Rows[i];
+         real_t b = x(i);
+         for ( ; row != NULL; row = row->Prev)
+         {
+            yp[row->Column] += std::pow(std::abs(row->Value), p) * b;
+         }
+      }
+      return;
+   }
+
+   EnsureMultTranspose();
+   if (At)
+   {
+      At->PowAbsMult(p, x, y);
+   }
+   else
+   {
+      for (int i = 0; i < height; i++)
+      {
+         const real_t xi = x[i];
+         const int end = I[i+1];
+         for (int j = I[i]; j < end; j++)
+         {
+            const int Jj = J[j];
+            y[Jj] += std::pow(std::abs(A[j]), p) * xi;
+         }
+      }
+   }
+}
+
 real_t SparseMatrix::InnerProduct(const Vector &x, const Vector &y) const
 {
    MFEM_ASSERT(x.Size() == Width(), "x.Size() = " << x.Size()

linalg/sparsemat.hpp

@@ -422,13 +422,22 @@ public:
    void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
 
    /// y = |A| * x, using entry-wise absolute values of matrix A
-   void AbsMult(const Vector &x, Vector &y) const;
+   void AbsMult(const Vector &x, Vector &y) const override;
 
    /// y = |At| * x, using entry-wise absolute values of the transpose of matrix A
    /** If the matrix is modified, call ResetTranspose() and optionally
        EnsureMultTranspose() to make sure this method uses the correct updated
        transpose. */
-   void AbsMultTranspose(const Vector &x, Vector &y) const;
+   void AbsMultTranspose(const Vector &x, Vector &y) const override;
+
+   /// y = |A|**p * x, using entry-wise absolute values of matrix A
+   void PowAbsMult(real_t p, const Vector &x, Vector &y) const;
+
+   /// y = |At|**p * x, using entry-wise absolute values of the transpose of matrix A
+   /** If the matrix is modified, call ResetTranspose() and optionally
+       EnsureMultTranspose() to make sure this method uses the correct updated
+       transpose. */
+   void PowAbsMultTranspose(real_t p, const Vector &x, Vector &y) const;
 
    /// Compute y^t A x
    real_t InnerProduct(const Vector &x, const Vector &y) const;

linalg/vector.cpp

@@ -310,6 +310,20 @@ void Vector::Reciprocal()
    mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = 1.0/y[i]; });
 }
 
+void Vector::PowerAbs(const real_t p)
+{
+   MFEM_ASSERT(p != 0.0, "requires p != 0.0");
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   if (p == 1.0)
+   {
+      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = std::abs(y[i]); });
+      return;
+   }
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = std::pow(std::abs(y[i]), p); });
+}
+
 void add(const Vector &v1, const Vector &v2, Vector &v)
 {
    MFEM_ASSERT(v.size == v1.size && v.size == v2.size,

linalg/vector.hpp

@@ -353,6 +353,9 @@ public:
    /// (*this)(i) = 1.0 / (*this)(i)
    void Reciprocal();
 
+   /// (*this)(i) = abs((*this)(i))^p
+   void PowerAbs(const real_t p);
+
    /// Swap the contents of two Vectors
    inline void Swap(Vector &other);
 

makefile

@@ -125,11 +125,11 @@ EXAMPLE_TEST_DIRS := examples
 
 MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff hooke \
- multidomain dpg hdiv-linear-solver spde
+ multidomain dpg hdiv-linear-solver spde diag-smoothers
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys shifted dpg)
+ toys shifted dpg diag-smoothers)
 
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 

miniapps/CMakeLists.txt

@@ -37,3 +37,4 @@ add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
+add_subdirectory(diag-smoothers)

miniapps/diag-smoothers/CMakeLists.txt

@@ -0,0 +1,131 @@
+# Copyright (c) 2010-2024, 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.
+
+if(MFEM_USE_MPI)
+  set(MESH_FILES_FROM_MESHING
+    icf.mesh
+    cube.mesh
+    # square01.mesh
+    # square01-tri.mesh
+    # stretched2D.mesh
+    # amr-quad-q2.mesh
+    # blade.mesh
+    # cube-tet.mesh
+  )
+  
+  set(MESH_FILES_FROM_GSLIB
+    triple-pt-1.mesh
+    triple-pt-2.mesh
+  )
+  
+  set(MESH_FILES_FROM_DATA
+    beam-tet.mesh
+    square-disc-p2.mesh
+    fichera-mixed-p2.mesh
+    amr-quad.mesh
+    ref-cube.mesh
+    ref-square.mesh
+  )
+  
+  # Add a target to copy the mesh files from the source directory; used by sample
+  # runs.
+  set(SRC_MESH_FILES)
+  foreach(MESH_FILE ${MESH_FILES_FROM_MESHING})
+    list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/../meshing/${MESH_FILE})
+  endforeach()
+  
+  foreach(MESH_FILE ${MESH_FILES_FROM_GSLIB})
+    list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/../gslib/${MESH_FILE})
+  endforeach()
+  
+  foreach(MESH_FILE ${MESH_FILES_FROM_DATA})
+    list(APPEND SRC_MESH_FILES ${CMAKE_CURRENT_SOURCE_DIR}/../../data/${MESH_FILE})
+  endforeach()
+  
+  add_custom_command(OUTPUT data_is_copied
+    COMMAND ${CMAKE_COMMAND} -E make_directory meshes
+    COMMAND ${CMAKE_COMMAND} -E copy_if_different ${SRC_MESH_FILES} ./meshes
+    COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
+    COMMENT "Copying meshing miniapps data files ...")
+  add_custom_target(copy_meshes DEPENDS data_is_copied)
+
+  add_custom_target(clean_meshes
+    COMMAND ${CMAKE_COMMAND} -E remove_directory meshes
+    COMMENT "Cleaning up meshes directory..."
+  )
+
+  add_custom_target(clean_common
+    COMMAND ${CMAKE_COMMAND} -E echo "Cleaning common directory..."
+    COMMAND make -C ${CMAKE_CURRENT_SOURCE_DIR}/../common clean
+    COMMENT "Cleaning common directory..."
+  )
+
+  set_property(DIRECTORY APPEND PROPERTY ADDITIONAL_MAKE_CLEAN_FILES meshes)
+  
+  # Define miniapps
+  set(MINIAPPS
+    lpq-jacobi
+    mg-lpq-jacobi
+    abs-l1-jacobi
+    mg-abs-l1-jacobi
+  )
+
+  # Add miniapps
+  foreach(APP ${MINIAPPS})
+    add_mfem_miniapp(${APP}
+      MAIN ${APP}.cpp
+      ${MFEM_MINIAPPS_COMMON_HEADERS}
+      EXTRA_SOURCES lpq-common.cpp
+      EXTRA_HEADERS lpq-common.hpp
+      LIBRARIES mfem-common mfem
+    )
+    add_dependencies(${APP} copy_meshes)
+  endforeach()
+
+  # Add the corresponding tests to the "test" target
+  if (MFEM_ENABLE_TESTING)
+    add_test(NAME lpq-jacobi_np${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:lpq-jacobi> -m meshes/cube.mesh -rs 1 -rp 1
+      -s 1 -i 1 -p 1.5 -q 0.75 -pc 1 -no-mon -no-vis
+      ${MPIEXEC_POSTFLAGS})
+
+    add_test(NAME abs-l1-jacobi_np${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:abs-l1-jacobi> -m meshes/cube.mesh -rs 1 -rp 1
+      -s 1 -i 1 -a 4 -pc 1 -no-mon -no-vis
+      ${MPIEXEC_POSTFLAGS})
+
+    add_test(NAME mg-lpq-jacobi_np${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:mg-lpq-jacobi> -m meshes/cube.mesh -rs 1 -rp 1
+      -ol 1 -gl 1 -s 1 -i 1 -p 1.5 -q 0.75 -no-mon -no-vis
+      ${MPIEXEC_POSTFLAGS})
+
+    add_test(NAME mg-abs-l1-jacobi_np${MFEM_MPI_NP}
+      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+      ${MPIEXEC_PREFLAGS}
+      $<TARGET_FILE:mg-abs-l1-jacobi> -m meshes/cube.mesh -rs 1 -rp 1
+      -ol 1 -gl 1 -s 1 -i 1 -a 4 -no-mon -no-vis
+      ${MPIEXEC_POSTFLAGS})
+  endif()
+
+  # Add custom target for cleaning everything
+  add_custom_target(clean_all
+    COMMAND ${CMAKE_MAKE_PROGRAM} clean
+    COMMAND ${CMAKE_COMMAND} -E remove_directory meshes
+    COMMAND ${CMAKE_MAKE_PROGRAM} -C ${CMAKE_CURRENT_SOURCE_DIR}/../common clean
+    COMMENT "Cleaning all build artifacts and meshes..."
+  )
+endif()

miniapps/diag-smoothers/README.md

@@ -0,0 +1,75 @@
+```
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               https://mfem.org
+```
+
+This directory contains some drivers reimplementing basic examples in MFEM,
+making use of the L(p,q)-Jacobi family of preconditioners/smoothers.
+
+Make sure you are familiar with the following examples:
+- `ex1p` [Laplace Problem](https://github.com/mfem/mfem/blob/master/examples/ex1p.cpp)
+- `ex2p` [Linear Elasticity](https://github.com/mfem/mfem/blob/master/examples/ex2p.cpp)
+- `ex3p` [Definite Maxwell Problem](https://github.com/mfem/mfem/blob/master/examples/ex3p.cpp)
+- `ex26p` [Multigrid Preconditioner](https://github.com/mfem/mfem/blob/master/examples/ex26p.cpp)
+
+The code has *four* drivers: `lpq-jacobi`, `mg-lpq-jacobi`, `abs-l1-jacobi`,
+and `mg-abs-l1-jacobi`. All these drivers have the capability to solve the following
+problems:
+- An L2-projection into a conforming H1-space.
+- A diffusion problem.
+- A linear elasticity problem.
+- A definite Maxwell probelm.
+
+For later reference, we say a smoother `M` is `A`-convergent if `M + M^T - A` is SPD,
+this is `(Ax,x) < (Mx, x) + (M^T x, x) = 2 (Mx,x)`. It suffices to find a constant
+`c < 2` such that `(Ax,x) < c(Mx,x)` to say that `M` is `A`-convergent.
+
+# L(p,q)-Jacobi preconditioners for fully assembled systems
+
+The driver `lpq-jacobi` allows the user solve the above mentioned problems,
+utilizing a fully assembled system. The matrix of the system is required to construct
+the L(p,q)-Jacobi preconditioners.
+
+The L(p,q)-Jacobi preconditioners can be described by `D_{p,q) = diag( D^{1 + q - p}
+|A|^{p} D^{-q} 1 )`, where `1` is the constant vector, `D` is the diagonal of `A`,
+and the operations are understood *entrywise*. *All of this matrices are SPD*.
+
+This code allows the uses to chose a problem (integrator) `mass, diffusion, elasticity, maxwell`,
+a mesh, a solver `cg, sli`, and a Kershaw transformation to define the system. The user can
+modify the frequency of the solution, the type of peconditioner `none, global, element`,
+the order of the preconditioner (in case of using L(p,q)-Jacobi) `p_order, q_order`, the
+polynomial degree of the underlying FES, the number of refinements (in serial and in parallel),
+tolerance and maximum number of iterations for the solver, compuation on device, get `.csv`
+outputs with the monitor option, and visualization to GLVis.
+
+# Absolute-value L(1)-Jacobi preconditioner for different assembly levels
+
+The driver `lpq-jacobi` allows the user solve the above mentioned problems,
+utilizing a different types of assembly levels. Our interest lies on `AssemblyLevel::PARTIAL`.
+As the FEM operator has the structure `A = P^T G^T B^T D B G P`
+(see [this](https://mfem.org/performance/)), and the standard L(1)-Jacobi can be writen as
+`D_1 = diag( |A|1 ). A triangle inequality implies
+`D_{abs} = diag( |P^T| |G^T| |B^T| |D| |B| |G| |P| 1 )` is also `A`-convergent.
+
+The MFEM interface allows to make use of `AbsMult` with the purpose of unwrap a composition
+(of different kind) of operators as their absolute-value application. Similar run-time
+options are available.
+
+# Multigrid wrappers
+
+The drivers `mg-lpq-jacobi` and `mg-abs-l1-jacobi` are basically the multigrid counterparts
+of the previously mentioned drivers. The wrapper (akin to
+[`ex26p`](https://github.com/mfem/mfem/blob/master/examples/ex26p.cpp))
+allow the user to do geometric refinement or order refinement.
+
+# Caveat: elementwise L(p,q)-Jacobi
+
+A small interface is implemented for performing an element-based L(p,q)-Jacobi.
+The function is defined in `lpq-common.xpp` and calls the element matrices,
+compute the element L(p,q)-preconditioner, and then accumulates them into one operator.

miniapps/diag-smoothers/abs-l1-jacobi.cpp

@@ -0,0 +1,450 @@
+// Copyright (c) 2010-2024, 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.
+//
+//                 -----------------------------------------
+//                   Absolute L(1)-Jacobi smoothers miniapp
+//                 -----------------------------------------
+//
+// This miniapp illustrates the implementation of an (slightly generalized)
+// absolute-L(1) Jacobi preconditioner. This preconditioner is tested in different
+// settings. We use Stationary Linear Iterations and Preconditioned Conjugate
+// Gradient as the main solvers.
+// We consider a H1-mass matrix, a diffusion matrix, a elasticity system, and a
+// definite Maxwell system.
+//
+// The preconditioner can be defined at run-time. Similarly, the mesh can be
+// modified by a Kershaw transformation at run-time. Relative tolerance and
+// maximum number of iterations can be modified as well.
+//
+// Compile with: make l1-partial
+//
+// Sample runs: mpirun -np 4 ./l1-partial
+//              mpirun -np 4 ./l1-partial -s 1 -i 3
+//              mpirun -np 4 ./l1-partial -m meshes/icf.mesh -f 0.5
+//              mpirun -np 4 ./l1-partial -rs 2 -rp 0
+//              mpirun -np 4 ./l1-partial -t 1e5 -ni 100 -vis
+//              mpirun -np 4 ./l1-partial -m meshes/beam-tet.mesh -Ky 0.5 -Kz 0.5
+
+#include "lpq-common.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace lpq_common;
+
+int main(int argc, char *argv[])
+{
+   /// 1. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   /// 2. Parse command line options.
+   string mesh_file = "meshes/cube.mesh";
+   // System properties
+   int order = 1;
+   SolverType solver_type = cg;
+   IntegratorType integrator_type = diffusion;
+   LpqType pc_type = global;
+   int assembly_type_int = 4;
+   AssemblyLevel assembly_type;
+   // Number of refinements
+   int refine_serial = 0;
+   int refine_parallel = 0;
+   // Preconditioner parameters
+   real_t p_order = 1.0;
+   real_t q_order = 0.0;
+   // Solver parameters
+   real_t rel_tol = 1e-10;
+   real_t max_iter = 3000;
+   // Kershaw Transformation
+   real_t eps_y = 0.0;
+   real_t eps_z = 0.0;
+   // Other options
+   string device_config = "cpu";
+   bool use_monitor = false;
+   bool visualization = false;
+
+   // Construct argument parser
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree)");
+   args.AddOption((int*)&solver_type, "-s", "--solver",
+                  "Solvers to be considered:"
+                  "\n\t0: Stationary Linear Iteration"
+                  "\n\t1: Preconditioned Conjugate Gradient");
+   args.AddOption((int*)&integrator_type, "-i", "--integrator",
+                  "Integrators to be considered:"
+                  "\n\t0: MassIntegrator"
+                  "\n\t1: DiffusionIntegrator"
+                  "\n\t2: ElasticityIntegrator"
+                  "\n\t3: CurlCurlIntegrator + VectorFEMassIntegrator");
+   args.AddOption(&assembly_type_int, "-a", "--assembly",
+                  "Assembly level to be considered:"
+                  "\n\t0: LEGACY"
+                  "\n\t1: LEGACYFULL (Deprecated)"
+                  "\n\t2: FULL"
+                  "\n\t3: ELEMENT"
+                  "\n\t4: PARTIAL"
+                  "\n\t5: NONE");
+   args.AddOption((int*)&pc_type, "-pc", "--preconditioner",
+                  "Preconditioners to be considered:"
+                  "\n\t0: No preconditioner"
+                  "\n\t1: L(p,q)-Jacobi preconditioner"
+                  "\n\t2: Element L(p,q)-Jacobi preconditioner");
+   args.AddOption(&refine_serial, "-rs", "--refine-serial",
+                  "Number of serial refinements");
+   args.AddOption(&refine_parallel, "-rp", "--refine-parallel",
+                  "Number of parallel refinements");
+   args.AddOption(&p_order, "-p", "--p-order",
+                  "P-order for L(p,q)-Jacobi preconditioner");
+   args.AddOption(&q_order, "-q", "--q-order",
+                  "Q-order for L(p,q)-Jacobi preconditioner");
+   args.AddOption(&rel_tol, "-t", "--tolerance",
+                  "Relative tolerance for the iterative solver");
+   args.AddOption(&max_iter, "-ni", "--iterations",
+                  "Maximum number of iterations");
+   args.AddOption(&eps_y, "-Ky", "--Kershaw-y",
+                  "Kershaw transform factor, eps_y in (0,1]");
+   args.AddOption(&eps_z, "-Kz", "--Kershaw-z",
+                  "Kershaw transform factor, eps_z in (0,1]");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&use_monitor, "-mon", "--monitor", "-no-mon",
+                  "--no-monitor",
+                  "Enable or disable Data Monitor.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   MFEM_ASSERT(p_order > 0.0, "p needs to be positive");
+   MFEM_ASSERT((0 <= solver_type) && (solver_type < num_solvers), "");
+   MFEM_ASSERT((0 <= integrator_type) && (integrator_type < num_integrators), "");
+   MFEM_ASSERT((0 <= assembly_type_int) && (assembly_type_int < 6), "");
+   MFEM_ASSERT((0 <= pc_type) && (pc_type < num_lpq_pc), "");
+   MFEM_ASSERT((0.0 <= eps_y) && (eps_y <= 1.0), "eps_y in [0,1]");
+   MFEM_ASSERT((0.0 <= eps_z) && (eps_z <= 1.0), "eps_z in [0,1]");
+
+   kappa = freq * M_PI;
+
+   ostringstream file_name;
+   if (use_monitor)
+   {
+      file_name << "ABS-"
+                << "O" << order
+                << "I" << (int) integrator_type
+                << "S" << (int) solver_type
+                << "A" << assembly_type_int
+                << ".csv";
+   }
+
+   switch (assembly_type_int)
+   {
+      case 0:
+         assembly_type = AssemblyLevel::LEGACY;
+         break;
+      case 1:
+         assembly_type = AssemblyLevel::LEGACYFULL;
+         break;
+      case 2:
+         assembly_type = AssemblyLevel::FULL;
+         break;
+      case 3:
+         assembly_type = AssemblyLevel::ELEMENT;
+         break;
+      case 4:
+         assembly_type = AssemblyLevel::PARTIAL;
+         break;
+      case 5:
+         assembly_type = AssemblyLevel::NONE;
+         break;
+      default:
+         MFEM_ABORT("Unsupported option!");
+   }
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   /// 3. Read the serial mesh from the given mesh file.
+   ///    For convinience, the meshes are available in
+   ///    ./meshes, and the number of serial and parallel
+   ///    refinements are user-defined.
+   Mesh *serial_mesh = new Mesh(mesh_file);
+   for (int ls = 0; ls < refine_serial; ls++)
+   {
+      serial_mesh->UniformRefinement();
+   }
+
+   /// 4. Define a parallel mesh by a partitioning of the serial mesh.
+   ///    Number of parallel refinements given by the user.
+   ParMesh *mesh = new ParMesh(MPI_COMM_WORLD, *serial_mesh);
+   delete serial_mesh;
+   for (int lp = 0; lp < refine_parallel; lp++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   dim = mesh->Dimension();
+   space_dim = mesh->SpaceDimension();
+
+   bool cond_z = (dim < 3)?true:(eps_z != 0); // lazy check
+   if (eps_y != 0.0 && cond_z)
+   {
+      if (dim < 3) { eps_z = 0.0; }
+      common::KershawTransformation kershawT(dim, eps_y, eps_z);
+      mesh->Transform(kershawT);
+   }
+
+   /// 5. Define a finite element space on the mesh. We use different spaces
+   ///    and collections for different systems.
+   ///    - H1-conforming Lagrange elements for the H1-mass matrix and the
+   ///      diffusion problem.
+   /// TODO(Gabriel): Elasticity not implemented yet for partial assembly
+   ///    - Vector H1-conforming Lagrange elements for the elasticity problem.
+   ///    - H(curl)-conforming Nedelec elements for the definite Maxwell problem.
+   FiniteElementCollection *fec;
+   ParFiniteElementSpace *fespace;
+
+   switch (integrator_type)
+   {
+      case mass: case diffusion:
+         fec = new H1_FECollection(order, dim);
+         fespace = new ParFiniteElementSpace(mesh, fec);
+         break;
+      case elasticity:
+         fec = new H1_FECollection(order, dim);
+         fespace = new ParFiniteElementSpace(mesh, fec, dim);
+         break;
+      case maxwell:
+         fec = new ND_FECollection(order, dim);
+         fespace = new ParFiniteElementSpace(mesh, fec);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check FiniteElementCollection");
+   }
+
+   HYPRE_BigInt sys_size = fespace->GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      mfem::out << "Number of unknowns: " << sys_size << endl;
+   }
+
+   /// 6. Extract the list of the essential boundary DoFs. We mark all boundary
+   ///    attibutes as essential. Then we get the list of essential DoFs.
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   /// 7. Define the linear system. Set up the bilinear form a(.,.) and the
+   ///    linear form b(.). The current implemented systems are the following:
+   ///    - (u,v), i.e., L2-projection.
+   ///    - (grad(u), grad(v)), i.e., diffusion operator.
+   ///    - (div(u), div(v)) + (e(u),e(v)), i.e., elasticity operator.
+   ///    - (curl(u), curl(v)) + (u,v), i.e., definite Maxwell operator.
+   ///    The linear form has the standard form (f,v).
+   ///    Define the matrices and vectors associated to the forms, and project
+   ///    the required boundary data into the GridFunction solution.
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   ParLinearForm *b = new ParLinearForm(fespace);
+
+   // These pointers are owned by the forms
+   LinearFormIntegrator *lfi = nullptr;
+   BilinearFormIntegrator *bfi = nullptr;
+   // Required for a static_cast
+   SumIntegrator *sum_bfi = nullptr;
+
+   // These pointers are not owned by the integrators
+   FunctionCoefficient *scalar_u = nullptr;
+   FunctionCoefficient *scalar_f = nullptr;
+   VectorFunctionCoefficient *vector_u = nullptr;
+   VectorFunctionCoefficient *vector_f = nullptr;
+
+   ConstantCoefficient one(1.0);
+
+   // These variables will define the linear system
+   ParGridFunction x(fespace), y(fespace);
+   OperatorPtr A;
+   Vector B, X;
+
+   x = 0.0;
+
+   switch (integrator_type)
+   {
+      case mass:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         lfi = new DomainLFIntegrator(*scalar_u);
+         bfi = new MassIntegrator(one);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case diffusion:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         scalar_f = new FunctionCoefficient(diffusion_source);
+         lfi = new DomainLFIntegrator(*scalar_f);
+         bfi = new DiffusionIntegrator(one);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case elasticity:
+         vector_u = new VectorFunctionCoefficient(space_dim, elasticity_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, elasticity_source);
+         lfi = new VectorDomainLFIntegrator(*vector_f);
+         bfi = new ElasticityIntegrator(one, one);
+         x.ProjectBdrCoefficient(*vector_u, ess_bdr);
+         break;
+      case maxwell:
+         vector_u = new VectorFunctionCoefficient(space_dim, maxwell_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, maxwell_source);
+         lfi = new VectorFEDomainLFIntegrator(*vector_f);
+         bfi = new SumIntegrator();
+         sum_bfi = static_cast<SumIntegrator*>(bfi);
+         sum_bfi->AddIntegrator(new CurlCurlIntegrator(one));
+         sum_bfi->AddIntegrator(new VectorFEMassIntegrator(one));
+         x.ProjectBdrCoefficientTangent(*vector_u, ess_bdr);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check ParLinearForm");
+   }
+
+   a->SetAssemblyLevel(assembly_type);
+   a->AddDomainIntegrator(bfi);
+   a->Assemble();
+
+   b->AddDomainIntegrator(lfi);
+   b->Assemble();
+
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   /// 8. Construct the preconditioner. Uses AbsMult to construct an appoximation
+   ///    of the diagonal of the matrix.
+
+   Solver *jacobi = nullptr;
+   Vector ones(fespace->GetTrueVSize());
+   Vector diag(fespace->GetTrueVSize());
+
+   switch (pc_type)
+   {
+      case none:
+         break;
+      case global:
+         if (Mpi::Root())
+         {
+            mfem::out << "WARNING: p_order and q_order are "
+                      << "ignored for Abs-Value L1-Jacobi!"
+                      << endl;
+         }
+         ones = 1.0;
+         A->AbsMult(ones, diag);
+         jacobi = new OperatorJacobiSmoother(diag, ess_tdof_list);
+         break;
+      case element:
+         AssembleElementLpqJacobiDiag(*a, p_order, q_order, diag);
+         jacobi = new OperatorJacobiSmoother(diag, ess_tdof_list);
+         break;
+      default:
+         mfem_error("Invalid preconditioner type!");
+   }
+
+   /// 9. Construct the solver. The implemented solvers are the following:
+   ///    - Stationary Linear Iteration
+   ///    - Preconditioned Conjugate Gradient
+   ///    Then, solve the system with the used-selected solver.
+   Solver *solver = nullptr;
+   DataMonitor *monitor = nullptr;
+
+   switch (solver_type)
+   {
+      case sli:
+         solver = new SLISolver(MPI_COMM_WORLD);
+         break;
+      case cg:
+         solver = new CGSolver(MPI_COMM_WORLD);
+         break;
+      default:
+         mfem_error("Invalid solver type!");
+   }
+   solver->SetOperator(*A);
+
+   IterativeSolver *it_solver = dynamic_cast<IterativeSolver *>(solver);
+   if (it_solver)
+   {
+      it_solver->SetRelTol(rel_tol);
+      it_solver->SetMaxIter(max_iter);
+      it_solver->SetPrintLevel(1);
+      if (use_monitor)
+      {
+         monitor = new DataMonitor(file_name.str(), MONITOR_DIGITS);
+         it_solver->SetMonitor(*monitor);
+      }
+      if (jacobi)
+      {
+         it_solver->SetPreconditioner(*jacobi);
+      }
+   }
+
+   solver->Mult(B, X);
+
+   /// 10. Recover the solution x as a grid function. Send the data by socket
+   ///     to a GLVis server.
+   a->RecoverFEMSolution(X, *b, x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   /// 11. Compute and print the L^2 norm of the error, print elapsed times
+   {
+      real_t error = 0.0;
+      switch (integrator_type)
+      {
+         case mass: case diffusion:
+            error = x.ComputeL2Error(*scalar_u);
+            break;
+         case elasticity: case maxwell:
+            error = x.ComputeL2Error(*vector_u);
+            break;
+         default:
+            mfem_error("Invalid integrator type! Check ComputeL2Error");
+      }
+      if (Mpi::Root())
+      {
+         mfem::out << "\n|| u_h - u ||_{L^2} = " << error << "\n" << endl;
+      }
+   }
+
+   /// 12. Free the memory used
+   delete solver;
+   if (jacobi) { delete jacobi; }
+   delete a;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete mesh;
+   if (monitor) { delete monitor; }
+   if (scalar_u) { delete scalar_u; }
+   if (scalar_f) { delete scalar_f; }
+   if (vector_u) { delete vector_u; }
+   if (vector_f) { delete vector_f; }
+
+   return 0;
+}

miniapps/diag-smoothers/lpq-common.cpp

@@ -0,0 +1,408 @@
+// Copyright (c) 2010-2024, 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 "lpq-common.hpp"
+
+using namespace std;
+using namespace mfem;
+
+namespace lpq_common
+{
+
+int MONITOR_DIGITS = 20;
+int MG_MAX_ITER = 10;
+real_t MG_REL_TOL = std::sqrt(1e-10);
+
+int dim = 0;
+int space_dim = 0;
+real_t freq = 1.0;
+real_t kappa = 1.0;
+
+// Custom monitor that prints a csv-formatted file
+DataMonitor::DataMonitor(string file_name, int ndigits)
+   : os(file_name),
+     precision(ndigits)
+{
+   if (Mpi::Root())
+   {
+      mfem::out << "Saving iterations into: " << file_name << endl;
+   }
+   os << "it,res,sol" << endl;
+   os << fixed << setprecision(precision);
+}
+
+void DataMonitor::MonitorResidual(int it, real_t norm, const Vector &x,
+                                  bool final)
+{
+   os << it << "," << norm << ",";
+}
+
+void DataMonitor::MonitorSolution(int it, real_t norm, const Vector &x,
+                                  bool final)
+{
+   os << norm << endl;
+}
+
+// L(p,q) general geometric multigrid method, derived from GeometricMultigrid
+LpqGeometricMultigrid::LpqGeometricMultigrid(
+   ParFiniteElementSpaceHierarchy& fes_hierarchy,
+   Array<int>& ess_bdr,
+   IntegratorType it,
+   SolverType st,
+   real_t p_order,
+   real_t q_order)
+   : GeometricMultigrid(fes_hierarchy, ess_bdr),
+     integrator_type(it),
+     solver_type(st),
+     p_order(p_order),
+     q_order(q_order),
+     coarse_pc(nullptr),
+     one(1.0)
+{
+   ConstructCoarseOperatorAndSolver(fes_hierarchy.GetFESpaceAtLevel(0));
+   for (int l = 1; l < fes_hierarchy.GetNumLevels(); ++l)
+   {
+      ConstructOperatorAndSmoother(fes_hierarchy.GetFESpaceAtLevel(l), l);
+   }
+}
+
+void LpqGeometricMultigrid::ConstructCoarseOperatorAndSolver(
+   ParFiniteElementSpace& coarse_fespace)
+{
+   ConstructBilinearForm(coarse_fespace);
+
+   HypreParMatrix* coarse_mat = new HypreParMatrix();
+   bfs[0]->FormSystemMatrix(*essentialTrueDofs[0], *coarse_mat);
+
+   Solver* coarse_solver = nullptr;
+   switch (solver_type)
+   {
+      case sli:
+         coarse_solver = new SLISolver(MPI_COMM_WORLD);
+         break;
+      case cg:
+         coarse_solver = new CGSolver(MPI_COMM_WORLD);
+         break;
+      default:
+         mfem_error("Invalid solver type!");
+   }
+
+   coarse_pc = new OperatorLpqJacobiSmoother(*coarse_mat, *essentialTrueDofs[0],
+                                             p_order, q_order);
+
+   IterativeSolver *it_solver = dynamic_cast<IterativeSolver*>(coarse_solver);
+   if (it_solver)
+   {
+      it_solver->SetRelTol(MG_REL_TOL);
+      it_solver->SetMaxIter(MG_MAX_ITER);
+      it_solver->SetPrintLevel(-1);
+      it_solver->SetPreconditioner(*coarse_pc);
+   }
+   coarse_solver->SetOperator(*coarse_mat);
+   AddLevel(coarse_mat, coarse_solver, true, true);
+}
+
+void LpqGeometricMultigrid::ConstructOperatorAndSmoother(
+   ParFiniteElementSpace& fespace, int level)
+{
+   const Array<int> &ess_tdof_list = *essentialTrueDofs[level];
+   ConstructBilinearForm(fespace);
+
+   auto level_mat = new HypreParMatrix();
+   bfs.Last()->FormSystemMatrix(ess_tdof_list, *level_mat);
+
+   Solver* smoother = new OperatorLpqJacobiSmoother(*level_mat,
+                                                    ess_tdof_list,
+                                                    p_order,
+                                                    q_order);
+
+   AddLevel(level_mat, smoother, true, true);
+}
+
+
+void LpqGeometricMultigrid::ConstructBilinearForm(ParFiniteElementSpace&
+                                                  fespace)
+{
+   ParBilinearForm* form = new ParBilinearForm(&fespace);
+   switch (integrator_type)
+   {
+      case mass:
+         form->AddDomainIntegrator(new MassIntegrator);
+         break;
+      case diffusion:
+         form->AddDomainIntegrator(new DiffusionIntegrator);
+         break;
+      case elasticity:
+         form->AddDomainIntegrator(new ElasticityIntegrator(one, one));
+         break;
+      case maxwell:
+         form->AddDomainIntegrator(new CurlCurlIntegrator(one));
+         form->AddDomainIntegrator(new VectorFEMassIntegrator(one));
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check ParBilinearForm");
+   }
+   form->Assemble();
+   bfs.Append(form);
+}
+
+// Abs-L(1) general geometric multigrid method, derived from GeometricMultigrid
+AbsL1GeometricMultigrid::AbsL1GeometricMultigrid(
+   ParFiniteElementSpaceHierarchy& fes_hierarchy,
+   Array<int>& ess_bdr,
+   IntegratorType it,
+   SolverType st,
+   AssemblyLevel al)
+   : GeometricMultigrid(fes_hierarchy, ess_bdr),
+     integrator_type(it),
+     solver_type(st),
+     assembly_level(al),
+     coarse_pc(nullptr),
+     one(1.0)
+{
+   // BilinearForm::FormSystemMatrix does not handle the ownership of A_l.
+   // GeometricMultigrid owns the forms, and deletes them.
+   mg_owned = !(AssemblyLevel::LEGACY == assembly_level);
+
+   ConstructCoarseOperatorAndSolver(fes_hierarchy.GetFESpaceAtLevel(0));
+   for (int l = 1; l < fes_hierarchy.GetNumLevels(); ++l)
+   {
+      ConstructOperatorAndSmoother(fes_hierarchy.GetFESpaceAtLevel(l), l);
+   }
+}
+
+void AbsL1GeometricMultigrid::ConstructCoarseOperatorAndSolver(
+   ParFiniteElementSpace& coarse_fespace)
+{
+   ConstructBilinearForm(coarse_fespace);
+
+   OperatorPtr coarse_mat;
+   coarse_mat.SetType(Operator::ANY_TYPE);
+   bfs[0]->FormSystemMatrix(*essentialTrueDofs[0], coarse_mat);
+   coarse_mat.SetOperatorOwner(false);
+
+   // Create smoother
+   Vector local_ones(coarse_mat->Height());
+   Vector result(coarse_mat->Height());
+
+   local_ones = 1.0;
+   coarse_mat->AbsMult(local_ones, result);
+
+   coarse_pc = new OperatorJacobiSmoother(result, *essentialTrueDofs[0]);
+
+   Solver* coarse_solver = nullptr;
+   switch (solver_type)
+   {
+      case sli:
+         coarse_solver = new SLISolver(MPI_COMM_WORLD);
+         break;
+      case cg:
+         coarse_solver = new CGSolver(MPI_COMM_WORLD);
+         break;
+      default:
+         mfem_error("Invalid solver type!");
+   }
+   coarse_solver->SetOperator(*coarse_mat);
+
+   IterativeSolver *it_solver = dynamic_cast<IterativeSolver*>(coarse_solver);
+   if (it_solver)
+   {
+      it_solver->SetRelTol(MG_REL_TOL);
+      it_solver->SetMaxIter(MG_MAX_ITER);
+      it_solver->SetPrintLevel(-1);
+      it_solver->SetPreconditioner(*coarse_pc);
+   }
+
+   AddLevel(coarse_mat.Ptr(), coarse_solver, mg_owned, true);
+}
+
+void AbsL1GeometricMultigrid::ConstructOperatorAndSmoother(
+   ParFiniteElementSpace& fespace, int level)
+{
+   const Array<int> &ess_tdof_list = *essentialTrueDofs[level];
+   ConstructBilinearForm(fespace);
+
+   OperatorPtr level_mat;
+   level_mat.SetType(Operator::ANY_TYPE);
+   bfs.Last()->FormSystemMatrix(ess_tdof_list, level_mat);
+   level_mat.SetOperatorOwner(false);
+
+   // Create smoother
+   Vector local_ones(level_mat->Height());
+   Vector result(level_mat->Height());
+
+   local_ones = 1.0;
+   level_mat->AbsMult(local_ones, result);
+
+   Solver* smoother = new OperatorJacobiSmoother(result, ess_tdof_list);
+
+   AddLevel(level_mat.Ptr(), smoother, mg_owned, true);
+}
+
+
+void AbsL1GeometricMultigrid::ConstructBilinearForm(ParFiniteElementSpace&
+                                                    fespace)
+{
+   ParBilinearForm* form = new ParBilinearForm(&fespace);
+   form->SetAssemblyLevel(assembly_level);
+   switch (integrator_type)
+   {
+      case mass:
+         form->AddDomainIntegrator(new MassIntegrator(one));
+         break;
+      case diffusion:
+         form->AddDomainIntegrator(new DiffusionIntegrator(one));
+         break;
+      case elasticity:
+         form->AddDomainIntegrator(new ElasticityIntegrator(one, one));
+         break;
+      case maxwell:
+         form->AddDomainIntegrator(new CurlCurlIntegrator(one));
+         form->AddDomainIntegrator(new VectorFEMassIntegrator(one));
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check ParBilinearForm");
+   }
+   form->Assemble();
+   bfs.Append(form);
+}
+
+
+void AssembleElementLpqJacobiDiag(ParBilinearForm& form, real_t p, real_t q,
+                                  Vector& diag)
+{
+   ParBilinearForm temp_form(form.ParFESpace());
+   temp_form.AllocateMatrix();
+   for (int i = 0; i < form.ParFESpace()->GetNE(); ++i)
+   {
+      DenseMatrix emat_i;
+      form.ComputeElementMatrix(i, emat_i);
+      Vector right(emat_i.Height());
+      Vector temp(emat_i.Height());
+      Vector left(emat_i.Height());
+
+      right = 1.0;
+      if (q!=0.0)
+      {
+         emat_i.GetDiag(right);
+         right.PowerAbs(-q);
+      }
+
+      emat_i.PowAbsMult(p, right, temp);
+
+      left = temp;
+      if (1.0 + q - p!= 0.0)
+      {
+         emat_i.GetDiag(left);
+         left.PowerAbs(1.0 + q - p);
+         left *= temp;
+      }
+
+      DenseMatrix temp_emat_i;
+      temp_emat_i.Diag(left.GetData(), left.Size());
+      temp_form.AssembleElementMatrix(i, temp_emat_i, 1);
+   }
+   temp_form.Finalize();
+   auto mat = temp_form.ParallelAssemble();
+   mat->AssembleDiagonal(diag);
+   delete mat;
+}
+
+real_t diffusion_solution(const Vector &x)
+{
+   if (dim == 3)
+   {
+      return sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2)) + 1.0;
+   }
+   else
+   {
+      return sin(kappa * x(0)) * sin(kappa * x(1)) + 1.0;
+   }
+}
+
+real_t diffusion_source(const Vector &x)
+{
+   if (dim == 3)
+   {
+      return dim * kappa * kappa * sin(kappa * x(0)) * sin(kappa * x(1)) * sin(
+                kappa * x(2));
+   }
+   else
+   {
+      return dim * kappa * kappa * sin(kappa * x(0)) * sin(kappa * x(1));
+   }
+}
+
+void elasticity_solution(const Vector &x, Vector &u)
+{
+   if (dim == 3)
+   {
+      u(0) = sin(kappa * x(0));
+      u(1) = sin(kappa * x(1));
+      u(2) = sin(kappa * x(2));
+   }
+   else
+   {
+      u(0) = sin(kappa * x(0));
+      u(1) = sin(kappa * x(1));
+      if (x.Size() == 3) { u(2) = 0.0; }
+   }
+}
+
+void elasticity_source(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = 3.0 * kappa * kappa * sin(kappa * x(0));
+      f(1) = 3.0 * kappa * kappa * sin(kappa * x(1));
+      f(2) = 3.0 * kappa * kappa * sin(kappa * x(2));
+   }
+   else
+   {
+      f(0) = 3.0 * kappa * kappa * sin(kappa * x(0));
+      f(1) = 3.0 * kappa * kappa * sin(kappa * x(1));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}
+
+void maxwell_solution(const Vector &x, Vector &u)
+{
+   if (dim == 3)
+   {
+      u(0) = sin(kappa * x(1));
+      u(1) = sin(kappa * x(2));
+      u(2) = sin(kappa * x(0));
+   }
+   else
+   {
+      u(0) = sin(kappa * x(1));
+      u(1) = sin(kappa * x(0));
+      if (x.Size() == 3) { u(2) = 0.0; }
+   }
+}
+
+void maxwell_source(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
+      f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
+   }
+   else
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}
+
+} // end namespace lpq_common

miniapps/diag-smoothers/lpq-common.hpp

@@ -0,0 +1,159 @@
+// Copyright (c) 2010-2024, 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_LPQ_COMMON_HPP
+#define MFEM_LPQ_COMMON_HPP
+
+#include "mfem.hpp"
+#include "miniapps/common/mfem-common.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+namespace lpq_common
+{
+
+extern int MONITOR_DIGITS;
+extern int MG_MAX_ITER;
+extern real_t MG_REL_TOL;
+
+extern int dim;
+extern int space_dim;
+extern real_t freq;
+extern real_t kappa;
+
+// Enumerator for the different solvers to implement
+enum SolverType
+{
+   sli,
+   cg,
+   num_solvers,  // last
+};
+
+// Enumerator for the different integrators to implement
+enum IntegratorType
+{
+   mass,
+   diffusion,
+   elasticity,
+   maxwell,
+   num_integrators,  // last
+};
+
+// Enumerator for the different types of preconditioners
+enum LpqType
+{
+   none,
+   global,
+   element,
+   num_lpq_pc,  // last
+};
+
+/// @brief Custom monitor that prints a csv-formatted file
+class DataMonitor : public IterativeSolverMonitor
+{
+private:
+   ofstream os;
+   int precision;
+public:
+   DataMonitor(string file_name, int ndigits);
+   void MonitorResidual(int it, real_t norm, const Vector &x, bool final);
+   void MonitorSolution(int it, real_t norm, const Vector &x, bool final);
+};
+
+/// @brief L(p,q)-Jacobi custom general geometric multigrid method.
+///
+/// Intermediate levels use L(p,q)-Jacobi preconditioner. Coarsest level uses a
+/// used-selected solver with an L(p,q)-Jacobi preconditioner. Assumes that
+/// the forms will be fully assembled.
+class LpqGeometricMultigrid : public GeometricMultigrid
+{
+public:
+   LpqGeometricMultigrid(ParFiniteElementSpaceHierarchy& fes_hierarchy,
+                         Array<int>& ess_bdr,
+                         IntegratorType it,
+                         SolverType st,
+                         real_t p_order,
+                         real_t q_order);
+
+   ~LpqGeometricMultigrid() { delete coarse_pc; }
+
+private:
+   IntegratorType integrator_type;
+   SolverType solver_type;
+   real_t p_order;
+   real_t q_order;
+   OperatorLpqJacobiSmoother* coarse_pc;
+   ConstantCoefficient one;
+
+   void ConstructCoarseOperatorAndSolver(ParFiniteElementSpace& coarse_fespace);
+
+   void ConstructOperatorAndSmoother(ParFiniteElementSpace& fespace, int level);
+
+   void ConstructBilinearForm(ParFiniteElementSpace& fespace);
+
+};
+
+/// @brief Abs-L(1)-Jacobi custom general geometric multigrid method.
+///
+/// Intermediate levels use Abs-L(1)-Jacobi preconditioner by applying the
+/// level matrix to the constant vector one. These are wrapped by an
+/// OperatorJacobiSmoother. Coarsest level uses a used-selected solver
+/// with an Abs-L(1)-Jacobi smoother. The assembly level is user-defined.
+///
+/// @warning The construction of the smoother is based on the application of
+/// AbsMult, which usually unfolds component-wise. E.g., if `A = B C`, then
+/// `|A|x = |B|(|C| x)`.
+class AbsL1GeometricMultigrid : public GeometricMultigrid
+{
+public:
+   AbsL1GeometricMultigrid(ParFiniteElementSpaceHierarchy& fes_hierarchy,
+                           Array<int>& ess_bdr,
+                           IntegratorType it,
+                           SolverType st,
+                           AssemblyLevel al);
+
+   ~AbsL1GeometricMultigrid() { delete coarse_pc; }
+
+   bool GetOwnershipLevelOperators() const { return mg_owned; }
+
+private:
+   IntegratorType integrator_type;
+   SolverType solver_type;
+   AssemblyLevel assembly_level;
+   bool mg_owned;
+   OperatorJacobiSmoother* coarse_pc;
+   ConstantCoefficient one;
+
+   void ConstructCoarseOperatorAndSolver(ParFiniteElementSpace& coarse_fespace);
+
+   void ConstructOperatorAndSmoother(ParFiniteElementSpace& fespace, int level);
+
+   void ConstructBilinearForm(ParFiniteElementSpace& fespace);
+
+};
+
+void AssembleElementLpqJacobiDiag(ParBilinearForm& form, real_t p, real_t q,
+                                  Vector& diag);
+
+real_t diffusion_solution(const Vector &x);
+real_t diffusion_source(const Vector &x);
+
+void elasticity_solution(const Vector &x, Vector &u);
+void elasticity_source(const Vector &x, Vector &f);
+
+void maxwell_solution(const Vector &x, Vector &u);
+void maxwell_source(const Vector &x, Vector &f);
+
+} // namespace lpq_jacobi
+#endif // MFEM_LPQ_COMMON_HPP

miniapps/diag-smoothers/lpq-jacobi.cpp

@@ -0,0 +1,412 @@
+// Copyright (c) 2010-2024, 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.
+//
+//                     ---------------------------------
+//                      L(p,q)-Jacobi smoothers miniapp
+//                     ---------------------------------
+//
+// This miniapp illustrates the use of a family of smoothers and preconditioners
+// of the L(p,q)-Jacobi family. These preconditioners are tested in different
+// settings. We use Stationary Linear Iterations and Preconditioned Conjugate
+// Gradient as the main solvers. We consider a H1-mass matrix, a diffusion matrix,
+// a elasticity system, and a definite Maxwell system.
+//
+// The preconditioner can be defined at run-time. Similarly, the mesh can be
+// modified by a Kershaw transformation at run-time. Relative tolerance and
+// maximum number of iterations can be modified as well.
+//
+// There is an analogous driver with a multigrid method implemented (cf. ex26(p)).
+//
+// Compile with: make lpq-jacobi
+//
+// Sample runs: mpirun -np 4 ./lpq-jacobi
+//              mpirun -np 4 ./lpq-jacobi -p 2.0 -q 1.0
+//              mpirun -np 4 ./lpq-jacobi -s 1 -i 3
+//              mpirun -np 4 ./lpq-jacobi -m meshes/icf.mesh -f 0.5
+//              mpirun -np 4 ./lpq-jacobi -rs 2 -rp 0
+//              mpirun -np 4 ./lpq-jacobi -t 1e5 -ni 100 -vis
+//              mpirun -np 4 ./lpq-jacobi -m meshes/beam-tet.mesh -Ky 0.5 -Kz 0.5
+
+#include "lpq-common.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace lpq_common;
+
+
+int main(int argc, char *argv[])
+{
+   /// 1. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   /// 2. Parse command line options.
+   string mesh_file = "meshes/cube.mesh";
+   // System properties
+   int order = 1;
+   SolverType solver_type = cg;
+   IntegratorType integrator_type = diffusion;
+   LpqType pc_type = global;
+   // Number of refinements
+   int refine_serial = 0;
+   int refine_parallel = 0;
+   // Preconditioner parameters
+   real_t p_order = 1.0;
+   real_t q_order = 0.0;
+   // Solver parameters
+   real_t rel_tol = 1e-10;
+   real_t max_iter = 3000;
+   // Kershaw Transformation
+   real_t eps_y = 0.0;
+   real_t eps_z = 0.0;
+   // Other options
+   string device_config = "cpu";
+   bool use_monitor = false;
+   bool visualization = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree)");
+   args.AddOption((int*)&solver_type, "-s", "--solver",
+                  "Solvers to be considered:"
+                  "\n\t0: Stationary Linear Iteration"
+                  "\n\t1: Preconditioned Conjugate Gradient");
+   args.AddOption((int*)&integrator_type, "-i", "--integrator",
+                  "Integrators to be considered:"
+                  "\n\t0: MassIntegrator"
+                  "\n\t1: DiffusionIntegrator"
+                  "\n\t2: ElasticityIntegrator"
+                  "\n\t3: CurlCurlIntegrator + VectorFEMassIntegrator");
+   args.AddOption((int*)&pc_type, "-pc", "--preconditioner",
+                  "Preconditioners to be considered:"
+                  "\n\t0: No preconditioner"
+                  "\n\t1: L(p,q)-Jacobi preconditioner"
+                  "\n\t2: Element L(p,q)-Jacobi preconditioner");
+   args.AddOption(&refine_serial, "-rs", "--refine-serial",
+                  "Number of serial refinements");
+   args.AddOption(&refine_parallel, "-rp", "--refine-parallel",
+                  "Number of parallel refinements");
+   args.AddOption(&p_order, "-p", "--p-order",
+                  "P-order for L(p,q)-Jacobi preconditioner");
+   args.AddOption(&q_order, "-q", "--q-order",
+                  "Q-order for L(p,q)-Jacobi preconditioner");
+   args.AddOption(&rel_tol, "-t", "--tolerance",
+                  "Relative tolerance for the iterative solver");
+   args.AddOption(&max_iter, "-ni", "--iterations",
+                  "Maximum number of iterations");
+   args.AddOption(&eps_y, "-Ky", "--Kershaw-y",
+                  "Kershaw transform factor, eps_y in (0,1]");
+   args.AddOption(&eps_z, "-Kz", "--Kershaw-z",
+                  "Kershaw transform factor, eps_z in (0,1]");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&use_monitor, "-mon", "--monitor", "-no-mon",
+                  "--no-monitor",
+                  "Enable or disable Data Monitor.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   MFEM_ASSERT(p_order > 0.0, "p needs to be positive");
+   MFEM_ASSERT((0 <= solver_type) && (solver_type < num_solvers), "");
+   MFEM_ASSERT((0 <= integrator_type) && (integrator_type < num_integrators), "");
+   MFEM_ASSERT((0 <= pc_type) && (pc_type < num_lpq_pc), "");
+   MFEM_ASSERT((0.0 <= eps_y) && (eps_y <= 1.0), "eps_y in [0,1]");
+   MFEM_ASSERT((0.0 <= eps_z) && (eps_z <= 1.0), "eps_z in [0,1]");
+
+   kappa = freq * M_PI;
+
+   ostringstream file_name;
+   if (use_monitor)
+   {
+      file_name << "LPQ-"
+                << "O" << order
+                << "I" << (int) integrator_type
+                << "S" << (int) solver_type
+                << fixed << setprecision(4)
+                << "P" << (int) (p_order * 1000)
+                << "Q" << (int) (q_order * 1000)
+                << ".csv";
+   }
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   /// 3. Read the serial mesh from the given mesh file.
+   ///    For convinience, the meshes are available in
+   ///    ./meshes, and the number of serial and parallel
+   ///    refinements are user-defined.
+   Mesh *serial_mesh = new Mesh(mesh_file);
+   for (int ls = 0; ls < refine_serial; ls++)
+   {
+      serial_mesh->UniformRefinement();
+   }
+
+   /// 4. Define a parallel mesh by a partitioning of the serial mesh.
+   ///    Number of parallel refinements given by the user. If defined,
+   ///    apply Kershaw transformation.
+   ParMesh *mesh = new ParMesh(MPI_COMM_WORLD, *serial_mesh);
+   delete serial_mesh;
+   for (int lp = 0; lp < refine_parallel; lp++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   dim = mesh->Dimension();
+   space_dim = mesh->SpaceDimension();
+
+   bool cond_z = (dim < 3)?true:(eps_z != 0.0); // lazy check
+   if (eps_y != 0.0 && cond_z)
+   {
+      if (dim < 3) { eps_z = 0.0; }
+      common::KershawTransformation kershawT(dim, eps_y, eps_z);
+      mesh->Transform(kershawT);
+   }
+
+   /// 5. Define a finite element space on the mesh. We use different spaces
+   ///    and collections for different systems.
+   ///    - H1-conforming Lagrange elements for the H1-mass matrix and the
+   ///      diffusion problem.
+   ///    - Vector H1-conforming Lagrange elements for the elasticity problem.
+   ///    - H(curl)-conforming Nedelec elements for the definite Maxwell problem.
+   FiniteElementCollection *fec;
+   ParFiniteElementSpace *fespace;
+   switch (integrator_type)
+   {
+      case mass: case diffusion:
+         fec = new H1_FECollection(order, dim);
+         fespace = new ParFiniteElementSpace(mesh, fec);
+         break;
+      case elasticity:
+         fec = new H1_FECollection(order, dim);
+         fespace = new ParFiniteElementSpace(mesh, fec, dim);
+         break;
+      case maxwell:
+         fec = new ND_FECollection(order, dim);
+         fespace = new ParFiniteElementSpace(mesh, fec);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check FiniteElementCollection");
+   }
+
+   HYPRE_BigInt sys_size = fespace->GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      mfem::out << "Number of unknowns: " << sys_size << endl;
+   }
+
+   /// 6. Extract the list of the essential boundary DoFs. We mark all boundary
+   ///    attibutes as essential. Then we get the list of essential DoFs.
+   Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   /// 7. Define the linear system. Set up the bilinear form a(.,.) and the
+   ///    linear form b(.). The current implemented systems are the following:
+   ///    - (u,v), i.e., L2-projection.
+   ///    - (grad(u), grad(v)), i.e., diffusion operator.
+   ///    - (div(u), div(v)) + (e(u),e(v)), i.e., elasticity operator.
+   ///    - (curl(u), curl(v)) + (u,v), i.e., definite Maxwell operator.
+   ///    The linear form has the standard form (f,v).
+   ///    Define the matrices and vectors associated to the forms, and project
+   ///    the required boundary data into the GridFunction solution.
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   ParLinearForm *b = new ParLinearForm(fespace);
+
+   // These pointers are owned by the forms
+   LinearFormIntegrator *lfi = nullptr;
+   BilinearFormIntegrator *bfi = nullptr;
+   // Required for a static_cast
+   SumIntegrator *sum_bfi = nullptr;
+
+   // These pointers are not owned by the integrators
+   FunctionCoefficient *scalar_u = nullptr;
+   FunctionCoefficient *scalar_f = nullptr;
+   VectorFunctionCoefficient *vector_u = nullptr;
+   VectorFunctionCoefficient *vector_f = nullptr;
+
+   ConstantCoefficient one(1.0);
+
+   // These variables will define the linear system
+   ParGridFunction x(fespace);
+   HypreParMatrix A;
+   Vector B, X;
+
+   x = 0.0;
+
+   switch (integrator_type)
+   {
+      case mass:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         lfi = new DomainLFIntegrator(*scalar_u);
+         bfi = new MassIntegrator(one);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case diffusion:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         scalar_f = new FunctionCoefficient(diffusion_source);
+         lfi = new DomainLFIntegrator(*scalar_f);
+         bfi = new DiffusionIntegrator(one);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case elasticity:
+         vector_u = new VectorFunctionCoefficient(space_dim, elasticity_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, elasticity_source);
+         lfi = new VectorDomainLFIntegrator(*vector_f);
+         bfi = new ElasticityIntegrator(one, one);
+         x.ProjectBdrCoefficient(*vector_u, ess_bdr);
+         break;
+      case maxwell:
+         vector_u = new VectorFunctionCoefficient(space_dim, maxwell_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, maxwell_source);
+         lfi = new VectorFEDomainLFIntegrator(*vector_f);
+         bfi = new SumIntegrator();
+         sum_bfi = static_cast<SumIntegrator*>(bfi);
+         sum_bfi->AddIntegrator(new CurlCurlIntegrator(one));
+         sum_bfi->AddIntegrator(new VectorFEMassIntegrator(one));
+         x.ProjectBdrCoefficientTangent(*vector_u, ess_bdr);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check ParLinearForm");
+   }
+
+   a->AddDomainIntegrator(bfi);
+   a->Assemble();
+   b->AddDomainIntegrator(lfi);
+   b->Assemble();
+
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   /// 8. Construct the preconditioner. User-inputs define the p_order and q_order
+   ///    of the L(p,q)-Jacobi type smoother.
+   // D_{p,q} = diag( D^{1+q-p} |A|^p D^{-q} 1) , where D = diag(A)
+
+   Solver *lpq_jacobi = nullptr;
+   real_t bound = 0.0;
+   Vector diag(fespace->GlobalTrueVSize());
+
+   switch (pc_type)
+   {
+      case none:
+         break;
+      case global:
+         lpq_jacobi = new OperatorLpqJacobiSmoother(A, ess_tdof_list, p_order,
+                                                    q_order);
+         bound = static_cast<OperatorLpqJacobiSmoother*>
+                 (lpq_jacobi)->CheckSpectralBoundConstant();
+         break;
+      case element:
+         AssembleElementLpqJacobiDiag(*a, p_order, q_order, diag);
+         lpq_jacobi = new OperatorJacobiSmoother(diag, ess_tdof_list);
+         break;
+      default:
+         mfem_error("Invalid preconditioner type!");
+   }
+
+   /// 9. Construct the solver. The implemented solvers are the following:
+   ///    - Stationary Linear Iteration
+   ///    - Preconditioned Conjugate Gradient
+   ///    Then, solve the system with the used-selected solver.
+   Solver *solver = nullptr;
+   DataMonitor *monitor = nullptr;
+
+   switch (solver_type)
+   {
+      case sli:
+         solver = new SLISolver(MPI_COMM_WORLD);
+         break;
+      case cg:
+         solver = new CGSolver(MPI_COMM_WORLD);
+         break;
+      default:
+         mfem_error("Invalid solver type!");
+   }
+   solver->SetOperator(A);
+
+   IterativeSolver *it_solver = dynamic_cast<IterativeSolver *>(solver);
+   if (it_solver)
+   {
+      it_solver->SetRelTol(rel_tol);
+      it_solver->SetMaxIter(max_iter);
+      it_solver->SetPrintLevel(1);
+      if (use_monitor)
+      {
+         monitor = new DataMonitor(file_name.str(), MONITOR_DIGITS);
+         it_solver->SetMonitor(*monitor);
+      }
+      if (lpq_jacobi)
+      {
+         it_solver->SetPreconditioner(*lpq_jacobi);
+      }
+   }
+
+   solver->Mult(B, X);
+
+   /// 10. Recover the solution x as a grid function. Send the data by socket
+   ///     to a GLVis server.
+   a->RecoverFEMSolution(X, *b, x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   /// 11. Compute and print the L^2 norm of the error
+   {
+      real_t error = 0.0;
+      switch (integrator_type)
+      {
+         case mass: case diffusion:
+            error = x.ComputeL2Error(*scalar_u);
+            break;
+         case elasticity: case maxwell:
+            error = x.ComputeL2Error(*vector_u);
+            break;
+         default:
+            mfem_error("Invalid integrator type! Check ComputeL2Error");
+      }
+      if (Mpi::Root())
+      {
+         mfem::out << "\n|| u_h - u ||_{L^2} = " << error << "\n" << endl;
+         if (pc_type==global) { mfem::out << "Spectral bound is: " << bound << endl; }
+      }
+   }
+
+   /// 12. Free the memory used
+   delete solver;
+   if (lpq_jacobi) { delete lpq_jacobi; }
+   delete a;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete mesh;
+   if (monitor) { delete monitor; }
+   if (scalar_u) { delete scalar_u; }
+   if (scalar_f) { delete scalar_f; }
+   if (vector_u) { delete vector_u; }
+   if (vector_f) { delete vector_f; }
+
+   return 0;
+}

miniapps/diag-smoothers/makefile

@@ -0,0 +1,140 @@
+# Copyright (c) 2010-2024, 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/lpq-jacobi/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+
+-include $(DEFAULTS_MK)
+-include $(CONFIG_MK)
+
+LPQ_COMMON_SRC = lpq-common.cpp
+LPQ_COMMON_OBJ = $(LPQ_COMMON_SRC:.cpp=.o)
+
+PAR_MINIAPPS = lpq-jacobi mg-lpq-jacobi abs-l1-jacobi mg-abs-l1-jacobi
+
+MINIAPPS = $(if $(MFEM_USE_MPI:NO=),$(PAR_MINIAPPS),)
+
+COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
+COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+# Mesh directories and files
+# Copy meshes into ./meshes
+SRC_MESHING_DIR = ../meshing
+SRC_GSLIB_DIR = ../gslib
+SRC_DATA_DIR = ../../data
+MESHES_DEST_DIR = ./meshes
+
+MESH_FILES_FROM_MESHING = icf.mesh cube.mesh
+MESH_FILES_FROM_GSLIB = triple-pt-1.mesh triple-pt-2.mesh
+MESH_FILES_FROM_DATA = beam-tet.mesh square-disc-p2.mesh fichera-mixed-p2.mesh \
+					   amr-quad.mesh ref-cube.mesh ref-square.mesh
+
+SRC_MESH_FILES = $(addprefix $(SRC_MESHING_DIR)/, $(MESH_FILES_FROM_MESHING)) \
+                 $(addprefix $(SRC_GSLIB_DIR)/, $(MESH_FILES_FROM_GSLIB)) \
+                 $(addprefix $(SRC_DATA_DIR)/, $(MESH_FILES_FROM_DATA))
+
+# Phony targets
+.PHONY: all lib-common copy_meshes clean clean-build clean-exec test
+
+# Main targets
+all: copy_meshes $(MINIAPPS)
+
+# Build rules
+$(MINIAPPS): %: %.cpp $(LPQ_COMMON_OBJ) lpq-common.hpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_FLAGS) $< $(LPQ_COMMON_OBJ) -o $@ $(COMMON_LIB) $(MFEM_LIBS)
+
+$(LPQ_COMMON_OBJ): %.o: %.cpp
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+# Library and mesh targets
+lib-common:
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
+
+copy_meshes: $(MESHES_DEST_DIR) $(addprefix $(MESHES_DEST_DIR)/, $(notdir $(SRC_MESH_FILES)))
+	@echo "Copying meshing miniapps data files ..."
+
+# Create the destination directory if it doesn't exist
+$(MESHES_DEST_DIR):
+	mkdir -p $(MESHES_DEST_DIR)
+
+$(MESHES_DEST_DIR)/%: $(SRC_MESHING_DIR)/% | $(MESHES_DEST_DIR)
+	cp $< $@
+
+$(MESHES_DEST_DIR)/%: $(SRC_GSLIB_DIR)/% | $(MESHES_DEST_DIR)
+	cp $< $@
+
+$(MESHES_DEST_DIR)/%: $(SRC_DATA_DIR)/% | $(MESHES_DEST_DIR)
+	cp $< $@
+
+# Test targets
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+
+lpq-jacobi-test-par: lpq-jacobi
+	@$(call mfem-test,$<, $(RUN_MPI), Lpq-Jqcobi miniapp, \
+		-m meshes/cube.mesh \
+		-rs 1 -rp 1 \
+		-s 1 -i 1 \
+		-p 1.5 -q 0.75 \
+		-pc 1 \
+		-no-mon)
+
+abs-l1-jacobi-test-par: abs-l1-jacobi
+	@$(call mfem-test,$<, $(RUN_MPI), Abs-value-L1-Jqcobi miniapp, \
+		-m meshes/cube.mesh \
+		-rs 1 -rp 1 \
+		-s 1 -i 1 -a 4 \
+		-pc 1 \
+		-no-mon)
+
+mg-lpq-jacobi-test-par: mg-lpq-jacobi
+	@$(call mfem-test,$<, $(RUN_MPI), MG Lpq-Jqcobi miniapp, \
+		-m meshes/cube.mesh \
+		-rs 1 -rp 1 \
+		-ol 1 -gl 1 \
+		-s 1 -i 1 \
+		-p 1.5 -q 0.75 \
+		-no-mon)
+
+mg-abs-l1-jacobi-test-par: mg-abs-l1-jacobi
+	@$(call mfem-test,$<, $(RUN_MPI), MG Abs-value-L1-Jqcobi miniapp, \
+		-m meshes/cube.mesh \
+		-rs 1 -rp 1\
+		-ol 1 -gl 1 \
+		-s 1 -i 1 -a 4 \
+		-no-mon)
+
+# Clean targets
+clean: clean-build clean-exec
+	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common clean
+
+clean-build:
+	rm -f *.o *~ $(MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	rm -rf $(MESHES_DEST_DIR)
+	rm -f *.csv
+
+# Error handling
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)

miniapps/diag-smoothers/mg-abs-l1-jacobi.cpp

@@ -0,0 +1,444 @@
+// Copyright (c) 2010-2024, 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.
+//
+//                   --------------------------------------
+//                    MG Abs L(1)-Jacobi smoothers miniapp
+//                   --------------------------------------
+//
+// (See lpq-jacobi.cpp first)
+//
+// This miniapp illustrates the use of an absolute value L(1)-Jacobi smoother.
+// We use a multigrid approach (cf. ex26(p)). The global solver and
+// the coarse level solver are user-selected. The current options are SLI and
+// PCG. The intermediate levels are directy smoothed with the absolute value
+// L(1)-Jacobi preconditioner. The systems to solve correspond to a mass matrix,
+// a difussion system, an elasticity system, and a definite Maxwell system.
+//
+// The preconditioner can be defined at run-time. Similarly, the mesh can be
+// modified by a Kershaw transformation at run-time. Relative tolerance and
+// maximum number of iterations can be modified as well.
+//
+// Compile with: make mg-lpq-jacobi
+//
+// Sample runs: mpirun -np 4 ./mg-abs-l1-jacobi
+//              mpirun -np 4 ./mg-abs-l1-jacobi -s 1 -i 3
+//              mpirun -np 4 ./mg-abs-l1-jacobi -m meshes/icf.mesh -f 0.5
+//              mpirun -np 4 ./mg-abs-l1-jacobi -rs 2 -rp 0
+//              mpirun -np 4 ./mg-abs-l1-jacobi -t 1e5 -ni 100 -vis
+//              mpirun -np 4 ./mg-abs-l1-jacobi -m meshes/beam-tet.mesh -Ky 0.5 -Kz 0.5
+
+#include "lpq-common.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace lpq_common;
+
+int main(int argc, char *argv[])
+{
+   /// 1. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   /// 2. Parse command line options.
+   string mesh_file = "meshes/cube.mesh";
+   // System properties
+   int order = 1;
+   SolverType solver_type = cg;
+   IntegratorType integrator_type = diffusion;
+   int assembly_type_int = 4;
+   AssemblyLevel assembly_type;
+   // Number of refinements
+   int refine_serial = 0;
+   int refine_parallel = 0;
+   // Number of geometric and order levels
+   int geometric_levels = 1;
+   int order_levels = 1;
+   // Solver parameters
+   real_t rel_tol = 1e-10;
+   real_t max_iter = 3000;
+   // Kershaw Transformation
+   real_t eps_y = 0.0;
+   real_t eps_z = 0.0;
+   // Other options
+   string device_config = "cpu";
+   bool use_monitor = false;
+   bool visualization = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree)");
+   args.AddOption(&geometric_levels, "-gl", "--geometric-levels",
+                  "Number of geometric refinements (levels) done prior to order refinements.");
+   args.AddOption(&order_levels, "-ol", "--order-levels",
+                  "Number of order refinements (levels). "
+                  "Finest level in the hierarchy has order 2^{or}.");
+   args.AddOption((int*)&solver_type, "-s", "--solver",
+                  "Solvers to be considered:"
+                  "\n\t0: Stationary Linear Iteration"
+                  "\n\t1: Preconditioned Conjugate Gradient");
+   args.AddOption((int*)&integrator_type, "-i", "--integrator",
+                  "Integrators to be considered:"
+                  "\n\t0: MassIntegrator"
+                  "\n\t1: DiffusionIntegrator"
+                  "\n\t2: ElasticityIntegrator"
+                  "\n\t3: CurlCurlIntegrator + VectorFEMassIntegrator");
+   args.AddOption(&assembly_type_int, "-a", "--assembly",
+                  "Assembly level to be considered:"
+                  "\n\t0: LEGACY"
+                  "\n\t1: LEGACYFULL (Deprecated)"
+                  "\n\t2: FULL"
+                  "\n\t3: ELEMENT"
+                  "\n\t4: PARTIAL"
+                  "\n\t5: NONE");
+   args.AddOption(&refine_serial, "-rs", "--refine-serial",
+                  "Number of serial refinements");
+   args.AddOption(&refine_parallel, "-rp", "--refine-parallel",
+                  "Number of parallel refinements");
+   args.AddOption(&rel_tol, "-t", "--tolerance",
+                  "Relative tolerance for the iterative solver");
+   args.AddOption(&max_iter, "-ni", "--iterations",
+                  "Maximum number of iterations");
+   args.AddOption(&eps_y, "-Ky", "--Kershaw-y",
+                  "Kershaw transform factor, eps_y in (0,1]");
+   args.AddOption(&eps_z, "-Kz", "--Kershaw-z",
+                  "Kershaw transform factor, eps_z in (0,1]");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&use_monitor, "-mon", "--monitor", "-no-mon",
+                  "--no-monitor",
+                  "Enable or disable Data Monitor.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   MFEM_ASSERT((0 <= solver_type) && (solver_type < num_solvers), "");
+   MFEM_ASSERT((0 <= integrator_type) && (integrator_type < num_integrators), "");
+   MFEM_ASSERT((0 <= assembly_type_int) && (assembly_type_int < 6), "");
+   MFEM_ASSERT(geometric_levels >= 0, "geometric_level needs to be non-negative");
+   MFEM_ASSERT(order_levels >= 0, "order_level needs to be non-negative");
+   MFEM_ASSERT((0.0 <= eps_y) && (eps_y <= 1.0), "eps_y in [0,1]");
+   MFEM_ASSERT((0.0 <= eps_z) && (eps_z <= 1.0), "eps_z in [0,1]");
+
+   kappa = freq * M_PI;
+
+   ostringstream file_name;
+   if (use_monitor)
+   {
+      file_name << "MGABS-"
+                << "G" << geometric_levels
+                << "O" << order_levels
+                << "O" << order
+                << "I" << (int) integrator_type
+                << "S" << (int) solver_type
+                << "A" << assembly_type_int
+                << ".csv";
+   }
+
+   string assembly_description;
+   switch (assembly_type_int)
+   {
+      case 0:
+         assembly_type = AssemblyLevel::LEGACY;
+         assembly_description = "Using Legacy type of assembly level...";
+         break;
+      case 1:
+         assembly_type = AssemblyLevel::LEGACYFULL;
+         assembly_description =
+            "Using Legacy Full type of assembly level... (Deprecated)";
+         break;
+      case 2:
+         assembly_type = AssemblyLevel::FULL;
+         assembly_description = "Using Full type of assembly level...";
+         break;
+      case 3:
+         assembly_type = AssemblyLevel::ELEMENT;
+         assembly_description = "Using Element type of assembly level...";
+         break;
+      case 4:
+         assembly_type = AssemblyLevel::PARTIAL;
+         assembly_description = "Using Partial type of assembly level...";
+         break;
+      case 5:
+         assembly_type = AssemblyLevel::NONE;
+         assembly_description = "Using matrix-free type of assembly level...";
+         break;
+      default:
+         MFEM_ABORT("Unsupported option!");
+   }
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   /// 3. Read the serial mesh from the given mesh file.
+   ///    For convinience, the meshes are available in
+   ///    ./meshes, and the number of serial and parallel
+   ///    refinements are user-defined.
+   Mesh *serial_mesh = new Mesh(mesh_file);
+   for (int ls = 0; ls < refine_serial; ls++)
+   {
+      serial_mesh->UniformRefinement();
+   }
+
+   /// 4. Define a parallel mesh by a partitioning of the serial mesh.
+   ///    Number of parallel refinements given by the user. If defined,
+   ///    apply Kershaw transformation.
+   ParMesh *mesh = new ParMesh(MPI_COMM_WORLD, *serial_mesh);
+   delete serial_mesh;
+   for (int lp = 0; lp < refine_parallel; lp++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   dim = mesh->Dimension();
+   space_dim = mesh->SpaceDimension();
+
+   bool cond_z = (dim < 3)?true:(eps_z != 0.0); // lazy check
+   if (eps_y != 0.0 && cond_z)
+   {
+      if (dim < 3) { eps_z = 0.0; }
+      common::KershawTransformation kershawT(dim, eps_y, eps_z);
+      mesh->Transform(kershawT);
+   }
+
+   /// 5. Define a finite element space on the mesh. We use different spaces
+   ///    and collections for different systems.
+   ///    - H1-conforming Lagrange elements for the H1-mass matrix and the
+   ///      diffusion problem.
+   ///    - Vector H1-conforming Lagrange elements for the elasticity problem.
+   ///    - H(curl)-conforming Nedelec elements for the definite Maxwell problem.
+   FiniteElementCollection *fec;
+   ParFiniteElementSpace *coarse_fes;
+   switch (integrator_type)
+   {
+      case mass: case diffusion:
+         fec = new H1_FECollection(order, dim);
+         coarse_fes = new ParFiniteElementSpace(mesh, fec);
+         break;
+      case elasticity:
+         fec = new H1_FECollection(order, dim);
+         coarse_fes= new ParFiniteElementSpace(mesh, fec, dim);
+         break;
+      case maxwell:
+         fec = new ND_FECollection(order, dim);
+         coarse_fes= new ParFiniteElementSpace(mesh, fec);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check FiniteElementCollection");
+   }
+
+   if (order > 1)
+   {
+      if (Mpi::Root())
+      {
+         mfem::out << "Warning! Polynomial order provided. "
+                   << "Ignoring order level..." << endl;
+      }
+      order_levels = 0;
+   }
+
+   /// 6. Define a finite element space hierarchy for the multigrid solver.
+   ///    Define a FEC array for the order-refinement levels. Add the refinements
+   ///    to the hierarchy.
+   Array<FiniteElementCollection*> fec_array;
+   fec_array.Append(fec);
+   // Transfer ownership of mesh and coarse_fes to fes_hierarchy
+   ParFiniteElementSpaceHierarchy* fes_hierarchy = new
+   ParFiniteElementSpaceHierarchy(mesh, coarse_fes, true, true);
+
+   for (int lg = 0; lg < geometric_levels; ++lg)
+   {
+      fes_hierarchy->AddUniformlyRefinedLevel();
+   }
+   for (int lo = 0; lo < order_levels; ++lo)
+   {
+      switch (integrator_type)
+      {
+         case mass: case diffusion: case elasticity:
+            fec_array.Append(new H1_FECollection(std::pow(2, lo + 1), dim));
+            break;
+         case maxwell:
+            fec_array.Append(new ND_FECollection(std::pow(2, lo + 1), dim));
+            break;
+         default:
+            mfem_error("Invalid integrator type! "
+                       "Check FiniteElementCollection for order refinements...");
+      }
+      fes_hierarchy->AddOrderRefinedLevel(fec_array.Last());
+   }
+
+   HYPRE_BigInt sys_size = fes_hierarchy->GetFinestFESpace().GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      mfem::out << "Number of unknowns: " << sys_size << endl;
+      mfem::out << assembly_description << endl;
+   }
+
+   /// 7. Extract the list of the essential boundary DoFs. We mark all boundary
+   ///    attibutes as essential. LpqGeometricMultigrid will determine
+   ///    the DoFs per level.
+   Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   if (mesh->bdr_attributes.Size()) { ess_bdr = 1; }
+
+   /// 8. Define the linear system. Set up the linear form b(.).
+   ///    The linear form has the standar form (f,v).
+   ///    Define the matrices and vectors associated to the forms, and project
+   ///    the required boundary data into the GridFunction solution.
+   ParLinearForm *b = new ParLinearForm(&fes_hierarchy->GetFinestFESpace());
+   LinearFormIntegrator *lfi = nullptr;
+
+   // These pointers are not owned by the integrators
+   FunctionCoefficient *scalar_u = nullptr;
+   FunctionCoefficient *scalar_f = nullptr;
+   VectorFunctionCoefficient *vector_u = nullptr;
+   VectorFunctionCoefficient *vector_f = nullptr;
+
+   ConstantCoefficient one(1.0);
+
+   // These variables will define the linear system
+   ParGridFunction x(&fes_hierarchy->GetFinestFESpace());
+   OperatorPtr A;
+   Vector B, X;
+
+   x = 0.0;
+
+   switch (integrator_type)
+   {
+      case mass:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         lfi = new DomainLFIntegrator(*scalar_u);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case diffusion:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         scalar_f = new FunctionCoefficient(diffusion_source);
+         lfi = new DomainLFIntegrator(*scalar_f);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case elasticity:
+         vector_u = new VectorFunctionCoefficient(space_dim, elasticity_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, elasticity_source);
+         lfi = new VectorDomainLFIntegrator(*vector_f);
+         x.ProjectBdrCoefficient(*vector_u, ess_bdr);
+         break;
+      case maxwell:
+         vector_u = new VectorFunctionCoefficient(space_dim, maxwell_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, maxwell_source);
+         lfi = new VectorFEDomainLFIntegrator(*vector_f);
+         x.ProjectBdrCoefficientTangent(*vector_u, ess_bdr);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check ParLinearForm");
+   }
+   b->AddDomainIntegrator(lfi);
+   b->Assemble();
+
+   /// 9. Define a geometric multigrid solver. The bilinear form
+   ///    a(.,.) is assembled internally. Set up the type of cycles
+   ///    and form the linear system.
+   auto mg = new AbsL1GeometricMultigrid(*fes_hierarchy,
+                                         ess_bdr,
+                                         integrator_type,
+                                         solver_type,
+                                         assembly_type);
+   mg->SetCycleType(Multigrid::CycleType::VCYCLE, 1, 1);
+   mg->FormFineLinearSystem(x, *b, A, X, B);
+
+   A.SetOperatorOwner(mg->GetOwnershipLevelOperators());
+
+   Solver *solver = nullptr;
+   DataMonitor *monitor = nullptr;
+
+   switch (solver_type)
+   {
+      case sli:
+         solver = new SLISolver(MPI_COMM_WORLD);
+         break;
+      case cg:
+         solver = new CGSolver(MPI_COMM_WORLD);
+         break;
+      default:
+         mfem_error("Invalid solver type!");
+   }
+   solver->SetOperator(*A.Ptr());
+
+   IterativeSolver *it_solver = dynamic_cast<IterativeSolver*>(solver);
+   if (it_solver)
+   {
+      it_solver->SetRelTol(rel_tol);
+      it_solver->SetMaxIter(max_iter);
+      it_solver->SetPrintLevel(1);
+      it_solver->SetPreconditioner(*mg);
+      if (use_monitor)
+      {
+         monitor = new DataMonitor(file_name.str(), MONITOR_DIGITS);
+         it_solver->SetMonitor(*monitor);
+      }
+   }
+
+   solver->Mult(B, X);
+
+   /// 10. Recover the solution x as a grid function. Send the data by socket
+   ///     to a GLVis server.
+   mg->RecoverFineFEMSolution(X, *b, x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *fes_hierarchy->GetFinestFESpace().GetParMesh()
+               << x << flush;
+   }
+
+   /// 11. Compute and print the L^2 norm of the error
+   {
+      real_t error = 0.0;
+      switch (integrator_type)
+      {
+         case mass: case diffusion:
+            error = x.ComputeL2Error(*scalar_u);
+            break;
+         case elasticity: case maxwell:
+            error = x.ComputeL2Error(*vector_u);
+            break;
+         default:
+            mfem_error("Invalid integrator type! Check ComputeL2Error");
+      }
+      if (Mpi::Root())
+      {
+         mfem::out << "\n|| u_h - u ||_{L^2} = " << error << "\n" << endl;
+      }
+   }
+
+   /// 12. Free the memory used
+   delete mg;
+   delete solver;
+   delete b;
+   if (monitor) { delete monitor; }
+   if (scalar_u) { delete scalar_u; }
+   if (scalar_f) { delete scalar_f; }
+   if (vector_u) { delete vector_u; }
+   if (vector_f) { delete vector_f; }
+   for (int level = 0; level < fec_array.Size(); ++level)
+   {
+      delete fec_array[level];
+   }
+   delete fes_hierarchy;
+
+   return 0;
+}

miniapps/diag-smoothers/mg-lpq-jacobi.cpp

@@ -0,0 +1,411 @@
+// Copyright (c) 2010-2024, 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.
+//
+//                   --------------------------------------
+//                     MG L(p,q)-Jacobi smoothers miniapp
+//                   --------------------------------------
+//
+// (See lpq-jacobi.cpp first)
+//
+// This miniapp illustrates the use of a family of smoothers and preconditioners
+// of the L(p,q)-Jacobi family. These preconditioners are tested in different
+// settings. We use a multigrid approach (cf. ex26(p)). The global solver and
+// the coarse level solver are user-selected. The current options are SLI and
+// PCG. The intermediate levels are directy smoothed with the L(p,q)-Jacobi
+// preconditioner. The systems to solve correspond to a mass matrix, a difussion
+// system, an elasticity system, and a definite Maxwell system.
+//
+// The preconditioner can be defined at run-time. Similarly, the mesh can be
+// modified by a Kershaw transformation at run-time. Relative tolerance and
+// maximum number of iterations can be modified as well.
+//
+// Compile with: make mg-lpq-jacobi
+//
+// Sample runs: mpirun -np 4 ./mg-lpq-jacobi
+//              mpirun -np 4 ./mg-lpq-jacobi -p 2.0 -q 1.0
+//              mpirun -np 4 ./mg-lpq-jacobi -s 1 -i 3
+//              mpirun -np 4 ./mg-lpq-jacobi -m meshes/icf.mesh -f 0.5
+//              mpirun -np 4 ./mg-lpq-jacobi -rs 2 -rp 0
+//              mpirun -np 4 ./mg-lpq-jacobi -t 1e5 -ni 100 -vis
+//              mpirun -np 4 ./mg-lpq-jacobi -m meshes/beam-tet.mesh -Ky 0.5 -Kz 0.5
+
+#include "lpq-common.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace lpq_common;
+
+int main(int argc, char *argv[])
+{
+   /// 1. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   /// 2. Parse command line options.
+   string mesh_file = "meshes/cube.mesh";
+   // System properties
+   int order = 1;
+   SolverType solver_type = cg;
+   IntegratorType integrator_type = diffusion;
+   // Number of refinements
+   int refine_serial = 0;
+   int refine_parallel = 0;
+   // Number of geometric and order levels
+   int geometric_levels = 1;
+   int order_levels = 1;
+   // Preconditioner parameters
+   real_t p_order = 1.0;
+   real_t q_order = 0.0;
+   // Solver parameters
+   real_t rel_tol = 1e-10;
+   real_t max_iter = 3000;
+   // Kershaw Transformation
+   real_t eps_y = 0.0;
+   real_t eps_z = 0.0;
+   // Other options
+   string device_config = "cpu";
+   bool use_monitor = false;
+   bool visualization = false;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree)");
+   args.AddOption(&geometric_levels, "-gl", "--geometric-levels",
+                  "Number of geometric refinements (levels) done prior to order refinements.");
+   args.AddOption(&order_levels, "-ol", "--order-levels",
+                  "Number of order refinements (levels). "
+                  "Finest level in the hierarchy has order 2^{or}.");
+   args.AddOption((int*)&solver_type, "-s", "--solver",
+                  "Solvers to be considered:"
+                  "\n\t0: Stationary Linear Iteration"
+                  "\n\t1: Preconditioned Conjugate Gradient");
+   args.AddOption((int*)&integrator_type, "-i", "--integrator",
+                  "Integrators to be considered:"
+                  "\n\t0: MassIntegrator"
+                  "\n\t1: DiffusionIntegrator"
+                  "\n\t2: ElasticityIntegrator"
+                  "\n\t3: CurlCurlIntegrator + VectorFEMassIntegrator");
+   args.AddOption(&refine_serial, "-rs", "--refine-serial",
+                  "Number of serial refinements");
+   args.AddOption(&refine_parallel, "-rp", "--refine-parallel",
+                  "Number of parallel refinements");
+   args.AddOption(&p_order, "-p", "--p-order",
+                  "P-order for L(p,q)-Jacobi preconditioner");
+   args.AddOption(&q_order, "-q", "--q-order",
+                  "Q-order for L(p,q)-Jacobi preconditioner");
+   args.AddOption(&rel_tol, "-t", "--tolerance",
+                  "Relative tolerance for the iterative solver");
+   args.AddOption(&max_iter, "-ni", "--iterations",
+                  "Maximum number of iterations");
+   args.AddOption(&eps_y, "-Ky", "--Kershaw-y",
+                  "Kershaw transform factor, eps_y in (0,1]");
+   args.AddOption(&eps_z, "-Kz", "--Kershaw-z",
+                  "Kershaw transform factor, eps_z in (0,1]");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&use_monitor, "-mon", "--monitor", "-no-mon",
+                  "--no-monitor",
+                  "Enable or disable Data Monitor.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   MFEM_ASSERT(p_order > 0.0, "p needs to be positive");
+   MFEM_ASSERT((0 <= solver_type) && (solver_type < num_solvers), "");
+   MFEM_ASSERT((0 <= integrator_type) && (integrator_type < num_integrators), "");
+   MFEM_ASSERT(geometric_levels >= 0, "geometric_level needs to be non-negative");
+   MFEM_ASSERT(order_levels >= 0, "order_level needs to be non-negative");
+   MFEM_ASSERT((0.0 <= eps_y) && (eps_y <= 1.0), "eps_y in [0,1]");
+   MFEM_ASSERT((0.0 <= eps_z) && (eps_z <= 1.0), "eps_z in [0,1]");
+
+   kappa = freq * M_PI;
+
+   ostringstream file_name;
+   if (use_monitor)
+   {
+      file_name << "MGLPQ-"
+                << "G" << geometric_levels
+                << "O" << order_levels
+                << "O" << order
+                << "I" << (int) integrator_type
+                << "S" << (int) solver_type
+                << fixed << setprecision(4)
+                << "P" << (int) (p_order * 1000)
+                << "Q" << (int) (q_order * 1000)
+                << ".csv";
+   }
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   /// 3. Read the serial mesh from the given mesh file.
+   ///    For convinience, the meshes are available in
+   ///    ./meshes, and the number of serial and parallel
+   ///    refinements are user-defined.
+   Mesh *serial_mesh = new Mesh(mesh_file);
+   for (int ls = 0; ls < refine_serial; ls++)
+   {
+      serial_mesh->UniformRefinement();
+   }
+
+   /// 4. Define a parallel mesh by a partitioning of the serial mesh.
+   ///    Number of parallel refinements given by the user. If defined,
+   ///    apply Kershaw transformation.
+   ParMesh *mesh = new ParMesh(MPI_COMM_WORLD, *serial_mesh);
+   delete serial_mesh;
+   for (int lp = 0; lp < refine_parallel; lp++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   dim = mesh->Dimension();
+   space_dim = mesh->SpaceDimension();
+
+   bool cond_z = (dim < 3)?true:(eps_z != 0.0); // lazy check
+   if (eps_y != 0.0 && cond_z)
+   {
+      if (dim < 3) { eps_z = 0.0; }
+      common::KershawTransformation kershawT(dim, eps_y, eps_z);
+      mesh->Transform(kershawT);
+   }
+
+   /// 5. Define a finite element space on the mesh. We use different spaces
+   ///    and collections for different systems.
+   ///    - H1-conforming Lagrange elements for the H1-mass matrix and the
+   ///      diffusion problem.
+   ///    - Vector H1-conforming Lagrange elements for the elasticity problem.
+   ///    - H(curl)-conforming Nedelec elements for the definite Maxwell problem.
+   FiniteElementCollection *fec = nullptr;
+   ParFiniteElementSpace *coarse_fes = nullptr;
+   switch (integrator_type)
+   {
+      case mass: case diffusion:
+         fec = new H1_FECollection(order, dim);
+         coarse_fes = new ParFiniteElementSpace(mesh, fec);
+         break;
+      case elasticity:
+         fec = new H1_FECollection(order, dim);
+         coarse_fes= new ParFiniteElementSpace(mesh, fec, dim);
+         break;
+      case maxwell:
+         fec = new ND_FECollection(order, dim);
+         coarse_fes= new ParFiniteElementSpace(mesh, fec);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check FiniteElementCollection");
+   }
+
+   if (order > 1)
+   {
+      if (Mpi::Root())
+      {
+         mfem::out << "Warning! Polynomial order provided. "
+                   << "Ignoring order level..." << endl;
+      }
+      order_levels = 0;
+   }
+
+   /// 6. Define a finite element space hierarchy for the multigrid solver.
+   ///    Define a FEC array for the order-refinement levels. Add the refinements
+   ///    to the hierarchy.
+   Array<FiniteElementCollection*> fec_array;
+   fec_array.Append(fec);
+   // Transfer ownership of mesh and coarse_fes to fes_hierarchy
+   ParFiniteElementSpaceHierarchy* fes_hierarchy = new
+   ParFiniteElementSpaceHierarchy(mesh, coarse_fes, true, true);
+
+   for (int lg = 0; lg < geometric_levels; ++lg)
+   {
+      fes_hierarchy->AddUniformlyRefinedLevel();
+   }
+   for (int lo = 0; lo < order_levels; ++lo)
+   {
+      switch (integrator_type)
+      {
+         case mass: case diffusion: case elasticity:
+            fec_array.Append(new H1_FECollection(std::pow(2, lo + 1), dim));
+            break;
+         case maxwell:
+            fec_array.Append(new ND_FECollection(std::pow(2, lo + 1), dim));
+            break;
+         default:
+            mfem_error("Invalid integrator type! "
+                       "Check FiniteElementCollection for order refinements...");
+      }
+      fes_hierarchy->AddOrderRefinedLevel(fec_array.Last());
+   }
+
+   HYPRE_BigInt sys_size = fes_hierarchy->GetFinestFESpace().GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      mfem::out << "Number of unknowns: " << sys_size << endl;
+   }
+
+   /// 7. Extract the list of the essential boundary DoFs. We mark all boundary
+   ///    attibutes as essential. LpqGeometricMultigrid will determine
+   ///    the DoFs per level.
+   Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   if (mesh->bdr_attributes.Size()) { ess_bdr = 1; }
+
+   /// 8. Define the linear system. Set up the linear form b(.).
+   ///    The linear form has the standar form (f,v).
+   ///    Define the matrices and vectors associated to the forms, and project
+   ///    the required boundary data into the GridFunction solution.
+   ParLinearForm *b = new ParLinearForm(&fes_hierarchy->GetFinestFESpace());
+   LinearFormIntegrator *lfi = nullptr;
+
+   // These pointers are not owned by the integrators
+   FunctionCoefficient *scalar_u = nullptr;
+   FunctionCoefficient *scalar_f = nullptr;
+   VectorFunctionCoefficient *vector_u = nullptr;
+   VectorFunctionCoefficient *vector_f = nullptr;
+
+   ConstantCoefficient one(1.0);
+
+   // These variables will define the linear system
+   ParGridFunction x(&fes_hierarchy->GetFinestFESpace());
+   OperatorPtr A;
+   Vector B, X;
+
+   x = 0.0;
+
+   switch (integrator_type)
+   {
+      case mass:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         lfi = new DomainLFIntegrator(*scalar_u);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case diffusion:
+         scalar_u = new FunctionCoefficient(diffusion_solution);
+         scalar_f = new FunctionCoefficient(diffusion_source);
+         lfi = new DomainLFIntegrator(*scalar_f);
+         x.ProjectBdrCoefficient(*scalar_u, ess_bdr);
+         break;
+      case elasticity:
+         vector_u = new VectorFunctionCoefficient(space_dim, elasticity_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, elasticity_source);
+         lfi = new VectorDomainLFIntegrator(*vector_f);
+         x.ProjectBdrCoefficient(*vector_u, ess_bdr);
+         break;
+      case maxwell:
+         vector_u = new VectorFunctionCoefficient(space_dim, maxwell_solution);
+         vector_f = new VectorFunctionCoefficient(space_dim, maxwell_source);
+         lfi = new VectorFEDomainLFIntegrator(*vector_f);
+         x.ProjectBdrCoefficientTangent(*vector_u, ess_bdr);
+         break;
+      default:
+         mfem_error("Invalid integrator type! Check ParLinearForm");
+   }
+   b->AddDomainIntegrator(lfi);
+   b->Assemble();
+
+   /// 9. Define a geometric multigrid solver. The bilinear form
+   ///    a(.,.) is assembled internally. Set up the type of cycles
+   ///    and form the linear system.
+   auto mg = new LpqGeometricMultigrid(*fes_hierarchy,
+                                       ess_bdr,
+                                       integrator_type,
+                                       solver_type,
+                                       p_order,
+                                       q_order);
+   mg->SetCycleType(Multigrid::CycleType::VCYCLE, 1, 1);
+   mg->FormFineLinearSystem(x, *b, A, X, B);
+
+   Solver *solver = nullptr;
+   DataMonitor *monitor = nullptr;
+
+   switch (solver_type)
+   {
+      case sli:
+         solver = new SLISolver(MPI_COMM_WORLD);
+         break;
+      case cg:
+         solver = new CGSolver(MPI_COMM_WORLD);
+         break;
+      default:
+         mfem_error("Invalid solver type!");
+   }
+   solver->SetOperator(*A.Ptr());
+
+   IterativeSolver *it_solver = dynamic_cast<IterativeSolver*>(solver);
+   if (it_solver)
+   {
+      it_solver->SetRelTol(rel_tol);
+      it_solver->SetMaxIter(max_iter);
+      it_solver->SetPrintLevel(1);
+      it_solver->SetPreconditioner(*mg);
+      if (use_monitor)
+      {
+         monitor = new DataMonitor(file_name.str(), MONITOR_DIGITS);
+         it_solver->SetMonitor(*monitor);
+      }
+   }
+
+   solver->Mult(B, X);
+
+   /// 10. Recover the solution x as a grid function. Send the data by socket
+   ///     to a GLVis server.
+   mg->RecoverFineFEMSolution(X, *b, x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *fes_hierarchy->GetFinestFESpace().GetParMesh()
+               << x << flush;
+   }
+
+   /// 11. Compute and print the L^2 norm of the error
+   {
+      real_t error = 0.0;
+      switch (integrator_type)
+      {
+         case mass: case diffusion:
+            error = x.ComputeL2Error(*scalar_u);
+            break;
+         case elasticity: case maxwell:
+            error = x.ComputeL2Error(*vector_u);
+            break;
+         default:
+            mfem_error("Invalid integrator type! Check ComputeL2Error");
+      }
+      if (Mpi::Root())
+      {
+         mfem::out << "\n|| u_h - u ||_{L^2} = " << error << "\n" << endl;
+      }
+   }
+
+   /// 12. Free the memory used
+   delete mg;
+   delete solver;
+   delete b;
+   if (monitor) { delete monitor; }
+   if (scalar_u) { delete scalar_u; }
+   if (scalar_f) { delete scalar_f; }
+   if (vector_u) { delete vector_u; }
+   if (vector_f) { delete vector_f; }
+   for (int level = 0; level < fec_array.Size(); ++level)
+   {
+      delete fec_array[level];
+   }
+   delete fes_hierarchy;
+
+   return 0;
+}

tests/unit/linalg/test_matrix_hypre.cpp

@@ -243,6 +243,93 @@ TEST_CASE("HypreParMatrixAbsMult",  "[Parallel], [HypreParMatrixAbsMult]")
    }
 }
 
+TEST_CASE("HypreParMatrixPowAbsMult",  "[Parallel], [HypreParMatrixPowAbsMult]")
+{
+   int rank;
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+   int dim = 2;
+   int ne = 4;
+   double power = 2.0;
+   for (int order = 1; order <= 3; ++order)
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(
+                     ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
+      ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
+      mesh.Clear();
+      FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
+      FiniteElementCollection *l2_coll(new L2_FECollection(order, dim));
+      ParFiniteElementSpace R_space(pmesh, hdiv_coll);
+      ParFiniteElementSpace W_space(pmesh, l2_coll);
+
+      int n = R_space.GetTrueVSize();
+      int m = W_space.GetTrueVSize();
+      ParMixedBilinearForm a(&R_space, &W_space);
+      a.AddDomainIntegrator(new VectorFEDivergenceIntegrator);
+      a.Assemble();
+      a.Finalize();
+
+      HypreParMatrix *A = a.ParallelAssemble();
+      HypreParMatrix *A_pow_abs = new HypreParMatrix(*A);
+
+      hypre_ParCSRMatrix * AparCSR = *A_pow_abs;
+      A_pow_abs->HypreReadWrite();
+
+      int nnzd = AparCSR->diag->num_nonzeros;
+      real_t *d_diag_data = AparCSR->diag->data;
+      mfem::hypre_forall(nnzd, [=] MFEM_HOST_DEVICE (int i)
+      {
+         d_diag_data[i] = std::pow(fabs(d_diag_data[i]), power);
+      });
+
+      int nnzoffd = AparCSR->offd->num_nonzeros;
+      real_t *d_offd_data = AparCSR->offd->data;
+      mfem::hypre_forall(nnzoffd, [=] MFEM_HOST_DEVICE (int i)
+      {
+         d_offd_data[i] = std::pow(fabs(d_offd_data[i]), power);
+      });
+
+      Vector X0(n), X1(n);
+      Vector Y0(m), Y1(m);
+
+      X0.Randomize();
+      Y0.Randomize(1);
+      Y1.Randomize(1);
+      A->PowAbsMult(power,3.4,X0,-2.3,Y0);
+      A_pow_abs->Mult(3.4,X0,-2.3,Y1);
+
+      Y1 -= Y0;
+      double error = Y1.Norml2();
+
+      mfem::out << "Testing PowAbsMult:   order: " << order
+                << ", error norm on rank "
+                << rank << ": " << error << std::endl;
+
+      REQUIRE(error == MFEM_Approx(0.0));
+
+      MPI_Barrier(MPI_COMM_WORLD);
+
+      Y0.Randomize();
+      X0.Randomize(1);
+      X1.Randomize(1);
+      A->PowAbsMultTranspose(power,3.4,Y0,-2.3,X0);
+      A_pow_abs->MultTranspose(3.4,Y0,-2.3,X1);
+      X1 -= X0;
+
+      error = X1.Norml1();
+      mfem::out << "Testing PowAbsMultT:  order: " << order
+                << ", error norm on rank "
+                << rank << ": " << error << std::endl;
+
+      REQUIRE(error == MFEM_Approx(0.0));
+
+      delete A;
+      delete A_pow_abs;
+      delete hdiv_coll;
+      delete l2_coll;
+      delete pmesh;
+   }
+}
+
 #endif // MFEM_USE_MPI
 
 } // namespace mfem

tests/unit/linalg/test_matrix_sparse.cpp

@@ -78,6 +78,69 @@ TEST_CASE("SparseMatrixAbsMult", "[SparseMatrixAbsMult]")
    }
 }
 
+TEST_CASE("SparseMatrixPowAbsMult", "[SparseMatrixPowAbsMult]")
+{
+   int dim = 2;
+   int ne = 4;
+   real_t power = 1.5;
+   for (int order = 1; order <= 3; ++order)
+   {
+      CAPTURE(order);
+
+      Mesh mesh = Mesh::MakeCartesian2D(
+                     ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
+      FiniteElementCollection *hdiv_coll(new RT_FECollection(order, dim));
+      FiniteElementCollection *l2_coll(new L2_FECollection(order, dim));
+      FiniteElementSpace R_space(&mesh, hdiv_coll);
+      FiniteElementSpace W_space(&mesh, l2_coll);
+
+      int n = R_space.GetTrueVSize();
+      int m = W_space.GetTrueVSize();
+      MixedBilinearForm a(&R_space, &W_space);
+      a.AddDomainIntegrator(new VectorFEDivergenceIntegrator);
+      a.Assemble();
+      a.Finalize();
+
+      SparseMatrix &A = a.SpMat();
+      SparseMatrix *Aabs = new SparseMatrix(A);
+
+      int nnz = Aabs->NumNonZeroElems();
+      for (int j = 0; j < nnz; j++)
+      {
+         Aabs->GetData()[j] = std::pow(fabs(Aabs->GetData()[j]), power);
+      }
+
+      Vector X0(n), X1(n);
+      Vector Y0(m), Y1(m);
+
+      X0.Randomize();
+      Y0.Randomize(1);
+      Y1.Randomize(1);
+      A.PowAbsMult(power,X0,Y0);
+      Aabs->Mult(X0,Y1);
+
+      Y1 -= Y0;
+      double error = Y1.Norml2();
+
+      REQUIRE(error == MFEM_Approx(0.0));
+
+      Y0.Randomize();
+      X0.Randomize(1);
+      X1.Randomize(1);
+      A.PowAbsMultTranspose(power,Y0,X0);
+      Aabs->MultTranspose(Y0,X1);
+      X1 -= X0;
+
+      error = X1.Norml2();
+
+      REQUIRE(error == MFEM_Approx(0.0));
+
+      delete Aabs;
+      delete hdiv_coll;
+      delete l2_coll;
+   }
+}
+
 TEST_CASE("SparseMatrix printing", "[SparseMatrix]")
 {
    // Create a test sparse matrix and print it using different methods