From 58cebc616737a0ffaa36baa1de7649c561ffda2c Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 21 Apr 2025 08:35:29 -0400 Subject: [PATCH] feat(poly): preconditioner is now being computed --- src/poly/solver/private/polySolver.cpp | 41 +++++++++++--------------- src/poly/solver/public/polySolver.h | 8 ++++- src/poly/utils/private/operator.cpp | 4 +++ src/poly/utils/public/operator.h | 2 ++ 4 files changed, 31 insertions(+), 24 deletions(-) diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index adb23ed..1ea037e 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -68,6 +68,7 @@ namespace laneEmden { } } + PolySolver::PolySolver(double n, double order) { // --- Check the polytropic index --- @@ -188,8 +189,7 @@ void PolySolver::solve() const { mfem::Vector zero_rhs(block_offsets.Last()); zero_rhs = 0.0; - - mfem::NewtonSolver newtonSolver = setupNewtonSolver(); + solverBundle sb = setupNewtonSolver(); // EMB 2025: Calling Mult gets the gradient of the operator for the NewtonSolver // This then is set as the operator for the solver for NewtonSolver @@ -199,12 +199,10 @@ void PolySolver::solve() const { // with updating the preconditioner at every newton step as the // changes to the jacobian are automatically propagated through the // solving chain. This is at least true with MFEM 4.8-rc0 - newtonSolver.Mult(zero_rhs, state_vector); + sb.newton.Mult(zero_rhs, state_vector); // --- Save and view the solution --- saveAndViewSolution(state_vector); - - } SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const { @@ -335,32 +333,29 @@ void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); } -mfem::NewtonSolver PolySolver::setupNewtonSolver() const { +solverBundle PolySolver::setupNewtonSolver() const { // --- Load configuration parameters --- double newtonRelTol, newtonAbsTol, gmresRelTol, gmresAbsTol; int newtonMaxIter, newtonPrintLevel, gmresMaxIter, gmresPrintLevel; LoadSolverUserParams(newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel, gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel); + solverBundle solver; + solver.solver.SetRelTol(gmresRelTol); + solver.solver.SetAbsTol(gmresAbsTol); + solver.solver.SetMaxIter(gmresMaxIter); + solver.solver.SetPrintLevel(gmresPrintLevel); + + solver.solver.SetPreconditioner(m_polytropOperator->GetPreconditioner()); // --- Set up the Newton solver --- - mfem::NewtonSolver newtonSolver; - newtonSolver.SetRelTol(newtonRelTol); - newtonSolver.SetAbsTol(newtonAbsTol); - newtonSolver.SetMaxIter(newtonMaxIter); - newtonSolver.SetPrintLevel(newtonPrintLevel); - newtonSolver.SetOperator(*m_polytropOperator); + solver.newton.SetRelTol(newtonRelTol); + solver.newton.SetAbsTol(newtonAbsTol); + solver.newton.SetMaxIter(newtonMaxIter); + solver.newton.SetPrintLevel(newtonPrintLevel); + solver.newton.SetOperator(*m_polytropOperator); // --- Created the linear solver which is used to invert the jacobian --- - mfem::GMRESSolver gmresSolver; - gmresSolver.SetRelTol(gmresRelTol); - gmresSolver.SetAbsTol(gmresAbsTol); - gmresSolver.SetMaxIter(gmresMaxIter); - gmresSolver.SetPrintLevel(gmresPrintLevel); + solver.newton.SetSolver(solver.solver); - // build_preconditioner(); - - - newtonSolver.SetSolver(gmresSolver); - - return newtonSolver; + return solver; } diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 081e087..b0a1a32 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -19,6 +19,12 @@ namespace laneEmden { double thetaSeriesExpansion(double xi, double n, int order); } +// Struct to persist lifetime of the linear and nonlinear solvers +struct solverBundle { + mfem::GMRESSolver solver; + mfem::NewtonSolver newton; +}; + class PolySolver { public: // Public methods PolySolver(double n, double order); @@ -58,7 +64,7 @@ private: // Private methods std::pair, mfem::Array> findCenterElement() const; void setInitialGuess() const; void saveAndViewSolution(const mfem::BlockVector& state_vector) const; - mfem::NewtonSolver setupNewtonSolver() const; + solverBundle setupNewtonSolver() const; void setupOperator() const; void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel, diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp index a21624f..618c820 100644 --- a/src/poly/utils/private/operator.cpp +++ b/src/poly/utils/private/operator.cpp @@ -125,6 +125,10 @@ const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const { return *m_jacobian; } +mfem::BlockDiagonalPreconditioner& PolytropeOperator::GetPreconditioner() const { + return *m_schurPreconditioner; +} + void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const { if (!m_isFinalized) { MFEM_ABORT("PolytropeOperator::Mult called before finalize"); diff --git a/src/poly/utils/public/operator.h b/src/poly/utils/public/operator.h index 3e2c16b..42f5797 100644 --- a/src/poly/utils/public/operator.h +++ b/src/poly/utils/public/operator.h @@ -35,6 +35,8 @@ public: const mfem::BlockOperator &GetJacobianOperator() const; + mfem::BlockDiagonalPreconditioner &GetPreconditioner() const; + private: Probe::LogManager& m_logManager = Probe::LogManager::getInstance();