fix(poly): added stabilization term

least squares stabalization term seems to have stabilized polytrope and mostly resolve the overshoot mode (in more non linear cases (n > 2) the mode does reapear; however, it is much less pronounced

src/poly/solver/private/polySolver.cpp

@@ -114,24 +114,25 @@ void PolySolver::assembleBlockSystem() {
 
     const std::unique_ptr<formBundle> forms = buildIndividualForms(blockOffsets);
 
-    const double penalty_param = m_config.get<double>("Poly::Solver::ZeroDerivativePenalty", 1.0);
+    // const double penalty_param = m_config.get<double>("Poly::Solver::ZeroDerivativePenalty", 1.0);
-    mfem::Array<int> thetaCenterDofs, phiCenterDofs;
+    // mfem::Array<int> thetaCenterDofs, phiCenterDofs;
-    std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
+    // std::tie(thetaCenterDofs, phiCenterDofs) = findCenterElement();
-    mfem::SparseMatrix& D_mat = forms->D->SpMat();
+    // mfem::SparseMatrix& D_mat = forms->D->SpMat();
-
+    //
-    for (int i = 0; i < phiCenterDofs.Size(); ++i)
+    // for (int i = 0; i < phiCenterDofs.Size(); ++i)
-    {
+    // {
-        const int dof_idx = phiCenterDofs[i];
+    //     const int dof_idx = phiCenterDofs[i];
-        if (dof_idx >= 0 && dof_idx < D_mat.Height()) {
+    //     if (dof_idx >= 0 && dof_idx < D_mat.Height()) {
-            D_mat(dof_idx, dof_idx) += penalty_param;
+    //         D_mat(dof_idx, dof_idx) += penalty_param;
-        }
+    //     }
-    }
+    // }
 
     // --- Build the BlockOperator ---
     m_polytropOperator = std::make_unique<PolytropeOperator>(
         std::move(forms->M),
         std::move(forms->Q),
         std::move(forms->D),
+        std::move(forms->S),
         std::move(forms->f),
         blockOffsets);
 }
@@ -152,26 +153,26 @@ std::unique_ptr<formBundle> PolySolver::buildIndividualForms(const mfem::Array<i
         std::make_unique<mfem::MixedBilinearForm>(m_fePhi.get(), m_feTheta.get()),
         std::make_unique<mfem::MixedBilinearForm>(m_feTheta.get(), m_fePhi.get()),
         std::make_unique<mfem::BilinearForm>(m_fePhi.get()),
+        std::make_unique<mfem::BilinearForm>(m_feTheta.get()),
         std::make_unique<mfem::NonlinearForm>(m_feTheta.get())
     );
 
-    // --- -M negation -> M ---
-    mfem::Vector negOneVec(m_mesh.SpaceDimension());
-    negOneVec = -1.0;
-
-    m_negationCoeff = std::make_unique<mfem::VectorConstantCoefficient>(negOneVec);
-
     // --- Add the integrators to the forms ---
-    forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator());
+    forms->M->AddDomainIntegrator(new mfem::MixedVectorWeakDivergenceIntegrator()); // M ∫∇ψ^θ·N^φ dV
-    forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator());
+    forms->Q->AddDomainIntegrator(new mfem::MixedVectorGradientIntegrator());       // Q ∫ψ^φ·∇N^θ dV
-    forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator());
+    forms->D->AddDomainIntegrator(new mfem::VectorFEMassIntegrator());              // D ∫ψ^φ·N^φ dV
+    forms->S->AddDomainIntegrator(new mfem::DiffusionIntegrator());                 // S ∫∇ψ^θ·∇N^θ dV
+    forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex));
 
     // --- Assemble and Finalize the forms ---
     assembleAndFinalizeForm(forms->M);
     assembleAndFinalizeForm(forms->Q);
     assembleAndFinalizeForm(forms->D);
+    assembleAndFinalizeForm(forms->S);
+
+    // Note: The NonlinearForm does not need to / cannot be finalized, as it is not a matrix form. Rather, the operator
+    //       will evaluate the nonlinear form during the solve phase.
 
-    forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex));
 
     return forms;
 }
@@ -179,7 +180,7 @@ std::unique_ptr<formBundle> PolySolver::buildIndividualForms(const mfem::Array<i
 void PolySolver::assembleAndFinalizeForm(auto &f) {
     // This constructs / ensures the matrix representation for each form
     // Assemble => Computes the local element matrices across the domain. Adds these to the global matrix . Adds these to the global matrix.
-    // Finalize => Builds the sparsity pattern and allows the SparseMatrix representation to be extracted.
+    // Finalize => Builds the sparsity pattern and allows the SparseMatrix representation to be extracted (CSR).
     f->Assemble();
     f->Finalize();
 }

src/poly/solver/public/polySolver.h

@@ -49,6 +49,7 @@ struct formBundle {
     std::unique_ptr<mfem::MixedBilinearForm> M; //<-- M ∫∇ψ^θ·N^φ dV
     std::unique_ptr<mfem::MixedBilinearForm> Q; //<-- Q ∫ψ^φ·∇N^θ dV
     std::unique_ptr<mfem::BilinearForm> D; //<-- D ∫ψ^φ·N^φ dV
+    std::unique_ptr<mfem::BilinearForm> S; //<-- S ∫∇ψ^θ·∇N^θ dV
     std::unique_ptr<mfem::NonlinearForm> f; //<-- f(θ) ∫ψ^θ·θ^n dV
 };
 
@@ -84,8 +85,6 @@ private: // Private Attributes
 
     std::unique_ptr<mfem::OperatorJacobiSmoother> m_prec;
 
-    std::unique_ptr<mfem::VectorConstantCoefficient> m_negationCoeff;
-
 
 private: // Private methods
     PolySolver(mfem::Mesh& mesh, double n, double order);

src/poly/utils/private/polytropeOperator.cpp

@@ -134,10 +134,10 @@ PolytropeOperator::PolytropeOperator(
   std::unique_ptr<mfem::MixedBilinearForm> M,
   std::unique_ptr<mfem::MixedBilinearForm> Q,
   std::unique_ptr<mfem::BilinearForm> D,
+  std::unique_ptr<mfem::BilinearForm> S,
   std::unique_ptr<mfem::NonlinearForm> f,
   const mfem::Array<int> &blockOffsets) :
 
-  // TODO: Need to update this so that the size is that of the reduced system operator
   mfem::Operator(blockOffsets.Last()), // Initialize the base class with the total size of the block offset vector
   m_blockOffsets(blockOffsets),
   m_jacobian(nullptr) {
@@ -145,6 +145,7 @@ PolytropeOperator::PolytropeOperator(
     m_M = std::move(M);
     m_Q = std::move(Q);
     m_D = std::move(D);
+    m_S = std::move(S);
     m_f = std::move(f);
 
     // Use Gauss-Seidel smoother to approximate the inverse of the matrix
@@ -181,21 +182,33 @@ void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) con
     mfem::Vector Mphi_term(theta_size);
     mfem::Vector Dphi_term(phi_size);
     mfem::Vector Qtheta_term(phi_size);
+    mfem::Vector Stheta_term(theta_size);
+
 
     // Calculate R0 and R1 terms
-    // R0 = f(θ) + Mɸ
+    // R0 = f(θ) + (1+c)Mɸ + cSθ
-    // R1 = Dɸ - Qθ
+    // R1 = (1+c)Dɸ - (1+c)Qθ
 
     MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult");
 
-    m_f->Mult(x_theta, f_term);
+    m_f->Mult(x_theta, f_term);      // f(θ)
-    m_M->Mult(x_phi, Mphi_term);
+    m_M->Mult(x_phi, Mphi_term);     // Mɸ
-    m_D->Mult(x_phi, Dphi_term);
+    m_D->Mult(x_phi, Dphi_term);     // Dɸ
-    m_Q->Mult(x_theta, Qtheta_term);
+    m_Q->Mult(x_theta, Qtheta_term); // Qθ
+    m_S->Mult(x_theta, Stheta_term); // Sθ
 
-    add(f_term, Mphi_term, y_R0);
+    // PERF: these multiplications may be able to be refactored into the matrices so they only need to be done once.
-    subtract(Dphi_term, Qtheta_term, y_R1);
+    Stheta_term *= m_stabilizationCoefficient;            // cSθ
+    Qtheta_term *= m_IncrementedStabilizationCoefficient; // (1+c)Qθ
+    Mphi_term *= m_IncrementedStabilizationCoefficient;   // (1+c)Mɸ
+    Dphi_term *= m_IncrementedStabilizationCoefficient;   // (1+c)Dɸ
 
+    mfem::Vector y_R0_temp(theta_size);
+    add(f_term, Mphi_term, y_R0_temp);      // R0 = f(θ) + (1+c)Mɸ
+    add(y_R0_temp, Stheta_term, y_R0);      // R0 = f(θ) + (1+c)Mɸ + cSθ
+    subtract(Dphi_term, Qtheta_term, y_R1); // R1 = (1+c)Dɸ - (1+c)Qθ
+
+    // --- Scatter the residual vector y into the free dofs ---
     yFree.SetSize(m_reducedBlockOffsets.Last());
     MFEM_ASSERT(m_freeDofs.Size() == m_reducedBlockOffsets.Last(), "PolytropeOperator::Mult: Size of free dofs does not match reduced block offsets size.");
     for (int i = 0, j = 0; i < y.Size(); ++i) {
@@ -207,7 +220,6 @@ void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) con
 }
 
 mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const {
-    //TODO: This now needs to be updated to deal with the reduced system size
     if (!m_isFinalized) {
         MFEM_ABORT("PolytropeOperator::GetGradient called before finalize");
     }
@@ -219,7 +231,7 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const
     // PERF: There are a lot of copies and loops here, probably performance could be gained by flattering some of these.
     auto &grad = m_f->GetGradient(x_theta);
     // updatePreconditioner(grad);
-    const auto gradMatrix = dynamic_cast<mfem::SparseMatrix*>(&grad);
+    const auto gradMatrix = dynamic_cast<mfem::SparseMatrix*>(&grad); // ∂f(θ)/∂θ
 
     if (gradMatrix == nullptr) {
         MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix.");
@@ -231,7 +243,11 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const
             m_theta_ess_tdofs.first,
             m_theta_ess_tdofs.first
             )
-        );
+        ); // ∂f(θ)/∂θ (now reduced to only the free DOFs)
+
+
+    // TODO: Confirm this actually does what I want it to do
+    *m_gradReduced += *m_ScaledSReduced; // ∂f(θ)/∂θ + cS (reduced to only the free DOFs)
 
     m_jacobian->SetBlock(0, 0, m_gradReduced.get());
 
@@ -346,29 +362,73 @@ void PolytropeOperator::populate_free_dof_array() {
 
 // PolytropeOperator: Private Helper Methods - Construction and Setup
 void PolytropeOperator::construct_matrix_representations() {
+    // --- Construct the sparse matrix representations of M, Q, D, and S ---
     m_Mmat = std::make_unique<mfem::SparseMatrix>(m_M->SpMat());
     m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
     m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
+    m_Smat = std::make_unique<mfem::SparseMatrix>(m_S->SpMat());
 
-    m_MReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::build_reduced_matrix(*m_Mmat, m_phi_ess_tdofs.first, m_theta_ess_tdofs.first));
+    // --- Reduce the matrices to only the free DOFs ---
-    m_QReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::build_reduced_matrix(*m_Qmat, m_theta_ess_tdofs.first, m_phi_ess_tdofs.first));
+    m_MReduced = std::make_unique<mfem::SparseMatrix>(
-    m_DReduced = std::make_unique<mfem::SparseMatrix>(serif::utilities::build_reduced_matrix(*m_Dmat, m_phi_ess_tdofs.first, m_phi_ess_tdofs.first));
+        serif::utilities::build_reduced_matrix(
+            *m_Mmat,
+            m_phi_ess_tdofs.first,
+            m_theta_ess_tdofs.first)
+    );
+    m_QReduced = std::make_unique<mfem::SparseMatrix>(
+        serif::utilities::build_reduced_matrix(
+            *m_Qmat,
+            m_theta_ess_tdofs.first,
+            m_phi_ess_tdofs.first)
+    );
+    m_DReduced = std::make_unique<mfem::SparseMatrix>(
+        serif::utilities::build_reduced_matrix(
+            *m_Dmat,
+            m_phi_ess_tdofs.first,
+            m_phi_ess_tdofs.first)
+    );
+    m_SReduced = std::make_unique<mfem::SparseMatrix>(
+        serif::utilities::build_reduced_matrix(
+            *m_Smat,
+            m_theta_ess_tdofs.first,
+            m_theta_ess_tdofs.first)
+    );
 
-    m_negQ_mat = std::make_unique<mfem::ScaledOperator>(m_QReduced.get(), -1.0);
+    // --- Scale the reduced matrices by the stabilization coefficients ---
+    mfem::SparseMatrix MScaledTemp(*m_MReduced); // Create a temporary copy of the M matrix for scaling
+    MScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the M matrix by the incremented stabilization coefficient
+    m_MScaledReduced = std::make_unique<mfem::SparseMatrix>(MScaledTemp); // Store the scaled M matrix so that it persists
+
+    mfem::SparseMatrix QScaledTemp(*m_QReduced);
+    QScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the Q matrix by the incremented stabilization coefficient
+    m_QScaledReduced = std::make_unique<mfem::SparseMatrix>(QScaledTemp); // Store the scaled Q matrix so that it persists
+
+    mfem::SparseMatrix DRScaledTemp(*m_DReduced);
+    DRScaledTemp *= m_IncrementedStabilizationCoefficient; // Scale the D matrix by the incremented stabilization coefficient
+    m_DScaledReduced = std::make_unique<mfem::SparseMatrix>(DRScaledTemp); // Store the scaled D matrix so that it persists
+
+    // Scale the S matrix by the stabilization coefficient (so that the allocations only need to be done once)
+    mfem::SparseMatrix SScaledTemp(*m_SReduced); // Create a temporary copy of the S matrix for scaling
+    SScaledTemp *= m_stabilizationCoefficient; // Scale the S matrix by the stabilization coefficient
+    m_ScaledSReduced = std::make_unique<mfem::SparseMatrix>(SScaledTemp); // Store the scaled S matrix so that it persists
+
+    // --- Create the scaled operator for -(1+c)Q ---
+    m_negQ_mat = std::make_unique<mfem::ScaledOperator>(m_QScaledReduced.get(), -1.0);
 }
 
 void PolytropeOperator::construct_reduced_block_offsets() {
     m_reducedBlockOffsets.SetSize(3);
     m_reducedBlockOffsets[0] = 0; // R0 block (theta)
     m_reducedBlockOffsets[1] = m_MReduced->Height(); // R1 block (theta)
-    m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi)
+    m_reducedBlockOffsets[2] = m_QReduced->Height() + m_reducedBlockOffsets[1]; // R2 block (phi + theta)
 }
 
 void PolytropeOperator::construct_jacobian_constant_terms() {
     m_jacobian = std::make_unique<mfem::BlockOperator>(m_reducedBlockOffsets);
-    m_jacobian->SetBlock(0, 1, m_MReduced.get());     //<-  M (constant)
+
-    m_jacobian->SetBlock(1, 0, m_negQ_mat.get());     //<- -Q (constant)
+    m_jacobian->SetBlock(0, 1, m_MScaledReduced.get());     //<-  (1+c)M (constant, reduced to only free DOFs)
-    m_jacobian->SetBlock(1, 1, m_DReduced.get());     //<-  D (constant)
+    m_jacobian->SetBlock(1, 0, m_negQ_mat.get());           //<- -(1+c)Q (constant, reduced to only free DOFs)
+    m_jacobian->SetBlock(1, 1, m_DScaledReduced.get());     //<-  (1+c)D (constant, reduced to only free DOFs)
 }
 
 void PolytropeOperator::scatter_boundary_conditions() {

src/poly/utils/public/polytropeOperator.h

@@ -152,6 +152,7 @@ public:
         std::unique_ptr<mfem::MixedBilinearForm> M,
         std::unique_ptr<mfem::MixedBilinearForm> Q,
         std::unique_ptr<mfem::BilinearForm> D,
+        std::unique_ptr<mfem::BilinearForm> S,
         std::unique_ptr<mfem::NonlinearForm> f,
         const mfem::Array<int> &blockOffsets);
 
@@ -286,19 +287,33 @@ private:
     std::unique_ptr<mfem::MixedBilinearForm> m_M; ///< Bilinear form M, coupling θ and φ.
     std::unique_ptr<mfem::MixedBilinearForm> m_Q; ///< Bilinear form Q, coupling φ and θ.
     std::unique_ptr<mfem::BilinearForm> m_D;      ///< Bilinear form D, acting on φ.
+    std::unique_ptr<mfem::BilinearForm> m_S;
     std::unique_ptr<mfem::NonlinearForm> m_f;     ///< Nonlinear form f, acting on θ.
 
     // --- Full Matrix Representations (owned, derived from forms) ---
     std::unique_ptr<mfem::SparseMatrix> m_Mmat;   ///< Sparse matrix representation of M.
     std::unique_ptr<mfem::SparseMatrix> m_Qmat;   ///< Sparse matrix representation of Q.
     std::unique_ptr<mfem::SparseMatrix> m_Dmat;   ///< Sparse matrix representation of D.
+    std::unique_ptr<mfem::SparseMatrix> m_Smat;
 
     // --- Reduced Matrix Representations (owned, after eliminating essential DOFs) ---
     std::unique_ptr<mfem::SparseMatrix> m_MReduced; ///< Reduced M matrix (free DOFs only).
     std::unique_ptr<mfem::SparseMatrix> m_QReduced; ///< Reduced Q matrix (free DOFs only).
     std::unique_ptr<mfem::SparseMatrix> m_DReduced; ///< Reduced D matrix (free DOFs only).
+    std::unique_ptr<mfem::SparseMatrix> m_SReduced; ///< Reduced S matrix (free DOFs only, used for least squares stabilization).
     mutable std::unique_ptr<mfem::SparseMatrix> m_gradReduced; ///< Reduced gradient operator (G) for the nonlinear part f(θ).
 
+    // --- Scaled Reduced Matrix Representations (owned, after eliminating essential DOFs and scaling by stabilization coefficients) ---
+    std::unique_ptr<mfem::SparseMatrix> m_MScaledReduced; ///< Scaled M matrix (free DOFs only, scaled by stabilization coefficient).
+    std::unique_ptr<mfem::SparseMatrix> m_QScaledReduced; ///< Scaled Q matrix (free DOFs only, scaled by stabilization coefficient).
+    std::unique_ptr<mfem::SparseMatrix> m_DScaledReduced; ///< Scaled D matrix (free DOFs only, scaled by stabilization coefficient).
+    std::unique_ptr<mfem::SparseMatrix> 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 φ.
+    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 ---
     mutable mfem::Vector m_state;                 ///< Full state vector [θ, φ]^T, including essential DOFs.
     mfem::Array<int> m_freeDofs;                  ///< Array of indices for free (non-essential) DOFs.