feat(poly): locking phi surface flux and fixed phi boundary condition application

src/poly/solver/private/polySolver.cpp

@@ -156,7 +156,7 @@ std::unique_ptr<formBundle> PolySolver::buildIndividualForms(const mfem::Array<i
 
     forms->f->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(m_polytropicIndex));
 
-    return std::move(forms);
+    return forms;
 }
 
 void PolySolver::assembleAndFinalizeForm(auto &f) {
@@ -264,11 +264,11 @@ void PolySolver::setInitialGuess() const {
         [this](const mfem::Vector &x) {
             const double r = x.Norml2();
             const double radius = Probe::getMeshRadius(*m_mesh);
-            const double u = 1/radius;
+            // const double u = 1/radius;
 
             // return (-1.0/radius) * r + 1;
             // return -std::pow((u*r), 2)+1.0; // The series expansion is a better guess; however, this is cheaper and ensures that the value at the surface is very close to zero in a way that the series expansion does not
-            return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 5);
+            return laneEmden::thetaSeriesExpansion(r, m_polytropicIndex, 10);
         }
     );
 
@@ -399,7 +399,8 @@ solverBundle PolySolver::setupNewtonSolver() const {
     solver.solver.SetMaxIter(gmresMaxIter);
     solver.solver.SetPrintLevel(gmresPrintLevel);
 
-    solver.solver.SetPreconditioner(m_polytropOperator->GetPreconditioner());
+    // Preconditioner turned off because the polytrope operator seems *very* well conditioned without it
+    // solver.solver.SetPreconditioner(m_polytropOperator->GetPreconditioner());
     // --- Set up the Newton solver ---
     solver.newton.SetRelTol(newtonRelTol);
     solver.newton.SetAbsTol(newtonAbsTol);

src/poly/utils/meson.build

@@ -30,6 +30,7 @@ dependencies = [
   quill_dep,
   config_dep,
   types_dep,
+  mfemanalysis_dep,
 ]
 
 libpolyutils = static_library('polyutils',

src/poly/utils/private/operator.cpp

@@ -21,70 +21,13 @@
 #include "operator.h"
 #include "4DSTARTypes.h"
 #include "mfem.hpp"
+#include "mfem_smout.h"
 #include <memory>
 
+#include "../../../../utils/debugUtils/MFEMAnalysisUtils/MFEMAnalysis-cpp/src/include/mfem_smout.h"
+
 static int s_PolyOperatorMultCalledCount = 0;
-
+static int s_PolyOperatorGradCalledCount = 0;
-void writeDenseMatrixToCSV(const std::string &filename, int precision, const mfem::DenseMatrix *mat) {
-    if (!mat) {
-        throw std::runtime_error("The operator is not a SparseMatrix.");
-    }
-
-    std::ofstream outfile(filename);
-    if (!outfile.is_open()) {
-        throw std::runtime_error("Failed to open file: " + filename);
-    }
-
-
-    int height = mat->Height();
-    int width = mat->Width();
-
-    // Set precision for floating-point output
-    outfile << std::fixed << std::setprecision(precision);
-
-    for (int i = 0; i < width; i++) {
-        outfile << i;
-        if (i < width - 1) {
-            outfile << ",";
-        }
-        else {
-            outfile << "\n";
-        }
-    }
-
-    // Iterate through rows
-    for (int i = 0; i < height; ++i) {
-        for (int j = 0; j < width; ++j) {
-            outfile << mat->Elem(i, j);
-            if (j < width - 1) {
-                outfile << ",";
-            }
-        }
-        outfile << std::endl;
-    }
-
-    outfile.close();
-}
-
-/**
- * @brief Writes the dense representation of an MFEM Operator (if it's a SparseMatrix) to a CSV file.
- *
- * @param op The MFEM Operator to write.
- * @param filename The name of the output CSV file.
- * @param precision Number of decimal places for floating-point values.
- */
- void writeOperatorToCSV(const mfem::Operator &op,
-    const std::string &filename,
-    const int precision = 6) // Add precision argument
-{
-    // Attempt to cast the Operator to a SparseMatrix
-    const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&op);
-    if (!sparse_mat) {
-        throw std::runtime_error("The operator is not a SparseMatrix.");
-    }
-    const mfem::DenseMatrix *mat = sparse_mat->ToDenseMatrix();
-    writeDenseMatrixToCSV(filename, precision, mat);
-}
 
 void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::SparseMatrix>& invMat, const std::string& name="matrix") {
     // PERF: This likely can be made much more efficient and will probably be called in tight loops, a good
@@ -138,14 +81,32 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
         return; // do nothing if already finalized
     }
 
-    m_negM_op = std::make_unique<mfem::ScaledOperator>(m_M.get(), -1.0);
+    m_Mmat = std::make_unique<mfem::SparseMatrix>(m_M->SpMat());
-    m_negQ_op = std::make_unique<mfem::ScaledOperator>(m_Q.get(), -1.0);
+    m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
+    m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
+
+    for (const auto& dof: m_theta_ess_tdofs.first) {
+        m_Mmat->EliminateRow(dof);
+        m_Qmat->EliminateCol(dof);
+    }
+
+    for (const auto& dof: m_phi_ess_tdofs.first) {
+        m_Mmat->EliminateCol(dof);
+        m_Qmat->EliminateRow(dof);
+        m_Dmat->EliminateRowCol(dof);
+    }
+
+    // m_negM_op = std::make_unique<mfem::ScaledOperator>(m_Mmat.get(), -1.0);
+    // m_negQ_op = std::make_unique<mfem::ScaledOperator>(m_Qmat.get(), -1.0);
+
+    m_Mmat->operator*=(-1.0);
+    m_Qmat->operator*=(-1.0);
 
     // Set up the constant parts of the jacobian now
     m_jacobian = std::make_unique<mfem::BlockOperator>(m_blockOffsets);
-    m_jacobian->SetBlock(0, 1, m_negM_op.get()); //<- -M (constant)
+    m_jacobian->SetBlock(0, 1, m_Mmat.get()); //<- -M (constant)
-    m_jacobian->SetBlock(1, 0, m_negQ_op.get()); //<- -Q (constant)
+    m_jacobian->SetBlock(1, 0, m_Qmat.get()); //<- -Q (constant)
-    m_jacobian->SetBlock(1, 1, m_D.get());       //<-  D (constant)
+    m_jacobian->SetBlock(1, 1, m_Dmat.get());       //<-  D (constant)
 
     // Build the initial preconditioner based on some initial guess
     const auto &grad = m_f->GetGradient(initTheta);
@@ -214,28 +175,20 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     // -- Apply essential boundary conditions --
     for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
         if (int idx = m_theta_ess_tdofs.first[i]; idx >= 0 && idx < y_R0.Size()) {
-            const double &targetValue = m_theta_ess_tdofs.second[i];
+            // const double &targetValue = m_theta_ess_tdofs.second[i];
             // y_block.GetBlock(0)[idx] = targetValue - x_theta(idx); // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
-            y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form
+            y_block.GetBlock(0)[idx] = 0; // Zero out the essential theta dofs in the bi-linear form // TODO Check if this is double zeroing (i.e if they were already removed maybe I am removing something that should not be removed here)
+        }
+    }
+    std::cout << "θ Norm: " << y_block.GetBlock(0).Norml2() << std::endl;
+
+    for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
+        if (int idx = m_phi_ess_tdofs.first[i]; idx >= 0 && idx < y_R1.Size()) {
+            y_block.GetBlock(1)[idx] = 0; // Zero out the essential phi dofs in the bi-linear form // TODO Check if this is double zeroing (i.e if they were already removed maybe I am removing something that should not be removed here)
         }
     }
 
-    // TODO: Are the residuals for φ being calculated correctly?
+    std::cout << "φ Norm: " << y_block.GetBlock(1).Norml2() << std::endl;
-
-    std::ofstream outfileθ("PolyOperatorMultCallTheta_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc);
-    outfileθ << "dof,R_θ\n";
-    for (int i = 0; i < y_R0.Size(); i++) {
-        outfileθ << i << "," << y_R0(i) << "\n";
-    }
-    outfileθ.close();
-    std::ofstream outfileφ("PolyOperatorMultCallPhi_" + std::to_string(s_PolyOperatorMultCalledCount) + ".csv", std::ios::out | std::ios::trunc);
-    outfileφ << "dof,R_ɸ\n";
-    for (int i = 0; i < y_R1.Size(); i++) {
-        outfileφ << i << "," << y_R1(i) << "\n";
-    }
-    outfileφ.close();
-    ++s_PolyOperatorMultCalledCount;
-
 }
 
 
@@ -270,6 +223,7 @@ void PolytropeOperator::updateInverseSchurCompliment() const {
     // ⎣0     S^-1⎦
     m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get());
     m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get());
+
 }
 
 void PolytropeOperator::updatePreconditioner(const mfem::Operator &grad) const {
@@ -287,11 +241,10 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
     const mfem::Vector& x_theta = x_block.GetBlock(0);
 
     auto &grad = m_f->GetGradient(x_theta);
-    updatePreconditioner(grad);
+    // auto *gradPtr = &grad;
+    // updatePreconditioner(grad);
 
     m_jacobian->SetBlock(0, 0, &grad);
-    // The other blocks are set up in finalize since they are constant. Only J00 depends on the current state of theta
-
     return *m_jacobian;
 }
 
@@ -305,26 +258,6 @@ void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess
     } else {
         MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs");
     }
-
-    if (m_M) {
-        m_M->EliminateTestEssentialBC(theta_ess_tdofs.first);
-        m_M->EliminateTrialEssentialBC(phi_ess_tdofs.first);
-    } else {
-        MFEM_ABORT("m_M is null in PolytropeOperator::SetEssentialTrueDofs");
-    }
-
-    if (m_Q) {
-        m_Q->EliminateTrialEssentialBC(theta_ess_tdofs.first);
-        m_Q->EliminateTestEssentialBC(phi_ess_tdofs.first);
-    } else {
-        MFEM_ABORT("m_Q is null in PolytropeOperator::SetEssentialTrueDofs");
-    }
-
-    if (m_D) {
-        m_D->EliminateEssentialBC(phi_ess_tdofs.first);
-    } else {
-        MFEM_ABORT("m_D is null in PolytropeOperator::SetEssentialTrueDofs");
-    }
 }
 
 void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) {

src/poly/utils/public/operator.h

@@ -70,6 +70,7 @@ private:
     SSE::MFEMArrayPair m_theta_ess_tdofs;
     SSE::MFEMArrayPair m_phi_ess_tdofs;
 
+    std::unique_ptr<mfem::SparseMatrix> m_Mmat, m_Qmat, m_Dmat;
     std::unique_ptr<mfem::ScaledOperator> m_negM_op;
     std::unique_ptr<mfem::ScaledOperator> m_negQ_op;
     mutable std::unique_ptr<mfem::BlockOperator> m_jacobian;