Merge pull request #220 from awslabs/sjg/estimator-solve-dev

Add `"EstimatorMG"` option for estimation mass matrix linear solve

docs/src/config/solver.md

@@ -337,6 +337,9 @@ directory specified by [`config["Problem"]["Output"]`](problem.md#config%5B%22Pr
     "PCSide": <string>,
     "DivFreeTol": <float>,
     "DivFreeMaxIts": <float>,
+    "EstimatorTol": <float>,
+    "EstimatorMaxIts": <float>,
+    "EstimatorMG": <bool>,
     "GSOrthogonalization": <string>
 }
 ```
@@ -445,9 +448,12 @@ the eigenmode simulation type.
 `"EstimatorTol" [1e-6]` :  Relative tolerance for flux projection used in the
 error estimate calculation.
 
-`"EstimatorMaxIts" [1000]` :  Maximum number of iterations for flux projection use in the
+`"EstimatorMaxIts" [10000]` :  Maximum number of iterations for flux projection use in the
 error estimate calculation.
 
+`"EstimatorMG" [false]` :  Set to true in order to enable multigrid preconditioner with AMG
+coarse solve for the error estimate linear solver, instead of just Jacobi.
+
 `"GSOrthogonalization" ["MGS"]` :  Gram-Schmidt variant used to explicitly orthogonalize
 vectors in Krylov subspace methods or other parts of the code.
 

palace/drivers/drivensolver.cpp

@@ -142,8 +142,8 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
 
   // Initialize structures for storing and reducing the results of error estimation.
   CurlFluxErrorEstimator<ComplexVector> estimator(
-      spaceop.GetMaterialOp(), spaceop.GetNDSpace(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0);
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
 
   // Main frequency sweep loop.
@@ -242,8 +242,8 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
 
   // Initialize structures for storing and reducing the results of error estimation.
   CurlFluxErrorEstimator<ComplexVector> estimator(
-      spaceop.GetMaterialOp(), spaceop.GetNDSpace(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0);
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
 
   // Configure the PROM operator which performs the parameter space sampling and basis

palace/drivers/eigensolver.cpp

@@ -273,8 +273,8 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   // Calculate and record the error indicators.
   Mpi::Print("\nComputing solution error estimates\n");
   CurlFluxErrorEstimator<ComplexVector> estimator(
-      spaceop.GetMaterialOp(), spaceop.GetNDSpace(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0);
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
   if (!KM)
   {

palace/drivers/electrostaticsolver.cpp

@@ -44,6 +44,13 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   Vector RHS(K->Height());
   std::vector<Vector> V(nstep);
 
+  // Initialize structures for storing and reducing the results of error estimation.
+  GradFluxErrorEstimator estimator(
+      laplaceop.GetMaterialOp(), laplaceop.GetH1Space(), laplaceop.GetRTSpaces(),
+      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
+      iodata.solver.linear.estimator_mg);
+  ErrorIndicator indicator;
+
   // Main loop over terminal boundaries.
   Mpi::Print("\nComputing electrostatic fields for {:d} terminal boundar{}\n", nstep,
              (nstep > 1) ? "ies" : "y");
@@ -65,19 +72,24 @@ ElectrostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
                linalg::Norml2(laplaceop.GetComm(), V[step]),
                linalg::Norml2(laplaceop.GetComm(), RHS));
 
+    // Calculate and record the error indicators.
+    Mpi::Print(" Updating solution error estimates\n");
+    estimator.AddErrorIndicator(V[step], indicator);
+
     // Next terminal.
     step++;
   }
 
   // Postprocess the capacitance matrix from the computed field solutions.
   BlockTimer bt1(Timer::POSTPRO);
   SaveMetadata(ksp);
-  return {Postprocess(laplaceop, postop, V), laplaceop.GlobalTrueVSize()};
+  Postprocess(laplaceop, postop, V, indicator);
+  return {indicator, laplaceop.GlobalTrueVSize()};
 }
 
-ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
-                                                PostOperator &postop,
-                                                const std::vector<Vector> &V) const
+void ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop, PostOperator &postop,
+                                      const std::vector<Vector> &V,
+                                      const ErrorIndicator &indicator) const
 {
   // Postprocess the Maxwell capacitance matrix. See p. 97 of the COMSOL AC/DC Module manual
   // for the associated formulas based on the electric field energy based on a unit voltage
@@ -90,18 +102,6 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
   int nstep = static_cast<int>(terminal_sources.size());
   mfem::DenseMatrix C(nstep), Cm(nstep);
   Vector E(Grad.Height()), Vij(Grad.Width());
-
-  // Calculate and record the error indicators.
-  Mpi::Print("\nComputing solution error estimates\n");
-  GradFluxErrorEstimator estimator(laplaceop.GetMaterialOp(), laplaceop.GetH1Space(),
-                                   iodata.solver.linear.estimator_tol,
-                                   iodata.solver.linear.estimator_max_it, 0);
-  ErrorIndicator indicator;
-  for (int i = 0; i < nstep; i++)
-  {
-    estimator.AddErrorIndicator(V[i], indicator);
-  }
-
   int i = 0;
   for (const auto &[idx, data] : terminal_sources)
   {
@@ -145,7 +145,7 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
       else if (j > i)
       {
         linalg::AXPBYPCZ(1.0, V[i], 1.0, V[j], 0.0, Vij);
-        postop.SetVGridFunction(Vij);
+        postop.SetVGridFunction(Vij, false);
         double Ue = postop.GetEFieldEnergy();
         C(i, j) = Ue - 0.5 * (C(i, i) + C(j, j));
         Cm(i, j) = -C(i, j);
@@ -156,7 +156,6 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
   mfem::DenseMatrix Cinv(C);
   Cinv.Invert();  // In-place, uses LAPACK (when available) and should be cheap
   PostprocessTerminals(terminal_sources, C, Cinv, Cm);
-  return indicator;
 }
 
 void ElectrostaticSolver::PostprocessTerminals(

palace/drivers/electrostaticsolver.hpp

@@ -35,8 +35,8 @@ class Timer;
 class ElectrostaticSolver : public BaseSolver
 {
 private:
-  ErrorIndicator Postprocess(LaplaceOperator &laplaceop, PostOperator &postop,
-                             const std::vector<Vector> &V) const;
+  void Postprocess(LaplaceOperator &laplaceop, PostOperator &postop,
+                   const std::vector<Vector> &V, const ErrorIndicator &indicator) const;
 
   void PostprocessTerminals(const std::map<int, mfem::Array<int>> &terminal_sources,
                             const mfem::DenseMatrix &C, const mfem::DenseMatrix &Cinv,

palace/drivers/magnetostaticsolver.cpp

@@ -44,6 +44,13 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   Vector RHS(K->Height());
   std::vector<Vector> A(nstep);
 
+  // Initialize structures for storing and reducing the results of error estimation.
+  CurlFluxErrorEstimator<Vector> estimator(
+      curlcurlop.GetMaterialOp(), curlcurlop.GetNDSpaces(),
+      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
+      iodata.solver.linear.estimator_mg);
+  ErrorIndicator indicator;
+
   // Main loop over current source boundaries.
   Mpi::Print("\nComputing magnetostatic fields for {:d} source boundar{}\n", nstep,
              (nstep > 1) ? "ies" : "y");
@@ -67,19 +74,24 @@ MagnetostaticSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
                linalg::Norml2(curlcurlop.GetComm(), A[step]),
                linalg::Norml2(curlcurlop.GetComm(), RHS));
 
+    // Calculate and record the error indicators.
+    Mpi::Print(" Updating solution error estimates\n");
+    estimator.AddErrorIndicator(A[step], indicator);
+
     // Next source.
     step++;
   }
 
   // Postprocess the capacitance matrix from the computed field solutions.
   BlockTimer bt1(Timer::POSTPRO);
   SaveMetadata(ksp);
-  return {Postprocess(curlcurlop, postop, A), curlcurlop.GlobalTrueVSize()};
+  Postprocess(curlcurlop, postop, A, indicator);
+  return {indicator, curlcurlop.GlobalTrueVSize()};
 }
 
-ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
-                                                PostOperator &postop,
-                                                const std::vector<Vector> &A) const
+void MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop, PostOperator &postop,
+                                      const std::vector<Vector> &A,
+                                      const ErrorIndicator &indicator) const
 {
   // Postprocess the Maxwell inductance matrix. See p. 97 of the COMSOL AC/DC Module manual
   // for the associated formulas based on the magnetic field energy based on a current
@@ -94,18 +106,6 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
   mfem::DenseMatrix M(nstep), Mm(nstep);
   Vector B(Curl.Height()), Aij(Curl.Width());
   Vector Iinc(nstep);
-
-  // Calculate and record the error indicators.
-  Mpi::Print("\nComputing solution error estimates\n");
-  CurlFluxErrorEstimator<Vector> estimator(
-      curlcurlop.GetMaterialOp(), curlcurlop.GetNDSpace(),
-      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0);
-  ErrorIndicator indicator;
-  for (int i = 0; i < nstep; i++)
-  {
-    estimator.AddErrorIndicator(A[i], indicator);
-  }
-
   int i = 0;
   for (const auto &[idx, data] : surf_j_op)
   {
@@ -153,7 +153,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
       else if (j > i)
       {
         linalg::AXPBYPCZ(1.0, A[i], 1.0, A[j], 0.0, Aij);
-        postop.SetAGridFunction(Aij);
+        postop.SetAGridFunction(Aij, false);
         double Um = postop.GetHFieldEnergy();
         M(i, j) = Um / (Iinc(i) * Iinc(j)) -
                   0.5 * (M(i, i) * Iinc(i) / Iinc(j) + M(j, j) * Iinc(j) / Iinc(i));
@@ -165,7 +165,6 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
   mfem::DenseMatrix Minv(M);
   Minv.Invert();  // In-place, uses LAPACK (when available) and should be cheap
   PostprocessTerminals(surf_j_op, M, Minv, Mm);
-  return indicator;
 }
 
 void MagnetostaticSolver::PostprocessTerminals(const SurfaceCurrentOperator &surf_j_op,

palace/drivers/magnetostaticsolver.hpp

@@ -33,8 +33,8 @@ class Timer;
 class MagnetostaticSolver : public BaseSolver
 {
 private:
-  ErrorIndicator Postprocess(CurlCurlOperator &curlcurlop, PostOperator &postop,
-                             const std::vector<Vector> &A) const;
+  void Postprocess(CurlCurlOperator &curlcurlop, PostOperator &postop,
+                   const std::vector<Vector> &A, const ErrorIndicator &indicator) const;
 
   void PostprocessTerminals(const SurfaceCurrentOperator &surf_j_op,
                             const mfem::DenseMatrix &M, const mfem::DenseMatrix &Minv,

palace/drivers/transientsolver.cpp

@@ -79,9 +79,9 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   Mpi::Print("\n");
 
   // Initialize structures for storing and reducing the results of error estimation.
-  CurlFluxErrorEstimator<Vector> estimator(spaceop.GetMaterialOp(), spaceop.GetNDSpace(),
-                                           iodata.solver.linear.estimator_tol,
-                                           iodata.solver.linear.estimator_max_it, 0);
+  CurlFluxErrorEstimator<Vector> estimator(
+      spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.linear.estimator_mg);
   ErrorIndicator indicator;
 
   // Main time integration loop.

palace/linalg/amg.cpp

@@ -9,21 +9,24 @@ namespace palace
 BoomerAmgSolver::BoomerAmgSolver(int cycle_it, int smooth_it, int print)
   : mfem::HypreBoomerAMG()
 {
-  SetPrintLevel((print > 1) ? print - 1 : 0);
-  SetMaxIter(cycle_it);
-  SetTol(0.0);
+  HYPRE_BoomerAMGSetPrintLevel(*this, (print > 1) ? print - 1 : 0);
+  HYPRE_BoomerAMGSetMaxIter(*this, cycle_it);
+  HYPRE_BoomerAMGSetTol(*this, 0.0);
 
   // Set additional BoomerAMG options.
   int agg_levels = 1;  // Number of aggressive coarsening levels
   double theta = 0.5;  // AMG strength parameter = 0.25 is 2D optimal (0.5-0.8 for 3D)
+  int relax_type = 8;  // 8 = l1-symm. GS, 13 = l1-GS, 18 = l1-Jacobi, 16 = Chebyshev
   if (mfem::Device::Allows(mfem::Backend::DEVICE_MASK))
   {
     // Modify options for GPU-supported features.
     agg_levels = 0;
+    relax_type = 18;
   }
 
-  SetAggressiveCoarsening(agg_levels);
-  SetStrengthThresh(theta);
+  HYPRE_BoomerAMGSetAggNumLevels(*this, agg_levels);
+  HYPRE_BoomerAMGSetStrongThreshold(*this, theta);
+  HYPRE_BoomerAMGSetRelaxType(*this, relax_type);
   HYPRE_BoomerAMGSetNumSweeps(*this, smooth_it);
 
   // int coarse_relax_type = 8;  // l1-symm. GS (inexact coarse solve)

palace/linalg/chebyshev.cpp

@@ -13,16 +13,18 @@ namespace
 
 double GetLambdaMax(MPI_Comm comm, const Operator &A, const Vector &dinv)
 {
+  // Assumes A SPD (diag(A) > 0) to use Hermitian eigenvalue solver.
   DiagonalOperator Dinv(dinv);
   ProductOperator DinvA(Dinv, A);
-  return linalg::SpectralNorm(comm, DinvA, false);
+  return linalg::SpectralNorm(comm, DinvA, true);
 }
 
 double GetLambdaMax(MPI_Comm comm, const ComplexOperator &A, const ComplexVector &dinv)
 {
+  // Assumes A SPD (diag(A) > 0) to use Hermitian eigenvalue solver.
   ComplexDiagonalOperator Dinv(dinv);
   ComplexProductOperator DinvA(Dinv, A);
-  return linalg::SpectralNorm(comm, DinvA, false);
+  return linalg::SpectralNorm(comm, DinvA, A.IsReal());
 }
 
 template <bool Transpose = false>

palace/linalg/chebyshev.hpp

@@ -41,7 +41,7 @@ class ChebyshevSmoother : public Solver<OperType>
   double lambda_max, sf_max;
 
   // Temporary vector for smoother application.
-  mutable VecType d;
+  mutable VecType d, r;
 
 public:
   ChebyshevSmoother(MPI_Comm comm, int smooth_it, int poly_order, double sf_max);
@@ -50,14 +50,22 @@ class ChebyshevSmoother : public Solver<OperType>
 
   void Mult(const VecType &x, VecType &y) const override
   {
-    MFEM_ABORT("ChebyshevSmoother implements Mult2 using an additional preallocated "
-               "temporary vector!");
+    if (r.Size() != y.Size())
+    {
+      r.SetSize(y.Size());
+      r.UseDevice(true);
+    }
+    Mult2(x, y, r);
   }
 
   void MultTranspose(const VecType &x, VecType &y) const override
   {
-    MFEM_ABORT("ChebyshevSmoother implements MultTranspose2 using an additional "
-               "preallocated temporary vector!");
+    if (r.Size() != y.Size())
+    {
+      r.SetSize(y.Size());
+      r.UseDevice(true);
+    }
+    MultTranspose2(x, y, r);
   }
 
   void Mult2(const VecType &x, VecType &y, VecType &r) const override;
@@ -97,7 +105,7 @@ class ChebyshevSmoother1stKind : public Solver<OperType>
   double theta, delta, sf_max, sf_min;
 
   // Temporary vector for smoother application.
-  mutable VecType d;
+  mutable VecType d, r;
 
 public:
   ChebyshevSmoother1stKind(MPI_Comm comm, int smooth_it, int poly_order, double sf_max,
@@ -107,14 +115,22 @@ class ChebyshevSmoother1stKind : public Solver<OperType>
 
   void Mult(const VecType &x, VecType &y) const override
   {
-    MFEM_ABORT("ChebyshevSmoother1stKind implements Mult2 using an additional preallocated "
-               "temporary vector!");
+    if (r.Size() != y.Size())
+    {
+      r.SetSize(y.Size());
+      r.UseDevice(true);
+    }
+    Mult2(x, y, r);
   }
 
   void MultTranspose(const VecType &x, VecType &y) const override
   {
-    MFEM_ABORT("ChebyshevSmoother1stKind implements MultTranspose2 using an additional "
-               "preallocated temporary vector!");
+    if (r.Size() != y.Size())
+    {
+      r.SetSize(y.Size());
+      r.UseDevice(true);
+    }
+    MultTranspose2(x, y, r);
   }
 
   void Mult2(const VecType &x, VecType &y, VecType &r) const override;