diff --git a/src/poly/solver/public/polySolver.h b/src/poly/solver/public/polySolver.h index 3e140fd..8d0007e 100644 --- a/src/poly/solver/public/polySolver.h +++ b/src/poly/solver/public/polySolver.h @@ -32,161 +32,579 @@ #include "probe.h" #include "quill/Logger.h" - +/** + * @brief Namespace for Lane-Emden equation related utility functions. + * + * Provides functions to compute coefficients and evaluate the series expansion + * solution to the Lane-Emden equation, which describes the structure of a + * spherically symmetric polytropic star. + * The Lane-Emden equation is given by: + * \f[ + * \frac{1}{\xi^2} \frac{d}{d\xi} \left( \xi^2 \frac{d\theta}{d\xi} \right) = -\theta^n + * \f] + * where \f$\xi\f$ is a dimensionless radius and \f$\theta\f$ is related to the density, + * and \f$n\f$ is the polytropic index. + */ namespace laneEmden { + /** + * @brief Computes the coefficient \f$a_k\f$ for the Lane-Emden series expansion. + * + * The series solution for \f$\theta(\xi)\f$ is given by \f$\theta(\xi) = \sum_{k=0}^{\infty} a_k \xi^k\f$. + * The coefficients \f$a_k\f$ are determined by substituting the series into the Lane-Emden equation. + * Specifically, \f$a_0 = 1\f$, \f$a_1 = 0\f$, and for \f$k \ge 2\f$, + * \f$a_k = -\frac{c_{k-2,n}}{k(k+1)}\f$. + * + * @param k The index of the coefficient. + * @param n The polytropic index. + * @return The value of the coefficient \f$a_k\f$. + * @see c(const int m, const double n) + */ double a (const int k, const double n); + + /** + * @brief Computes the auxiliary coefficient \f$c_{m,n}\f$ used in determining \f$a_k\f$. + * + * The term \f$\theta^n\f$ in the Lane-Emden equation can also be expanded as a series + * \f$\theta^n(\xi) = \sum_{m=0}^{\infty} c_{m,n} \xi^m\f$. + * The coefficients \f$c_{m,n}\f$ are related to \f$a_k\f$ by: + * \f$c_{0,n} = a_0^n\f$ + * \f$c_{m,n} = \frac{1}{m a_0} \sum_{j=1}^{m} (j n - m + j) a_j c_{m-j,n}\f$ for \f$m > 0\f$. + * + * @param m The index of the coefficient. + * @param n The polytropic index. + * @return The value of the coefficient \f$c_{m,n}\f$. + * @see a(const int k, const double n) + */ double c(const int m, const double n); + + /** + * @brief Computes the Lane-Emden function \f$\theta(\xi)\f$ using a series expansion. + * + * Evaluates the series \f$\theta(\xi) = \sum_{k=0}^{\text{order}-1} a_k \xi^k\f$ up to a specified order. + * This provides an approximate solution to the Lane-Emden equation, particularly accurate for small \f$\xi\f$. + * + * @param xi The dimensionless radius \f$\xi\f$. + * @param n The polytropic index. + * @param order The number of terms to include in the series expansion. + * @return The approximate value of \f$\theta(\xi)\f$. + * + * @example + * @code + * double xi = 0.5; + * double n = 1.5; + * int series_order = 10; + * double theta_val = laneEmden::thetaSeriesExpansion(xi, n, series_order); + * // theta_val will be an approximation of the Lane-Emden function at xi=0.5 for n=1.5 + * @endcode + */ double thetaSeriesExpansion(const double xi, const double n, const int order); } -// Struct to persist lifetime of the linear and nonlinear solvers +/** + * @brief Structure to manage the lifetime of MFEM solver objects. + * + * This structure ensures that the `mfem::GMRESSolver` outlives the `mfem::NewtonSolver` + * that might use it. The `mfem::NewtonSolver` often takes a pointer or reference to a + * linear solver, and if the linear solver is destroyed first, it can lead to dangling + * pointers and crashes. + * + * @note The order of declaration of members is crucial: `solver` must be declared + * before `newton` to ensure it is constructed first and destroyed last. + */ struct solverBundle { - mfem::GMRESSolver solver; // Must be first so it lives longer than the newton solver - mfem::NewtonSolver newton; // Must be second so that when it is destroyed the solver is still alive preventing a double delete + mfem::GMRESSolver solver; ///< The linear solver (e.g., GMRES). Must be declared first. + mfem::NewtonSolver newton; ///< The nonlinear solver (e.g., Newton). Must be declared second. }; +/** + * @brief Structure to hold the various bilinear and nonlinear forms for the polytrope problem. + * + * This structure bundles the MFEM forms that represent the discretized weak formulation + * of the mixed variable polytropic equations. + * The system being solved is typically: + * \f[ + * \begin{aligned} + * \boldsymbol{\phi} + \nabla \theta &= \mathbf{0} \\ + * \nabla \cdot \boldsymbol{\phi} - \theta^n &= 0 + * \end{aligned} + * \f] + * After integration by parts and discretization, these lead to matrix equations involving + * the forms stored here. + */ struct formBundle { - std::unique_ptr M; //<-- M ∫∇ψ^θ·N^φ dV - std::unique_ptr Q; //<-- Q ∫ψ^φ·∇N^θ dV - std::unique_ptr D; //<-- D ∫ψ^φ·N^φ dV - std::unique_ptr S; //<-- S ∫∇ψ^θ·∇N^θ dV - std::unique_ptr f; //<-- f(θ) ∫ψ^θ·θ^n dV + /** + * @brief Mixed bilinear form \f$ M(\psi^\theta, N^\phi) = \int_\Omega \nabla\psi^\theta \cdot N^\phi \,dV \f$. + * This form arises from the term \f$\nabla \theta\f$ in the first equation, tested with a vector test function \f$N^\phi\f$. + * It couples the scalar field \f$\theta\f$ (potential) with the vector field \f$\boldsymbol{\phi}\f$ (flux). + */ + std::unique_ptr M; + + /** + * @brief Mixed bilinear form \f$ Q(\psi^\phi, N^\theta) = \int_\Omega \psi^\phi \cdot \nabla N^\theta \,dV \f$. + * This form arises from the term \f$\boldsymbol{\phi}\f$ in the first equation, tested with a scalar test function \f$N^\theta\f$, + * after integration by parts of the \f$\nabla \theta\f$ term. + * It can also be seen as related to the transpose of a discrete gradient operator. + */ + std::unique_ptr Q; + + /** + * @brief Bilinear form \f$ D(\psi^\phi, N^\phi) = \int_\Omega \psi^\phi \cdot N^\phi \,dV \f$. + * This is a mass matrix for the vector field \f$\boldsymbol{\phi}\f$. It arises from the \f$\boldsymbol{\phi}\f$ term + * in the first equation when tested with a vector test function \f$N^\phi\f$. + */ + std::unique_ptr D; + + /** + * @brief Bilinear form \f$ S(\psi^\theta, N^\theta) = \int_\Omega \nabla\psi^\theta \cdot \nabla N^\theta \,dV \f$. + * This is a stiffness matrix (Laplacian) for the scalar field \f$\theta\f$. It is used for stabilization terms + * or alternative formulations. + */ + std::unique_ptr S; + + /** + * @brief Nonlinear form \f$ f(\theta)(\psi^\theta) = \int_\Omega \psi^\theta \cdot \theta^n \,dV \f$. + * This form arises from the nonlinear source term \f$\theta^n\f$ in the second equation, tested with a scalar + * test function \f$\psi^\theta\f$. + */ + std::unique_ptr f; }; +/** + * @brief Solves the Lane-Emden equation for a polytropic star using a mixed finite element method. + * + * This class sets up and solves the system of equations describing a polytropic + * stellar model. The Lane-Emden equation, in its second-order form, is: + * \f[ + * \frac{1}{\xi^2} \frac{d}{d\xi} \left( \xi^2 \frac{d\theta}{d\xi} \right) = -\theta^n + * \f] + * where \f$\xi\f$ is a dimensionless radius, \f$\theta\f$ is a dimensionless temperature/potential + * related to density (\f$\rho = \rho_c \theta^n\f$), and \f$n\f$ is the polytropic index. + * + * This solver uses a mixed formulation by introducing \f$\boldsymbol{\phi} = -\nabla \theta\f$: + * \f[ + * \begin{aligned} + * \boldsymbol{\phi} + \nabla \theta &= \mathbf{0} \\ + * \nabla \cdot \boldsymbol{\phi} - \theta^n &= 0 + * \end{aligned} + * \f] + * These equations are discretized using H1 finite elements for \f$\theta\f$ and + * RT (Raviart-Thomas) finite elements for \f$\boldsymbol{\phi}\f$. The resulting nonlinear + * system is solved using a Newton method, with the `PolytropeOperator` class + * defining the system residual and Jacobian. + * + * Boundary conditions typically include \f$\theta(\xi_c) = 1\f$ at the center (\f$\xi_c \approx 0\f$), + * \f$\theta(\xi_1) = 0\f$ at the surface (\f$\xi_1\f$ is the first root of \f$\theta\f$), and + * \f$\boldsymbol{\phi} \cdot \mathbf{n} = -\frac{d\theta}{d\xi}\f$ at the surface, which is related to the + * surface gravity. At the center, \f$\boldsymbol{\phi}(\xi_c) = \mathbf{0}\f$ due to symmetry. + */ class PolySolver final{ public: // Public methods + /** + * @brief Constructs a PolySolver instance. + * + * Initializes the solver for a given polytropic index \f$n\f$ and finite element order. + * This constructor internally calls `prepareMesh` to load and scale the mesh. + * + * @param n The polytropic index. Must be \f$0 \le n < 5\f$. + * @param order The polynomial order of the H1 finite elements for \f$\theta\f$. + * The RT elements for \f$\boldsymbol{\phi}\f$ will be of order `order-1`. + * @throw std::runtime_error If \f$n \ge 5\f$ or \f$n < 0\f$ (via `prepareMesh`). + * + * @example + * @code + * try { + * PolySolver solver(1.5, 2); // Polytropic index n=1.5, FE order 2 for theta + * solver.solve(); + * } catch (const std::exception& e) { + * std::cerr << "Error: " << e.what() << std::endl; + * } + * @endcode + */ PolySolver(const double n, const double order); + + /** + * @brief Destructor for PolySolver. + * + * Defaulted, handles cleanup of owned resources like `std::unique_ptr` members. + */ ~PolySolver(); + /** + * @brief Solves the polytropic system. + * + * This method orchestrates the solution process: + * 1. Sets the initial guess for \f$\theta\f$ and \f$\boldsymbol{\phi}\f$. + * 2. Applies essential boundary conditions to the `PolytropeOperator`. + * 3. Finalizes the `PolytropeOperator` with the initial state. + * 4. Sets up the Newton solver (including the GMRES linear solver for Jacobian systems). + * 5. Runs the Newton iteration to find the solution for the free DOFs. + * 6. Reconstructs the full solution vector from the free DOFs. + * 7. Saves and/or views the solution based on configuration. + * + * @throw std::runtime_error or `mfem::ErrorException` if the solver fails to converge or encounters numerical issues. + */ void solve() const; + /** + * @brief Gets the polytropic index \f$n\f$. + * @return The polytropic index. + */ double getN() const { return m_polytropicIndex; } + + /** + * @brief Gets the finite element order used for the \f$\theta\f$ variable. + * @return The polynomial order of the H1 finite elements. + */ double getOrder() const { return m_feOrder; } + + /** + * @brief Gets a reference to the computational mesh. + * @return A reference to the `mfem::Mesh` object. + */ mfem::Mesh& getMesh() const { return m_mesh; } + + /** + * @brief Gets a reference to the solution grid function for \f$\theta\f$. + * @return A reference to the `mfem::GridFunction` storing the \f$\theta\f$ solution. + * @note The solution is populated after `solve()` has been successfully called. + */ mfem::GridFunction& getSolution() const { return *m_theta; } private: // Private Attributes - Config& m_config; - Probe::LogManager& m_logManager; - quill::Logger* m_logger; - double m_polytropicIndex, m_feOrder; + // --- Configuration and Logging --- + Config& m_config; ///< Reference to the global configuration manager instance. + Probe::LogManager& m_logManager; ///< Reference to the global log manager instance. + quill::Logger* m_logger; ///< Pointer to the specific logger instance for this class. - mfem::Mesh& m_mesh; - std::unique_ptr m_fecH1; - std::unique_ptr m_fecRT; + // --- Physical and Discretization Parameters --- + double m_polytropicIndex; ///< The polytropic index \f$n\f$. + double m_feOrder; ///< The polynomial order for H1 elements (\f$\theta\f$). RT elements for \f$\boldsymbol{\phi}\f$ are of order `m_feOrder - 1`. - std::unique_ptr m_feTheta; - std::unique_ptr m_fePhi; + // --- MFEM Core Objects --- + mfem::Mesh& m_mesh; ///< Reference to the computational mesh (owned by ResourceManager). + std::unique_ptr m_fecH1; ///< Finite Element Collection for \f$\theta\f$ (H1 elements). + std::unique_ptr m_fecRT; ///< Finite Element Collection for \f$\boldsymbol{\phi}\f$ (Raviart-Thomas elements). - std::unique_ptr m_theta; - std::unique_ptr m_phi; + std::unique_ptr m_feTheta; ///< Finite Element Space for \f$\theta\f$. + std::unique_ptr m_fePhi; ///< Finite Element Space for \f$\boldsymbol{\phi}\f$. - std::unique_ptr m_polytropOperator; - - std::unique_ptr m_prec; + // --- Solution Vectors --- + std::unique_ptr m_theta; ///< Grid function for the scalar potential \f$\theta\f$. + std::unique_ptr m_phi; ///< Grid function for the vector flux \f$\boldsymbol{\phi}\f$. + // --- Operator and Solver Components --- + std::unique_ptr m_polytropOperator; ///< The main nonlinear operator for the mixed system. + std::unique_ptr m_prec; ///< Preconditioner (currently seems unused in `polySolver.cpp`). private: // Private methods + /** + * @brief Private constructor that takes an existing mesh. + * + * Initializes FE collections, spaces, grid functions, and assembles the block system. + * This is called by the public constructor after `prepareMesh`. + * + * @param mesh A reference to an initialized `mfem::Mesh`. + * @param n The polytropic index. + * @param order The polynomial order for H1 finite elements. + */ PolySolver(mfem::Mesh& mesh, double n, double order); + /** + * @brief Prepares the mesh for the simulation. + * + * Loads a generic sphere mesh from the `ResourceManager` and scales it + * to the dimensionless radius \f$\xi_1(n)\f$ (the first zero of the Lane-Emden function + * for the given polytropic index \f$n\f$). + * + * @param n The polytropic index. + * @return A reference to the prepared `mfem::Mesh` (owned by `ResourceManager`). + * @throw std::runtime_error If \f$n \ge 5.0\f$ or \f$n < 0.0\f$, as \f$\xi_1(n)\f$ is typically + * undefined or infinite outside this range for physical polytropes. + * @throw std::runtime_error If the mesh resource "mesh:polySphere" cannot be found or loaded. + * + * @details The scaling factor \f$\xi_1(n)\f$ is obtained from `polycoeff::x1(n)`. + * The mesh is expected to be a unit sphere initially. + */ static mfem::Mesh& prepareMesh(double n); - + /** + * @brief Assembles the block system operator `m_polytropOperator`. + * + * This involves: + * 1. Computing block offsets for \f$\theta\f$ and \f$\boldsymbol{\phi}\f$ variables. + * 2. Building individual bilinear and nonlinear forms (M, Q, D, S, f) using `buildIndividualForms`. + * 3. Constructing the `PolytropeOperator` with these forms and offsets. + */ void assembleBlockSystem(); /** - * @breif Compute the block offsets for the operator. These are the offsets that define which dofs belong to which variable. + * @brief Compute the block offsets for the operator. These are the offsets that define which dofs belong to which variable. * * @details - * * Create the block offsets. These define the start of each block in the combined vector. - * Block offsets will be [0, thetaDofs, thetaDofs + phiDofs]. + * Block offsets will be \f$[0, N_\theta, N_\theta + N_\phi]\f$, where \f$N_\theta\f$ is the number of + * degrees of freedom for \f$\theta\f$ and \f$N_\phi\f$ is for \f$\boldsymbol{\phi}\f$. * The interpretation of this is that each block tells the operator where in the flattened (1D) vector * the degrees of freedom or coefficients for that free parameter start and end. I.e. - * we know that in any flattened vector will have a size thetaDofs + phiDofs. The theta dofs will span - * from blockOffsets[0] -> blockOffsets[1] and the phiDofs will span from blockOffsets[1] -> blockOffsets[2]. + * we know that in any flattened vector will have a size \f$N_\theta + N_\phi\f$. The \f$\theta\f$ dofs will span + * from `blockOffsets[0]` to `blockOffsets[1]` and the \f$\boldsymbol{\phi}\f$ dofs will span from `blockOffsets[1]` to `blockOffsets[2]`. * - * This is the same for matrices only in 2D (rows and columns) + * This is the same for matrices only in 2D (rows and columns). * * The key point here is that this is fundamentally an accounting structure, it is here to keep track of what * parts of vectors and matrices belong to which variable. * - * Also note that we use VSize rather than Size. Size referees to the number of true dofs. That is the dofs which - * still are present in the system after eliminating boundary conditions. This is the wrong size to use if we are - * trying to account for the true size of the system. VSize on the other hand refers to the total number of dofs. + * Also note that we use `VSize` rather than `Size`. `Size` refers to the number of true dofs (after eliminating + * boundary conditions). `VSize` refers to the total number of dofs before BC elimination, which is needed here. * - * @return blockOffsets The offsets for the blocks in the operator + * @return `mfem::Array` The offsets for the blocks in the operator. * - * @pre m_feTheta and m_fePhi must be valid, populated FiniteElementSpaces. + * @pre `m_feTheta` and `m_fePhi` must be valid, populated `mfem::FiniteElementSpace` objects. */ mfem::Array computeBlockOffsets() const; /** - * @breif Build the individual forms for the block operator (M, Q, D, and f) + * @brief Build the individual forms for the block operator (M, Q, D, S, and f). * - * @param blockOffsets The offsets for the blocks in the operator - * @return forms The forms for the block operator + * @param blockOffsets The offsets for the blocks in the operator, computed by `computeBlockOffsets`. + * @return A `std::unique_ptr` containing the assembled forms. * - * @note These forms are build exactly how they are defined in the derivation. This means that Mform -> M not -M and Qform -> Q not -Q. + * @note These forms are built exactly how they are defined in the derivation. + * For example, `Mform` corresponds to \f$+M\f$, not \f$-M\f$. + * The `PolytropeOperator` handles any necessary sign changes for the final system assembly. * * @details - * Computes the block offsets - * \f$\{0,\;|\theta|,\;|\theta|+|\phi|\}\f$, then builds and finalizes - * the MixedBilinearForms \c Mform, \c Qform and the BilinearForm \c Dform, - * plus the NonlinearForm \c fform. Finally, these are handed off to - * \c PolytropeOperator along with the offsets. - * - * The discretized weak form is + * This method initializes and assembles the discrete forms corresponding to the weak formulation + * of the mixed polytropic system: * \f[ - * R(X) - * = \begin{pmatrix} - * f(\theta) - M\,\theta \\[6pt] - * D\,\theta - Q\,\phi + * M = \int_\Omega \nabla\psi^\theta \cdot N^\phi \,dV \quad (\text{from } \nabla\theta \text{ term}) + * \f] + * \f[ + * Q = \int_\Omega \psi^\phi \cdot \nabla N^\theta \,dV \quad (\text{from } \boldsymbol{\phi} \text{ term, related to } -M^T) + * \f] + * \f[ + * D = \int_\Omega \psi^\phi \cdot N^\phi \,dV \quad (\text{mass matrix for } \boldsymbol{\phi}) + * \f] + * \f[ + * S = \int_\Omega \nabla\psi^\theta \cdot \nabla N^\theta \,dV \quad (\text{stiffness matrix for } \theta \text{, for stabilization}) + * \f] + * \f[ + * f(\theta)(\psi^\theta) = \int_\Omega \psi^\theta \cdot \theta^n \,dV \quad (\text{from nonlinear term } \theta^n) + * \f] + * The `PolytropeOperator` will then use these to form a system like: + * \f[ + * R(X) = \begin{pmatrix} + * \text{nonlin_op}(\theta) + M\,\boldsymbol{\phi} \\ + * D\,\boldsymbol{\phi} - Q\,\theta * \end{pmatrix} - * = \mathbf{0}, - * \f] - * with - * \f[ - * M = \int \nabla\psi^\theta \;\cdot\; N^\phi \,dV,\quad - * D = \int \psi^\phi \;\cdot\; N^\phi \,dV, - * \quad - * Q = \int \psi^\phi \;\cdot\; \nabla N^\theta \,dV, - * \quad - * f(\theta) = \int \psi^\theta \;\cdot\; \theta^n \,dV. + * = \mathbf{0} * \f] + * (The exact structure depends on the `PolytropeOperator`'s internal assembly, which might involve stabilization terms using S). * - * @note MFEM’s MixedVectorWeakDivergenceIntegrator implements - * \f$ -\nabla\!\cdot \f$ so we supply a –1 coefficient to make - * `Mform` represent the +M from the derivation. The single negation - * in `PolytropeOperator` then restores the final block sign. - * - - * - * @pre \c m_feTheta and \c m_fePhi must be valid, populated FiniteElementSpaces. - * @post \c m_polytropOperator is constructed with assembled forms and offsets. - * + * @pre `m_feTheta` and `m_fePhi` must be valid, populated `mfem::FiniteElementSpace` objects. + * `m_polytropicIndex` must be set. + * @post The returned `formBundle` contains unique pointers to assembled (and finalized where appropriate) MFEM forms. */ std::unique_ptr buildIndividualForms(const mfem::Array& blockOffsets); /** - * @brief Assemble and finalize the form (Must be a form that can be assembled and finalized) + * @brief Assemble and finalize a given MFEM form. * - * @param f form which is to be assembled and finalized + * This template function calls `Assemble()` and `Finalize()` on the provided form. + * `Assemble()` computes local element matrices and adds them to the global matrix. + * `Finalize()` builds the sparsity pattern and allows the `SparseMatrix` representation to be extracted. * - * @pre f is a valid form that can be assembled and finalized (Such as Bilinear or MixedBilinearForm) - * @post f is assembled and finalized + * @tparam FormType The type of the MFEM form (e.g., `mfem::BilinearForm`, `mfem::MixedBilinearForm`). + * @param f A reference to the form to be assembled and finalized. + * + * @pre `f` must be a valid form object that supports `Assemble()` and `Finalize()` methods. + * @post `f` is assembled and finalized, and its sparse matrix representation is available. + * + * @note Nonlinear forms like `mfem::NonlinearForm` are typically not "finalized" in this sense, + * as they don't directly produce a sparse matrix but are evaluated. This function is + * intended for linear forms. */ static void assembleAndFinalizeForm(auto &f); + + /** + * @brief Computes the essential true degrees of freedom (DOFs) and their values for boundary conditions. + * + * This method determines the DOFs that correspond to: + * 1. \f$\theta = 1\f$ at the center of the star. + * 2. \f$\theta = 0\f$ on the surface of the star (boundary attribute 1). + * 3. \f$\boldsymbol{\phi} \cdot \mathbf{n} = \text{surface\_flux}\f$ on the surface, where the flux is derived from \f$\theta'(\xi_1)\f$. + * (This is applied to the normal component of \f$\boldsymbol{\phi}\f$). + * 4. Potentially \f$\boldsymbol{\phi} = \mathbf{0}\f$ at the center (though the current implementation in `polySolver.cpp` + * seems to set components of \f$\boldsymbol{\phi}\f$ to zero at DOFs associated with the center element(s)). + * + * @return An `SSE::MFEMArrayPairSet` containing two pairs: + * - The first pair is for \f$\theta\f$: (`mfem::Array` of essential TDOF indices, `mfem::Array` of corresponding values). + * - The second pair is for \f$\boldsymbol{\phi}\f$: (`mfem::Array` of essential TDOF indices, `mfem::Array` of corresponding values). + * + * @details "True DOFs" (tdof) in MFEM refer to the actual degrees of freedom in the global system, + * as opposed to local DOFs within an element. Essential boundary conditions fix the values of these DOFs. + * The center condition for \f$\theta\f$ is applied to DOFs identified by `findCenterElement()`. + * Surface conditions are applied to boundary attributes marked as essential (typically attribute 1). + * The surface flux for \f$\boldsymbol{\phi}\f$ is obtained from `polycoeff::thetaSurfaceFlux(m_polytropicIndex)`. + * + * @pre `m_mesh`, `m_feTheta`, `m_fePhi` must be initialized. `m_polytropicIndex` must be set. + * @throw `mfem::ErrorException` (via `MFEM_ABORT`) if `findCenterElement()` fails to locate the center vertex. + */ SSE::MFEMArrayPairSet getEssentialTrueDof() const; + + /** + * @brief Finds the degrees of freedom (DOFs) associated with the geometric center (origin) of the mesh. + * + * This method identifies: + * 1. The mesh vertex located at coordinates (0,0,0). + * 2. The H1 DOFs (\f$\theta\f$) associated with this center vertex. + * 3. The RT DOFs (\f$\boldsymbol{\phi}\f$) associated with mesh elements that share this center vertex. + * + * These DOFs are used to apply boundary conditions at the center of the polytrope, + * such as \f$\theta(\xi_c)=1\f$ and \f$\boldsymbol{\phi}(\xi_c)=\mathbf{0}\f$. + * + * @return A `std::pair` of `mfem::Array`: + * - `first`: Array of \f$\theta\f$ (H1) DOF indices at the center. + * - `second`: Array of \f$\boldsymbol{\phi}\f$ (RT) DOF indices associated with central elements. + * + * @throw `mfem::ErrorException` (via `MFEM_ABORT`) if no vertex is found at the origin (within tolerance). + * + * @details For RT elements, DOFs are typically associated with faces (or edges in 2D). This method collects + * DOFs from all elements connected to the center vertex. For H1 elements, DOFs can be directly + * associated with vertices. + * + * @pre `m_mesh`, `m_feTheta`, `m_fePhi` must be initialized. + */ std::pair, mfem::Array> findCenterElement() const; + + /** + * @brief Sets the initial guess for the solution variables \f$\theta\f$ and \f$\boldsymbol{\phi}\f$. + * + * - For \f$\theta\f$: + * - The interior initial guess is based on the Lane-Emden series expansion: + * \f$\theta(\xi) = \sum a_k \xi^k\f$ (using `laneEmden::thetaSeriesExpansion`). + * - The boundary condition \f$\theta(\xi_1) = 0\f$ is projected onto the surface. + * - For \f$\boldsymbol{\phi}\f$: + * - The initial guess is \f$\boldsymbol{\phi} = -\nabla \theta_{\text{init}}\f$, where \f$\theta_{\text{init}}\f$ is the initial guess for \f$\theta\f$. + * - The boundary condition for the normal component \f$\boldsymbol{\phi} \cdot \mathbf{n}\f$ at the surface is projected. + * This value is \f$\theta'(\xi_1)\f$, obtained from `polycoeff::thetaSurfaceFlux`. + * - \f$\boldsymbol{\phi}\f$ components corresponding to central DOFs (from `findCenterElement`) are set to 0. + * + * The method uses `mfem::GridFunction::ProjectCoefficient` and `ProjectBdrCoefficient` (or `ProjectBdrCoefficientNormal`) + * for these projections. + * + * @note The projection of boundary conditions for \f$\boldsymbol{\phi}\f$ (RT elements) might not be exact due to + * H(div) continuity requirements. Inhomogeneous BC enforcement in `PolytropeOperator` handles this more robustly. + * + * @pre `m_theta`, `m_phi`, `m_mesh`, `m_feTheta`, `m_fePhi` must be initialized. + * `m_polytropicIndex` must be set. + * Configuration settings for viewing initial guess should be loaded if desired. + * @post `m_theta` and `m_phi` grid functions are populated with the initial guess. + */ void setInitialGuess() const; + + /** + * @brief Saves the 1D radial solution and optionally displays the 2D/3D solution using GLVis. + * + * Extracts the \f$\theta\f$ and \f$\boldsymbol{\phi}\f$ components from the converged `state_vector`. + * - If configured (`Poly:Output:1D:Save`), it extracts a 1D solution along a specified ray + * (defined by co-latitude and longitude) using `Probe::getRaySolution` and saves it to a CSV file. + * - If configured (`Poly:Output:View`), it displays \f$\theta\f$ and \f$\boldsymbol{\phi}\f$ using `Probe::glVisView`. + * + * @param state_vector The full solution vector (block vector containing \f$\theta\f$ and \f$\boldsymbol{\phi}\f$) + * obtained from the Newton solver, typically after reconstruction by `PolytropeOperator`. + * + * @pre `m_polytropOperator`, `m_feTheta`, `m_fePhi` must be initialized. + * Configuration settings for output and viewing must be loaded. + * @post Solution data is saved and/or visualized according to configuration. + */ void saveAndViewSolution(const mfem::BlockVector& state_vector) const; + + /** + * @brief Sets up the Newton solver and its associated linear solver (GMRES). + * + * 1. Loads solver parameters (tolerances, max iterations, print levels) for both Newton + * and GMRES from the configuration file (via `LoadSolverUserParams`). + * 2. Creates a `solverBundle` to manage the lifetimes of `mfem::GMRESSolver` and `mfem::NewtonSolver`. + * 3. Configures the GMRES solver with its parameters. + * 4. Configures the Newton solver with its parameters, sets the `PolytropeOperator` as the + * nonlinear system operator, and sets the configured GMRES solver as the linear solver + * for inverting Jacobians. + * + * @return A `solverBundle` struct containing the configured Newton and GMRES solvers. + * The ownership of the solvers within the bundle is managed by the bundle itself. + * + * @pre `m_polytropOperator` must be initialized and finalized. + * Configuration settings for solver parameters must be available. + * @post A fully configured `solverBundle` is returned, ready for `newton.Mult()`. + */ solverBundle setupNewtonSolver() const; + + /** + * @brief Sets the essential true DOFs on the `PolytropeOperator`. + * + * Calls `getEssentialTrueDof()` to compute the boundary condition DOFs and values, + * and then passes them to `m_polytropOperator->set_essential_true_dofs()`. + * This step is crucial for the `PolytropeOperator` to correctly handle + * boundary conditions when forming reduced systems and applying residuals/Jacobians. + * + * @pre `m_polytropOperator` must be initialized. + * `m_mesh`, `m_feTheta`, `m_fePhi`, `m_polytropicIndex` must be set (for `getEssentialTrueDof`). + * @post Essential boundary conditions are registered with the `m_polytropOperator`. + * The `PolytropeOperator` will likely be marked as not finalized after this call. + */ void setOperatorEssentialTrueDofs() const; + /** + * @brief Loads Newton and GMRES solver parameters from the configuration. + * + * Retrieves relative tolerance, absolute tolerance, maximum iterations, and print level + * for both the Newton solver and the GMRES linear solver from the `Config` instance. + * Default values are used if specific configuration keys are not found. + * + * Configuration keys are typically prefixed with `Poly:Solver:Newton:` and `Poly:Solver:GMRES:`. + * Example keys: `Poly:Solver:Newton:RelTol`, `Poly:Solver:GMRES:MaxIter`. + * + * @param[out] newtonRelTol Relative tolerance for Newton solver. + * @param[out] newtonAbsTol Absolute tolerance for Newton solver. + * @param[out] newtonMaxIter Maximum iterations for Newton solver. + * @param[out] newtonPrintLevel Print level for Newton solver. + * @param[out] gmresRelTol Relative tolerance for GMRES solver. + * @param[out] gmresAbsTol Absolute tolerance for GMRES solver. + * @param[out] gmresMaxIter Maximum iterations for GMRES solver. + * @param[out] gmresPrintLevel Print level for GMRES solver. + * + * @pre `m_config` must be a valid reference to a `Config` object. + * `m_logger` should be initialized if debug logging is enabled. + * @post Output parameters are populated with values from configuration or defaults. + * Solver parameters are logged at DEBUG level. + */ void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel, double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const; + /** + * @brief Utility function to get and save the coordinates of degrees of freedom for a finite element space. + * + * For each element in the mesh: + * 1. Gets the DOFs associated with that element from the given `fes`. + * 2. Gets the physical coordinates of the element's center (or a reference point). + * 3. Writes the DOF indices, mesh radius, element center radius, and element center x,y,z coordinates + * to the specified output file in CSV format. + * + * This can be useful for debugging or analyzing the distribution of DOFs. + * + * @param fes The `mfem::FiniteElementSpace` for which to get DOF coordinates. + * @param filename The name of the output CSV file. + * + * @pre `fes` must be a valid, initialized `mfem::FiniteElementSpace`. + * The mesh associated with `fes` must be valid. + * @post A CSV file named `filename` is created/truncated and populated with DOF coordinate information. + * + * @note For DOFs shared by multiple elements, this function might list them multiple times, + * associated with each element they belong to. The output format lists all DOFs for an element + * on one line, pipe-separated, followed by coordinate data for that element. + */ static void GetDofCoordinates(const mfem::FiniteElementSpace &fes, const std::string& filename); }; + diff --git a/src/poly/utils/public/polytropeOperator.h b/src/poly/utils/public/polytropeOperator.h index 828214b..b8ba482 100644 --- a/src/poly/utils/public/polytropeOperator.h +++ b/src/poly/utils/public/polytropeOperator.h @@ -310,8 +310,8 @@ private: std::unique_ptr m_ScaledSReduced; ///< Scaled S matrix (free DOFs only, scaled by stabilization coefficient). - // --- Stabilization Coefficient --- (Perhapses this should be a user parameter...) - static constexpr double m_stabilizationCoefficient = -1.0; ///< Stabilization coefficient for the system, used to more tightly couple ∇θ and φ. + // --- Stabilization Coefficient --- (Perhapses this should be a user parameter...) // TODO: Sort out why this is negative (sign error?) + static constexpr double m_stabilizationCoefficient = -2.0; ///< Stabilization coefficient for the system, used to more tightly couple ∇θ and φ. static constexpr double m_IncrementedStabilizationCoefficient = 1.0 + m_stabilizationCoefficient; ///< 1 + Stabilization coefficient for the system, used to more tightly couple ∇θ and φ. // --- State Vectors and DOF Management ---