diff --git a/src/poly/solver/private/polySolver.cpp b/src/poly/solver/private/polySolver.cpp index 5ae30f1..29c221c 100644 --- a/src/poly/solver/private/polySolver.cpp +++ b/src/poly/solver/private/polySolver.cpp @@ -194,6 +194,9 @@ void PolySolver::solve() const { const solverBundle sb = setupNewtonSolver(); sb.newton.Mult(zero_rhs, state_vector); + // TODO: Need some system here to reconstruct the full system vector from what MFEM's Newton Solver + // will return, since that will end up being the vector for the reduced system. + // --- Save and view an approximate 1D solution --- saveAndViewSolution(state_vector); } diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 11c4071..8b14be6 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -190,4 +190,4 @@ private: // Private methods static void GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename); -}; \ No newline at end of file +}; diff --git a/src/poly/utils/private/polytropeOperator.cpp b/src/poly/utils/private/polytropeOperator.cpp index 6eab1e3..1be8402 100644 --- a/src/poly/utils/private/polytropeOperator.cpp +++ b/src/poly/utils/private/polytropeOperator.cpp @@ -168,6 +168,11 @@ int PolytropeOperator::get_reduced_system_size() const { return m_reducedBlockOffsets.Last(); } +const mfem::Vector &PolytropeOperator::reconstruct_full_state_vector(const mfem::Vector &reducedState) const { + m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector + return m_state; +} + void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const { if (!m_isFinalized) { MFEM_ABORT("PolytropeOperator::Mult called before finalize"); diff --git a/src/poly/utils/public/polytropeOperator.h b/src/poly/utils/public/polytropeOperator.h index 6b905ee..0714c64 100644 --- a/src/poly/utils/public/polytropeOperator.h +++ b/src/poly/utils/public/polytropeOperator.h @@ -26,45 +26,128 @@ #include "probe.h" +/** + * @brief Represents the Schur complement operator used in the solution process. + * + * This class computes S = D - Q * GradInv * M. + */ class SchurCompliment final : public mfem::Operator { public: + /** + * @brief Constructs a SchurCompliment operator. + * @param QOp The Q matrix operator. + * @param DOp The D matrix operator. + * @param MOp The M matrix operator. + * @param GradInvOp The inverse of the gradient operator. + */ SchurCompliment(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp); + + /** + * @brief Constructs a SchurCompliment operator without the inverse gradient initially. + * The inverse gradient must be set later using updateInverseNonlinearJacobian. + * @param QOp The Q matrix operator. + * @param DOp The D matrix operator. + * @param MOp The M matrix operator. + */ SchurCompliment(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp); + + /** + * @brief Destructor. + */ ~SchurCompliment() override = default; + + /** + * @brief Applies the Schur complement operator: y = (D - Q * GradInv * M) * x. + * @param x The input vector. + * @param y The output vector. + */ void Mult(const mfem::Vector &x, mfem::Vector &y) const override; - void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver & - GradInvOp); + + /** + * @brief Sets all operators for the Schur complement. + * @param QOp The Q matrix operator. + * @param DOp The D matrix operator. + * @param MOp The M matrix operator. + * @param GradInvOp The inverse of the gradient operator. + */ + void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp); + + /** + * @brief Updates the inverse of the nonlinear Jacobian (GradInv). + * @param gradInv The new inverse gradient solver. + */ void updateInverseNonlinearJacobian(const mfem::Solver &gradInv); private: + /** + * @brief Updates the constant matrix terms (Q, D, M). + * @param QOp The Q matrix operator. + * @param DOp The D matrix operator. + * @param MOp The M matrix operator. + */ void updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp); private: + // Pointers to external operators (not owned by this class) + const mfem::SparseMatrix* m_QOp = nullptr; ///< Pointer to the Q matrix operator. + const mfem::SparseMatrix* m_DOp = nullptr; ///< Pointer to the D matrix operator. + const mfem::SparseMatrix* m_MOp = nullptr; ///< Pointer to the M matrix operator. + const mfem::Solver* m_GradInvOp = nullptr; ///< Pointer to the inverse of the gradient operator. - // Note that these are not owned by this class - const mfem::SparseMatrix* m_QOp = nullptr; - const mfem::SparseMatrix* m_DOp = nullptr; - const mfem::SparseMatrix* m_MOp = nullptr; - const mfem::Solver* m_GradInvOp = nullptr; - int m_nPhi = 0; - int m_nTheta = 0; + // Dimensions + int m_nPhi = 0; ///< Dimension related to the phi variable (typically number of rows/cols of D). + int m_nTheta = 0; ///< Dimension related to the theta variable (typically number of rows of M). - mutable std::unique_ptr m_matrixForm; + // Owned resources + mutable std::unique_ptr m_matrixForm; ///< Optional: if a matrix representation is ever explicitly formed. }; +/** + * @brief Provides an approximate inverse of the SchurCompliment operator using GMRES. + */ class GMRESInverter final : public mfem::Operator { public: + /** + * @brief Constructs a GMRESInverter. + * @param op The SchurCompliment operator to invert. + */ explicit GMRESInverter(const SchurCompliment& op); + + /** + * @brief Destructor. + */ ~GMRESInverter() override = default; + + /** + * @brief Applies the GMRES solver to approximate op^-1 * x. + * @param x The input vector. + * @param y The output vector (approximation of op^-1 * x). + */ void Mult(const mfem::Vector &x, mfem::Vector &y) const override; private: - const SchurCompliment& m_op; - mfem::GMRESSolver m_solver; + const SchurCompliment& m_op; ///< Reference to the SchurCompliment operator to be inverted. + mfem::GMRESSolver m_solver; ///< GMRES solver instance. }; +/** + * @brief Represents the coupled nonlinear operator for the polytropic system. + * + * This operator defines the system F(X) = 0, where X = [θ, φ]^T, + * and F(X) = [ f(θ) + Mφ ] + * [ Dφ - Qθ ]. + * It handles essential boundary conditions and the construction of the Jacobian. + */ class PolytropeOperator final : public mfem::Operator { public: + /** + * @brief Constructs the PolytropeOperator. + * @param M The M bilinear form (coupling θ and φ). + * @param Q The Q bilinear form (coupling φ and θ). + * @param D The D bilinear form (acting on φ). + * @param f The nonlinear form f(θ). + * @param blockOffsets Offsets defining the blocks for θ and φ variables. + */ PolytropeOperator( std::unique_ptr M, std::unique_ptr Q, @@ -72,78 +155,197 @@ public: std::unique_ptr f, const mfem::Array &blockOffsets); - void populate_free_dof_array(); - - + /** + * @brief Destructor. + */ ~PolytropeOperator() override = default; + /** + * @brief Applies the PolytropeOperator: y = F(x). + * This computes the residual of the nonlinear system. + * @param xFree The vector of free (non-essential) degrees of freedom. + * @param yFree The resulting residual vector corresponding to free DOFs. + */ void Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const override; + /** + * @brief Computes the Jacobian of the PolytropeOperator at a given state xFree. + * The Jacobian is of the form: + * J = [ G M ] + * [ -Q D ] + * where G is the gradient of f(θ). + * @param xFree The vector of free DOFs at which to evaluate the gradient. + * @return A reference to the Jacobian operator. + */ mfem::Operator& GetGradient(const mfem::Vector &xFree) const override; - void set_essential_true_dofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs); - void set_essential_true_dofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set); - - SSE::MFEMArrayPairSet get_essential_true_dofs() const; - - bool isFinalized() const { return m_isFinalized; } - + /** + * @brief Finalizes the operator setup. + * This must be called after setting essential true DOFs and before using Mult or GetGradient. + * It constructs reduced matrices, block offsets, and populates free DOFs. + * @param initTheta Initial guess for θ, used to evaluate the nonlinear form if necessary during setup. + */ void finalize(const mfem::Vector &initTheta); + /** + * @brief Checks if the operator has been finalized. + * @return True if finalize() has been successfully called, false otherwise. + */ + bool isFinalized() const { return m_isFinalized; } + + /** + * @brief Sets the essential true degrees of freedom for both θ and φ variables. + * Marks the operator as not finalized. + * @param theta_ess_tdofs Pair of arrays: (indices of essential DOFs for θ, values at these DOFs). + * @param phi_ess_tdofs Pair of arrays: (indices of essential DOFs for φ, values at these DOFs). + */ + void set_essential_true_dofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs); + + /** + * @brief Sets the essential true degrees of freedom using a pair of pairs. + * Marks the operator as not finalized. + * @param ess_tdof_pair_set A pair containing the essential TDOF pairs for theta and phi. + */ + void set_essential_true_dofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set); + + /** + * @brief Gets the currently set essential true degrees of freedom. + * @return A pair containing the essential TDOF pairs for theta and phi. + */ + SSE::MFEMArrayPairSet get_essential_true_dofs() const; + + /** + * @brief Gets the block offsets for the full (unreduced) system. + * @return Constant reference to the array of block offsets. + */ const mfem::Array& get_block_offsets() const { return m_blockOffsets; } + + /** + * @brief Gets the block offsets for the reduced system (after eliminating essential DOFs). + * @return Constant reference to the array of reduced block offsets. + */ const mfem::Array& get_reduced_block_offsets() const {return m_reducedBlockOffsets; } + /** + * @brief Gets the Jacobian operator. + * Asserts that the operator is finalized and the Jacobian has been computed. + * @return Constant reference to the block Jacobian operator. + */ const mfem::BlockOperator &get_jacobian_operator() const; + /** + * @brief Gets the block diagonal preconditioner for the Schur complement system. + * Asserts that the operator is finalized and the preconditioner has been computed. + * @return Reference to the block diagonal preconditioner. + */ mfem::BlockDiagonalPreconditioner &get_preconditioner() const; + /** + * @brief Gets the size of the reduced system (number of free DOFs). + * Asserts that the operator is finalized. + * @return The total number of free degrees of freedom. + */ int get_reduced_system_size() const; + /** + * @brief Reconstructs the full state vector (including essential DOFs) from a reduced state vector (free DOFs). + * @param reducedState The vector containing only the free degrees of freedom. + * @return Constant reference to the internal full state vector, updated with the reducedState. + */ + [[nodiscard]] const mfem::Vector &reconstruct_full_state_vector(const mfem::Vector &reducedState) const; + + /** + * @brief Populates the internal array of free (non-essential) degree of freedom indices. + * This is called during finalize(). + */ + void populate_free_dof_array(); private: - Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); - quill::Logger* m_logger = m_logManager.getLogger("log"); - std::unique_ptr m_M; - std::unique_ptr m_Q; - std::unique_ptr m_D; - std::unique_ptr m_f; + // --- Logging --- + Probe::LogManager& m_logManager = Probe::LogManager::getInstance(); ///< Reference to the global log manager. + quill::Logger* m_logger = m_logManager.getLogger("log"); ///< Pointer to the specific logger instance. - std::unique_ptr m_Mmat; - std::unique_ptr m_Qmat; - std::unique_ptr m_Dmat; + // --- Input Bilinear/Nonlinear Forms (owned) --- + std::unique_ptr m_M; ///< Bilinear form M, coupling θ and φ. + std::unique_ptr m_Q; ///< Bilinear form Q, coupling φ and θ. + std::unique_ptr m_D; ///< Bilinear form D, acting on φ. + std::unique_ptr m_f; ///< Nonlinear form f, acting on θ. - mutable mfem::Vector m_state; - mfem::Array m_freeDofs; + // --- Full Matrix Representations (owned, derived from forms) --- + std::unique_ptr m_Mmat; ///< Sparse matrix representation of M. + std::unique_ptr m_Qmat; ///< Sparse matrix representation of Q. + std::unique_ptr m_Dmat; ///< Sparse matrix representation of D. - // These are used to store the reduced system matrices after essential dofs have been eliminated - std::unique_ptr m_MReduced; - std::unique_ptr m_QReduced; - std::unique_ptr m_DReduced; + // --- Reduced Matrix Representations (owned, after eliminating essential DOFs) --- + std::unique_ptr m_MReduced; ///< Reduced M matrix (free DOFs only). + std::unique_ptr m_QReduced; ///< Reduced Q matrix (free DOFs only). + std::unique_ptr m_DReduced; ///< Reduced D matrix (free DOFs only). - const mfem::Array m_blockOffsets; - mfem::Array m_reducedBlockOffsets; + // --- State Vectors and DOF Management --- + mutable mfem::Vector m_state; ///< Full state vector [θ, φ]^T, including essential DOFs. + mfem::Array m_freeDofs; ///< Array of indices for free (non-essential) DOFs. - SSE::MFEMArrayPair m_theta_ess_tdofs; - SSE::MFEMArrayPair m_phi_ess_tdofs; + // --- Block Offsets --- + const mfem::Array m_blockOffsets; ///< Block offsets for the full system [0, size(θ), size(θ)+size(φ)]. + mfem::Array m_reducedBlockOffsets; ///< Block offsets for the reduced system (free DOFs). - std::unique_ptr m_negQ_mat; - mutable std::unique_ptr m_jacobian; + // --- Essential Boundary Conditions --- + SSE::MFEMArrayPair m_theta_ess_tdofs; ///< Essential true DOFs for θ (indices and values). + SSE::MFEMArrayPair m_phi_ess_tdofs; ///< Essential true DOFs for φ (indices and values). - mutable std::unique_ptr m_schurCompliment; - mutable std::unique_ptr m_invSchurCompliment; - mutable std::unique_ptr m_invNonlinearJacobian; + // --- Jacobian and Preconditioner Components (owned, mutable) --- + std::unique_ptr m_negQ_mat; ///< Scaled operator for -Q_reduced. + mutable std::unique_ptr m_jacobian; ///< Jacobian operator J = [G M; -Q D]_reduced. + mutable std::unique_ptr m_schurCompliment; ///< Schur complement S = D_reduced - Q_reduced * G_inv_reduced * M_reduced. + mutable std::unique_ptr m_invSchurCompliment; ///< Approximate inverse of the Schur complement. + mutable std::unique_ptr m_invNonlinearJacobian; ///< Solver for the inverse of the G block (gradient of f(θ)_reduced). + mutable std::unique_ptr m_schurPreconditioner; ///< Block diagonal preconditioner for the system. - mutable std::unique_ptr m_schurPreconditioner; - - bool m_isFinalized = false; + // --- State Flags --- + bool m_isFinalized = false; ///< Flag indicating if finalize() has been called. private: - void update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const; - void update_inverse_schur_compliment() const; - void update_preconditioner(const mfem::Operator &grad) const; - void scatter_boundary_conditions(); + /** + * @brief Constructs the sparse matrix representations of M, Q, and D, and their reduced forms. + * Called during finalize(). + */ void construct_matrix_representations(); - void construct_reduced_block_offsets(); - void construct_jacobian_constant_terms(); -}; + /** + * @brief Constructs the block offsets for the reduced system. + * Called during finalize(). + */ + void construct_reduced_block_offsets(); + + /** + * @brief Constructs the constant terms of the Jacobian operator (M, -Q, D). + * The (0,0) block (gradient of f) is set in GetGradient. + * Called during finalize(). + */ + void construct_jacobian_constant_terms(); + + /** + * @brief Scatters the values of essential boundary conditions into the full state vector. + * Called during finalize(). + */ + void scatter_boundary_conditions(); + + /** + * @brief Updates the solver for the inverse of the nonlinear Jacobian block (G_00). + * @param grad The gradient operator (G_00) of the nonlinear part f(θ). + */ + void update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const; + + /** + * @brief Updates the inverse Schur complement operator and its components. + * This is typically called after the nonlinear Jacobian part has been updated. + */ + void update_inverse_schur_compliment() const; + + /** + * @brief Updates the preconditioner components. + * This involves updating the inverse nonlinear Jacobian and then the inverse Schur complement. + * @param grad The gradient operator (G_00) of the nonlinear part f(θ). + */ + void update_preconditioner(const mfem::Operator &grad) const; +};