From e7ad73c0f94a8a6f40336a29da641cd67c7aa4ec Mon Sep 17 00:00:00 2001
From: Emily Boudreaux <thomas@boudreauxmail.com>
Date: Mon, 21 Apr 2025 07:33:39 -0400
Subject: [PATCH] refactor(jacobiInvert): moved all jacobi inverting logic into
 dedicated function

---
 src/poly/utils/private/operator.cpp | 81 +++++++++--------------------
 1 file changed, 26 insertions(+), 55 deletions(-)

diff --git a/src/poly/utils/private/operator.cpp b/src/poly/utils/private/operator.cpp
index 7d42500..5432d18 100644
--- a/src/poly/utils/private/operator.cpp
+++ b/src/poly/utils/private/operator.cpp
@@ -125,50 +125,6 @@ const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const {
      return *m_jacobian;
 }
 
-mfem::Vector PolytropeOperator::GetJ00Diag() const {
-}
-
-mfem::Vector PolytropeOperator::GetJ01Diag() const {
-     if (!m_isFinalized) {
-         MFEM_ABORT("PolytropeOperator::Get01Diag -> GetJ01Diag called before finalization of PolytropeOperator");
-     }
-     if (m_Mmat == nullptr) {
-         MFEM_ABORT("PolytropeOperator::Get01Diag -> M sparse matrix has not been initialized before GetJ01Diag() call.");
-     }
-
-     mfem::Vector J01diag;
-     m_Mmat->GetDiag(J01diag);
-     J01diag *= -1;
-     return J01diag;
-}
-
-mfem::Vector PolytropeOperator::GetJ10diag() const {
-     if (!m_isFinalized) {
-         MFEM_ABORT("PolytropeOperator::Get10Diag -> GetJ10Diag called before finalization of PolytropeOperator");
-     }
-     if (m_Qmat == nullptr) {
-         MFEM_ABORT("PolytropeOperator::Get10Diag -> Q sparse matrix has not been initialized before GetJ10Diag() call.");
-     }
-     mfem::Vector J10diag;
-     m_Qmat->GetDiag(J10diag);
-     J10diag *= -1;
-     return J10diag;
-}
-
-mfem::Vector PolytropeOperator::GetJ11diag() const {
-     if (!m_isFinalized) {
-         MFEM_ABORT("PolytropeOperator::Get11Diag -> GetJ11Diag called before finalization of PolytropeOperator");
-     }
-     if (m_Dmat == nullptr) {
-         MFEM_ABORT("PolytropeOperator::Get11Diag -> D sparse matrix has not been initialized before GetJ11Diag() call.");
-     }
-     mfem::Vector J11diag;
-     m_Dmat->GetDiag(J11diag);
-     J11diag *= -1;
-     return J11diag;
-}
-
-
 void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     if (!m_isFinalized) {
         MFEM_ABORT("PolytropeOperator::Mult called before finalize");
@@ -209,7 +165,6 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     subtract(Dphi_term, Qtheta_term, y_R1);
 
 
-
     // -- Apply essential boundary conditions --
     for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
         int idx = m_theta_ess_tdofs.first[i];
@@ -232,16 +187,35 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
 
 }
 
+// TODO: I was *very* stupid and I accidently deleted a lot of the code
+//       which I had written to find and use the preconditioner. This needs
+//       to be reimplemented. Once that is working you can get back to
+//       Trying to understand the multimodal hump in the residuals vector.
+//       There is a jupyter notebook about this. I was thinking that this was
+//       perhapse related to the non consistent application of boundary conditions.
+void approxJacobiInvert(const mfem::SparseMatrix& mat, std::unique_ptr<mfem::SparseMatrix>& invMat, std::string name="matrix") {
+    // TODO: This likely can be made much more efficient and will probably be called in tight loops, a good
+    // place for some easy optimization might be here.
+    // Get the diagonal of the matrix
+    mfem::Vector diag;
+    mat.GetDiag(diag);
+
+    // Invert the diagonal
+    for (int i = 0; i < diag.Size(); i++) {
+        if (diag(i) == 0.0) {
+            MFEM_ABORT("Diagonal element (" + std::to_string(i) +") in " + name + " is zero, cannot invert.");
+        }
+        diag(i) = 1.0 / diag(i);
+    }
+
+    invMat = std::make_unique<mfem::SparseMatrix>(diag);
+}
+
 void PolytropeOperator::updateInverseNonlinearJacobian(const mfem::Operator &grad) const {
     if (const auto *sparse_mat = dynamic_cast<const mfem::SparseMatrix*>(&grad); sparse_mat != nullptr) {
-        mfem::Vector gradDiag;
-        sparse_mat->GetDiag(gradDiag);
-        for (int i = 0; i < gradDiag.Size(); i++) {
-            gradDiag(i) = 1.0/gradDiag(i); // Invert the diagonals of the jacobian
-        }
-        m_invNonlinearJacobian = std::make_unique<mfem::SparseMatrix>(gradDiag);
+        approxJacobiInvert(*sparse_mat, m_invNonlinearJacobian, "Nonlinear Jacobian");
     } else {
-        MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear Jacobian");
+        MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear jacobian where nonlinear jacobian is not trivially castable to a sparse matrix");
     }
 }
 
@@ -258,9 +232,6 @@ void PolytropeOperator::updateInverseSchurCompliment() const {
 
     schurCompliment->Add(-1, *temp_QxdfInvxM); // D - Q * df^{-1} * M
 
-    writeOperatorToCSV(*schurCompliment, "/Users/tboudreaux/Desktop/SchursCompliment.csv");
-    std::cout << "Wrote" << std::endl;
-
     // Note: EMB (April 9, 2025) Where I left for the day, so I don't pull my hair out tomorrow:
     /*
      * Need to calculate the inverse of Schur's compliment to precondition the jacobian