Squashing commits on AIR

.gitignore

@@ -6,6 +6,8 @@
 # Also ignore OS-specific files like .DS_Store on Mac
 # ------------------------------------------------------------------------------
 
+*DS_Store
+
 # Object and library files
 *.o
 /libmfem.*
@@ -24,6 +26,7 @@ CMakeFiles/
 /deps.mk
 config/_config.hpp
 config/config.mk
+config/user.cmake
 config/sample-runs-build.log
 doc/CodeDocumentation.conf
 doc/CodeDocumentation.html

config/defaults.mk

@@ -103,9 +103,9 @@ MFEM_MPI_NP = 4
 # config.hpp. The values below are the defaults for generating the actual values
 # in config.mk and config.hpp.
 
-MFEM_USE_MPI           = NO
+MFEM_USE_MPI           = YES
 MFEM_USE_METIS         = $(MFEM_USE_MPI)
-MFEM_USE_METIS_5       = NO
+MFEM_USE_METIS_5       = YES
 MFEM_DEBUG             = NO
 MFEM_USE_EXCEPTIONS    = NO
 MFEM_USE_GZSTREAM      = NO

examples/CMakeLists.txt

@@ -28,6 +28,7 @@ list(APPEND ALL_EXE_SRCS
   ex19.cpp
   ex20.cpp
   ex21.cpp
+  ex23.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -53,6 +54,7 @@ if (MFEM_USE_MPI)
     ex19p.cpp
     ex20p.cpp
     ex21p.cpp
+    ex23p.cpp
     )
 endif()
 

examples/ex23.cpp

@@ -0,0 +1,553 @@
+//                               MFEM Example 23
+//
+// Compile with: make ex23
+//
+// Sample runs:
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -r 2 -dt 0.005
+//    ex23 -m ../data/periodic-square.mesh -p 0 -r 2 -dt 0.01
+//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.01
+//    ex23 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex23 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex23 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
+//    ex23 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex23 -m ../data/star-mixed.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9 -d 0.05
+//    ex23 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
+//    ex23 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation
+//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
+//               D is a diffusion coefficient,
+//               v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of implicit
+//               ODE time integrators, the definition of periodic boundary
+//               conditions through periodic meshes, as well as the use of GLVis
+//               for persistent visualization of a time-evolving solution. The
+//               saving of time-dependent data files for external visualization
+//               with VisIt (visit.llnl.gov) is also illustrated.
+//
+//               This example is a merger of examples 9 and 14.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// 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);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the implicit or explicit solve for du/dt. */
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   SparseMatrix *A;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   DSmoother *A_prec;
+   GMRESSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   FE_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~FE_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   int ode_solver_type = 3;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   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(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
+                  "\t   11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // Implicit L-stable methods
+      case 1:  ode_solver = new BackwardEulerSolver; break;
+      case 2:  ode_solver = new SDIRK23Solver(2); break;
+      case 3:  ode_solver = new SDIRK33Solver; break;
+      // Explicit methods
+      case 11: ode_solver = new ForwardEulerSolver; break;
+      case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 13: ode_solver = new RK3SSPSolver; break;
+      case 14: ode_solver = new RK4Solver; break;
+      case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. 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 < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the parallel discontinuous DG finite element space on the
+   //    parallel refined mesh of the given polynomial order.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient diff_coef(d_coef);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm s(&fes);
+   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
+                                                         kappa));
+   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   BilinearForm k(&fes);
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k.AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k.AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m.Assemble(skip_zeros);
+   m.Finalize(skip_zeros);
+   s.Assemble(skip_zeros);
+   s.Finalize(skip_zeros);
+   k.Assemble(skip_zeros);
+   k.Finalize(skip_zeros);
+   b.Assemble();
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   {
+      ofstream omesh("ex23.mesh");
+      omesh.precision(precision);
+      mesh.Print(omesh);
+      ofstream osol("ex23-init.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example23", &mesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example23", &mesh);
+         dc->SetPrecision(precision);
+      }
+      dc->RegisterField("solution", &u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   FE_Evolution adv(m.SpMat(), s.SpMat(), k.SpMat(), b);
+
+   double t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      double dt_real = min(dt, t_final - t);
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+
+         if (visualization)
+         {
+            sout << "solution\n" << mesh << u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   // 9. Save the final solution in parallel. This output can be viewed later
+   //    using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
+   {
+      ofstream osol("ex23-final.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete dc;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                           SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     M_prec(M),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      SparseMatrix * SK = Add(1.0, S, -1.0, K);
+      A = Add(1.0, M, _dt, *SK);
+      delete SK;
+      dt = _dt;
+
+      A_prec = new DSmoother(*A);
+      A_solver = new GMRESSolver;
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void FE_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+
+   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// 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 0:
+      {
+         // 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 3:
+      {
+         // 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)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/ex23p.cpp

@@ -0,0 +1,673 @@
+//                       MFEM Example 23 - Parallel Version
+//
+// Compile with: make ex23p
+//
+// Sample runs:
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -dt 0.005
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 1 -dt 0.005 -tf 9
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 1 -dt 0.005 -tf 9
+//    mpirun -np 4 ex23p -m ../data/amr-quad.mesh -p 1 -rp 1 -dt 0.002 -tf 9
+//    mpirun -np 4 ex23p -m ../data/star-q3.mesh -p 1 -rp 1 -dt 0.004 -tf 9
+//    mpirun -np 4 ex23p -m ../data/star-mixed.mesh -p 1 -rp 1 -dt 0.004 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 1 -rp 1 -dt 0.005 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 2 -rp 1 -dt 0.005 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 3 -rp 1 -dt 0.005 -tf 9 -d 0.05
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 3 -rp 2 -dt 0.0025 -tf 9 -vs 20
+//    mpirun -np 4 ex23p -m ../data/periodic-cube.mesh -p 0 -o 2 -rp 1 -dt 0.01 -tf 8
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation
+//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
+//               D is a diffusion coefficient,
+//               v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of implicit
+//               ODE time integrators, the definition of periodic boundary
+//               conditions through periodic meshes, as well as the use of GLVis
+//               for persistent visualization of a time-evolving solution. The
+//               saving of time-dependent data files for external visualization
+//               with VisIt (visit.llnl.gov) is also illustrated.
+//
+//               This example is a merger of examples 9 and 14.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem, use_gmres;
+bool use_AIR;
+
+// 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);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+struct AIR_parameters {
+   double distanceR;
+   std::string prerelax;
+   std::string postrelax;
+   int interp_type;
+   int relax_type;
+   int coarsen_type;
+   double strength_tolC;
+   double strength_tolR;
+   double filter_tolR;
+   double filterA_tol;
+};
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the implicit or explicit solve for du/dt. */
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   HypreParMatrix *A;
+   HypreParMatrix A_s;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   HypreBoomerAMG *AMG_solver;
+   HypreGMRES *GMRES_solver;
+   double dt;
+   int blocksize;
+
+   mutable Vector z;
+
+public:
+   FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                const Vector &_b, int order);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~FE_Evolution() { delete GMRES_solver; delete AMG_solver; delete A; }
+};
+
+
+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 = 0;
+   use_gmres = true;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 0;
+   int order = 3;
+   int ode_solver_type = 3;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   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(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
+                  "\t   11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 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();
+
+   // 4. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // Implicit L-stable methods
+      case 1:  ode_solver = new BackwardEulerSolver; break;
+      case 2:  ode_solver = new SDIRK23Solver(2); break;
+      case 3:  ode_solver = new SDIRK33Solver; break;
+      // Explicit methods
+      case 11: ode_solver = new ForwardEulerSolver; break;
+      case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 13: ode_solver = new RK3SSPSolver; break;
+      case 14: ode_solver = new RK4Solver; break;
+      case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         delete mesh;
+         return 3;
+   }
+
+   // 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();
+   }
+   if (mesh->NURBSext)
+   {
+      mesh->SetCurvature(max(order, 1));
+   }
+   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();
+   }
+
+   // Get mesh size, compare with time step and diffusion coefficient. Use
+   // classical AMG for diffusion-dominated problems, and nonsymmetric AMG
+   // based on approximate ideal restriction (AIR) for advection dominated.
+   double h_min, h_max, k_min, k_max;
+   pmesh->GetCharacteristics(h_min, h_max, k_min, k_max);
+   if (dt > d_coef*h_max) use_AIR = true;
+   else use_AIR = true;
+   cout << "ratio = " << d_coef*h_max / dt << "\n";
+
+   // 7. Define the parallel discontinuous DG finite element space on the
+   //    parallel refined mesh of the given polynomial order.
+   DG_FECollection fec(order, dim);
+   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 parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient diff_coef(d_coef);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   ParBilinearForm *m = new ParBilinearForm(fes);
+   m->AddDomainIntegrator(new MassIntegrator);
+
+   ParBilinearForm *s = new ParBilinearForm(fes);
+   s->AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
+                                                          kappa));
+   s->AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   ParBilinearForm *k = new ParBilinearForm(fes);
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k->AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k->AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+
+   ParLinearForm *b = new ParLinearForm(fes);
+   b->AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m->Assemble(skip_zeros);
+   m->Finalize(skip_zeros);
+   s->Assemble(skip_zeros);
+   s->Finalize(skip_zeros);
+   k->Assemble(skip_zeros);
+   k->Finalize(skip_zeros);
+   b->Assemble();
+
+   HypreParMatrix *M = m->ParallelAssemble();
+   HypreParMatrix *S = s->ParallelAssemble();
+   HypreParMatrix *K = k->ParallelAssemble();
+   HypreParVector *B = b->ParallelAssemble();
+
+   // 9. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   ParGridFunction *u = new ParGridFunction(fes);
+   u->ProjectCoefficient(u0);
+   HypreParVector *U = u->GetTrueDofs();
+
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "ex23-mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex23-init." << setfill('0') << setw(6) << myid;
+      ofstream omesh(mesh_name.str().c_str());
+      omesh.precision(precision);
+      pmesh->Print(omesh);
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example23-Parallel", pmesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example23-Parallel", pmesh);
+         dc->SetPrecision(precision);
+         // To save the mesh using MFEM's parallel mesh format:
+         // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
+      }
+      dc->RegisterField("solution", u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         if (myid == 0)
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+         visualization = false;
+         if (myid == 0)
+         {
+            cout << "GLVis visualization disabled.\n";
+         }
+      }
+      else
+      {
+         sout << "parallel " << num_procs << " " << myid << "\n";
+         sout.precision(precision);
+         sout << "solution\n" << *pmesh << *u;
+         sout << "pause\n";
+         sout << flush;
+         if (myid == 0)
+            cout << "GLVis visualization paused."
+                 << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 10. Define the time-dependent evolution operator describing the ODE
+   //     right-hand side, and perform time-integration (looping over the time
+   //     iterations, ti, with a time-step dt).
+   FE_Evolution adv(*M, *S, *K, *B, order);
+
+   double t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      double dt_real = min(dt, t_final - t);
+      ode_solver->Step(*U, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         if (myid == 0)
+         {
+            cout << "time step: " << ti << ", time: " << t << endl;
+         }
+
+         // 11. Extract the parallel grid function corresponding to the finite
+         //     element approximation U (the local solution on each processor).
+         *u = *U;
+
+         if (visualization)
+         {
+            sout << "parallel " << num_procs << " " << myid << "\n";
+            sout << "solution\n" << *pmesh << *u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   // 12. Save the final solution in parallel. This output can be viewed later
+   //     using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
+   {
+      *u = *U;
+      ostringstream sol_name;
+      sol_name << "ex23-final." << setfill('0') << setw(6) << myid;
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // 13. Free the used memory.
+   delete U;
+   delete u;
+   delete B;
+   delete b;
+   delete K;
+   delete k;
+   delete S;
+   delete s;
+   delete M;
+   delete m;
+   delete fes;
+   delete pmesh;
+   delete ode_solver;
+   delete dc;
+
+   MPI_Finalize();
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                           HypreParMatrix &_K, const Vector &_b, int order)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), b(_b), GMRES_solver(NULL), AMG_solver(NULL),
+     M_prec(M), M_solver(M.GetComm()), A(NULL),
+     dt(-1.0), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+
+   // DG block size given by (FEorder+1)^2 on square meshes.
+   blocksize = (order+1)*(order+1);
+}
+
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+
+void FE_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   if ((fabs(dt - _dt) > 1e-4 * _dt) || !A)
+   {
+      delete GMRES_solver;
+      delete AMG_solver;
+      delete A;
+
+      dt = _dt;
+      HypreParMatrix *SK = Add(1.0, S, -1.0, K);
+      A = Add(1.0, M, dt, *SK);
+      delete SK;
+
+      BlockInvScal(A, &A_s, NULL, NULL, blocksize, 0);
+
+      int print_level = 1;
+      AMG_solver = new HypreBoomerAMG(A_s);
+      AMG_solver->SetMaxLevels(50);      
+      if (use_AIR) {
+         AMG_solver->SetLAIROptions(1.5, "", "FFC", 0.1, 0.01, 0.0,
+                                    100, 3, 0.0, 10, -1, 1);
+                                // 100, 3, 0.0, 6, -1, 1);
+      }
+      else {
+         AMG_solver->SetInterpolation(0);
+         AMG_solver->SetCoarsening(6);
+         AMG_solver->SetAggressiveCoarsening(1);
+      }
+
+      if (use_gmres) {
+         GMRES_solver = new HypreGMRES(A_s);
+         GMRES_solver->SetTol(1e-12);
+         GMRES_solver->SetMaxIter(100);
+         GMRES_solver->SetPrintLevel(print_level);
+         GMRES_solver->SetPreconditioner(*AMG_solver);
+         GMRES_solver->iterative_mode = false;
+      }
+      else {
+         AMG_solver->SetPrintLevel(print_level);
+         AMG_solver->SetTol(1e-12);
+         AMG_solver->SetMaxIter(100);
+      }
+   }
+
+   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+
+   // Scale block inverse to right hand side
+   HypreParVector b_s;
+   BlockInvScal(A, NULL, &z, &b_s, blocksize, 2);
+   GMRES_solver->Mult(b_s, y);
+
+}
+
+// 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 0:
+      {
+         // 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 3:
+      {
+         // 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)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/ex9p.cpp

@@ -37,7 +37,7 @@ using namespace mfem;
 
 // Choice for the problem setup. The fluid velocity, initial condition and
 // inflow boundary condition are chosen based on this parameter.
-int problem;
+int problem, trisolve, use_gmres;
 
 // Velocity coefficient
 void velocity_function(const Vector &x, Vector &v);
@@ -51,6 +51,19 @@ double inflow_function(const Vector &x);
 // Mesh bounding box
 Vector bb_min, bb_max;
 
+struct AIR_parameters {
+   double distanceR;
+   std::string prerelax;
+   std::string postrelax;
+   int interp_type;
+   int relax_type;
+   int coarsen_type;
+   double strength_tolC;
+   double strength_tolR;
+   double filter_tolR;
+   double filterA_tol;
+};
+
 
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
@@ -60,19 +73,33 @@ Vector bb_min, bb_max;
 class FE_Evolution : public TimeDependentOperator
 {
 private:
-   HypreParMatrix &M, &K;
+   HypreParMatrix &M, &K, *A, A_s;
    const Vector &b;
    HypreSmoother M_prec;
    CGSolver M_solver;
 
+   // Preconditioner/solvers for A
+   HypreBoomerAMG *AMG_solver;
+   HypreGMRES     *GMRES_solver;
+   HypreTriSolve *preconditioner;
+   AIR_parameters &AIR;
+
+   double dt;
+   int blocksize;
+
    mutable Vector z;
 
 public:
-   FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_K, const Vector &_b);
+   FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_K, const Vector &_b,
+                int order, AIR_parameters &_AIR);
 
+   /** Solve the Backward-Euler equation: d = f(x + dt*d, t+dt), where u_t = f(x,t).
+       This is the only requirement for high-order SDIRK implicit integration.*/
+   virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
    virtual void Mult(const Vector &x, Vector &y) const;
 
-   virtual ~FE_Evolution() { }
+   virtual ~FE_Evolution();
+
 };
 
 
@@ -86,21 +113,28 @@ int main(int argc, char *argv[])
 
    // 2. Parse command-line options.
    problem = 0;
+   use_gmres = 0;
+   trisolve = 0;
    const char *mesh_file = "../data/periodic-hexagon.mesh";
    int ser_ref_levels = 2;
    int par_ref_levels = 0;
    int order = 3;
-   int ode_solver_type = 4;
+   int ode_solver_type = 3;
    double t_final = 10.0;
    double dt = 0.01;
    bool visualization = true;
    bool visit = false;
    bool binary = false;
    int vis_steps = 5;
+   int basis_type = 1;
 
    int precision = 8;
    cout.precision(precision);
 
+    AIR_parameters AIR = {1, "", "FA", 100, 10, 10, 0.1, 0.01, 0.0, 1e-4};
+    const char* temp_prerelax = "";
+    const char* temp_postrelax = "FA";
+
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
                   "Mesh file to use.");
@@ -119,6 +153,32 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption(&basis_type, "-b", "--basis-type",
+                  "DG finite element basis type. 0 for G-Leg, 1 for G-Lob.");
+   args.AddOption(&use_gmres, "-gmres", "--use-gmres",
+                  "Boolean to use GMRES as solver (default with AIR preconditioning).");
+   args.AddOption(&trisolve, "-trisolve", "--precond-trisolve",
+                  "Precondition GMRES with an on-processor triangular solve.");
+   args.AddOption(&(AIR.distanceR), "-Ad", "--AIR-distance",
+                  "Distance restriction neighborhood for AIR.");
+   args.AddOption(&(AIR.interp_type), "-Ai", "--AIR-interpolation",
+                  "Index for hypre interpolation routine.");
+   args.AddOption(&(AIR.coarsen_type), "-Ac", "--AIR-coarsen_type",
+                  "Index for hypre coarsening routine.");
+   args.AddOption(&(AIR.strength_tolC), "-AsC", "--AIR-strengthC",
+                  "Theta value determining strong connections for AIR (coarsen_type).");
+   args.AddOption(&(AIR.strength_tolR), "-AsR", "--AIR-strengthR",
+                  "Theta value determining strong connections for AIR (restriction).");
+   args.AddOption(&(AIR.filter_tolR), "-AfR", "--AIR-filterR",
+                  "Theta value eliminating small entries in restriction (after building).");
+   args.AddOption(&(AIR.filterA_tol), "-Af", "--AIR-filter",
+                  "Theta value to eliminate small connections in AIR hierarchy. Use -1 to specify O(h).");
+   args.AddOption(&(AIR.relax_type), "-Ar", "--AIR-relaxation",
+                  "Index for hypre relaxation routine.");
+   args.AddOption(&temp_prerelax, "-Ar1", "--AIR-prerelax",
+                  "String denoting prerelaxation scheme; e.g., FCC.");
+   args.AddOption(&temp_postrelax, "-Ar2", "--AIR-postrelax",
+                  "String denoting postrelaxation scheme; e.g., FFC.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -131,6 +191,9 @@ int main(int argc, char *argv[])
    args.AddOption(&vis_steps, "-vs", "--visualization-steps",
                   "Visualize every n-th timestep.");
    args.Parse();
+   AIR.prerelax = std::string(temp_prerelax);
+   AIR.postrelax = std::string(temp_postrelax);
+   if (trisolve) use_gmres = 1;
    if (!args.Good())
    {
       if (myid == 0)
@@ -155,11 +218,21 @@ int main(int argc, char *argv[])
    ODESolver *ode_solver = NULL;
    switch (ode_solver_type)
    {
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(1.0); break;
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
-      case 6: ode_solver = new RK6Solver; break;
+      // Implicit L-stable methods
+      case 1:  ode_solver = new BackwardEulerSolver; break;
+      case 2:  ode_solver = new SDIRK23Solver(2); break;
+      case 3:  ode_solver = new SDIRK33Solver; break;
+      // Explicit methods
+      case 11: ode_solver = new ForwardEulerSolver; break;
+      case 12: ode_solver = new RK2Solver(1.0); break;
+      case 13: ode_solver = new RK3SSPSolver; break;
+      case 14: ode_solver = new RK4Solver; break;
+      case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      case 16: ode_solver = new RK6Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
       default:
          if (myid == 0)
          {
@@ -195,8 +268,10 @@ int main(int argc, char *argv[])
    }
 
    // 7. Define the parallel discontinuous DG finite element space on the
-   //    parallel refined mesh of the given polynomial order.
-   DG_FECollection fec(order, dim);
+   //    parallel refined mesh of the given polynomial order. Basis_type=1
+   //    gives Gauss-Lobatto quadrature points, which are preferable for
+   //    nonsymmetric AMG implict solves. 
+   DG_FECollection fec(order, dim, basis_type);
    ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
 
    HYPRE_Int global_vSize = fes->GlobalTrueVSize();
@@ -314,7 +389,7 @@ int main(int argc, char *argv[])
    // 10. Define the time-dependent evolution operator describing the ODE
    //     right-hand side, and perform time-integration (looping over the time
    //     iterations, ti, with a time-step dt).
-   FE_Evolution adv(*M, *K, *B);
+   FE_Evolution adv(*M, *K, *B, order, AIR);
 
    double t = 0.0;
    adv.SetTime(t);
@@ -387,9 +462,11 @@ int main(int argc, char *argv[])
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_K,
-                           const Vector &_b)
-   : TimeDependentOperator(_M.Height()),
-     M(_M), K(_K), b(_b), M_solver(M.GetComm()), z(_M.Height())
+                           const Vector &_b, int order,
+                           AIR_parameters &_AIR)
+   : TimeDependentOperator(_M.Height()), A(NULL), AMG_solver(NULL),
+     GMRES_solver(NULL), preconditioner(NULL), M(_M), K(_K), b(_b),
+     M_solver(M.GetComm()), z(_M.Height()), AIR(_AIR)
 {
    M_prec.SetType(HypreSmoother::Jacobi);
    M_solver.SetPreconditioner(M_prec);
@@ -400,8 +477,23 @@ FE_Evolution::FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_K,
    M_solver.SetAbsTol(0.0);
    M_solver.SetMaxIter(100);
    M_solver.SetPrintLevel(0);
+
+   // DG block size given by (FEorder+1)^2 on square meshes.
+   blocksize = (order+1)*(order+1);
+   dt = -1;
 }
 
+
+FE_Evolution::~FE_Evolution()
+{
+   BlockInvScal(NULL, NULL, NULL, NULL, 0, -1);
+   if (A) delete A;
+   if (AMG_solver) delete AMG_solver;
+   if (GMRES_solver) delete GMRES_solver;
+   if (preconditioner) delete preconditioner;
+}
+
+
 void FE_Evolution::Mult(const Vector &x, Vector &y) const
 {
    // y = M^{-1} (K x + b)
@@ -411,6 +503,74 @@ void FE_Evolution::Mult(const Vector &x, Vector &y) const
 }
 
 
+// Solve the equation:
+//    u_t = M^{-1}(Ku + b), 
+// by solving associated linear system
+//    (M - dt*K) d = K*u + b
+void FE_Evolution::ImplicitSolve(const double dt_, const Vector &u, Vector &du_dt)
+{
+   // if A is NULL or dt has changed since A was built, rebuild matrix and solver.
+   if ( (fabs(dt - dt_) > 1e-4 * dt) || !A ) {
+      delete GMRES_solver;
+      delete AMG_solver;
+      delete preconditioner;
+      delete A;
+
+      dt = dt_;
+      A = HypreParMatrixAdd(1.0, M, -1.0*dt, K);
+
+      // Scale A by block-diagonal inverse
+      BlockInvScal(A, &A_s, NULL, NULL, blocksize, 0);
+
+      int print_level = 1;
+      if (!trisolve) {
+         AMG_solver = new HypreBoomerAMG(A_s);
+         AMG_solver->SetLAIROptions(AIR.distanceR, AIR.prerelax, AIR.postrelax,
+                                    AIR.strength_tolC, AIR.strength_tolR, AIR.filter_tolR,
+                                    AIR.interp_type, AIR.relax_type, AIR.filterA_tol,
+                                    AIR.coarsen_type, -1, 1);  
+         AMG_solver->SetMaxLevels(50);
+         if (use_gmres) {
+            GMRES_solver = new HypreGMRES(A_s);
+            GMRES_solver->SetTol(1e-12);
+            GMRES_solver->SetMaxIter(100);
+            GMRES_solver->SetPrintLevel(print_level);
+            GMRES_solver->SetPreconditioner(*AMG_solver);
+            GMRES_solver->iterative_mode = false;
+         }
+         else {
+            AMG_solver->SetPrintLevel(print_level);
+            AMG_solver->SetTol(1e-12);
+            AMG_solver->SetMaxIter(100);
+         }
+      }
+      else {
+         preconditioner = new HypreTriSolve();
+         GMRES_solver = new HypreGMRES(A_s);
+         GMRES_solver->SetTol(1e-12);
+         GMRES_solver->SetMaxIter(100);
+         GMRES_solver->SetPrintLevel(print_level);
+         GMRES_solver->SetPreconditioner(*preconditioner);
+         GMRES_solver->SetZeroInintialIterate();
+         GMRES_solver->iterative_mode = false;
+      }
+   }
+
+   K.Mult(u, z);
+   z += b;
+
+   // scale the rhs and solve system 
+   HypreParVector z_s;
+   BlockInvScal(A, NULL, &z, &z_s, blocksize, 2);
+   if (use_gmres){
+      GMRES_solver->Mult(z_s, du_dt);
+   }
+   else {
+      AMG_solver->Mult(z_s, du_dt);
+   }
+}
+
+
 // Velocity coefficient
 void velocity_function(const Vector &x, Vector &v)
 {

examples/makefile

@@ -22,9 +22,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex21
+  ex18 ex19 ex20 ex21 ex23
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p MFEM_adv ex23TRp exETRp
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -117,7 +117,7 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -rf Example5* Example9* Example15* Example16*
+	@rm -rf Example5* Example9* Example15* Example16* Example23*
 	@rm -f sphere_refined.* sol.* sol_u.* sol_p.*
 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
 	@rm -f deformed.* velocity.* elastic_energy.* mode_*
@@ -125,4 +125,9 @@ clean-exec:
 	@rm -f vortex-mesh.* vortex.mesh vortex-?-init.* vortex-?-final.*
 	@rm -f deformation.* pressure.*
 	@rm -f ex20.dat ex20p_?????.dat gnuplot_ex20.inp gnuplot_ex20p.inp
+<<<<<<< HEAD
+	@rm -f ex22*.mesh ex22*.sol ex22p_*.*
+	@rm -f ex23.mesh ex23-mesh.* ex23-init.* ex23-final.*
+=======
 	@rm -f ex21*.mesh ex21*.sol ex21p_*.*
+>>>>>>> f8a3a379d13841c6e63b9d8fbc868aa325b8afc0

linalg/hypre.cpp

@@ -172,6 +172,13 @@ HypreParVector::HypreParVector(ParFiniteElementSpace *pfes)
    _SetDataAndSize_();
    own_ParVector = 1;
 }
+  
+void HypreParVector::WrapHypreParVector(hypre_ParVector *y)
+{
+   x = y;
+   _SetDataAndSize_();
+   own_ParVector = 0;
+}
 
 Vector * HypreParVector::GlobalVector() const
 {
@@ -974,17 +981,26 @@ static void MakeWrapper(const hypre_CSRMatrix *mat, SparseMatrix &wrapper)
    wrapper.Swap(tmp);
 }
 
+
 void HypreParMatrix::GetDiag(SparseMatrix &diag) const
 {
    MakeWrapper(A->diag, diag);
 }
 
+
 void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const
 {
    MakeWrapper(A->offd, offd);
    cmap = A->col_map_offd;
 }
 
+
+void HypreParMatrix::GetProcRows(SparseMatrix &colCSRMat)
+{
+   MakeWrapper(hypre_MergeDiagAndOffd(A), colCSRMat);
+}
+
+
 void HypreParMatrix::GetBlocks(Array2D<HypreParMatrix*> &blocks,
                                bool interleaved_rows,
                                bool interleaved_cols) const
@@ -1545,6 +1561,40 @@ void HypreParMatrix::Destroy()
    }
 }
 
+/* job =  0, extract block diagonal of A and scale A into C
+ * job =  1, job 0 + scale b into d
+ * job =  2, use A to scale b only
+ */
+int BlockInvScal(const HypreParMatrix *A, HypreParMatrix *C,
+                 const Vector *b, HypreParVector *d, int block, int job)
+{
+   if (0 == job || 1 == job)
+   {
+      hypre_ParCSRMatrix *C_hypre;
+      hypre_ParcsrBdiagInvScal(*A, block, &C_hypre);
+      /* XXX: FIXME drop in BdiagInvScal */
+      hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
+      (*C).WrapHypreParCSRMatrix(C_hypre);
+   }
+
+   if (1 == job || 2 == job)
+   {
+      HypreParVector *b_Hypre = new HypreParVector(A->GetComm(), A->GetGlobalNumRows(), 
+                                                   b->GetData(), A->GetRowStarts());
+      hypre_ParVector *d_hypre;
+      hypre_ParvecBdiagInvScal(*b_Hypre, block, &d_hypre, *A);
+      
+      delete b_Hypre;
+
+      d->WrapHypreParVector(d_hypre);
+      d->SetOwnership(true);
+
+      return 0;
+   }
+
+   return -1;
+}
+
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
                     double beta,  const HypreParMatrix &B)
 {
@@ -1562,6 +1612,18 @@ HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
    return C;
 }
 
+HypreParMatrix *HypreParMatrixAdd(double alpha, const HypreParMatrix &A,
+                                  double beta,  const HypreParMatrix &B)
+{
+   hypre_ParCSRMatrix *C_hypre;
+
+   hypre_ParcsrAdd(alpha, A, beta, B, &C_hypre);
+
+   HypreParMatrix *C = new HypreParMatrix(C_hypre);
+
+   return C;
+}
+
 HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B)
 {
    hypre_ParCSRMatrix * ab;
@@ -2363,6 +2425,12 @@ void HypreGMRES::SetTol(double tol)
    HYPRE_GMRESSetTol(gmres_solver, tol);
 }
 
+void HypreGMRES::SetAbsTol(double tol)
+{
+   HYPRE_GMRESSetTol(gmres_solver, 0.0);
+   HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
+}
+
 void HypreGMRES::SetMaxIter(int max_iter)
 {
    HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
@@ -2537,6 +2605,75 @@ HypreBoomerAMG::HypreBoomerAMG(HypreParMatrix &A) : HypreSolver(&A)
    SetDefaultOptions();
 }
 
+void HypreBoomerAMG::Mult(const HypreParVector &b, HypreParVector &x) const
+{
+   int myid;
+   HYPRE_Int time_index = 0;
+   HYPRE_Int num_iterations;
+   double final_res_norm;
+   MPI_Comm comm;
+   HYPRE_Int print_level;
+
+   HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
+
+   HYPRE_ParCSRMatrixGetComm(*A, &comm);
+
+   if (!setup_called)
+   {
+      if (print_level > 0)
+      {
+         time_index = hypre_InitializeTiming("BoomerAMG Setup");
+         hypre_BeginTiming(time_index);
+      }
+
+      HYPRE_BoomerAMGSetup(amg_precond, *A, b, x);
+      setup_called = 1;
+
+      if (print_level > 0)
+      {
+         hypre_EndTiming(time_index);
+         hypre_PrintTiming("Setup phase times", comm);
+         hypre_FinalizeTiming(time_index);
+         hypre_ClearTiming();
+      }
+   }
+
+   if (print_level > 0)
+   {
+      time_index = hypre_InitializeTiming("BoomerAMG Solve");
+      hypre_BeginTiming(time_index);
+   }
+
+   if (!iterative_mode)
+   {
+      x = 0.0;
+   }
+
+   HYPRE_BoomerAMGSolve(amg_precond, *A, b, x);
+
+   if (print_level > 0)
+   {
+      hypre_EndTiming(time_index);
+      hypre_PrintTiming("Solve phase times", comm);
+      hypre_FinalizeTiming(time_index);
+      hypre_ClearTiming();
+
+      HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_iterations);
+      HYPRE_BoomerAMGGetFinalRelativeResidualNorm(amg_precond,
+                                                    &final_res_norm);
+
+      MPI_Comm_rank(comm, &myid);
+
+      if (myid == 0)
+      {
+         mfem::out << "BoomerAMG Iterations = " << num_iterations << endl
+                   << "Final Relative Residual Norm = " << final_res_norm
+                   << endl;
+      }
+   }
+}
+
+
 void HypreBoomerAMG::SetDefaultOptions()
 {
    // AMG coarsening options:
@@ -2784,6 +2921,234 @@ void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace)
    error_mode = IGNORE_HYPRE_ERRORS;
 }
 
+void HypreBoomerAMG::SetCoord(int coord_dim, float *coord)
+{
+   HYPRE_BoomerAMGSetPlotGrids (amg_precond, 1);
+   //HYPRE_BoomerAMGSetPlotFileName (amg_precond, plot_file_name);
+   HYPRE_BoomerAMGSetCoordDim (amg_precond, coord_dim);
+   HYPRE_BoomerAMGSetCoordinates (amg_precond, coord);
+}
+
+
+void HypreBoomerAMG::SetLAIROptions(int distance,
+                                    std::string prerelax,
+                                    std::string postrelax, 
+                                    double strength_tolC,
+                                    double strength_tolR,
+                                    double filter_tolR,
+                                    int interp_type, 
+                                    int relax_type,
+                                    double filterA_tol, 
+                                    int splitting,
+                                    int blksize,
+                                    int Sabs)
+{
+   int ns_down, ns_up, ns_coarse;
+   if (distance > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+      std::string F("F");
+      std::string C("C");
+      std::string A("A");
+
+      // Array to store relaxation scheme and pass to Hypre
+      int **grid_relax_points = (int **) malloc(4*sizeof(int *));
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = (int *) malloc(sizeof(int)*ns_down);
+      grid_relax_points[2] = (int *) malloc(sizeof(int)*ns_up);
+      grid_relax_points[3] = (int *) malloc(sizeof(int));
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme 
+      for(unsigned int i = 0; i<ns_down; i++) {
+         if (prerelax.compare(i,1,F) == 0) {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax.compare(i,1,C) == 0) {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax.compare(i,1,A) == 0) {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme 
+      for(unsigned int i = 0; i<ns_up; i++) {
+         if (postrelax.compare(i,1,F) == 0) {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax.compare(i,1,C) == 0) {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax.compare(i,1,A) == 0) {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, distance);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+   
+   //HYPRE_BoomerAMGSetMaxRowSum(amg_precond, 0.8);
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   if (blksize > 0)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+      //HYPRE_BoomerAMGSetNodalLevels(amg_precond, 1);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, splitting);
+   
+   /* does not support aggressive coarsening */
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+   
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+   
+   if (distance > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+      HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (distance > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      /* type = -1: drop based on row inf-norm */
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+
+   //HYPRE_BoomerAMGSetMaxCoarseSize(amg_precond, 1000);
+}
+
+
+void HypreBoomerAMG::SetNAIROptions(int neumann_degree,
+                                    std::string prerelax,
+                                    std::string postrelax, 
+                                    double strength_tolC,
+                                    double strength_tolR,
+                                    double filter_tolR,
+                                    int interp_type, 
+                                    int relax_type,
+                                    double filterA_tol, 
+                                    int splitting,
+                                    int blksize,
+                                    int Sabs)
+{
+   int ns_down, ns_up, ns_coarse;
+   if (neumann_degree > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+      std::string F("F");
+      std::string C("C");
+      std::string A("A");
+
+      // Array to store relaxation scheme and pass to Hypre
+      int **grid_relax_points = (int **) malloc(4*sizeof(int *));
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = (int *) malloc(sizeof(int)*ns_down);
+      grid_relax_points[2] = (int *) malloc(sizeof(int)*ns_up);
+      grid_relax_points[3] = (int *) malloc(sizeof(int));
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme 
+      for(unsigned int i = 0; i<ns_down; i++) {
+         if (prerelax.compare(i,1,F) == 0) {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax.compare(i,1,C) == 0) {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax.compare(i,1,A) == 0) {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme 
+      for(unsigned int i = 0; i<ns_up; i++) {
+         if (postrelax.compare(i,1,F) == 0) {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax.compare(i,1,C) == 0) {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax.compare(i,1,A) == 0) {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, 3+neumann_degree);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+   
+   //HYPRE_BoomerAMGSetMaxRowSum(amg_precond, 0.8);
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   if (blksize > 0)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+      //HYPRE_BoomerAMGSetNodalLevels(amg_precond, 1);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, splitting);
+   
+   /* does not support aggressive coarsening */
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+   
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+   
+   if (neumann_degree > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (neumann_degree > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      /* type = -1: drop based on row inf-norm */
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+
+   //HYPRE_BoomerAMGSetMaxCoarseSize(amg_precond, 1000);
+}
+
+
 HypreBoomerAMG::~HypreBoomerAMG()
 {
    for (int i = 0; i < rbms.Size(); i++)
@@ -3796,4 +4161,4 @@ HypreAME::StealEigenvectors()
 
 }
 
-#endif
+#endif

linalg/hypre.hpp

@@ -84,6 +84,9 @@ private:
    inline void _SetDataAndSize_();
 
 public:
+
+  HypreParVector() {}
+  
    /** @brief Creates vector with given global size and parallel partitioning of
        the rows/columns given by @a col. */
    /** @anchor hypre_partitioning_descr
@@ -116,9 +119,9 @@ public:
    /// MPI communicator
    MPI_Comm GetComm() { return x->comm; }
 
-   /// Returns the parallel row/column partitioning
-   /** See @ref hypre_partitioning_descr "here" for a description of the
-       partitioning array. */
+   void WrapHypreParVector(hypre_ParVector *y);
+   
+   /// Returns the row partitioning
    inline HYPRE_Int *Partitioning() { return x->partitioning; }
 
    /// Returns the global number of rows
@@ -234,22 +237,25 @@ public:
    /// An empty matrix to be used as a reference to an existing matrix
    HypreParMatrix();
 
-   /// Converts hypre's format to HypreParMatrix
-   /** If @a owner is false, ownership of @a a is not transferred */
-   explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
+   void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner = true)
    {
-      Init();
       A = a;
       if (!owner) { ParCSROwner = 0; }
       height = GetNumRows();
       width = GetNumCols();
    }
 
-   /// Creates block-diagonal square parallel matrix.
-   /** Diagonal is given by @a diag which must be in CSR format (finalized). The
-       new HypreParMatrix does not take ownership of any of the input arrays.
-       See @ref hypre_partitioning_descr "here" for a description of the row
-       partitioning array @a row_starts.
+   /// Converts hypre's format to HypreParMatrix
+   /** If @a owner is false, ownership of @a a is not transferred */
+   explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
+   {
+      Init();
+      WrapHypreParCSRMatrix(a, owner);
+   }
+
+   /** Creates block-diagonal square parallel matrix. Diagonal is given by diag
+       which must be in CSR format (finalized). The new HypreParMatrix does not
+       take ownership of any of the input arrays.
 
        @warning The ordering of the columns in each row in @a *diag may be
        changed by this constructor to ensure that the first entry in each row is
@@ -393,6 +399,8 @@ public:
    void GetDiag(SparseMatrix &diag) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;
+   /// Get on-processor rows as CSR matrix.
+   void GetProcRows(SparseMatrix &colCSRMat);
 
    /** Split the matrix into M x N equally sized blocks of parallel matrices.
        The size of 'blocks' must already be set to M x N. */
@@ -541,12 +549,18 @@ public:
    Type GetType() const { return Hypre_ParCSR; }
 };
 
+int BlockInvScal(const HypreParMatrix *A, HypreParMatrix *C,
+                 const Vector *b, HypreParVector *d, int block, int job);
+
 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
     and `B` use the same row and column partitions and the same `col_map_offd`
     arrays. */
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
                     double beta,  const HypreParMatrix &B);
 
+HypreParMatrix *HypreParMatrixAdd(double alpha, const HypreParMatrix &A,
+                                  double beta,  const HypreParMatrix &B);
+
 /// Returns the matrix A * B
 HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B);
 /// Returns the matrix A + B
@@ -627,7 +641,7 @@ public:
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
    enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
-               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
+               GS = 6, TS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
 
    HypreSmoother();
@@ -729,6 +743,25 @@ public:
    virtual ~HypreSolver();
 };
 
+
+/// Abstract class for hypre's solvers and preconditioners
+class HypreTriSolve : public HypreSolver
+{
+public:
+   HypreTriSolve() : HypreSolver() { }
+   explicit HypreTriSolve(HypreParMatrix &A) : HypreSolver(&A) { }
+   virtual operator HYPRE_Solver() const { return NULL; }
+
+   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSetup; }
+   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSolve; }
+
+   HypreParMatrix* GetData() { return A; }
+   virtual ~HypreTriSolve() { }
+};
+
+
 /// PCG solver in hypre
 class HyprePCG : public HypreSolver
 {
@@ -788,6 +821,7 @@ public:
    HypreGMRES(HypreParMatrix &_A);
 
    void SetTol(double tol);
+   void SetAbsTol(double tol);
    void SetMaxIter(int max_iter);
    void SetKDim(int dim);
    void SetLogging(int logging);
@@ -938,9 +972,79 @@ public:
        As with SetSystemsOptions(), this solver assumes Ordering::byVDIM. */
    void SetElasticityOptions(ParFiniteElementSpace *fespace);
 
+   /* distance parameter takes on values {1,2,15} for lAIR, meaning R is built using
+   distance 1 neighbors, distance two neighbors, or distance two on processor and
+   distance 1 off processor (i.e., distance 1.5 --> 15).        */
+   void SetLAIROptions(int distance=15,  std::string prerelax="",
+                      std::string postrelax="FFC", double strength_tol=0.1,
+                      double strength_tolR=0.01, double filter_tolR=0.0,
+                      int interp_type=100, int relax_type=3, double filterA_tol=0.0, 
+                      int splitting=6, int blksize=0, int Sabs=0);
+
+   void SetNAIROptions(int neumann_degree=2,  std::string prerelax="A",
+                      std::string postrelax="F", double strength_tol=0.1,
+                      double strength_tolR=0.01, double filter_tolR=0.0,
+                      int interp_type=100, int relax_type=10, double filterA_tol=0.0, 
+                      int splitting=6, int blksize=0, int Sabs=0);
+
+   void SetCoord(int dim, float *coord);
+   
    void SetPrintLevel(int print_level)
    { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
 
+   void SetMaxIter(int max_iter)
+   { HYPRE_BoomerAMGSetMaxIter(amg_precond, max_iter); }
+
+   void SetMaxLevels(int max_levels)
+   { HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); }
+
+   void SetTol(double tol)
+   { HYPRE_BoomerAMGSetTol(amg_precond, tol); }
+
+   void SetStrengthThresh(double strength)
+   { HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength); }
+
+   void SetStrengthThreshR(double strengthR)
+   { HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strengthR); }
+
+   void SetFilterThreshR(double filterR)
+   { HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filterR); }
+
+   void SetInterpolation(int interp_type)
+   { HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); }
+
+   void SetRestriction(int restrict_type)
+   { HYPRE_BoomerAMGSetRestriction(amg_precond, restrict_type); }
+
+   void SetCoarsening(int coarsen_type)
+   { HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); }
+
+   void SetRelaxType(int relax_type)
+   { HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); }
+
+   void SetRelaxCycle(int prerelax, int postrelax)
+   { HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, prerelax,  1);
+     HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, postrelax, 2); }
+
+   void GetNumIterations(int &num_it)
+   { HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_it); }
+
+   void SetCycleType(int cycle_type)
+   { HYPRE_BoomerAMGSetCycleType(amg_precond, cycle_type); }
+
+   void SetNodal(int blocksize)
+   { HYPRE_BoomerAMGSetNumFunctions(amg_precond, blocksize);
+     HYPRE_BoomerAMGSetNodal(amg_precond, 1); }
+
+   void SetAggressiveCoarsening(int num_levels)
+   { HYPRE_BoomerAMGSetAggNumLevels(amg_precond, num_levels); }
+
+   void SetTriangular()
+   { HYPRE_BoomerAMGSetIsTriangular(amg_precond, 1); }
+
+   void SetGMRESSwitchR(int gmres_switch)
+   { HYPRE_BoomerAMGSetGMRESSwitchR(amg_precond, gmres_switch); }
+
    /// The typecast to HYPRE_Solver returns the internal amg_precond
    virtual operator HYPRE_Solver() const { return amg_precond; }
 
@@ -949,6 +1053,9 @@ public:
    virtual HYPRE_PtrToParSolverFcn SolveFcn() const
    { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
 
+   virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
+   using HypreSolver::Mult;
+
    virtual ~HypreBoomerAMG();
 };