Formatting.

Grouping Mesh methods in Doxygen [mesh-group-doc-dev]

CHANGELOG

@@ -41,6 +41,8 @@ Meshing improvements
 - Added new methods in the Mesh class to set and get attributes on NURBS patches
   and patch boundaries.
 
+- Added HIP support to the SUNDIALS interface.
+
 - TMOP improvement: added asymptotically-balanced compound metrics 90, 94, 328,
   338. Added the tmop-metric-magnitude tool for tracking how metrics change
   under geometric perturbations.
@@ -54,6 +56,9 @@ Discretization improvements
 
 - Added support for p-refined meshes in FindPointsGSLIB.
 
+- Added radial basis functions and reproducing kernels as alternatives to 
+  polynomial functions. These are accessed through the LocalKernelFECollection.
+
 Linear and nonlinear solvers
 ----------------------------
 - Updated interface to MUMPS direct solver to support multiple right-hand

CMakeLists.txt

@@ -317,6 +317,9 @@ if (MFEM_USE_SUNDIALS)
   if (MFEM_USE_CUDA)
     list(APPEND SUNDIALS_COMPONENTS NVector_Cuda)
   endif()
+  if (MFEM_USE_HIP)
+    list(APPEND SUNDIALS_COMPONENTS NVector_Hip)
+  endif()
   find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
 endif()
 

INSTALL

@@ -628,9 +628,13 @@ The specific libraries and their options are:
   both MPI and hypre.
   If MFEM_USE_CUDA is enabled, we expect that SUNDIALS is built with support
   for CUDA.
-  URL: http://computation.llnl.gov/projects/sundials/sundials-software
+  If MFEM_USE_HIP is enabled, we expect that SUNDIALS is built with support
+  for HIP.
+  URL: http://computing.llnl.gov/projects/sundials/sundials-software
   Options: SUNDIALS_OPT, SUNDIALS_LIB.
-  Versions: SUNDIALS >= 5.0.0, SUNDIALS >= 5.4.0 for CUDA support.
+  Versions: SUNDIALS >= 5.0.0,
+            SUNDIALS >= 5.4.0 for CUDA support, and
+            SUNDIALS >= 5.7.0 for HIP support.
 
 - SuiteSparse (optional), used when MFEM_USE_SUITESPARSE = YES.
   URL: http://faculty.cse.tamu.edu/davis/suitesparse.html

config/cmake/MFEMConfig.cmake.in

@@ -55,6 +55,8 @@ set(MFEM_USE_SIMD @MFEM_USE_SIMD@)
 set(MFEM_USE_ADIOS2 @MFEM_USE_ADIOS2@)
 set(MFEM_USE_MOONOLITH @MFEM_USE_MOONOLITH@)
 set(MFEM_USE_CODIPACK @MFEM_USE_CODIPACK@)
+set(MFEM_USE_MKL_CPARDISO @MFEM_USE_MKL_CPARDISO@)
+set(MFEM_USE_MKL_PARDISO @MFEM_USE_MKL_PARDISO@)
 set(MFEM_USE_ADFORWARD @MFEM_USE_ADFORWARD@)
 set(MFEM_USE_CALIPER @MFEM_USE_CALIPER@)
 set(MFEM_USE_ALGOIM @MFEM_USE_ALGOIM@)

config/cmake/config.hpp.in

@@ -186,6 +186,9 @@
 // Enable interface to the MKL CPardiso library.
 #cmakedefine MFEM_USE_MKL_CPARDISO
 
+// Enable interface to the MKL Pardiso library.
+#cmakedefine MFEM_USE_MKL_PARDISO
+
 // Use forward mode for automatic differentiation.
 #cmakedefine MFEM_USE_ADFORWARD
 

config/cmake/modules/FindSUNDIALS.cmake

@@ -22,8 +22,8 @@ mfem_find_package(SUNDIALS SUNDIALS SUNDIALS_DIR
     "include" nvector/nvector_serial.h "lib" sundials_nvecserial
   ADD_COMPONENT NVector_Cuda
     "include" nvector/nvector_cuda.h "lib" sundials_nveccuda
-  ADD_COMPONENT NVector_ParHyp
-    "include" nvector/nvector_parhyp.h "lib" sundials_nvecparhyp
+  ADD_COMPONENT NVector_Hip
+    "include" nvector/nvector_hip.h "lib" sundials_nvechip
   ADD_COMPONENT NVector_Parallel
     "include" nvector/nvector_parallel.h "lib" sundials_nvecparallel
   ADD_COMPONENT NVector_MPIPlusX

config/defaults.mk

@@ -267,6 +267,9 @@ endif
 ifeq ($(MFEM_USE_CUDA),YES)
    SUNDIALS_LIB += -lsundials_nveccuda
 endif
+ifeq ($(MFEM_USE_HIP),YES)
+   SUNDIALS_LIB += -lsundials_nvechip
+endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES
 

examples/ex14.cpp

@@ -18,6 +18,8 @@
 //               ex14 -m ../data/amr-quad.mesh -r 3
 //               ex14 -m ../data/amr-hex.mesh
 //               ex14 -m ../data/fichera-amr.mesh
+//               ex14 -m ../data/inline-quad.mesh -r 2 -rk
+//               ex14 -m ../data/star.mesh -r 0 -o 3 -rk
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
@@ -43,10 +45,11 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/star.mesh";
    int ref_levels = -1;
    int order = 1;
+   bool rk = false;
    double sigma = -1.0;
    double kappa = -1.0;
    double eta = 0.0;
-   bool visualization = 1;
+   bool visualization = true;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -54,7 +57,9 @@ int main(int argc, char *argv[])
    args.AddOption(&ref_levels, "-r", "--refine",
                   "Number of times to refine the mesh uniformly, -1 for auto.");
    args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) >= 0.");
+                  "Finite element order (polynomial degree) >= 0 OR reproducing kernel order.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&sigma, "-s", "--sigma",
                   "One of the three DG penalty parameters, typically +1/-1."
                   " See the documentation of class DGDiffusionIntegrator.");
@@ -75,6 +80,10 @@ int main(int argc, char *argv[])
    {
       kappa = (order+1)*(order+1);
    }
+   if (rk && sigma < 0.0)
+   {
+      sigma = 1.0;
+   }
    args.PrintOptions(cout);
 
    // 2. Read the mesh from the given mesh file. We can handle triangular,
@@ -104,7 +113,15 @@ int main(int argc, char *argv[])
 
    // 4. Define a finite element space on the mesh. Here we use discontinuous
    //    finite elements of the specified order >= 0.
-   FiniteElementCollection *fec = new DG_FECollection(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 4, 6, order, 1.01 + order, 0.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
    cout << "Number of unknowns: " << fespace->GetVSize() << endl;
 

examples/ex14p.cpp

@@ -17,6 +17,8 @@
 //               mpirun -np 4 ex14p -m ../data/inline-segment.mesh -rs 5
 //               mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 3
 //               mpirun -np 4 ex14p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex14p -m ../data/star.mesh -rs 1 -o 2 -rk
+//               mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 2 -rk
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
@@ -81,6 +83,7 @@ int main(int argc, char *argv[])
    int ser_ref_levels = -1;
    int par_ref_levels = 2;
    int order = 1;
+   bool rk = false;
    double sigma = -1.0;
    double kappa = -1.0;
    double eta = 0.0;
@@ -96,6 +99,8 @@ int main(int argc, char *argv[])
                   "Number of times to refine the mesh uniformly in parallel.");
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&sigma, "-s", "--sigma",
                   "One of the three DG penalty parameters, typically +1/-1."
                   " See the documentation of class DGDiffusionIntegrator.");
@@ -119,6 +124,10 @@ int main(int argc, char *argv[])
    {
       kappa = (order+1)*(order+1);
    }
+   if (rk && sigma < 0.0)
+   {
+      sigma = 1.0;
+   }
    if (myid == 0)
    {
       args.PrintOptions(cout);
@@ -163,7 +172,15 @@ int main(int argc, char *argv[])
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use discontinuous finite elements of the specified order >= 0.
-   FiniteElementCollection *fec = new DG_FECollection(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 4, 6, order, 1.01 + order, 0.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
    HYPRE_BigInt size = fespace->GlobalTrueVSize();
    if (myid == 0)

examples/ex18.cpp

@@ -9,6 +9,8 @@
 //       ex18 -p 1 -r 0 -o 5 -s 6
 //       ex18 -p 2 -r 1 -o 1 -s 3
 //       ex18 -p 2 -r 0 -o 3 -s 3
+//       ex18 -p 1 -r 1 -o 3 -s 4 -rk
+//       ex18 -p 2 -r 0 -o 3 -s 3 -rk
 //
 // Description:  This example code solves the compressible Euler system of
 //               equations, a model nonlinear hyperbolic PDE, with a
@@ -65,6 +67,7 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/periodic-square.mesh";
    int ref_levels = 1;
    int order = 3;
+   bool rk = false;
    int ode_solver_type = 4;
    double t_final = 2.0;
    double dt = -0.01;
@@ -84,6 +87,8 @@ int main(int argc, char *argv[])
                   "Number of times to refine the mesh uniformly.");
    args.AddOption(&order, "-o", "--order",
                   "Order (degree) of the finite elements.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                   "ODE solver: 1 - Forward Euler,\n\t"
                   "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
@@ -139,13 +144,21 @@ int main(int argc, char *argv[])
 
    // 5. Define the discontinuous DG finite element space of the given
    //    polynomial order on the refined mesh.
-   DG_FECollection fec(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 5, 7, order, 2.01 + order, 1.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    // Finite element space for a scalar (thermodynamic quantity)
-   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace fes(&mesh, fec);
    // Finite element space for a mesh-dim vector quantity (momentum)
-   FiniteElementSpace dfes(&mesh, &fec, dim, Ordering::byNODES);
+   FiniteElementSpace dfes(&mesh, fec, dim, Ordering::byNODES);
    // Finite element space for all variables together (total thermodynamic state)
-   FiniteElementSpace vfes(&mesh, &fec, num_equation, Ordering::byNODES);
+   FiniteElementSpace vfes(&mesh, fec, num_equation, Ordering::byNODES);
 
    // This example depends on this ordering of the space.
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
@@ -304,6 +317,7 @@ int main(int argc, char *argv[])
 
    // Free the used memory.
    delete ode_solver;
+   delete fec;
 
    return 0;
 }

examples/ex18p.cpp

@@ -9,6 +9,8 @@
 //       mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 5 -s 6
 //       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 1 -s 3
 //       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 3 -s 3
+//       mpirun -np 4 ex18p -p 1 -rs 1 -rp 0 -o 3 -s 4 -rk
+//       mpirun -np 4 ex18p -p 2 -rs 1 -rp 0 -o 3 -s 3 -rk
 //
 // Description:  This example code solves the compressible Euler system of
 //               equations, a model nonlinear hyperbolic PDE, with a
@@ -70,6 +72,7 @@ int main(int argc, char *argv[])
    int ser_ref_levels = 0;
    int par_ref_levels = 1;
    int order = 3;
+   bool rk = false;
    int ode_solver_type = 4;
    double t_final = 2.0;
    double dt = -0.01;
@@ -93,6 +96,8 @@ int main(int argc, char *argv[])
                   " partitioning.");
    args.AddOption(&order, "-o", "--order",
                   "Order (degree) of the finite elements.");
+   args.AddOption(&rk, "-rk", "--rk", "-no-rk", "--no-rk",
+                  "Use reproducing kernel functions");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                   "ODE solver: 1 - Forward Euler,\n\t"
                   "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
@@ -161,13 +166,21 @@ int main(int argc, char *argv[])
 
    // 7. Define the discontinuous DG finite element space of the given
    //    polynomial order on the refined mesh.
-   DG_FECollection fec(order, dim);
+   FiniteElementCollection *fec;
+   if (rk)
+   {
+      fec = new LocalKernelFECollection(dim, 5, 7, order, 2.01 + order, 1.0);
+   }
+   else
+   {
+      fec = new DG_FECollection(order, dim);
+   }
    // Finite element space for a scalar (thermodynamic quantity)
-   ParFiniteElementSpace fes(&pmesh, &fec);
+   ParFiniteElementSpace fes(&pmesh, fec);
    // Finite element space for a mesh-dim vector quantity (momentum)
-   ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
+   ParFiniteElementSpace dfes(&pmesh, fec, dim, Ordering::byNODES);
    // Finite element space for all variables together (total thermodynamic state)
-   ParFiniteElementSpace vfes(&pmesh, &fec, num_equation, Ordering::byNODES);
+   ParFiniteElementSpace vfes(&pmesh, fec, num_equation, Ordering::byNODES);
 
    // This example depends on this ordering of the space.
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
@@ -375,6 +388,7 @@ int main(int argc, char *argv[])
 
    // Free the used memory.
    delete ode_solver;
+   delete fec;
 
    return 0;
 }

examples/ex9p.cpp

@@ -536,8 +536,10 @@ int main(int argc, char *argv[])
       if (!sout)
       {
          if (Mpi::Root())
+         {
             cout << "Unable to connect to GLVis server at "
                  << vishost << ':' << visport << endl;
+         }
          visualization = false;
          if (Mpi::Root())
          {
@@ -552,8 +554,10 @@ int main(int argc, char *argv[])
          sout << "pause\n";
          sout << flush;
          if (Mpi::Root())
+         {
             cout << "GLVis visualization paused."
                  << " Press space (in the GLVis window) to resume it.\n";
+         }
       }
    }
 

examples/sundials/CMakeLists.txt

@@ -68,11 +68,43 @@ if (MFEM_ENABLE_TESTING)
       add_test(NAME ${TEST_NAME}_ser
         COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
     else()
-      add_test(NAME ${TEST_NAME}_np=4
+      add_test(NAME ${TEST_NAME}_np=${MFEM_MPI_NP}
         COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
         ${MPIEXEC_PREFLAGS}
         $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
         ${MPIEXEC_POSTFLAGS})
     endif()
   endforeach()
-endif()
+
+  # Add CUDA/HIP tests.
+  set(DEVICE_EXAMPLES
+      # serial examples with device support:
+      ex9
+      # parallel examples with device support:
+      ex9p)
+  set(MFEM_TEST_DEVICE)
+  if (MFEM_USE_CUDA)
+    set(MFEM_TEST_DEVICE "cuda")
+  elseif (MFEM_USE_HIP)
+    set(MFEM_TEST_DEVICE "hip")
+  endif()
+  if (MFEM_TEST_DEVICE)
+    foreach(TEST_NAME ${DEVICE_EXAMPLES})
+      string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
+
+      set(THIS_TEST_OPTIONS "-no-vis" "-d" "${MFEM_TEST_DEVICE}")
+      list(APPEND THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
+
+      if (NOT (${TEST_NAME} MATCHES ".*p$"))
+        add_test(NAME ${PFX}${TEST_NAME}_${MFEM_TEST_DEVICE}_ser
+          COMMAND ${PFX}${TEST_NAME} ${THIS_TEST_OPTIONS})
+      else()
+        add_test(NAME ${PFX}${TEST_NAME}_${MFEM_TEST_DEVICE}_np=${MFEM_MPI_NP}
+          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+          ${MPIEXEC_PREFLAGS}
+          $<TARGET_FILE:${PFX}${TEST_NAME}> ${THIS_TEST_OPTIONS}
+          ${MPIEXEC_POSTFLAGS})
+      endif()
+    endforeach()
+  endif(MFEM_TEST_DEVICE)
+endif(MFEM_ENABLE_TESTING)

examples/sundials/README

@@ -12,8 +12,7 @@ use of MFEM features based on the SUNDIALS suite of time integration and
 non-linear solvers.
 
 To build these examples, make sure that MFEM is configured with the option
-"MFEM_USE_SUNDIALS = YES", see the top-level INSTALL file for details (version
-2.7 or higher of SUNDIALS is required).
+"MFEM_USE_SUNDIALS = YES", see the top-level INSTALL file for details.
 
 We recommend comparing the original example codes with the corresponding files
 in the current directory.

examples/sundials/ex9.cpp

@@ -280,15 +280,16 @@ int main(int argc, char *argv[])
       k.SetAssemblyLevel(AssemblyLevel::FULL);
    }
    m.AddDomainIntegrator(new MassIntegrator);
-   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
    k.AddInteriorFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
    k.AddBdrFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
 
    LinearForm b(&fes);
    b.AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
+      new BoundaryFlowIntegrator(inflow, velocity, alpha));
 
    m.Assemble();
    int skip_zeros = 0;

examples/sundials/ex9p.cpp

@@ -63,6 +63,66 @@ double inflow_function(const Vector &x);
 // Mesh bounding box
 Vector bb_min, bb_max;
 
+// Type of preconditioner for implicit time integrator
+enum class PrecType : int
+{
+   ILU = 0,
+   AIR = 1
+};
+
+#if MFEM_HYPRE_VERSION >= 21800
+// Algebraic multigrid preconditioner for advective problems based on
+// approximate ideal restriction (AIR). Most effective when matrix is
+// first scaled by DG block inverse, and AIR applied to scaled matrix.
+// See https://doi.org/10.1137/17M1144350.
+class AIR_prec : public Solver
+{
+private:
+   const HypreParMatrix *A;
+   // Copy of A scaled by block-diagonal inverse
+   HypreParMatrix A_s;
+
+   HypreBoomerAMG *AIR_solver;
+   int blocksize;
+
+public:
+   AIR_prec(int blocksize_) : AIR_solver(NULL), blocksize(blocksize_) { }
+
+   void SetOperator(const Operator &op)
+   {
+      width = op.Width();
+      height = op.Height();
+
+      A = dynamic_cast<const HypreParMatrix *>(&op);
+      MFEM_VERIFY(A != NULL, "AIR_prec requires a HypreParMatrix.")
+
+      // Scale A by block-diagonal inverse
+      BlockInverseScale(A, &A_s, NULL, NULL, blocksize,
+                        BlockInverseScaleJob::MATRIX_ONLY);
+      delete AIR_solver;
+      AIR_solver = new HypreBoomerAMG(A_s);
+      AIR_solver->SetAdvectiveOptions(1, "", "FA");
+      AIR_solver->SetPrintLevel(0);
+      AIR_solver->SetMaxLevels(50);
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      // Scale the rhs by block inverse and solve system
+      HypreParVector z_s;
+      BlockInverseScale(A, NULL, &x, &z_s, blocksize,
+                        BlockInverseScaleJob::RHS_ONLY);
+      AIR_solver->Mult(z_s, y);
+   }
+
+   ~AIR_prec()
+   {
+      delete AIR_solver;
+   }
+};
+#endif
+
+
 class DG_Solver : public Solver
 {
 private:
@@ -70,24 +130,37 @@ private:
    SparseMatrix M_diag;
    HypreParMatrix *A;
    GMRESSolver linear_solver;
-   BlockILU prec;
+   Solver *prec;
    double dt;
 public:
-   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
+             PrecType prec_type)
       : M(M_),
         K(K_),
         A(NULL),
         linear_solver(M.GetComm()),
-        prec(fes.GetFE(0)->GetDof(),
-             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {
+      int block_size = fes.GetFE(0)->GetDof();
+      if (prec_type == PrecType::ILU)
+      {
+         prec = new BlockILU(block_size,
+                             BlockILU::Reordering::MINIMUM_DISCARDED_FILL);
+      }
+      else if (prec_type == PrecType::AIR)
+      {
+#if MFEM_HYPRE_VERSION >= 21800
+         prec = new AIR_prec(block_size);
+#else
+         MFEM_ABORT("Must have MFEM_HYPRE_VERSION >= 21800 to use AIR.\n");
+#endif
+      }
       linear_solver.iterative_mode = false;
       linear_solver.SetRelTol(1e-9);
       linear_solver.SetAbsTol(0.0);
       linear_solver.SetMaxIter(100);
       linear_solver.SetPrintLevel(0);
-      linear_solver.SetPreconditioner(prec);
+      linear_solver.SetPreconditioner(*prec);
 
       M.GetDiag(M_diag);
    }
@@ -120,10 +193,12 @@ public:
 
    ~DG_Solver()
    {
+      delete prec;
       delete A;
    }
 };
 
+
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -141,7 +216,8 @@ private:
    mutable Vector z;
 
 public:
-   FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, const Vector &b_);
+   FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, const Vector &b_,
+                PrecType prec_type);
 
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
@@ -178,6 +254,11 @@ int main(int argc, char *argv[])
    bool adios2 = false;
    bool binary = false;
    int vis_steps = 5;
+#if MFEM_HYPRE_VERSION >= 21800
+   PrecType prec_type = PrecType::AIR;
+#else
+   PrecType prec_type = PrecType::ILU;
+#endif
 
    // Relative and absolute tolerances for CVODE and ARKODE.
    const double reltol = 1e-2, abstol = 1e-2;
@@ -218,6 +299,8 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
+                  "implicit solves. 0 for ILU, 1 for pAIR-AMG.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -238,13 +321,13 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      if (myid == 0)
+      if (Mpi::Root())
       {
          args.PrintUsage(cout);
       }
       return 1;
    }
-   if (myid == 0)
+   if (Mpi::Root())
    {
       args.PrintOptions(cout);
    }
@@ -252,7 +335,7 @@ int main(int argc, char *argv[])
    // check for valid ODE solver option
    if (ode_solver_type < 1 || ode_solver_type > 9)
    {
-      if (myid == 0)
+      if (Mpi::Root())
       {
          cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
       }
@@ -260,7 +343,7 @@ int main(int argc, char *argv[])
    }
 
    Device device(device_config);
-   if (myid == 0) { device.Print(); }
+   if (Mpi::Root()) { device.Print(); }
 
    // 3. Read the serial mesh from the given mesh file on all processors. We can
    //    handle geometrically periodic meshes in this code.
@@ -297,7 +380,7 @@ int main(int argc, char *argv[])
    ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
 
    HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
-   if (myid == 0)
+   if (Mpi::Root())
    {
       cout << "Number of unknowns: " << global_vSize << endl;
    }
@@ -328,15 +411,16 @@ int main(int argc, char *argv[])
    }
 
    m->AddDomainIntegrator(new MassIntegrator);
-   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   constexpr double alpha = -1.0;
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
    k->AddInteriorFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
    k->AddBdrFaceIntegrator(
-      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
 
    ParLinearForm *b = new ParLinearForm(fes);
    b->AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
+      new BoundaryFlowIntegrator(inflow, velocity, alpha));
 
    int skip_zeros = 0;
    m->Assemble();
@@ -435,11 +519,13 @@ int main(int argc, char *argv[])
       sout.open(vishost, visport);
       if (!sout)
       {
-         if (myid == 0)
+         if (Mpi::Root())
+         {
             cout << "Unable to connect to GLVis server at "
                  << vishost << ':' << visport << endl;
+         }
          visualization = false;
-         if (myid == 0)
+         if (Mpi::Root())
          {
             cout << "GLVis visualization disabled.\n";
          }
@@ -451,15 +537,17 @@ int main(int argc, char *argv[])
          sout << "solution\n" << *pmesh << *u;
          sout << "pause\n";
          sout << flush;
-         if (myid == 0)
+         if (Mpi::Root())
+         {
             cout << "GLVis visualization paused."
                  << " Press space (in the GLVis window) to resume it.\n";
+         }
       }
    }
 
    // 9. Define the time-dependent evolution operator describing the ODE
    //    right-hand side, and define the ODE solver used for time integration.
-   FE_Evolution adv(*m, *k, *B);
+   FE_Evolution adv(*m, *k, *B, prec_type);
 
    double t = 0.0;
    adv.SetTime(t);
@@ -511,7 +599,7 @@ int main(int argc, char *argv[])
 
       if (done || ti % vis_steps == 0)
       {
-         if (myid == 0)
+         if (Mpi::Root())
          {
             cout << "time step: " << ti << ", time: " << t << endl;
             if (cvode) { cvode->PrintInfo(); }
@@ -590,7 +678,7 @@ int main(int argc, char *argv[])
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
-                           const Vector &b_)
+                           const Vector &b_, PrecType prec_type)
    : TimeDependentOperator(M_.Height()),
      b(b_),
      M_solver(M_.ParFESpace()->GetComm()),
@@ -617,7 +705,7 @@ FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
       HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
       M_prec = hypre_prec;
 
-      dg_solver = new DG_Solver(M_mat, K_mat, *M_.FESpace());
+      dg_solver = new DG_Solver(M_mat, K_mat, *M_.FESpace(), prec_type);
    }
    else
    {
@@ -633,6 +721,10 @@ FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
    M_solver.SetPrintLevel(0);
 }
 
+// Solve the equation:
+//    u_t = M^{-1}(Ku + b),
+// by solving associated linear system
+//    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
 {
    K->Mult(x, z);

examples/sundials/makefile

@@ -23,6 +23,8 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex9 ex10 ex16
 PAR_EXAMPLES = ex9p ex10p ex16p
+SEQ_DEVICE_EXAMPLES = ex9
+PAR_DEVICE_EXAMPLES = ex9p
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
 else
@@ -54,10 +56,22 @@ include $(MFEM_TEST_MK)
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 SERIAL_NAME := Serial SUNDIALS example
 PARALLEL_NAME := Parallel SUNDIALS example
+SERIAL_CUDA_NAME := Serial SUNDIALS CUDA example
+PARALLEL_CUDA_NAME := Parallel SUNDIALS CUDA example
+SERIAL_HIP_NAME := Serial SUNDIALS HIP example
+PARALLEL_HIP_NAME := Parallel SUNDIALS HIP example
 %-test-par: %
 	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME))
 %-test-seq: %
 	@$(call mfem-test,$<,, $(SERIAL_NAME))
+%-test-par-cuda: %
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_CUDA_NAME),-d cuda)
+%-test-seq-cuda: %
+	@$(call mfem-test,$<,, $(SERIAL_CUDA_NAME),-d cuda)
+%-test-par-hip: %
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_HIP_NAME),-d hip)
+%-test-seq-hip: %
+	@$(call mfem-test,$<,, $(SERIAL_HIP_NAME),-d hip)
 
 # Testing: Specific execution options:
 # Example 9: test CVODE with CV_ADAMS (non-stiff implicit) time stepping
@@ -68,6 +82,16 @@ ex9-test-seq: ex9
 	@$(call mfem-test,$<,, $(SERIAL_NAME),$(EX9_ARGS))
 ex9p-test-par: ex9p
 	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_NAME),$(EX9P_ARGS))
+ex9-test-seq-cuda: ex9
+	@$(call mfem-test,$<,, $(SERIAL_CUDA_NAME),-d cuda $(EX9_ARGS))
+ex9p-test-par-cuda: ex9p
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_CUDA_NAME),-d cuda \
+	   $(EX9P_ARGS))
+ex9-test-seq-hip: ex9
+	@$(call mfem-test,$<,, $(SERIAL_HIP_NAME),-d hip $(EX9_ARGS))
+ex9p-test-par-hip: ex9p
+	@$(call mfem-test,$<, $(RUN_MPI), $(PARALLEL_HIP_NAME),-d hip \
+	   $(EX9P_ARGS))
 # Example 10: test CVODE with CV_BDF (stiff implicit) time stepping
 EX10_COMMON_ARGS := -m ../../data/beam-quad.mesh -o 2 -s 5 -dt 0.15 -tf 6 -vs 10
 EX10_ARGS  := $(EX10_COMMON_ARGS) -r  2

fem/CMakeLists.txt

@@ -68,6 +68,7 @@ set(SRCS
   fe/fe_nd.cpp
   fe/fe_nurbs.cpp
   fe/fe_pos.cpp
+  fe/fe_rbf.cpp
   fe/fe_rt.cpp
   fe/fe_ser.cpp
   fe_coll.cpp
@@ -175,6 +176,7 @@ set(HDRS
   fe/fe_nd.hpp
   fe/fe_nurbs.hpp
   fe/fe_pos.hpp
+  fe/fe_rbf.hpp
   fe/fe_rt.hpp
   fe/fe_ser.hpp
   fe_coll.hpp

fem/fe.hpp

@@ -26,5 +26,6 @@
 #include "fe/fe_nurbs.hpp"
 #include "fe/fe_pos.hpp"
 #include "fe/fe_ser.hpp"
+#include "fe/fe_rbf.hpp"
 
 #endif

fem/fe/fe_base.hpp

@@ -1026,9 +1026,6 @@ private:
 
    static Array2D<int> binom;
 
-   static void CalcMono(const int p, const double x, double *u);
-   static void CalcMono(const int p, const double x, double *u, double *d);
-
    static void CalcChebyshev(const int p, const double x, double *u);
    static void CalcChebyshev(const int p, const double x, double *u, double *d);
    static void CalcChebyshev(const int p, const double x, double *u, double *d,
@@ -1165,6 +1162,9 @@ public:
    static void CalcLegendre(const int p, const double x, double *u);
    static void CalcLegendre(const int p, const double x, double *u, double *d);
 
+   static void CalcMono(const int p, const double x, double *u);
+   static void CalcMono(const int p, const double x, double *u, double *d);
+
    ~Poly_1D();
 };
 

fem/fe/fe_rbf.hpp

@@ -0,0 +1,596 @@
+// Copyright (c) 2010-2023, 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_FE_RBF
+#define MFEM_FE_RBF
+
+#include "fe_base.hpp"
+
+namespace mfem
+{
+
+/** @brief Pure virtual class for dimensionless radial basis functions (RBFs).
+    Many RBFs are shaped like a Gaussian and are used here as an alternative
+    to polynomials in RBF and RK elements. The input for these
+    is a dimensionless radius r = ||x|| / h, where ||x|| is a distance and
+    h is the smoothing parameter, which controls the gradient of the RBF.
+    For reference: https://doi.org/10.1017/S0962492900000015 */
+class RBFKernel
+{
+public:
+   static const double GlobalRadius; // functions with r>=GR are considered global
+   RBFKernel() { };
+   virtual ~RBFKernel() { }
+
+   /// Evaluate the dimensionless RBF
+   virtual double BaseFunction(double r) const = 0;
+
+   /// Evaluate the derivative of the dimensionless RBF with respect to r
+   virtual double BaseDerivative(double r) const = 0;
+
+   /// Evaluate the second derivative of the dimensionless RBF with respect to r
+   virtual double BaseDerivative2(double r) const = 0;
+
+   /// The support radius, outside of which the function is zero if the function has compact support
+   virtual double Radius() const { return GlobalRadius; }
+
+   /// Does function have compact support?
+   virtual bool CompactSupport() const { return false; }
+
+   /** #brief This normalizes the smoothing parameter h such that h doesn't
+       need to be changed based on the choice of basis function */
+   virtual double HNorm() const = 0;
+};
+
+/// Gaussian RBF, exp(-r^2)
+class GaussianRBF : public RBFKernel
+{
+   // hNorm minimizes integral of Gaussian minus Wendland kernel over r=0,1
+   static const double hNorm;
+public:
+   GaussianRBF() { };
+   virtual ~GaussianRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double HNorm() const { return hNorm; }
+};
+
+/// Multiquadric RBF, sqrt(1+r^2)
+class MultiquadricRBF : public RBFKernel
+{
+   // Same as inverse multiquadric
+   static const double hNorm;
+public:
+   MultiquadricRBF() { };
+   virtual ~MultiquadricRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double HNorm() const { return hNorm; }
+};
+
+/// Inverse multiquadric RBF, 1/sqrt(1+r^2)
+class InvMultiquadricRBF : public RBFKernel
+{
+   // hNorm minimizes integral of Gaussian minus InvMQ kernel over r=0,0.5
+   static const double hNorm;
+public:
+   InvMultiquadricRBF() { };
+   virtual ~InvMultiquadricRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double HNorm() const { return hNorm; }
+};
+
+/** @brief Identitcal to the Gaussian RBF, but subtracted by a factor
+    such that the function is exactly zero at the chosen radius */
+class CompactGaussianRBF : public RBFKernel
+{
+   static const double hNorm;
+   const double radius;
+   double multK, shiftK;
+
+public:
+   CompactGaussianRBF(const double rad = 5.0);
+   virtual ~CompactGaussianRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return hNorm; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Identical to the Gaussian, but truncated (set to zero) at the chosen radius
+class TruncatedGaussianRBF : public RBFKernel
+{
+   static const double hNorm;
+   const double radius;
+public:
+   TruncatedGaussianRBF(const double rad = 5.0)
+      : radius(rad) { }
+   virtual ~TruncatedGaussianRBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return hNorm; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Wendland 11 RBF, (1-r)^3 * (1+3r) if r < 1
+class Wendland11RBF : public RBFKernel
+{
+   static const double radius;
+
+public:
+   Wendland11RBF() { }
+   virtual ~Wendland11RBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Wendland 31 RBF, (1-r)^4 * (1+4r) if r < 1
+class Wendland31RBF : public RBFKernel
+{
+   static const double radius;
+
+public:
+   Wendland31RBF() { };
+   virtual ~Wendland31RBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Wendland 33 RBF, (1-r)^8 * (1+8r+25r^2+32r^3) if r < 1
+class Wendland33RBF : public RBFKernel
+{
+   static const double radius;
+
+public:
+   Wendland33RBF() { };
+   virtual ~Wendland33RBF() { }
+
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+/// Class for storing and creating the various RBFs
+class RBFType
+{
+public:
+   /// Represent each type of function for input/output
+   enum
+   {
+      Gaussian = 0,
+      Multiquadric = 1,
+      InvMultiquadric = 2,
+      TruncatedGaussian = 3,
+      CompactGaussian = 4,
+      Wendland11 = 5,
+      Wendland31 = 6,
+      Wendland33 = 7,
+      NumRBFTypes = 8
+   };
+
+   /// Return the requested RBF
+   static RBFKernel *GetRBF(const int rbfType)
+   {
+      switch (rbfType)
+      {
+         case RBFType::Gaussian:
+            return new GaussianRBF();
+         case RBFType::Multiquadric:
+            return new MultiquadricRBF();
+         case RBFType::InvMultiquadric:
+            return new InvMultiquadricRBF();
+         case RBFType::TruncatedGaussian:
+            return new TruncatedGaussianRBF();
+         case RBFType::CompactGaussian:
+            return new CompactGaussianRBF();
+         case RBFType::Wendland11:
+            return new Wendland11RBF();
+         case RBFType::Wendland31:
+            return new Wendland31RBF();
+         case RBFType::Wendland33:
+            return new Wendland33RBF();
+      }
+      MFEM_ABORT("unknown RBF type");
+      return NULL;
+   }
+
+   /// Abort if rbfType is invalid
+   static int Check(const int rbfType)
+   {
+      MFEM_VERIFY(0 <= rbfType && rbfType < NumRBFTypes,
+                  "unknown RBF type: " << rbfType);
+      return rbfType;
+   }
+
+   /// Convert rbf int to identifier for storage
+   static char GetChar(const int rbfType)
+   {
+      static const char ident[] = { 'G', 'M', 'I',
+                                    'T', 'C',
+                                    '1', '3', '6'
+                                  };
+      return ident[Check(rbfType)];
+   }
+
+   /// Convert identifier to rbf int
+   static int GetType(const char rbfIdent)
+   {
+      switch (rbfIdent)
+      {
+         case 'G': return Gaussian;
+         case 'M': return Multiquadric;
+         case 'I': return InvMultiquadric;
+         case 'T': return TruncatedGaussian;
+         case 'C': return CompactGaussian;
+         case '1': return Wendland11;
+         case '3': return Wendland31;
+         case '6': return Wendland33;
+      }
+      MFEM_ABORT("unknown RBF identifier: " << rbfIdent);
+      return -1;
+   }
+};
+
+/// Dimensionless distance metrics, whose output is the input for RBFs
+class DistanceMetric
+{
+protected:
+   int dim;
+public:
+   /** @brief Create a distance metric
+       @param D    Reference space dimension */
+   DistanceMetric(int D) { dim = D; }
+   virtual ~DistanceMetric() { }
+
+   /// Set the reference dimension
+   virtual void SetDim(int D) { dim = D; }
+
+   /// Given a vector of length D, return a physical distance
+   virtual void Distance(const Vector &x,
+                         double &r) const = 0;
+
+   /** @brief Given a vector of length D, return the gradient of the
+       distance with respect to the original coordinates */
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const = 0;
+
+   /** @brief Given a vector of length D, return the Hessian of the
+       distance with respect to the original coordinates */
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const = 0;
+
+   /// Create an Lp distance metric for the requested dimension and norm
+   static DistanceMetric *GetDistance(int dim, int pnorm);
+};
+
+/// Dimensionless distance with r = |x| + |y| + ...
+class L1Distance : public DistanceMetric
+{
+public:
+   L1Distance(int D) : DistanceMetric(D) { };
+   virtual ~L1Distance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+/// Dimensionless distance with r = (x^2 + y^2 + ...)^(1/2)
+class L2Distance : public DistanceMetric
+{
+public:
+   L2Distance(int D) : DistanceMetric(D) { };
+   virtual ~L2Distance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+/// Dimensionless distance with r = (x^p + y^p + ...)^(1/p)
+class LpDistance : public DistanceMetric
+{
+   const int p;
+   const double pinv;
+public:
+   LpDistance(int D, int pnorm)
+      : DistanceMetric(D),
+        p(pnorm),
+        pinv(1. / static_cast<double>(p))
+   { };
+   virtual ~LpDistance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+
+/** @brief Pure virtual class for a finite element with radial basis functions
+    instead of polynomials inside each element */
+class KernelFiniteElement : public ScalarFiniteElement
+{
+private:
+   // Choose whether to interpolate or project when Project is called
+   bool interpolate = false;
+public:
+   KernelFiniteElement(int D, Geometry::Type G, int Do, int O, int F)
+      : ScalarFiniteElement(D, G, Do, O, F) { }
+   virtual ~KernelFiniteElement() { }
+
+   /// Converts integration rule to vector
+   virtual void IntRuleToVec(const IntegrationPoint &ip,
+                             Vector &vec) const;
+
+   /// Is base RBF compact?
+   virtual bool IsCompact() const = 0;
+
+   /// Return base kernel
+   virtual const RBFKernel *Kernel() const = 0;
+
+   /** @brief Return whether shape function is a tensor product, used in providing indices for compact support */
+   virtual bool TensorIndexed() const { return false; }
+
+   /** @brief Get range of indices (start and end for each dimension) that
+       that are nonzero for compact support for the given integration point */
+   virtual void GetTensorIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const
+   { MFEM_ABORT("GetTensorIndices(...)"); }
+
+   /** @brief Return total number of points in each dimension for tensor-
+       indexed points */
+   virtual void GetTensorNumPoints(int (&tNumPoints)[3]) const
+   { MFEM_ABORT("GetTensorNumPoints(...)"); }
+
+   using FiniteElement::Project;
+
+   virtual void Project(Coefficient &coeff, ElementTransformation &Trans,
+                        Vector &dofs) const;
+
+   virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans,
+                        Vector &dofs) const;
+
+   virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
+                        DenseMatrix &I) const;
+
+   virtual void GetLocalInterpolation(ElementTransformation &Trans,
+                                      DenseMatrix &I) const
+   { ScalarLocalInterpolation(Trans, I, *this); }
+
+   virtual void GetTransferMatrix(const FiniteElement &fe,
+                                  ElementTransformation &Trans,
+                                  DenseMatrix &I) const
+   { CheckScalarFE(fe).ScalarLocalInterpolation(Trans, I, *this); }
+};
+
+/** @brief Finite element using base radial basis functions without
+    polynomial corrections. */
+class RBFFiniteElement : public KernelFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable double r_scr, df_scr, ddf_scr;
+   mutable Vector x_scr, y_scr, dy_scr, dr_scr;
+   mutable DenseMatrix ddr_scr;
+   mutable int cInd[3][2];
+#endif
+   bool isCompact;
+   int dimPoints[3];
+   int numPointsD;
+   double delta; // Distance between points
+   double h;
+   double hPhys; // Shape parameter times distance between points times HNorm
+   double hPhysInv; // Inverse hPhys
+   double radPhys; // Radius adjusted by h
+   double faceFactor;
+   const RBFKernel *rbf;
+   const DistanceMetric *distance;
+   void InitializeGeometry();
+
+   // Get the dimensionless distance from x to the center of the RBF indexed i
+   virtual void DistanceVec(const int i,
+                            const Vector &x,
+                            Vector &y) const;
+
+public:
+   /** @brief Construct RBFFiniteElement
+       @param D           Reference space dimension
+       @param numPointsD  Number of points across the element in each dimension
+       @param rbfType     Type of radial basis function, from RBFType
+       @param distNorm    Norm used for distance, usually 2 = Euclidean distance
+       @param intOrder    Number of integration points per RBF point in each dimension
+       @param h           Shape parameter, approximately equal to the number of points in the support radius in one dimension
+       @param faceFactor  1.0 = points end on face, 0.0 = points end at dx/2 from face
+   */
+   RBFFiniteElement(const int D,
+                    const int numPointsD,
+                    const int rbfType,
+                    const int distNorm,
+                    const int intOrder,
+                    const double h,
+                    const double faceFactor);
+   virtual ~RBFFiniteElement() { delete rbf; delete distance; }
+
+   virtual bool TensorIndexed() const { return true; }
+   virtual void GetCompactIndices(const Vector &ip,
+                                  int (&indices)[3][2]) const;
+   virtual void GetGlobalIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const;
+   virtual void GetTensorIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const;
+   virtual void GetTensorNumPoints(int (&tNumPoints)[3]) const
+   {
+      tNumPoints[0] = dimPoints[0];
+      tNumPoints[1] = dimPoints[1];
+      tNumPoints[2] = dimPoints[2];
+   }
+
+   virtual bool IsCompact() const { return isCompact; }
+   virtual const RBFKernel *Kernel() const { return rbf; }
+
+   virtual void CalcShape(const IntegrationPoint &ip,
+                          Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &hess) const;
+};
+
+/** @brief Reproducing kernel finite element, which includes polynomial
+    corrections to the standard radial basis function finite element
+    to guarantee a chosen order of accuracy */
+class RKFiniteElement : public KernelFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable double f_scr;
+   mutable Vector x_scr, y_scr, g_scr, c_scr, s_scr, p_scr, df_scr;
+   mutable DenseMatrix q_scr, dq_scr, M_scr;
+   mutable Vector dc_scr[3], dp_scr[3];
+   mutable DenseMatrix dM_scr[3];
+   mutable DenseMatrixInverse Minv_scr;
+   mutable int cInd[3][2];
+   mutable int dimPoints[3];
+#endif
+   int polyOrd, numPoly, numPoly1d;
+   KernelFiniteElement *baseFE;
+
+   // Get the vector of polynomials for the corrections, evaluated at x
+   virtual void GetPoly(const Vector &x,
+                        Vector &p) const;
+   virtual void GetDPoly(const Vector &x,
+                         Vector &p,
+                         Vector (&dp)[3]) const;
+
+   /* Helper functions that return pieces of the RK evaluation
+      W_{RK,i} = P_i^T C_i W_{RBF,i}, where
+      P_i = [1, x, y, z, ...] is the polynomial vector evaluated at the point i,
+      M_i = P_i P_i^T W_{RBF,i} is a matrix used in calculating the corrections,
+      C_i = M_i^{-1} G are the RK corrections,
+      G = [1, 0, 0, ...] is a convencience vector */
+   virtual void GetG(Vector &g) const;
+   virtual void GetM(const Vector &baseShape,
+                     const IntegrationPoint &ip,
+                     DenseMatrix &M) const;
+   virtual void GetDM(const Vector &baseShape,
+                      const DenseMatrix &baseDeriv,
+                      const IntegrationPoint &ip,
+                      DenseMatrix &M,
+                      DenseMatrix (&dM)[3]) const;
+   virtual void AddToM(const Vector &p,
+                       const double &f,
+                       DenseMatrix &M) const;
+   virtual void AddToDM(const Vector &p,
+                        const Vector (&dp)[3],
+                        const double &f,
+                        const Vector &df,
+                        DenseMatrix (&dM)[3]) const;
+
+   // Given the corrections and base values, calculate the RK value at the ip
+   virtual void CalculateValues(const Vector &c,
+                                const Vector &baseShape,
+                                const IntegrationPoint &ip,
+                                Vector &shape) const;
+   virtual void CalculateDValues(const Vector &c,
+                                 const Vector (&dc)[3],
+                                 const Vector &baseShape,
+                                 const DenseMatrix &baseDShape,
+                                 const IntegrationPoint &ip,
+                                 DenseMatrix &dshape) const;
+
+   // The corrections make the distance dimensionless, so no smoothing parameter
+   virtual void DistanceVec(const int i,
+                            const Vector &x,
+                            Vector &y) const;
+
+public:
+   /** @brief Construct RBFFiniteElement
+       @param D           Reference space dimension
+       @param numPointsD  Number of points across the element in each dimension
+       @param rbfType     Type of radial basis function, from RBFType
+       @param distNorm    Norm used for distance, usually 2 = Euclidean distance
+       @param order       Order of polynomial correction, >= 0
+       @param intOrder    Number of integration points per RBF point in each dimension
+       @param h           Shape parameter, approximately equal to the number of points in the support radius in one dimension
+       @param faceFactor  1.0 = points end on face, 0.0 = points end at dx/2 from face
+   */
+   RKFiniteElement(const int D,
+                   const int numPointsD,
+                   const int rbfType,
+                   const int distNorm,
+                   const int order,
+                   const int intOrder,
+                   const double h,
+                   const double faceFactor);
+   virtual ~RKFiniteElement() { delete baseFE; }
+
+   virtual bool IsCompact() const { return baseFE->IsCompact(); }
+   virtual const RBFKernel *Kernel() const { return baseFE->Kernel(); }
+
+   static int GetNumPoly(int polyOrd, int dim);
+
+   virtual void CalcShape(const IntegrationPoint &ip,
+                          Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+};
+
+} // end namespace mfem
+
+#endif

fem/fe_coll.cpp

@@ -347,6 +347,36 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
          fec = new NURBSFECollection();
       }
    }
+   else if (!strncmp(name, "RBF", 3) || !strncmp(name, "RK", 2))
+   {
+      // Example: RK4_G_2_V_2D_0020_1.00_4.01
+      // (RK order 4, Gaussian, L2 dist, Value map, 2 dimensions,
+      // 20 points across element, face factor 1.0, smoothing length of 4.01)
+      const int dim = atoi(name + 10);
+      const int numPoints = atoi(name + 13);
+      const double h = atof(name + 23);
+      const double faceFactor = atof(name + 18);
+      const int rbfType = RBFType::GetType(name[4]);
+      const int distNorm = atoi(name + 6);
+      const int mapType = (name[8] == 'V'
+                           ? FiniteElement::VALUE
+                           : FiniteElement::INTEGRAL);
+      const int intOrder = 2; // keep this fixed for now
+
+      if (!strncmp(name, "RK", 2))
+      {
+         int order = atoi(name + 2);
+         fec = new LocalKernelFECollection(dim, numPoints, rbfType,
+                                           order, h, faceFactor,
+                                           intOrder, distNorm, mapType);
+      }
+      else
+      {
+         fec = new LocalKernelFECollection(dim, numPoints, rbfType,
+                                           -1, h, faceFactor,
+                                           intOrder, distNorm, mapType);
+      }
+   }
    else
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
@@ -3510,4 +3540,184 @@ FiniteElementCollection *NURBSFECollection::GetTraceCollection() const
    return NULL;
 }
 
+LocalKernelFECollection::LocalKernelFECollection(const int dim,
+                                                 const int numPointsD,
+                                                 const int rbfType,
+                                                 const int order,
+                                                 const double h,
+                                                 const double faceFactor,
+                                                 const int intOrder,
+                                                 const int distNorm,
+                                                 const int mapType)
+{
+   const char *mapStr = NULL;
+   switch (mapType)
+   {
+      case FiniteElement::VALUE:    mapStr = "V"; break;
+      case FiniteElement::INTEGRAL: mapStr = "I"; break;
+      default:
+         MFEM_ABORT("invalid mapType: " << mapType);
+   }
+   if (order == -1)
+   {
+      snprintf(d_name, 32, "RBF_%c_%d_%s_%dD_%04d_%1.2f_%.2f",
+               (int)RBFType::GetChar(rbfType), distNorm,
+               mapStr, dim, numPointsD, faceFactor, h);
+   }
+   else if (order >= 0)
+   {
+      snprintf(d_name, 32, "RK%d_%c_%d_%s_%dD_%04d_%1.2f_%.2f", order,
+               (int)RBFType::GetChar(rbfType), distNorm,
+               mapStr, dim, numPointsD, faceFactor, h);
+   }
+   else
+   {
+      MFEM_ABORT("invalid order: " << order);
+   }
+
+   for (int g = 0; g < Geometry::NumGeom; ++g)
+   {
+      L2_Elements[g] = NULL;
+      Tr_Elements[g] = NULL;
+   }
+   for (int i = 0; i < 2; i++)
+   {
+      SegDofOrd[i] = NULL;
+   }
+   OtherDofOrd = NULL;
+
+   if (dim == 0)
+   {
+      L2_Elements[Geometry::POINT] = new PointFiniteElement;
+   }
+   else if (dim == 1)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::SEGMENT]
+            = new RBFFiniteElement(1, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+      }
+      else
+      {
+         L2_Elements[Geometry::SEGMENT]
+            = new RKFiniteElement(1, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+      }
+      L2_Elements[Geometry::SEGMENT]->SetMapType(mapType);
+      Tr_Elements[Geometry::POINT] = new PointFiniteElement;
+   }
+   else if (dim == 2)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::SQUARE]
+            = new RBFFiniteElement(2, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SEGMENT]
+            = new RBFFiniteElement(1, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+      }
+      else
+      {
+         L2_Elements[Geometry::SQUARE]
+            = new RKFiniteElement(2, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SEGMENT]
+            = new RKFiniteElement(1, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+      }
+      L2_Elements[Geometry::SQUARE]->SetMapType(mapType);
+   }
+   else if (dim == 3)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::CUBE]
+            = new RBFFiniteElement(3, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SQUARE]
+            = new RBFFiniteElement(2, numPointsD, rbfType, distNorm,
+                                   intOrder, h, faceFactor);
+      }
+      else
+      {
+         L2_Elements[Geometry::CUBE]
+            = new RKFiniteElement(3, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+         Tr_Elements[Geometry::SQUARE]
+            = new RKFiniteElement(2, numPointsD, rbfType, distNorm,
+                                  order, intOrder, h, faceFactor);
+      }
+      L2_Elements[Geometry::CUBE]->SetMapType(mapType);
+   }
+
+   if (dim == 1)
+   {
+      SegDofOrd[0] = new int[2*numPointsD];
+      SegDofOrd[1] = SegDofOrd[0] + numPointsD;
+      for (int i = 0; i < numPointsD; ++i)
+      {
+         SegDofOrd[0][i] = i;
+         SegDofOrd[1][i] = numPointsD - i - 1;
+      }
+   }
+   else
+   {
+      const int geomType = TensorBasisElement::GetTensorProductGeometry(dim);
+      const int dof = L2_Elements[geomType]->GetDof();
+      OtherDofOrd = new int[dof];
+      for (int i = 0; i < dof; ++i)
+      {
+         OtherDofOrd[i] = i;
+      }
+   }
 }
+
+LocalKernelFECollection::~LocalKernelFECollection()
+{
+   delete [] OtherDofOrd;
+   delete [] SegDofOrd[0];
+   for (int i = 0; i < Geometry::NumGeom; ++i)
+   {
+      delete L2_Elements[i];
+   }
+}
+
+const FiniteElement *
+LocalKernelFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   return L2_Elements[GeomType];
+}
+
+const FiniteElement *
+LocalKernelFECollection::TraceFiniteElementForGeometry(Geometry::Type GeomType)
+const
+{
+   return Tr_Elements[GeomType];
+}
+
+int LocalKernelFECollection::DofForGeometry(Geometry::Type GeomType) const
+{
+   if (L2_Elements[GeomType])
+   {
+      return L2_Elements[GeomType]->GetDof();
+   }
+   return 0;
+}
+
+const int *LocalKernelFECollection::DofOrderForOrientation(
+   Geometry::Type GeomType,
+   int Or) const
+{
+   if (GeomType == Geometry::SEGMENT)
+   {
+      return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
+   }
+   else
+   {
+      return (Or == 0) ? OtherDofOrd : NULL;
+   }
+}
+
+} // namespace mfem

fem/fe_coll.hpp

@@ -1305,6 +1305,64 @@ public:
    virtual int GetContType() const { return DISCONTINUOUS; }
 };
 
-}
+/** @brief Radial basis function collection.
+    Instead of using polynomials as a basis, these methods use
+    Gaussian-like functions centered at chosen points. Reproducing
+    kernels augment the radial basis functions so that the resultant
+    basis can exactly represent functions of a chosen polynomial degree. */
+class LocalKernelFECollection : public FiniteElementCollection
+{
+private:
+   int maxDim;
+   char d_name[32];
+   ScalarFiniteElement *Tr_Elements[Geometry::NumGeom];
+   ScalarFiniteElement *L2_Elements[Geometry::NumGeom];
+   int *SegDofOrd[2]; // for rotating segment dofs in 1D
+   int *OtherDofOrd;
+
+public:
+   /** @brief Create a radial basis function (RBF) or reproducing kernel (RK)
+       finite element collection.
+       @param D           Reference space dimension
+       @param numPointsD  Number of points across the element in each dimension
+       @param rbfType     Type of radial basis function, from RBFType
+       @param order       Order of polynomial correction, >= 0 for RK or -1 for
+                          RBF
+       @param h           Shape parameter, approximately equal to the number of
+                          points in the support radius in one dimension, should
+                          generally be > order
+       @param distNorm    Norm used for distance, usually 2 = Euclidean distance
+       @param intOrder    Number of integration points per RBF point in each
+                          dimension
+       @param faceFactor  1.0 = points end on face, 0.0 = points end at dx/2
+                          from face
+   */
+   LocalKernelFECollection(const int D,
+                           const int numPointsD,
+                           const int rbfType,
+                           const int order,
+                           const double h,
+                           const double faceFactor = 0.0,
+                           const int intOrder = 2,
+                           const int distNorm = 2,
+                           const int mapType = FiniteElement::VALUE);
+   virtual ~LocalKernelFECollection();
+
+   virtual const FiniteElement *
+   FiniteElementForGeometry(Geometry::Type GeomType) const;
+
+   virtual const FiniteElement *
+   TraceFiniteElementForGeometry(Geometry::Type GeomType) const;
+
+   virtual int DofForGeometry(Geometry::Type GeomType) const;
+
+   virtual const int * DofOrderForOrientation(Geometry::Type GeomType,
+                                              int Or) const;
+
+   virtual const char * Name() const { return d_name; }
+   virtual int GetContType() const { return DISCONTINUOUS; }
+};
+
+} // namespace mfem
 
 #endif

general/mem_manager.cpp

@@ -851,7 +851,8 @@ void *MemoryManager::Register_(void *ptr, void *h_tmp, size_t bytes,
 
 void MemoryManager::Register2_(void *h_ptr, void *d_ptr, size_t bytes,
                                MemoryType h_mt, MemoryType d_mt,
-                               bool own, bool alias, unsigned &flags)
+                               bool own, bool alias, unsigned &flags,
+                               unsigned valid_flags)
 {
    MFEM_CONTRACT_VAR(alias);
    MFEM_ASSERT(exists, "Internal error!");
@@ -871,7 +872,7 @@ void MemoryManager::Register2_(void *h_ptr, void *d_ptr, size_t bytes,
    mm.InsertDevice(d_ptr, h_ptr, bytes, h_mt, d_mt);
    flags = (own ? flags | (Mem::OWNS_HOST | Mem::OWNS_DEVICE) :
             flags & ~(Mem::OWNS_HOST | Mem::OWNS_DEVICE)) |
-           Mem::VALID_HOST;
+           valid_flags;
 
    CheckHostMemoryType_(h_mt, h_ptr, alias);
 }

general/mem_manager.hpp

@@ -367,7 +367,8 @@ public:
 
    /** Wrap an externally pair of allocated pointers, @a h_ptr and @a d_ptr,
        of the given host MemoryType @a h_mt. */
-   /** The new memory object will have the device MemoryType set as valid.
+   /** The new memory object will have the device MemoryType set as valid unless
+       specified otherwise by the parameters @a valid_host and @a valid_device.
 
        The given @a h_ptr and @a d_ptr must be allocated appropriately for the
        given host MemoryType and its dual device MemoryType as defined by
@@ -376,13 +377,18 @@ public:
        The parameter @a own determines whether both @a h_ptr and @a d_ptr will
        be deleted when the method Delete() is called.
 
+       The parameters @a valid_host and @a valid_device determine which
+       pointers, host and/or device, will be marked as valid; at least one of
+       the two parameters must be set to true.
+
        @note Ownership can also be controlled by using the following methods:
          - ClearOwnerFlags,
          - SetHostPtrOwner,
          - SetDevicePtrOwner.
 
        @note The current memory is NOT deleted by this method. */
-   inline void Wrap(T *h_ptr, T *d_ptr, int size, MemoryType h_mt, bool own);
+   inline void Wrap(T *h_ptr, T *d_ptr, int size, MemoryType h_mt, bool own,
+                    bool valid_host = false, bool valid_device = true);
 
    /// Create a memory object that points inside the memory object @a base.
    /** The new Memory object uses the same MemoryType(s) as @a base.
@@ -645,7 +651,8 @@ private: // Static methods used by the Memory<T> class
    /// Register a pair of external host and device pointers
    static void Register2_(void *h_ptr, void *d_ptr, size_t bytes,
                           MemoryType h_mt, MemoryType d_mt,
-                          bool own, bool alias, unsigned &flags);
+                          bool own, bool alias, unsigned &flags,
+                          unsigned valid_flags);
 
    /// Register an alias. Note: base_h_ptr may be an alias.
    static void Alias_(void *base_h_ptr, size_t offset, size_t bytes,
@@ -958,17 +965,20 @@ inline void Memory<T>::Wrap(T *ptr, int size, MemoryType mt, bool own)
 }
 
 template <typename T>
-inline void Memory<T>::Wrap(T *ptr, T *d_ptr, int size, MemoryType mt, bool own)
+inline void Memory<T>::Wrap(T *h_ptr_, T *d_ptr, int size, MemoryType h_mt_,
+                            bool own, bool valid_host, bool valid_device)
 {
-   h_mt = mt;
+   h_mt = h_mt_;
    flags = 0;
-   h_ptr = ptr;
+   h_ptr = h_ptr_;
    capacity = size;
    MFEM_ASSERT(IsHostMemory(h_mt),"");
+   MFEM_ASSERT(valid_host || valid_device,"");
    const size_t bytes = size*sizeof(T);
    const MemoryType d_mt = MemoryManager::GetDualMemoryType(h_mt);
    MemoryManager::Register2_(h_ptr, d_ptr, bytes, h_mt, d_mt,
-                             own, false, flags);
+                             own, false, flags,
+                             valid_host*VALID_HOST|valid_device*VALID_DEVICE);
 }
 
 template <typename T>

linalg/sundials.cpp

@@ -20,8 +20,10 @@
 
 // SUNDIALS vectors
 #include <nvector/nvector_serial.h>
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA)
 #include <nvector/nvector_cuda.h>
+#elif defined(MFEM_USE_HIP)
+#include <nvector/nvector_hip.h>
 #endif
 #ifdef MFEM_USE_MPI
 #include <nvector/nvector_mpiplusx.h>
@@ -35,6 +37,14 @@
 // Access SUNDIALS object's content pointer
 #define GET_CONTENT(X) ( X->content )
 
+#if defined(MFEM_USE_CUDA)
+#define SUN_Hip_OR_Cuda(X) X##_Cuda
+#define SUN_HIP_OR_CUDA(X) X##_CUDA
+#elif defined(MFEM_USE_HIP)
+#define SUN_Hip_OR_Cuda(X) X##_Hip
+#define SUN_HIP_OR_CUDA(X) X##_HIP
+#endif
+
 using namespace std;
 
 #if (SUNDIALS_VERSION_MAJOR < 6)
@@ -112,16 +122,16 @@ MFEM_DEPRECATED N_Vector N_VNewEmpty_Parallel(MPI_Comm comm,
 
 #endif // MFEM_USE_MPI
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
 
 /// (DEPRECATED) Wrapper function for backwards compatibility with SUNDIALS
 /// version < 6
-MFEM_DEPRECATED N_Vector N_VNewWithMemHelp_Cuda(sunindextype length,
-                                                booleantype use_managed_mem,
-                                                SUNMemoryHelper helper,
-                                                SUNContext)
+MFEM_DEPRECATED N_Vector SUN_Hip_OR_Cuda(N_VNewWithMemHelp)(sunindextype length,
+                                                            booleantype use_managed_mem,
+                                                            SUNMemoryHelper helper,
+                                                            SUNContext)
 {
-   return N_VNewWithMemHelp_Cuda(length, use_managed_mem, helper);
+   return SUN_Hip_OR_Cuda(N_VNewWithMemHelp)(length, use_managed_mem, helper);
 }
 
 /// (DEPRECATED) Wrapper function for backwards compatibility with SUNDIALS
@@ -131,9 +141,9 @@ MFEM_DEPRECATED SUNMemoryHelper SUNMemoryHelper_NewEmpty(SUNContext)
    return SUNMemoryHelper_NewEmpty();
 }
 
-#endif // MFEM_USE_CUDA
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
 
-#if defined(MFEM_USE_MPI) && defined(MFEM_USE_CUDA)
+#if defined(MFEM_USE_MPI) && (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
 
 /// (DEPRECATED) Wrapper function for backwards compatibility with SUNDIALS
 /// version < 6
@@ -143,7 +153,7 @@ MFEM_DEPRECATED N_Vector N_VMake_MPIPlusX(MPI_Comm comm, N_Vector local_vector,
    return N_VMake_MPIPlusX(comm, local_vector);
 }
 
-#endif // MFEM_USE_MPI && MFEM_USE_CUDA
+#endif // MFEM_USE_MPI && (MFEM_USE_CUDA || MFEM_USE_HIP)
 
 #endif // SUNDIALS_VERSION_MAJOR < 6
 
@@ -206,7 +216,7 @@ Sundials::~Sundials()
 
 #endif // SUNDIALS_VERSION_MAJOR >= 6
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
 SundialsMemHelper::SundialsMemHelper(SUNContext context)
 {
    /* Allocate helper */
@@ -215,8 +225,8 @@ SundialsMemHelper::SundialsMemHelper(SUNContext context)
    /* Set the ops */
    h->ops->alloc     = SundialsMemHelper_Alloc;
    h->ops->dealloc   = SundialsMemHelper_Dealloc;
-   h->ops->copy      = SUNMemoryHelper_Copy_Cuda;
-   h->ops->copyasync = SUNMemoryHelper_CopyAsync_Cuda;
+   h->ops->copy      = SUN_Hip_OR_Cuda(SUNMemoryHelper_Copy);
+   h->ops->copyasync = SUN_Hip_OR_Cuda(SUNMemoryHelper_CopyAsync);
 }
 
 SundialsMemHelper::SundialsMemHelper(SundialsMemHelper&& that_helper)
@@ -240,25 +250,25 @@ int SundialsMemHelper::SundialsMemHelper_Alloc(SUNMemoryHelper helper,
 #endif
                                               )
 {
-   int length = memsize/sizeof(double);
    SUNMemory sunmem = SUNMemoryNewEmpty();
 
    sunmem->ptr = NULL;
    sunmem->own = SUNTRUE;
 
+   // memsize is the number of bytes to allocate, so we use Memory<char>
    if (mem_type == SUNMEMTYPE_HOST)
    {
-      Memory<double> mem(length, Device::GetHostMemoryType());
+      Memory<char> mem(memsize, Device::GetHostMemoryType());
       mem.SetHostPtrOwner(false);
-      sunmem->ptr  = mfem::HostReadWrite(mem, length);
+      sunmem->ptr  = mfem::HostReadWrite(mem, memsize);
       sunmem->type = SUNMEMTYPE_HOST;
       mem.Delete();
    }
    else if (mem_type == SUNMEMTYPE_DEVICE || mem_type == SUNMEMTYPE_UVM)
    {
-      Memory<double> mem(length, Device::GetDeviceMemoryType());
+      Memory<char> mem(memsize, Device::GetDeviceMemoryType());
       mem.SetDevicePtrOwner(false);
-      sunmem->ptr  = mfem::ReadWrite(mem, length);
+      sunmem->ptr  = mfem::ReadWrite(mem, memsize);
       sunmem->type = mem_type;
       mem.Delete();
    }
@@ -283,14 +293,14 @@ int SundialsMemHelper::SundialsMemHelper_Dealloc(SUNMemoryHelper helper,
    {
       if (sunmem->type == SUNMEMTYPE_HOST)
       {
-         Memory<double> mem(static_cast<double*>(sunmem->ptr), 1,
-                            Device::GetHostMemoryType(), true);
+         Memory<char> mem(static_cast<char*>(sunmem->ptr), 1,
+                          Device::GetHostMemoryType(), true);
          mem.Delete();
       }
       else if (sunmem->type == SUNMEMTYPE_DEVICE || sunmem->type == SUNMEMTYPE_UVM)
       {
-         Memory<double> mem(static_cast<double*>(sunmem->ptr), 1,
-                            Device::GetDeviceMemoryType(), true);
+         Memory<char> mem(static_cast<char*>(sunmem->ptr), 1,
+                          Device::GetDeviceMemoryType(), true);
          mem.Delete();
       }
       else
@@ -303,7 +313,7 @@ int SundialsMemHelper::SundialsMemHelper_Dealloc(SUNMemoryHelper helper,
    return 0;
 }
 
-#endif // MFEM_USE_CUDA
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
 
 
 // ---------------------------------------------------------------------------
@@ -329,12 +339,13 @@ void SundialsNVector::_SetNvecDataAndSize_(long glob_size)
          NV_LENGTH_S(local_x) = size;
          break;
       }
-#ifdef MFEM_USE_CUDA
-      case SUNDIALS_NVEC_CUDA:
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+      case SUN_HIP_OR_CUDA(SUNDIALS_NVEC):
       {
-         N_VSetHostArrayPointer_Cuda(HostReadWrite(), local_x);
-         N_VSetDeviceArrayPointer_Cuda(ReadWrite(), local_x);
-         static_cast<N_VectorContent_Cuda>(GET_CONTENT(local_x))->length = size;
+         SUN_Hip_OR_Cuda(N_VSetHostArrayPointer)(HostReadWrite(), local_x);
+         SUN_Hip_OR_Cuda(N_VSetDeviceArrayPointer)(ReadWrite(), local_x);
+         static_cast<SUN_Hip_OR_Cuda(N_VectorContent)>(GET_CONTENT(
+                                                          local_x))->length = size;
          break;
       }
 #endif
@@ -403,14 +414,14 @@ void SundialsNVector::_SetDataAndSize_()
          if (known) { data.ClearOwnerFlags(); }
          break;
       }
-#ifdef MFEM_USE_CUDA
-      case SUNDIALS_NVEC_CUDA:
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+      case SUN_HIP_OR_CUDA(SUNDIALS_NVEC):
       {
-         double *h_ptr = N_VGetHostArrayPointer_Cuda(local_x);
-         double *d_ptr = N_VGetDeviceArrayPointer_Cuda(local_x);
+         double *h_ptr = SUN_Hip_OR_Cuda(N_VGetHostArrayPointer)(local_x);
+         double *d_ptr = SUN_Hip_OR_Cuda(N_VGetDeviceArrayPointer)(local_x);
          const bool known = mm.IsKnown(h_ptr);
-         size = N_VGetLength_Cuda(local_x);
-         data.Wrap(h_ptr, d_ptr, size, Device::GetHostMemoryType(), false);
+         size = SUN_Hip_OR_Cuda(N_VGetLength)(local_x);
+         data.Wrap(h_ptr, d_ptr, size, Device::GetHostMemoryType(), false, false, true);
          if (known) { data.ClearOwnerFlags(); }
          UseDevice(true);
          break;
@@ -525,11 +536,12 @@ void SundialsNVector::SetDataAndSize(double *d, int s, long glob_size)
 N_Vector SundialsNVector::MakeNVector(bool use_device)
 {
    N_Vector x;
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_device)
    {
-      x = N_VNewWithMemHelp_Cuda(0, UseManagedMemory(), Sundials::GetMemHelper(),
-                                 Sundials::GetContext());
+      x = SUN_Hip_OR_Cuda(N_VNewWithMemHelp)(0, UseManagedMemory(),
+                                             Sundials::GetMemHelper(),
+                                             Sundials::GetContext());
    }
    else
    {
@@ -555,12 +567,13 @@ N_Vector SundialsNVector::MakeNVector(MPI_Comm comm, bool use_device)
    }
    else
    {
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
       if (use_device)
       {
-         x = N_VMake_MPIPlusX(comm, N_VNewWithMemHelp_Cuda(0, UseManagedMemory(),
-                                                           Sundials::GetMemHelper(),
-                                                           Sundials::GetContext()),
+         x = N_VMake_MPIPlusX(comm, SUN_Hip_OR_Cuda(N_VNewWithMemHelp)(0,
+                                                                       UseManagedMemory(),
+                                                                       Sundials::GetMemHelper(),
+                                                                       Sundials::GetContext()),
                               Sundials::GetContext());
       }
       else
@@ -569,7 +582,7 @@ N_Vector SundialsNVector::MakeNVector(MPI_Comm comm, bool use_device)
       }
 #else
       x = N_VNewEmpty_Parallel(comm, 0, 0, Sundials::GetContext());
-#endif // MFEM_USE_CUDA
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
    }
 
    MFEM_VERIFY(x, "Error in SundialsNVector::MakeNVector.");

linalg/sundials.hpp

@@ -32,13 +32,24 @@
 #if defined(MFEM_USE_CUDA) && ((SUNDIALS_VERSION_MAJOR == 5) && (SUNDIALS_VERSION_MINOR < 4))
 #error MFEM requires SUNDIALS version 5.4.0 or newer when MFEM_USE_CUDA=TRUE!
 #endif
+#if defined(MFEM_USE_HIP) && ((SUNDIALS_VERSION_MAJOR == 5) && (SUNDIALS_VERSION_MINOR < 7))
+#error MFEM requires SUNDIALS version 5.7.0 or newer when MFEM_USE_HIP=TRUE!
+#endif
+#if defined(MFEM_USE_CUDA) && !defined(SUNDIALS_NVECTOR_CUDA)
+#error MFEM_USE_CUDA=TRUE requires SUNDIALS to be built with CUDA support
+#endif
+#if defined(MFEM_USE_HIP) && !defined(SUNDIALS_NVECTOR_HIP)
+#error MFEM_USE_HIP=TRUE requires SUNDIALS to be built with HIP support
+#endif
 #include <sundials/sundials_matrix.h>
 #include <sundials/sundials_linearsolver.h>
 #include <arkode/arkode_arkstep.h>
 #include <cvodes/cvodes.h>
 #include <kinsol/kinsol.h>
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA)
 #include <sunmemory/sunmemory_cuda.h>
+#elif defined(MFEM_USE_HIP)
+#include <sunmemory/sunmemory_hip.h>
 #endif
 
 #include <functional>
@@ -62,10 +73,10 @@ using SUNContext = void*;
 namespace mfem
 {
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
 
 // ---------------------------------------------------------------------------
-// SUNMemory interface class (used when CUDA is enabled)
+// SUNMemory interface class (used when CUDA or HIP is enabled)
 // ---------------------------------------------------------------------------
 class SundialsMemHelper
 {
@@ -113,10 +124,10 @@ public:
 
 };
 
-#else // MFEM_USE_CUDA
+#else // MFEM_USE_CUDA || MFEM_USE_HIP
 
 // ---------------------------------------------------------------------------
-// Dummy SUNMemory interface class (used when CUDA is not enabled)
+// Dummy SUNMemory interface class (used when CUDA or HIP is not enabled)
 // ---------------------------------------------------------------------------
 class SundialsMemHelper
 {
@@ -130,7 +141,7 @@ public:
    }
 };
 
-#endif // MFEM_USE_CUDA
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
 
 
 /// Singleton class for SUNContext and SundialsMemHelper objects
@@ -290,17 +301,17 @@ public:
 #endif
 
    /// Create a N_Vector.
-   /** @param[in] use_device  If true, use the SUNDIALS CUDA N_Vector. */
+   /** @param[in] use_device  If true, use the SUNDIALS CUDA or HIP N_Vector. */
    static N_Vector MakeNVector(bool use_device);
 
 #ifdef MFEM_USE_MPI
    /// Create a parallel N_Vector.
    /** @param[in] comm  The MPI communicator to use.
-       @param[in] use_device  If true, use the SUNDIALS CUDA N_Vector. */
+       @param[in] use_device  If true, use the SUNDIALS CUDA or HIP N_Vector. */
    static N_Vector MakeNVector(MPI_Comm comm, bool use_device);
 #endif
 
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    static bool UseManagedMemory()
    {
       return Device::GetDeviceMemoryType() == MemoryType::MANAGED;

linalg/vector.cpp

@@ -15,13 +15,6 @@
 #include "vector.hpp"
 #include "../general/forall.hpp"
 
-#if defined(MFEM_USE_SUNDIALS)
-#include "sundials.hpp"
-#if defined(MFEM_USE_MPI)
-#include <nvector/nvector_parallel.h>
-#endif
-#endif
-
 #ifdef MFEM_USE_OPENMP
 #include <omp.h>
 #endif

linalg/vector.hpp

@@ -19,9 +19,7 @@
 #include "../general/globals.hpp"
 #include "../general/mem_manager.hpp"
 #include "../general/device.hpp"
-#ifdef MFEM_USE_SUNDIALS
-#include <nvector/nvector_serial.h>
-#endif
+
 #include <cmath>
 #include <iostream>
 #include <limits>

mesh/mesh.cpp

@@ -7943,7 +7943,7 @@ void Mesh::GetNode(int i, double *coord) const
       FiniteElementSpace *fes = Nodes->FESpace();
       for (int j = 0; j < spaceDim; j++)
       {
-         coord[j] = (*Nodes)(fes->DofToVDof(i, j));
+         coord[j] = AsConst(*Nodes)(fes->DofToVDof(i, j));
       }
    }
    else