Merge remote-tracking branch 'origin/main' into hughcars/amr-solve-es…

…timate-mark-refine-dev

docs/src/config/solver.md

@@ -89,10 +89,6 @@ Thus, this object is only relevant for [`config["Problem"]["Type"]: "Magnetostat
 `"Linear"` :  Top-level object for configuring the linear solver employed by all simulation
 types.
 
-### Advanced solver options
-
-  - `"PartialAssemblyInterpolators" [true]`
-
 ## `solver["Eigenmode"]`
 
 ```json
@@ -413,8 +409,10 @@ 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" [4]` :  Order of polynomial smoothing for geometric multigrid
-preconditioning (when `"UseMultigrid"` is `true`).
+`"MGSmoothOrder" [0]` :  Order of polynomial smoothing for geometric multigrid
+preconditioning (when `"UseMultigrid"` is `true`). A value less than 1 defaults to twice
+the solution order given in [`config["Solver"]["Order"]`]
+(problem.md#config%5B%22Solver%22%5D) or 4, whichever is larger.
 
 `"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
@@ -454,7 +452,6 @@ vectors in Krylov subspace methods or other parts of the code.
 ### Advanced linear solver options
 
   - `"InitialGuess" [true]`
-  - `"MGLegacyTransfer" [false]`
   - `"MGAuxiliarySmoother" [true]`
   - `"MGSmoothEigScaleMax" [1.0]`
   - `"MGSmoothEigScaleMin" [0.0]`

palace/drivers/basesolver.cpp

@@ -10,7 +10,7 @@
 #include <nlohmann/json.hpp>
 #include "drivers/transientsolver.hpp"
 #include "fem/errorindicator.hpp"
-#include "linalg/errorestimator.hpp"
+#include "fem/fespace.hpp"
 #include "linalg/ksp.hpp"
 #include "models/domainpostoperator.hpp"
 #include "models/postoperator.hpp"
@@ -216,17 +216,17 @@ void BaseSolver::SolveEstimateMarkRefine(
   }
   Mpi::Print("\nFinal Error Indicator: {:.3e}, Size: {}\n", indicators.Norml2(comm), ntdof);
 }
-void BaseSolver::SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const
+void BaseSolver::SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const
 {
   if (post_dir.length() == 0)
   {
     return;
   }
-  const mfem::ParFiniteElementSpace &fespace = fespaces.GetFinestFESpace();
+  const auto &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++)
+  for (std::size_t l = 0; l < fespaces.GetNumLevels(); l++)
   {
     ndofs[l] = fespaces.GetFESpaceAtLevel(l).GlobalTrueVSize();
   }

palace/drivers/basesolver.hpp

@@ -12,7 +12,6 @@
 namespace mfem
 {
 
-class ParFiniteElementSpaceHierarchy;
 class ParMesh;
 
 }  // namespace mfem
@@ -21,6 +20,7 @@ namespace palace
 {
 
 class ErrorIndicator;
+class FiniteElementSpaceHierarchy;
 class IoData;
 class PostOperator;
 class Timer;
@@ -95,7 +95,7 @@ class BaseSolver
   void SolveEstimateMarkRefine(std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) const;
 
   // These methods write different simulation metadata to a JSON file in post_dir.
-  void SaveMetadata(const mfem::ParFiniteElementSpaceHierarchy &fespaces) const;
+  void SaveMetadata(const FiniteElementSpaceHierarchy &fespaces) const;
   template <typename SolverType>
   void SaveMetadata(const SolverType &ksp) const;
   void SaveMetadata(const Timer &timer) const;

palace/drivers/drivensolver.cpp

@@ -117,7 +117,7 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
   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>();
+  const auto &Curl = spaceop.GetCurlMatrix();
 
   // 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
@@ -133,15 +133,14 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
 
   // Set up RHS vector for the incident field at port boundaries, and the vector for the
   // first frequency step.
-  ComplexVector RHS(Curl->Width()), E(Curl->Width()), B(Curl->Height());
+  ComplexVector RHS(Curl.Width()), E(Curl.Width()), B(Curl.Height());
   E = 0.0;
   B = 0.0;
 
   // Initialize structures for storing and reducing the results of error estimation.
   CurlFluxErrorEstimator<ComplexVector> estimator(
       spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold,
-      iodata.solver.pa_discrete_interp);
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold);
   ErrorIndicator indicator;
 
   // Main frequency sweep loop.
@@ -177,7 +176,8 @@ ErrorIndicator DrivenSolver::SweepUniform(SpaceOperator &spaceop, PostOperator &
     // PostOperator for all postprocessing operations.
     BlockTimer bt2(Timer::POSTPRO);
     double E_elec = 0.0, E_mag = 0.0;
-    Curl->Mult(E, B);
+    Curl.Mult(E.Real(), B.Real());
+    Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);
@@ -245,16 +245,15 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
 
   // Allocate negative curl matrix for postprocessing the B-field and vectors for the
   // high-dimensional field solution.
-  auto Curl = spaceop.GetCurlMatrix<ComplexOperator>();
-  ComplexVector E(Curl->Width()), B(Curl->Height());
+  const auto &Curl = spaceop.GetCurlMatrix();
+  ComplexVector E(Curl.Width()), B(Curl.Height());
   E = 0.0;
   B = 0.0;
 
   // Initialize structures for storing and reducing the results of error estimation.
   CurlFluxErrorEstimator<ComplexVector> estimator(
       spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold,
-      iodata.solver.pa_discrete_interp);
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold);
   ErrorIndicator indicator;
 
   // Configure the PROM operator which performs the parameter space sampling and basis
@@ -338,7 +337,8 @@ ErrorIndicator DrivenSolver::SweepAdaptive(SpaceOperator &spaceop, PostOperator
     // PostOperator for all postprocessing operations.
     BlockTimer bt4(Timer::POSTPRO);
     double E_elec = 0.0, E_mag = 0.0;
-    Curl->Mult(E, B);
+    Curl.Mult(E.Real(), B.Real());
+    Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);

palace/drivers/eigensolver.cpp

@@ -35,12 +35,12 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
   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>();
+  const auto &Curl = spaceop.GetCurlMatrix();
   SaveMetadata(spaceop.GetNDSpaces());
 
   // Configure objects for postprocessing.
   PostOperator postop(iodata, spaceop, "eigenmode");
-  ComplexVector E(Curl->Width()), B(Curl->Height());
+  ComplexVector E(Curl.Width()), B(Curl.Height());
 
   // Define and configure the eigensolver to solve the eigenvalue problem:
   //         (K + λ C + λ² M) u = 0    or    K u = -λ² M u
@@ -166,7 +166,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
         spaceop.GetMaterialOp(), spaceop.GetNDSpace(), spaceop.GetH1Spaces(),
         spaceop.GetAuxBdrTDofLists(), iodata.solver.linear.divfree_tol,
         iodata.solver.linear.divfree_max_it, divfree_verbose,
-        iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp);
+        iodata.solver.pa_order_threshold);
     eigen->SetDivFreeProjector(*divfree);
   }
 
@@ -192,9 +192,10 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
     eigen->SetInitialSpace(v0);  // Copies the vector
 
     // Debug
-    // auto Grad = spaceop.GetGradMatrix<ComplexOperator>();
+    // const auto &Grad = spaceop.GetGradMatrix();
     // ComplexVector r0(Grad->Width());
-    // Grad->MultTranspose(v0, r0);
+    // Grad.MultTranspose(v0.Real(), r0.Real());
+    // Grad.MultTranspose(v0.Imag(), r0.Imag());
     // r0.Print();
   }
 
@@ -260,8 +261,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
   // Calculate and record the error indicators.
   CurlFluxErrorEstimator<ComplexVector> estimator(
       spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold,
-      iodata.solver.pa_discrete_interp);
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold);
   ErrorIndicator indicator;
   for (int i = 0; i < iodata.solver.eigenmode.n; i++)
   {
@@ -296,7 +296,8 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh) cons
     // Compute B = -1/(iω) ∇ x E on the true dofs, and set the internal GridFunctions in
     // PostOperator for all postprocessing operations.
     eigen->GetEigenvector(i, E);
-    Curl->Mult(E, B);
+    Curl.Mult(E.Real(), B.Real());
+    Curl.Mult(E.Imag(), B.Imag());
     B *= -1.0 / (1i * omega);
     postop.SetEGridFunction(E);
     postop.SetBGridFunction(B);

palace/drivers/electrostaticsolver.cpp

@@ -86,21 +86,21 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
   // charges from the prescribed voltage to get C directly as:
   //         Q_i = ∫ ρ dV = ∫ ∇ ⋅ (ε E) dV = ∫ (ε E) ⋅ n dS
   // and C_ij = Q_i/V_j. The energy formulation avoids having to locally integrate E = -∇V.
-  auto Grad = laplaceop.GetGradMatrix();
+  const auto &Grad = laplaceop.GetGradMatrix();
   const std::map<int, mfem::Array<int>> &terminal_sources = laplaceop.GetSources();
   int nstep = static_cast<int>(terminal_sources.size());
   mfem::DenseMatrix C(nstep), Cm(nstep);
-  Vector E(Grad->Height()), Vij(Grad->Width());
+  Vector E(Grad.Height()), Vij(Grad.Width());
   if (iodata.solver.electrostatic.n_post > 0)
   {
     Mpi::Print("\n");
   }
 
   // Calculate and record the error indicators.
-  GradFluxErrorEstimator estimator(
-      laplaceop.GetMaterialOp(), laplaceop.GetH1Spaces(),
-      iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
-      iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp);
+  GradFluxErrorEstimator estimator(laplaceop.GetMaterialOp(), laplaceop.GetH1Spaces(),
+                                   iodata.solver.linear.estimator_tol,
+                                   iodata.solver.linear.estimator_max_it, 0,
+                                   iodata.solver.pa_order_threshold);
   ErrorIndicator indicator;
   for (int i = 0; i < nstep; i++)
   {
@@ -113,7 +113,7 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
     // Compute E = -∇V on the true dofs, and set the internal GridFunctions in PostOperator
     // for all postprocessing operations.
     E = 0.0;
-    Grad->AddMult(V[i], E, -1.0);
+    Grad.AddMult(V[i], E, -1.0);
     postop.SetEGridFunction(E);
     postop.SetVGridFunction(V[i]);
     double Ue = postop.GetEFieldEnergy();
@@ -151,7 +151,7 @@ ErrorIndicator ElectrostaticSolver::Postprocess(LaplaceOperator &laplaceop,
       {
         linalg::AXPBYPCZ(1.0, V[i], 1.0, V[j], 0.0, Vij);
         E = 0.0;
-        Grad->AddMult(Vij, E, -1.0);
+        Grad.AddMult(Vij, E, -1.0);
         postop.SetEGridFunction(E);
         double Ue = postop.GetEFieldEnergy();
         C(i, j) = Ue - 0.5 * (C(i, i) + C(j, j));

palace/drivers/magnetostaticsolver.cpp

@@ -88,11 +88,11 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
   //                         Φ_i = ∫ B ⋅ n_j dS
   // and M_ij = Φ_i/I_j. The energy formulation avoids having to locally integrate B =
   // ∇ x A.
-  auto Curl = curlcurlop.GetCurlMatrix();
+  const auto &Curl = curlcurlop.GetCurlMatrix();
   const SurfaceCurrentOperator &surf_j_op = curlcurlop.GetSurfaceCurrentOp();
   int nstep = static_cast<int>(surf_j_op.Size());
   mfem::DenseMatrix M(nstep), Mm(nstep);
-  Vector B(Curl->Height()), Aij(Curl->Width());
+  Vector B(Curl.Height()), Aij(Curl.Width());
   Vector Iinc(nstep);
   if (iodata.solver.magnetostatic.n_post > 0)
   {
@@ -103,7 +103,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
   CurlFluxErrorEstimator<Vector> estimator(
       curlcurlop.GetMaterialOp(), curlcurlop.GetNDSpaces(),
       iodata.solver.linear.estimator_tol, iodata.solver.linear.estimator_max_it, 0,
-      iodata.solver.pa_order_threshold, iodata.solver.pa_discrete_interp);
+      iodata.solver.pa_order_threshold);
   ErrorIndicator indicator;
   for (int i = 0; i < nstep; i++)
   {
@@ -120,7 +120,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
 
     // Compute B = ∇ x A on the true dofs, and set the internal GridFunctions in
     // PostOperator for all postprocessing operations.
-    Curl->Mult(A[i], B);
+    Curl.Mult(A[i], B);
     postop.SetBGridFunction(B);
     postop.SetAGridFunction(A[i]);
     double Um = postop.GetHFieldEnergy();
@@ -157,7 +157,7 @@ ErrorIndicator MagnetostaticSolver::Postprocess(CurlCurlOperator &curlcurlop,
       else if (j > i)
       {
         linalg::AXPBYPCZ(1.0, A[i], 1.0, A[j], 0.0, Aij);
-        Curl->Mult(Aij, B);
+        Curl.Mult(Aij, B);
         postop.SetBGridFunction(B);
         double Um = postop.GetHFieldEnergy();
         M(i, j) = Um / (Iinc(i) * Iinc(j)) -

palace/drivers/transientsolver.cpp

@@ -80,8 +80,7 @@ TransientSolver::Solve(const std::vector<std::unique_ptr<mfem::ParMesh>> &mesh)
   // Initialize structures for storing and reducing the results of error estimation.
   CurlFluxErrorEstimator<Vector> estimator(
       spaceop.GetMaterialOp(), spaceop.GetNDSpaces(), iodata.solver.linear.estimator_tol,
-      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold,
-      iodata.solver.pa_discrete_interp);
+      iodata.solver.linear.estimator_max_it, 0, iodata.solver.pa_order_threshold);
   ErrorIndicator indicator;
 
   // Main time integration loop.