Reducing the two nested GMRES solvers to just one. Also implemented a complex version of AMS.

(this commit also improves some documentation in auxiliary.hpp)

CHANGELOG

@@ -11,6 +11,9 @@
 Version 4.2.1 (development)
 ===========================
 
+- Added high-order matrix-free auxiliary Maxwell solver for H(curl) problems,
+  as described in Barker and Kolev 2020 (https://doi.org/10.1002/nla.2348).
+
 - Added matrix-free GPU-enabled implementations of GradientInterpolator and
   IdentityInterpolator.
 

examples/ex1p_complex.cpp

@@ -0,0 +1,435 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p_complex
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+double E_exact(const Vector &);
+void gradE_exact(const Vector &, Vector &);
+double f_exact(const Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+
+#define COMPLEX_VERSION
+#define NEUMANN
+
+const double omega = 1.4;
+const double eps = 1.0e-8;
+
+
+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/beam-tet.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (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);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 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 1,000 elements.
+   {
+      int ref_levels = (int)floor(log(1000./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. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   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;
+
+#ifndef NEUMANN
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+#endif
+
+   //const double imscale = 0.0;
+   const double imscale = -omega;
+
+   Coefficient *im = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);  // im part
+
+   FunctionCoefficient E_coef(E_exact);
+   VectorFunctionCoefficient grad_E(sdim, gradE_exact);
+   ProductCoefficient omegaE(imscale, E_coef);
+
+   // 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
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   FunctionCoefficient f(f_exact);
+#ifdef COMPLEX_VERSION
+   ParComplexLinearForm *b = new ParComplexLinearForm(fespace);
+   b->AddDomainIntegrator(new DomainLFIntegrator(f), NULL);
+#ifdef NEUMANN
+   b->AddBoundaryIntegrator(NULL, new BoundaryNormalLFIntegrator(grad_E));
+#endif
+   b->AddBoundaryIntegrator(NULL, new BoundaryLFIntegrator(omegaE)); // im part
+#endif
+
+   b->Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x by projecting the exact
+   //     solution. Note that only values from the boundary edges will be used
+   //     when eliminating the non-homogeneous boundary condition to modify the
+   //     r.h.s. vector b.
+   /*
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+   */
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   ConstantCoefficient E_im(0.0);
+   //x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
+   x.ProjectCoefficient(E_coef, E_im);
+#endif
+
+   // 11. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+   Coefficient *epscoef = new ConstantCoefficient(eps);
+   Coefficient *imabs = new ConstantCoefficient(fabs(imscale));  // im part
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new DiffusionIntegrator(*muinv), NULL);
+   a->AddDomainIntegrator(new MassIntegrator(*epscoef), NULL);
+   a->AddBoundaryIntegrator(NULL, new MassIntegrator(*im));  // im part
+#endif
+
+   // 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);
+
+   ParBilinearForm a_Re(fespace);
+   a_Re.AddDomainIntegrator(new DiffusionIntegrator(*muinv));
+   a_Re.AddDomainIntegrator(new MassIntegrator(*epscoef));
+
+   if (pa) { a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a_Re.Assemble();
+
+   OperatorPtr A_Re;
+   a_Re.FormSystemMatrix(ess_tdof_list, A_Re);
+
+   ParBilinearForm a_Im(fespace);
+   a_Im.AddBoundaryIntegrator(new MassIntegrator(*imabs));
+   a_Im.Assemble();
+
+   OperatorPtr A_Im;
+   a_Im.FormSystemMatrix(ess_tdof_list, A_Im);
+
+   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
+   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
+   //     partial assembly case).
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = fespace->GetTrueVSize();
+   offsets[2] = fespace->GetTrueVSize();
+   offsets.PartialSum();
+
+   //OperatorJacobiSmoother massJacobi(a_Im, ess_tdof_list);
+
+   StopWatch sw;
+   sw.Clear();
+   sw.Start();
+
+   if (pa) // Jacobi preconditioning in partial assembly mode
+   {
+      MFEM_VERIFY(false, "TODO");
+      //OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(1000);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      //cg.SetPreconditioner(Jacobi);
+      cg.Mult(B, X);
+   }
+   else
+   {
+      if (myid == 0)
+      {
+         cout << "Size of linear system: "
+              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+      }
+
+      HypreBoomerAMG amg(*A_Re.As<HypreParMatrix>());
+
+#ifdef COMPLEX_VERSION
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, &amg);
+      BlockDP.SetDiagonalBlock(1, &amg);
+
+      Complex_PMHSS PMHSS(A_Re.Ptr(), A_Im.Ptr(), &BlockDP, NULL, 1.0);
+
+      ComplexOperator AspdComplex(A_Re.Ptr(), A_Im.Ptr(), false, false);
+
+      GMRESSolver PMHSSgmres(MPI_COMM_WORLD);
+      PMHSSgmres.SetPrintLevel(1);
+      PMHSSgmres.SetKDim(100);
+      PMHSSgmres.SetMaxIter(100);
+      PMHSSgmres.SetRelTol(1e-6);
+      PMHSSgmres.SetAbsTol(0.0);
+      PMHSSgmres.SetOperator(AspdComplex);
+      PMHSSgmres.SetPreconditioner(PMHSS);
+
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      //gmres.SetPreconditioner(BlockDP);
+      gmres.SetPreconditioner(PMHSS);
+      //gmres.SetPreconditioner(PMHSSgmres);
+#else
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      gmres.SetPreconditioner(ams);
+#endif
+
+      gmres.Mult(B, X);
+   }
+
+   sw.Stop();
+   mfem::out << "Total solve time " <<sw.RealTime() << endl;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+#ifdef COMPLEX_VERSION
+      double err = x.real().ComputeL2Error(E_coef);
+#else
+      double err = x.ComputeL2Error(E_coef);
+#endif
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 16. 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);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+#ifdef COMPLEX_VERSION
+      x.real().Save(sol_ofs);
+#else
+      x.Save(sol_ofs);
+#endif
+   }
+
+   // 17. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+#ifdef COMPLEX_VERSION
+      sol_sock << "solution\n" << *pmesh << x.real() << flush;
+#else
+      sol_sock << "solution\n" << *pmesh << x << flush;
+#endif
+   }
+
+   // 18. Free the used memory.
+   delete a;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+#define VERSION_COS
+
+double E_exact(const Vector &x)
+{
+   if (dim == 3)
+   {
+#ifdef VERSION_COS
+      return cos(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+#else
+      return sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+#endif
+   }
+   else
+   {
+      return 0.0;
+   }
+}
+
+void gradE_exact(const Vector &x, Vector &grad)
+{
+   if (dim == 3)
+   {
+#ifdef VERSION_COS
+      grad(0) = -kappa * sin(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+      grad(1) = -kappa * sin(kappa * x(1)) * cos(kappa * x(0)) * cos(kappa * x(2));
+      grad(2) = -kappa * sin(kappa * x(2)) * cos(kappa * x(0)) * cos(kappa * x(1));
+#else
+      grad(0) = kappa * cos(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+      grad(1) = kappa * cos(kappa * x(1)) * sin(kappa * x(0)) * sin(kappa * x(2));
+      grad(2) = kappa * cos(kappa * x(2)) * sin(kappa * x(0)) * sin(kappa * x(1));
+#endif
+   }
+   else
+   {
+      MFEM_VERIFY(false, "");
+   }
+}
+
+// (grad u, grad v) + eps (u, v) = <grad u . n, v> - (div grad u, v) + eps (u, v)
+double f_exact(const Vector &x)
+{
+   if (dim == 3)
+   {
+      const double c = 3.0 * kappa * kappa;
+#ifdef VERSION_COS
+      return (eps + c) * cos(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+#else
+      return (eps + c) * sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+#endif
+   }
+   else
+   {
+      return 0.0;
+   }
+}

examples/ex3p.cpp

@@ -69,7 +69,10 @@ int main(int argc, char *argv[])
    bool static_cond = false;
    bool pa = false;
    const char *device_config = "cpu";
-   bool visualization = 1;
+   bool visualization = true;
+#ifdef MFEM_USE_AMGX
+   bool useAmgX = false;
+#endif
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -87,6 +90,11 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+#ifdef MFEM_USE_AMGX
+   args.AddOption(&useAmgX, "-amgx", "--useAmgX", "-no-amgx",
+                  "--no-useAmgX",
+                  "Enable or disable AmgX in MatrixFreeAMS.");
+#endif
 
    args.Parse();
    if (!args.Good())
@@ -159,9 +167,10 @@ int main(int argc, char *argv[])
    //    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;
+   Array<int> ess_bdr;
    if (pmesh->bdr_attributes.Size())
    {
-      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
       ess_bdr = 1;
       fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
@@ -205,20 +214,20 @@ int main(int argc, char *argv[])
    Vector B, X;
    a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
 
-   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
-   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
-   //     partial assembly case).
-
-   if (pa) // Jacobi preconditioning in partial assembly mode
+   // 13. Solve the system AX=B using PCG with an AMS preconditioner.
+   if (pa)
    {
-      OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
-
+#ifdef MFEM_USE_AMGX
+      MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr, useAmgX);
+#else
+      MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr);
+#endif
       CGSolver cg(MPI_COMM_WORLD);
       cg.SetRelTol(1e-12);
       cg.SetMaxIter(1000);
       cg.SetPrintLevel(1);
       cg.SetOperator(*A);
-      cg.SetPreconditioner(Jacobi);
+      cg.SetPreconditioner(ams);
       cg.Mult(B, X);
    }
    else

examples/ex3p_complex.cpp

@@ -0,0 +1,517 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p_complex
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+void E_exact(const Vector &, Vector &);
+void curlE_exact(const Vector &, Vector &);
+void f_exact(const Vector &, Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+
+#define COMPLEX_VERSION
+#define NEUMANN
+#define INDEFINITE
+
+const double omega = 1.4;
+
+
+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/beam-tet.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (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);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 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 1,000 elements.
+   {
+      int ref_levels = (int)floor(log(1000./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. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   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;
+   Array<int> ess_bdr;
+   ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+   ess_bdr = 0;
+
+#ifndef NEUMANN
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+#endif
+
+   //const double imscale = 0.0;
+   const double imscale = omega;
+
+   Coefficient *im = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);  // im part
+
+   VectorFunctionCoefficient E_Re(sdim, E_exact);
+   VectorFunctionCoefficient curlE_Re(sdim, curlE_exact);
+
+   ScalarVectorProductCoefficient omegaE(imscale, E_Re);  // im part
+   //ScalarVectorProductCoefficient omegaE(0.0, E_Re);  // im part
+
+   // 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
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(sdim, f_exact);
+#ifdef COMPLEX_VERSION
+   ParComplexLinearForm *b = new ParComplexLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f), NULL);
+   b->AddBoundaryIntegrator(NULL,
+                            new VectorFEDomainLFIntegrator(omegaE));  // im part
+#else
+   // Real version
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+#endif
+
+#ifdef NEUMANN
+   b->AddBoundaryIntegrator(new VectorFEBoundaryTangentLFIntegrator(curlE_Re),
+                            NULL);
+#endif
+
+   b->Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x by projecting the exact
+   //     solution. Note that only values from the boundary edges will be used
+   //     when eliminating the non-homogeneous boundary condition to modify the
+   //     r.h.s. vector b.
+   /*
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+   */
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   Vector zero(sdim);
+   zero = 0.0;
+   VectorConstantCoefficient E_Im(zero);
+   //x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
+   x.ProjectCoefficient(E_Re, E_Im);
+#else
+   ParGridFunction x(fespace);
+   x = 0.0;
+   x.ProjectCoefficient(E_Re);
+#endif
+
+   // 11. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+#ifdef INDEFINITE
+   Coefficient *sigma = new ConstantCoefficient(
+      -omega*omega); // indefinite -, definite +
+#else
+   Coefficient *sigma = new ConstantCoefficient(
+      omega*omega); // indefinite -, definite +
+#endif
+   Coefficient *abssigma = new ConstantCoefficient(omega*omega);
+   Coefficient *imabs = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *imabs = new ConstantCoefficient(0.0);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv), NULL);
+   //a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), new VectorFEMassIntegrator(*im));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), NULL);
+   a->AddBoundaryIntegrator(NULL, new VectorFEMassIntegrator(*im));  // im part
+#else
+   // Real version
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   //a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), new VectorFEMassIntegrator(*im));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+#endif
+
+   // 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);
+
+   ParBilinearForm a_Re(fespace);
+   a_Re.AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a_Re.AddDomainIntegrator(new VectorFEMassIntegrator(*abssigma));
+
+   //if (pa) { a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   a_Re.Assemble();
+
+   OperatorPtr A_Re;
+   a_Re.FormSystemMatrix(ess_tdof_list, A_Re);
+
+   ParBilinearForm a_Im(fespace);
+   a_Im.AddBoundaryIntegrator(new VectorFEMassIntegrator(*imabs));
+   a_Im.Assemble();
+
+   OperatorPtr A_Im;
+   a_Im.FormSystemMatrix(ess_tdof_list, A_Im);
+
+   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
+   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
+   //     partial assembly case).
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = fespace->GetTrueVSize();
+   offsets[2] = fespace->GetTrueVSize();
+   offsets.PartialSum();
+
+   //OperatorJacobiSmoother massJacobi(a_Im, ess_tdof_list);
+
+   StopWatch sw;
+   sw.Clear();
+   sw.Start();
+
+   if (pa) // Jacobi preconditioning in partial assembly mode
+   {
+      MFEM_VERIFY(false, "TODO");
+      //OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(1000);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      //cg.SetPreconditioner(Jacobi);
+      cg.Mult(B, X);
+   }
+   else
+   {
+      if (myid == 0)
+      {
+         cout << "Size of linear system: "
+              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+      }
+
+      //HypreAMS ams(*A_Re.As<HypreParMatrix>(), fespace);
+
+      // One option is to use the standard real-valued MatrixFreeAMS to precondition
+      // the real part of the complex system in the PMHSS preconditioner (BlockDiagonalPreconditioner).
+      // Another option is to use complex MatrixFreeAMS to precondition the
+      // complex system without PMHSS and without a BlockDiagonalPreconditioner.
+      //#define COMPLEX_AMS
+
+#ifdef MFEM_USE_AMGX
+      bool useAmgX = false;
+      cout << "Built with AMGX, using AMGX " << useAmgX << endl;
+      MatrixFreeAMS ams(a_Re, *A_Re, *fespace, muinv, abssigma, im, imabs, NULL,
+                        ess_bdr, useAmgX);
+      MatrixFreeAMS ams(a_Re, *A_Re, *fespace, muinv, abssigma, NULL, NULL, ess_bdr,
+                        useAmgX);
+#ifdef COMPLEX_AMS
+      MFEM_VERIFY(false, "TODO");
+#endif
+
+#else
+      cout << "Not built with AMGX" << endl;
+#ifdef COMPLEX_AMS
+      MatrixFreeAMS ams(a_Re, *A_Re, A.Ptr(), *fespace, muinv, abssigma, im, imabs,
+                        NULL, ess_bdr);
+#else
+      MatrixFreeAMS ams(a_Re, *A_Re, NULL, *fespace, muinv, abssigma, NULL, NULL,
+                        NULL, ess_bdr);
+#endif
+#endif
+
+#ifdef COMPLEX_VERSION
+
+#ifdef COMPLEX_AMS
+      //MFEM_VERIFY(false, "TODO");
+#else
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, &ams);
+      BlockDP.SetDiagonalBlock(1, &ams);
+
+      /*
+      BlockDiagonalPreconditioner BlockDP_Im(offsets);
+      BlockDP_Im.SetDiagonalBlock(0, &massJacobi); // TODO: this won't work if it has zeros on diagonal
+      BlockDP_Im.SetDiagonalBlock(1, &massJacobi);
+      */
+
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, &BlockDP_Im);
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, NULL, 2.0 * omega);
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, NULL, omega);
+      Complex_PMHSS PMHSS(A_Re.Ptr(), A_Im.Ptr(), &BlockDP, NULL, 1.0);
+
+      ComplexOperator AspdComplex(A_Re.Ptr(), A_Im.Ptr(), false, false);
+
+      GMRESSolver PMHSSgmres(MPI_COMM_WORLD);
+      PMHSSgmres.SetPrintLevel(1);
+      PMHSSgmres.SetKDim(100);
+      PMHSSgmres.SetMaxIter(100);
+      PMHSSgmres.SetRelTol(1e-6);
+      PMHSSgmres.SetAbsTol(0.0);
+      PMHSSgmres.SetOperator(AspdComplex);
+      PMHSSgmres.SetPreconditioner(PMHSS);
+#endif
+
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      //gmres.SetPreconditioner(BlockDP);
+#ifdef COMPLEX_AMS
+      //MFEM_VERIFY(false, "TODO");
+      gmres.SetPreconditioner(ams);
+#else
+      gmres.SetPreconditioner(PMHSS);
+      //gmres.SetPreconditioner(PMHSSgmres);
+#endif
+
+#else
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      gmres.SetPreconditioner(ams);
+#endif
+
+      gmres.Mult(B, X);
+   }
+
+   sw.Stop();
+   mfem::out << "Total solve time " <<sw.RealTime() << endl;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+#ifdef COMPLEX_VERSION
+      double err = x.real().ComputeL2Error(E_Re);
+#else
+      double err = x.ComputeL2Error(E_Re);
+#endif
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 16. 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);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+#ifdef COMPLEX_VERSION
+      x.real().Save(sol_ofs);
+#else
+      x.Save(sol_ofs);
+#endif
+   }
+
+   // 17. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+#ifdef COMPLEX_VERSION
+      sol_sock << "solution\n" << *pmesh << x.real() << flush;
+#else
+      sol_sock << "solution\n" << *pmesh << x << flush;
+#endif
+   }
+
+   // 18. Free the used memory.
+   delete a;
+   delete sigma;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+
+void E_exact(const Vector &x, Vector &E)
+{
+   if (dim == 3)
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(2));
+      E(2) = sin(kappa * x(0));
+   }
+   else
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(0));
+      if (x.Size() == 3) { E(2) = 0.0; }
+   }
+}
+
+void curlE_exact(const Vector &x, Vector &curl)
+{
+   if (dim == 3)
+   {
+      curl(0) = kappa * cos(kappa * x(2));
+      curl(1) = kappa * cos(kappa * x(0));
+      curl(2) = kappa * cos(kappa * x(1));
+   }
+   else
+   {
+      MFEM_VERIFY(false, "");
+   }
+}
+
+void f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      // indefinite -m, definite +m
+      const double c = kappa * kappa;
+#ifdef INDEFINITE
+      const double m = -omega * omega;
+#else
+      const double m = omega * omega;
+#endif
+      f(0) = (c + m) * sin(kappa * x(1));
+      f(1) = (c + m) * sin(kappa * x(2));
+      f(2) = (c + m) * sin(kappa * x(0));
+   }
+   else
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}

linalg/CMakeLists.txt

@@ -10,6 +10,7 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SRCS
+  auxiliary.cpp
   blockmatrix.cpp
   blockoperator.cpp
   blockvector.cpp
@@ -27,6 +28,7 @@ list(APPEND SRCS
   )
 
 list(APPEND HDRS
+  auxiliary.hpp
   blockmatrix.hpp
   blockoperator.hpp
   blockvector.hpp

linalg/auxiliary.hpp

@@ -0,0 +1,285 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_AUXILIARY
+#define MFEM_AUXILIARY
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "../general/tic_toc.hpp"
+#include "solvers.hpp"
+
+namespace mfem
+{
+
+// forward declarations
+class Coefficient;
+class MatrixCoefficient;
+class ParMesh;
+class ParBilinearForm;
+class ParDiscreteLinearOperator;
+
+/** @brief Auxiliary space solvers for MatrixFreeAMS preconditioner
+
+    Given an operator A and a transfer G, this will create a solver
+    that approximates (G^T A G)^{-1}. Used for two different
+    auxiliary spaces in the AMS cycle.
+
+    The produced solver is based on a low-order refined discretization
+    for the high-order H1 problem. */
+class MatrixFreeAuxiliarySpace : public Solver
+{
+public:
+   /** @brief Pi space constructor
+
+       In the AMS framework this auxiliary space has two coefficients.
+
+       @param mesh_lor       Low-order refined auxiliary mesh
+       @param alpha_coeff    coefficient on curl-curl term (1 if null)
+       @param beta_coeff     coefficient on mass term (1 if null)
+       @param beta_mcoeff    matrix coefficient on mass term
+       @param ess_bdr        attributes for essential boundaries
+       @param curlcurl_oper  High-order operator for the system
+       @param pi             Intentity interpolation operator
+       @param useAmgX_       Use AmgX instead of hypre for auxiliary solves
+       @param cg_iterations  number of CG iterations used to invert
+                             auxiliary system, choosing 0 means to use a
+                             single V-cycle
+   */
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* alpha_coeff, Coefficient* beta_coeff,
+      MatrixCoefficient* beta_mcoeff,
+      Array<int>& ess_bdr, Operator& curlcurl_oper, Operator& pi,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 0);
+
+   // Complex Pi space constructor
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* alpha_coeff, Coefficient* beta_coeff,
+      Coefficient* beta_imag, Coefficient* abs_beta_imag,
+      MatrixCoefficient* beta_mcoeff,
+      Array<int>& ess_bdr, Operator& curlcurl_oper, Operator *oper_complex,
+      Operator& pi,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 0);
+
+   /** @brief G space constructor
+
+       This has one coefficient in the AMS framework.
+
+       @param mesh_lor       Low-order refined auxiliary mesh
+       @param beta_coeff     coefficient on mass term (1 if null)
+       @param beta_mcoeff    matrix coefficient on mass term
+       @param ess_bdr        attributes for essential boundaries
+       @param curlcurl_oper  High-order operator for the system
+       @param g              Gradient interpolation operator
+       @param useAmgX_       Use AmgX instead of hypre for auxiliary solves
+       @param cg_iterations  number of CG iterations used to invert
+                             auxiliary system, choosing 0 means to
+                             use a single V-cycle
+   */
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* beta_coeff,
+      MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+      Operator& curlcurl_oper, Operator& g,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 1);
+
+   // Complex G space constructor
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* beta_coeff, Coefficient* beta_imag,
+      Coefficient* abs_beta_imag,
+      MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+      Operator& curlcurl_oper, Operator *oper_complex, Operator& g,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 1);
+
+   ~MatrixFreeAuxiliarySpace();
+
+   void Mult(const Vector& x, Vector& y) const;
+
+   void SetOperator(const Operator& op) {}
+
+private:
+   /** @brief Helper routine for constructors.
+
+       @param system_dimension is passed to HypreBoomerAMG::SetSystemsOptions
+   */
+   void SetupAMG(int system_dimension);
+   void SetupVCycle();
+
+   /// inner_cg_iterations > 99 applies an exact solve here
+   void SetupCG(Operator& curlcurl_oper, Operator& conn,
+                int inner_cg_iterations);
+
+   void SetupGMRES(Operator& curlcurl_oper, Operator& conn);
+
+   void SetupPMHSS();
+
+   MPI_Comm comm;
+   Array<int> ess_tdof_list;
+   HypreParMatrix * aspacematrix;
+   HypreParMatrix * aspacematrix_complex;
+   HypreParMatrix * aspacematrix_imag;
+   Solver * aspacepc;
+   Operator* matfree;
+   CGSolver* cg;
+   GMRESSolver* gmres;
+   GMRESSolver* gmres_PMHSS;
+   Operator* aspacewrapper;
+#ifdef MFEM_USE_AMGX
+   const bool useAmgX;
+#endif
+   mutable int inner_aux_iterations;
+
+   const bool imagBdry;
+
+   Complex_PMHSS *PMHSS = NULL;
+
+   Array<int> offsets;
+   Array<int> offsets_nd;
+   BlockDiagonalPreconditioner *BlockDP;
+
+   BlockOperator *conn_block;
+};
+
+
+/** @brief Perform AMS cycle with generic Operator objects.
+
+    Most users should use MatrixFreeAMS, which wraps this. */
+class GeneralAMS : public Solver
+{
+public:
+   /** @brief Constructor.
+
+       Most of these arguments just need a Mult() operation,
+       but pi and g also require MultTranspose() */
+   GeneralAMS(const Operator& curlcurl_op_,
+              Operator *oper_complex,
+              const Operator& pi_,
+              const Operator& gradient_,
+              const Operator& pispacesolver_,
+              const Operator& gspacesolver_,
+              const Operator& smoother_,
+              const Array<int>& ess_tdof_list_);
+   virtual ~GeneralAMS();
+
+   /// in principle this should set A_ = op;
+   void SetOperator(const Operator &op) {}
+
+   virtual void Mult(const Vector& x, Vector& y) const;
+
+private:
+   const Operator& curlcurl_op;
+   Operator *oper_complex;
+   const Operator& pi;
+   const Operator& gradient;
+   const Operator& pispacesolver;
+   const Operator& gspacesolver;
+   const Operator& smoother;
+   const Array<int> ess_tdof_list;
+
+   void FormResidual(const Vector& rhs, const Vector& x,
+                     Vector& residual) const;
+};
+
+
+/** @brief An auxiliary Maxwell solver for a high-order curl-curl
+    system without high-order assembly.
+
+    The auxiliary space solves are done using a low-order refined approach,
+    but all the interpolation operators, residuals, etc. are done in a
+    matrix-free manner.
+
+    See Barker and Kolev, Matrix-free preconditioning for high-order H(curl)
+    discretizations (https://doi.org/10.1002/nla.2348) */
+class MatrixFreeAMS : public Solver
+{
+public:
+   /** @brief Construct matrix-free AMS preconditioner
+
+       @param aform        BilinearForm for curl-curl problem, generally will
+                           have a CurlCurlIntegrator and possibly a
+                           VectorFEMassIntegrator.
+       @param oper         Operator to precondition.
+       @param nd_fespace   Underlying Nedelec finite element space.
+       @param alpha_coeff  coefficient on curl-curl term in Maxwell problem
+                           (can be null, in which case constant 1 is assumed)
+       @param beta_coeff   (scalar) coefficient on mass term in Maxwell problem
+       @param beta_mcoeff  (matrix) coefficient on mass term
+       @param ess_bdr      boundary *attributes* that are marked essential. In
+                           contrast to other MFEM cases, these are *attributes*
+                           not dofs, because we need to apply these boundary
+                           conditions to different bilinear forms.
+       @param useAmgX      use AmgX (instead of hypre) for LOR problems
+       @param inner_pi_its number of CG iterations on auxiliary pi space,
+                           may need more for difficult coefficients
+       @param inner_g_its  number of CG iterations on auxiliary g space,
+                           may need more for difficult coefficients
+       @param nd_smoother  optional user-provided smoother for Nedelec space,
+                           this object takes ownership and will delete.
+   */
+   MatrixFreeAMS(ParBilinearForm& aform, Operator& oper, Operator *oper_complex,
+                 ParFiniteElementSpace& nd_fespace, Coefficient* alpha_coeff,
+                 Coefficient* beta_coeff, Coefficient* beta_imag,
+                 Coefficient* abs_beta_imag, MatrixCoefficient* beta_mcoeff,
+                 Array<int>& ess_bdr,
+#ifdef MFEM_USE_AMGX
+                 bool useAmgX = false,
+#endif
+                 int inner_pi_its = 0, int inner_g_its = 1,
+                 Solver* nd_smoother = NULL);
+
+   ~MatrixFreeAMS();
+
+   void SetOperator(const Operator &op) {}
+
+   void Mult(const Vector& x, Vector& y) const { general_ams->Mult(x, y); }
+
+private:
+   GeneralAMS * general_ams;
+
+   Solver * smoother;
+   ParDiscreteLinearOperator * pa_grad;
+   OperatorPtr Gradient;
+   ParDiscreteLinearOperator * pa_interp;
+   OperatorPtr Pi;
+
+   Solver * Gspacesolver;
+   Solver * Pispacesolver;
+
+   ParFiniteElementSpace * h1_fespace;
+   ParFiniteElementSpace * h1_fespace_d;
+
+   Array<int> offsets_nd;
+   Array<int> offsets_vector;
+   Array<int> offsets_scalar;
+
+   BlockOperator *Pi_block;
+   BlockOperator *Gradient_block;
+   BlockOperator *smoother_block;
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif

linalg/linalg.hpp

@@ -29,6 +29,7 @@
 #include "solvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
+#include "auxiliary.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/operator.hpp

@@ -696,6 +696,42 @@ public:
    { A_.Mult(x, y); y *= a_; }
 };
 
+/// General sum operator: x -> A(x)+B(x)
+class SumOperator : public Operator
+{
+   const Operator *A, *B;
+   bool ownA, ownB;
+   mutable Vector z, w;
+   double cA, cB;
+
+public:
+   SumOperator(const Operator *A_, const Operator *B_,
+               bool ownA_, bool ownB_, double cA_, double cB_)
+      : Operator(A_->Height(), B_->Width()),
+        A(A_), B(B_), ownA(ownA_), ownB(ownB_), z(A_->Height()), w(A_->Width()),
+        cA(cA_), cB(cB_)
+   {
+      MFEM_VERIFY(A->Width() == B->Width() && A->Height() == B->Height(),
+                  "incompatible Operators: A->Width() = " << A->Width()
+                  << ", B->Height() = " << B->Height());
+
+      z.UseDevice(true);
+      w.UseDevice(true);
+   }
+
+   ~SumOperator()
+   {
+      if (ownA) { delete A; }
+      if (ownB) { delete B; }
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   { B->Mult(x, z); A->Mult(x, y); y *= cA; z *= cB; y += z;}
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { B->MultTranspose(x, w); A->MultTranspose(x, y); y *= cA; w *= cB; y += w;}
+
+};
 
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */

linalg/solvers.hpp

@@ -852,6 +852,169 @@ public:
 
 #endif // MFEM_USE_SUITESPARSE
 
+class Complex_PMHSS : public Solver
+{
+public:
+   Complex_PMHSS(Operator *Re, Operator *Im, Solver *prec_Re, Solver *prec_Im,
+                 double a_)
+      : Solver(2*Re->Height()), a(a_), A(Re, Im, false, false),
+        A_Re(Re, NULL, false, false),
+        A_Im(Im, NULL, false, false), u(2*Re->Height()), rhs(2*Re->Height()),
+        n(Re->Height())
+   {
+      MFEM_VERIFY(Re->Height() == Im->Height() && Re->Height() == Re->Width() &&
+                  Im->Height() == Im->Width(), "");
+      MFEM_VERIFY(this->Height() == A.Height(), "");
+
+      // Create CG solver for real operator aV + A_Re in complex space.
+
+      V = useIdentityV ? (Operator*) new IdentityOperator(this->Height()) :
+          (Operator*) &A_Re;
+
+      // In the case V = A_Re, it is faster to use a scaled operator than a SumOperator
+      Operator *sumOpRe = useIdentityV ? (Operator*) new SumOperator(V, &A_Re, false,
+                                                                     false, a, 1.0)
+                          : (Operator*) new ScaledOperator(&A_Re, a + 1.0);
+
+      SumOperator *sumOpIm = new SumOperator(V, &A_Im, false, false, a, 1.0);
+
+      CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
+      cg->SetRelTol(1e-6);
+      cg->SetMaxIter(1000);
+      cg->SetPrintLevel(0);
+      cg->SetOperator(*sumOpRe);
+      cg->SetPreconditioner(*prec_Re);
+      cg->iterative_mode = false;
+
+      SRe = cg;
+
+      CGSolver *cgi = new CGSolver(MPI_COMM_WORLD);
+      cgi->SetRelTol(1e-6);
+      cgi->SetMaxIter(1000);
+      cgi->SetPrintLevel(0);
+      cgi->SetOperator(*sumOpIm);
+      if (prec_Im && useIdentityV) { cgi->SetPreconditioner(*prec_Im); }
+      if (!useIdentityV) { cgi->SetPreconditioner(*prec_Re); }
+      cgi->iterative_mode = false;
+
+      /*
+      // For negative definite imaginary part, but then PMHSS does not work?
+      MINRESSolver *cgi = new MINRESSolver(MPI_COMM_WORLD);
+      cgi->SetRelTol(1e-12);
+      cgi->SetMaxIter(1000);
+      cgi->SetPrintLevel(0);
+      cgi->SetOperator(*sumOpIm);
+      if (prec_Im) cgi->SetPreconditioner(*prec_Im);
+      */
+
+      SIm = cgi;
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      MFEM_VERIFY(false, "Don't call SetOperator");
+   }
+
+   void ComputeResidual(const Vector &b, const Vector &sol, Vector &res) const
+   {
+      A.Mult(sol, res);
+      res -= b;
+   }
+
+   void Mult(const Vector &x, Vector &y) const
+   {
+      if (!(x.Size() == Height() && y.Size() == Height()))
+      {
+         std::cout << "bug";
+      }
+
+      MFEM_VERIFY(x.Size() == Height() && y.Size() == Height(), "");
+
+      const double initNorm = x.Norml2();
+      mfem::out << "MHSS RHS norm " << initNorm << '\n';
+
+      // With V = I, use modified HSS (MHSS) from Bai, Benzi, Chen 2010.
+      y = 0.0;
+
+      for (int it=0; it<maxiter; ++it)
+      {
+         // Solve (aI + Re) u = (aI - i Im) y + x
+
+         if (it == 0)
+         {
+            // Optimize the first iteration, when the initial guess is y=0.
+            SRe->Mult(x, u);
+         }
+         else
+         {
+            A_Im.Mult(y, u);  // u = Im y
+            // Set rhs = -i Im y = -i u
+            for (int j=0; j<n; ++j)
+            {
+               rhs[j] = u[n+j];
+               rhs[n+j] = -u[j];
+            }
+
+            rhs += x;
+
+            V->Mult(y, u);
+            rhs.Add(a, u);
+
+            SRe->Mult(rhs, u);
+         }
+
+         // Solve (aI + Im) y = (aI + i Re) u - i x
+
+         A_Re.Mult(u, y);  // y = Re u
+         // Set rhs = i (Re u - x) = i (y - x)
+         for (int j=0; j<n; ++j)
+         {
+            rhs[j] = -(y[n+j] - x[n+j]);
+            rhs[n+j] = y[j] - x[j];
+         }
+
+         if (useIdentityV)
+         {
+            //V->Mult(u, y);
+            //rhs.Add(a, y);
+            rhs.Add(a, u);
+         }
+         else
+         {
+            // Using V = A_Re
+            rhs.Add(a, y);
+         }
+
+         SIm->Mult(rhs, y);
+
+         ComputeResidual(x, y, rhs);
+         const double resNorm = rhs.Norml2();
+         mfem::out << "MHSS iter " << it << " residual norm " << resNorm << '\n';
+
+         if (resNorm / initNorm < tol)
+         {
+            mfem::out << "MHSS converged\n";
+            break;
+         }
+      }
+   }
+
+private:
+   const double a;
+   const int maxiter = 1;
+   ComplexOperator A, A_Re, A_Im;
+   mutable Vector u, rhs;
+   const int n;
+
+   const double tol = 1.0e-8;
+
+   const bool useIdentityV = false;
+   Operator *V = NULL;
+
+   Solver *SRe = NULL;
+   Solver *SIm = NULL;
+};
+
 }
 
 #endif // MFEM_SOLVERS