updting io benchmark

examples/CMakeLists.txt

@@ -64,6 +64,8 @@ if (MFEM_USE_MPI)
     ex25p.cpp
     ex26p.cpp
     ex27p.cpp
+    pa_oper.cpp
+    io_benchmark.cpp
     )
 endif()
 

examples/ex1p.cpp

@@ -70,6 +70,7 @@ int main(int argc, char *argv[])
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   int nfiles = 1;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -86,6 +87,7 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&nfiles, "-nf", "--num-files", "Number of files to write.");
    args.Parse();
    if (!args.Good())
    {
@@ -158,7 +160,7 @@ int main(int argc, char *argv[])
    {
       fec = new H1_FECollection(order = 1, dim);
    }
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec, 1, 0);
    HYPRE_Int size = fespace->GlobalTrueVSize();
    if (myid == 0)
    {
@@ -237,20 +239,52 @@ int main(int argc, char *argv[])
    //     local finite element solution on each processor.
    a->RecoverFEMSolution(X, *b, x);
 
+   std::string filename("nranks_");
+   filename += to_string(num_procs);
+   filename += ".gf";
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      x.Save(filename.c_str(), nfiles);
+      double t2 = MPI_Wtime();
+
+      if (myid == 0)
+      {
+         err << "elapsed write time: " << t2 - t1 << endl;
+      }
+   }
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      ParGridFunction new_x(fespace, filename.c_str());
+      double t2 = MPI_Wtime();
+      if (myid == 0)
+      {
+         err << "elapsed read time: " << t2 - t1 << endl;
+      }
+      // new_x -= x;
+      // out << "GF difference: " << new_x.Norml1() << endl;
+   }
+
    // 15. Save the refined mesh and the solution in parallel. This output can
    //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
    {
       ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
+      //mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << num_procs << setfill('0') << setw(6) << myid;
 
+      //ofstream mesh_ofs(mesh_name.str().c_str());
+      //mesh_ofs.precision(8);
+      //pmesh->Print(mesh_ofs);
+      double t1 = MPI_Wtime();
       ofstream sol_ofs(sol_name.str().c_str());
       sol_ofs.precision(8);
       x.Save(sol_ofs);
+      double t2 = MPI_Wtime();
+      if (myid == 0)
+      {
+         err << t2 - t1 << endl;
+      }
    }
 
    // 16. Send the solution by socket to a GLVis server.

examples/ex2p.cpp

@@ -274,6 +274,13 @@ int main(int argc, char *argv[])
       pmesh->SetNodalFESpace(fespace);
    }
 
+   {
+      x.Save("ex2p.gf", 1);
+      ParGridFunction new_x(fespace, "ex2p.gf");
+      new_x -= x;
+      out << "GF difference: " << new_x.Norml1() << endl;
+   }
+
    // 16. Save in parallel the displaced mesh and the inverted solution (which
    //     gives the backward displacements to the original grid). This output
    //     can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".

examples/io_benchmark.cpp

@@ -0,0 +1,292 @@
+//                       MFEM Example 1 - Parallel Version
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../data/star.mesh
+//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Device sample runs:
+//               mpirun -np 4 ex1p -pa -d cuda
+//               mpirun -np 4 ex1p -pa -d occa-cuda
+//               mpirun -np 4 ex1p -pa -d raja-omp
+//               mpirun -np 4 ex1p -pa -d ceed-cpu
+//               mpirun -np 4 ex1p -pa -d ceed-cuda
+//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mpi.h"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   // const char *mesh_file = "../data/star.mesh";
+   const char *mesh_file = "../data/square-disc.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = false;
+   int nfiles = 1;
+   // const char *out_file = "0_0.gf";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   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.AddOption(&nfiles, "-nf", "--num-files", "Number of files to write.");
+   // args.AddOption(&out_file, "-o", "--outfile",
+   //                "Name of file to write.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 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 (myid == 0) { 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);
+   int 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. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (pmesh->GetNodes())
+   {
+      fec = pmesh->GetNodes()->OwnFEC();
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec, 1, 0);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm *b = new ParLinearForm(fespace);
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x with initial guess of zero,
+   //     which satisfies the boundary conditions.
+   ParGridFunction x(fespace);
+   x = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the Diffusion
+   //     domain integrator.
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   // 13. Solve the linear system A X = B.
+   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
+   //     * With partial assembly, use Jacobi smoothing, for now.
+   Solver *prec = NULL;
+   if (pa)
+   {
+      if (UsesTensorBasis(*fespace))
+      {
+         prec = new OperatorJacobiSmoother(*a, ess_tdof_list);
+      }
+   }
+   else
+   {
+      prec = new HypreBoomerAMG;
+   }
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   if (prec) { cg.SetPreconditioner(*prec); }
+   cg.SetOperator(*A);
+   cg.Mult(B, X);
+   delete prec;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   std::string filename = to_string(num_procs) + "_" + to_string(nfiles) + "_";
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      x.Save(filename.c_str(), nfiles);
+      double t2 = MPI_Wtime();
+
+      double write_time = t2 - t1;
+      double average_write_time;
+      MPI_Reduce(&write_time, &average_write_time, 1,
+               MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+      if (myid == 0)
+      {
+         std::cout << "Average write time: " << average_write_time / num_procs << " for "
+                   << nfiles << " files and "  << num_procs << " ranks\n";
+      }
+   }
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      ParGridFunction temp_gf(fespace, filename.c_str());
+      double t2 = MPI_Wtime();
+
+      double read_time = t2 - t1;
+      double average_read_time;
+      MPI_Reduce(&read_time, &average_read_time, 1,
+                 MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+      if (myid == 0)
+      {
+         std::cout << "Average read time: " << average_read_time / num_procs << " for "
+                   << nfiles << " files and "  << num_procs << " ranks\n";
+      }
+   }
+
+   // 17. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   if (order > 0) { delete fec; }
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

examples/pa_oper.cpp

@@ -0,0 +1,906 @@
+//                                MFEM Example 9
+//
+// Compile with: make serial_nogpu
+//
+// Description:  This code solves the time-dependent advection-diffusion
+//               equation:
+//               \frac(\partial u}{\partial t}
+//                                 = \mathbf{a} \cdot \Nabla u - \nu \Nabla^2 u
+//               where a is a given advection velocity, \nu is the diffusion
+//               parameter, and u0(x) = u(0,x) is a given initial condition.
+//
+//               The demonstrates explicit time marching with H1 elements of
+//               arbitrary order. Periodic boundary conditions are used through
+//               periodic meshes. GLVis can be used for visualization of a
+//               time-evolving solution.
+
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+#include "mfem.hpp"
+#include "mpi.h"
+
+using namespace std;
+using namespace mfem;
+
+/** A time-dependent operator for the right-hand side of the ODE. The weak
+    form of du/dt = -a.grad(u) + nu Delta(u) is M du/dt = K u + b, where M and
+    K are the mass and advection-diffusion 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 AdvectionDiffusionEvolution : public mfem::TimeDependentOperator
+{
+public:
+   /// \param[in] M - bilinear form for mass matrix
+   /// \param[in] K - bilinear form for stiffness matrix
+   /// \param[in] b - load vector
+   AdvectionDiffusionEvolution(mfem::BilinearForm &M, mfem::BilinearForm &K,
+                               const mfem::Vector &b);
+
+   /// Perform the action of the operator: y = k = f(x, t), where k solves
+   /// Compute k = M^-1(Kx + l)
+   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Solve the implicit equation: k = f(x + dt k, t), for the unknown k at
+   /// the current time t. 
+   void ImplicitSolve(const double dt, const mfem::Vector &x,
+                      mfem::Vector &k) override;
+
+   virtual ~AdvectionDiffusionEvolution();
+
+private:
+   mfem::BilinearForm &M, &K;
+   const mfem::Vector &b;
+   /// solver for inverting mass matrix for explicit time-marching
+   std::unique_ptr<mfem::Solver> M_prec;
+   mfem::CGSolver M_solver;
+   /// solver for implicit time-marching
+   mfem::GSSmoother prec;
+   mfem::GMRESSolver linear_solver;
+   mfem::NewtonSolver newton;
+
+   mutable mfem::Vector z;
+
+   /// pointer-to-implementation idiom
+   /// Hides implementation details of this operator
+   class SystemOperator;
+   /// Operator that combines the linear spatial discretization with
+   /// the load vector into one operator used for implicit solves
+   std::unique_ptr<SystemOperator> combined_oper;
+
+   /// sets the state and dt for the combined operator
+   /// \param[in] dt - time increment
+   /// \param[in] x - the current state
+   void setOperParameters(double dt, const mfem::Vector *x);
+
+};
+
+class PAJacobianOperator : public mfem::Operator
+{
+public:
+   PAJacobianOperator(mfem::ParBilinearForm &_mass,
+                      mfem::ParBilinearForm &_stiff);
+
+   /// Compute r = J@k = M@k + dt*K@k
+   /// \param[in] k - dx/dt
+   /// \param[out] r - J@k = M@k + dt*K@k
+   void Mult(const mfem::Vector &k, mfem::Vector &r) const override;
+      
+   /// Set current dt values - needed to compute action of Jacobian.
+   void setParameters(double dt);
+
+private:
+   mfem::ParBilinearForm &mass;
+   mfem::ParBilinearForm &stiff;
+
+   double dt;
+};
+
+class ParSystemOperator : public mfem::Operator
+{
+public:
+   /// Nonlinear operator of the form that combines the mass, res, stiff,
+   /// and load elements for implicit/explicit ODE integration
+   /// \param[in] ess_bdr - array of boundaries attributes marked essential
+   /// \param[in] mass - bilinear form for mass matrix (not owned)
+   /// \param[in] res - nonlinear residual operator (not owned)
+   /// \param[in] stiff - bilinear form for stiffness matrix (not owned)
+   /// \param[in] load - load vector (not owned)
+   /// \param[in] a - used to move the spatial residual to the rhs
+   ParSystemOperator(mfem::ParBilinearForm &_mass,
+                     mfem::ParBilinearForm &_stiff);
+
+   /// Compute r = M@k + K@(x+dt*k)
+   /// (with `@` denoting matrix-vector multiplication)
+   /// \param[in] k - dx/dt 
+   /// \param[out] r - the residual
+   /// \note the signs on each operator must be accounted for elsewhere
+   void Mult(const mfem::Vector &k, mfem::Vector &r) const override;
+
+   /// Compute J = M + dt * K
+   /// \param[in] k - dx/dt 
+   mfem::Operator &GetGradient(const mfem::Vector &k) const override;
+
+   /// Set current dt and x values - needed to compute action and Jacobian.
+   void setParameters(double _dt, const mfem::Vector *_x);
+
+   ~ParSystemOperator();
+
+private:
+   mfem::ParBilinearForm &mass;
+   mfem::ParBilinearForm &stiff;
+   mutable mfem::HypreParMatrix *jacobian, *stiff_jacobian;
+
+   double dt;
+   const mfem::Vector *x;
+
+   mutable mfem::Vector work, work2;
+
+   std::unique_ptr<PAJacobianOperator> pa_jac;
+
+};
+
+/** A time-dependent operator for the right-hand side of the ODE. The weak
+    form of du/dt = -a.grad(u) + nu Delta(u) is M du/dt = K u + b, where M and
+    K are the mass and advection-diffusion 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 ParAdvectionDiffusionEvolution : public mfem::TimeDependentOperator
+{
+public:
+   /// \param[in] M - parallel bilinear form for mass matrix
+   /// \param[in] K - parallel bilinear form for stiffness matrix
+   ParAdvectionDiffusionEvolution(mfem::ParBilinearForm &M,
+                                  mfem::ParBilinearForm &K);
+
+   /// Perform the action of the operator: y = k = f(x, t), where k solves
+   /// Compute k = M^-1(Kx + l)
+   void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+   /// Solve the implicit equation: k = f(x + dt k, t), for the unknown k at
+   /// the current time t. 
+   void ImplicitSolve(const double dt, const mfem::Vector &x,
+                      mfem::Vector &k) override;
+
+   virtual ~ParAdvectionDiffusionEvolution();
+
+private:
+   mfem::OperatorHandle M_;
+   mfem::ParBilinearForm &M, &K;
+   /// solver for inverting mass matrix for explicit time-marching
+   std::unique_ptr<mfem::Solver> M_prec;
+   mfem::CGSolver M_solver;
+   /// solver for implicit time-marching
+   mfem::Solver *prec;
+   mfem::GMRESSolver linear_solver;
+   mfem::NewtonSolver newton;
+
+   mfem::Vector diag;
+   mutable mfem::Vector z, work, work2;
+
+   /// pointer-to-implementation idiom
+   /// Hides implementation details of this operator
+   /// Operator that combines the linear spatial discretization with
+   /// the load vector into one operator used for implicit solves
+   std::unique_ptr<ParSystemOperator> combined_oper;
+
+   /// sets the state and dt for the combined operator
+   /// \param[in] dt - time increment
+   /// \param[in] x - the current state
+   void setOperParameters(double dt, const mfem::Vector *x);
+
+};
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &X, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &X);
+
+// Inflow boundary condition
+double inflow_function(const Vector &X, const double t);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   problem = 3;
+   const char *mesh_file = "../data/periodic-square.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order = 3;
+   const char *device_config = "cpu";
+   int ode_solver_type = 22;
+   double t_final = 3 * 2*M_PI;
+   double dt = 0.01;
+   bool glvis = false;
+   bool paraview = false;
+   int vis_steps = 5;
+
+   double nu_val = 0.001;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   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(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&glvis, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&nu_val, "-nu", "--nu-value",
+                  "Value for \nu, the parameter that controls diffusion.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      std::cout << "Num ranks: " << num_procs << "\n";
+      args.PrintOptions(cout);
+   }
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 3. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 5. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter. If the mesh is of NURBS type, we convert it
+   //    to a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 7. Define the finite element space of the given
+   //    polynomial order on the refined mesh.
+   H1_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
+
+   HYPRE_Int global_vSize = fes->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << global_vSize << endl;
+   }
+
+   // 8. Set up and assemble the bilinear and linear forms corresponding to the
+   //    CG discretization.
+   /// negative to move the diffusion terms to the right side
+   ConstantCoefficient nu(-nu_val);
+   ConstantCoefficient one(1.0);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   ParBilinearForm *m_pa = new ParBilinearForm(fes);
+   ParBilinearForm *k_pa = new ParBilinearForm(fes);
+   m_pa->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   k_pa->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+   /// create mass matrix
+   m_pa->AddDomainIntegrator(new MassIntegrator(one));
+   /// add advection terms to stiffness matrix
+   k_pa->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   /// add diffusion terms to stiffness matrix
+   k_pa->AddDomainIntegrator(new DiffusionIntegrator(nu));
+
+   m_pa->Assemble();
+   int skip_zeros = 0;
+   k_pa->Assemble(skip_zeros);
+   m_pa->Finalize();
+   k_pa->Finalize(skip_zeros);
+
+   ParBilinearForm *m = new ParBilinearForm(fes);
+   ParBilinearForm *k = new ParBilinearForm(fes);
+   /// create mass matrix
+   m->AddDomainIntegrator(new MassIntegrator);
+   /// add advection terms to stiffness matrix
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   /// add diffusion terms to stiffness matrix
+   k->AddDomainIntegrator(new DiffusionIntegrator(nu));
+
+   m->Assemble();
+   k->Assemble(skip_zeros);
+   m->Finalize();
+   k->Finalize(skip_zeros);
+
+
+   ParGridFunction *u = new ParGridFunction(fes);
+   u->UseDevice(true);
+   u->ProjectCoefficient(u0);
+
+
+   HypreParVector *U = u->GetTrueDofs();
+
+   ParSystemOperator pso(*m, *k);
+   ParSystemOperator pso_pa(*m_pa, *k_pa);
+
+   pso.setParameters(dt, U);
+   pso_pa.setParameters(dt, U);
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_r(U->Size());
+   double t1 = MPI_Wtime();
+   pso.Mult(*U, pso_r);
+   double t2 = MPI_Wtime();
+   double fa_mult_time = t2 - t1;
+   double average_fa_mult_time;
+   MPI_Reduce(&fa_mult_time, &average_fa_mult_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "FA Mult time: " << average_fa_mult_time / num_procs << endl;
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_pa_r(U->Size());
+   double t3 = MPI_Wtime();
+   pso_pa.Mult(*U, pso_pa_r);
+   double t4 = MPI_Wtime();
+   double pa_mult_time = t4 - t3;
+   double average_pa_mult_time;
+   MPI_Reduce(&pa_mult_time, &average_pa_mult_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "FA Mult time: " << average_pa_mult_time / num_procs << endl;
+   
+   double local_mult_speedup = (t2-t1) / (t4-t3);
+   double global_mult_speedup;
+   MPI_Reduce(&local_mult_speedup, &global_mult_speedup, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+
+   if (myid == 0)
+      std::cout << "PA mult speedup: " << global_mult_speedup / num_procs << endl;
+
+   mfem::Vector diff_r(pso_pa_r);
+   diff_r -= pso_r;
+   // std::cout << "r diff: " << diff_r.Norml2() << std::endl;
+
+   mfem::Operator &pso_jac = pso.GetGradient(*U);
+   mfem::Operator &pso_pa_jac = pso_pa.GetGradient(*U);
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_jac_r(U->Size());
+   double t5 = MPI_Wtime();
+   pso_jac.Mult(*U, pso_jac_r);
+   double t6 = MPI_Wtime();
+   double fa_jac_mult_time = t6-t5;
+   double average_fa_jac_time;
+   MPI_Reduce(&fa_jac_mult_time, &average_fa_jac_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "FA Jac Mult time: " << average_fa_jac_time / num_procs << endl;
+
+   MPI_Barrier(MPI_COMM_WORLD);
+   mfem::Vector pso_pa_jac_r(U->Size());
+   double t7 = MPI_Wtime();
+   pso_pa_jac.Mult(*U, pso_pa_jac_r);
+   double t8 = MPI_Wtime();
+   double pa_jac_mult_time = t8-t7;
+   double average_pa_jac_time;
+   MPI_Reduce(&pa_jac_mult_time, &average_pa_jac_time, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "PA Jac Mult time: " << average_pa_jac_time / num_procs << endl;
+
+   double local_jac_speedup = (t6-t5) / (t8-t7);
+   double global_jac_speedup;
+   MPI_Reduce(&local_jac_speedup, &global_jac_speedup, 1,
+              MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+   if (myid == 0)
+      std::cout << "PA Jac mult speedup: " << global_jac_speedup / num_procs << endl;
+
+   // 13. Free the used memory.
+   delete U;
+   delete u;
+   delete k;
+   delete m;
+   delete fes;
+   delete pmesh;
+
+   MPI_Finalize();
+   return 0;
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 3:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 0:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x, const double t)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}
+
+class AdvectionDiffusionEvolution::SystemOperator : public mfem::Operator
+{
+public:
+   /// Nonlinear operator of the form that combines the mass, res, stiff,
+   /// and load elements for implicit/explicit ODE integration
+   /// \param[in] mass - bilinear form for mass matrix (not owned)
+   /// \param[in] res - nonlinear residual operator (not owned)
+   /// \param[in] stiff - bilinear form for stiffness matrix (not owned)
+   /// \param[in] load - load vector (not owned)
+   /// \param[in] a - used to move the spatial residual to the rhs
+   SystemOperator(BilinearForm &_mass, BilinearForm &_stiff,
+                  const mfem::Vector &b)
+      : Operator(_mass.Height()), mass(_mass), stiff(_stiff),
+        load(b), Jacobian(NULL), dt(0.0), x(NULL), work(height)
+      { }
+
+   /// Compute r = M@k + K@(x+dt*k) + l
+   /// (with `@` denoting matrix-vector multiplication)
+   /// \param[in] k - dx/dt 
+   /// \param[out] r - the residual
+   /// \note the signs on each operator must be accounted for elsewhere
+   void Mult(const mfem::Vector &k, mfem::Vector &r) const override
+   {
+      /// work = x+dt*k = x+dt*dx/dt = x+dx
+      add(1.0, *x, dt, k, work);
+      r = 0.0;
+      stiff.AddMult(work, r);
+      r += load;
+      mass.AddMult(k, r, -1.0);
+   }
+
+   /// Compute J = M + dt * K
+   /// \param[in] k - dx/dt 
+   mfem::Operator &GetGradient(const mfem::Vector &k) const override
+   {
+      delete Jacobian;
+      Jacobian = Add(-1.0, mass.SpMat(), dt, stiff.SpMat());
+      return *Jacobian;
+   }
+
+   /// Set current dt and x values - needed to compute action and Jacobian.
+   void setParameters(double _dt, const mfem::Vector *_x)
+   {
+      dt = _dt;
+      x = _x;
+   };
+
+   ~SystemOperator() {delete Jacobian;};
+
+private:
+   BilinearForm &mass;
+   BilinearForm &stiff;
+   const mfem::Vector &load;
+   mutable mfem::SparseMatrix *Jacobian;
+
+   double dt;
+   const mfem::Vector *x;
+
+   mutable mfem::Vector work, work2;
+
+};
+
+AdvectionDiffusionEvolution::AdvectionDiffusionEvolution(
+      BilinearForm &_M, BilinearForm &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()), M(_M), K(_K), b(_b),
+     z(_M.Height())
+{
+   bool pa = M.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   Array<int> ess_tdof_list;
+   if (pa)
+   {
+      M_prec.reset(new OperatorJacobiSmoother(M, ess_tdof_list));
+      M_solver.SetOperator(M);
+   }
+   else
+   {
+      M_prec.reset(new DSmoother(M.SpMat()));
+      M_solver.SetOperator(M.SpMat());
+   }
+
+   combined_oper.reset(new SystemOperator(_M, _K, _b));
+
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+
+   linear_solver.iterative_mode = true;
+   linear_solver.SetRelTol(1e-12);
+   linear_solver.SetAbsTol(0.0);
+   linear_solver.SetMaxIter(100);
+   linear_solver.SetPrintLevel(0);
+   linear_solver.SetPreconditioner(prec);
+
+   newton.iterative_mode = false;
+   newton.SetRelTol(1e-9);
+   newton.SetAbsTol(0.0);
+   newton.SetMaxIter(100);
+   newton.SetPrintLevel(-1);
+   newton.SetSolver(linear_solver);
+   newton.SetOperator(*combined_oper);
+}
+
+void AdvectionDiffusionEvolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void AdvectionDiffusionEvolution::ImplicitSolve(const double dt,
+                                                const Vector &x,
+                                                Vector &k)
+{
+   setOperParameters(dt, &x);
+   Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
+   newton.Mult(zero, k);
+   MFEM_VERIFY(newton.GetConverged(), "Newton solver did not converge!");
+}
+
+void AdvectionDiffusionEvolution::setOperParameters(double dt,
+                                                    const mfem::Vector *x)
+{
+   combined_oper->setParameters(dt, x);
+}
+
+AdvectionDiffusionEvolution::~AdvectionDiffusionEvolution() {}
+
+
+PAJacobianOperator::PAJacobianOperator(ParBilinearForm &_mass, ParBilinearForm &_stiff)
+ : Operator(_mass.ParFESpace()->GetTrueVSize()), mass(_mass), stiff(_stiff),
+   dt(0.0) { }
+
+
+void PAJacobianOperator::Mult(const mfem::Vector &k, mfem::Vector &r) const
+{
+   r.UseDevice(true);
+   r = 0.0;
+   stiff.TrueAddMult(k, r, dt);
+   mass.TrueAddMult(k, r, -1.0);
+}
+   
+void PAJacobianOperator::setParameters(const double _dt)
+{
+   dt = _dt;
+};
+
+ParSystemOperator::ParSystemOperator(ParBilinearForm &_mass, ParBilinearForm &_stiff)
+   : Operator(_mass.ParFESpace()->GetTrueVSize()), mass(_mass), stiff(_stiff),
+   jacobian(NULL), stiff_jacobian(NULL), dt(0.0), x(NULL),
+   work(height)
+{
+   pa_jac.reset(new PAJacobianOperator(mass, stiff));
+}
+
+/// Compute r = M@k + K@(x+dt*k)
+/// (with `@` denoting matrix-vector multiplication)
+/// \param[in] k - dx/dt 
+/// \param[out] r - the residual
+/// \note the signs on each operator must be accounted for elsewhere
+void ParSystemOperator::Mult(const mfem::Vector &k, mfem::Vector &r) const
+{
+   r = 0.0;
+   work.UseDevice(true);
+   work = 0.0;
+   /// work = x+dt*k = x+dt*dx/dt = x+dx
+   if (x)
+   {
+      add(1.0, *x, dt, k, work);
+   }
+
+   stiff.TrueAddMult(work, r);
+   mass.TrueAddMult(k, r, -1.0);
+}
+
+/// Compute J = M + dt * K
+/// \param[in] k - dx/dt 
+mfem::Operator &ParSystemOperator::GetGradient(const mfem::Vector &k) const
+{
+   bool mass_pa = mass.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   bool stiff_pa = stiff.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   
+   if (mass_pa && stiff_pa)
+   {
+      return *pa_jac.get();
+   }
+   else
+   {
+      delete stiff_jacobian;
+      delete jacobian;
+      jacobian = mass.ParallelAssemble();
+      *jacobian *= -1.0; //alpha;
+      stiff_jacobian = stiff.ParallelAssemble();
+      jacobian->Add(dt, *stiff_jacobian);
+      return *jacobian;
+   }
+}
+
+/// Set current dt and x values - needed to compute action and Jacobian.
+void ParSystemOperator::setParameters(const double _dt, const mfem::Vector *_x)
+{
+   dt = _dt;
+   x = _x;
+   pa_jac->setParameters(_dt);
+};
+
+ParSystemOperator::~ParSystemOperator()
+{
+   delete jacobian;
+   delete stiff_jacobian;
+};
+
+ParAdvectionDiffusionEvolution::ParAdvectionDiffusionEvolution(
+      ParBilinearForm &_M, ParBilinearForm &_K)
+   : TimeDependentOperator(_M.ParFESpace()->GetTrueVSize()), M(_M), K(_K), z(_M.Height())
+{
+   bool mass_pa = M.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+   bool stiff_pa = K.GetAssemblyLevel() == AssemblyLevel::PARTIAL;
+
+   Array<int> ess_tdof_list;
+   M_solver = CGSolver(MPI_COMM_WORLD);
+   if (mass_pa)
+   {
+      M_prec.reset(new OperatorJacobiSmoother(M, ess_tdof_list));
+      M_solver.SetOperator(M);
+   }
+   else
+   {
+      M_.Reset(_M.ParallelAssemble(), true);
+
+      // M_prec.reset(new HypreSmoother());
+      // M_solver.SetOperator(M.As<HypreParMatrix>());
+      HypreParMatrix &M_mat = *M_.As<HypreParMatrix>();
+      // HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
+      M_prec.reset(new HypreSmoother(M_mat, HypreSmoother::Jacobi));
+   }
+
+   combined_oper.reset(new ParSystemOperator(_M, _K));
+
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+
+   if (mass_pa && stiff_pa)
+   {
+      diag.UseDevice(true);
+      diag.SetSize(M.ParFESpace()->GetTrueVSize());
+      diag = 0.0;
+      work.UseDevice(true);
+      work2.UseDevice(true);
+      work.SetSize(M.ParFESpace()->GetTrueVSize());
+      work2.SetSize(M.ParFESpace()->GetTrueVSize());
+      work = 0.0;
+      work2 = 0.0;
+      M.AssembleDiagonal(work);
+
+      ParBilinearForm k(M.ParFESpace());
+      ConstantCoefficient nu(-0.01);
+      k.AddDomainIntegrator(new mfem::DiffusionIntegrator(nu));
+      k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      k.Assemble(0);
+      k.Finalize(0);
+      k.AssembleDiagonal(work2);
+
+      double dt = 0.1;
+      add(-1.0, work, dt, work2, diag);
+
+      prec = new OperatorChebyshevSmoother(combined_oper.get(), diag,
+                                           ess_tdof_list, 5,
+                                           M.ParFESpace()->GetComm());
+   }
+   else
+   {
+      prec = new HypreSmoother();
+   }
+   
+   linear_solver = GMRESSolver(MPI_COMM_WORLD);
+   linear_solver.iterative_mode = true;
+   linear_solver.SetRelTol(1e-12);
+   linear_solver.SetAbsTol(0.0);
+   linear_solver.SetMaxIter(2000);
+   linear_solver.SetPrintLevel(0);
+   linear_solver.SetPreconditioner(*prec);
+   linear_solver.SetKDim(2000);
+
+   newton.iterative_mode = true;
+   newton.SetRelTol(1e-9);
+   newton.SetAbsTol(0.0);
+   newton.SetMaxIter(10);
+   newton.SetPrintLevel(-1);
+   newton.SetSolver(linear_solver);
+   newton.SetOperator(*combined_oper);
+}
+
+void ParAdvectionDiffusionEvolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   M_solver.Mult(z, y);
+}
+
+void ParAdvectionDiffusionEvolution::ImplicitSolve(const double dt,
+                                                   const Vector &x,
+                                                   Vector &k)
+{
+   setOperParameters(dt, &x);
+   Vector zero; // empty vector is interpreted as zero r.h.s. by NewtonSolver
+   newton.Mult(zero, k);
+   MFEM_VERIFY(newton.GetConverged(), "Newton solver did not converge!");
+}
+
+void ParAdvectionDiffusionEvolution::setOperParameters(const double dt,
+                                                       const mfem::Vector *x)
+{
+   combined_oper->setParameters(dt, x);
+}
+
+ParAdvectionDiffusionEvolution::~ParAdvectionDiffusionEvolution() {delete prec;}

fem/pgridfunc.cpp

@@ -16,6 +16,7 @@
 #include "fem.hpp"
 #include <iostream>
 #include <limits>
+#include <string>
 #include "../general/forall.hpp"
 using namespace std;
 
@@ -78,6 +79,229 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
    fes = pfes;
 }
 
+ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf,
+                                 const char *_filename)
+   : GridFunction(pf), pfes(pf)
+{
+   MPI_Comm fes_comm;
+   int fes_rank, n_fes_ranks;
+   fes_comm = pfes->GetComm();
+   MPI_Comm_size(fes_comm, &n_fes_ranks);
+   MPI_Comm_rank(fes_comm, &fes_rank);
+
+   std::string filename(_filename);
+   std::string file_prefix;
+   std::string file_ext;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         file_prefix = (filename.substr(0, i));
+         file_ext = (filename.substr(i, filename.length() - i));
+      }
+   }
+
+   int nfiles = 1;
+   if (fes_rank == 0)
+   {
+      int n_rfes_ranks;
+      int tmp[2];
+      std::string mpi_filename;
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         mpi_filename = file_prefix + to_string(0) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + to_string(0);
+      }
+      
+      MPI_File fh;
+      MPI_File_open(MPI_COMM_SELF, mpi_filename.c_str(), MPI_MODE_RDONLY,
+                                                      MPI_INFO_NULL, &fh);
+      MPI_File_read_at(fh, 0, tmp, 2, MPI_INT, MPI_STATUS_IGNORE);
+      MPI_File_close(&fh);
+      
+      n_rfes_ranks = tmp[0];
+      nfiles = tmp[1];
+
+      MFEM_ASSERT(n_fes_ranks == n_rfes_ranks,
+         "ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf,"
+               " const char *_filename):\n"
+               "\tThe number of MPI ranks used to save the GridFunction is\n"
+               "\tnot the same as the number used to load it!");
+   }
+   MPI_Bcast(&nfiles, 1, MPI_INT, 0, fes_comm);
+
+   int color = fes_rank * nfiles / n_fes_ranks;
+   
+   MPI_Comm file_comm;
+   MPI_Comm_split(fes_comm, color, fes_rank, &file_comm);
+
+   int file_rank, n_file_ranks;
+   MPI_Comm_size(file_comm, &n_file_ranks);
+   MPI_Comm_rank(file_comm, &file_rank);
+
+   std::string mpi_filename;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos) {
+         mpi_filename = file_prefix + std::to_string(color) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + std::to_string(color);
+      }
+   }
+
+   MPI_File fh;
+   MPI_File_open(file_comm, mpi_filename.c_str(), MPI_MODE_RDONLY,
+                                                  MPI_INFO_NULL, &fh);
+
+   int *dof_counts = new int[5*n_file_ranks];
+   int **nv = new int*[n_file_ranks];
+   int **nvdofs = new int*[n_file_ranks];
+   int **nedofs = new int*[n_file_ranks];
+   int **nfdofs = new int*[n_file_ranks];
+   int **nrdofs = new int*[n_file_ranks];
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      nv[i]     = &dof_counts[i*5+0];
+      nvdofs[i] = &dof_counts[i*5+1];
+      nedofs[i] = &dof_counts[i*5+2];
+      nfdofs[i] = &dof_counts[i*5+3];
+      nrdofs[i] = &dof_counts[i*5+4];
+   }
+
+   *nv[file_rank] = pfes->GetVSize();
+   *nvdofs[file_rank] = pfes->GetNVDofs();
+   *nedofs[file_rank] = pfes->GetNEDofs();
+   *nfdofs[file_rank] = pfes->GetNFDofs();
+
+   int vdim = pfes->GetVDim();
+   *nrdofs[file_rank] = *nv[file_rank] / vdim - *nvdofs[file_rank] -
+                                       *nedofs[file_rank] - *nfdofs[file_rank];
+
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &dof_counts[0], 5,
+                 MPI_INT, file_comm);
+
+   double *data_ = HostWrite();
+
+   MPI_Offset header_offset = 0;
+   header_offset += 2 * sizeof(int);
+   MPI_Offset v_offset, e_offset, f_offset, r_offset;
+
+   int total_vdofs = 0, total_edofs = 0, total_fdofs = 0, total_rdofs = 0;
+   int total_scalar_dofs = 0;
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      total_vdofs += *nvdofs[i];
+      total_edofs += *nedofs[i];
+      total_fdofs += *nfdofs[i];
+      total_rdofs += *nrdofs[i];
+      total_scalar_dofs += *nv[i];
+   }
+
+   total_scalar_dofs /= vdim;
+
+   if (pfes->GetOrdering() == Ordering::byNODES)
+   {
+      for (int d = 0; d < vdim; ++d)
+      {
+         int v_data_offset = 0 + *nv[file_rank] * d / vdim ;
+         int e_data_offset = v_data_offset + *nvdofs[file_rank];
+         int f_data_offset = e_data_offset + *nedofs[file_rank];
+         int r_data_offset = f_data_offset + *nfdofs[file_rank];
+
+         v_offset = header_offset;
+         e_offset = header_offset;
+         f_offset = header_offset;
+         r_offset = header_offset;
+         
+         v_offset += total_scalar_dofs * d * sizeof(double);
+         e_offset += (total_vdofs + total_scalar_dofs * d) * sizeof(double);
+         f_offset += (total_vdofs + total_edofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+         r_offset += (total_vdofs + total_edofs + total_fdofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+
+
+         for (int i = 0; i < file_rank; ++i)
+         {
+            v_offset += *nvdofs[i] * sizeof(double);
+            e_offset += *nedofs[i] * sizeof(double);
+            f_offset += *nfdofs[i] * sizeof(double);
+            r_offset += *nrdofs[i] * sizeof(double);
+         }
+
+         MPI_File_read_at_all(fh, v_offset, &data_[v_data_offset],
+                              *nvdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, e_offset, &data_[e_data_offset],
+                              *nedofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, f_offset, &data_[f_data_offset],
+                              *nfdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, r_offset, &data_[r_data_offset],
+                              *nrdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+      }
+   }
+   else
+   {
+      v_offset = header_offset;
+      e_offset = v_offset + total_vdofs * vdim * sizeof(double);
+      f_offset = e_offset + total_edofs * vdim * sizeof(double);
+      r_offset = f_offset + total_fdofs * vdim * sizeof(double);
+
+      for (int i = 0; i < file_rank; ++i)
+      {
+         v_offset += *nvdofs[i] * sizeof(double) * vdim;
+         e_offset += *nedofs[i] * sizeof(double) * vdim;
+         f_offset += *nfdofs[i] * sizeof(double) * vdim;
+         r_offset += *nrdofs[i] * sizeof(double) * vdim;
+      }
+
+      int v_data_offset = 0;
+      int e_data_offset = v_data_offset + *nvdofs[file_rank] * vdim;
+      int f_data_offset = e_data_offset + *nedofs[file_rank] * vdim;
+      int r_data_offset = f_data_offset + *nfdofs[file_rank] * vdim;
+      
+      MPI_File_read_at_all(fh, v_offset, &data_[v_data_offset],
+                           *nvdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, e_offset, &data_[e_data_offset],
+                           *nedofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, f_offset, &data_[f_data_offset],
+                           *nfdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, r_offset, &data_[r_data_offset],
+                           *nrdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+   }
+
+   MPI_File_close(&fh);
+   MPI_Comm_free(&file_comm);
+
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   delete[] dof_counts;
+   delete[] nv;
+   delete[] nvdofs;
+   delete[] nedofs;
+   delete[] nfdofs;
+   delete[] nrdofs;
+}
+
+
 void ParGridFunction::Update()
 {
    face_nbr_data.Destroy();
@@ -518,6 +742,201 @@ void ParGridFunction::Save(adios2stream &out,
 }
 #endif
 
+void ParGridFunction::Save(const char *_filename, const int nfiles)
+{
+   MPI_Comm fes_comm;
+   int fes_rank, n_fes_ranks;
+   fes_comm = pfes->GetComm();
+
+   MPI_Comm_size(fes_comm, &n_fes_ranks);
+   MPI_Comm_rank(fes_comm, &fes_rank);
+
+   int color = fes_rank * nfiles / n_fes_ranks;
+   
+   MPI_Comm file_comm;
+   MPI_Comm_split(fes_comm, color, fes_rank, &file_comm);
+
+   int file_rank, n_file_ranks;
+   MPI_Comm_size(file_comm, &n_file_ranks);
+   MPI_Comm_rank(file_comm, &file_rank);
+
+   std::string filename(_filename);
+   std::string file_prefix;
+   std::string file_ext;
+   std::string mpi_filename;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         file_prefix = (filename.substr(0, i));
+         file_ext = (filename.substr(i, filename.length() - i));
+         mpi_filename = file_prefix + std::to_string(color) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + std::to_string(color);
+      }
+   }
+
+   MPI_File fh;
+   MPI_File_open(file_comm, mpi_filename.c_str(), MPI_MODE_CREATE |
+                                                  MPI_MODE_WRONLY,
+                                                  MPI_INFO_NULL, &fh);
+
+   int *dof_counts = new int[5*n_file_ranks];
+   int **nv = new int*[n_file_ranks];
+   int **nvdofs = new int*[n_file_ranks];
+   int **nedofs = new int*[n_file_ranks];
+   int **nfdofs = new int*[n_file_ranks];
+   int **nrdofs = new int*[n_file_ranks];
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      nv[i]     = &dof_counts[i*5+0];
+      nvdofs[i] = &dof_counts[i*5+1];
+      nedofs[i] = &dof_counts[i*5+2];
+      nfdofs[i] = &dof_counts[i*5+3];
+      nrdofs[i] = &dof_counts[i*5+4];
+   }
+
+   *nv[file_rank] = pfes->GetVSize();
+   *nvdofs[file_rank] = pfes->GetNVDofs();
+   *nedofs[file_rank] = pfes->GetNEDofs();
+   *nfdofs[file_rank] = pfes->GetNFDofs();
+
+   int vdim = pfes->GetVDim();
+   *nrdofs[file_rank] = *nv[file_rank] / vdim - *nvdofs[file_rank] -
+                                       *nedofs[file_rank] - *nfdofs[file_rank];
+
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &dof_counts[0], 5,
+                 MPI_INT, file_comm);
+
+   double *data_  = const_cast<double*>(HostRead());
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   MPI_Offset header_offset = 0;
+
+   if (file_rank == 0)
+   {
+      int tmp[] = {n_fes_ranks, nfiles};
+      MPI_File_write_at(fh, header_offset, &tmp, 2, MPI_INT,
+                        MPI_STATUS_IGNORE);
+   }
+
+   header_offset += 2 * sizeof(int);
+
+   MPI_Offset v_offset, e_offset, f_offset, r_offset;
+
+   int total_vdofs = 0, total_edofs = 0, total_fdofs = 0, total_rdofs = 0;
+   int total_scalar_dofs = 0;
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      total_vdofs += *nvdofs[i];
+      total_edofs += *nedofs[i];
+      total_fdofs += *nfdofs[i];
+      total_rdofs += *nrdofs[i];
+      total_scalar_dofs += *nv[i];
+   }
+
+   total_scalar_dofs /= vdim;
+
+   if (pfes->GetOrdering() == Ordering::byNODES)
+   {
+      for (int d = 0; d < vdim; ++d)
+      {
+         int v_data_offset = 0 + *nv[file_rank] * d / vdim ;
+         int e_data_offset = v_data_offset + *nvdofs[file_rank];
+         int f_data_offset = e_data_offset + *nedofs[file_rank];
+         int r_data_offset = f_data_offset + *nfdofs[file_rank];
+
+         v_offset = header_offset;
+         e_offset = header_offset;
+         f_offset = header_offset;
+         r_offset = header_offset;
+         
+         v_offset += total_scalar_dofs * d * sizeof(double);
+         e_offset += (total_vdofs + total_scalar_dofs * d) * sizeof(double);
+         f_offset += (total_vdofs + total_edofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+         r_offset += (total_vdofs + total_edofs + total_fdofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+
+         for (int i = 0; i < file_rank; ++i)
+         {
+            v_offset += *nvdofs[i] * sizeof(double);
+            e_offset += *nedofs[i] * sizeof(double);
+            f_offset += *nfdofs[i] * sizeof(double);
+            r_offset += *nrdofs[i] * sizeof(double);
+         }
+
+         MPI_File_write_at_all(fh, v_offset, &data_[v_data_offset],
+                               *nvdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, e_offset, &data_[e_data_offset],
+                               *nedofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, f_offset, &data_[f_data_offset],
+                               *nfdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, r_offset, &data_[r_data_offset],
+                               *nrdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+      }
+   }
+   else
+   {
+      v_offset = header_offset;
+      e_offset = v_offset + total_vdofs * vdim * sizeof(double);
+      f_offset = e_offset + total_edofs * vdim * sizeof(double);
+      r_offset = f_offset + total_fdofs * vdim * sizeof(double);
+
+      for (int i = 0; i < file_rank; ++i)
+      {
+         v_offset += *nvdofs[i] * sizeof(double) * vdim;
+         e_offset += *nedofs[i] * sizeof(double) * vdim;
+         f_offset += *nfdofs[i] * sizeof(double) * vdim;
+         r_offset += *nrdofs[i] * sizeof(double) * vdim;
+      }
+
+      int v_data_offset = 0;
+      int e_data_offset = v_data_offset + *nvdofs[file_rank] * vdim;
+      int f_data_offset = e_data_offset + *nedofs[file_rank] * vdim;
+      int r_data_offset = f_data_offset + *nfdofs[file_rank] * vdim;
+      
+      MPI_File_write_at_all(fh, v_offset, &data_[v_data_offset],
+                            *nvdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, e_offset, &data_[e_data_offset],
+                            *nedofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, f_offset, &data_[f_data_offset],
+                            *nfdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, r_offset, &data_[r_data_offset],
+                            *nrdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+   }
+
+   MPI_File_close(&fh);
+   MPI_Comm_free(&file_comm);
+   
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   delete[] dof_counts;
+   delete[] nv;
+   delete[] nvdofs;
+   delete[] nedofs;
+   delete[] nfdofs;
+   delete[] nrdofs;
+}
+
 void ParGridFunction::SaveAsOne(std::ostream &out)
 {
    int i, p;

fem/pgridfunc.hpp

@@ -83,6 +83,13 @@ public:
        constructed. The new ParGridFunction assumes ownership of both. */
    ParGridFunction(ParMesh *pmesh, std::istream &input);
 
+   /// Construct a ParGridFunction by loading a ParGridFunction saved using
+   /// ParGridFunction::Save(char *filename, int nfiles).
+   /** The parallel space @a *pf and the space used by the GridFunction saved
+       in @a *filename should match. The number of ranks used when loading the
+       ParGridFunction must be the same as when it was saved. */
+   ParGridFunction(ParFiniteElementSpace *pf, const char *filename);
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use ParFiniteElementSpace%s
        that have the same size.
@@ -324,6 +331,20 @@ public:
       const adios2stream::data_type type = adios2stream::data_type::point_data) const;
 #endif
 
+   /** Save the local grid functions to n number of files, where each file will
+       contain the grid functions from potentially multiple ranks. This is
+       similar to the syncIO approach from "Fu, Jing, et al. 'Scalable parallel
+       I/O alternatives for massively parallel partitioned solver systems.'
+       2010 IEEE International Symposium on Parallel & Distributed Processing,
+       Workshops and Phd Forum (IPDPSW). IEEE, 2010." 
+       @param[in] filename - filename for output files with extension
+       @param[in] nfiles - number of files to write using MPI-IO
+       @note - takes into account the signs of the local dofs.
+       @note - writes a binary file without the FESpace header; the saved file
+       should only be loaded by the accompanying constructor:
+       ParGridFunction(ParFiniteElementSpace *pf, const char *filename) */
+   void Save(const char *filename, const int nfiles = 1);
+
    /// Merge the local grid functions
    void SaveAsOne(std::ostream &out = mfem::out);