Update examples/sundials/ex16.cpp

examples/ex_sconv_p.cpp

@@ -1,503 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static double freq = 0.5, kappa;
-static int dim;
-
-double u_func(const Vector &);
-
-enum SCA_TYPE {INVALID_SCA_TYPE = -1,
-               H1_TYPE = 0,
-               L2_TYPE,
-               L2I_TYPE,
-               NUM_SCA_TYPES
-              };
-enum CONV_TYPE {INVALID_CONV_TYPE = -1,
-                PROJECTION = 0,
-                INTERPOLATION_OP,
-                SOLVE,
-                SOLVE_W_DBC,
-                NUM_CONV_TYPES
-               };
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
-                                   FiniteElementCollection &fec);
-string GetTypeName(SCA_TYPE type);
-string GetConvTypeName(CONV_TYPE type);
-string GetConvTypeShortName(CONV_TYPE type);
-
-void Projection(const ParGridFunction &v0, ParGridFunction &v1);
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1);
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c);
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   MPI_Session mpi(argc, argv);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ser_ref_levels = 0;
-   int par_ref_levels = 0;
-   int order0 = 1;
-   int order1 = 1;
-   int type0 = 0;
-   int type1 = 1;
-   int conv_type = -1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order0, "-o0", "--initial-order",
-                  "Finite element order (polynomial degree) "
-                  "for initial field.");
-   args.AddOption(&order1, "-o1", "--final-order",
-                  "Finite element order (polynomial degree) "
-                  "for final field.");
-   args.AddOption(&type0, "-t0", "--initial-type",
-                  "Set the basis type for the initial field: "
-                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
-   args.AddOption(&type1, "-t1", "--final-type",
-                  "Set the basis type for the final field: "
-                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
-   args.AddOption(&conv_type, "-c", "--conversion-type",
-                  "Set the conversion scheme: "
-                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
-                  "3-Least Squares with BC, -1 loop over all.");
-   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
-                  " solution.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      if (mpi.Root()) { args.PrintUsage(cout); }
-      return 1;
-   }
-   if (mpi.Root()) { args.PrintOptions(cout); }
-   kappa = freq * M_PI;
-
-   // 3. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   if (mpi.Root()) { device.Print(); }
-
-   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   dim = mesh->Dimension();
-
-   // 5. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
-   //    specified on the command line with -rs).
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution (1 time by
-   //    default, or specified on the command line with -rp). Once the parallel
-   //    mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh.UniformRefinement();
-   }
-
-   FunctionCoefficient uCoef(u_func);
-
-   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
-
-   if (mpi.Root())
-   {
-      cout << "L2 Errors:" << endl;
-   }
-   int t0a = (type0 == -1) ? 0 : type0;
-   int t0b = (type0 == -1) ? NUM_SCA_TYPES : (type0+1);
-   for (int t0 = t0a; t0 < t0b; t0++)
-   {
-      FiniteElementCollection *fec0 = GetFECollection((SCA_TYPE)t0, order0);
-      ParFiniteElementSpace   *fes0 = GetFESpace((SCA_TYPE)t0, pmesh, *fec0);
-
-      ParGridFunction x0(fes0);
-      x0.ProjectCoefficient(uCoef);
-
-      double err0 = x0.ComputeL2Error(uCoef);
-      if (mpi.Root())
-      {
-         cout << "Initial " << GetTypeName((SCA_TYPE)t0)
-              << ": \t\t" << err0 << endl;
-      }
-
-      // nn. Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")";
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock0(vishost, visport);
-         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
-                   << pmesh.GetMyRank() << '\n';
-         sol_sock0.precision(8);
-         sol_sock0 << "solution\n" << pmesh << x0
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
-                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
-                   << flush;
-      }
-
-      int t1a = (type1 == -1) ? 0 : type1;
-      int t1b = (type1 == -1) ? NUM_SCA_TYPES : (type1+1);
-      for (int t1 = t1a; t1 < t1b; t1++)
-      {
-         FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, order1);
-         ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, pmesh, *fec1);
-
-         ParGridFunction y1(fes1);
-
-         if (mpi.Root())
-         {
-            cout << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")"
-                 << " -> "
-                 << GetTypeName((SCA_TYPE)t1) << "(" << order1 << ")"
-                 << ":" << endl;
-         }
-
-         int c01a = (conv_type == -1) ? 0 : conv_type;
-         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
-         for (int c01 = c01a; c01 < c01b; c01++)
-         {
-            string cmnt = "";
-
-            switch ((CONV_TYPE)c01)
-            {
-               case PROJECTION:
-                  Projection(x0, y1);
-                  break;
-               case INTERPOLATION_OP:
-                  cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
-                         "(should match projection)" : "(not expected to succeed)";
-                  InterpolationOp(x0, y1);
-                  break;
-               case SOLVE:
-                  LeastSquares((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1);
-                  break;
-               case SOLVE_W_DBC:
-                  LeastSquaresBC((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1, uCoef);
-                  break;
-               default:
-                  y1 = 0.0;
-            }
-
-            double err1 = y1.ComputeL2Error(uCoef);
-            cout << GetConvTypeName((CONV_TYPE)c01)
-                 << "\t\t" << err1 << "\t" << cmnt << endl;
-
-            if (visualization)
-            {
-               ostringstream oss;
-               oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")" << " --"
-                   << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
-                   << GetTypeName((SCA_TYPE)t1)<< "(" << order1 << ")";
-               char vishost[] = "localhost";
-               int  visport   = 19916;
-               socketstream sol_sock1(vishost, visport);
-               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
-                         << pmesh.GetMyRank() << '\n';
-               sol_sock1.precision(8);
-               sol_sock1 << "solution\n" << pmesh << y1
-                         << "window_title '" << oss.str() << "'"
-                         << "window_geometry "
-                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
-                                  Ws * (t0 - t0a))
-                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
-                         << " " << Ww << " " << Wh
-                         << flush;
-            }
-         }
-         if (mpi.Root())
-         {
-            cout << endl;
-         }
-
-         delete fes1;
-         delete fec1;
-      }
-
-      delete fes0;
-      delete fec0;
-
-      if (t0 < t0b - 1)
-      {
-         char c;
-         if (mpi.Root())
-         {
-            cout << "press (q)uit or (c)ontinue --> " << flush;
-            cin >> c;
-         }
-         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-         if (c != 'c')
-         {
-            break;
-         }
-      }
-      if (mpi.Root())
-      {
-         cout << endl;
-      }
-   }
-
-   return 0;
-}
-
-double u_func(const Vector &x)
-{
-   double kx = kappa * x[0];
-   double ky = kappa * x[1];
-   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
-
-   // Add the gradient of a scalar function
-   return cos(kx) * cos(ky) * cos(kz);
-}
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return new H1_FECollection(p, dim);
-      case L2_TYPE:
-         return new L2_FECollection(p-1, dim);
-      case L2I_TYPE:
-         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
-                                    FiniteElement::INTEGRAL);
-      default:
-         return NULL;
-   }
-}
-
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
-                                   ParMesh &pmesh,
-                                   FiniteElementCollection &fec)
-{
-   return new ParFiniteElementSpace(&pmesh, &fec);
-}
-
-string GetTypeName(SCA_TYPE type)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return "    H1";
-      case L2_TYPE:
-         return "    L2";
-      case L2I_TYPE:
-         return "    L2I";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Projection            ";
-      case INTERPOLATION_OP:
-         return "Interpolation Operator";
-      case SOLVE:
-         return "Least Squares         ";
-      case SOLVE_W_DBC:
-         return "Least Squares with BC ";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeShortName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Proj";
-      case INTERPOLATION_OP:
-         return "Interp";
-      case SOLVE:
-         return "LS";
-      case SOLVE_W_DBC:
-         return "LSwBC";
-      default:
-         return "--";
-   }
-}
-
-/** Perform a naive projection from one scalar field to another.
-
-   This scheme simply evaluates v0 at the interpolation points of v1.
-
-   If v0 has reduced continuity compared to v1 this can produce
-   results that depend on the order in which the elements are
-   traversed.
-
-   Suitable conversions:
-   H1 -> L2
-   H1 -> DG (same as L2)
- */
-void Projection(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   GridFunctionCoefficient v0Coef(&v0);
-   v1.ProjectCoefficient(v0Coef);
-}
-
-/** In theory this interpolation scheme should be equivalent to projection.
-
-    Building an interpolastion matrix could lead to computational
-    efficiency compared to simple projection if the operator will be
-    used several times.
-
-    Unfortunately this is broken for several combinations of source
-    and target fields.
-*/
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
-   op.AddDomainInterpolator(new IdentityInterpolator);
-   op.Assemble();
-   op.Finalize();
-
-   op.Mult(v0, v1);
-}
-
-/** Compute a least-squares best fit using the target basis functions.
-
-    This scheme is more difficult to setup and more computationally
-    expensive but the results can be significantly better than simple
-    projections.
-*/
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   HypreParMatrix * M = m.ParallelAssemble();
-
-   HypreDiagScale diag(*M);
-   HyprePCG pcg(*M);
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   Vector B, X;
-   b.ParallelAssemble(B);
-
-   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-
-   delete M;
-}
-
-/** Compute a least-squares best fit with boundary conditions.
-
-    This scheme is virtually identical to the previous one but it
-    makes use of boundary values, when available, to improve the
-    accuracy.  This scheme can produce significantly better results
-    when the normal derivative of the field is large near the
-    boundary.  This is particularly true when the field is
-    under-resolved near the boundary.
-*/
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
-   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
-      ess_bdr = 1;
-      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   if (t1 == H1_TYPE)
-   {
-      v1.ProjectBdrCoefficient(c, ess_bdr);
-   }
-
-   OperatorPtr M;
-   Vector B, X;
-   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
-
-   HypreDiagScale diag(*M.As<HypreParMatrix>());
-   HyprePCG pcg(*M.As<HypreParMatrix>());
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-}

examples/ex_vconv_p.cpp

@@ -1,586 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static double freq = 0.5, kappa;
-static int dim;
-
-void u_func(const Vector &, Vector &);
-
-enum VEC_TYPE {INVALID_VEC_TYPE = -1,
-               H1V_TYPE = 0,
-               ND_TYPE,
-               RT_TYPE,
-               L2V_TYPE,
-               NUM_VEC_TYPES
-              };
-enum CONV_TYPE {INVALID_CONV_TYPE = -1,
-                PROJECTION = 0,
-                INTERPOLATION_OP,
-                SOLVE,
-                SOLVE_W_DBC,
-                NUM_CONV_TYPES
-               };
-
-FiniteElementCollection * GetFECollection(VEC_TYPE type, int p);
-ParFiniteElementSpace * GetFESpace(VEC_TYPE type, ParMesh &pmesh,
-                                   FiniteElementCollection &fec);
-string GetTypeName(VEC_TYPE type);
-string GetConvTypeName(CONV_TYPE type);
-
-void Projection(const ParGridFunction &v0, ParGridFunction &v1);
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
-void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
-                  VEC_TYPE t1, ParGridFunction &v1);
-void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
-                    VEC_TYPE t1, ParGridFunction &v1,
-                    VectorCoefficient &vc);
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   MPI_Session mpi(argc, argv);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ser_ref_levels = 0;
-   int par_ref_levels = 0;
-   int order0 = 1;
-   int order1 = 1;
-   int type0 = 0;
-   int type1 = 1;
-   int conv_type = -1;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order0, "-o0", "--initial-order",
-                  "Finite element order (polynomial degree) "
-                  "for initial field.");
-   args.AddOption(&order1, "-o1", "--final-order",
-                  "Finite element order (polynomial degree) "
-                  "for final field.");
-   args.AddOption(&type0, "-t0", "--initial-type",
-                  "Set the basis type for the initial field: "
-                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
-   args.AddOption(&type1, "-t1", "--final-type",
-                  "Set the basis type for the final field: "
-                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
-   args.AddOption(&conv_type, "-c", "--conversion-type",
-                  "Set the conversion scheme: "
-                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
-                  "3-Least Squares with BC, -1 loop over all.");
-   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
-                  " solution.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      if (mpi.Root()) { args.PrintUsage(cout); }
-      return 1;
-   }
-   if (mpi.Root()) { args.PrintOptions(cout); }
-   kappa = freq * M_PI;
-
-   // 3. Enable hardware devices such as GPUs, and programming models such as
-   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
-   Device device(device_config);
-   if (mpi.Root()) { device.Print(); }
-
-   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   dim = mesh->Dimension();
-
-   // 5. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
-   //    specified on the command line with -rs).
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution (1 time by
-   //    default, or specified on the command line with -rp). Once the parallel
-   //    mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh.UniformRefinement();
-   }
-
-   VectorFunctionCoefficient uCoef(dim, u_func);
-
-   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
-
-   if (mpi.Root())
-   {
-      cout << "L2 Errors:" << endl;
-   }
-   int t0a = (type0 == -1) ? 0 : type0;
-   int t0b = (type0 == -1) ? NUM_VEC_TYPES : (type0+1);
-   for (int t0 = t0a; t0 < t0b; t0++)
-   {
-      FiniteElementCollection *fec0 = GetFECollection((VEC_TYPE)t0, order0);
-      ParFiniteElementSpace   *fes0 = GetFESpace((VEC_TYPE)t0, pmesh, *fec0);
-
-      ParGridFunction x0(fes0);
-      x0.ProjectCoefficient(uCoef);
-
-      double err0 = x0.ComputeL2Error(uCoef);
-      if (mpi.Root())
-      {
-         cout << "Initial " << GetTypeName((VEC_TYPE)t0)
-              << ": \t\t" << err0 << endl;
-      }
-
-      // nn. Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << GetTypeName((VEC_TYPE)t0);
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock0(vishost, visport);
-         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
-                   << pmesh.GetMyRank() << '\n';
-         sol_sock0.precision(8);
-         sol_sock0 << "solution\n" << pmesh << x0
-                   << "keys vvv "
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
-                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
-                   << flush;
-      }
-
-      int t1a = (type1 == -1) ? 0 : type1;
-      int t1b = (type1 == -1) ? NUM_VEC_TYPES : (type1+1);
-      for (int t1 = t1a; t1 < t1b; t1++)
-      {
-         FiniteElementCollection *fec1 = GetFECollection((VEC_TYPE)t1, order1);
-         ParFiniteElementSpace   *fes1 = GetFESpace((VEC_TYPE)t1, pmesh, *fec1);
-
-         ParGridFunction x1(fes1);
-
-         if (mpi.Root())
-         {
-            cout << GetTypeName((VEC_TYPE)t0) << " -> "
-                 << GetTypeName((VEC_TYPE)t1) << ":" << endl;
-         }
-
-         int c01a = (conv_type == -1) ? 0 : conv_type;
-         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
-         for (int c01 = c01a; c01 < c01b; c01++)
-         {
-            switch ((CONV_TYPE)c01)
-            {
-               case PROJECTION:
-                  Projection(x0, x1);
-                  break;
-               case INTERPOLATION_OP:
-                  // InterpolationOp(x0, x1);
-                  x1 = 0.0;
-                  break;
-               case SOLVE:
-                  LeastSquares((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1);
-                  break;
-               case SOLVE_W_DBC:
-                  LeastSquaresBC((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1, uCoef);
-                  break;
-               default:
-                  x1 = 0.0;
-            }
-
-            double err1 = x1.ComputeL2Error(uCoef);
-            cout << GetConvTypeName((CONV_TYPE)c01)
-                 << "\t\t" << err1 << endl;
-
-            if (visualization)
-            {
-               ostringstream oss;
-               oss << GetTypeName((VEC_TYPE)t0) << " --" << c01 << "--> "
-                   << GetTypeName((VEC_TYPE)t1);
-               char vishost[] = "localhost";
-               int  visport   = 19916;
-               socketstream sol_sock1(vishost, visport);
-               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
-                         << pmesh.GetMyRank() << '\n';
-               sol_sock1.precision(8);
-               sol_sock1 << "solution\n" << pmesh << x1
-                         << "keys vvv "
-                         << "window_title '" << oss.str() << "'"
-                         << "window_geometry "
-                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
-                                  Ws * (t0 - t0a))
-                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
-                         << " " << Ww << " " << Wh
-                         << flush;
-            }
-         }
-         if (mpi.Root())
-         {
-            cout << endl;
-         }
-
-         delete fes1;
-         delete fec1;
-      }
-
-      delete fes0;
-      delete fec0;
-
-      if (t0 < t0b - 1)
-      {
-         char c;
-         if (mpi.Root())
-         {
-            cout << "press (q)uit or (c)ontinue --> " << flush;
-            cin >> c;
-         }
-         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-         if (c != 'c')
-         {
-            break;
-         }
-      }
-      if (mpi.Root())
-      {
-         cout << endl;
-      }
-   }
-
-   return 0;
-}
-
-void u_func(const Vector &x, Vector &u)
-{
-   u.SetSize(dim);
-
-   double kx = kappa * x[0];
-   double ky = kappa * x[1];
-   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
-
-   // Add the gradient of a scalar function
-   u(0) = sin(kx) * cos(ky);
-   u(1) = cos(kx) * sin(ky);
-   if (dim == 3)
-   {
-      u(0) *= cos(kz);
-      u(1) *= cos(kz);
-      u(2) = cos(kx) * cos(ky) * sin(kz);
-   }
-
-   // Add the curl of a vector function
-   u(0) -= cos(kx) * sin(ky);
-   u(1) += sin(kx) * cos(ky);
-   if (dim == 3)
-   {
-      u(0) += cos(kx) * sin(kz);
-      u(1) -= cos(ky) * sin(kz);
-      u(2) += (sin(ky) - sin(kx)) * cos(kz);
-   }
-}
-
-FiniteElementCollection * GetFECollection(VEC_TYPE type, int p)
-{
-   switch (type)
-   {
-      case H1V_TYPE:
-         return new H1_FECollection(p, dim);
-      case ND_TYPE:
-         return new ND_FECollection(p, dim);
-      case RT_TYPE:
-         return new RT_FECollection(p-1, dim);
-      case L2V_TYPE:
-         return new L2_FECollection(p-1, dim);
-      default:
-         return NULL;
-   }
-}
-
-ParFiniteElementSpace * GetFESpace(VEC_TYPE type,
-                                   ParMesh &pmesh,
-                                   FiniteElementCollection &fec)
-{
-   switch (type)
-   {
-      case H1V_TYPE:
-      case L2V_TYPE:
-         return new ParFiniteElementSpace(&pmesh, &fec, dim);
-      case ND_TYPE:
-      case RT_TYPE:
-         return new ParFiniteElementSpace(&pmesh, &fec);
-      default:
-         return NULL;
-   }
-}
-
-string GetTypeName(VEC_TYPE type)
-{
-   switch (type)
-   {
-      case H1V_TYPE:
-         return "    H1V";
-      case ND_TYPE:
-         return "H(Curl)";
-      case RT_TYPE:
-         return " H(Div)";
-      case L2V_TYPE:
-         return "    L2V";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Projection            ";
-      case INTERPOLATION_OP:
-         return "Interpolation Operator";
-      case SOLVE:
-         return "Least Squares         ";
-      case SOLVE_W_DBC:
-         return "Least Squares with BC ";
-      default:
-         return "--";
-   }
-}
-
-/** Perform a naive projection from one vector field to another.
-
-   This scheme simply evaluates v0 at the interpolation points of v1.
-
-   If v0 has reduced continuity compared to v1 this can produce
-   results that depend on the order in which the elements are
-   traversed.
-
-   Suitable conversions:
-   H1V     -> H(Curl), H(Div), or L2V
-   H(Curl) -> L2V
-   H(Div)  -> L2V
-
- */
-void Projection(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   VectorGridFunctionCoefficient v0Coef(&v0);
-   v1.ProjectCoefficient(v0Coef);
-}
-
-/** In theory this interpolation scheme should be equivalent to projection.
-
-    Building an interpolastion matrix could lead to computational
-    efficiency compared to simple projection if the operator will be
-    used several times.
-
-    Unfortunately this is broken for several combinations of source
-    and target fields.
-*/
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
-   op.AddDomainInterpolator(new IdentityInterpolator);
-   op.Assemble();
-   op.Finalize();
-
-   op.Mult(v0, v1);
-}
-
-/** Compute a least-squares best fit using the target basis functions.
-
-    This scheme is more difficult to setup and more computationally
-    expensive but the results can be significantly better than simple
-    projections.
-*/
-void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
-                  VEC_TYPE t1, ParGridFunction &v1)
-{
-   bool trans = false;
-   ParFiniteElementSpace *fes0, *fes1;
-   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
-       (t1 == ND_TYPE || t1 == RT_TYPE))
-   {
-      fes0 = v1.ParFESpace();
-      fes1 = v0.ParFESpace();
-      trans = true;
-   }
-   else
-   {
-      fes0 = v0.ParFESpace();
-      fes1 = v1.ParFESpace();
-   }
-   ParMixedBilinearForm op(fes0, fes1);
-   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      op.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      op.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   if (trans)
-   {
-      op.MultTranspose(v0, b);
-   }
-   else
-   {
-      op.Mult(v0, b);
-   }
-
-   ParBilinearForm m(v1.ParFESpace());
-   if (t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      m.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      m.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   m.Assemble();
-   m.Finalize();
-
-   HypreParMatrix * M = m.ParallelAssemble();
-
-   HypreDiagScale diag(*M);
-   HyprePCG pcg(*M);
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   Vector B, X;
-   b.ParallelAssemble(B);
-
-   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-
-   delete M;
-}
-
-/** Compute a least-squares best fit with boundary conditions.
-
-    This scheme is virtually identical to the previous one but it
-    makes use of boundary values, when available, to improve the
-    accuracy.  This scheme can produce significantly better results
-    when the normal derivative of the field is large near the
-    boundary.  This is particularly true when the field is
-    under-resolved near the boundary.
-*/
-void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
-                    VEC_TYPE t1, ParGridFunction &v1,
-                    VectorCoefficient &vc)
-{
-   bool trans = false;
-   ParFiniteElementSpace *fes0, *fes1;
-   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
-       (t1 == ND_TYPE || t1 == RT_TYPE))
-   {
-      fes0 = v1.ParFESpace();
-      fes1 = v0.ParFESpace();
-      trans = true;
-   }
-   else
-   {
-      fes0 = v0.ParFESpace();
-      fes1 = v1.ParFESpace();
-   }
-   ParMixedBilinearForm op(fes0, fes1);
-   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      op.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      op.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   if (trans)
-   {
-      op.MultTranspose(v0, b);
-   }
-   else
-   {
-      op.Mult(v0, b);
-   }
-
-   ParBilinearForm m(v1.ParFESpace());
-   if (t1 == ND_TYPE || t1 == RT_TYPE)
-   {
-      m.AddDomainIntegrator(new VectorFEMassIntegrator);
-   }
-   else
-   {
-      m.AddDomainIntegrator(new VectorMassIntegrator);
-   }
-   m.Assemble();
-   m.Finalize();
-
-   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
-   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
-      ess_bdr = 1;
-      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   if (t1 == H1V_TYPE)
-   {
-      v1.ProjectBdrCoefficient(vc, ess_bdr);
-   }
-   if (t1 == ND_TYPE)
-   {
-      v1.ProjectBdrCoefficientTangent(vc, ess_bdr);
-   }
-   else if (t1 == RT_TYPE)
-   {
-      v1.ProjectBdrCoefficientNormal(vc, ess_bdr);
-   }
-
-   OperatorPtr M;
-   Vector B, X;
-   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
-
-   HypreDiagScale diag(*M.As<HypreParMatrix>());
-   HyprePCG pcg(*M.As<HypreParMatrix>());
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-}

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

miniapps/tools/gf-conv.cpp

@@ -1,577 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-static int dim;
-
-enum SCA_TYPE {INVALID_SCA_TYPE = -1,
-               H1_TYPE = 0,
-               L2_TYPE,
-               L2I_TYPE,
-               NUM_SCA_TYPES
-              };
-enum CONV_TYPE {INVALID_CONV_TYPE = -1,
-                PROJECTION = 0,
-                INTERPOLATION_OP,
-                SOLVE,
-                SOLVE_W_DBC,
-                NUM_CONV_TYPES
-               };
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
-                                   FiniteElementCollection &fec);
-void parseFieldNames(const char * field_name_c_str,
-                     vector<string> &field_names);
-
-string GetTypeName(SCA_TYPE type);
-string GetTypeShortName(SCA_TYPE type);
-string GetConvTypeName(CONV_TYPE type);
-string GetConvTypeShortName(CONV_TYPE type);
-
-void Projection(const ParGridFunction &v0, ParGridFunction &v1);
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1);
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c);
-
-int main(int argc, char *argv[])
-{
-#ifdef MFEM_USE_MPI
-   Mpi::Init();
-   if (!Mpi::Root()) { mfem::out.Disable(); mfem::err.Disable(); }
-   Hypre::Init();
-#endif
-
-   // Parse command-line options.
-   const char *coll_name = NULL;
-   int cycle = 0;
-
-   const char *field_name_c_str = "ALL";
-
-   Array<int> orders;
-   Array<int> types;
-   Array<int> conv_types;
-   bool static_cond = false;
-   bool pa = false;
-   const char *device_config = "cpu";
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&coll_name, "-r", "--root-file",
-                  "Set the VisIt data collection root file prefix.", true);
-   args.AddOption(&cycle, "-c", "--cycle", "Set the cycle index to read.");
-   args.AddOption(&field_name_c_str, "-fn", "--field-names",
-                  "List of field names to get values from.");
-   args.AddOption(&orders, "-o", "--final-order",
-                  "Finite element orders for each final field "
-                  "(an array of integers for multiple fields).");
-   args.AddOption(&types, "-t", "--final-type",
-                  "Set the basis type for the final fields: "
-                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
-   args.AddOption(&conv_types, "-ct", "--conversion-type",
-                  "Set the conversion schemes: "
-                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
-                  "3-Least Squares with BC, -1 loop over all.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
-                  "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(mfem::out);
-      return 1;
-   }
-   args.PrintOptions(mfem::out);
-
-#ifdef MFEM_USE_MPI
-   VisItDataCollection dc(MPI_COMM_WORLD, coll_name);
-#else
-   VisItDataCollection dc(coll_name);
-#endif
-
-   dc.Load(cycle);
-
-   if (dc.Error() != DataCollection::NO_ERROR)
-   {
-      mfem::out << "Error loading VisIt data collection: " << coll_name << endl;
-      return 1;
-   }
-
-   dim = dc.GetMesh()->Dimension();
-   int spaceDim = dc.GetMesh()->SpaceDimension();
-
-   mfem::out << endl;
-   mfem::out << "Collection Name:    " << dc.GetCollectionName() << endl;
-   mfem::out << "Manifold Dimension: " << dim << endl;
-   mfem::out << "Space Dimension:    " << spaceDim << endl;
-   mfem::out << "Cycle:              " << dc.GetCycle() << endl;
-   mfem::out << "Time:               " << dc.GetTime() << endl;
-   mfem::out << "Time Step:          " << dc.GetTimeStep() << endl;
-   mfem::out << endl;
-
-   typedef DataCollection::FieldMapType fields_t;
-   const fields_t &fields = dc.GetFieldMap();
-   // Print the names of all fields.
-   mfem::out << "fields: [ ";
-   for (fields_t::const_iterator it = fields.begin(); it != fields.end(); ++it)
-   {
-      if (it != fields.begin()) { mfem::out << ", "; }
-      mfem::out << it->first;
-   }
-   mfem::out << " ]" << endl;
-
-   // Parsing desired field names
-   vector<string> field_names;
-   parseFieldNames(field_name_c_str, field_names);
-
-   if (field_names.size() == 1)
-   {
-      if (field_names[0] == "ALL")
-      {
-         fields_t::const_iterator it = fields.begin();
-         field_names[0] = it->first; it++;
-         for ( ; it != fields.end(); ++it)
-         {
-            field_names.push_back(it->first);
-         }
-      }
-   }
-
-   if (orders.Size() < field_names.size())
-   {
-      int size = orders.Size();
-      int order = (size > 0) ? orders[0] : 1;
-
-      orders.SetSize(field_names.size());
-      for (int i=size; i < field_names.size(); i++)
-      {
-         orders[i] = order;
-      }
-   }
-
-   if (types.Size() < field_names.size())
-   {
-      int size = types.Size();
-      int type = (size > 0) ? types[0] : 0;
-
-      types.SetSize(field_names.size());
-      for (int i=size; i < field_names.size(); i++)
-      {
-         types[i] = type;
-      }
-   }
-
-   if (conv_types.Size() < field_names.size())
-   {
-      int size = conv_types.Size();
-      int type = (size > 0) ? conv_types[0] : 0;
-
-      conv_types.SetSize(field_names.size());
-      for (int i=size; i < field_names.size(); i++)
-      {
-         conv_types[i] = type;
-      }
-   }
-
-   // Print field names to be extracted
-   mfem::out << "Extracting fields: ";
-   for (int i=0; i < field_names.size(); i++)
-   {
-      mfem::out << " \"" << field_names[i] << "\"";
-   }
-   mfem::out << endl;
-
-#ifdef MFEM_USE_MPI
-   ParMesh *mesh = dynamic_cast<ParMesh*>(dc.GetMesh());
-#else
-   Mesh *mesh = dc.GetMesh();
-#endif
-   if (mesh == NULL)
-   {
-      mfem::out << "Problem with mesh\n";
-      return 1;
-   }
-
-   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
-
-   // Loop over all requested fields.
-   for (int i=0; i < field_names.size(); i++)
-   {
-#ifdef MFEM_USE_MPI
-      ParGridFunction *x0 = dc.GetParField(field_names[i]);
-#else
-      GridFunction *x0 = dc.GetField(field_names[i]);
-#endif
-      if (x0 == NULL)
-      {
-         mfem::out << "Problem with x0 for field \"" << field_names[i] << "\"\n";
-         continue;
-      }
-
-      int t0 = 0;
-
-      // nn. Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << field_names[i];
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock0(vishost, visport);
-#ifdef MFEM_USE_MPI
-         sol_sock0 << "parallel " << mesh->GetNRanks() << ' '
-                   << mesh->GetMyRank() << '\n';
-#endif
-         sol_sock0.precision(8);
-         sol_sock0 << "solution\n" << *mesh << *x0
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << Ws * (t0) << " " << Ws * (t0) << " "
-                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
-                   << flush;
-      }
-
-      int t1 = types[i];
-      FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, orders[i]);
-      ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, *mesh, *fec1);
-
-      ParGridFunction *y1 = new ParGridFunction(fes1);
-
-      mfem::out << GetTypeName((SCA_TYPE)t1) << "(" << orders[i] << ")"
-                << ":" << endl;
-
-      int c01 = conv_types[i];
-      string cmnt = "";
-
-      switch ((CONV_TYPE)c01)
-      {
-         case PROJECTION:
-            Projection(*x0, *y1);
-            break;
-         case INTERPOLATION_OP:
-            cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
-                   "(should match projection)" : "(not expected to succeed)";
-            InterpolationOp(*x0, *y1);
-            break;
-         case SOLVE:
-            LeastSquares((SCA_TYPE)t0, *x0, (SCA_TYPE)t1, *y1);
-            break;
-         default:
-            *y1 = 0.0;
-      }
-
-      {
-         ostringstream oss;
-         oss << field_names[i] << "_" << GetConvTypeShortName((CONV_TYPE)c01)
-             << "_" << GetTypeShortName((SCA_TYPE)t1) << "_o" << orders[i];
-
-         dc.RegisterField(oss.str(), y1);
-      }
-
-      if (visualization)
-      {
-         ostringstream oss;
-         oss << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
-             << GetTypeName((SCA_TYPE)t1)<< "(" << orders[i] << ")";
-         char vishost[] = "localhost";
-         int  visport   = 19916;
-         socketstream sol_sock1(vishost, visport);
-#ifdef MFEM_USE_MPI
-         sol_sock1 << "parallel " << mesh->GetNRanks() << ' '
-                   << mesh->GetMyRank() << '\n';
-#endif
-         sol_sock1.precision(8);
-         sol_sock1 << "solution\n" << *mesh << y1
-                   << "window_title '" << oss.str() << "'"
-                   << "window_geometry "
-                   << (int)((Ww + Fw) * (1.5 + c01) +
-                            Ws * (t0))
-                   << " " << (Wh + Fh) * (t1) + Ws * (t0)
-                   << " " << Ww << " " << Wh
-                   << flush;
-      }
-      mfem::out << endl;
-
-      // delete fes1;
-      // delete fec1;
-   }
-   dc.Save();
-
-   return 0;
-}
-
-FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return new H1_FECollection(p, dim);
-      case L2_TYPE:
-         return new L2_FECollection(p-1, dim);
-      case L2I_TYPE:
-         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
-                                    FiniteElement::INTEGRAL);
-      default:
-         return NULL;
-   }
-}
-
-ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
-                                   ParMesh &pmesh,
-                                   FiniteElementCollection &fec)
-{
-   return new ParFiniteElementSpace(&pmesh, &fec);
-}
-
-string GetTypeName(SCA_TYPE type)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return "    H1";
-      case L2_TYPE:
-         return "    L2";
-      case L2I_TYPE:
-         return "    L2I";
-      default:
-         return "--";
-   }
-}
-
-string GetTypeShortName(SCA_TYPE type)
-{
-   switch (type)
-   {
-      case H1_TYPE:
-         return "H1";
-      case L2_TYPE:
-         return "L2";
-      case L2I_TYPE:
-         return "L2I";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Projection            ";
-      case INTERPOLATION_OP:
-         return "Interpolation Operator";
-      case SOLVE:
-         return "Least Squares         ";
-      case SOLVE_W_DBC:
-         return "Least Squares with BC ";
-      default:
-         return "--";
-   }
-}
-
-string GetConvTypeShortName(CONV_TYPE type)
-{
-   switch (type)
-   {
-      case PROJECTION:
-         return "Proj";
-      case INTERPOLATION_OP:
-         return "Interp";
-      case SOLVE:
-         return "LS";
-      case SOLVE_W_DBC:
-         return "LSwBC";
-      default:
-         return "--";
-   }
-}
-
-void parseFieldNames(const char * field_name_c_str, vector<string> &field_names)
-{
-   string field_name_str(field_name_c_str);
-   string field_name;
-
-   for (string::iterator it=field_name_str.begin();
-        it!=field_name_str.end(); it++)
-   {
-      if (*it == '\\')
-      {
-         it++;
-         field_name.push_back(*it);
-      }
-      else if (*it == ' ')
-      {
-         if (!field_name.empty())
-         {
-            field_names.push_back(field_name);
-         }
-         field_name.clear();
-      }
-      else if (it == field_name_str.end() - 1)
-      {
-         field_name.push_back(*it);
-         field_names.push_back(field_name);
-      }
-      else
-      {
-         field_name.push_back(*it);
-      }
-   }
-   if (field_names.size() == 0)
-   {
-      field_names.push_back("ALL");
-   }
-}
-
-/** Perform a naive projection from one scalar field to another.
-
-   This scheme simply evaluates v0 at the interpolation points of v1.
-
-   If v0 has reduced continuity compared to v1 this can produce
-   results that depend on the order in which the elements are
-   traversed.
-
-   Suitable conversions:
-   H1 -> L2
-   H1 -> DG (same as L2)
- */
-void Projection(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   GridFunctionCoefficient v0Coef(&v0);
-   v1.ProjectCoefficient(v0Coef);
-}
-
-/** In theory this interpolation scheme should be equivalent to projection.
-
-    Building an interpolastion matrix could lead to computational
-    efficiency compared to simple projection if the operator will be
-    used several times.
-
-    Unfortunately this is broken for several combinations of source
-    and target fields.
-*/
-void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
-{
-   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
-   op.AddDomainInterpolator(new IdentityInterpolator);
-   op.Assemble();
-   op.Finalize();
-
-   op.Mult(v0, v1);
-}
-
-/** Compute a least-squares best fit using the target basis functions.
-
-    This scheme is more difficult to setup and more computationally
-    expensive but the results can be significantly better than simple
-    projections.
-*/
-void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
-                  SCA_TYPE t1, ParGridFunction &v1)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   HypreParMatrix * M = m.ParallelAssemble();
-
-   HypreDiagScale diag(*M);
-   HyprePCG pcg(*M);
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   Vector B, X;
-   b.ParallelAssemble(B);
-
-   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-
-   delete M;
-}
-
-/** Compute a least-squares best fit with boundary conditions.
-
-    This scheme is virtually identical to the previous one but it
-    makes use of boundary values, when available, to improve the
-    accuracy.  This scheme can produce significantly better results
-    when the normal derivative of the field is large near the
-    boundary.  This is particularly true when the field is
-    under-resolved near the boundary.
-*/
-void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
-                    SCA_TYPE t1, ParGridFunction &v1,
-                    Coefficient &c)
-{
-   ParFiniteElementSpace *fes0, *fes1;
-   fes0 = v0.ParFESpace();
-   fes1 = v1.ParFESpace();
-
-   ParMixedBilinearForm op(fes0, fes1);
-   op.AddDomainIntegrator(new MassIntegrator);
-   op.Assemble();
-   op.Finalize();
-
-   ParLinearForm b(v1.ParFESpace());
-   op.Mult(v0, b);
-
-   ParBilinearForm m(v1.ParFESpace());
-   m.AddDomainIntegrator(new MassIntegrator);
-   m.Assemble();
-   m.Finalize();
-
-   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
-   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
-      ess_bdr = 1;
-      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   if (t1 == H1_TYPE)
-   {
-      v1.ProjectBdrCoefficient(c, ess_bdr);
-   }
-
-   OperatorPtr M;
-   Vector B, X;
-   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
-
-   HypreDiagScale diag(*M.As<HypreParMatrix>());
-   HyprePCG pcg(*M.As<HypreParMatrix>());
-   pcg.SetPreconditioner(diag);
-   pcg.SetTol(1e-12);
-   pcg.SetMaxIter(1000);
-
-   pcg.Mult(B, X);
-   v1.Distribute(X);
-}