refactor(serif): refactored entire codebase into serif and sub namespaces

src/polytrope/utils/meson.build

@@ -0,0 +1,50 @@
+# ***********************************************************************
+#
+#   Copyright (C) 2025 -- The 4D-STAR Collaboration
+#   File Author: Emily Boudreaux
+#   Last Modified: March 19, 2025
+#
+#   4DSSE is free software; you can use it and/or modify
+#   it under the terms and restrictions the GNU General Library Public
+#   License version 3 (GPLv3) as published by the Free Software Foundation.
+#
+#   4DSSE is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+#   See the GNU Library General Public License for more details.
+#
+#   You should have received a copy of the GNU Library General Public License
+#   along with this software; if not, write to the Free Software
+#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# *********************************************************************** #
+polyutils_sources = files(
+  'private/integrators.cpp',
+  'private/polytropeOperator.cpp',
+  'private/utilities.cpp',
+)
+
+dependencies = [
+  mfem_dep,
+  macros_dep,
+  probe_dep,
+  quill_dep,
+  config_dep,
+  types_dep,
+  mfemanalysis_dep,
+]
+
+libpolyutils = static_library('polyutils',
+  polyutils_sources,
+  include_directories : include_directories('./public'),
+  cpp_args: ['-fvisibility=default'],
+  dependencies: dependencies,
+  install: true
+)
+
+polyutils_dep = declare_dependency(
+  include_directories : include_directories('./public'),
+  link_with : libpolyutils,
+  sources : polyutils_sources,
+  dependencies : dependencies
+)

src/polytrope/utils/private/integrators.cpp

@@ -0,0 +1,146 @@
+/* ***********************************************************************
+//
+//   Copyright (C) 2025 -- The 4D-STAR Collaboration
+//   File Author: Emily Boudreaux
+//   Last Modified: April 21, 2025
+//
+//   4DSSE is free software; you can use it and/or modify
+//   it under the terms and restrictions the GNU General Library Public
+//   License version 3 (GPLv3) as published by the Free Software Foundation.
+//
+//   4DSSE is distributed in the hope that it will be useful,
+//   but WITHOUT ANY WARRANTY; without even the implied warranty of
+//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+//   See the GNU Library General Public License for more details.
+//
+//   You should have received a copy of the GNU Library General Public License
+//   along with this software; if not, write to the Free Software
+//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// *********************************************************************** */
+#include "mfem.hpp"
+#include <cmath>
+
+#include "integrators.h"
+#include "config.h"
+#include <string>
+
+namespace serif {
+namespace polytrope {
+namespace polyMFEMUtils {
+  NonlinearPowerIntegrator::NonlinearPowerIntegrator(const double n) :
+    m_polytropicIndex(n),
+    m_epsilon(serif::config::Config::getInstance().get<double>("Poly:Solver:Epsilon", 1.0e-8)) {
+
+    if (m_polytropicIndex < 0.0) {
+      throw std::invalid_argument("Polytropic index must be non-negative.");
+    }
+    if (m_polytropicIndex > 5.0) {
+      throw std::invalid_argument("Polytropic index must be less than 5.0.");
+    }
+    if (m_epsilon <= 0.0) {
+      throw std::invalid_argument("Epsilon (Poly:Solver:Epsilon) must be positive non-zero.");
+    }
+  }
+
+  void NonlinearPowerIntegrator::AssembleElementVector(
+    const mfem::FiniteElement &el,
+    mfem::ElementTransformation &Trans,
+    const mfem::Vector &elfun,
+    mfem::Vector &elvect) {
+      
+      const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
+      int dof = el.GetDof();
+      elvect.SetSize(dof);
+      elvect = 0.0;
+
+      mfem::Vector shape(dof);
+      mfem::Vector physCoord;
+      for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
+        mfem::IntegrationPoint ip = ir->IntPoint(iqp);
+        Trans.SetIntPoint(&ip);
+        const double weight = ip.weight * Trans.Weight();
+  
+        el.CalcShape(ip, shape);
+  
+        double u_val = 0.0;
+        for (int j = 0; j < dof; j++) {
+          u_val += elfun(j) * shape(j);
+        }
+
+        double u_nl;
+        Trans.Transform(ip, physCoord);
+        const double r = physCoord.Norml2();
+        std::ofstream outFile("r.dat", std::ios::app);
+        outFile << r << '\n';
+        outFile.close();
+        if (r > m_regularizationRadius) {
+          if (u_val < m_epsilon) {
+            u_nl = fmod(u_val, m_polytropicIndex, m_epsilon);
+          } else {
+            u_nl = std::pow(u_val, m_polytropicIndex);
+          }
+        } else {
+          u_nl = 1.0 - m_polytropicIndex * m_regularizationCoeff * std::pow(r, 2);
+        }
+  
+        for (int i = 0; i < dof; i++){
+          elvect(i) += shape(i) * u_nl * weight;
+        }
+      }
+    }
+  
+
+    void NonlinearPowerIntegrator::AssembleElementGrad (
+      const mfem::FiniteElement &el,
+      mfem::ElementTransformation &Trans,
+      const mfem::Vector &elfun,
+      mfem::DenseMatrix &elmat) {
+    
+        const mfem::IntegrationRule *ir = &mfem::IntRules.Get(el.GetGeomType(), 2 * el.GetOrder() + 3);
+        const int dof = el.GetDof();
+        elmat.SetSize(dof);
+        elmat = 0.0;
+        mfem::Vector shape(dof);
+        mfem::DenseMatrix dshape(dof, 3);
+        mfem::DenseMatrix invJ(3, 3);
+        mfem::Vector physCoord;
+
+        for (int iqp = 0; iqp < ir->GetNPoints(); iqp++) {
+          const mfem::IntegrationPoint &ip = ir->IntPoint(iqp);
+          Trans.SetIntPoint(&ip);
+          const double weight = ip.weight * Trans.Weight();
+          Trans.Transform(ip, physCoord);
+          double r = physCoord.Norml2();
+    
+          el.CalcShape(ip, shape);
+
+          double u_val = 0.0;
+    
+          for (int j = 0; j < dof; j++) {
+            u_val += elfun(j) * shape(j);
+          }
+
+          double d_u_nl;
+          if (r > m_regularizationRadius) {
+            if (u_val < m_epsilon) {
+              d_u_nl = dfmod(m_epsilon, m_polytropicIndex);
+            } else {
+              d_u_nl = m_polytropicIndex * std::pow(u_val, m_polytropicIndex - 1.0);
+            }
+          } else {
+            d_u_nl = 0.0;
+          }
+
+          for (int i = 0; i < dof; i++) {
+            for (int j = 0; j < dof; j++) {
+              elmat(i, j) += shape(i) * d_u_nl * shape(j) * weight;
+            }
+          }
+
+    }
+  }
+
+} // namespace polyMFEMUtils
+} // namespace polytrope
+} // namespace serif

src/polytrope/utils/private/polytropeOperator.cpp

@@ -0,0 +1,493 @@
+/* ***********************************************************************
+//
+//   Copyright (C) 2025 -- The 4D-STAR Collaboration
+//   File Author: Emily Boudreaux
+//   Last Modified: April 21, 2025
+//
+//   4DSSE is free software; you can use it and/or modify
+//   it under the terms and restrictions the GNU General Library Public
+//   License version 3 (GPLv3) as published by the Free Software Foundation.
+//
+//   4DSSE is distributed in the hope that it will be useful,
+//   but WITHOUT ANY WARRANTY; without even the implied warranty of
+//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+//   See the GNU Library General Public License for more details.
+//
+//   You should have received a copy of the GNU Library General Public License
+//   along with this software; if not, write to the Free Software
+//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// *********************************************************************** */
+#include "polytropeOperator.h"
+#include "4DSTARTypes.h"
+#include "utilities.h"
+
+#include "mfem.hpp"
+#include "mfem_smout.h"
+#include <memory>
+
+#include "config.h"
+
+namespace serif {
+namespace polytrope {
+    // --- SchurCompliment Class Implementation ---
+
+    // SchurCompliment: Constructors
+SchurCompliment::SchurCompliment(
+    const mfem::SparseMatrix &QOp,
+    const mfem::SparseMatrix &DOp,
+    const mfem::SparseMatrix &MOp,
+    const mfem::Solver &GradInvOp ) :
+mfem::Operator(DOp.Height(), DOp.Width()) {
+    // Initialize sizes
+    SetOperator(QOp, DOp, MOp, GradInvOp);
+    m_nPhi = m_DOp->Height();
+    m_nTheta = m_MOp->Height();
+}
+
+SchurCompliment::SchurCompliment(
+    const mfem::SparseMatrix &QOp,
+    const mfem::SparseMatrix &DOp,
+    const mfem::SparseMatrix &MOp) :
+mfem::Operator(DOp.Height(), DOp.Width())
+{
+    updateConstantTerms(QOp, DOp, MOp);
+    m_nPhi = m_DOp->Height();
+    m_nTheta = m_MOp->Height();
+}
+
+// SchurCompliment: Public Interface Methods
+void SchurCompliment::Mult(const mfem::Vector &x, mfem::Vector &y) const {
+    // Check that the input vector is the correct size
+    if (x.Size() != m_nPhi) {
+        MFEM_ABORT("Input vector x has size " + std::to_string(x.Size()) + ", expected " + std::to_string(m_nPhi));
+    }
+    if (y.Size() != m_nPhi) {
+        MFEM_ABORT("Output vector y has size " + std::to_string(y.Size()) + ", expected " + std::to_string(m_nPhi));
+    }
+
+    // Check that the operators are set
+    if (m_QOp == nullptr) {
+        MFEM_ABORT("QOp is null in SchurCompliment::Mult");
+    }
+    if (m_DOp == nullptr) {
+        MFEM_ABORT("DOp is null in SchurCompliment::Mult");
+    }
+    if (m_MOp == nullptr) {
+        MFEM_ABORT("MOp is null in SchurCompliment::Mult");
+    }
+    if (m_GradInvOp == nullptr) {
+        MFEM_ABORT("GradInvOp is null in SchurCompliment::Mult");
+    }
+
+    mfem::Vector v1(m_nTheta); // M * x
+    m_MOp -> Mult(x, v1); // M * x
+
+    mfem::Vector v2(m_nTheta); // GradInv * M * x
+    m_GradInvOp -> Mult(v1, v2); // GradInv * M * x
+
+    mfem::Vector v3(m_nPhi); // Q * GradInv * M * x
+    m_QOp -> Mult(v2, v3); // Q * GradInv * M * x
+
+    mfem::Vector v4(m_nPhi); // D * x
+    m_DOp -> Mult(x, v4); // D * x
+
+    subtract(v4, v3, y); // (D - Q * GradInv * M) * x
+}
+
+void SchurCompliment::SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp) {
+    updateConstantTerms(QOp, DOp, MOp);
+    updateInverseNonlinearJacobian(GradInvOp);
+}
+
+void SchurCompliment::updateInverseNonlinearJacobian(const mfem::Solver &gradInv) {
+    m_GradInvOp = &gradInv;
+}
+
+// SchurCompliment: Private Helper Methods
+void SchurCompliment::updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp) {
+    m_QOp = &QOp;
+    m_DOp = &DOp;
+    m_MOp = &MOp;
+}
+
+
+// --- GMRESInverter Class Implementation ---
+
+// GMRESInverter: Constructor
+GMRESInverter::GMRESInverter(const SchurCompliment &op) :
+mfem::Operator(op.Height(), op.Width()),
+m_op(op) {
+    m_solver.SetOperator(m_op);
+    m_solver.SetMaxIter(100);
+    m_solver.SetRelTol(1e-1);
+    m_solver.SetAbsTol(1e-1);
+}
+
+// GMRESInverter: Public Interface Methods
+void GMRESInverter::Mult(const mfem::Vector &x, mfem::Vector &y) const {
+    m_solver.Mult(x, y); // Approximates m_op^-1 * x
+}
+
+
+// --- PolytropeOperator Class Implementation ---
+    // PolytropeOperator: Constructor
+    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) {
+
+        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
+        // t = 0 (symmetric Gauss-Seidel), 1 (forward Gauss-Seidel), 2 (backward Gauss-Seidel)
+        // iterations = 3
+        m_invNonlinearJacobian = std::make_unique<mfem::GSSmoother>(0, 3);
+    }
+
+// PolytropeOperator: Core Operator Overrides
+void PolytropeOperator::Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const {
+    if (!m_isFinalized) {
+        MFEM_ABORT("PolytropeOperator::Mult called before finalize");
+    }
+
+    // TODO: confirm that the vectors xFree and m_freeDofs are always parallel
+    m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector
+    mfem::Vector y;
+    y.SetSize(m_blockOffsets.Last());
+
+    // -- Create BlockVector views for input x and output y
+    mfem::BlockVector x_block(const_cast<mfem::Vector&>(m_state), m_blockOffsets);
+    mfem::BlockVector y_block(y, m_blockOffsets);
+
+    // -- Get Vector views for individual blocks
+    const mfem::Vector &x_theta = x_block.GetBlock(0);
+    const mfem::Vector &x_phi = x_block.GetBlock(1);
+    mfem::Vector &y_R0 = y_block.GetBlock(0); // Residual Block 0 (theta)
+    mfem::Vector &y_R1 = y_block.GetBlock(1); // Residual Block 1 (phi)
+
+    int theta_size = m_blockOffsets[1] - m_blockOffsets[0];
+    int phi_size = m_blockOffsets[2] - m_blockOffsets[1];
+
+    mfem::Vector f_term(theta_size);
+    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(θ) + (1+c)Mɸ + cSθ
+    // 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);      // f(θ)
+    m_M->Mult(x_phi, Mphi_term);     // Mɸ
+    m_D->Mult(x_phi, Dphi_term);     // Dɸ
+    m_Q->Mult(x_theta, Qtheta_term); // Qθ
+    m_S->Mult(x_theta, Stheta_term); // Sθ
+
+    // PERF: these multiplications may be able to be refactored into the matrices so they only need to be done once.
+    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) {
+        if (m_freeDofs.Find(i) != -1) {
+            yFree[j] = y[i];
+            j++;
+        }
+    }
+}
+
+mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &xFree) const {
+    if (!m_isFinalized) {
+        MFEM_ABORT("PolytropeOperator::GetGradient called before finalize");
+    }
+    m_state.SetSubVector(m_freeDofs, xFree); // Scatter the free dofs from the input vector xFree into the state vector
+    // --- Get the gradient of f ---
+    mfem::BlockVector x_block(const_cast<mfem::Vector&>(m_state), m_blockOffsets);
+    const mfem::Vector& x_theta = x_block.GetBlock(0);
+
+    // 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); // ∂f(θ)/∂θ
+
+    if (gradMatrix == nullptr) {
+        MFEM_ABORT("PolytropeOperator::GetGradient: Gradient is not a SparseMatrix.");
+    }
+
+    m_gradReduced = std::make_unique<mfem::SparseMatrix> (
+        serif::utilities::build_reduced_matrix(
+            *gradMatrix,
+            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());
+
+    return *m_jacobian;
+}
+
+
+// PolytropeOperator: Setup and Finalization
+void PolytropeOperator::finalize(const mfem::Vector &initTheta) {
+    using serif::utilities::build_reduced_matrix;
+
+    if (m_isFinalized) {
+        return;
+    }
+
+    // These functions must be called in this order since they depend on each others post state
+    // TODO: Refactor this so that either there are explicit checks to make sure the order is correct or make
+    //       them pure functions
+    construct_matrix_representations();
+    construct_reduced_block_offsets();
+    construct_jacobian_constant_terms();
+    scatter_boundary_conditions();
+    populate_free_dof_array();
+
+    // Override the size based on the reduced system
+    height = m_reducedBlockOffsets.Last();
+    width = m_reducedBlockOffsets.Last();
+
+    m_isFinalized = true;
+
+}
+
+// PolytropeOperator: Essential True DOF Management
+void PolytropeOperator::set_essential_true_dofs(const serif::types::MFEMArrayPair& theta_ess_tdofs, const serif::types::MFEMArrayPair& phi_ess_tdofs) {
+    m_isFinalized = false;
+    m_theta_ess_tdofs = theta_ess_tdofs;
+    m_phi_ess_tdofs = phi_ess_tdofs;
+
+    if (m_f) {
+        m_f->SetEssentialTrueDofs(theta_ess_tdofs.first);
+    } else {
+        MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs");
+    }
+}
+
+void PolytropeOperator::set_essential_true_dofs(const serif::types::MFEMArrayPairSet& ess_tdof_pair_set) {
+    set_essential_true_dofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second);
+}
+
+serif::types::MFEMArrayPairSet PolytropeOperator::get_essential_true_dofs() const {
+     return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs);
+}
+
+// PolytropeOperator: Getter Methods
+const mfem::BlockOperator &PolytropeOperator::get_jacobian_operator() const {
+    if (m_jacobian == nullptr) {
+        MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call.");
+    }
+    if (!m_isFinalized) {
+        MFEM_ABORT("PolytropeOperator not finalized prior to call to GetJacobianOperator().");
+    }
+
+    return *m_jacobian;
+}
+
+mfem::BlockDiagonalPreconditioner& PolytropeOperator::get_preconditioner() const {
+    if (m_schurPreconditioner == nullptr) {
+        MFEM_ABORT("Schur preconditioner has not been initialized before GetPreconditioner() call.");
+    }
+    if (!m_isFinalized) {
+        MFEM_ABORT("PolytropeOperator not finalized prior to call to GetPreconditioner().");
+    }
+    return *m_schurPreconditioner;
+}
+
+int PolytropeOperator::get_reduced_system_size() const {
+    if (!m_isFinalized) {
+        MFEM_ABORT("PolytropeOperator not finalized prior to call to GetReducedSystemSize().");
+    }
+    return m_reducedBlockOffsets.Last();
+}
+
+// PolytropeOperator: State Reconstruction
+const mfem::Vector &PolytropeOperator::reconstruct_full_state_vector(const mfem::Vector &reducedState) const {
+    m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector
+    return m_state;
+}
+
+const mfem::BlockVector PolytropeOperator::reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const {
+    m_state.SetSubVector(m_freeDofs, reducedState); // Scatter the reduced state vector into the full state vector
+    mfem::BlockVector x_block(m_state, m_blockOffsets);
+    return x_block;
+}
+
+// PolytropeOperator: DOF Population
+void PolytropeOperator::populate_free_dof_array() {
+    m_freeDofs.SetSize(0);
+    for (int i = 0; i < m_blockOffsets.Last(); i++) {
+        const int thetaSearchIndex = i;
+        const int phiSearchIndex = i - m_blockOffsets[1];
+        if (phiSearchIndex < 0){
+            if (m_theta_ess_tdofs.first.Find(thetaSearchIndex) == -1) {
+                m_freeDofs.Append(i);
+            }
+        } else {
+            if (m_phi_ess_tdofs.first.Find(phiSearchIndex) == -1) {
+                m_freeDofs.Append(i);
+            }
+        }
+    }
+}
+
+// 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());
+
+    // --- Reduce the matrices to only the free DOFs ---
+    m_MReduced = std::make_unique<mfem::SparseMatrix>(
+        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)
+    );
+
+    // --- 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 + theta)
+}
+
+void PolytropeOperator::construct_jacobian_constant_terms() {
+    m_jacobian = std::make_unique<mfem::BlockOperator>(m_reducedBlockOffsets);
+
+    m_jacobian->SetBlock(0, 1, m_MScaledReduced.get());     //<-  (1+c)M (constant, reduced to only free DOFs)
+    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() {
+    mfem::Vector thetaStateValues(m_theta_ess_tdofs.first.Size());
+    for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
+        thetaStateValues[i] = m_theta_ess_tdofs.second[i];
+    }
+    mfem::Vector phiStateValues(m_phi_ess_tdofs.first.Size());
+    for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
+        phiStateValues[i] = m_phi_ess_tdofs.second[i]; // TODO: figure out if this needs to be normalized
+    }
+
+    mfem::Array<int> phiDofIndices(m_phi_ess_tdofs.first.Size());
+    for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
+        phiDofIndices[i] = m_phi_ess_tdofs.first[i] + m_blockOffsets[1];
+    }
+
+    m_state.SetSize(m_blockOffsets.Last());
+    m_state = 0.0;
+    m_state.SetSubVector(m_theta_ess_tdofs.first, thetaStateValues);
+    m_state.SetSubVector(phiDofIndices, phiStateValues);
+
+}
+
+// PolytropeOperator: Private Helper Methods - Jacobian and Preconditioner Updates
+void PolytropeOperator::update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const {
+    m_invNonlinearJacobian->SetOperator(grad);
+}
+
+void PolytropeOperator::update_inverse_schur_compliment() const {
+    // TODO: This entire function could probably be refactored out
+    if (!m_isFinalized) {
+        MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before finalize");
+    }
+    if (m_invNonlinearJacobian == nullptr) {
+        MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian");
+    }
+    if (m_schurCompliment == nullptr) {
+        MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseSchurCompliment");
+    }
+
+    m_schurCompliment->updateInverseNonlinearJacobian(*m_invNonlinearJacobian);
+
+    if (m_schurPreconditioner == nullptr) {
+        m_schurPreconditioner = std::make_unique<mfem::BlockDiagonalPreconditioner>(m_blockOffsets);
+    }
+
+    m_schurPreconditioner->SetDiagonalBlock(0, m_invNonlinearJacobian.get());
+    m_schurPreconditioner->SetDiagonalBlock(1, m_invSchurCompliment.get());
+
+}
+
+void PolytropeOperator::update_preconditioner(const mfem::Operator &grad) const {
+    update_inverse_nonlinear_jacobian(grad);
+    update_inverse_schur_compliment();
+}
+
+} // namespace polytrope
+} // namespace serif

src/polytrope/utils/private/utilities.cpp

@@ -0,0 +1,132 @@
+#include "utilities.h"
+#include "mfem.hpp"
+#include <memory>
+
+namespace serif::utilities {
+    mfem::SparseMatrix build_reduced_matrix(
+        const mfem::SparseMatrix& matrix,
+        const mfem::Array<int>& trialEssentialDofs,
+        const mfem::Array<int>& testEssentialDofs
+    ) {
+        int M_orig = matrix.Height();
+        int N_orig = matrix.Width();
+
+        // 1. Create boolean lookup tables for eliminated rows/columns
+        // These tables help quickly check if an original row/column index is eliminated.
+        mfem::Array<bool> row_is_eliminated(M_orig);
+        row_is_eliminated = false; // Initialize all to false (no rows eliminated yet)
+        for (int i = 0; i < testEssentialDofs.Size(); ++i) {
+            int r_idx = testEssentialDofs[i];
+            if (r_idx >= 0 && r_idx < M_orig) { // Check for valid index
+                row_is_eliminated[r_idx] = true;
+            }
+        }
+
+        mfem::Array<bool> col_is_eliminated(N_orig);
+        col_is_eliminated = false; // Initialize all to false (no columns eliminated yet)
+        for (int i = 0; i < trialEssentialDofs.Size(); ++i) {
+            int c_idx = trialEssentialDofs[i];
+            if (c_idx >= 0 && c_idx < N_orig) { // Check for valid index
+                col_is_eliminated[c_idx] = true;
+            }
+        }
+
+        // 2. Create mappings from old (original) indices to new indices
+        // Also, count the number of rows and columns in the new, reduced matrix.
+        mfem::Array<int> old_row_to_new_row(M_orig);
+        int num_new_rows = 0;
+        for (int i = 0; i < M_orig; ++i) {
+            if (row_is_eliminated[i]) {
+                old_row_to_new_row[i] = -1; // Mark as eliminated (no corresponding new row)
+            } else {
+                old_row_to_new_row[i] = num_new_rows++; // Assign the next available new row index
+            }
+        }
+
+        mfem::Array<int> old_col_to_new_col(N_orig);
+        int num_new_cols = 0;
+        for (int i = 0; i < N_orig; ++i) {
+            if (col_is_eliminated[i]) {
+                old_col_to_new_col[i] = -1; // Mark as eliminated (no corresponding new column)
+            } else {
+                old_col_to_new_col[i] = num_new_cols++; // Assign the next available new column index
+            }
+        }
+
+        // 3. Create the new sparse matrix with the calculated reduced dimensions.
+        // It's initially empty (no non-zero entries).
+        mfem::SparseMatrix A_new(num_new_rows, num_new_cols);
+
+        // 4. Iterate through the non-zero entries of the original matrix (A_orig).
+        // A_orig is typically stored in Compressed Sparse Row (CSR) format.
+        // GetI() returns row pointers, GetJ() returns column indices, GetData() returns values.
+        const int* I_orig = matrix.GetI();
+        const int* J_orig = matrix.GetJ();
+        const double* V_orig = matrix.GetData(); // Assuming double, can be templated if needed
+
+        for (int r_orig = 0; r_orig < M_orig; ++r_orig) {
+            // If the original row is marked for elimination, skip all its entries.
+            if (row_is_eliminated[r_orig]) {
+                continue;
+            }
+
+            // Get the new row index for the current original row.
+            int r_new = old_row_to_new_row[r_orig];
+
+            // Iterate through non-zero entries in the current original row r_orig.
+            // I_orig[r_orig] is the start index in J_orig and V_orig for row r_orig.
+            // I_orig[r_orig+1]-1 is the end index.
+            for (int k = I_orig[r_orig]; k < I_orig[r_orig + 1]; ++k) {
+                int c_orig = J_orig[k]; // Original column index of the non-zero entry
+                double val = V_orig[k]; // Value of the non-zero entry
+
+                if (col_is_eliminated[c_orig]) {
+                    continue;
+                }
+
+                int c_new = old_col_to_new_col[c_orig];
+
+                A_new.Add(r_new, c_new, val);
+            }
+        }
+
+        A_new.Finalize();
+
+        return A_new;
+    }
+
+    mfem::Vector build_dof_identification_vector(const mfem::Array<int>& allDofs, const::mfem::Array<int>& highlightDofs) {
+        mfem::Vector v(allDofs.Size());
+        v = 0.0; // Initialize the vector to zero
+        v.SetSubVector(highlightDofs, 1.0); // Set the highlighted dofs to 1.0
+        return v;
+    }
+
+    mfem::GridFunction compute_curl(mfem::GridFunction& phi_gf) {
+        mfem::Mesh* mesh = phi_gf.FESpace()->GetMesh();
+        const int dim = mesh->Dimension();
+        // Match the polynomial order of the original RT space for consistency.
+        const int order = phi_gf.FESpace()->GetOrder(0);
+        mfem::Vector curl_mag_vec;
+
+        for (int ne = 0; ne < mesh->GetNE(); ++ne) {
+            if (mesh->GetElementType(ne) != mfem::Element::TRIANGLE &&
+                mesh->GetElementType(ne) != mfem::Element::QUADRILATERAL &&
+                mesh->GetElementType(ne) != mfem::Element::TETRAHEDRON &&
+                mesh->GetElementType(ne) != mfem::Element::HEXAHEDRON) {
+                throw std::invalid_argument("Mesh element type not supported for curl computation.");
+            }
+            mfem::IsoparametricTransformation T;
+            mesh->GetElementTransformation(ne, &T);
+            phi_gf.GetCurl(T, curl_mag_vec);
+            std::cout << "HERE" << std::endl;
+        }
+        mfem::L2_FECollection fac(order, dim);
+        mfem::FiniteElementSpace fs(mesh, &fac);
+        mfem::GridFunction curl_gf(&fs);
+        curl_gf = 0.0;
+        return curl_gf;
+
+    }
+
+}

src/polytrope/utils/public/integrators.h

@@ -0,0 +1,107 @@
+/* ***********************************************************************
+//
+//   Copyright (C) 2025 -- The 4D-STAR Collaboration
+//   File Author: Emily Boudreaux
+//   Last Modified: April 21, 2025
+//
+//   4DSSE is free software; you can use it and/or modify
+//   it under the terms and restrictions the GNU General Library Public
+//   License version 3 (GPLv3) as published by the Free Software Foundation.
+//
+//   4DSSE is distributed in the hope that it will be useful,
+//   but WITHOUT ANY WARRANTY; without even the implied warranty of
+//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+//   See the GNU Library General Public License for more details.
+//
+//   You should have received a copy of the GNU Library General Public License
+//   along with this software; if not, write to the Free Software
+//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// *********************************************************************** */
+#pragma once
+
+#include "mfem.hpp"
+#include <string>
+#include "config.h"
+#include "probe.h"
+
+
+/**
+ * @file integrators.h
+ * @brief A collection of utilities for working with MFEM and solving the lane-emden equation.
+ */
+
+namespace serif {
+namespace polytrope {
+/**
+ * @namespace polyMFEMUtils
+ * @brief A namespace for utilities for working with MFEM and solving the lane-emden equation.
+ */
+namespace polyMFEMUtils {
+  /**
+   * @brief A class for nonlinear power integrator.
+   */
+  class NonlinearPowerIntegrator: public mfem::NonlinearFormIntegrator {
+  public:
+    /**
+     * @brief Constructor for NonlinearPowerIntegrator.
+     *
+     * @param coeff The function coefficient.
+     * @param n The polytropic index.
+     */
+    NonlinearPowerIntegrator(double n);
+
+    /**
+     * @brief Assembles the element vector.
+     *
+     * @param el The finite element.
+     * @param Trans The element transformation.
+     * @param elfun The element function.
+     * @param elvect The element vector to be assembled.
+     */
+     virtual void AssembleElementVector(const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::Vector &elvect) override;
+    /**
+     * @brief Assembles the element gradient.
+     *
+     * @param el The finite element.
+     * @param Trans The element transformation.
+     * @param elfun The element function.
+     * @param elmat The element matrix to be assembled.
+     */
+     virtual void AssembleElementGrad (const mfem::FiniteElement &el, mfem::ElementTransformation &Trans, const mfem::Vector &elfun, mfem::DenseMatrix &elmat) override;
+  private:
+    serif::config::Config& m_config = serif::config::Config::getInstance();
+    serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance();
+    quill::Logger* m_logger = m_logManager.getLogger("log");
+    double m_polytropicIndex;
+    double m_epsilon;
+    static constexpr double m_regularizationRadius = 0.15; ///< Regularization radius for the epsilon function, used to avoid singularities in the power law.
+    static constexpr double m_regularizationCoeff = 1.0/6.0; ///< Coefficient for the regularization term, used to ensure smoothness in the power law.
+  };
+
+  inline double dfmod(const double epsilon, const double n) {
+    if (n == 0.0) {
+      return 0.0;
+    }
+    if (n == 1.0) {
+      return 1.0;
+    }
+    return n * std::pow(epsilon, n - 1.0);
+  }
+
+  inline double fmod(const double theta, const double n, const double epsilon) {
+    if (n == 0.0) {
+      return 1.0;
+    }
+    // For n != 0
+    const double y_prime_at_epsilon = dfmod(epsilon, n); // Uses the robust dfmod
+    const double y_at_epsilon = std::pow(epsilon, n);     // epsilon^n
+
+    // f_mod(theta) = y_at_epsilon + y_prime_at_epsilon * (theta - epsilon)
+    return y_at_epsilon + y_prime_at_epsilon * (theta - epsilon);
+  }
+
+
+} // namespace polyMFEMUtils
+} // namespace polytrope
+} // namespace serif

src/polytrope/utils/public/polytropeOperator.h

@@ -0,0 +1,399 @@
+/* ***********************************************************************
+//
+//   Copyright (C) 2025 -- The 4D-STAR Collaboration
+//   File Author: Emily Boudreaux
+//   Last Modified: April 21, 2025
+//
+//   4DSSE is free software; you can use it and/or modify
+//   it under the terms and restrictions the GNU General Library Public
+//   License version 3 (GPLv3) as published by the Free Software Foundation.
+//
+//   4DSSE is distributed in the hope that it will be useful,
+//   but WITHOUT ANY WARRANTY; without even the implied warranty of
+//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+//   See the GNU Library General Public License for more details.
+//
+//   You should have received a copy of the GNU Library General Public License
+//   along with this software; if not, write to the Free Software
+//   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// *********************************************************************** */
+#pragma once
+
+#include "mfem.hpp"
+#include "4DSTARTypes.h"
+#include <memory>
+
+#include "probe.h"
+
+namespace serif {
+namespace polytrope {
+
+/**
+ * @brief Represents the Schur complement operator used in the solution process.
+ *
+ * This class computes S = D - Q * GradInv * M.
+ */
+class SchurCompliment final : public mfem::Operator {
+public:
+    /**
+     * @brief Constructs a SchurCompliment operator.
+     * @param QOp The Q matrix operator.
+     * @param DOp The D matrix operator.
+     * @param MOp The M matrix operator.
+     * @param GradInvOp The inverse of the gradient operator.
+     */
+    SchurCompliment(
+        const mfem::SparseMatrix &QOp,
+        const mfem::SparseMatrix &DOp,
+        const mfem::SparseMatrix &MOp,
+        const mfem::Solver &GradInvOp
+    );
+
+    /**
+     * @brief Constructs a SchurCompliment operator without the inverse gradient initially.
+     * The inverse gradient must be set later using updateInverseNonlinearJacobian.
+     * @param QOp The Q matrix operator.
+     * @param DOp The D matrix operator.
+     * @param MOp The M matrix operator.
+     */
+    SchurCompliment(
+        const mfem::SparseMatrix &QOp,
+        const mfem::SparseMatrix &DOp,
+        const mfem::SparseMatrix &MOp
+    );
+
+    /**
+     * @brief Destructor.
+     */
+    ~SchurCompliment() override = default;
+
+    /**
+     * @brief Applies the Schur complement operator: y = (D - Q * GradInv * M) * x.
+     * @param x The input vector.
+     * @param y The output vector.
+     */
+    void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+    /**
+     * @brief Sets all operators for the Schur complement.
+     * @param QOp The Q matrix operator.
+     * @param DOp The D matrix operator.
+     * @param MOp The M matrix operator.
+     * @param GradInvOp The inverse of the gradient operator.
+     */
+    void SetOperator(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp, const mfem::Solver &GradInvOp);
+
+    /**
+     * @brief Updates the inverse of the nonlinear Jacobian (GradInv).
+     * @param gradInv The new inverse gradient solver.
+     */
+    void updateInverseNonlinearJacobian(const mfem::Solver &gradInv);
+
+private:
+    /**
+     * @brief Updates the constant matrix terms (Q, D, M).
+     * @param QOp The Q matrix operator.
+     * @param DOp The D matrix operator.
+     * @param MOp The M matrix operator.
+     */
+    void updateConstantTerms(const mfem::SparseMatrix &QOp, const mfem::SparseMatrix &DOp, const mfem::SparseMatrix &MOp);
+
+private:
+    // Pointers to external operators (not owned by this class)
+    const mfem::SparseMatrix* m_QOp = nullptr;         ///< Pointer to the Q matrix operator.
+    const mfem::SparseMatrix* m_DOp = nullptr;         ///< Pointer to the D matrix operator.
+    const mfem::SparseMatrix* m_MOp = nullptr;         ///< Pointer to the M matrix operator.
+    const mfem::Solver* m_GradInvOp = nullptr;         ///< Pointer to the inverse of the gradient operator.
+
+    // Dimensions
+    int m_nPhi = 0;                                    ///< Dimension related to the phi variable (typically number of rows/cols of D).
+    int m_nTheta = 0;                                  ///< Dimension related to the theta variable (typically number of rows of M).
+
+    // Owned resources
+    mutable std::unique_ptr<mfem::SparseMatrix> m_matrixForm; ///< Optional: if a matrix representation is ever explicitly formed.
+};
+
+/**
+ * @brief Provides an approximate inverse of the SchurCompliment operator using GMRES.
+ */
+class GMRESInverter final : public mfem::Operator {
+public:
+    /**
+     * @brief Constructs a GMRESInverter.
+     * @param op The SchurCompliment operator to invert.
+     */
+    explicit GMRESInverter(const SchurCompliment& op);
+
+    /**
+     * @brief Destructor.
+     */
+    ~GMRESInverter() override = default;
+
+    /**
+     * @brief Applies the GMRES solver to approximate op^-1 * x.
+     * @param x The input vector.
+     * @param y The output vector (approximation of op^-1 * x).
+     */
+    void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+private:
+    const SchurCompliment& m_op; ///< Reference to the SchurCompliment operator to be inverted.
+    mfem::GMRESSolver m_solver;   ///< GMRES solver instance.
+};
+
+/**
+ * @brief Represents the coupled nonlinear operator for the polytropic system.
+ *
+ * This operator defines the system F(X) = 0, where X = [θ, φ]^T,
+ * and F(X) = [ f(θ) + Mφ ]
+ *            [ Dφ  - Qθ ].
+ * It handles essential boundary conditions and the construction of the Jacobian.
+ */
+class PolytropeOperator final : public mfem::Operator {
+public:
+    /**
+     * @brief Constructs the PolytropeOperator.
+     * @param M The M bilinear form (coupling θ and φ).
+     * @param Q The Q bilinear form (coupling φ and θ).
+     * @param D The D bilinear form (acting on φ).
+     * @param f The nonlinear form f(θ).
+     * @param blockOffsets Offsets defining the blocks for θ and φ variables.
+     */
+    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);
+
+    /**
+     * @brief Destructor.
+     */
+    ~PolytropeOperator() override = default;
+
+    /**
+     * @brief Applies the PolytropeOperator: y = F(x).
+     * This computes the residual of the nonlinear system.
+     * @param xFree The vector of free (non-essential) degrees of freedom.
+     * @param yFree The resulting residual vector corresponding to free DOFs.
+     */
+    void Mult(const mfem::Vector &xFree, mfem::Vector &yFree) const override;
+
+    /**
+     * @brief Computes the Jacobian of the PolytropeOperator at a given state xFree.
+     * The Jacobian is of the form:
+     * J = [ G  M ]
+     *     [ -Q D ]
+     * where G is the gradient of f(θ).
+     * @param xFree The vector of free DOFs at which to evaluate the gradient.
+     * @return A reference to the Jacobian operator.
+     */
+    mfem::Operator& GetGradient(const mfem::Vector &xFree) const override;
+
+    /**
+     * @brief Finalizes the operator setup.
+     * This must be called after setting essential true DOFs and before using Mult or GetGradient.
+     * It constructs reduced matrices, block offsets, and populates free DOFs.
+     * @param initTheta Initial guess for θ, used to evaluate the nonlinear form if necessary during setup.
+     */
+    void finalize(const mfem::Vector &initTheta);
+
+    /**
+     * @brief Checks if the operator has been finalized.
+     * @return True if finalize() has been successfully called, false otherwise.
+     */
+    bool isFinalized() const { return m_isFinalized; }
+
+    /**
+     * @brief Sets the essential true degrees of freedom for both θ and φ variables.
+     * Marks the operator as not finalized.
+     * @param theta_ess_tdofs Pair of arrays: (indices of essential DOFs for θ, values at these DOFs).
+     * @param phi_ess_tdofs Pair of arrays: (indices of essential DOFs for φ, values at these DOFs).
+     */
+    void set_essential_true_dofs(const serif::types::MFEMArrayPair& theta_ess_tdofs, const serif::types::MFEMArrayPair& phi_ess_tdofs);
+
+    /**
+     * @brief Sets the essential true degrees of freedom using a pair of pairs.
+     * Marks the operator as not finalized.
+     * @param ess_tdof_pair_set A pair containing the essential TDOF pairs for theta and phi.
+     */
+    void set_essential_true_dofs(const serif::types::MFEMArrayPairSet& ess_tdof_pair_set);
+
+
+    /**
+     * @brief Reconstructs the full state vector (including essential DOFs) from a reduced state vector (free DOFs).
+     * @param reducedState The vector containing only the free degrees of freedom.
+     * @return Constant reference to the internal full state vector, updated with the reducedState.
+     */
+    [[nodiscard]] const mfem::Vector &reconstruct_full_state_vector(const mfem::Vector &reducedState) const;
+
+    /**
+     * @breif Reconstruct the full state vector (including essential DOFs) from a reduced state vector (free DOFs) as well as the block offsets.
+     * @param reducedState The vector containing only the free degrees of freedom.
+     * @return Constant reference to the internal full state vector, updated with the reducedState as a block vector.
+     */
+    [[nodiscard]] const mfem::BlockVector reconstruct_full_block_state_vector(const mfem::Vector &reducedState) const;
+
+    /**
+     * @brief Populates the internal array of free (non-essential) degree of freedom indices.
+     * This is called during finalize().
+     */
+    void populate_free_dof_array();
+
+    /// --- Getters for key internal state and operators ---
+    /**
+     * @brief Gets the Jacobian operator.
+     * Asserts that the operator is finalized and the Jacobian has been computed.
+     * @return Constant reference to the block Jacobian operator.
+     */
+    const mfem::BlockOperator &get_jacobian_operator() const;
+
+    /**
+     * @brief Gets the block diagonal preconditioner for the Schur complement system.
+     * Asserts that the operator is finalized and the preconditioner has been computed.
+     * @return Reference to the block diagonal preconditioner.
+     */
+    mfem::BlockDiagonalPreconditioner &get_preconditioner() const;
+
+
+    /// --- Getters for information on internal state ---
+    /**
+     * @brief Gets the full state vector, including essential DOFs.
+     * @return Constant reference to the internal state vector.
+     */
+    const mfem::Array<int>& get_free_dofs() const { return m_freeDofs; } ///< Getter for the free DOFs array.
+
+    /**
+     * @brief Gets the size of the reduced system (number of free DOFs).
+     * Asserts that the operator is finalized.
+     * @return The total number of free degrees of freedom.
+     */
+    int get_reduced_system_size() const;
+
+    /**
+     * @brief Gets the currently set essential true degrees of freedom.
+     * @return A pair containing the essential TDOF pairs for theta and phi.
+     */
+    serif::types::MFEMArrayPairSet get_essential_true_dofs() const;
+
+    /**
+     * @brief Gets the block offsets for the full (unreduced) system.
+     * @return Constant reference to the array of block offsets.
+     */
+    const mfem::Array<int>& get_block_offsets() const { return m_blockOffsets; }
+
+    /**
+     * @brief Gets the block offsets for the reduced system (after eliminating essential DOFs).
+     * @return Constant reference to the array of reduced block offsets.
+     */
+    const mfem::Array<int>& get_reduced_block_offsets() const {return m_reducedBlockOffsets; }
+
+private:
+    // --- Logging ---
+    serif::probe::LogManager& m_logManager = serif::probe::LogManager::getInstance(); ///< Reference to the global log manager.
+    quill::Logger* m_logger = m_logManager.getLogger("log");            ///< Pointer to the specific logger instance.
+
+    // --- Input Bilinear/Nonlinear Forms (owned) ---
+    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...) // 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 ---
+    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.
+
+    // --- Block Offsets ---
+    const mfem::Array<int> m_blockOffsets;        ///< Block offsets for the full system [0, size(θ), size(θ)+size(φ)].
+    mfem::Array<int> m_reducedBlockOffsets;       ///< Block offsets for the reduced system (free DOFs).
+
+    // --- Essential Boundary Conditions ---
+    serif::types::MFEMArrayPair m_theta_ess_tdofs;         ///< Essential true DOFs for θ (indices and values).
+    serif::types::MFEMArrayPair m_phi_ess_tdofs;           ///< Essential true DOFs for φ (indices and values).
+
+    // --- Jacobian and Preconditioner Components (owned, mutable) ---
+    std::unique_ptr<mfem::ScaledOperator> m_negQ_mat; ///< Scaled operator for -Q_reduced.
+    mutable std::unique_ptr<mfem::BlockOperator> m_jacobian; ///< Jacobian operator J = [G M; -Q D]_reduced.
+    mutable std::unique_ptr<SchurCompliment> m_schurCompliment; ///< Schur complement S = D_reduced - Q_reduced * G_inv_reduced * M_reduced.
+    mutable std::unique_ptr<GMRESInverter> m_invSchurCompliment; ///< Approximate inverse of the Schur complement.
+    mutable std::unique_ptr<mfem::Solver> m_invNonlinearJacobian; ///< Solver for the inverse of the G block (gradient of f(θ)_reduced).
+    mutable std::unique_ptr<mfem::BlockDiagonalPreconditioner> m_schurPreconditioner; ///< Block diagonal preconditioner for the system.
+
+    // --- State Flags ---
+    bool m_isFinalized = false;                   ///< Flag indicating if finalize() has been called.
+
+private:
+    /**
+     * @brief Constructs the sparse matrix representations of M, Q, and D, and their reduced forms.
+     * Called during finalize().
+     */
+    void construct_matrix_representations();
+
+    /**
+     * @brief Constructs the block offsets for the reduced system.
+     * Called during finalize().
+     */
+    void construct_reduced_block_offsets();
+
+    /**
+     * @brief Constructs the constant terms of the Jacobian operator (M, -Q, D).
+     * The (0,0) block (gradient of f) is set in GetGradient.
+     * Called during finalize().
+     */
+    void construct_jacobian_constant_terms();
+
+    /**
+     * @brief Scatters the values of essential boundary conditions into the full state vector.
+     * Called during finalize().
+     */
+    void scatter_boundary_conditions();
+
+    /**
+     * @brief Updates the solver for the inverse of the nonlinear Jacobian block (G_00).
+     * @param grad The gradient operator (G_00) of the nonlinear part f(θ).
+     */
+    void update_inverse_nonlinear_jacobian(const mfem::Operator &grad) const;
+
+    /**
+     * @brief Updates the inverse Schur complement operator and its components.
+     * This is typically called after the nonlinear Jacobian part has been updated.
+     */
+    void update_inverse_schur_compliment() const;
+
+    /**
+     * @brief Updates the preconditioner components.
+     * This involves updating the inverse nonlinear Jacobian and then the inverse Schur complement.
+     * @param grad The gradient operator (G_00) of the nonlinear part f(θ).
+     */
+    void update_preconditioner(const mfem::Operator &grad) const;
+};
+
+} // namespace polytrope
+} // namespace serif

src/polytrope/utils/public/utilities.h

@@ -0,0 +1,62 @@
+#pragma once
+
+#include "mfem.hpp"
+
+namespace serif::utilities {
+    [[nodiscard]] mfem::SparseMatrix build_reduced_matrix(
+        const mfem::SparseMatrix& matrix,
+        const mfem::Array<int>& trialEssentialDofs,
+        const mfem::Array<int>& testEssentialDofs
+    );
+
+    /**
+     * @brief Generate a vector of 1s and 0s where 1 elemetns cooresponds to queried dofs. Useful for degugging
+     * @param allDofs array, counding from 0, of all dofs in the system
+     * @param highlightDofs the dofs that you want to identify
+     * @return
+     *
+     * *Example Usage:*
+     * One could use this to identify, for example, which dofs are being identified as the central dofs
+     * @code
+     * ...
+     * mfem::Array<int> phiDofs, thetaDofs;
+     * phiDofs.SetSize(m_fePhi->GetNDofs());
+     * thetaDofs.SetSize(m_feTheta->GetNDofs());
+     * const mfem::Vector phiHighlightVector = serif::utilities::build_dof_identification_vector(phiDofs, phiCenterDofs);
+     * const mfem::Vector thetaHighlightVector = serif::utilities::build_dof_identification_vector(thetaDofs, thetaCenterDofs);
+     * Probe::glVisView(
+     *     const_cast<mfem::Vector&>(phiHighlightVector),
+     *     *m_fePhi,
+     *     "Phi Center Dofs"
+     *     );
+     * Probe::glVisView(
+     *     const_cast<mfem::Vector&>(thetaHighlightVector),
+     *     *m_feTheta,
+     *     "Theta Center Dofs"
+     *     );
+     * @endcode
+     */
+    mfem::Vector build_dof_identification_vector(
+        const mfem::Array<int>& allDofs,
+        const::mfem::Array<int>& highlightDofs
+        );
+
+    /**
+     * @brief Computes the curl of a given H(div) grid function (e.g., from an RT space).
+     *
+     * This function is crucial for diagnosing spurious, non-physical modes in mixed FEM
+     * formulations where the curl of a gradient field is expected to be zero.
+     *
+     * @param phi_gf The GridFunction representing the vector field (e.g., φ). It is expected
+     * to be in an H(div)-conforming space like Raviart-Thomas (RT).
+     * @return A std::pair containing two new grid functions:
+     * - pair.first: A unique_ptr to the vector curl field (∇ × φ). This field will
+     * be in an H(curl)-conforming Nedelec (ND) space.
+     * - pair.second: A unique_ptr to the scalar magnitude of the curl (||∇ × φ||).
+     * This field will be in an L2 space.
+     *
+     * @note The returned unique_ptrs manage the lifetime of the new GridFunctions and their
+     * associated FiniteElementSpaces and FECollections, preventing memory leaks.
+     */
+    mfem::GridFunction compute_curl(mfem::GridFunction& phi_gf);
+}