feat(poly): major work on preconditioner for block form of lane emden equation

working on a "smart" schur compliment preconditioner for the block form of the lane emden equation. Currently this is stub and should not be considered usable

docs/derivations/laneEmdenVariationalBlockForm.aux

@@ -1,9 +1,11 @@
 \relax 
 \@writefile{toc}{\contentsline {section}{\numberline {1}Continuous Variational Form}{1}{}\protected@file@percent }
 \@writefile{toc}{\contentsline {section}{\numberline {2}Discritized Variational Form}{1}{}\protected@file@percent }
-\@writefile{toc}{\contentsline {subsection}{\numberline {2.1}A Few Quick Notes}{3}{}\protected@file@percent }
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.1}The Jacobian}{3}{}\protected@file@percent }
-\@writefile{toc}{\contentsline {section}{\numberline {3}Representation in FEM}{3}{}\protected@file@percent }
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.2}Preconditioner}{4}{}\protected@file@percent }
-\@writefile{toc}{\contentsline {subsection}{\numberline {3.1}MFEM Integrators}{4}{}\protected@file@percent }
+\@writefile{toc}{\contentsline {subsection}{\numberline {2.3}A Few Quick Notes}{5}{}\protected@file@percent }
-\@writefile{lot}{\contentsline {table}{\numberline {1}{\ignorespaces Selection of MFEM Mixed Bilinear Form Integrators}}{4}{}\protected@file@percent }
+\@writefile{toc}{\contentsline {section}{\numberline {3}Representation in FEM}{5}{}\protected@file@percent }
-\newlabel{tab:mfem_mixed_integrators}{{1}{4}{}{table.1}{}}
+\@writefile{toc}{\contentsline {subsection}{\numberline {3.1}MFEM Integrators}{5}{}\protected@file@percent }
-\gdef \@abspage@last{4}
+\@writefile{lot}{\contentsline {table}{\numberline {1}{\ignorespaces Selection of MFEM Mixed Bilinear Form Integrators}}{6}{}\protected@file@percent }
+\newlabel{tab:mfem_mixed_integrators}{{1}{6}{}{table.1}{}}
+\gdef \@abspage@last{6}

docs/derivations/laneEmdenVariationalBlockForm.tex

@@ -122,6 +122,114 @@ We can then set this up as a matrix operation
     \end{bmatrix}
 \end{align}
 
+From this form we can easily see that the residual matrix is
+
+\begin{align}
+    R &= \begin{bmatrix}
+        f(\bar{\theta}) - M\bar{\phi} \\
+        D\bar{\phi} - Q\bar{\theta}
+    \end{bmatrix}
+\end{align}
+
+\subsection{The Jacobian}
+We need to define the Jacobian of this system of equations so that we can use it
+in our Newton-Raphson method. Generally the Jacobian is the matrix of partial derivitives wrt. the state vector. We will let our state vector, $X$, be
+\begin{align}
+    X = \begin{bmatrix}
+        \bar{\theta} \\
+        \bar{\phi}
+    \end{bmatrix}
+\end{align}
+
+So then the Jacobian is
+\begin{align}
+    J &= \begin{bmatrix}
+        \frac{\partial}{\partial \theta}\left(f(\theta) - M\phi\right) & \frac{\partial}{\partial \phi}\left(f(\theta) - M\phi\right) \\
+        \frac{\partial}{\partial \theta}\left(D\phi - Q\theta\right) & \frac{\partial}{\partial \phi}\left(D\phi - Q\theta\right)
+    \end{bmatrix} \\
+    J &= \begin{bmatrix}
+        \frac{df}{d\theta} - \phi\frac{\partial M}{\partial \theta} & -M-\phi\frac{\partial M}{\partial \phi} \\
+        -Q - \theta\frac{\partial Q}{\partial \theta} & D + \phi\frac{\partial D}{\partial \phi} - \theta\frac{\partial Q}{\partial \phi}
+    \end{bmatrix}
+\end{align}
+
+Finally, we know that the matrices $M$, $D$, and $Q$ are constant with respect to $\theta$ and $\phi$. Therefore, we can drop the partial derivatives with respect to $\theta$ and $\phi$ from the Jacobian. This gives us
+\begin{align}
+  \mathbf{J} &= \begin{bmatrix}
+        \frac{df}{d\theta} & -M \\
+        -Q & D
+    \end{bmatrix}
+\end{align}
+
+\noindent In a fully assembled, distritized form this will look like
+
+\begin{align}
+  \mathbf{J} = \begin{bmatrix} \frac{df}{d\theta}_{00}          & \dots  & \frac{df}{d\theta}_{0n_{\theta}}          & -M_{00}                  & \dots            & -M_{0n_{\phi}} \\
+  \vdots          & \ddots &                          & \vdots                   & \ddots           & \\
+    \frac{df}{d\theta}_{n_{\theta}0} &        & \frac{df}{d\theta}_{n_{\theta}n_{\theta}} & -M_{n_{\theta}0}         &                  & -M_{n_{\theta}n_{\phi}} \\
+  -Q_{00}         & \dots  & -Q_{0n_{\theta}}         & D_{00}                   & \dots            & D_{0n_{\phi}} \\
+  \vdots          & \ddots &                          & \vdots                   & \ddots           & \\
+  -Q_{n_{\phi}0}  &        & -Q_{n_{\phi}n_{\theta}}  & D_{n_{\phi}0}            &                  & D_{n_{\phi}n_{\phi}}
+  \end{bmatrix}
+\end{align}
+
+\noindent Where $N_{dof}^{\theta} = n_{\theta}$ is the number of degrees of freedom on $\theta$, and $N_{dof}^{\phi} = n_{\phi}$ is the number of degrees of freedom on $\phi$. Note how the Jacobian is a matrix of size $\left(N_{dof}^{\theta} + N_{dof}^{\phi} \times N_{dof}^{\theta} + N_{dof}^{\phi}\right)$
+
+\subsection{Preconditioner}
+Due to the eventual size of these matrices we would like to be able to solve
+each step in this using a memory efficiet approach. Krylov solvers, such as
+GMRES, allow for matrix free iterative solutions (as long the concept of
+multiplication is defined). However, for these systems to be well formed for
+such solvers it is useful for us to use a preconditioner. However, this is a
+somewhat strongly coupled system where we cannot simply use the inverse
+diagonals of the matrix as a preconditioner. Instead, to encode the coupling we
+will use Schur's Compliment. Each Newton iteration we solve the equation
+\begin{align}
+  \mathbf{J}\Delta \vec{x} = \vec{b}
+\end{align}
+If we expand this out
+\begin{align}
+  \begin{bmatrix} \mathbf{\dot{f}} & -\mathbf{M} \\
+    -\mathbf{Q} & \mathbf{D}
+  \end{bmatrix}\begin{bmatrix} \theta \\
+    \phi
+    \end{bmatrix} = \begin{bmatrix} b_{0} \\
+  b_{1}\end{bmatrix}
+\end{align}
+We can pull out the first equation from this system
+\begin{align}
+  \mathbf{\dot{f}}\theta - \mathbf{M}\phi &= b_{0} \\
+  \theta &= \mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi
+\end{align}
+Then if we pull out the second equation from the system
+\begin{align}
+  -\mathbf{Q}\theta + \mathbf{D}\phi &= b_{1} \\
+  -\mathbf{Q}\left(\mathbf{\dot{f}}^{-1}b_{0} + \mathbf{\dot{f}}^{-1}\mathbf{M}\phi\right) + \mathbf{D}\phi &= b_{1}
+\end{align}
+rearanging terms a bit
+\begin{align}
+  -\mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\phi + \mathbf{D}\phi &= b_{1} \\
+  \left(\mathbf{D} - \mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}\right)\phi &= b_{1} + \mathbf{Q}\mathbf{\dot{f}}^{-1}b_{0}
+\end{align}
+The term $\mathbf{D}-\mathbf{Q}\mathbf{\dot{f}}^{-1}\mathbf{M}$ is Schur's Compliment for this system, and we will represent this by the symbol $\mathbf{\tilde{S}}$. We can use Schur's Compilment to precondition our equation if we let the preconditioner be of the form
+\begin{align}
+  \mathbf{P} = \begin{bmatrix} \mathbf{\dot{f}}^{-1} & 0 \\
+    0 & \mathbf{\tilde{S}}^{-1}
+  \end{bmatrix}
+\end{align}
+So then the preconditioned equation which can be more easily solved by some Krylov solver (such as GMRES) is
+\begin{align}
+  \mathbf{P}\mathbf{J}\Delta \vec{x} = \mathbf{P}\vec{b}
+\end{align}
+
+It is easy to see here that for this system to be solvable / well defined both
+$\mathbf{\tilde{S}}$ and $\mathbf{\dot{f}}$ need to be invertable
+matrices. They are both easily shown to be square (with $\mathbf{\tilde{S}}$
+having a size $\left(N_{dof}^{\phi}\times N_{dof}^{\phi}\right)$ and
+$\mathbf{\dot{f}}$ having a size $\left(N_{dof}^{\theta}\times
+N_{dof}^{\theta}\right)$).
+
+
 \subsection{A Few Quick Notes}
 A few notes on the dimensions of $\mathbf{M}$, $\mathbf{Q}$, $\mathbf{D}$, and $f(\bar{\theta})$.
 \begin{itemize}

src/meson.build

@@ -3,6 +3,7 @@
 # as there are dependencies which exist between them.
 
 # Utility Libraries
+subdir('types')
 subdir('misc')
 subdir('config')
 subdir('probe')

src/poly/solver/meson.build

@@ -35,7 +35,8 @@ dependencies = [
   probe_dep,
   quill_dep,
   config_dep,
-  resourceManager_dep
+  resourceManager_dep,
+  types_dep,
 ]
 
 libPolySolver = static_library('polySolver',

src/poly/solver/private/polySolver.cpp

@@ -21,34 +21,33 @@
 #include "mfem.hpp"
 
 #include <memory>
-#include <string>
 #include <stdexcept>
+#include <string>
 #include <utility>
 
 #include "polySolver.h"
-#include "integrators.h"
+#include "4DSTARTypes.h"
-#include "polyCoeff.h"
 #include "config.h"
+#include "integrators.h"
+#include "operator.h"
+#include "polyCoeff.h"
 #include "probe.h"
 #include "resourceManager.h"
 #include "resourceManagerTypes.h"
-#include "operator.h"
-
-#include "debug.h"
 
 #include "quill/LogMacros.h"
 
 
 namespace laneEmden {
 
-    double a (int k, double n) {
+    double a (int k, double n) { // NOLINT(*-no-recursion)
         if ( k == 0 ) { return 1; }
         if ( k == 1 ) { return 0; }
         else { return -(c(k-2, n)/(std::pow(k, 2)+k)); }
 
     }
 
-    double c(int m, double n) {
+    double c(int m, double n) { // NOLINT(*-no-recursion)
         if ( m == 0 ) { return std::pow(a(0, n), n); }
         else {
             double termOne = 1.0/(m*a(0, n));
@@ -60,7 +59,7 @@ namespace laneEmden {
         }
     }
 
-    double thetaSerieseExpansion(double xi, double n, int order) {
+    double thetaSeriesExpansion(double xi, double n, int order) {
         double acc = 0;
         for (int k = 0; k < order; k++) {
             acc += a(k, n) * std::pow(xi, k);
@@ -101,7 +100,7 @@ PolySolver::PolySolver(double n, double order) {
 
 }
 
-PolySolver::~PolySolver() {}
+PolySolver::~PolySolver() = default;
 
 void PolySolver::assembleBlockSystem() {
 
@@ -121,20 +120,8 @@ void PolySolver::assembleBlockSystem() {
     blockOffsets[2] = feSpaces[1]->GetVSize();
     blockOffsets.PartialSum();
 
-    // Coefficients
-    mfem::ConstantCoefficient negOneCoeff(-1.0);
-    mfem::ConstantCoefficient oneCoeff(1.0);
-
-    mfem::Vector negOneVec(mfem::Vector(3));
-    mfem::Vector oneVec(mfem::Vector(3));
-    negOneVec = -1.0;
-    oneVec = 1.0;
-
-    mfem::VectorConstantCoefficient negOneVCoeff(negOneVec);
-    mfem::VectorConstantCoefficient oneVCoeff(oneVec);
-
     // Add integrators to block form one by one
-    // We add integrators cooresponding to each term in the weak form
+    // We add integrators corresponding to each term in the weak form
     // The block form of the residual matrix
     // ⎡ 0 -M ⎤ ⎡ θ ⎤ + ⎡f(θ)⎤ = ⎡ 0 ⎤ = R(X)
     // ⎣ -Q D ⎦ ⎣ Φ ⎦   ⎣  0 ⎦   ⎣ 0 ⎦
@@ -171,6 +158,7 @@ void PolySolver::assembleBlockSystem() {
 
     // --- Assemble the NonlinearForm (f) ---
     auto fform = std::make_unique<mfem::NonlinearForm>(m_feTheta.get());
+    mfem::ConstantCoefficient oneCoeff(1.0);
     fform->AddDomainIntegrator(new polyMFEMUtils::NonlinearPowerIntegrator(oneCoeff, m_polytropicIndex));
     // TODO: Add essential boundary conditions to the nonlinear form
 
@@ -185,7 +173,7 @@ void PolySolver::assembleBlockSystem() {
 
 }
 
-void PolySolver::solve(){
+void PolySolver::solve() const {
     // --- Set the initial guess for the solution ---
     setInitialGuess();
 
@@ -203,6 +191,14 @@ void PolySolver::solve(){
 
     mfem::NewtonSolver newtonSolver = setupNewtonSolver();
 
+    // EMB 2025: Calling Mult gets the gradient of the operator for the NewtonSolver
+    // This then is set as the operator for the solver for NewtonSolver
+    // The solver (assuming it is an iterative solver) then sets the
+    // operator for its preconditioner to this.
+    // What this means is that there is no need to manually deal
+    // with updating the preconditioner at every newton step as the
+    // changes to the jacobian are automatically propagated through the
+    // solving chain. This is at least true with MFEM 4.8-rc0
     newtonSolver.Mult(zero_rhs, state_vector);
 
     // --- Save and view the solution ---
@@ -211,25 +207,40 @@ void PolySolver::solve(){
 
 }
 
-std::pair<mfem::Array<int>, mfem::Array<int>> PolySolver::getEssentialTrueDof() {
+SSE::MFEMArrayPairSet PolySolver::getEssentialTrueDof() const {
     mfem::Array<int> theta_ess_tdof_list;
     mfem::Array<int> phi_ess_tdof_list;
 
-    mfem::Array<int> centerDofs = findCenterElement();
+    mfem::Array<int> thetaCenterDofs, phiCenterDofs;
+    mfem::Array<double> thetaCenterVals, phiCenterVals;
+    std::tie(thetaCenterDofs, phi_ess_tdof_list) = findCenterElement();
+    thetaCenterVals.SetSize(thetaCenterDofs.Size());
+    phiCenterVals.SetSize(phi_ess_tdof_list.Size());
 
-    phi_ess_tdof_list.Append(centerDofs);
+    thetaCenterVals = 1.0;
+    phiCenterVals = 0.0;
 
     mfem::Array<int> ess_brd(m_mesh->bdr_attributes.Max());
     ess_brd = 1;
-    m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list);
-    // combine the essential dofs with the center dofs
-    theta_ess_tdof_list.Append(centerDofs);
 
-    return std::make_pair(theta_ess_tdof_list, phi_ess_tdof_list);
+    mfem::Array<double> thetaSurfaceVals;
+    m_feTheta->GetEssentialTrueDofs(ess_brd, theta_ess_tdof_list);
+    thetaSurfaceVals.SetSize(theta_ess_tdof_list.Size());
+    thetaSurfaceVals = 0.0;
+    // combine the essential dofs with the center dofs
+    theta_ess_tdof_list.Append(thetaCenterDofs);
+    thetaSurfaceVals.Append(thetaCenterVals);
+
+    SSE::MFEMArrayPair thetaPair = std::make_pair(theta_ess_tdof_list, thetaSurfaceVals);
+    SSE::MFEMArrayPair phiPair = std::make_pair(phi_ess_tdof_list, phiCenterVals);
+    SSE::MFEMArrayPairSet pairSet = std::make_pair(thetaPair, phiPair);
+
+    return pairSet;
 }
 
-mfem::Array<int> PolySolver::findCenterElement() {
+std::pair<mfem::Array<int>, mfem::Array<int>> PolySolver::findCenterElement() const {
-    mfem::Array<int> centerDofs;
+    mfem::Array<int> thetaCenterDofs;
+    mfem::Array<int> phiCenterDofs;
     mfem::DenseMatrix centerPoint(m_mesh->SpaceDimension(), 1);
     centerPoint(0, 0) = 0.0;
     centerPoint(1, 0) = 0.0;
@@ -241,12 +252,14 @@ mfem::Array<int> PolySolver::findCenterElement() {
     mfem::Array<int> tempDofs;
     for (int i = 0; i < elementIDs.Size(); i++) {
         m_feTheta->GetElementDofs(elementIDs[i], tempDofs);
-        centerDofs.Append(tempDofs);
+        thetaCenterDofs.Append(tempDofs);
+        m_fePhi->GetElementDofs(elementIDs[i], tempDofs);
+        phiCenterDofs.Append(tempDofs);
     }
-    return centerDofs;
+    return std::make_pair(thetaCenterDofs, phiCenterDofs);
 }
 
-void PolySolver::setInitialGuess() {
+void PolySolver::setInitialGuess() const {
     // --- Set the initial guess for the solution ---
     mfem::FunctionCoefficient thetaInitGuess (
         [this](const mfem::Vector &x) {
@@ -261,39 +274,38 @@ void PolySolver::setInitialGuess() {
         [this](const mfem::Vector &x, mfem::Vector &v) {
             double radius = Probe::getMeshRadius(*m_mesh);
             double u = -1/std::pow(radius,2);
-            v(0) = 2*std::abs(x(0))*u;
+            v(0) = 2*x(0)*u;
-            v(1) = 2*std::abs(x(1))*u;
+            v(1) = 2*x(1)*u;
-            v(2) = 2*std::abs(x(2))*u;
+            v(2) = 2*x(2)*u;
         }
     );
     m_theta->ProjectCoefficient(thetaInitGuess);
     m_phi->ProjectCoefficient(phiInitGuess);
     if (m_config.get<bool>("Poly:Solver:ViewInitialGuess", false)) {
-        Probe::glVisView(*m_theta, *m_mesh, "initialGuess");
+        Probe::glVisView(*m_theta, *m_mesh, "θ init");
+        Probe::glVisView(*m_phi, *m_mesh, "ɸ init");
     }
 
 }
 
-void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) {
+void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) const {
     mfem::BlockVector x_block(const_cast<mfem::BlockVector&>(state_vector), m_polytropOperator->GetBlockOffsets());
     mfem::Vector& x_theta = x_block.GetBlock(0);
     mfem::Vector& x_phi = x_block.GetBlock(1);
 
-    bool doView = m_config.get<bool>("Poly:Output:View", false);
+    if (m_config.get<bool>("Poly:Output:View", false)) {
-    if (doView) {
         Probe::glVisView(x_theta, *m_feTheta, "θ Solution");
         Probe::glVisView(x_phi, *m_fePhi, "ɸ Solution");
     }
 
     // --- Extract the Solution ---
-    bool write11DSolution = m_config.get<bool>("Poly:Output:1D:Save", true);
+    if (bool write11DSolution = m_config.get<bool>("Poly:Output:1D:Save", true)) {
-    if (write11DSolution) {
+        auto solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
-        std::string solutionPath = m_config.get<std::string>("Poly:Output:1D:Path", "polytropeSolution_1D.csv");
+        auto derivSolPath = "d" + solutionPath;
-        std::string derivSolPath = "d" + solutionPath;
 
-        double rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
+        auto rayCoLatitude = m_config.get<double>("Poly:Output:1D:RayCoLatitude", 0.0);
-        double rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
+        auto rayLongitude = m_config.get<double>("Poly:Output:1D:RayLongitude", 0.0);
-        int raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
+        auto raySamples = m_config.get<int>("Poly:Output:1D:RaySamples", 100);
 
         std::vector rayDirection = {rayCoLatitude, rayLongitude};
 
@@ -302,10 +314,10 @@ void PolySolver::saveAndViewSolution(const mfem::BlockVector& state_vector) {
     }
 }
 
-void PolySolver::setupOperator() {
+void PolySolver::setupOperator() const {
-    mfem::Array<int> theta_ess_tdof_list, phi_ess_tdof_list;
+
-    std::tie(theta_ess_tdof_list, phi_ess_tdof_list) = getEssentialTrueDof();
+    SSE::MFEMArrayPairSet ess_tdof_pair_set = getEssentialTrueDof();
-    m_polytropOperator->SetEssentialTrueDofs(theta_ess_tdof_list, phi_ess_tdof_list);
+    m_polytropOperator->SetEssentialTrueDofs(ess_tdof_pair_set);
 
     // -- Finalize the operator --
     m_polytropOperator->finalize();
@@ -316,20 +328,52 @@ void PolySolver::setupOperator() {
     }
 }
 
-mfem::NewtonSolver PolySolver::setupNewtonSolver(){
+void PolySolver::LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel, double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const {
-    // --- Load configuration parameters ---
+    newtonRelTol = m_config.get<double>("Poly:Solver:Newton:RelTol", 1e-7);
-    double newtonRelTol = m_config.get<double>("Poly:Solver:Newton:RelTol", 1e-7);
+    newtonAbsTol = m_config.get<double>("Poly:Solver:Newton:AbsTol", 1e-7);
-    double newtonAbsTol = m_config.get<double>("Poly:Solver:Newton:AbsTol", 1e-7);
+    newtonMaxIter = m_config.get<int>("Poly:Solver:Newton:MaxIter", 200);
-    int newtonMaxIter = m_config.get<int>("Poly:Solver:Newton:MaxIter", 200);
+    newtonPrintLevel = m_config.get<int>("Poly:Solver:Newton:PrintLevel", 1);
-    int newtonPrintLevel = m_config.get<int>("Poly:Solver:Newton:PrintLevel", 1);
 
-    double gmresRelTol = m_config.get<double>("Poly:Solver:GMRES:RelTol", 1e-10);
+    gmresRelTol = m_config.get<double>("Poly:Solver:GMRES:RelTol", 1e-10);
-    double gmresAbsTol = m_config.get<double>("Poly:Solver:GMRES:AbsTol", 1e-12);
+    gmresAbsTol = m_config.get<double>("Poly:Solver:GMRES:AbsTol", 1e-12);
-    int gmresMaxIter = m_config.get<int>("Poly:Solver:GMRES:MaxIter", 2000);
+    gmresMaxIter = m_config.get<int>("Poly:Solver:GMRES:MaxIter", 2000);
-    int gmresPrintLevel = m_config.get<int>("Poly:Solver:GMRES:PrintLevel", 0);
+    gmresPrintLevel = m_config.get<int>("Poly:Solver:GMRES:PrintLevel", 0);
 
     LOG_DEBUG(m_logger, "Newton Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel);
     LOG_DEBUG(m_logger, "GMRES Solver (relTol: {:0.2E}, absTol: {:0.2E}, maxIter: {}, printLevel: {})", gmresRelTol, gmresAbsTol, gmresMaxIter, gmresPrintLevel);
+}
+
+mfem::BlockDiagonalPreconditioner PolySolver::build_preconditioner() const {
+    // --- Set up the preconditioners. The non-linear form will use a Chebyshev Preconditioner while the linear terms will use a simpler Jacobi preconditioner ---
+    mfem::BlockDiagonalPreconditioner prec(m_polytropOperator->GetBlockOffsets());
+    const mfem::BlockOperator &jacobian = m_polytropOperator->GetJacobianOperator();
+
+    // Get all the blocks. J00 -> Non-linear form (df(θ)/dθ), J01-> -M, J10 -> -Q, J11 -> D
+    const mfem::Operator& J00 = jacobian.GetBlock(0, 0);
+    mfem::Vector J00diag;
+
+    J00.AssembleDiagonal(J00diag);
+
+    SSE::MFEMArrayPairSet ess_tdof_pair_set = m_polytropOperator->GetEssentialTrueDofs();
+
+    // TODO: This order may need to be tuned (EMB)
+    // --- ess_tdof_pair_set.first -> (theta dof ids, theta dof vals).first -> theta dof ids
+    // --- ess_tdof_pair_set.second -> (phi dof ids, phi dof vals).first -> phi dof ids
+    mfem::OperatorChebyshevSmoother J00Prec(J00, J00diag, ess_tdof_pair_set.first.first, 2);
+    mfem::OperatorJacobiSmoother J11Prec(m_polytropOperator->GetJ11diag(), ess_tdof_pair_set.second.first);
+
+    prec.SetDiagonalBlock(0, &J00Prec);
+    prec.SetDiagonalBlock(1, &J11Prec);
+
+    return prec;
+}
+
+mfem::NewtonSolver PolySolver::setupNewtonSolver() const {
+    // --- Load configuration parameters ---
+    double newtonRelTol, newtonAbsTol, gmresRelTol, gmresAbsTol;
+    int newtonMaxIter, newtonPrintLevel, gmresMaxIter, gmresPrintLevel;
+    LoadSolverUserParams(newtonRelTol, newtonAbsTol, newtonMaxIter, newtonPrintLevel, gmresRelTol, gmresAbsTol,
+                         gmresMaxIter, gmresPrintLevel);
 
     // --- Set up the Newton solver ---
     mfem::NewtonSolver newtonSolver;
@@ -338,13 +382,18 @@ mfem::NewtonSolver PolySolver::setupNewtonSolver(){
     newtonSolver.SetMaxIter(newtonMaxIter);
     newtonSolver.SetPrintLevel(newtonPrintLevel);
     newtonSolver.SetOperator(*m_polytropOperator);
+
+    // --- Created the linear solver which is used to invert the jacobian ---
     mfem::GMRESSolver gmresSolver;
     gmresSolver.SetRelTol(gmresRelTol);
     gmresSolver.SetAbsTol(gmresAbsTol);
     gmresSolver.SetMaxIter(gmresMaxIter);
     gmresSolver.SetPrintLevel(gmresPrintLevel);
+
+    // build_preconditioner();
+
+
     newtonSolver.SetSolver(gmresSolver);
-    // newtonSolver.SetAdaptiveLinRtol();
 
     return newtonSolver;
 }

src/poly/solver/public/polySolver.h

@@ -1,21 +1,22 @@
 #ifndef POLYSOLVER_H
 #define POLYSOLVER_H
 
-#include "linalg/solvers.hpp"
 #include "mfem.hpp"
 #include <memory>
 #include <utility>
 
 #include "integrators.h"
+#include "4DSTARTypes.h"
 #include "operator.h"
 #include "config.h"
 #include "probe.h"
 #include "quill/Logger.h"
 
+
 namespace laneEmden {
     double a (int k, double n);
     double c(int m, double n);
-    double thetaSerieseExpansion(double xi, double n, int order);
+    double thetaSeriesExpansion(double xi, double n, int order);
 }
 
 class PolySolver {
@@ -23,7 +24,7 @@ public: // Public methods
     PolySolver(double n, double order);
     ~PolySolver();
 
-    void solve();
+    void solve() const;
 
     double getN() { return m_polytropicIndex; }
     double getOrder() { return m_feOrder; }
@@ -48,17 +49,22 @@ private: // Private Attributes
 
     std::unique_ptr<PolytropeOperator> m_polytropOperator;
 
+    std::unique_ptr<mfem::OperatorJacobiSmoother> m_prec;
+
 
 private: // Private methods
     void assembleBlockSystem();
-    std::pair<mfem::Array<int>, mfem::Array<int>> getEssentialTrueDof();
+    SSE::MFEMArrayPairSet getEssentialTrueDof() const;
-    mfem::Array<int> findCenterElement();
+    std::pair<mfem::Array<int>, mfem::Array<int>> findCenterElement() const;
-    void setInitialGuess();
+    void setInitialGuess() const;
-    void saveAndViewSolution(const mfem::BlockVector& state_vector);
+    void saveAndViewSolution(const mfem::BlockVector& state_vector) const;
-    mfem::NewtonSolver setupNewtonSolver();
+    mfem::NewtonSolver setupNewtonSolver() const;
-    void setupOperator();
+    void setupOperator() const;
 
+    void LoadSolverUserParams(double &newtonRelTol, double &newtonAbsTol, int &newtonMaxIter, int &newtonPrintLevel,
+                              double &gmresRelTol, double &gmresAbsTol, int &gmresMaxIter, int &gmresPrintLevel) const;
 
+    mfem::BlockDiagonalPreconditioner build_preconditioner() const;
 };
 
 #endif // POLYSOLVER_H

src/poly/utils/meson.build

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

src/poly/utils/private/operator.cpp

@@ -1,9 +1,71 @@
 #include "operator.h"
+#include "4DSTARTypes.h"
+#include "linalg/densemat.hpp"
+#include "linalg/sparsemat.hpp"
 #include "mfem.hpp"
 #include "linalg/vector.hpp"
 #include <memory>
 
-#include "debug.h"
+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,
+    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);
+}
 
 PolytropeOperator::PolytropeOperator(
 
@@ -31,12 +93,81 @@ void PolytropeOperator::finalize() {
     m_Qmat = std::make_unique<mfem::SparseMatrix>(m_Q->SpMat());
     m_Dmat = std::make_unique<mfem::SparseMatrix>(m_D->SpMat());
 
+    for (const int thetaDof : m_theta_ess_tdofs.first) {
+        m_Mmat->EliminateRow(thetaDof);
+        m_Qmat->EliminateCol(thetaDof);
+    }
+    // These are commented out because they theoretically are wrong (need to think more about how to apply essential dofs to a vector div field)
+    // for (const int phiDof : m_phi_ess_tdofs.first) {
+    //     if (phiDof >=0 && phiDof < m_Dmat->Height()) {
+    //         m_Dmat->EliminateRowCol(phiDof);
+    //         m_Qmat->EliminateRow(phiDof);
+    //         m_Mmat->EliminateCol(phiDof);
+    //     }
+    // }
+
     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);
 
+    // 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(1, 0, m_negQ_op.get()); // -Q (constant)
+    m_jacobian->SetBlock(1, 1, m_Dmat.get()); // D (constant)
+
     m_isFinalized = true;
 }
 
+const mfem::BlockOperator &PolytropeOperator::GetJacobianOperator() const {
+     if (m_jacobian == nullptr) {
+         MFEM_ABORT("Jacobian has not been initialized before GetJacobianOperator() call.");
+     }
+     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) {
@@ -60,7 +191,7 @@ void PolytropeOperator::Mult(const mfem::Vector &x, mfem::Vector &y) const {
     mfem::Vector Dphi_term(phi_size);
     mfem::Vector Qtheta_term(phi_size);
 
-    // Caucluate R0 and R1 terms
+    // Calculate R0 and R1 terms
     // R0 = f(θ) - Mɸ
     // R1 = Dɸ - Qθ
 
@@ -78,22 +209,75 @@ 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_tofs.Size(); i++) {
+    for (int i = 0; i < m_theta_ess_tdofs.first.Size(); i++) {
-        int idx = m_theta_ess_tofs[i];
+        int idx = m_theta_ess_tdofs.first[i];
         if (idx >= 0 && idx < y_R0.Size()) {
-            y_block.GetBlock(0)[idx] = 0.0; // Zero out the essential theta dofs in the bilinear form
+            const double &targetValue = m_theta_ess_tdofs.second[i];
+            y_block.GetBlock(0)[idx] = targetValue - x_theta(idx); // 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
         }
     }
 
-    for (int i = 0; i < m_phi_ess_tofs.Size(); i++) {
+    // TODO look into how the true dof -> vector component works
-        int idx = m_phi_ess_tofs[i];
+    // for (int i = 0; i < m_phi_ess_tdofs.first.Size(); i++) {
-        if (idx >= 0 && idx < y_R1.Size()) {
+    //     int idx = m_phi_ess_tdofs.first[i];
-            y_block.GetBlock(1)[idx] = 0.0; // Zero out the essential phi dofs in the bilinear form
+    //     if (idx >= 0 && idx < y_R1.Size()) {
+    //         // const double &targetValue = m_phi_ess_tdofs.second[i];
+    //         // y_block.GetBlock(1)[idx] = targetValue - x_phi(idx); // Zero out the essential phi dofs in the bi-linear form
+    //         y_block.GetBlock(1)[idx] = 0; // Zero out the essential phi dofs in the bi-linear form
+    //     }
+    // }
+
+}
+
+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);
+    } else {
+        MFEM_ABORT("PolytropeOperator::GetGradient called on nonlinear Jacobian");
     }
 }
 
+void PolytropeOperator::updateInverseSchurCompliment() const {
+    // TODO Add a flag in to make sure this tracks in parallel (i.e. every time the non linear jacobian inverse is updated set the flag to true and then check if the flag is true here and if so do work (if not throw error). then at the end of this function set it to false.
+    if (m_invNonlinearJacobian == nullptr) {
+        MFEM_ABORT("PolytropeOperator::updateInverseSchurCompliment called before updateInverseNonlinearJacobian");
     }
+    mfem::SparseMatrix* schurCompliment(m_Dmat.get()); // This is now a copy of D
+
+    // Calculate S = D - Q df^{-1} M
+    mfem::SparseMatrix* temp_QxdfInv = mfem::Mult(*m_Qmat, *m_invNonlinearJacobian); // Q * df^{-1}
+    mfem::SparseMatrix* temp_QxdfInvxM = mfem::Mult(*temp_QxdfInv, *m_Mmat); // Q * df^{-1} * M
+
+    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
+     * I think I have this close; however, when running there is a seg fault when converting to a DenseMatrix for
+     *   writing out to CSV. This bodes poorly for the underlying memory structure.
+     *
+     * Also, there is a lot of math happening in these tight loops, probably need to think about that and
+     *   make it more efficient.
+     */
+
+    // TODO Also I need to actually invert the schur compliment since so far I have just found its non inverted state
+    // Should probably add an energy check to make sure it is roughly diagonal
+
+    // I did this manually in python for the non linear jacobian and the energy along the diagonal was ~99.999 % of the
+    // total energy indicating that it is very strongly diagonal.
+}
+
 mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
     if (!m_isFinalized) {
         MFEM_ABORT("PolytropeOperator::GetGradient called before finalize");
@@ -102,31 +286,34 @@ mfem::Operator& PolytropeOperator::GetGradient(const mfem::Vector &x) const {
     mfem::BlockVector x_block(const_cast<mfem::Vector&>(x), m_blockOffsets);
     const mfem::Vector& x_theta = x_block.GetBlock(0);
 
-    mfem::Operator& J00 = m_f->GetGradient(x_theta);
+    auto &grad = m_f->GetGradient(x_theta);
-    if (m_jacobian == nullptr) {
+    updateInverseNonlinearJacobian(grad);
-        m_jacobian = std::make_unique<mfem::BlockOperator>(m_blockOffsets);
+    updateInverseSchurCompliment();
-        m_jacobian->SetBlock(0, 0, &J00); // df/dθ (state-dependent)
+
-        m_jacobian->SetBlock(0, 1, m_negM_op.get()); // -M (constant)
+    m_jacobian->SetBlock(0, 0, &grad);
-        m_jacobian->SetBlock(1, 0, m_negQ_op.get()); // -Q (constant)
+    // The other blocks are set up in finalize since they are constant. Only J00 depends on the current state of theta
-        m_jacobian->SetBlock(1, 1, m_Dmat.get()); // D (constant)
-    } else {
-        // The Jacobian already exists, we only need to update the first block
-        // since the other blocks have a constant derivitive (they are linear)
-        m_jacobian->SetBlock(0, 0, &J00);
-    }
 
     return *m_jacobian;
 }
 
-void PolytropeOperator::SetEssentialTrueDofs(const mfem::Array<int> &theta_ess_tofs,
+void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs) {
-                                             const mfem::Array<int> &phi_ess_tofs) {
     m_isFinalized = false;
-    m_theta_ess_tofs = theta_ess_tofs;
+    m_theta_ess_tdofs = theta_ess_tdofs;
-    m_phi_ess_tofs = phi_ess_tofs;
+    m_phi_ess_tdofs = phi_ess_tdofs;
 
     if (m_f) {
-        m_f->SetEssentialTrueDofs(theta_ess_tofs);
+        m_f->SetEssentialTrueDofs(theta_ess_tdofs.first);
+
+        // This should be zeroing out the row; however, I am getting a segfault
     } else {
         MFEM_ABORT("m_f is null in PolytropeOperator::SetEssentialTrueDofs");
     }
 }
+
+void PolytropeOperator::SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set) {
+    SetEssentialTrueDofs(ess_tdof_pair_set.first, ess_tdof_pair_set.second);
+}
+
+SSE::MFEMArrayPairSet PolytropeOperator::GetEssentialTrueDofs() const {
+     return std::make_pair(m_theta_ess_tdofs, m_phi_ess_tdofs);
+}

src/poly/utils/public/operator.h

@@ -2,8 +2,11 @@
 #define POLY_UTILS_OPERATOR_H
 
 #include "mfem.hpp"
+#include "4DSTARTypes.h"
 #include <memory>
 
+#include "probe.h"
+
 class PolytropeOperator : public mfem::Operator {
 public:
     PolytropeOperator(
@@ -15,10 +18,14 @@ public:
     ~PolytropeOperator() override = default;
 
     void Mult(const mfem::Vector &x, mfem::Vector &y) const override;
+
+
     mfem::Operator& GetGradient(const mfem::Vector &x) const override;
 
-    void SetEssentialTrueDofs(const mfem::Array<int> &theta_ess_tofs,
+    void SetEssentialTrueDofs(const SSE::MFEMArrayPair& theta_ess_tdofs, const SSE::MFEMArrayPair& phi_ess_tdofs);
-                            const mfem::Array<int> &phi_ess_tofs);
+    void SetEssentialTrueDofs(const SSE::MFEMArrayPairSet& ess_tdof_pair_set);
+
+    SSE::MFEMArrayPairSet GetEssentialTrueDofs() const;
 
     bool isFinalized() const { return m_isFinalized; }
 
@@ -26,7 +33,18 @@ public:
 
     const mfem::Array<int>& GetBlockOffsets() const { return m_blockOffsets; }
 
+    const mfem::BlockOperator &GetJacobianOperator() const;
+
+    // Get the diagonals of -M, -Q, and D. I use J01, J10, and J11 for those here since the diagonals
+    // are not just from the matrixes, but are the scaled versions in the jacobian specificially (-1*M & -1*Q)
+    mfem::Vector GetJ00Diag() const;
+    mfem::Vector GetJ01Diag() const;
+    mfem::Vector GetJ10diag() const;
+    mfem::Vector GetJ11diag() const;
+
 private:
+    Probe::LogManager& m_logManager = Probe::LogManager::getInstance();
+    quill::Logger* m_logger = m_logManager.getLogger("log");
     std::unique_ptr<mfem::MixedBilinearForm> m_M;
     std::unique_ptr<mfem::MixedBilinearForm> m_Q;
     std::unique_ptr<mfem::BilinearForm> m_D;
@@ -34,19 +52,39 @@ private:
 
     const mfem::Array<int> m_blockOffsets;
 
-    mfem::Array<int> m_theta_ess_tofs;
+    SSE::MFEMArrayPair m_theta_ess_tdofs;
-    mfem::Array<int> m_phi_ess_tofs;
+    SSE::MFEMArrayPair m_phi_ess_tdofs;
 
     std::unique_ptr<mfem::SparseMatrix> m_Mmat;
     std::unique_ptr<mfem::SparseMatrix> m_Qmat;
     std::unique_ptr<mfem::SparseMatrix> 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;
 
+    // TODO I think these need to be calculated in the GetGradient every time since they will always change
+    mutable std::unique_ptr<mfem::SparseMatrix> m_invSchurCompliment;
+    mutable std::unique_ptr<mfem::SparseMatrix> m_invNonlinearJacobian;
+
+    /*
+     * The schur preconditioner has the form
+     *
+     * | df/drtheta^-1 0    |
+     * | 0             S^-1 |
+     *
+     * Where S is the Schur compliment of the system
+     *
+     */
+
+    // TODO: I have not combined these parts yet and they need to be combined
+    mutable std::unique_ptr<mfem::SparseMatrix> m_schurPreconditioner;
+
     bool m_isFinalized = false;
+
+private:
+    void updateInverseNonlinearJacobian(const mfem::Operator &grad) const;
+    void updateInverseSchurCompliment() const;
 };
 
 

src/resource/private/resourceManagerTypes.cpp

@@ -22,7 +22,7 @@ std::string getFirstSegment(const std::string& input) {
 Resource createResource(const std::string& type, const std::string& path) {
     static const std::unordered_map<std::string, std::function<Resource(const std::string&)>> factoryMap = {
         {"opac", [](const std::string& p) { return Resource(
-            std::make_unique<OpatIO>(p));
+            std::make_unique<opat::OPAT>(opat::readOPAT(p)));
          }},
         {"mesh", [](const std::string& p) { return Resource(
             std::make_unique<MeshIO>(p));

src/resource/public/resourceManagerTypes.h

@@ -30,7 +30,7 @@
  * @endcode
  */
 using Resource = std::variant<
-                    std::unique_ptr<OpatIO>,
+                    std::unique_ptr<opat::OPAT>,
                     std::unique_ptr<MeshIO>,
                     std::unique_ptr<EosIO>>;
 

src/types/meson.build

@@ -0,0 +1 @@
+types_dep = declare_dependency(include_directories: include_directories('public'))

src/types/public/4DSTARTypes.h

@@ -0,0 +1,13 @@
+#ifndef _4DSTAR_TYPES_H
+#define _4DSTAR_TYPES_H
+
+#include <utility>
+#include "mfem.hpp"
+
+// TODO : Need a better namespace name for these types
+namespace SSE {
+    typedef std::pair<mfem::Array<int>, mfem::Array<double>> MFEMArrayPair;
+    typedef std::pair<MFEMArrayPair, MFEMArrayPair> MFEMArrayPairSet;
+}
+
+#endif // _4DSTAR_TYPES_H