Complex-valued operators for multigrid preconditioning

docs/src/config/solver.md

@@ -316,6 +316,7 @@ directory specified by [`config["Problem"]["Output"]`]
     "MGCycleIts": <int>,
     "MGSmoothIts": <int>,
     "MGSmoothOrder": <int>,
+    "PCMatReal": <bool>,
     "PCMatShifted": <bool>,
     "PCSide": <string>,
     "DivFreeTol": <float>,
@@ -394,13 +395,17 @@ multigrid preconditioners (when `"UseMultigrid"` is `true` or `"Type"` is `"AMS"
 `"MGSmoothIts" [1]` :  Number of pre- and post-smooth iterations used for multigrid
 preconditioners (when `"UseMultigrid"` is `true` or `"Type"` is `"AMS"` or `"BoomerAMG"`).
 
-`"MGSmoothOrder" [3]` :  Order of polynomial smoothing for geometric multigrid
+`"MGSmoothOrder" [4]` :  Order of polynomial smoothing for geometric multigrid
 preconditioning (when `"UseMultigrid"` is `true`).
 
+`"PCMatReal" [false]` :  When set to `true`, constructs the preconditioner for frequency
+domain problems using a real-valued approximation of the system matrix. This is always
+performed for the coarsest multigrid level regardless of the setting of `"PCMatReal"`.
+
 `"PCMatShifted" [false]` :  When set to `true`, constructs the preconditioner for frequency
-domain problems using a real SPD approximation of the system matrix, which can help
-performance at high frequencies (relative to the lowest nonzero eigenfrequencies of the
-model).
+domain problems using a positive definite approximation of the system matrix by flipping
+the sign for the mass matrix contribution, which can help performance at high frequencies
+(relative to the lowest nonzero eigenfrequencies of the model).
 
 `"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.
@@ -427,6 +432,9 @@ vectors in Krylov subspace methods or other parts of the code.
   - `"InitialGuess" [true]`
   - `"MGLegacyTransfer" [false]`
   - `"MGAuxiliarySmoother" [true]`
+  - `"MGSmoothEigScaleMax" [1.0]`
+  - `"MGSmoothEigScaleMin" [0.0]`
+  - `"MGSmoothChebyshev4th" [true]`
   - `"PCLowOrderRefined" [false]`
   - `"ColumnOrdering" ["Default"]` :  `"METIS"`, `"ParMETIS"`,`"Scotch"`, `"PTScotch"`,
     `"Default"`

palace/drivers/basesolver.cpp

@@ -82,20 +82,26 @@ BaseSolver::BaseSolver(const IoData &iodata, bool root, int size, int num_thread
   }
 }
 
-void BaseSolver::SaveMetadata(const mfem::ParFiniteElementSpace &fespace) const
+void BaseSolver::SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const
 {
   if (post_dir.length() == 0)
   {
     return;
   }
-  const HYPRE_BigInt ndof = fespace.GlobalTrueVSize();
+  const mfem::ParFiniteElementSpace &fespace = fespaces.GetFinestFESpace();
   HYPRE_BigInt ne = fespace.GetParMesh()->GetNE();
   Mpi::GlobalSum(1, &ne, fespace.GetComm());
+  std::vector<HYPRE_BigInt> ndofs(fespaces.GetNumLevels());
+  for (int l = 0; l < fespaces.GetNumLevels(); l++)
+  {
+    ndofs[l] = fespaces.GetFESpaceAtLevel(l).GlobalTrueVSize();
+  }
   if (root)
   {
     json meta = LoadMetadata(post_dir);
     meta["Problem"]["MeshElements"] = ne;
-    meta["Problem"]["DegreesOfFreedom"] = ndof;
+    meta["Problem"]["DegreesOfFreedom"] = ndofs.back();
+    meta["Problem"]["MultigridDegreesOfFreedom"] = ndofs;
     WriteMetadata(post_dir, meta);
   }
 }

palace/drivers/basesolver.hpp

@@ -12,7 +12,7 @@
 namespace mfem
 {
 
-class ParFiniteElementSpace;
+class ParFiniteElementSpaceHierarchy;
 class ParMesh;
 
 }  // namespace mfem
@@ -77,7 +77,7 @@ class BaseSolver
                      Timer &timer) const = 0;
 
   // These methods write different simulation metadata to a JSON file in post_dir.
-  void SaveMetadata(const mfem::ParFiniteElementSpace &fespace) const;
+  void SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const;
   template <typename SolverType>
   void SaveMetadata(const SolverType &ksp) const;
   void SaveMetadata(const Timer &timer) const;

palace/drivers/drivensolver.cpp

@@ -42,7 +42,7 @@ void DrivenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
                  "Reverting to uniform sweep!\n");
     adaptive = false;
   }
-  SaveMetadata(spaceop.GetNDSpace());
+  SaveMetadata(spaceop.GetNDSpaces());
 
   // Frequencies will be sampled uniformly in the frequency domain. Index sets are for
   // computing things like S-parameters in postprocessing.
@@ -115,11 +115,11 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
   // Because the Dirichlet BC is always homogenous, no special elimination is required on
   // the RHS. Assemble the linear system for the initial frequency (so we can call
   // KspSolver::SetOperators). Compute everything at the first frequency step.
-  auto K = spaceop.GetComplexStiffnessMatrix(Operator::DIAG_ONE);
-  auto C = spaceop.GetComplexDampingMatrix(Operator::DIAG_ZERO);
-  auto M = spaceop.GetComplexMassMatrix(Operator::DIAG_ZERO);
-  auto A2 = spaceop.GetComplexExtraSystemMatrix(omega0, Operator::DIAG_ZERO);
-  auto Curl = spaceop.GetComplexCurlMatrix();
+  auto K = spaceop.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
+  auto C = spaceop.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
+  auto M = spaceop.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
+  auto A2 = spaceop.GetExtraSystemMatrix<ComplexOperator>(omega0, Operator::DIAG_ZERO);
+  auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();
 
   // Set up the linear solver and set operators for the first frequency step. The
   // preconditioner for the complex linear system is constructed from a real approximation
@@ -154,7 +154,7 @@ void DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &postop, in
     if (step > step0)
     {
       // Update frequency-dependent excitation and operators.
-      A2 = spaceop.GetComplexExtraSystemMatrix(omega, Operator::DIAG_ZERO);
+      A2 = spaceop.GetExtraSystemMatrix<ComplexOperator>(omega, Operator::DIAG_ZERO);
       A = spaceop.GetSystemMatrix(std::complex<double>(1.0, 0.0), 1i * omega,
                                   std::complex<double>(-omega * omega, 0.0), K.get(),
                                   C.get(), M.get(), A2.get());
@@ -236,7 +236,7 @@ void DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator &postop, i
 
   // Allocate negative curl matrix for postprocessing the B-field and vectors for the
   // high-dimensional field solution.
-  auto Curl = spaceop.GetComplexCurlMatrix();
+  auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();
   ComplexVector E(Curl->Width()), B(Curl->Height());
   E = 0.0;
   B = 0.0;

palace/drivers/eigensolver.cpp

@@ -30,11 +30,11 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
   // computational range. The damping matrix may be nullptr.
   timer.Lap();
   SpaceOperator spaceop(iodata, mesh);
-  auto K = spaceop.GetComplexStiffnessMatrix(Operator::DIAG_ONE);
-  auto C = spaceop.GetComplexDampingMatrix(Operator::DIAG_ZERO);
-  auto M = spaceop.GetComplexMassMatrix(Operator::DIAG_ZERO);
-  auto Curl = spaceop.GetComplexCurlMatrix();
-  SaveMetadata(spaceop.GetNDSpace());
+  auto K = spaceop.GetStiffnessMatrix<ComplexOperator>(Operator::DIAG_ONE);
+  auto C = spaceop.GetDampingMatrix<ComplexOperator>(Operator::DIAG_ZERO);
+  auto M = spaceop.GetMassMatrix<ComplexOperator>(Operator::DIAG_ZERO);
+  auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();
+  SaveMetadata(spaceop.GetNDSpaces());
 
   // Configure objects for postprocessing.
   PostOperator postop(iodata, spaceop, "eigenmode");
@@ -190,7 +190,7 @@ void EigenSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
     eigen->SetInitialSpace(v0);  // Copies the vector
 
     // Debug
-    // auto Grad = spaceop.GetComplexGradMatrix();
+    // auto Grad = spaceop.GetGradMatrix<ComplexOperator>();
     // ComplexVector r0(Grad->Width());
     // Grad->MultTranspose(v0, r0);
     // r0.Print();

palace/drivers/electrostaticsolver.cpp

@@ -25,7 +25,7 @@ void ElectrostaticSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mes
   timer.Lap();
   LaplaceOperator laplaceop(iodata, mesh);
   auto K = laplaceop.GetStiffnessMatrix();
-  SaveMetadata(laplaceop.GetH1Space());
+  SaveMetadata(laplaceop.GetH1Spaces());
 
   // Set up the linear solver.
   KspSolver ksp(iodata, laplaceop.GetH1Spaces());

palace/drivers/magnetostaticsolver.cpp

@@ -25,7 +25,7 @@ void MagnetostaticSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mes
   timer.Lap();
   CurlCurlOperator curlcurlop(iodata, mesh);
   auto K = curlcurlop.GetStiffnessMatrix();
-  SaveMetadata(curlcurlop.GetNDSpace());
+  SaveMetadata(curlcurlop.GetNDSpaces());
 
   // Set up the linear solver.
   KspSolver ksp(iodata, curlcurlop.GetNDSpaces(), &curlcurlop.GetH1Spaces());

palace/drivers/transientsolver.cpp

@@ -37,7 +37,7 @@ void TransientSolver::Solve(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh,
     delta_t = std::min(delta_t, 0.95 * dt_max);
   }
   int nstep = GetNumSteps(0.0, iodata.solver.transient.max_t, delta_t);
-  SaveMetadata(spaceop.GetNDSpace());
+  SaveMetadata(spaceop.GetNDSpaces());
 
   // Time stepping is uniform in the time domain. Index sets are for computing things like
   // port voltages and currents in postprocessing.