Complex-valued coarse preconditioner

docs/src/config/solver.md

@@ -322,6 +322,7 @@ under the directory specified by
     "MGSmoothOrder": <int>,
     "PCMatReal": <bool>,
     "PCMatShifted": <bool>,
+    "ComplexCoarseSolve": <bool>,
     "PCSide": <string>,
     "DivFreeTol": <float>,
     "DivFreeMaxIts": <float>,
@@ -420,6 +421,9 @@ domain problems using a positive definite approximation of the system matrix by
 the sign for the mass matrix contribution, which can help performance at high frequencies
 (relative to the lowest nonzero eigenfrequencies of the model).
 
+`"ComplexCoarseSolve" [true]` : When set to `true`, the coarse-level solver uses the true
+complex-valued system matrix. When set to `false`, the real-valued approximation is used.
+
 `"PCSide" ["Default"]` :  Side for preconditioning. Not all options are available for all
 iterative solver choices, and the default choice depends on the iterative solver used.
 

palace/linalg/ksp.cpp

@@ -114,7 +114,7 @@ std::unique_ptr<IterativeSolver<OperType>> ConfigureKrylovSolver(MPI_Comm comm,
 }
 
 template <typename OperType, typename T, typename... U>
-auto MakeWrapperSolver(U &&...args)
+auto MakeWrapperSolver(bool complex_coarse, U &&...args)
 {
   // Sparse direct solver types copy the input matrix, so there is no need to save the
   // parallel assembled operator.
@@ -131,7 +131,7 @@ auto MakeWrapperSolver(U &&...args)
 #endif
                                     false);
   return std::make_unique<MfemWrapperSolver<OperType>>(
-      std::make_unique<T>(std::forward<U>(args)...), save_assembled);
+      std::make_unique<T>(std::forward<U>(args)...), save_assembled, complex_coarse);
 }
 
 template <typename OperType>
@@ -144,6 +144,7 @@ ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata,
   std::unique_ptr<Solver<OperType>> pc;
   const auto type = iodata.solver.linear.type;
   const int print = iodata.problem.verbose - 1;
+  const bool complex_coarse = iodata.solver.linear.complex_coarse_solve;
   switch (type)
   {
     case config::LinearSolverData::Type::AMS:
@@ -152,40 +153,42 @@ ConfigurePreconditionerSolver(MPI_Comm comm, const IoData &iodata,
       MFEM_VERIFY(aux_fespaces, "AMS solver relies on both primary space "
                                 "and auxiliary spaces for construction!");
       pc = MakeWrapperSolver<OperType, HypreAmsSolver>(
-          iodata, fespaces.GetNumLevels() > 1, fespaces.GetFESpaceAtLevel(0),
-          aux_fespaces->GetFESpaceAtLevel(0), print);
+          complex_coarse, iodata, fespaces.GetNumLevels() > 1,
+          fespaces.GetFESpaceAtLevel(0), aux_fespaces->GetFESpaceAtLevel(0), print);
       break;
     case config::LinearSolverData::Type::BOOMER_AMG:
-      pc = MakeWrapperSolver<OperType, BoomerAmgSolver>(iodata, fespaces.GetNumLevels() > 1,
-                                                        print);
+      pc = MakeWrapperSolver<OperType, BoomerAmgSolver>(complex_coarse, iodata,
+                                                        fespaces.GetNumLevels() > 1, print);
       break;
     case config::LinearSolverData::Type::SUPERLU:
 #if defined(MFEM_USE_SUPERLU)
-      pc = MakeWrapperSolver<OperType, SuperLUSolver>(comm, iodata, print);
+      pc = MakeWrapperSolver<OperType, SuperLUSolver>(complex_coarse, comm, iodata, print);
 #else
       MFEM_ABORT("Solver was not built with SuperLU_DIST support, please choose a "
                  "different solver!");
 #endif
       break;
     case config::LinearSolverData::Type::STRUMPACK:
 #if defined(MFEM_USE_STRUMPACK)
-      pc = MakeWrapperSolver<OperType, StrumpackSolver>(comm, iodata, print);
+      pc =
+          MakeWrapperSolver<OperType, StrumpackSolver>(complex_coarse, comm, iodata, print);
 #else
       MFEM_ABORT("Solver was not built with STRUMPACK support, please choose a "
                  "different solver!");
 #endif
       break;
     case config::LinearSolverData::Type::STRUMPACK_MP:
 #if defined(MFEM_USE_STRUMPACK)
-      pc = MakeWrapperSolver<OperType, StrumpackMixedPrecisionSolver>(comm, iodata, print);
+      pc = MakeWrapperSolver<OperType, StrumpackMixedPrecisionSolver>(complex_coarse, comm,
+                                                                      iodata, print);
 #else
       MFEM_ABORT("Solver was not built with STRUMPACK support, please choose a "
                  "different solver!");
 #endif
       break;
     case config::LinearSolverData::Type::MUMPS:
 #if defined(MFEM_USE_MUMPS)
-      pc = MakeWrapperSolver<OperType, MumpsSolver>(comm, iodata, print);
+      pc = MakeWrapperSolver<OperType, MumpsSolver>(complex_coarse, comm, iodata, print);
 #else
       MFEM_ABORT(
           "Solver was not built with MUMPS support, please choose a different solver!");

palace/linalg/solver.cpp

@@ -52,7 +52,34 @@ void MfemWrapperSolver<ComplexOperator>::SetOperator(const ComplexOperator &op)
   }
   if (hAr && hAi)
   {
-    A.reset(mfem::Add(1.0, *hAr, 1.0, *hAi));
+    if (complex_matrix)
+    {
+      // A = [Ar, -Ai]
+      //     [Ai,  Ar]
+      mfem::Array2D<const mfem::HypreParMatrix *> blocks(2, 2);
+      mfem::Array2D<double> block_coeffs(2, 2);
+      blocks(0, 0) = hAr;
+      blocks(0, 1) = hAi;
+      blocks(1, 0) = hAi;
+      blocks(1, 1) = hAr;
+      block_coeffs(0, 0) = 1.0;
+      block_coeffs(0, 1) = -1.0;
+      block_coeffs(1, 0) = 1.0;
+      block_coeffs(1, 1) = 1.0;
+      A.reset(mfem::HypreParMatrixFromBlocks(blocks, &block_coeffs));
+      idx1.SetSize(op.Width());
+      idx2.SetSize(op.Width());
+      for (int i = 0; i < op.Width(); i++)
+      {
+        idx1[i] = i;
+        idx2[i] = i + op.Width();
+      }
+    }
+    else
+    {
+      // A = Ar + Ai.
+      A.reset(mfem::Add(1.0, *hAr, 1.0, *hAi));
+    }
     if (PtAPr)
     {
       PtAPr->StealParallelAssemble();
@@ -101,13 +128,32 @@ template <>
 void MfemWrapperSolver<ComplexOperator>::Mult(const ComplexVector &x,
                                               ComplexVector &y) const
 {
-  mfem::Array<const Vector *> X(2);
-  mfem::Array<Vector *> Y(2);
-  X[0] = &x.Real();
-  X[1] = &x.Imag();
-  Y[0] = &y.Real();
-  Y[1] = &y.Imag();
-  pc->ArrayMult(X, Y);
+  if (pc->Height() == x.Size())
+  {
+    mfem::Array<const Vector *> X(2);
+    mfem::Array<Vector *> Y(2);
+    X[0] = &x.Real();
+    X[1] = &x.Imag();
+    Y[0] = &y.Real();
+    Y[1] = &y.Imag();
+    pc->ArrayMult(X, Y);
+  }
+  else
+  {
+    Vector X(2 * x.Size()), Y(2 * y.Size()), yr, yi;
+    X.UseDevice(true);
+    Y.UseDevice(true);
+    yr.UseDevice(true);
+    yi.UseDevice(true);
+    X.SetSubVector(idx1, x.Real());
+    X.SetSubVector(idx2, x.Imag());
+    pc->Mult(X, Y);
+    Y.ReadWrite();
+    yr.MakeRef(Y, 0, y.Size());
+    yi.MakeRef(Y, y.Size(), y.Size());
+    y.Real() = yr;
+    y.Imag() = yi;
+  }
 }
 
 }  // namespace palace

palace/linalg/solver.hpp

@@ -76,17 +76,25 @@ class MfemWrapperSolver : public Solver<OperType>
   // The actual mfem::Solver.
   std::unique_ptr<mfem::Solver> pc;
 
-  // Real-valued system matrix A = Ar + Ai in parallel assembled form.
+  // System matrix A in parallel assembled form.
   std::unique_ptr<mfem::HypreParMatrix> A;
 
   // Whether or not to save the parallel assembled matrix after calling
   // mfem::Solver::SetOperator (some solvers copy their input).
   bool save_assembled;
 
+  // Whether to use the exact complex-valued system matrix or the real-valued
+  // approximation A = Ar + Ai.
+  bool complex_matrix = true;
+
+  // Indices of real and imaginary parts of the complex system RHS/solution.
+  mfem::Array<int> idx1, idx2;
+
 public:
-  MfemWrapperSolver(std::unique_ptr<mfem::Solver> &&pc, bool save_assembled = true)
+  MfemWrapperSolver(std::unique_ptr<mfem::Solver> &&pc, bool save_assembled = true,
+                    bool complex_matrix = true)
     : Solver<OperType>(pc->iterative_mode), pc(std::move(pc)),
-      save_assembled(save_assembled)
+      save_assembled(save_assembled), complex_matrix(complex_matrix)
   {
   }
 

palace/utils/configfile.cpp

@@ -1779,6 +1779,7 @@ void LinearSolverData::SetUp(json &solver)
   // Preconditioner-specific options.
   pc_mat_real = linear->value("PCMatReal", pc_mat_real);
   pc_mat_shifted = linear->value("PCMatShifted", pc_mat_shifted);
+  complex_coarse_solve = linear->value("ComplexCoarseSolve", complex_coarse_solve);
   pc_side_type = linear->value("PCSide", pc_side_type);
   sym_fact_type = linear->value("ColumnOrdering", sym_fact_type);
   strumpack_compression_type =
@@ -1821,6 +1822,7 @@ void LinearSolverData::SetUp(json &solver)
 
   linear->erase("PCMatReal");
   linear->erase("PCMatShifted");
+  linear->erase("ComplexCoarseSolve");
   linear->erase("PCSide");
   linear->erase("ColumnOrdering");
   linear->erase("STRUMPACKCompressionType");
@@ -1865,6 +1867,7 @@ void LinearSolverData::SetUp(json &solver)
 
     std::cout << "PCMatReal: " << pc_mat_real << '\n';
     std::cout << "PCMatShifted: " << pc_mat_shifted << '\n';
+    std::cout << "ComplexCoarseSolve: " << complex_coarse_solve << '\n';
     std::cout << "PCSide: " << pc_side_type << '\n';
     std::cout << "ColumnOrdering: " << sym_fact_type << '\n';
     std::cout << "STRUMPACKCompressionType: " << strumpack_compression_type << '\n';

palace/utils/configfile.hpp

@@ -873,6 +873,10 @@ struct LinearSolverData
   // (makes the preconditoner matrix SPD).
   int pc_mat_shifted = -1;
 
+  // For frequency domain applications, use the complex-valued system matrix in the sparse
+  // direct solver.
+  bool complex_coarse_solve = true;
+
   // Choose left or right preconditioning.
   enum class SideType
   {

scripts/schema/config/solver.json

@@ -117,7 +117,8 @@
         "MGSmoothEigScaleMin": { "type": "number", "minimum": 0 },
         "MGSmoothChebyshev4th": { "type": "boolean" },
         "PCMatReal": { "type": "boolean" },
-        "PCMatShifted": { "type": "boolean" },
+          "PCMatShifted": { "type": "boolean" },
+          "ComplexCoarseSolve": {"type": "boolean"},
         "PCSide": { "type": "string" },
         "ColumnOrdering": { "type": "string" },
         "STRUMPACKCompressionType": { "type": "string" },