diff --git a/docs/src/config/solver.md b/docs/src/config/solver.md index d278c96c4..2f0d05de6 100644 --- a/docs/src/config/solver.md +++ b/docs/src/config/solver.md @@ -316,6 +316,7 @@ directory specified by [`config["Problem"]["Output"]`] "MGCycleIts": , "MGSmoothIts": , "MGSmoothOrder": , + "PCMatReal": , "PCMatShifted": , "PCSide": , "DivFreeTol": , @@ -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"` diff --git a/palace/drivers/basesolver.cpp b/palace/drivers/basesolver.cpp index 89fa8433c..6d277b5cb 100644 --- a/palace/drivers/basesolver.cpp +++ b/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 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); } } diff --git a/palace/drivers/basesolver.hpp b/palace/drivers/basesolver.hpp index 2bbbdfa03..ca3e4ed14 100644 --- a/palace/drivers/basesolver.hpp +++ b/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 void SaveMetadata(const SolverType &ksp) const; void SaveMetadata(const Timer &timer) const; diff --git a/palace/drivers/drivensolver.cpp b/palace/drivers/drivensolver.cpp index ae6d3b341..d2ffd07b4 100644 --- a/palace/drivers/drivensolver.cpp +++ b/palace/drivers/drivensolver.cpp @@ -42,7 +42,7 @@ void DrivenSolver::Solve(std::vector> &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(Operator::DIAG_ONE); + auto C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); + auto M = spaceop.GetMassMatrix(Operator::DIAG_ZERO); + auto A2 = spaceop.GetExtraSystemMatrix(omega0, Operator::DIAG_ZERO); + 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 @@ -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(omega, Operator::DIAG_ZERO); A = spaceop.GetSystemMatrix(std::complex(1.0, 0.0), 1i * omega, std::complex(-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(); ComplexVector E(Curl->Width()), B(Curl->Height()); E = 0.0; B = 0.0; diff --git a/palace/drivers/eigensolver.cpp b/palace/drivers/eigensolver.cpp index ce38f9022..2ac19f0ee 100644 --- a/palace/drivers/eigensolver.cpp +++ b/palace/drivers/eigensolver.cpp @@ -30,11 +30,11 @@ void EigenSolver::Solve(std::vector> &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(Operator::DIAG_ONE); + auto C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); + auto M = spaceop.GetMassMatrix(Operator::DIAG_ZERO); + auto Curl = spaceop.GetCurlMatrix(); + SaveMetadata(spaceop.GetNDSpaces()); // Configure objects for postprocessing. PostOperator postop(iodata, spaceop, "eigenmode"); @@ -190,7 +190,7 @@ void EigenSolver::Solve(std::vector> &mesh, eigen->SetInitialSpace(v0); // Copies the vector // Debug - // auto Grad = spaceop.GetComplexGradMatrix(); + // auto Grad = spaceop.GetGradMatrix(); // ComplexVector r0(Grad->Width()); // Grad->MultTranspose(v0, r0); // r0.Print(); diff --git a/palace/drivers/electrostaticsolver.cpp b/palace/drivers/electrostaticsolver.cpp index fd4888f0e..7a83f4752 100644 --- a/palace/drivers/electrostaticsolver.cpp +++ b/palace/drivers/electrostaticsolver.cpp @@ -25,7 +25,7 @@ void ElectrostaticSolver::Solve(std::vector> &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()); diff --git a/palace/drivers/magnetostaticsolver.cpp b/palace/drivers/magnetostaticsolver.cpp index 38be4ae47..0ed721468 100644 --- a/palace/drivers/magnetostaticsolver.cpp +++ b/palace/drivers/magnetostaticsolver.cpp @@ -25,7 +25,7 @@ void MagnetostaticSolver::Solve(std::vector> &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()); diff --git a/palace/drivers/transientsolver.cpp b/palace/drivers/transientsolver.cpp index 270197b91..0ee395ed2 100644 --- a/palace/drivers/transientsolver.cpp +++ b/palace/drivers/transientsolver.cpp @@ -37,7 +37,7 @@ void TransientSolver::Solve(std::vector> &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. diff --git a/palace/linalg/chebyshev.cpp b/palace/linalg/chebyshev.cpp index 446c652b5..36a8f8bec 100644 --- a/palace/linalg/chebyshev.cpp +++ b/palace/linalg/chebyshev.cpp @@ -3,7 +3,6 @@ #include "chebyshev.hpp" -#include #include #include "linalg/rap.hpp" @@ -20,27 +19,22 @@ void GetInverseDiagonal(const ParOperator &A, Vector &dinv) dinv.Reciprocal(); } -void GetInverseDiagonal(const ComplexParOperator &A, Vector &dinv) +void GetInverseDiagonal(const ComplexParOperator &A, ComplexVector &dinv) { - MFEM_VERIFY(A.HasReal() && !A.HasImag(), - "ComplexOperator for ChebyshevSmoother must be real-valued for now!"); + MFEM_VERIFY(A.HasReal() || A.HasImag(), + "Invalid zero ComplexOperator for ChebyshevSmoother!"); dinv.SetSize(A.Height()); - A.Real()->AssembleDiagonal(dinv); + dinv.SetSize(A.Height()); + dinv = 0.0; + if (A.HasReal()) + { + A.Real()->AssembleDiagonal(dinv.Real()); + } + if (A.HasImag()) + { + A.Imag()->AssembleDiagonal(dinv.Imag()); + } dinv.Reciprocal(); - // MFEM_VERIFY(A.HasReal() || A.HasImag(), - // "Invalid zero ComplexOperator for ChebyshevSmoother!"); - // dinv.SetSize(A.Height()); - // dinv.SetSize(A.Height()); - // dinv = 0.0; - // if (A.HasReal()) - // { - // A.Real()->AssembleDiagonal(dinv.Real()); - // } - // if (A.HasImag()) - // { - // A.Imag()->AssembleDiagonal(dinv.Imag()); - // } - // dinv.Reciprocal(); } double GetLambdaMax(MPI_Comm comm, const Operator &A, const Vector &dinv) @@ -50,48 +44,13 @@ double GetLambdaMax(MPI_Comm comm, const Operator &A, const Vector &dinv) return linalg::SpectralNorm(comm, DinvA, false); } -double GetLambdaMax(MPI_Comm comm, const ComplexOperator &A, const Vector &dinv) +double GetLambdaMax(MPI_Comm comm, const ComplexOperator &A, const ComplexVector &dinv) { - MFEM_VERIFY(A.HasReal() && !A.HasImag(), - "ComplexOperator for ChebyshevSmoother must be real-valued for now!"); - DiagonalOperator Dinv(dinv); - ProductOperator DinvA(Dinv, *A.Real()); + ComplexDiagonalOperator Dinv(dinv); + ComplexProductOperator DinvA(Dinv, A); return linalg::SpectralNorm(comm, DinvA, false); } -} // namespace - -template -ChebyshevSmoother::ChebyshevSmoother(int smooth_it, int poly_order) - : Solver(), pc_it(smooth_it), order(poly_order), A(nullptr) -{ -} - -template -void ChebyshevSmoother::SetOperator(const OperType &op) -{ - using ParOperType = - typename std::conditional::value, - ComplexParOperator, ParOperator>::type; - - A = &op; - r.SetSize(op.Height()); - d.SetSize(op.Height()); - - const auto *PtAP = dynamic_cast(&op); - MFEM_VERIFY(PtAP, - "ChebyshevSmoother requires a ParOperator or ComplexParOperator operator!"); - GetInverseDiagonal(*PtAP, dinv); - - // Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. See - // mfem::OperatorChebyshevSmoother or Adams et al., Parallel multigrid smoothing: - // polynomial versus Gauss-Seidel, JCP (2003). - lambda_max = 1.01 * GetLambdaMax(PtAP->GetComm(), *A, dinv); -} - -namespace -{ - template inline void ApplyOp(const Operator &A, const Vector &x, Vector &y) { @@ -137,31 +96,28 @@ inline void ApplyOrder0(double sr, const Vector &dinv, const Vector &r, Vector & const int N = d.Size(); const auto *DI = dinv.Read(); const auto *R = r.Read(); - auto *D = d.ReadWrite(); + auto *D = d.Write(); mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { D[i] = sr * DI[i] * R[i]; }); } template -inline void ApplyOrder0(const double sr, const Vector &dinv, const ComplexVector &r, +inline void ApplyOrder0(const double sr, const ComplexVector &dinv, const ComplexVector &r, ComplexVector &d) { const int N = dinv.Size(); - // const auto *DIR = dinv.Real().Read(); - // const auto *DII = dinv.Imag().Read(); - const auto *DIR = dinv.Read(); + const auto *DIR = dinv.Real().Read(); + const auto *DII = dinv.Imag().Read(); const auto *RR = r.Real().Read(); const auto *RI = r.Imag().Read(); - auto *DR = d.Real().ReadWrite(); - auto *DI = d.Imag().ReadWrite(); + auto *DR = d.Real().Write(); + auto *DI = d.Imag().Write(); if constexpr (!Transpose) { mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { - // DR[i] = sr * (DIR[i] * RR[i] - DII[i] * RI[i]); - // DI[i] = sr * (DII[i] * RR[i] + DIR[i] * RI[i]); - DR[i] = sr * DIR[i] * RR[i]; - DI[i] = sr * DIR[i] * RI[i]; + DR[i] = sr * (DIR[i] * RR[i] - DII[i] * RI[i]); + DI[i] = sr * (DII[i] * RR[i] + DIR[i] * RI[i]); }); } else @@ -169,10 +125,8 @@ inline void ApplyOrder0(const double sr, const Vector &dinv, const ComplexVector mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { - // DR[i] = sr * (DIR[i] * RR[i] + DII[i] * RI[i]); - // DI[i] = sr * (-DII[i] * RR[i] + DIR[i] * RI[i]); - DR[i] = sr * DIR[i] * RR[i]; - DI[i] = sr * DIR[i] * RI[i]; + DR[i] = sr * (DIR[i] * RR[i] + DII[i] * RI[i]); + DI[i] = sr * (-DII[i] * RR[i] + DIR[i] * RI[i]); }); } } @@ -189,13 +143,12 @@ inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv, } template -inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv, +inline void ApplyOrderK(const double sd, const double sr, const ComplexVector &dinv, const ComplexVector &r, ComplexVector &d) { const int N = dinv.Size(); - // const auto *DIR = dinv.Real().Read(); - // const auto *DII = dinv.Imag().Read(); - const auto *DIR = dinv.Read(); + const auto *DIR = dinv.Real().Read(); + const auto *DII = dinv.Imag().Read(); const auto *RR = r.Real().Read(); const auto *RI = r.Imag().Read(); auto *DR = d.Real().ReadWrite(); @@ -205,10 +158,8 @@ inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv, mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { - // DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] - DII[i] * RI[i]); - // DI[i] = sd * DI[i] + sr * (DII[i] * RR[i] + DIR[i] * RI[i]); - DR[i] = sd * DR[i] + sr * DIR[i] * RR[i]; - DI[i] = sd * DI[i] + sr * DIR[i] * RI[i]; + DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] - DII[i] * RI[i]); + DI[i] = sd * DI[i] + sr * (DII[i] * RR[i] + DIR[i] * RI[i]); }); } else @@ -216,16 +167,45 @@ inline void ApplyOrderK(const double sd, const double sr, const Vector &dinv, mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { - // DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] + DII[i] * RI[i]); - // DI[i] = sd * DI[i] + sr * (-DII[i] * RR[i] + DIR[i] * RI[i]); - DR[i] = sd * DR[i] + sr * DIR[i] * RR[i]; - DI[i] = sd * DI[i] + sr * DIR[i] * RI[i]; + DR[i] = sd * DR[i] + sr * (DIR[i] * RR[i] + DII[i] * RI[i]); + DI[i] = sd * DI[i] + sr * (-DII[i] * RR[i] + DIR[i] * RI[i]); }); } } } // namespace +template +ChebyshevSmoother::ChebyshevSmoother(int smooth_it, int poly_order, double sf_max) + : Solver(), pc_it(smooth_it), order(poly_order), A(nullptr), lambda_max(0.0), + sf_max(sf_max) +{ + MFEM_VERIFY(order > 0, "Polynomial order for Chebyshev smoothing must be positive!"); +} + +template +void ChebyshevSmoother::SetOperator(const OperType &op) +{ + using ParOperType = + typename std::conditional::value, + ComplexParOperator, ParOperator>::type; + + A = &op; + r.SetSize(op.Height()); + d.SetSize(op.Height()); + this->height = op.Height(); + this->width = op.Width(); + + const auto *PtAP = dynamic_cast(&op); + MFEM_VERIFY(PtAP, + "ChebyshevSmoother requires a ParOperator or ComplexParOperator operator!"); + GetInverseDiagonal(*PtAP, dinv); + + // Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. See + // mfem::OperatorChebyshevSmoother or Adams et al. (2003). + lambda_max = sf_max * GetLambdaMax(PtAP->GetComm(), *A, dinv); +} + template void ChebyshevSmoother::Mult(const VecType &x, VecType &y) const { @@ -258,7 +238,84 @@ void ChebyshevSmoother::Mult(const VecType &x, VecType &y) const } } +template +ChebyshevSmoother1stKind::ChebyshevSmoother1stKind(int smooth_it, int poly_order, + double sf_max, double sf_min) + : Solver(), pc_it(smooth_it), order(poly_order), A(nullptr), theta(0.0), + sf_max(sf_max), sf_min(sf_min) +{ + MFEM_VERIFY(order > 0, "Polynomial order for Chebyshev smoothing must be positive!"); +} + +template +void ChebyshevSmoother1stKind::SetOperator(const OperType &op) +{ + using ParOperType = + typename std::conditional::value, + ComplexParOperator, ParOperator>::type; + + A = &op; + r.SetSize(op.Height()); + d.SetSize(op.Height()); + this->height = op.Height(); + this->width = op.Width(); + + const auto *PtAP = dynamic_cast(&op); + MFEM_VERIFY( + PtAP, + "ChebyshevSmoother1stKind requires a ParOperator or ComplexParOperator operator!"); + GetInverseDiagonal(*PtAP, dinv); + + // Set up Chebyshev coefficients using the computed maximum eigenvalue estimate. The + // optimized estimate of lambda_min comes from (2.24) of Phillips and Fischer (2022). + if (sf_min <= 0.0) + { + sf_min = 1.69 / (std::pow(order, 1.68) + 2.11 * order + 1.98); + } + const double lambda_max = sf_max * GetLambdaMax(PtAP->GetComm(), *A, dinv); + const double lambda_min = sf_min * lambda_max; + theta = 0.5 * (lambda_max + lambda_min); + delta = 0.5 * (lambda_max - lambda_min); +} + +template +void ChebyshevSmoother1stKind::Mult(const VecType &x, VecType &y) const +{ + // Apply smoother: y = y + p(A) (x - A y) . + for (int it = 0; it < pc_it; it++) + { + if (this->initial_guess || it > 0) + { + ApplyOp(*A, y, r); + linalg::AXPBY(1.0, x, -1.0, r); + } + else + { + r = x; + y = 0.0; + } + + // 1th-kind Chebyshev smoother, from Phillips and Fischer or Adams. + ApplyOrder0(1.0 / theta, dinv, r, d); + double rhop = delta / theta; + for (int k = 1; k < order; k++) + { + y += d; + ApplyOp(*A, d, r, -1.0); + const double rho = 1.0 / (2.0 * theta / delta - rhop); + const double sd = rho * rhop; + const double sr = 2.0 * rho / delta; + ApplyOrderK(sd, sr, dinv, r, d); + rhop = rho; + } + y += d; + } +} + template class ChebyshevSmoother; template class ChebyshevSmoother; +template class ChebyshevSmoother1stKind; +template class ChebyshevSmoother1stKind; + } // namespace palace diff --git a/palace/linalg/chebyshev.hpp b/palace/linalg/chebyshev.hpp index 47a8de3e3..56df15a44 100644 --- a/palace/linalg/chebyshev.hpp +++ b/palace/linalg/chebyshev.hpp @@ -14,9 +14,10 @@ namespace palace // // Matrix-free diagonally-scaled Chebyshev smoothing. This is largely the same as // mfem::OperatorChebyshevSmoother allows a nonzero initial guess and uses alternative -// methods to estimate the largest eigenvalue. See also Phillips and Fischer, Optimal -// Chebyshev smoothers and one-sided V-cycles, arXiv:2210.03179v1 (2022) for reference on -// the 4th-kind Chebyshev polynomial smoother. +// methods to estimate the largest eigenvalue. We use a smoother based on Chebyshev +// polynomials of the 4th-kind as proposed in recent work. +// Reference: Phillips and Fischer, Optimal Chebyshev smoothers and one-sided V-cycles, +// arXiv:2210.03179v1 (2022). // template class ChebyshevSmoother : public Solver @@ -31,16 +32,57 @@ class ChebyshevSmoother : public Solver const OperType *A; // Inverse diagonal scaling of the operator (real-valued for now). - Vector dinv; + VecType dinv; // Maximum operator eigenvalue for Chebyshev polynomial smoothing. - double lambda_max; + double lambda_max, sf_max; // Temporary vectors for smoother application. mutable VecType r, d; public: - ChebyshevSmoother(int smooth_it, int poly_order); + ChebyshevSmoother(int smooth_it, int poly_order, double sf_max); + + void SetOperator(const OperType &op) override; + + void Mult(const VecType &x, VecType &y) const override; + + void MultTranspose(const VecType &x, VecType &y) const override + { + Mult(x, y); // Assumes operator symmetry + } +}; + +// +// Matrix-free diagonally-scaled Chebyshev smoothing using standard 1st-kind Chebyshev +// polynomials. +// Reference: Adams et al., Parallel multigrid smoothing: polynomial versus Gauss–Seidel, +// JCP (2003). +// +template +class ChebyshevSmoother1stKind : public Solver +{ + using VecType = typename Solver::VecType; + +private: + // Number of smoother iterations and polynomial order. + const int pc_it, order; + + // System matrix (not owned). + const OperType *A; + + // Inverse diagonal scaling of the operator (real-valued for now). + VecType dinv; + + // Parameters depending on maximum and minimum operator eigenvalue estimates for Chebyshev + // polynomial smoothing. + double theta, delta, sf_max, sf_min; + + // Temporary vectors for smoother application. + mutable VecType r, d; + +public: + ChebyshevSmoother1stKind(int smooth_it, int poly_order, double sf_max, double sf_min); void SetOperator(const OperType &op) override; diff --git a/palace/linalg/distrelaxation.cpp b/palace/linalg/distrelaxation.cpp index 2012630e9..6228ca791 100644 --- a/palace/linalg/distrelaxation.cpp +++ b/palace/linalg/distrelaxation.cpp @@ -4,7 +4,6 @@ #include "distrelaxation.hpp" #include -#include #include "fem/multigrid.hpp" #include "linalg/chebyshev.hpp" #include "linalg/rap.hpp" @@ -15,7 +14,8 @@ namespace palace template DistRelaxationSmoother::DistRelaxationSmoother( mfem::ParFiniteElementSpace &nd_fespace, mfem::ParFiniteElementSpace &h1_fespace, - int smooth_it, int cheby_smooth_it, int cheby_order, int pa_order_threshold) + int smooth_it, int cheby_smooth_it, int cheby_order, double cheby_sf_max, + double cheby_sf_min, bool cheby_4th_kind, int pa_order_threshold) : Solver(), pc_it(smooth_it), A(nullptr), A_G(nullptr), dbc_tdof_list_G(nullptr) { // Construct discrete gradient matrix for the auxiliary space. @@ -29,8 +29,20 @@ DistRelaxationSmoother::DistRelaxationSmoother( } // Initialize smoothers. - B = std::make_unique>(cheby_smooth_it, cheby_order); - B_G = std::make_unique>(cheby_smooth_it, cheby_order); + if (cheby_4th_kind) + { + B = std::make_unique>(cheby_smooth_it, cheby_order, + cheby_sf_max); + B_G = std::make_unique>(cheby_smooth_it, cheby_order, + cheby_sf_max); + } + else + { + B = std::make_unique>(cheby_smooth_it, cheby_order, + cheby_sf_max, cheby_sf_min); + B_G = std::make_unique>(cheby_smooth_it, cheby_order, + cheby_sf_max, cheby_sf_min); + } B_G->SetInitialGuess(false); } @@ -47,16 +59,17 @@ void DistRelaxationSmoother::SetOperators(const OperType &op, "Invalid operator sizes for DistRelaxationSmoother!"); A = &op; A_G = &op_G; + r.SetSize(op.Height()); + x_G.SetSize(op_G.Height()); + y_G.SetSize(op_G.Height()); + this->height = op.Height(); + this->width = op.Width(); const auto *PtAP_G = dynamic_cast(&op_G); MFEM_VERIFY(PtAP_G, "ChebyshevSmoother requires a ParOperator or ComplexParOperator operator!"); dbc_tdof_list_G = PtAP_G->GetEssentialTrueDofs(); - r.SetSize(op.Height()); - x_G.SetSize(op_G.Height()); - y_G.SetSize(op_G.Height()); - // Set up smoothers for A and A_G. B->SetOperator(op); B_G->SetOperator(op_G); diff --git a/palace/linalg/distrelaxation.hpp b/palace/linalg/distrelaxation.hpp index 7c00f86d9..d988bdf20 100644 --- a/palace/linalg/distrelaxation.hpp +++ b/palace/linalg/distrelaxation.hpp @@ -53,7 +53,8 @@ class DistRelaxationSmoother : public Solver public: DistRelaxationSmoother(mfem::ParFiniteElementSpace &nd_fespace, mfem::ParFiniteElementSpace &h1_fespace, int smooth_it, - int cheby_smooth_it, int cheby_order, int pa_order_threshold); + int cheby_smooth_it, int cheby_order, double cheby_sf_max, + double cheby_sf_min, bool cheby_4th_kind, int pa_order_threshold); void SetOperator(const OperType &op) override { diff --git a/palace/linalg/divfree.cpp b/palace/linalg/divfree.cpp index 3453a8a0e..c1def4a56 100644 --- a/palace/linalg/divfree.cpp +++ b/palace/linalg/divfree.cpp @@ -68,7 +68,7 @@ DivFreeSolver::DivFreeSolver(const MaterialOperator &mat_op, auto amg = std::make_unique>(std::make_unique(1, 1, 0)); auto gmg = std::make_unique>( - std::move(amg), h1_fespaces, nullptr, 1, 1, 2, pa_order_threshold); + std::move(amg), h1_fespaces, nullptr, 1, 1, 2, 1.0, 0.0, true, pa_order_threshold); auto pcg = std::make_unique>(h1_fespaces.GetFinestFESpace().GetComm(), print); diff --git a/palace/linalg/gmg.cpp b/palace/linalg/gmg.cpp index 3bca69dab..b2ad28151 100644 --- a/palace/linalg/gmg.cpp +++ b/palace/linalg/gmg.cpp @@ -16,7 +16,8 @@ GeometricMultigridSolver::GeometricMultigridSolver( std::unique_ptr> &&coarse_solver, mfem::ParFiniteElementSpaceHierarchy &fespaces, mfem::ParFiniteElementSpaceHierarchy *aux_fespaces, int cycle_it, int smooth_it, - int cheby_order, int pa_order_threshold) + int cheby_order, double cheby_sf_max, double cheby_sf_min, bool cheby_4th_kind, + int pa_order_threshold) : Solver(), pc_it(cycle_it), A(fespaces.GetNumLevels()), P(fespaces.GetNumLevels() - 1), dbc_tdof_lists(fespaces.GetNumLevels() - 1), B(fespaces.GetNumLevels()), X(fespaces.GetNumLevels()), Y(fespaces.GetNumLevels()), @@ -43,13 +44,25 @@ GeometricMultigridSolver::GeometricMultigridSolver( { if (aux_fespaces) { + const int cheby_smooth_it = 1; B[l] = std::make_unique>( - fespaces.GetFESpaceAtLevel(l), aux_fespaces->GetFESpaceAtLevel(l), smooth_it, 1, - cheby_order, pa_order_threshold); + fespaces.GetFESpaceAtLevel(l), aux_fespaces->GetFESpaceAtLevel(l), smooth_it, + cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min, cheby_4th_kind, + pa_order_threshold); } else { - B[l] = std::make_unique>(smooth_it, cheby_order); + const int cheby_smooth_it = smooth_it; + if (cheby_4th_kind) + { + B[l] = std::make_unique>(cheby_smooth_it, cheby_order, + cheby_sf_max); + } + else + { + B[l] = std::make_unique>( + cheby_smooth_it, cheby_order, cheby_sf_max, cheby_sf_min); + } } } } @@ -106,6 +119,8 @@ void GeometricMultigridSolver::SetOperator(const OperType &op) Y[l].SetSize(A[l]->Height()); R[l].SetSize(A[l]->Height()); } + this->height = op.Height(); + this->width = op.Width(); } template diff --git a/palace/linalg/gmg.hpp b/palace/linalg/gmg.hpp index 1705594d5..0134ddd08 100644 --- a/palace/linalg/gmg.hpp +++ b/palace/linalg/gmg.hpp @@ -56,7 +56,9 @@ class GeometricMultigridSolver : public Solver GeometricMultigridSolver(std::unique_ptr> &&coarse_solver, mfem::ParFiniteElementSpaceHierarchy &fespaces, mfem::ParFiniteElementSpaceHierarchy *aux_fespaces, int cycle_it, - int smooth_it, int cheby_order, int pa_order_threshold); + int smooth_it, int cheby_order, double cheby_sf_max, + double cheby_sf_min, bool cheby_4th_kind, + int pa_order_threshold); GeometricMultigridSolver(const IoData &iodata, std::unique_ptr> &&coarse_solver, mfem::ParFiniteElementSpaceHierarchy &fespaces, @@ -64,7 +66,9 @@ class GeometricMultigridSolver : public Solver : GeometricMultigridSolver( std::move(coarse_solver), fespaces, aux_fespaces, iodata.solver.linear.mg_cycle_it, iodata.solver.linear.mg_smooth_it, - iodata.solver.linear.mg_smooth_order, iodata.solver.pa_order_threshold) + iodata.solver.linear.mg_smooth_order, iodata.solver.linear.mg_smooth_sf_max, + iodata.solver.linear.mg_smooth_sf_min, iodata.solver.linear.mg_smooth_cheby_4th, + iodata.solver.pa_order_threshold) { } diff --git a/palace/linalg/hcurl.cpp b/palace/linalg/hcurl.cpp index e4a45972f..e42e46127 100644 --- a/palace/linalg/hcurl.cpp +++ b/palace/linalg/hcurl.cpp @@ -79,7 +79,8 @@ WeightedHCurlNormSolver::WeightedHCurlNormSolver( nd_fespaces.GetFESpaceAtLevel(0), h1_fespaces.GetFESpaceAtLevel(0), 1, 1, 1, false, false, 0)); auto gmg = std::make_unique>( - std::move(ams), nd_fespaces, &h1_fespaces, 1, 1, 2, pa_order_threshold); + std::move(ams), nd_fespaces, &h1_fespaces, 1, 1, 2, 1.0, 0.0, true, + pa_order_threshold); auto pcg = std::make_unique>(nd_fespaces.GetFinestFESpace().GetComm(), print); diff --git a/palace/linalg/iterative.cpp b/palace/linalg/iterative.cpp index ccb6ff529..531af33f9 100644 --- a/palace/linalg/iterative.cpp +++ b/palace/linalg/iterative.cpp @@ -368,7 +368,7 @@ void CgSolver::Mult(const VecType &b, VecType &x) const Mpi::Print(comm, "{}{:{}d} KSP residual norm ||r||_B = {:.6e}\n", std::string(tab_width, ' '), it, int_width, res); } - if (print_opts.summary || (print_opts.warnings && !converged)) + if (print_opts.summary || (print_opts.warnings && eps > 0.0 && !converged)) { Mpi::Print(comm, "{}PCG solver {} in {:d} iteration{}", std::string(tab_width, ' '), converged ? "converged" : "did NOT converge", it, (it == 1) ? "" : "s"); @@ -650,7 +650,7 @@ void GmresSolver::Mult(const VecType &b, VecType &x) const Mpi::Print(comm, "{}{:{}d} (restart {:d}) KSP residual norm {:.6e}\n", std::string(tab_width, ' '), it, int_width, restart, beta); } - if (print_opts.summary || (print_opts.warnings && !converged)) + if (print_opts.summary || (print_opts.warnings && eps > 0.0 && !converged)) { Mpi::Print(comm, "{}GMRES solver {} in {:d} iteration{}", std::string(tab_width, ' '), converged ? "converged" : "did NOT converge", it, (it == 1) ? "" : "s"); @@ -814,7 +814,7 @@ void FgmresSolver::Mult(const VecType &b, VecType &x) const Mpi::Print(comm, "{}{:{}d} (restart {:d}) KSP residual norm {:.6e}\n", std::string(tab_width, ' '), it, int_width, restart, beta); } - if (print_opts.summary || (print_opts.warnings && !converged)) + if (print_opts.summary || (print_opts.warnings && eps > 0.0 && !converged)) { Mpi::Print(comm, "{}FGMRES solver {} in {:d} iteration{}", std::string(tab_width, ' '), converged ? "converged" : "did NOT converge", it, (it == 1) ? "" : "s"); diff --git a/palace/linalg/iterative.hpp b/palace/linalg/iterative.hpp index 1ab2b8c75..536632ebb 100644 --- a/palace/linalg/iterative.hpp +++ b/palace/linalg/iterative.hpp @@ -74,13 +74,18 @@ class IterativeSolver : public Solver } // Set the operator for the solver. - void SetOperator(const OperType &op) override { A = &op; } + void SetOperator(const OperType &op) override + { + A = &op; + this->height = op.Height(); + this->width = op.Width(); + } // Set the preconditioner for the solver. void SetPreconditioner(const Solver &pc) { B = &pc; } // Returns if the previous solve converged or not. - bool GetConverged() const { return converged; } + bool GetConverged() const { return converged && (rel_tol > 0.0 || abs_tol > 0.0); } // Returns the initial (absolute) residual for the previous solve. double GetInitialRes() const { return initial_res; } diff --git a/palace/linalg/jacobi.cpp b/palace/linalg/jacobi.cpp index f6f003f1a..9230eea87 100644 --- a/palace/linalg/jacobi.cpp +++ b/palace/linalg/jacobi.cpp @@ -4,28 +4,104 @@ #include "jacobi.hpp" #include +#include "linalg/rap.hpp" namespace palace { -void JacobiSmoother::SetOperator(const Operator &op) +namespace { - height = op.Height(); - width = op.Width(); - dinv.SetSize(height); - op.AssembleDiagonal(dinv); + +void GetInverseDiagonal(const ParOperator &A, Vector &dinv) +{ + dinv.SetSize(A.Height()); + A.AssembleDiagonal(dinv); + dinv.Reciprocal(); +} + +void GetInverseDiagonal(const ComplexParOperator &A, ComplexVector &dinv) +{ + MFEM_VERIFY(A.HasReal() || A.HasImag(), + "Invalid zero ComplexOperator for JacobiSmoother!"); + dinv.SetSize(A.Height()); + dinv = 0.0; + if (A.HasReal()) + { + A.Real()->AssembleDiagonal(dinv.Real()); + } + if (A.HasImag()) + { + A.Imag()->AssembleDiagonal(dinv.Imag()); + } dinv.Reciprocal(); } -void JacobiSmoother::Mult(const Vector &x, Vector &y) const +template +inline void Apply(const Vector &dinv, const Vector &x, Vector &y) { - MFEM_ASSERT(!iterative_mode, - "JacobiSmoother is not implemented for iterative_mode = true!"); - const int N = height; + const int N = dinv.Size(); const auto *DI = dinv.Read(); const auto *X = x.Read(); auto *Y = y.Write(); mfem::forall(N, [=] MFEM_HOST_DEVICE(int i) { Y[i] = DI[i] * X[i]; }); } +template +inline void Apply(const ComplexVector &dinv, const ComplexVector &x, ComplexVector &y) +{ + const int N = dinv.Size(); + const auto *DIR = dinv.Real().Read(); + const auto *DII = dinv.Imag().Read(); + const auto *XR = x.Real().Read(); + const auto *XI = x.Imag().Read(); + auto *YR = y.Real().Write(); + auto *YI = y.Imag().Write(); + if constexpr (!Transpose) + { + mfem::forall(N, + [=] MFEM_HOST_DEVICE(int i) + { + YR[i] = DIR[i] * XR[i] - DII[i] * XI[i]; + YI[i] = DII[i] * XR[i] + DIR[i] * XI[i]; + }); + } + else + { + mfem::forall(N, + [=] MFEM_HOST_DEVICE(int i) + { + YR[i] = DIR[i] * XR[i] + DII[i] * XI[i]; + YI[i] = -DII[i] * XR[i] + DIR[i] * XI[i]; + }); + } +} + +} // namespace + +template +void JacobiSmoother::SetOperator(const OperType &op) +{ + using ParOperType = + typename std::conditional::value, + ComplexParOperator, ParOperator>::type; + + this->height = op.Height(); + this->width = op.Width(); + + const auto *PtAP = dynamic_cast(&op); + MFEM_VERIFY(PtAP, + "JacobiSmoother requires a ParOperator or ComplexParOperator operator!"); + GetInverseDiagonal(*PtAP, dinv); +} + +template +void JacobiSmoother::Mult(const VecType &x, VecType &y) const +{ + MFEM_ASSERT(!this->initial_guess, "JacobiSmoother does not use initial guess!"); + Apply(dinv, x, y); +} + +template class JacobiSmoother; +template class JacobiSmoother; + } // namespace palace diff --git a/palace/linalg/jacobi.hpp b/palace/linalg/jacobi.hpp index 25e4735f0..3296c98d0 100644 --- a/palace/linalg/jacobi.hpp +++ b/palace/linalg/jacobi.hpp @@ -4,31 +4,34 @@ #ifndef PALACE_LINALG_JACOBI_SMOOTHER_HPP #define PALACE_LINALG_JACOBI_SMOOTHER_HPP -#include #include "linalg/operator.hpp" +#include "linalg/solver.hpp" #include "linalg/vector.hpp" namespace palace { // -// Simple Jacobi smoother using the diagonal vector from Operator::AssembleDiagonal(), +// Simple Jacobi smoother using the diagonal vector from OperType::AssembleDiagonal(), // which allows for (approximate) diagonal construction for matrix-free operators. // -class JacobiSmoother : public mfem::Solver +template +class JacobiSmoother : public Solver { + using VecType = typename Solver::VecType; + private: - // Inverse diagonal scaling of the operator. - Vector dinv; + // Inverse diagonal scaling of the operator (real-valued for now). + VecType dinv; public: - JacobiSmoother() : mfem::Solver() {} + JacobiSmoother() : Solver() {} - void SetOperator(const Operator &op) override; + void SetOperator(const OperType &op) override; - void Mult(const Vector &x, Vector &y) const override; + void Mult(const VecType &x, VecType &y) const override; - void MultTranspose(const Vector &x, Vector &y) const override { Mult(x, y); } + void MultTranspose(const VecType &x, VecType &y) const override { Mult(x, y); } }; } // namespace palace diff --git a/palace/linalg/solver.cpp b/palace/linalg/solver.cpp index c78a00c85..fcecbe39b 100644 --- a/palace/linalg/solver.cpp +++ b/palace/linalg/solver.cpp @@ -10,6 +10,8 @@ template <> void WrapperSolver::SetOperator(const Operator &op) { pc->SetOperator(op); + this->height = op.Height(); + this->width = op.Width(); } template <> @@ -19,6 +21,8 @@ void WrapperSolver::SetOperator(const ComplexOperator &op) "WrapperSolver::SetOperator requires an operator which is purely real for " "mfem::Solver!"); pc->SetOperator(*op.Real()); + this->height = op.Height(); + this->width = op.Width(); } template <> diff --git a/palace/linalg/solver.hpp b/palace/linalg/solver.hpp index 54f31a9a3..82755a1f2 100644 --- a/palace/linalg/solver.hpp +++ b/palace/linalg/solver.hpp @@ -19,7 +19,7 @@ namespace palace // Abstract base class for real-valued or complex-valued solvers. template -class Solver +class Solver : public OperType { static_assert(std::is_same::value || std::is_same::value, @@ -33,7 +33,7 @@ class Solver bool initial_guess; public: - Solver(bool initial_guess = false) : initial_guess(initial_guess) {} + Solver(bool initial_guess = false) : OperType(), initial_guess(initial_guess) {} virtual ~Solver() = default; // Configure whether or not to use an initial guess when applying the solver. @@ -42,11 +42,8 @@ class Solver // Set the operator associated with the solver, or update it if called repeatedly. virtual void SetOperator(const OperType &op) = 0; - // Apply the solver. - virtual void Mult(const VecType &x, VecType &y) const = 0; - // Apply the solver for the transpose problem. - virtual void MultTranspose(const VecType &x, VecType &y) const + void MultTranspose(const VecType &x, VecType &y) const override { MFEM_ABORT("MultTranspose() is not implemented for base class Solver!"); } diff --git a/palace/models/curlcurloperator.cpp b/palace/models/curlcurloperator.cpp index e1385417c..ba00d3a63 100644 --- a/palace/models/curlcurloperator.cpp +++ b/palace/models/curlcurloperator.cpp @@ -128,7 +128,8 @@ std::unique_ptr CurlCurlOperator::GetStiffnessMatrix() auto &nd_fespace_l = nd_fespaces.GetFESpaceAtLevel(l); if (print_hdr) { - Mpi::Print(" Level {:d}: {:d} unknowns", l, nd_fespace_l.GlobalTrueVSize()); + Mpi::Print(" Level {:d} (p = {:d}): {:d} unknowns", l, + nd_fespace_l.GetMaxElementOrder(), nd_fespace_l.GlobalTrueVSize()); } constexpr auto MatType = MaterialPropertyType::INV_PERMEABILITY; MaterialPropertyCoefficient muinv_func(mat_op); diff --git a/palace/models/laplaceoperator.cpp b/palace/models/laplaceoperator.cpp index 799c42b64..146f5cd94 100644 --- a/palace/models/laplaceoperator.cpp +++ b/palace/models/laplaceoperator.cpp @@ -150,7 +150,8 @@ std::unique_ptr LaplaceOperator::GetStiffnessMatrix() auto &h1_fespace_l = h1_fespaces.GetFESpaceAtLevel(l); if (print_hdr) { - Mpi::Print(" Level {:d}: {:d} unknowns", l, h1_fespace_l.GlobalTrueVSize()); + Mpi::Print(" Level {:d} (p = {:d}): {:d} unknowns", l, + h1_fespace_l.GetMaxElementOrder(), h1_fespace_l.GlobalTrueVSize()); } constexpr auto MatType = MaterialPropertyType::PERMITTIVITY_REAL; MaterialPropertyCoefficient epsilon_func(mat_op); diff --git a/palace/models/romoperator.cpp b/palace/models/romoperator.cpp index a41e184f2..37fbdf6c2 100644 --- a/palace/models/romoperator.cpp +++ b/palace/models/romoperator.cpp @@ -84,9 +84,9 @@ RomOperator::RomOperator(const IoData &iodata, SpaceOperator &spaceop) : spaceop // simply by setting diagonal entries of the system matrix for the corresponding dofs. // Because the Dirichlet BC is always homogenous, no special elimination is required on // the RHS. The damping matrix may be nullptr. - K = spaceop.GetComplexStiffnessMatrix(Operator::DIAG_ONE); - C = spaceop.GetComplexDampingMatrix(Operator::DIAG_ZERO); - M = spaceop.GetComplexMassMatrix(Operator::DIAG_ZERO); + K = spaceop.GetStiffnessMatrix(Operator::DIAG_ONE); + C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); + M = spaceop.GetMassMatrix(Operator::DIAG_ZERO); MFEM_VERIFY(K && M, "Invalid empty HDM matrices when constructing PROM!"); // Set up RHS vector (linear in frequency part) for the incident field at port boundaries, @@ -158,7 +158,7 @@ void RomOperator::SolveHDM(double omega, ComplexVector &e) { // Compute HDM solution at the given frequency. The system matrix, A = K + iω C - ω² M + // A2(ω) is built by summing the underlying operator contributions. - A2 = spaceop.GetComplexExtraSystemMatrix(omega, Operator::DIAG_ZERO); + A2 = spaceop.GetExtraSystemMatrix(omega, Operator::DIAG_ZERO); has_A2 = (A2 != nullptr); auto A = spaceop.GetSystemMatrix(std::complex(1.0, 0.0), 1i * omega, std::complex(-omega * omega, 0.0), K.get(), @@ -264,7 +264,7 @@ void RomOperator::AssemblePROM(double omega) // only nonzero on boundaries, will be empty if not needed. if (has_A2) { - A2 = spaceop.GetComplexExtraSystemMatrix(omega, Operator::DIAG_ZERO); + A2 = spaceop.GetExtraSystemMatrix(omega, Operator::DIAG_ZERO); ProjectMatInternal(spaceop.GetComm(), V, *A2, Ar, r, 0); } else @@ -368,40 +368,44 @@ double RomOperator::ComputeMaxError(int num_cand, double &omega_star) std::vector PC; if (Mpi::Root(spaceop.GetComm())) { - // Sample with uniform probability. - PC.reserve(num_cand); - std::sample(P_m_PS.begin(), P_m_PS.end(), std::back_inserter(PC), num_cand, engine); - -#if 0 - // Sample with weighted probability by distance from the set of already sampled - // points. - std::vector weights(P_m_PS.size()); - weights = static_cast(weights.Size()); - PC.reserve(num_cand); - for (auto sample : PS) + if constexpr (false) { - int i = std::distance(P_m_PS.begin(), P_m_PS.lower_bound(sample)); - int il = i-1; - while (il >= 0) + // Sample with weighted probability by distance from the set of already sampled + // points. + std::vector weights(P_m_PS.size()); + PC.reserve(num_cand); + for (auto sample : PS) { - weights[il] = std::min(weights[il], static_cast(i-il)); - il--; + int i = std::distance(P_m_PS.begin(), P_m_PS.lower_bound(sample)); + int il = i - 1; + while (il >= 0) + { + weights[il] = std::min(weights[il], static_cast(i - il)); + il--; + } + int iu = i; + while (iu < weights.size()) + { + weights[iu] = std::min(weights[iu], static_cast(1 + iu - i)); + iu++; + } } - int iu = i; - while (iu < weights.size()) + for (int i = 0; i < num_cand; i++) { - weights[iu] = std::min(weights[iu], static_cast(1+iu-i)); - iu++; + std::discrete_distribution dist(weights.begin(), weights.end()); + auto res = dist(engine); + auto it = P_m_PS.begin(); + std::advance(it, res); + PC.push_back(*it); + weights[res] = 0.0; // No replacement } } - for (int i = 0; i < num_cand; i++) + else { - std::discrete_distribution dist(weights.begin(), weights.end()); - int res = dist(engine); - PC.push_back(P_m_PS[res]); - weights[res] = 0.0; // No replacement + // Sample with uniform probability. + PC.reserve(num_cand); + std::sample(P_m_PS.begin(), P_m_PS.end(), std::back_inserter(PC), num_cand, engine); } -#endif } else { diff --git a/palace/models/spaceoperator.cpp b/palace/models/spaceoperator.cpp index d846c820c..22641ee56 100644 --- a/palace/models/spaceoperator.cpp +++ b/palace/models/spaceoperator.cpp @@ -73,12 +73,13 @@ mfem::Array SetUpBoundaryProperties(const IoData &iodata, const mfem::ParMe SpaceOperator::SpaceOperator(const IoData &iodata, const std::vector> &mesh) : pa_order_threshold(iodata.solver.pa_order_threshold), skip_zeros(0), - pc_lor(iodata.solver.linear.pc_mat_lor), - pc_shifted(iodata.solver.linear.pc_mat_shifted), print_hdr(true), print_prec_hdr(true), + pc_mat_real(iodata.solver.linear.pc_mat_real), + pc_mat_shifted(iodata.solver.linear.pc_mat_shifted), + pc_mat_lor(iodata.solver.linear.pc_mat_lor), print_hdr(true), print_prec_hdr(true), dbc_marker(SetUpBoundaryProperties(iodata, *mesh.back())), nd_fecs(utils::ConstructFECollections( iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels, - iodata.solver.linear.mg_coarsen_type, pc_lor)), + iodata.solver.linear.mg_coarsen_type, pc_mat_lor)), h1_fecs(utils::ConstructFECollections( iodata.solver.order, mesh.back()->Dimension(), iodata.solver.linear.mg_max_levels, iodata.solver.linear.mg_coarsen_type, false)), @@ -267,7 +268,8 @@ std::unique_ptr BuildAuxOperator(mfem::ParFiniteElementSpace &fespace, } // namespace -std::unique_ptr +template +std::unique_ptr SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy diag_policy) { PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); @@ -282,12 +284,22 @@ SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy diag_policy) auto k = BuildOperator(GetNDSpace(), &df, &f, (SumCoefficient *)nullptr, &fb, pa_order_threshold, skip_zeros); - auto K = std::make_unique(std::move(k), GetNDSpace()); - K->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return K; + if constexpr (std::is_same::value) + { + auto K = std::make_unique(std::move(k), nullptr, GetNDSpace()); + K->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return K; + } + else + { + auto K = std::make_unique(std::move(k), GetNDSpace()); + K->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return K; + } } -std::unique_ptr +template +std::unique_ptr SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy diag_policy) { PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); @@ -302,79 +314,32 @@ SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy diag_policy) auto c = BuildOperator(GetNDSpace(), (SumCoefficient *)nullptr, &f, (SumCoefficient *)nullptr, &fb, pa_order_threshold, skip_zeros); - auto C = std::make_unique(std::move(c), GetNDSpace()); - C->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return C; -} - -std::unique_ptr SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy diag_policy) -{ - PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); - const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); - SumMatrixCoefficient f(sdim), fb(sdim); - AddRealMassCoefficients(1.0, f); - AddRealMassBdrCoefficients(1.0, fb); - if (f.empty() && fb.empty()) - { - return {}; - } - - auto m = BuildOperator(GetNDSpace(), (SumCoefficient *)nullptr, &f, - (SumCoefficient *)nullptr, &fb, pa_order_threshold, skip_zeros); - auto M = std::make_unique(std::move(m), GetNDSpace()); - M->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return M; -} - -std::unique_ptr -SpaceOperator::GetComplexStiffnessMatrix(Operator::DiagonalPolicy diag_policy) -{ - PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); - const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); - SumMatrixCoefficient df(sdim), f(sdim), fb(sdim); - AddStiffnessCoefficients(1.0, df, f); - AddStiffnessBdrCoefficients(1.0, fb); - if (df.empty() && f.empty() && fb.empty()) + if constexpr (std::is_same::value) { - return {}; + auto C = std::make_unique(std::move(c), nullptr, GetNDSpace()); + C->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return C; } - - auto k = BuildOperator(GetNDSpace(), &df, &f, (SumCoefficient *)nullptr, &fb, - pa_order_threshold, skip_zeros); - auto K = std::make_unique(std::move(k), nullptr, GetNDSpace()); - K->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return K; -} - -std::unique_ptr -SpaceOperator::GetComplexDampingMatrix(Operator::DiagonalPolicy diag_policy) -{ - PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); - const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); - SumMatrixCoefficient f(sdim), fb(sdim); - AddDampingCoefficients(1.0, f); - AddDampingBdrCoefficients(1.0, fb); - if (f.empty() && fb.empty()) + else { - return {}; + auto C = std::make_unique(std::move(c), GetNDSpace()); + C->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return C; } - - auto c = BuildOperator(GetNDSpace(), (SumCoefficient *)nullptr, &f, - (SumCoefficient *)nullptr, &fb, pa_order_threshold, skip_zeros); - auto C = std::make_unique(std::move(c), nullptr, GetNDSpace()); - C->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return C; } -std::unique_ptr -SpaceOperator::GetComplexMassMatrix(Operator::DiagonalPolicy diag_policy) +template +std::unique_ptr SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy diag_policy) { PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); SumMatrixCoefficient fr(sdim), fi(sdim), fbr(sdim); AddRealMassCoefficients(1.0, fr); AddRealMassBdrCoefficients(1.0, fbr); - AddImagMassCoefficients(1.0, fi); + if constexpr (std::is_same::value) + { + AddImagMassCoefficients(1.0, fi); + } if (fr.empty() && fbr.empty() && fi.empty()) { return {}; @@ -392,14 +357,24 @@ SpaceOperator::GetComplexMassMatrix(Operator::DiagonalPolicy diag_policy) (SumCoefficient *)nullptr, (SumCoefficient *)nullptr, pa_order_threshold, skip_zeros); } - auto M = std::make_unique(std::move(mr), std::move(mi), GetNDSpace()); - M->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return M; + if constexpr (std::is_same::value) + { + auto M = + std::make_unique(std::move(mr), std::move(mi), GetNDSpace()); + M->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return M; + } + else + { + auto M = std::make_unique(std::move(mr), GetNDSpace()); + M->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return M; + } } -std::unique_ptr -SpaceOperator::GetComplexExtraSystemMatrix(double omega, - Operator::DiagonalPolicy diag_policy) +template +std::unique_ptr +SpaceOperator::GetExtraSystemMatrix(double omega, Operator::DiagonalPolicy diag_policy) { PrintHeader(GetH1Space(), GetNDSpace(), GetRTSpace(), pa_order_threshold, print_hdr); const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); @@ -422,9 +397,20 @@ SpaceOperator::GetComplexExtraSystemMatrix(double omega, ai = BuildOperator(GetNDSpace(), (SumCoefficient *)nullptr, (SumCoefficient *)nullptr, &dfbi, &fbi, pa_order_threshold, skip_zeros); } - auto A = std::make_unique(std::move(ar), std::move(ai), GetNDSpace()); - A->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); - return A; + if constexpr (std::is_same::value) + { + auto A = + std::make_unique(std::move(ar), std::move(ai), GetNDSpace()); + A->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return A; + } + else + { + MFEM_VERIFY(!ai, "Unexpected imaginary part in GetExtraSystemMatrix!"); + auto A = std::make_unique(std::move(ar), GetNDSpace()); + A->SetEssentialTrueDofs(nd_dbc_tdof_lists.back(), diag_policy); + return A; + } } namespace @@ -668,6 +654,8 @@ template std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, double a1, double a2, double a3) { + // XX TODO: Test complex PC matrix assembly for l == 0 if coarse solve supports it + // XX TODO: Handle complex coeff a0/a1/a2/a3 (like GetSystemMatrix) if (print_prec_hdr) { Mpi::Print("\nAssembling multigrid hierarchy:\n"); @@ -685,55 +673,55 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub auto &fespace_l = fespaces.GetFESpaceAtLevel(l); if (print_prec_hdr) { - Mpi::Print(" Level {:d}{}: {:d} unknowns", l, (s == 0) ? "" : " (auxiliary)", + Mpi::Print(" Level {:d}{} (p = {:d}): {:d} unknowns", l, + (s == 0) ? "" : " (auxiliary)", fespace_l.GetMaxElementOrder(), fespace_l.GlobalTrueVSize()); } const int sdim = GetNDSpace().GetParMesh()->SpaceDimension(); SumMatrixCoefficient dfr(sdim), fr(sdim), fi(sdim), fbr(sdim), fbi(sdim); SumCoefficient dfbr, dfbi; - // if (s > 0) - // { - - // // XX TODO: Test complex PC matrix assembly for s > 0 - // // (or s == 0 if coarse solve supports it) - // // XX TODO: Handle complex coeff a0/a1/a2 (like SumOperator) - - // AddStiffnessCoefficients(a0, dfr, fr); - // AddStiffnessBdrCoefficients(a0, fbr); - // AddDampingCoefficients(a1, fi); - // AddDampingBdrCoefficients(a1, fbi); - // AddRealMassCoefficients(pc_shifted ? std::abs(a2) : a2, fr); - // AddRealMassBdrCoefficients(pc_shifted ? std::abs(a2) : a2, fbr); - // AddImagMassCoefficients(a2, fi); - // AddExtraSystemBdrCoefficients(a3, dfbr, dfbi, fbr, fbi); - // } - // else + if (!std::is_same::value || pc_mat_real || l == 0) { + // Real-valued system matrix (approximation) for preconditioning. AddStiffnessCoefficients(a0, dfr, fr); AddStiffnessBdrCoefficients(a0, fbr); AddDampingCoefficients(a1, fr); AddDampingBdrCoefficients(a1, fbr); - AddAbsMassCoefficients(pc_shifted ? std::abs(a2) : a2, fr); - AddRealMassBdrCoefficients(pc_shifted ? std::abs(a2) : a2, fbr); + AddAbsMassCoefficients(pc_mat_shifted ? std::abs(a2) : a2, fr); + AddRealMassBdrCoefficients(pc_mat_shifted ? std::abs(a2) : a2, fbr); AddExtraSystemBdrCoefficients(a3, dfbr, dfbr, fbr, fbr); } + else + { + // Build preconditioner based on the actual complex-valued system matrix. + AddStiffnessCoefficients(a0, dfr, fr); + AddStiffnessBdrCoefficients(a0, fbr); + AddDampingCoefficients(a1, fi); + AddDampingBdrCoefficients(a1, fbi); + AddRealMassCoefficients(pc_mat_shifted ? std::abs(a2) : a2, fr); + AddRealMassBdrCoefficients(pc_mat_shifted ? std::abs(a2) : a2, fbr); + AddImagMassCoefficients(a2, fi); + AddExtraSystemBdrCoefficients(a3, dfbr, dfbi, fbr, fbi); + } std::unique_ptr br, bi; if (!dfr.empty() || !fr.empty() || !dfbr.empty() || !fbr.empty()) { - br = (s == 0) - ? BuildOperator(fespace_l, &dfr, &fr, &dfbr, &fbr, - (l > 0) ? pa_order_threshold : 100, skip_zeros, pc_lor) - : BuildAuxOperator(fespace_l, &fr, &fbr, - (l > 0) ? pa_order_threshold : 100, skip_zeros, pc_lor); + br = + (s == 0) + ? BuildOperator(fespace_l, &dfr, &fr, &dfbr, &fbr, + (l > 0) ? pa_order_threshold : 100, skip_zeros, pc_mat_lor) + : BuildAuxOperator(fespace_l, &fr, &fbr, (l > 0) ? pa_order_threshold : 100, + skip_zeros, pc_mat_lor); } if (!fi.empty() || !dfbi.empty() || !fbi.empty()) { - bi = (s == 0) - ? BuildOperator(fespace_l, (SumCoefficient *)nullptr, &fi, &dfbi, &fbi, - (l > 0) ? pa_order_threshold : 100, skip_zeros, pc_lor) - : BuildAuxOperator(fespace_l, &fi, &fbi, - (l > 0) ? pa_order_threshold : 100, skip_zeros, pc_lor); + bi = + (s == 0) + ? BuildOperator(fespace_l, (SumCoefficient *)nullptr, &fi, &dfbi, &fbi, + (l > 0) ? pa_order_threshold : 100, skip_zeros, pc_mat_lor) + : BuildAuxOperator(fespace_l, &fi, &fbi, (l > 0) ? pa_order_threshold : 100, + skip_zeros, pc_mat_lor); } if (print_prec_hdr) { @@ -741,7 +729,7 @@ std::unique_ptr SpaceOperator::GetPreconditionerMatrix(double a0, doub { HYPRE_BigInt nnz = br_spm->NumNonZeroElems(); Mpi::GlobalSum(1, &nnz, fespace_l.GetComm()); - Mpi::Print(", {:d} NNZ{}\n", nnz, pc_lor ? " (LOR)" : ""); + Mpi::Print(", {:d} NNZ{}\n", nnz, pc_mat_lor ? " (LOR)" : ""); } else { @@ -786,6 +774,7 @@ auto BuildGrad(mfem::ParFiniteElementSpace &h1_fespace, } // namespace +template <> std::unique_ptr SpaceOperator::GetCurlMatrix() { return std::make_unique( @@ -793,13 +782,15 @@ std::unique_ptr SpaceOperator::GetCurlMatrix() true); } -std::unique_ptr SpaceOperator::GetComplexCurlMatrix() +template <> +std::unique_ptr SpaceOperator::GetCurlMatrix() { return std::make_unique( BuildCurl(GetNDSpace(), GetRTSpace(), pa_order_threshold), nullptr, GetNDSpace(), GetRTSpace(), true); } +template <> std::unique_ptr SpaceOperator::GetGradMatrix() { return std::make_unique( @@ -807,7 +798,8 @@ std::unique_ptr SpaceOperator::GetGradMatrix() true); } -std::unique_ptr SpaceOperator::GetComplexGradMatrix() +template <> +std::unique_ptr SpaceOperator::GetGradMatrix() { return std::make_unique( BuildGrad(GetH1Space(), GetNDSpace(), pa_order_threshold), nullptr, GetH1Space(), @@ -1002,6 +994,25 @@ void SpaceOperator::GetRandomInitialVector(ComplexVector &v) linalg::SetSubVector(v, nd_dbc_tdof_lists.back(), 0.0); } +template std::unique_ptr + SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy); +template std::unique_ptr + SpaceOperator::GetStiffnessMatrix(Operator::DiagonalPolicy); + +template std::unique_ptr + SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy); +template std::unique_ptr + SpaceOperator::GetDampingMatrix(Operator::DiagonalPolicy); + +template std::unique_ptr SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy); +template std::unique_ptr + SpaceOperator::GetMassMatrix(Operator::DiagonalPolicy); + +template std::unique_ptr +SpaceOperator::GetExtraSystemMatrix(double, Operator::DiagonalPolicy); +template std::unique_ptr +SpaceOperator::GetExtraSystemMatrix(double, Operator::DiagonalPolicy); + template std::unique_ptr SpaceOperator::GetSystemMatrix(double, double, double, const Operator *, const Operator *, const Operator *, diff --git a/palace/models/spaceoperator.hpp b/palace/models/spaceoperator.hpp index 5e532d3a5..b68646503 100644 --- a/palace/models/spaceoperator.hpp +++ b/palace/models/spaceoperator.hpp @@ -32,8 +32,9 @@ class SpaceOperator private: const int pa_order_threshold; // Order above which to use partial assembly vs. full const int skip_zeros; // Skip zeros during full assembly of operators - const bool pc_lor; // Use low-order refined (LOR) space for preconditioner - const bool pc_shifted; // Use shifted mass matrix for preconditioner + const bool pc_mat_real; // Use real-valued matrix for preconditioner + const bool pc_mat_shifted; // Use shifted mass matrix for preconditioner + const bool pc_mat_lor; // Use low-order refined (LOR) space for preconditioner // Helper variables for log file printing. bool print_hdr, print_prec_hdr; @@ -129,17 +130,15 @@ class SpaceOperator // A = K + iω C - ω² (Mr + i Mi) + A2(ω) . // For time domain problems, any one of K, C, or M = Mr can be constructed. The argument // ω is required only for the constructing the "extra" matrix A2(ω). - std::unique_ptr GetStiffnessMatrix(Operator::DiagonalPolicy diag_policy); - std::unique_ptr GetDampingMatrix(Operator::DiagonalPolicy diag_policy); - std::unique_ptr GetMassMatrix(Operator::DiagonalPolicy diag_policy); - std::unique_ptr - GetComplexStiffnessMatrix(Operator::DiagonalPolicy diag_policy); - std::unique_ptr - GetComplexDampingMatrix(Operator::DiagonalPolicy diag_policy); - std::unique_ptr - GetComplexMassMatrix(Operator::DiagonalPolicy diag_policy); - std::unique_ptr - GetComplexExtraSystemMatrix(double omega, Operator::DiagonalPolicy diag_policy); + template + std::unique_ptr GetStiffnessMatrix(Operator::DiagonalPolicy diag_policy); + template + std::unique_ptr GetDampingMatrix(Operator::DiagonalPolicy diag_policy); + template + std::unique_ptr GetMassMatrix(Operator::DiagonalPolicy diag_policy); + template + std::unique_ptr GetExtraSystemMatrix(double omega, + Operator::DiagonalPolicy diag_policy); // Construct the complete frequency or time domain system matrix using the provided // stiffness, damping, mass, and extra matrices: @@ -171,10 +170,10 @@ class SpaceOperator // Construct and return the discrete curl or gradient matrices. The complex variants // return a matrix suitable for applying to complex-valued vectors. - std::unique_ptr GetCurlMatrix(); - std::unique_ptr GetComplexCurlMatrix(); - std::unique_ptr GetGradMatrix(); - std::unique_ptr GetComplexGradMatrix(); + template + std::unique_ptr GetCurlMatrix(); + template + std::unique_ptr GetGradMatrix(); // Assemble the right-hand side source term vector for an incident field or current source // applied on specified excited boundaries. The return value indicates whether or not the diff --git a/palace/models/timeoperator.cpp b/palace/models/timeoperator.cpp index 1758c81fb..ac7c77b04 100644 --- a/palace/models/timeoperator.cpp +++ b/palace/models/timeoperator.cpp @@ -52,9 +52,9 @@ class TimeDependentCurlCurlOperator : public mfem::SecondOrderTimeDependentOpera // handled simply by setting diagonal entries of the mass matrix for the corresponding // dofs. Because the Dirichlet BC is always homogenous, no special elimination is // required on the RHS. Diagonal entries are set in M (so M is non-singular). - K = spaceop.GetStiffnessMatrix(Operator::DIAG_ZERO); - C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); - M = spaceop.GetMassMatrix(Operator::DIAG_ONE); + K = spaceop.GetStiffnessMatrix(Operator::DIAG_ZERO); + C = spaceop.GetDampingMatrix(Operator::DIAG_ZERO); + M = spaceop.GetMassMatrix(Operator::DIAG_ONE); // Set up RHS vector for the current source term: -g'(t) J, where g(t) handles the time // dependence. @@ -67,8 +67,7 @@ class TimeDependentCurlCurlOperator : public mfem::SecondOrderTimeDependentOpera pcg->SetInitialGuess(iodata.solver.linear.initial_guess); pcg->SetRelTol(iodata.solver.linear.tol); pcg->SetMaxIter(iodata.solver.linear.max_it); - auto jac = - std::make_unique>(std::make_unique()); + auto jac = std::make_unique>(); kspM = std::make_unique(std::move(pcg), std::move(jac)); kspM->SetOperators(*M, *M); } @@ -144,7 +143,7 @@ TimeOperator::TimeOperator(const IoData &iodata, SpaceOperator &spaceop, std::function &djcoef) { // Construct discrete curl matrix for B-field time integration. - Curl = spaceop.GetCurlMatrix(); + Curl = spaceop.GetCurlMatrix(); // Allocate space for solution vectors. E.SetSize(Curl->Width()); @@ -201,13 +200,11 @@ double TimeOperator::GetMaxTimeStep() const // Solver for M⁻¹. constexpr double lin_tol = 1.0e-9; constexpr int max_lin_it = 500; - mfem::CGSolver pcg(comm); + CgSolver pcg(comm, 0); pcg.SetRelTol(lin_tol); pcg.SetMaxIter(max_lin_it); - pcg.SetPrintLevel(0); pcg.SetOperator(M); - - JacobiSmoother jac; + JacobiSmoother jac; jac.SetOperator(M); pcg.SetPreconditioner(jac); diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 5cb31f1b7..9ca8ef28e 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -1640,14 +1640,18 @@ void LinearSolverData::SetUp(json &solver) mg_max_levels = linear->value("MGMaxLevels", mg_max_levels); mg_coarsen_type = linear->value("MGCoarsenType", mg_coarsen_type); mg_legacy_transfer = linear->value("MGLegacyTransfer", mg_legacy_transfer); - mg_smooth_aux = linear->value("MGAuxiliarySmoother", mg_smooth_aux); mg_cycle_it = linear->value("MGCycleIts", mg_cycle_it); + mg_smooth_aux = linear->value("MGAuxiliarySmoother", mg_smooth_aux); mg_smooth_it = linear->value("MGSmoothIts", mg_smooth_it); mg_smooth_order = linear->value("MGSmoothOrder", mg_smooth_order); + mg_smooth_sf_max = linear->value("MGSmoothEigScaleMax", mg_smooth_sf_max); + mg_smooth_sf_min = linear->value("MGSmoothEigScaleMin", mg_smooth_sf_min); + mg_smooth_cheby_4th = linear->value("MGSmoothChebyshev4th", mg_smooth_cheby_4th); // Preconditioner-specific options - pc_mat_lor = linear->value("PCLowOrderRefined", pc_mat_lor); + pc_mat_real = linear->value("PCMatReal", pc_mat_real); pc_mat_shifted = linear->value("PCMatShifted", pc_mat_shifted); + pc_mat_lor = linear->value("PCLowOrderRefined", pc_mat_lor); pc_side_type = linear->value("PCSide", pc_side_type); sym_fact_type = linear->value("ColumnOrdering", sym_fact_type); strumpack_compression_type = @@ -1673,13 +1677,17 @@ void LinearSolverData::SetUp(json &solver) linear->erase("MGMaxLevels"); linear->erase("MGCoarsenType"); linear->erase("MGLegacyTransfer"); - linear->erase("MGAuxiliarySmoother"); linear->erase("MGCycleIts"); + linear->erase("MGAuxiliarySmoother"); linear->erase("MGSmoothIts"); linear->erase("MGSmoothOrder"); + linear->erase("MGSmoothEigScaleMax"); + linear->erase("MGSmoothEigScaleMin"); + linear->erase("MGSmoothChebyshev4th"); - linear->erase("PCLowOrderRefined"); + linear->erase("PCMatReal"); linear->erase("PCMatShifted"); + linear->erase("PCLowOrderRefined"); linear->erase("PCSide"); linear->erase("ColumnOrdering"); linear->erase("STRUMPACKCompressionType"); @@ -1706,13 +1714,17 @@ void LinearSolverData::SetUp(json &solver) // std::cout << "MGMaxLevels: " << mg_max_levels << '\n'; // std::cout << "MGCoarsenType: " << mg_coarsen_type << '\n'; // std::cout << "MGLegacyTransfer: " << mg_legacy_transfer << '\n'; - // std::cout << "MGAuxiliarySmoother: " << mg_smooth_aux << '\n'; // std::cout << "MGCycleIts: " << mg_cycle_it << '\n'; + // std::cout << "MGAuxiliarySmoother: " << mg_smooth_aux << '\n'; // std::cout << "MGSmoothIts: " << mg_smooth_it << '\n'; // std::cout << "MGSmoothOrder: " << mg_smooth_order << '\n'; + // std::cout << "MGSmoothEigScaleMax: " << mg_smooth_sf_max << '\n'; + // std::cout << "MGSmoothEigScaleMin: " << mg_smooth_sf_min << '\n'; + // std::cout << "MGSmoothChebyshev4th: " << mg_smooth_cheby_4th << '\n'; - // std::cout << "PCLowOrderRefined: " << pc_mat_lor << '\n'; + // std::cout << "PCMatReal: " << pc_mat_real << '\n'; // std::cout << "PCMatShifted: " << pc_mat_shifted << '\n'; + // std::cout << "PCLowOrderRefined: " << pc_mat_lor << '\n'; // std::cout << "PCSide: " << pc_side_type << '\n'; // std::cout << "ColumnOrdering: " << sym_fact_type << '\n'; // std::cout << "STRUMPACKCompressionType: " << strumpack_compression_type << '\n'; diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 0b13d6c1d..0496df948 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -731,13 +731,13 @@ struct LinearSolverData // transfer operators. bool mg_legacy_transfer = false; - // Use auxiliary space smoothers on geometric multigrid levels. - int mg_smooth_aux = -1; - // Number of iterations for preconditioners which support it. For multigrid, this is the // number of V-cycles per Krylov solver iteration. int mg_cycle_it = 1; + // Use auxiliary space smoothers on geometric multigrid levels. + int mg_smooth_aux = -1; + // Number of pre-/post-smoothing iterations at each geometric or algebraic multigrid // level. int mg_smooth_it = 1; @@ -745,14 +745,27 @@ struct LinearSolverData // Order of polynomial smoothing for geometric multigrid. int mg_smooth_order = 4; - // Enable low-order refined (LOR) preconditioner construction. Only available for meshes - // based on tensor-product elements. - bool pc_mat_lor = false; + // Safety factors for eigenvalue estimates associated with Chebyshev smoothing for + // geometric multigrid. + double mg_smooth_sf_max = 1.0; + double mg_smooth_sf_min = 0.0; + + // Smooth based on 4th-kind Chebyshev polynomials for geometric multigrid, otherwise + // use standard 1st-kind polynomials. + bool mg_smooth_cheby_4th = true; + + // For frequency domain applications, precondition linear systems with a real-valued + // approximation to the system matrix. + bool pc_mat_real = false; // For frequency domain applications, precondition linear systems with a shifted matrix // (makes the preconditoner matrix SPD). int pc_mat_shifted = -1; + // Enable low-order refined (LOR) preconditioner construction. Only available for meshes + // based on tensor-product elements. + bool pc_mat_lor = false; + // Choose left or right preconditioning. enum class SideType { diff --git a/scripts/schema/config/solver.json b/scripts/schema/config/solver.json index 498cb3ec4..bae4e1ae1 100644 --- a/scripts/schema/config/solver.json +++ b/scripts/schema/config/solver.json @@ -110,8 +110,12 @@ "MGCycleIts": { "type": "integer", "exclusiveMinimum": 0 }, "MGSmoothIts": { "type": "integer", "exclusiveMinimum": 0 }, "MGSmoothOrder": { "type": "integer", "exclusiveMinimum": 0 }, - "PCLowOrderRefined": { "type": "boolean" }, + "MGSmoothEigScaleMax": { "type": "number", "exclusiveMinimum": 0 }, + "MGSmoothEigScaleMin": { "type": "number", "minimum": 0 }, + "MGSmoothChebyshev4th": { "type": "boolean" }, + "PCMatReal": { "type": "boolean" }, "PCMatShifted": { "type": "boolean" }, + "PCLowOrderRefined": { "type": "boolean" }, "PCSide": { "type": "string" }, "ColumnOrdering": { "type": "string" }, "STRUMPACKCompressionType": { "type": "string" },