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;

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.
    /**