fix(operator): changed MULT (residuals) to deal with negated M properly

src/poly/utils/private/operator.cpp

@@ -26,9 +26,6 @@
 
 #include "../../../../utils/debugUtils/MFEMAnalysisUtils/MFEMAnalysis-cpp/src/include/mfem_smout.h"
 
-static int s_PolyOperatorMultCalledCount = 0;
-static int s_PolyOperatorGradCalledCount = 0;
-
 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
     // PERF:   place for some easy optimization might be here.
@@ -99,7 +96,8 @@ void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
     // 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);
+    // TODO Replace this with a scaled operator when I am done writing these out to disk
+    m_Mmat->operator*=(-1.0); // Sparse matrix only defines an inplace operator*= not an operator* so we need to do this
     m_Qmat->operator*=(-1.0);
 
     // Set up the constant parts of the jacobian now
@@ -165,7 +163,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     MFEM_ASSERT(m_f.get() != nullptr, "NonlinearForm m_f is null in PolytropeOperator::Mult");
 
     m_f->Mult(x_theta, f_term); // fixme: this may be the wrong way to assemble m_f?
-    m_M->Mult(x_phi, Mphi_term);
+    m_M->Mult(x_phi, Mphi_term); // Does the order of operations work out between this and the next subtract to make Mφ negative properly? I think so but double check
     m_D->Mult(x_phi, Dphi_term);
     m_Q->Mult(x_theta, Qtheta_term);
 
@@ -175,16 +173,16 @@ 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];
-            // 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 // 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)
+            const double &targetValue = m_theta_ess_tdofs.second[i];
+            y_block.GetBlock(0)[idx] = x_theta(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
         }
     }
     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)
+            const double &targetValue = m_phi_ess_tdofs.second[i];
+            y_block.GetBlock(1)[idx] = x_phi(idx) - targetValue; // inhomogenous essential bc. This is commented out since seems it results in dramatic instabilies arrising
         }
     }
 
@@ -205,14 +203,15 @@ void PolytropeOperator::updateInverseSchurCompliment() const {
     if (m_invNonlinearJacobian == nullptr) {
         MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian");
     }
-    mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // This is now a copy of D
+    mfem::SparseMatrix* schurCompliment(&m_D->SpMat()); // Represents S in the preconditioner, starts as a copy of D before being modified in place
 
     // Calculate S = D - Q df^{-1} M
-    mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(m_Q->SpMat(), *m_invNonlinearJacobian); // Q * df^{-1}
-    const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, m_M->SpMat()); // Q * df^{-1} * M
+    mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(*m_Qmat, *m_invNonlinearJacobian); // Q * df^{-1}
+    const mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, *m_Mmat); // Q * df^{-1} * M
 
     // PERF: Could potentially add some caching in here to only update the preconditioner when some condition has been met
-    schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M
+    schurCompliment->Add(1, *temp_QxdfInvxM); // D - Q * df^{-1} * M
+    saveSparseMatrixBinary(*schurCompliment, "schurCompliment.bin");
     approxJacobiInvert(*schurCompliment, m_invSchurCompliment, "Schur Compliment");
 
     if (m_schurPreconditioner == nullptr) {