Update examples/sundials/ex16.cpp

examples/sundials/ex16.cpp

@@ -28,10 +28,8 @@
 //
 //               The example demonstrates the use of nonlinear operators (the
 //               class ConductionOperator defining C(u)), as well as their
-//               implicit time integration. Note that implementing the method
-//               ConductionOperator::ImplicitSolve is the only requirement for
-//               high-order implicit (SDIRK) time integration. By default, this
-//               example uses the SUNDIALS ODE solvers from CVODE and ARKODE.
+//               implicit time integration. By default, this example uses the
+//               SUNDIALS ODE solvers from CVODE and ARKODE.
 //
 //               We recommend viewing examples 2, 9 and 10 before viewing this
 //               example.
@@ -51,15 +49,16 @@ using namespace mfem;
  *  and K(u) is the diffusion operator with diffusivity depending on u:
  *  (\kappa + \alpha u).
  *
- *  Class ConductionOperatorOperator represents the above ODE operator in the
- *  general form F(u, k, t) = G(u, t) where
+ *  Class ConductionOperator represents the above ODE operator as a
+ *  TimeDependentOperator for use with native MFEM integrators and CVODE
+ *  integrators, i.e., F(u, k, t) = G(u, t) with F(u, du/dt, t) = du/dt and
+ *  G(u, t) = -K(u) u
  *
- *    1. F(u, du/dt, t) = du/dt            (ODE is expressed in EXPLICIT form)
- *       G(u, t) = - inv(M) K(u) u
- *    2. F(u, du/dt, t) = M du/dt          (ODE is expressed in IMPLICIT form)
- *       G(u, t) = - K(u) u
+ *  Class ConductionOperator represents the above ODE operator as an
+ *  ARKStepODE for use with ARKODE integrators, i.e., either M du/dt = -K(u) u
+ *  (mass form) or du/dt = -inv(M) K(u) u (MFEM form)
  */
-class ConductionOperator : public TimeDependentOperator
+class ConductionOperator : public TimeDependentOperator, public ARKStepODE
 {
    FiniteElementSpace &fespace;
    Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
@@ -81,50 +80,90 @@ class ConductionOperator : public TimeDependentOperator
 
    mutable Vector z; // auxiliary vector
 
+   const bool use_mass_form;
+
 public:
 
    ConductionOperator(FiniteElementSpace &f, const real_t alpha,
                       const real_t kappa, const Vector &u,
-                      const Type &ode_expression_type);
+                      const bool use_mass_form);
 
    // Compute K(u_n) for use as an approximation in - K(u) u
    void SetConductionTensor(const Vector &u);
 
-   /** Compute G(u, t) as defined in the IMPLICIT expression form of the ODE
-    operator, i.e., @a v = - K(u_n) @a u. Note that K(u_n) is an
-      approximation to K(u). */
-   void ExplicitMult(const Vector &u, Vector &v) const override;
+   // ********* methods for MFEM native time integrators *********
 
-   /** Solve for k in F(u, k, t) = G(u, t) for either EXPLICIT or IMPLICIT
-       expression forms of the ODE operator, i.e., @a k = - inv(M) K(u_n) @a u.
+   /** Solve for k in F(u, k, t) = G(u, t), i.e., @a k = - inv(M) K(u_n) @a u.
        Note that K(u_n) is an approximation to K(u). */
    void Mult(const Vector &u, Vector &k) const override;
 
-   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t) for either EXPLICIT
-       or IMPLICIT expression forms of the ODE operator, i.e.,
-       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u . Note that K(u_n) is an
-       approximation to K(u). */
+   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t), i.e.,
+       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u .
+       Note that K(u_n) is an approximation to K(u). */
    void ImplicitSolve(const real_t gam, const Vector &u, Vector &k) override;
 
-   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F for
-       either EXPLICIT or IMPLICIT expression forms of the ODE operator, i.e.,
-       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
-       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
-   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
-                        real_t gam) override;
+      // ********* methods for ARKODE time integrators *********
+
+   // TODO: add comments
+   int ARKSize() const override;
+
+   // TODO: add comments
+   bool ARKInMassForm() const override;
+
+   // TODO: add comments
+   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const override;
+
+   // TODO: add comments
+   int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &fu,
+                        int jok, int *jcur, real_t gam) override;
 
    /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
        with the residual @a r providing either
-        1. @a r = G - F = inv(M) f(u) - k       (EXPLICIT expression form)
-        1. @a r = G - F = f(u) - M k            (IMPLICIT expression form)
+        1. @a r = G - F = inv(M) f(u) - k       (MFEM form)
+        1. @a r = G - F = f(u) - M k            (mass form)
        */
-   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
+   int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
 
-   int SUNMassSetup() override;
+   int ARKMassSetup(const real_t t) override;
 
-   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override;
+   int ARKMassSolve(const Vector &b, Vector &x, real_t tol) override;
 
-   int SUNMassMult(const Vector &x, Vector &v) override;
+   int ARKMassMult(const Vector &x, Vector &v) override;
+
+   // ********* methods for CVODE time integrators *********
+   // note these methods merely call the corresponding ARKStepODE methods until
+   // the CVODESolver is refactored to use specialized interface like ARKStepODE
+
+   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F, i.e.,
+       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
+       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
+   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
+                        real_t gam) override
+   {
+      return ARKImplicitSetup(u, 0.0, fu, jok, jcur, gam); // the ODE is autonomous
+   }
+
+   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
+       with the residual @a r providing @a r = G - F = inv(M) f(u) - k. */
+   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override
+   {
+      return ARKImplicitSolve(r, dk, tol);
+   }
+
+   int SUNMassSetup() override
+   {
+      return ARKMassSetup(0.0); // the ODE is autonomous
+   }
+
+   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override
+   {
+      return ARKMassSolve(b, x, tol);
+   }
+
+   int SUNMassMult(const Vector &x, Vector &v) override
+   {
+      return ARKMassMult(x, v);
+   }
 };
 
 real_t InitialTemperature(const Vector &x)
@@ -245,16 +284,7 @@ int main(int argc, char *argv[])
    u_gf.GetTrueDofs(u);
 
    // 6. Initialize the conduction ODE operator and the visualization.
-   ConductionOperator::Type ode_expression_type;
-   if (use_mass_solver)
-   {
-      ode_expression_type = ConductionOperator::Type::IMPLICIT;
-   }
-   else
-   {
-      ode_expression_type = ConductionOperator::Type::EXPLICIT;
-   }
-   ConductionOperator oper(fespace, alpha, kappa, u, ode_expression_type);
+   ConductionOperator oper(fespace, alpha, kappa, u, use_mass_solver);
 
    u_gf.SetFromTrueDofs(u);
    {
@@ -352,7 +382,7 @@ int main(int argc, char *argv[])
          }
          std::unique_ptr<ARKStepSolver> arkode(
             new ARKStepSolver(arkode_solver_type));
-         arkode->Init(oper);
+         arkode->Init(&oper);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
          if (ode_solver_type == 11 || ode_solver_type == 14)
@@ -445,9 +475,10 @@ int main(int argc, char *argv[])
 ConductionOperator::ConductionOperator(FiniteElementSpace &fes,
                                        const real_t alpha, const real_t kappa,
                                        const Vector &u,
-                                       const Type &ode_expression_type)
-   : TimeDependentOperator(fes.GetTrueVSize(), 0.0, ode_expression_type),
-     fespace(fes), M(&fespace), alpha(alpha), kappa(kappa), z(height)
+                                       const bool use_mass_form)
+   : TimeDependentOperator(fes.GetTrueVSize(), 0.0),
+     fespace(fes), M(&fespace), alpha(alpha), kappa(kappa), z(height),
+     use_mass_form(use_mass_form)
 {
    // specify a relative tolerance for all solves with MFEM integrators
    const real_t rel_tol = 1e-8;
@@ -474,6 +505,16 @@ ConductionOperator::ConductionOperator(FiniteElementSpace &fes,
    SetConductionTensor(u);
 }
 
+int ConductionOperator::ARKSize() const
+{
+   return z.Size();
+}
+
+bool ConductionOperator::ARKInMassForm() const
+{
+   return use_mass_form;
+}
+
 void ConductionOperator::SetConductionTensor(const Vector &u)
 {
    // Compute K(u_n).
@@ -491,17 +532,27 @@ void ConductionOperator::SetConductionTensor(const Vector &u)
    K->FormSystemMatrix(ess_tdof_list, Kmat);
 }
 
-void ConductionOperator::ExplicitMult(const Vector &u, Vector &v) const
+void ConductionOperator::ARKEvaluateRHS(const Vector &u, const real_t t,
+                                        Vector &result) const
 {
-   // Compute - K(u_n) u.
-   Kmat.Mult(u, v);
-   v.Neg();
+   if (use_mass_form) // compute -K(u_n) u.
+   {
+      Kmat.Mult(u, result);
+      result.Neg();
+   }
+   else // compute -inv(M) K(u_n) u
+   {
+      Kmat.Mult(u, z);
+      z.Neg();
+      M_solver.Mult(z, result);
+   }
 }
 
 void ConductionOperator::Mult(const Vector &u, Vector &k) const
 {
    // Compute - inv(M) K(u_n) u.
-   ExplicitMult(u, z);
+   Kmat.Mult(u, z);
+   z.Neg();
    M_solver.Mult(z, k);
 }
 
@@ -509,14 +560,16 @@ void ConductionOperator::ImplicitSolve(const real_t gam, const Vector &u,
                                        Vector &k)
 {
    // Solve for k in M k = - K(u_n) [u + gam*k].
-   ExplicitMult(u, z);
+   Kmat.Mult(u, z);
+   z.Neg();
    T = std::unique_ptr<SparseMatrix>(Add(1.0, Mmat, gam, Kmat));
    T_solver.SetOperator(*T);
    T_solver.Mult(z, k);
 }
 
-int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
-                                         int jok, int *jcur, real_t gam)
+int ConductionOperator::ARKImplicitSetup(const Vector &u, const real_t t,
+                                         const Vector &fu, int jok, int *jcur,
+                                         real_t gam)
 {
    // Compute T = M + gamma K(u_n).
    T = std::unique_ptr<SparseMatrix>(Add(1.0, Mmat, gam, Kmat));
@@ -525,22 +578,22 @@ int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
    return SUN_SUCCESS;
 }
 
-int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
+int ConductionOperator::ARKImplicitSolve(const Vector &r, Vector &dk,
                                          real_t tol)
 {
    // Solve the system [M + gamma K(u_n)] dk = - K(u_n) u - M k.
    // What value r is providing depends on the ODE expression form:
-   //   EXPLICIT form: r = -inv(M) K(u_n) u - k
-   //   IMPLICIT form: r = -K(u_n) u - M k
+   //   MFEM form: r = -inv(M) K(u_n) u - k
+   //   mass form: r = -K(u_n) u - M k
    T_solver.SetRelTol(tol);
-   if (isExplicit())
+   if (use_mass_form)
    {
-      Mmat.Mult(r, z);
-      T_solver.Mult(z, dk);
+      T_solver.Mult(r, dk);
    }
    else
    {
-      T_solver.Mult(r, dk);
+      Mmat.Mult(r, z);
+      T_solver.Mult(z, dk);
    }
    if (T_solver.GetConverged())
    {
@@ -552,13 +605,13 @@ int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
    }
 }
 
-int ConductionOperator::SUNMassSetup()
+int ConductionOperator::ARKMassSetup(const real_t t)
 {
    // Do nothing b/c mass solver was setup in constructor.
    return SUN_SUCCESS;
 }
 
-int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
+int ConductionOperator::ARKMassSolve(const Vector &b, Vector &x, real_t tol)
 {
    // Solve the system M x = b.
    M_solver.SetRelTol(tol);
@@ -573,7 +626,7 @@ int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
    }
 }
 
-int ConductionOperator::SUNMassMult(const Vector &x, Vector &v)
+int ConductionOperator::ARKMassMult(const Vector &x, Vector &v)
 {
    // Compute M x.
    Mmat.Mult(x, v);

examples/sundials/ex16p.cpp

@@ -29,10 +29,8 @@
 //
 //               The example demonstrates the use of nonlinear operators (the
 //               class ConductionOperator defining C(u)), as well as their
-//               implicit time integration. Note that implementing the method
-//               ConductionOperator::ImplicitSolve is the only requirement for
-//               high-order implicit (SDIRK) time integration. By default, this
-//               example uses the SUNDIALS ODE solvers from CVODE and ARKODE.
+//               implicit time integration. By default, this example uses the
+//               SUNDIALS ODE solvers from CVODE and ARKODE.
 //
 //               We recommend viewing examples 2, 9 and 10 before viewing this
 //               example.
@@ -52,15 +50,16 @@ using namespace mfem;
  *  and K(u) is the diffusion operator with diffusivity depending on u:
  *  (\kappa + \alpha u).
  *
- *  Class ConductionOperatorOperator represents the above ODE operator in the
- *  general form F(u, k, t) = G(u, t) where either
+ *  Class ConductionOperator represents the above ODE operator as a
+ *  TimeDependentOperator for use with native MFEM integrators and CVODE
+ *  integrators, i.e., F(u, k, t) = G(u, t) with F(u, du/dt, t) = du/dt and
+ *  G(u, t) = -K(u) u
  *
- *    1. F(u, du/dt, t) = du/dt            (ODE is expressed in EXPLICIT form)
- *       G(u, t) = - inv(M) K(u) u
- *    2. F(u, du/dt, t) = M du/dt          (ODE is expressed in IMPLICIT form)
- *       G(u, t) = - K(u) u
+ *  Class ConductionOperator represents the above ODE operator as an
+ *  ARKStepODE for use with ARKODE integrators, i.e., either M du/dt = -K(u) u
+ *  (mass form) or du/dt = -inv(M) K(u) u (MFEM form)
  */
-class ConductionOperator : public TimeDependentOperator
+class ConductionOperator : public TimeDependentOperator, public ARKStepODE
 {
    ParFiniteElementSpace &fespace;
    Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
@@ -82,50 +81,90 @@ class ConductionOperator : public TimeDependentOperator
 
    mutable Vector z; // auxiliary vector
 
+   const bool use_mass_form;
+
 public:
 
    ConductionOperator(ParFiniteElementSpace &f, const real_t alpha,
                       const real_t kappa, const Vector &u,
-                      const Type &ode_expression_type);
+                      const bool use_mass_form);
 
    // Compute K(u_n) for use as an approximation in - K(u) u
    void SetConductionTensor(const Vector &u);
 
-   /** Compute G(u, t) as defined in the IMPLICIT expression form of the ODE
-       operator, i.e., @a v = - K(u_n) @a u. Note that K(u_n) is an
-       approximation to K(u). */
-   void ExplicitMult(const Vector &u, Vector &v) const override;
+   // ********* methods for MFEM native time integrators *********
 
-   /** Solve for k in F(u, k, t) = G(u, t) for either EXPLICIT or IMPLICIT
-       expression forms of the ODE operator, i.e., @a k = - inv(M) K(u_n) @a u.
+   /** Solve for k in F(u, k, t) = G(u, t), i.e., @a k = - inv(M) K(u_n) @a u.
        Note that K(u_n) is an approximation to K(u). */
    void Mult(const Vector &u, Vector &k) const override;
 
-   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t) for either EXPLICIT
-       or IMPLICIT expression forms of the ODE operator, i.e.,
-       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u . Note that K(u_n) is an
-       approximation to K(u). */
+   /** Solve for k in F(u + gam*k, k, t) = G(u + gam*k, t), i.e.,
+       [ M + @a gam K(u_n) ] @a k = - K(u_n) @a u .
+       Note that K(u_n) is an approximation to K(u). */
    void ImplicitSolve(const real_t gam, const Vector &u, Vector &k) override;
 
-   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F for
-       either EXPLICIT or IMPLICIT expression forms of the ODE operator, i.e.,
-       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
-       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
-   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
-                        real_t gam) override;
+   // ********* methods for ARKODE time integrators *********
+
+   // TODO: add comments
+   int ARKSize() const override;
+
+   // TODO: add comments
+   bool ARKInMassForm() const override;
+
+   // TODO: add comments
+   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const override;
+
+   // TODO: add comments
+   int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &fu,
+                        int jok, int *jcur, real_t gam) override;
 
    /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
        with the residual @a r providing either
-        1. @a r = G - F = inv(M) f(u) - k       (EXPLICIT expression form)
-        1. @a r = G - F = f(u) - M k            (IMPLICIT expression form)
+        1. @a r = G - F = inv(M) f(u) - k       (MFEM form)
+        1. @a r = G - F = f(u) - M k            (mass form)
        */
-   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
+   int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
 
-   int SUNMassSetup() override;
+   int ARKMassSetup(const real_t t) override;
 
-   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override;
+   int ARKMassSolve(const Vector &b, Vector &x, real_t tol) override;
 
-   int SUNMassMult(const Vector &x, Vector &v) override;
+   int ARKMassMult(const Vector &x, Vector &v) override;
+
+   // ********* methods for CVODE time integrators *********
+   // note these methods merely call the corresponding ARKStepODE methods until
+   // the CVODESolver is refactored to use specialized interface like ARKStepODE
+
+   /** Setup to solve for dk in [dF/dk + gam*dF/du - gam*dG/du] dk = G - F, i.e.,
+       [M - @a gam Jf(u)] dk = G - F, where Jf(u) is an approximation of the
+       Jacobian of -K(u) u. The approximation chosen here is Jf(u) = -K(u_n). */
+   int SUNImplicitSetup(const Vector &u, const Vector &fu, int jok, int *jcur,
+                        real_t gam) override
+   {
+      return ARKImplicitSetup(u, 0.0, fu, jok, jcur, gam); // the ODE is autonomous
+   }
+
+   /** Solve for @a dk in the system in SUNImplicitSetup to the given tolerance,
+       with the residual @a r providing @a r = G - F = inv(M) f(u) - k. */
+   int SUNImplicitSolve(const Vector &r, Vector &dk, real_t tol) override
+   {
+      return ARKImplicitSolve(r, dk, tol);
+   }
+
+   int SUNMassSetup() override
+   {
+      return ARKMassSetup(0.0); // the ODE is autonomous
+   }
+
+   int SUNMassSolve(const Vector &b, Vector &x, real_t tol) override
+   {
+      return ARKMassSolve(b, x, tol);
+   }
+
+   int SUNMassMult(const Vector &x, Vector &v) override
+   {
+      return ARKMassMult(x, v);
+   }
 };
 
 real_t InitialTemperature(const Vector &x)
@@ -273,16 +312,7 @@ int main(int argc, char *argv[])
    u_gf.GetTrueDofs(u);
 
    // 8. Initialize the conduction ODE operator and the visualization.
-   ConductionOperator::Type ode_expression_type;
-   if (use_mass_solver)
-   {
-      ode_expression_type = ConductionOperator::Type::IMPLICIT;
-   }
-   else
-   {
-      ode_expression_type = ConductionOperator::Type::EXPLICIT;
-   }
-   ConductionOperator oper(fespace, alpha, kappa, u, ode_expression_type);
+   ConductionOperator oper(fespace, alpha, kappa, u, use_mass_solver);
 
    u_gf.SetFromTrueDofs(u);
    {
@@ -394,7 +424,7 @@ int main(int argc, char *argv[])
          }
          std::unique_ptr<ARKStepSolver> arkode(
             new ARKStepSolver(MPI_COMM_WORLD, arkode_solver_type));
-         arkode->Init(oper);
+         arkode->Init(&oper);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
          if (ode_solver_type == 11 || ode_solver_type == 14)
@@ -497,10 +527,11 @@ int main(int argc, char *argv[])
 ConductionOperator::ConductionOperator(ParFiniteElementSpace &fes,
                                        const real_t alpha, const real_t kappa,
                                        const Vector &u,
-                                       const Type &ode_expression_type)
-   : TimeDependentOperator(fes.GetTrueVSize(), 0.0, ode_expression_type),
+                                       const bool use_mass_form)
+   : TimeDependentOperator(fes.GetTrueVSize(), 0.0),
      fespace(fes), M(&fespace), alpha(alpha), kappa(kappa),
-     M_solver(fes.GetComm()), T_solver(fes.GetComm()), z(height)
+     M_solver(fes.GetComm()), T_solver(fes.GetComm()), z(height),
+     use_mass_form(use_mass_form)
 {
    // specify a relative tolerance for all solves with MFEM integrators
    const real_t rel_tol = 1e-8;
@@ -528,6 +559,16 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &fes,
    SetConductionTensor(u);
 }
 
+int ConductionOperator::ARKSize() const
+{
+   return z.Size();
+}
+
+bool ConductionOperator::ARKInMassForm() const
+{
+   return use_mass_form;
+}
+
 void ConductionOperator::SetConductionTensor(const Vector &u)
 {
    // Compute K(u_n).
@@ -545,17 +586,27 @@ void ConductionOperator::SetConductionTensor(const Vector &u)
    K->FormSystemMatrix(ess_tdof_list, Kmat);
 }
 
-void ConductionOperator::ExplicitMult(const Vector &u, Vector &v) const
+void ConductionOperator::ARKEvaluateRHS(const Vector &u, const real_t t,
+                                        Vector &result) const
 {
-   // Compute - K(u_n) u.
-   Kmat.Mult(u, v);
-   v.Neg();
+   if (use_mass_form) // compute -K(u_n) u.
+   {
+      Kmat.Mult(u, result);
+      result.Neg();
+   }
+   else // compute -inv(M) K(u_n) u
+   {
+      Kmat.Mult(u, z);
+      z.Neg();
+      M_solver.Mult(z, result);
+   }
 }
 
 void ConductionOperator::Mult(const Vector &u, Vector &k) const
 {
    // Compute - inv(M) K(u_n) u.
-   ExplicitMult(u, z);
+   Kmat.Mult(u, z);
+   z.Neg();
    M_solver.Mult(z, k);
 }
 
@@ -563,14 +614,16 @@ void ConductionOperator::ImplicitSolve(const real_t gam, const Vector &u,
                                        Vector &k)
 {
    // Solve for k in M k = - K(u_n) [u + gam*k].
-   ExplicitMult(u, z);
+   Kmat.Mult(u, z);
+   z.Neg();
    T = std::unique_ptr<HypreParMatrix>(Add(1.0, Mmat, gam, Kmat));
    T_solver.SetOperator(*T);
    T_solver.Mult(z, k);
 }
 
-int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
-                                         int jok, int *jcur, real_t gam)
+int ConductionOperator::ARKImplicitSetup(const Vector &u, const real_t t,
+                                         const Vector &fu, int jok, int *jcur,
+                                         real_t gam)
 {
    // Compute T = M + gamma K(u_n).
    T = std::unique_ptr<HypreParMatrix>(Add(1.0, Mmat, gam, Kmat));
@@ -579,22 +632,22 @@ int ConductionOperator::SUNImplicitSetup(const Vector &u, const Vector &fu,
    return SUN_SUCCESS;
 }
 
-int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
+int ConductionOperator::ARKImplicitSolve(const Vector &r, Vector &dk,
                                          real_t tol)
 {
    // Solve the system [M + gamma K(u_n)] dk = - K(u_n) u - M k.
    // What value r is providing depends on the ODE expression form:
-   //   EXPLICIT form: r = -inv(M) K(u_n) u - k
-   //   IMPLICIT form: r = -K(u_n) u - M k
+   //   MFEM form: r = -inv(M) K(u_n) u - k
+   //   mass form: r = -K(u_n) u - M k
    T_solver.SetRelTol(tol);
-   if (isExplicit())
+   if (use_mass_form)
    {
-      Mmat.Mult(r, z);
-      T_solver.Mult(z, dk);
+      T_solver.Mult(r, dk);
    }
    else
    {
-      T_solver.Mult(r, dk);
+      Mmat.Mult(r, z);
+      T_solver.Mult(z, dk);
    }
    if (T_solver.GetConverged())
    {
@@ -606,13 +659,13 @@ int ConductionOperator::SUNImplicitSolve(const Vector &r, Vector &dk,
    }
 }
 
-int ConductionOperator::SUNMassSetup()
+int ConductionOperator::ARKMassSetup(const real_t t)
 {
    // Do nothing b/c mass solver was setup in constructor.
    return SUN_SUCCESS;
 }
 
-int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
+int ConductionOperator::ARKMassSolve(const Vector &b, Vector &x, real_t tol)
 {
    // Solve the system M x = b.
    M_solver.SetRelTol(tol);
@@ -627,7 +680,7 @@ int ConductionOperator::SUNMassSolve(const Vector &b, Vector &x, real_t tol)
    }
 }
 
-int ConductionOperator::SUNMassMult(const Vector &x, Vector &v)
+int ConductionOperator::ARKMassMult(const Vector &x, Vector &v)
 {
    // Compute M x.
    Mmat.Mult(x, v);

examples/sundials/ex9.cpp

@@ -119,7 +119,7 @@ public:
     and advection matrices, and b describes the flow on the boundary. This can
     be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
     used to evaluate the right-hand side. */
-class FE_Evolution : public TimeDependentOperator
+class FE_Evolution : public TimeDependentOperator, public ARKStepODE
 {
 private:
    BilinearForm &M, &K;
@@ -133,9 +133,14 @@ private:
 public:
    FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
 
+   // TimeDependentOperator methods for MFEM native and CVODE time integrators
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
 
+   // ARKStepODE methods for ARKODE time integrators
+   int ARKSize() const override;
+   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector& result) const override;
+
    virtual ~FE_Evolution();
 };
 
@@ -404,14 +409,14 @@ int main(int argc, char *argv[])
          ode_solver = cvode; break;
       case 8:
          arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
-         arkode->Init(adv);
+         arkode->Init(&adv);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
          arkode->SetOrder(4);
          ode_solver = arkode; break;
       case 9:
          arkode = new ARKStepSolver(ARKStepSolver::EXPLICIT);
-         arkode->Init(adv);
+         arkode->Init(&adv);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
          arkode->SetERKTableNum(ARKODE_FEHLBERG_13_7_8);
@@ -520,6 +525,19 @@ void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
    dg_solver->Mult(z, k);
 }
 
+int FE_Evolution::ARKSize() const
+{
+   return z.Size();
+}
+
+void FE_Evolution::ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(u, z);
+   z += b;
+   M_solver.Mult(z, result);
+}
+
 FE_Evolution::~FE_Evolution()
 {
    delete M_prec;

examples/sundials/ex9p.cpp

@@ -206,7 +206,7 @@ public:
     and advection matrices, and b describes the flow on the boundary. This can
     be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
     used to evaluate the right-hand side. */
-class FE_Evolution : public TimeDependentOperator
+class FE_Evolution : public TimeDependentOperator, public ARKStepODE
 {
 private:
    OperatorHandle M, K;
@@ -221,9 +221,14 @@ public:
    FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, const Vector &b_,
                 PrecType prec_type);
 
+   // TimeDependentOperator methods for MFEM native and CVODE time integrators
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
 
+   // ARKStepODE methods for ARKODE time integrators
+   int ARKSize() const override;
+   void ARKEvaluateRHS(const Vector &u, const real_t t, Vector& result) const override;
+
    virtual ~FE_Evolution();
 };
 
@@ -575,7 +580,7 @@ int main(int argc, char *argv[])
       case 8:
       case 9:
          arkode = new ARKStepSolver(MPI_COMM_WORLD, ARKStepSolver::EXPLICIT);
-         arkode->Init(adv);
+         arkode->Init(&adv);
          arkode->SetSStolerances(reltol, abstol);
          arkode->SetMaxStep(dt);
          if (ode_solver_type == 9)
@@ -743,6 +748,19 @@ void FE_Evolution::Mult(const Vector &x, Vector &y) const
    M_solver.Mult(z, y);
 }
 
+int FE_Evolution::ARKSize() const
+{
+   return z.Size();
+}
+
+void FE_Evolution::ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const
+{
+   // y = M^{-1} (K x + b)
+   K->Mult(u, z);
+   z += b;
+   M_solver.Mult(z, result);
+}
+
 FE_Evolution::~FE_Evolution()
 {
    delete M_prec;

general/array.cpp

@@ -0,0 +1,214 @@
+// Copyright (c) 2010-2025, 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.
+
+
+// Abstract array data type
+
+#include "array.hpp"
+#include "../general/forall.hpp"
+#include <fstream>
+#include <type_traits>
+
+namespace mfem
+{
+
+template <class T>
+void Array<T>::Print(std::ostream &os, int width) const
+{
+   for (int i = 0; i < size; i++)
+   {
+      os << data[i];
+      if ( !((i+1) % width) || i+1 == size )
+      {
+         os << '\n';
+      }
+      else
+      {
+         os << " ";
+      }
+   }
+}
+
+template <class T>
+void Array<T>::Save(std::ostream &os, int fmt) const
+{
+   if (fmt == 0)
+   {
+      os << size << '\n';
+   }
+   for (int i = 0; i < size; i++)
+   {
+      os << operator[](i) << '\n';
+   }
+}
+
+template <class T>
+void Array<T>::Load(std::istream &in, int fmt)
+{
+   if (fmt == 0)
+   {
+      int new_size;
+      in >> new_size;
+      SetSize(new_size);
+   }
+   for (int i = 0; i < size; i++)
+   {
+      in >> operator[](i);
+   }
+}
+
+template <class T>
+T Array<T>::Max() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T max = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (max < operator[](i))
+      {
+         max = operator[](i);
+      }
+   }
+
+   return max;
+}
+
+template <class T>
+T Array<T>::Min() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T min = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (operator[](i) < min)
+      {
+         min = operator[](i);
+      }
+   }
+
+   return min;
+}
+
+// Partial Sum
+template <class T>
+void Array<T>::PartialSum()
+{
+   T sum = static_cast<T>(0);
+   for (int i = 0; i < size; i++)
+   {
+      sum+=operator[](i);
+      operator[](i) = sum;
+   }
+}
+
+template <class T>
+void Array<T>::Abs()
+{
+   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
+   const bool useDevice = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(useDevice);
+   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = std::abs(y[i]);
+   });
+}
+
+// Sum
+template <class T>
+T Array<T>::Sum() const
+{
+   T sum = static_cast<T>(0);
+   for (int i = 0; i < size; i++)
+   {
+      sum+=operator[](i);
+   }
+
+   return sum;
+}
+
+template <class T>
+int Array<T>::IsSorted() const
+{
+   T val_prev = operator[](0), val;
+   for (int i = 1; i < size; i++)
+   {
+      val=operator[](i);
+      if (val < val_prev)
+      {
+         return 0;
+      }
+      val_prev = val;
+   }
+
+   return 1;
+}
+
+template <class T>
+bool Array<T>::IsConstant() const
+{
+   if (size < 2) { return true; }
+   const T v0 = data[0];
+   for (int i = 1; i < size; i++)
+   {
+      if (data[i] != v0)
+      {
+         return false;
+      }
+   }
+
+   return true;
+}
+
+template <class T>
+void Array2D<T>::Load(const char *filename, int fmt)
+{
+   std::ifstream in;
+   in.open(filename, std::ifstream::in);
+   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
+   Load(in, fmt);
+   in.close();
+}
+
+template <class T>
+void Array2D<T>::Print(std::ostream &os, int width_)
+{
+   int height = this->NumRows();
+   int width  = this->NumCols();
+
+   for (int i = 0; i < height; i++)
+   {
+      os << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         os << (*this)(i,j);
+         if ( (j+1) == width_ || (j+1) % width_ == 0 )
+         {
+            os << '\n';
+         }
+         else
+         {
+            os << ' ';
+         }
+      }
+   }
+}
+
+template class Array<char>;
+template class Array<int>;
+template class Array<long long>;
+template class Array<real_t>;
+template class Array2D<int>;
+template class Array2D<real_t>;
+
+} // namespace mfem

general/array.hpp

@@ -16,13 +16,9 @@
 #include "mem_manager.hpp"
 #include "device.hpp"
 #include "error.hpp"
-#include "forall.hpp"
 #include "globals.hpp"
-#include "reducers.hpp"
-#include "scan.hpp"
 
 #include <iostream>
-#include <fstream>
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
@@ -139,8 +135,6 @@ public:
    /// Return the device flag of the Memory object used by the Array
    bool UseDevice() const { return data.UseDevice(); }
 
-   void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
-
    /// Return true if the data will be deleted by the Array
    inline bool OwnsData() const { return data.OwnsHostPtr(); }
 
@@ -281,11 +275,11 @@ public:
 
    /** @brief Find the maximal element in the array, using the comparison
        operator `<` for class T. */
-   inline T Max() const;
+   T Max() const;
 
    /** @brief Find the minimal element in the array, using the comparison
        operator `<` for class T. */
-   inline T Min() const;
+   T Min() const;
 
    /// Sorts the array in ascending order. This requires operator< to be defined for T.
    void Sort() { std::sort((T*)data, data + size); }
@@ -303,22 +297,22 @@ public:
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   inline int IsSorted() const;
+   int IsSorted() const;
 
    /// Does the Array have Size zero.
    bool IsEmpty() const { return Size() == 0; }
 
    /// Return true if all entries of the array are the same.
-   inline bool IsConstant() const;
+   bool IsConstant() const;
 
    /// Fill the entries of the array with the cumulative sum of the entries.
-   inline void PartialSum();
+   void PartialSum();
 
    /// Replace each entry of the array with its absolute value.
-   inline void Abs();
+   void Abs();
 
    /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   inline T Sum() const;
+   T Sum() const;
 
    /// Set all entries of the array to the provided constant.
    inline void operator=(const T &a);
@@ -803,14 +797,8 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
    SetSize(src.Size());
-
-   const bool use_dev = UseDevice() || src.UseDevice();
-   const auto x = src.Read(use_dev);
-   auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = x[i];
-   });
+   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+   return *this;
 }
 
 template <class T>
@@ -1026,24 +1014,19 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
    sa.SetSize(sa_size);
-   const bool use_dev = UseDevice() || sa.UseDevice();
-   const auto x = Read(use_dev);
-   auto y = sa.Write(use_dev);
-   mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
+   for (int i = 0; i < sa_size; i++)
    {
-      y[i] = x[offset + i];
-   });
+      sa[i] = (*this)[offset+i];
+   }
 }
 
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   const bool use_dev = UseDevice();
-   auto x = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   for (int i = 0; i < size; i++)
    {
-      x[i] = a;
-   });
+      data[i] = a;
+   }
 }
 
 template <class T>
@@ -1052,153 +1035,6 @@ inline void Array<T>::Assign(const T *p)
    data.CopyFromHost(p, Size());
 }
 
-template <class T>
-inline void Array<T>::Print(std::ostream &os, int width) const
-{
-   for (int i = 0; i < size; i++)
-   {
-      os << data[i];
-      if ( !((i+1) % width) || i+1 == size )
-      {
-         os << '\n';
-      }
-      else
-      {
-         os << " ";
-      }
-   }
-}
-
-template <class T>
-inline void Array<T>::Save(std::ostream &os, int fmt) const
-{
-   if (fmt == 0)
-   {
-      os << size << '\n';
-   }
-   for (int i = 0; i < size; i++)
-   {
-      os << operator[](i) << '\n';
-   }
-}
-
-template <class T>
-void Array<T>::Load(std::istream &in, int fmt)
-{
-   if (fmt == 0)
-   {
-      int new_size;
-      in >> new_size;
-      SetSize(new_size);
-   }
-   for (int i = 0; i < size; i++)
-   {
-      in >> operator[](i);
-   }
-}
-
-template <class T>
-inline T Array<T>::Max() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T max = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (max < operator[](i))
-      {
-         max = operator[](i);
-      }
-   }
-
-   return max;
-}
-
-template <class T>
-inline T Array<T>::Min() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T min = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (operator[](i) < min)
-      {
-         min = operator[](i);
-      }
-   }
-
-   return min;
-}
-
-// Partial Sum
-template <class T>
-inline void Array<T>::PartialSum()
-{
-   auto data_ptr = ReadWrite(UseDevice());
-   InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
-}
-
-template <class T>
-inline void Array<T>::Abs()
-{
-   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
-   const bool useDevice = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = std::abs(y[i]);
-   });
-}
-
-// Sum
-template <class T>
-inline T Array<T>::Sum() const
-{
-   T sum = static_cast<T>(0);
-   if (size > 0)
-   {
-      const auto m_data = Read(UseDevice());
-      reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
-      /*  */ SumReducer<T> {}, UseDevice());
-   }
-   return sum;
-}
-
-template <class T>
-inline int Array<T>::IsSorted() const
-{
-   T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
-   {
-      val=operator[](i);
-      if (val < val_prev)
-      {
-         return 0;
-      }
-      val_prev = val;
-   }
-
-   return 1;
-}
-
-template <class T>
-inline bool Array<T>::IsConstant() const
-{
-   if (size < 2) { return true; }
-   const T v0 = data[0];
-   for (int i = 1; i < size; i++)
-   {
-      if (data[i] != v0)
-      {
-         return false;
-      }
-   }
-
-   return true;
-}
-
 
 template <class T>
 inline const T &Array2D<T>::operator()(int i, int j) const
@@ -1238,40 +1074,6 @@ inline T *Array2D<T>::operator[](int i)
    return &array1d[i*N];
 }
 
-template <class T>
-void Array2D<T>::Load(const char *filename, int fmt)
-{
-   std::ifstream in;
-   in.open(filename, std::ifstream::in);
-   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
-   Load(in, fmt);
-   in.close();
-}
-
-template <class T>
-void Array2D<T>::Print(std::ostream &os, int width_)
-{
-   int height = this->NumRows();
-   int width  = this->NumCols();
-
-   for (int i = 0; i < height; i++)
-   {
-      os << "[row " << i << "]\n";
-      for (int j = 0; j < width; j++)
-      {
-         os << (*this)(i,j);
-         if ( (j+1) == width_ || (j+1) % width_ == 0 )
-         {
-            os << '\n';
-         }
-         else
-         {
-            os << ' ';
-         }
-      }
-   }
-}
-
 
 template <class T>
 inline void Swap(Array2D<T> &a, Array2D<T> &b)

general/reducers.hpp

@@ -12,6 +12,7 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP
 
+#include "array.hpp"
 #include "forall.hpp"
 
 #include <cmath>
@@ -513,33 +514,6 @@ template<class B, class R> struct reduction_kernel
       }
    }
 };
-
-template <class T>
-class ReductionWorkspace
-{
-   Memory<T> workspace;
-
-   static ReductionWorkspace &Instance()
-   {
-      static ReductionWorkspace instance;
-      return instance;
-   }
-
-   ~ReductionWorkspace() { workspace.Delete(); }
-
-public:
-   static T *Get(int num_blocks)
-   {
-      ReductionWorkspace &instance = Instance();
-      if (instance.workspace.Capacity() < num_blocks)
-      {
-         instance.workspace.Delete();
-         instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
-      }
-      return instance.workspace;
-   }
-};
-
 }
 
 /**
@@ -555,7 +529,8 @@ public:
  @tparam T value_type to operate on
  */
 template <class T, class B, class R>
-void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
+            Array<T> &workspace)
 {
    if (N == 0)
    {
@@ -592,7 +567,13 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
 
       red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
       // allocate res to fit block_size entries
-      auto work = internal::ReductionWorkspace<T>::Get(nblocks);
+      auto mt = workspace.GetMemory().GetMemoryType();
+      if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
+      {
+         mt = MemoryType::HOST_PINNED;
+      }
+      workspace.SetSize(nblocks, mt);
+      auto work = workspace.HostWrite();
       red.work = work;
       forall_2D(nblocks, block_size, 1, std::move(red));
       // wait for results

general/scan.hpp

@@ -28,37 +28,8 @@
 
 namespace mfem
 {
-
-namespace internal
-{
-class ScanWorkspace
-{
-   Memory<std::byte> workspace;
-   static ScanWorkspace &Instance()
-   {
-      static ScanWorkspace instance;
-      return instance;
-   }
-   ~ScanWorkspace() { workspace.Delete(); }
-public:
-   static std::byte *Get(int num_bytes)
-   {
-      ScanWorkspace &instance = Instance();
-      if (Size() < num_bytes)
-      {
-         instance.workspace.Delete();
-         instance.workspace.New(num_bytes);
-      }
-      return instance.workspace.Write(MemoryClass::DEVICE, Size());
-   }
-   static int Size()
-   {
-      return Instance().workspace.Capacity();
-   }
-};
-}
-
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
+/// std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -66,12 +37,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
+                       workspace.Write(), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -86,12 +57,11 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
       }
       // try allocating a larger buffer
       bytes = 0;
-      // get size of buffer
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
                         nullptr, bytes, d_in, d_out, num_items));
-      // resize buffer (in ScanWorkspace::Get) and try again
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
+                        workspace.Write(), bytes, d_in, d_out, num_items));
       return;
    }
 #endif
@@ -131,13 +101,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -154,9 +123,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, num_items));
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
       return;
    }
 #endif
@@ -195,13 +164,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
       if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       init_value, num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -218,9 +187,10 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        init_value, num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
       return;
    }
 #endif
@@ -243,7 +213,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 }
 
 /// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// std::plus<>{})
+/// workspace, std::plus<>{})
 template <class InputIt, class OutputIt, class T>
 void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
                    T init_value)

linalg/operator.hpp

@@ -578,7 +578,7 @@ public:
 
        Presently, this method is used by SUNDIALS ARKStep integrator, for more
        details, see the ARKode User Guide. */
-   virtual int SUNMassSetup();
+   MFEM_DEPRECATED virtual int SUNMassSetup();
 
    /** @brief Solve the mass matrix linear system  M @a x = @a b, where M is
        defined by the method SUNMassSetup().
@@ -591,7 +591,7 @@ public:
 
        Presently, this method is used by SUNDIALS ARKStep integrator, for more
        details, see the ARKode User Guide. */
-   virtual int SUNMassSolve(const Vector &b, Vector &x, real_t tol);
+   MFEM_DEPRECATED virtual int SUNMassSolve(const Vector &b, Vector &x, real_t tol);
 
    /** @brief Compute the mass matrix-vector product @a v = M @a x, where M is
        defined by the method SUNMassSetup().
@@ -603,7 +603,7 @@ public:
 
        Presently, this method is used by SUNDIALS ARKStep integrator, for more
        details, see the ARKode User Guide. */
-   virtual int SUNMassMult(const Vector &x, Vector &v);
+   MFEM_DEPRECATED virtual int SUNMassMult(const Vector &x, Vector &v);
 
    virtual ~TimeDependentOperator() { }
 };

linalg/sundials.cpp

@@ -1367,6 +1367,84 @@ CVODESSolver::~CVODESSolver()
 // ARKStep interface
 // ---------------------------------------------------------------------------
 
+ARKStepSolver::TimeDependentOperatorWrapper::TimeDependentOperatorWrapper(
+   TimeDependentOperator *f)
+{
+   tdo = f;
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKSize() const
+{
+   return tdo->Height();
+}
+
+bool ARKStepSolver::TimeDependentOperatorWrapper::ARKInMassForm() const
+{
+   return (tdo->isExplicit() == false);
+}
+
+void ARKStepSolver::TimeDependentOperatorWrapper::ARKSetEvalMode(
+   const ARKEvalMode new_eval_mode)
+{
+   if (new_eval_mode == NORMAL)
+      tdo->SetEvalMode(tdo->NORMAL);
+   else if (new_eval_mode == ADDITIVE_TERM_1)
+      tdo->SetEvalMode(tdo->ADDITIVE_TERM_1);
+   else if (new_eval_mode == ADDITIVE_TERM_2)
+      tdo->SetEvalMode(tdo->ADDITIVE_TERM_2);
+   else
+      mfem_error("Unrecognized evaluation mode.");
+}
+
+void ARKStepSolver::TimeDependentOperatorWrapper::ARKEvaluateRHS(
+   const Vector &u, const real_t t, Vector &result) const
+{
+   tdo->SetTime(t);
+   if (ARKInMassForm())
+      tdo->Mult(u, result);
+   else
+      tdo->ExplicitMult(u, result);
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKImplicitSetup(
+   const Vector &u, const real_t t, const Vector &v, int jok, int *jcur,
+   real_t gamma)
+{
+   tdo->SetTime(t);
+   return tdo->SUNImplicitSetup(u, v, jok, jcur, gamma);
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKImplicitSolve(
+   const Vector &r, Vector &dk, real_t tol)
+{
+   return tdo->SUNImplicitSolve(r, dk, tol);
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassSetup(const real_t t)
+{
+   tdo->SetTime(t);
+   return tdo->SUNMassSetup();
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassSolve(const Vector &b,
+   Vector &x, real_t tol)
+{
+   return tdo->SUNMassSolve(b, x, tol);
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassMult(const Vector &x,
+   const real_t t, Vector &v)
+{
+   tdo->SetTime(t);
+   return tdo->SUNMassMult(x, v);
+}
+
+int ARKStepSolver::TimeDependentOperatorWrapper::ARKMassMult(const Vector &x,
+   Vector &v)
+{
+   return tdo->SUNMassMult(x, v);
+}
+
 int ARKStepSolver::RHS1(sunrealtype t, const N_Vector y, N_Vector result,
                         void *user_data)
 {
@@ -1381,19 +1459,11 @@ int ARKStepSolver::RHS1(sunrealtype t, const N_Vector y, N_Vector result,
    // or fe(t, y) in one of
    //   1. y' = fe(t, y) + fi(t, y)
    //   2. M y' = fe(t, y) + fi(t, y)
-   self->f->SetTime(t);
    if (self->rk_type == IMEX)
    {
-      self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_1);
-   }
-   if (self->f->isExplicit()) // ODE is in form 1
-   {
-      self->f->Mult(mfem_y, mfem_result);
-   }
-   else // ODE is in form 2
-   {
-      self->f->ExplicitMult(mfem_y, mfem_result);
+      self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_1);
    }
+   self->f_arkstep->ARKEvaluateRHS(mfem_y, t, mfem_result);
 
    // Return success
    return (0);
@@ -1410,16 +1480,8 @@ int ARKStepSolver::RHS2(sunrealtype t, const N_Vector y, N_Vector result,
    // Compute fi(t, y) in one of
    //   1. y' = fe(t, y) + fi(t, y)       (ODE is expressed in EXPLICIT form)
    //   2. M y' = fe(t, y) + fi(y, t)     (ODE is expressed in IMPLICIT form)
-   self->f->SetTime(t);
-   self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
-   if (self->f->isExplicit())
-   {
-      self->f->Mult(mfem_y, mfem_result);
-   }
-   else
-   {
-      self->f->ExplicitMult(mfem_y, mfem_result);
-   }
+   self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_2);
+   self->f_arkstep->ARKEvaluateRHS(mfem_y, t, mfem_result);
 
    // Return success
    return (0);
@@ -1436,12 +1498,11 @@ int ARKStepSolver::LinSysSetup(sunrealtype t, N_Vector y, N_Vector fy,
    ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(A));
 
    // Compute the linear system
-   self->f->SetTime(t);
    if (self->rk_type == IMEX)
    {
-      self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
+      self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_2);
    }
-   return (self->f->SUNImplicitSetup(mfem_y, mfem_fy, jok, jcur, gamma));
+   return (self->f_arkstep->ARKImplicitSetup(mfem_y, t, mfem_fy, jok, jcur, gamma));
 }
 
 int ARKStepSolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
@@ -1454,9 +1515,9 @@ int ARKStepSolver::LinSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
    // Solve the linear system
    if (self->rk_type == IMEX)
    {
-      self->f->SetEvalMode(TimeDependentOperator::ADDITIVE_TERM_2);
+      self->f_arkstep->ARKSetEvalMode(ARKStepODE::ADDITIVE_TERM_2);
    }
-   return (self->f->SUNImplicitSolve(mfem_b, mfem_x, tol));
+   return (self->f_arkstep->ARKImplicitSolve(mfem_b, mfem_x, tol));
 }
 
 int ARKStepSolver::MassSysSetup(sunrealtype t, SUNMatrix M,
@@ -1465,8 +1526,7 @@ int ARKStepSolver::MassSysSetup(sunrealtype t, SUNMatrix M,
    ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(M));
 
    // Compute the mass matrix system
-   self->f->SetTime(t);
-   return (self->f->SUNMassSetup());
+   return (self->f_arkstep->ARKMassSetup(t));
 }
 
 int ARKStepSolver::MassSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
@@ -1477,7 +1537,7 @@ int ARKStepSolver::MassSysSolve(SUNLinearSolver LS, SUNMatrix, N_Vector x,
    ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(LS));
 
    // Solve the mass matrix system
-   return (self->f->SUNMassSolve(mfem_b, mfem_x, tol));
+   return (self->f_arkstep->ARKMassSolve(mfem_b, mfem_x, tol));
 }
 
 int ARKStepSolver::MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
@@ -1487,7 +1547,7 @@ int ARKStepSolver::MassMult1(SUNMatrix M, N_Vector x, N_Vector v)
    ARKStepSolver *self = static_cast<ARKStepSolver*>(GET_CONTENT(M));
 
    // Compute the mass matrix-vector product
-   return (self->f->SUNMassMult(mfem_x, mfem_v));
+   return (self->f_arkstep->ARKMassMult(mfem_x, mfem_v));
 }
 
 int ARKStepSolver::MassMult2(N_Vector x, N_Vector v, sunrealtype t,
@@ -1498,8 +1558,7 @@ int ARKStepSolver::MassMult2(N_Vector x, N_Vector v, sunrealtype t,
    ARKStepSolver *self = static_cast<ARKStepSolver*>(mtimes_data);
 
    // Compute the mass matrix-vector product
-   self->f->SetTime(t);
-   return (self->f->SUNMassMult(mfem_x, mfem_v));
+   return (self->f_arkstep->ARKMassMult(mfem_x, t, mfem_v));
 }
 
 ARKStepSolver::ARKStepSolver(Type type)
@@ -1518,13 +1577,12 @@ ARKStepSolver::ARKStepSolver(MPI_Comm comm, Type type)
 }
 #endif
 
-void ARKStepSolver::Init(TimeDependentOperator &f_)
+void ARKStepSolver::Init(ARKStepODE *f_ark_)
 {
-   // Initialize the base class
-   ODESolver::Init(f_);
+   f_arkstep = f_ark_;
 
    // Get the vector length
-   long local_size = f_.Height();
+   long local_size = f_arkstep->ARKSize();
 #ifdef MFEM_USE_MPI
    long global_size;
 #endif
@@ -1538,7 +1596,7 @@ void ARKStepSolver::Init(TimeDependentOperator &f_)
    }
 
    // Get current time
-   double t = f_.GetTime();
+   double t = f ? f->GetTime() : 0.0;
 
    if (sundials_mem)
    {
@@ -1617,6 +1675,12 @@ void ARKStepSolver::Init(TimeDependentOperator &f_)
    reinit = true;
 }
 
+void ARKStepSolver::Init(TimeDependentOperator &f_)
+{
+   f_tdo = std::make_unique<TimeDependentOperatorWrapper>(&f_);
+   Init(f_tdo.get());
+}
+
 void ARKStepSolver::Step(Vector &x, real_t &t, real_t &dt)
 {
    Y->MakeRef(x, 0, x.Size());
@@ -1709,6 +1773,9 @@ void ARKStepSolver::UseSundialsLinearSolver()
 
 void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
 {
+   // Check that the ODE is expressed in mass form
+   MFEM_VERIFY(f_arkstep->ARKInMassForm(), "ODE operator is not in mass form.")
+
    // Free any existing matrix and linear solver
    if (M != NULL)   { SUNMatDestroy(M); M = NULL; }
    if (LSM != NULL) { SUNLinSolFree(LSM); LSM = NULL; }
@@ -1739,13 +1806,13 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
    flag = MFEM_ARKode(SetMassFn)(sundials_mem, ARKStepSolver::MassSysSetup);
    MFEM_VERIFY(flag == ARK_SUCCESS,
                "error in " STR(MFEM_ARKode(SetMassFn)) "()");
-
-   // Check that the ODE is not expressed in EXPLICIT form
-   MFEM_VERIFY(!f->isExplicit(), "ODE operator is expressed in EXPLICIT form")
 }
 
 void ARKStepSolver::UseSundialsMassLinearSolver(int tdep)
 {
+   // Check that the ODE is expressed in mass form
+   MFEM_VERIFY(f_arkstep->ARKInMassForm(), "ODE operator is not in mass form.")
+
    // Free any existing matrix and linear solver
    if (M != NULL)   { SUNMatDestroy(A); M = NULL; }
    if (LSM != NULL) { SUNLinSolFree(LSM); LSM = NULL; }
@@ -1764,9 +1831,6 @@ void ARKStepSolver::UseSundialsMassLinearSolver(int tdep)
                                     ARKStepSolver::MassMult2, this);
    MFEM_VERIFY(flag == ARK_SUCCESS,
                "error in " STR(MFEM_ARKode(SetMassTimes)) "()");
-
-   // Check that the ODE is not expressed in EXPLICIT form
-   MFEM_VERIFY(!f->isExplicit(), "ODE operator is expressed in EXPLICIT form")
 }
 
 void ARKStepSolver::SetStepMode(int itask)

linalg/sundials.hpp

@@ -706,9 +706,130 @@ public:
 // Interface to ARKode's ARKStep module -- Additive Runge-Kutta methods
 // ---------------------------------------------------------------------------
 
+// Interface for defining ODE systems to be evolved using ARKStepSolver:
+//
+//  1)   du/dt = inv(M) f(u,t) ("MFEM" form)
+//  2) M dy/dt = f(u,t)        ("mass" form)
+//
+// where f(u,t) might be additively split, i.e., f(u,t) = f1(u,t) + f2(u,t)
+class ARKStepODE
+{
+
+public:
+
+   // the size of the ODE system
+   virtual int ARKSize() const = 0;
+
+   // return if the ODE system is of the form M du/dt = f(u,t), note the MFEM
+   // default is to use the form du/dt = int(M) f(u,t)
+   virtual bool ARKInMassForm() const { return false; };
+
+   // these flags are used by ARKStepSolver for switching between RK and ARK methods
+   enum ARKEvalMode
+   {  NORMAL,          // evaluate f(u,t)
+      ADDITIVE_TERM_1, // evaluate f1(u,t)
+      ADDITIVE_TERM_2  // evaluate f2(u,t)
+   };
+   virtual void ARKSetEvalMode(const ARKEvalMode new_eval_mode) {}
+
+   // evaluate either f(u,t) (mass form) or inv(M(t)) f(u,t) (MFEM form),
+   // which is necessary for solving ODEs with ERK or IMEX
+   virtual void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const
+   {
+      mfem_error("This function must be specified for ERK or IMEX methods.");
+   }
+
+  /** setup linear system for solving [M(t) - gamma Jf(u)] dk = f(u) - M(t) k,
+      which is necessary for solving ODEs with DIRK or IMEX methods
+      @param[in]  u     The state at which A(@a u,t) should be evaluated.
+      @param[in]  t     The time at which A(u,@a t) should be evaluated.
+      @param[in]  v     The value of inv(M) f(u,t) or f(u,t) for depending on form.
+      @param[in]  jok   Flag indicating if the Jacobian should be updated.
+      @param[out] jcur  Flag to signal if the Jacobian was updated.
+      @param[in]  gamma The scaled time step value.    */
+   virtual int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &v,
+                                int jok, int *jcur, real_t gamma)
+   {
+      mfem_error("This function must be specified for DIRK or IMEX methods.");
+   }
+
+  /** solve for dk in [M - gamma Jf(u)] dk = r, where r is either
+        inv(M) f(u,t) - k       (MFEM form)
+        f(u,t) - M k f(u) - M k (mass form)
+      when using DIRK or IMEX methods
+      @param[in]      r   inv(M) f(u,t) - k or f(u,t) - M k, depending on form.
+      @param[in,out]  dk  On input, the initial guess. On output, the solution.
+      @param[in]      tol Linear solve tolerance. */
+   virtual int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol)
+   {
+      mfem_error("This function must be specified for DIRK or IMEX methods.");
+   }
+
+  /** for mass form ODEs using an MFEM mass solver, setup the mass linear
+      system M(t) x = b
+      @param[in]  t  The time at which M(@a t) should be evaluated. */
+   virtual int ARKMassSetup(const real_t t)
+   {
+      mfem_error("This function must be specified to use MFEM mass solvers for mass form ODEs.");
+   }
+
+  /** for mass form ODEs using an MFEM mass solver, solve for x in M(t) x = b
+      @param[in]      b   The linear system right-hand side.
+      @param[in,out]  x   On input, the initial guess. On output, the solution.
+      @param[in]      tol Linear solve tolerance. */
+   virtual int ARKMassSolve(const Vector &b, Vector &x, real_t tol)
+   {
+      mfem_error("This function must be specified to use MFEM mass solver for mass form ODEs.");
+   }
+
+  /**  for mass form ODEs using an MFEM mass solver, evaluate M(t) x
+       @param[in]   x The vector to multiply.
+       @param[out]  v The result of the matrix-vector product. */
+   virtual int ARKMassMult(const Vector &x, Vector &v)
+   {
+      mfem_error("This function must be specified to use MFEM mass solver for mass form ODEs.");
+   }
+
+  /**  for mass form ODEs using a SUNDIALS mass solver, evaluate M(t) x
+       @param[in]   x The vector to multiply.
+       @param[in]   t The time at which M(@a t) should be evaluated.
+       @param[out]  v The result of the matrix-vector product. */
+   virtual int ARKMassMult(const Vector &x, const real_t t, Vector &v)
+   {
+      mfem_error("This function must be specified to use SUNDIALS mass solver for mass form ODEs.");
+   }
+
+};
+
 /// Interface to ARKode's ARKStep module -- additive Runge-Kutta methods.
 class ARKStepSolver : public ODESolver, public SundialsSolver
 {
+
+   // Wrapper class to provide backwards compatability with user code that
+   // derives from TimeDependentOperator instead of ARKStepODE
+   class TimeDependentOperatorWrapper : public ARKStepODE
+   {
+      TimeDependentOperator *tdo;
+
+   public:
+
+      TimeDependentOperatorWrapper(TimeDependentOperator *f);
+
+      int ARKSize() const override;
+      bool ARKInMassForm() const override;
+      void ARKSetEvalMode(const ARKEvalMode new_eval_mode) override;
+      void ARKEvaluateRHS(const Vector &u, const real_t t, Vector &result) const override;
+      int ARKImplicitSetup(const Vector &u, const real_t t, const Vector &v,
+         int jok, int *jcur, real_t gamma) override;
+      int ARKImplicitSolve(const Vector &r, Vector &dk, real_t tol) override;
+      int ARKMassSetup(const real_t t) override;
+      int ARKMassSolve(const Vector &b, Vector &x, real_t tol) override;
+      int ARKMassMult(const Vector &x, Vector &v) override;
+      int ARKMassMult(const Vector &x, const real_t t, Vector &v) override;
+
+   };
+
+
 public:
    /// Types of ARKODE solvers.
    enum Type
@@ -722,6 +843,8 @@ protected:
    Type rk_type;      ///< Runge-Kutta type.
    int step_mode;     ///< ARKStep step mode (ARK_NORMAL or ARK_ONE_STEP).
    bool use_implicit; ///< True for implicit or imex integration.
+   ARKStepODE* f_arkstep;
+   std::unique_ptr<TimeDependentOperatorWrapper> f_tdo; // for backwards compatibility
 
    /** @name Wrappers to compute the ODE RHS functions.
        RHS1 is explicit RHS and RHS2 the implicit RHS for IMEX integration. When
@@ -784,14 +907,19 @@ public:
        then ARKStepReInit() will be called in the next call to Step(). If the
        problem size has changed, the ARKStep memory is freed and realloced
        for the new problem size. */
-   /** @param[in] f_ The TimeDependentOperator that defines the ODE system
+   /** @param[in] f_ The ARKStepODE that defines the ODE system
 
        @note All other methods must be called after Init().
 
        @note If this method is called a second time with a different problem
        size, then any non-default user-set options will be lost and will need
        to be set again. */
-   void Init(TimeDependentOperator &f_) override;
+   void Init(ARKStepODE *f_ark_);
+
+   // This method is provided for backwards compatibility with classes that
+   // derive TimeDependentOperator instead of ARKStepODE; however, those classes
+   // should be migrated.
+   MFEM_DEPRECATED void Init(TimeDependentOperator &f_) override;
 
    /// Integrate the ODE with ARKode using the specified step mode.
    /**

linalg/vector.cpp

@@ -92,6 +92,18 @@ struct LpReducer
    }
 };
 
+static Array<real_t>& vector_workspace()
+{
+   static Array<real_t> instance;
+   return instance;
+}
+
+static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
+{
+   static Array<DevicePair<real_t, real_t>> instance;
+   return instance;
+}
+
 Vector::Vector(const Vector &v)
 {
    const int s = v.Size();
@@ -979,7 +991,7 @@ real_t Vector::Norml2() const
          }
       }
    },
-   L2Reducer{}, UseDevice());
+   L2Reducer{}, UseDevice(), Lpvector_workspace());
    // final answer
    return res.second * sqrt(res.first);
 }
@@ -994,7 +1006,7 @@ real_t Vector::Normlinf() const
    {
       r = fmax(r, fabs(m_data[i]));
    },
-   MaxReducer<real_t> {}, UseDevice());
+   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 
@@ -1008,7 +1020,7 @@ real_t Vector::Norml1() const
    {
       r += fabs(m_data[i]);
    },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 
@@ -1051,7 +1063,7 @@ real_t Vector::Normlp(real_t p) const
             }
          }
       },
-      LpReducer{p}, UseDevice());
+      LpReducer{p}, UseDevice(), Lpvector_workspace());
       // final answer
       return res.second * pow(res.first, 1.0 / p);
    } // end if p < infinity()
@@ -1084,7 +1096,7 @@ real_t Vector::operator*(const Vector &v) const
       {
          r += m_data[i] * v_data[i];
       },
-      SumReducer<real_t> {}, use_dev);
+      SumReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1155,7 +1167,7 @@ real_t Vector::Min() const
       {
          r = fmin(r, m_data[i]);
       },
-      MinReducer<real_t> {}, use_dev);
+      MinReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1201,7 +1213,7 @@ real_t Vector::Max() const
       {
          r = fmax(r, m_data[i]);
       },
-      MaxReducer<real_t> {}, use_dev);
+      MaxReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1236,7 +1248,7 @@ real_t Vector::Sum() const
    {
       r += m_data[i];
    },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 

tests/unit/general/test_reduction.cpp

@@ -22,6 +22,7 @@ using namespace mfem;
 
 TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -35,7 +36,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
       int res = 0;
       mfem::reduce(
       a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
-      SumReducer<int> {}, use_dev);
+      SumReducer<int> {}, use_dev, workspace);
       // correct for even-length summations
       int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
       CAPTURE(use_dev);
@@ -45,6 +46,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 {
+   Array<long long> workspace;
    Array<long long> a(64);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -62,7 +64,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev);
+         MultReducer<long long> {}, use_dev, workspace);
          long long expected = 0;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -74,7 +76,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev);
+         MultReducer<long long> {}, use_dev, workspace);
          long long expected = 21936950640377856;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -84,6 +86,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
    Array<unsigned> a(10);
    SECTION("{ Bit unset }")
    {
@@ -105,7 +108,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev);
+         BAndReducer<unsigned> {}, use_dev, workspace);
          CAPTURE(use_dev);
          REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
          REQUIRE((res & (1u << unset_bit)) == 0);
@@ -129,7 +132,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev);
+         BAndReducer<unsigned> {}, use_dev, workspace);
          CAPTURE(use_dev);
          REQUIRE(res == (1u << set_bit));
       }
@@ -138,6 +141,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
    Array<unsigned> a(0x210);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -153,7 +157,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
       mfem::reduce(
          a.Size(), res,
       [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
-      BOrReducer<unsigned> {}, use_dev);
+      BOrReducer<unsigned> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == 0x3ffu);
    }
@@ -161,6 +165,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -185,7 +190,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MinReducer<int> {}, use_dev);
+      MinReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == -10);
    }
@@ -193,6 +198,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -217,7 +223,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MaxReducer<int> {}, use_dev);
+      MaxReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == 999 - 10);
    }
@@ -225,6 +231,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 {
+   Array<DevicePair<int, int>> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -255,7 +262,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
             r.second = dptr[i];
          }
       },
-      MinMaxReducer<int> {}, use_dev);
+      MinMaxReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second == a.Size() - 11);
@@ -264,6 +271,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 {
+   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -289,7 +297,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMinReducer<double, int> {}, use_dev);
+      ArgMinReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second >= 0);
@@ -300,6 +308,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 {
+   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
 
    auto hptr = a.HostReadWrite();
@@ -328,7 +337,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMaxReducer<double, int> {}, use_dev);
+      ArgMaxReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == a.Size() - 11);
       REQUIRE(res.second >= 0);
@@ -339,6 +348,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
+   Array<MinMaxLocScalar<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -373,7 +383,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
             r.max_loc = i;
          }
       },
-      ArgMinMaxReducer<double, int> {}, use_dev);
+      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.min_val == -10);
       REQUIRE(res.min_loc >= 0);