feat(poly): preconditioner is now being computed

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;
 
-
+    solverBundle sb = setupNewtonSolver();
-    mfem::NewtonSolver newtonSolver = 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;
+    solver.newton.SetRelTol(newtonRelTol);
-    newtonSolver.SetRelTol(newtonRelTol);
+    solver.newton.SetAbsTol(newtonAbsTol);
-    newtonSolver.SetAbsTol(newtonAbsTol);
+    solver.newton.SetMaxIter(newtonMaxIter);
-    newtonSolver.SetMaxIter(newtonMaxIter);
+    solver.newton.SetPrintLevel(newtonPrintLevel);
-    newtonSolver.SetPrintLevel(newtonPrintLevel);
+    solver.newton.SetOperator(*m_polytropOperator);
-    newtonSolver.SetOperator(*m_polytropOperator);
 
     // --- Created the linear solver which is used to invert the jacobian ---
-    mfem::GMRESSolver gmresSolver;
+    solver.newton.SetSolver(solver.solver);
-    gmresSolver.SetRelTol(gmresRelTol);
-    gmresSolver.SetAbsTol(gmresAbsTol);
-    gmresSolver.SetMaxIter(gmresMaxIter);
-    gmresSolver.SetPrintLevel(gmresPrintLevel);
 
-    // build_preconditioner();
+    return solver;
-
-
-    newtonSolver.SetSolver(gmresSolver);
-
-    return newtonSolver;
 }

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<int>, mfem::Array<int>> 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,

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");

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();