diff --git a/EmilysRulesForSERiFCode.md b/EmilysRulesForSERiFCode.md
index 28024f61..3f63a3b7 100644
--- a/EmilysRulesForSERiFCode.md
+++ b/EmilysRulesForSERiFCode.md
@@ -37,4 +37,6 @@ for correctness.
 - Prefer more and smaller functions over fewer and larger functions.
 - Do not include the full mfem.hpp header. Its slow, I will send you an angry email. 
 - AI is a tool, its okay to use it, but use it as a tool not as a crutch. You must understand all the code you write.
-- The less something is related to physics the more you may consider AI. A build system configuration for example is a good candidate for AI. A physics function is not.
\ No newline at end of file
+- The less something is related to physics the more you may consider AI. A build system configuration for example is a good candidate for AI. A physics function is not.
+- Its not 1980, we are not working on dumb terminals connected to a PDP 11 with 80 charecters of width. Do not limit yourself to 80 charecter width lines if it limits the readability of that line
+- At the same time monitors do still have a finite width, as a line gets long consider if it makes sense to split it into multiple lines.
\ No newline at end of file
diff --git a/build-config/magicenum/meson.build b/build-config/magicenum/meson.build
index d68252b4..e13e0b54 100644
--- a/build-config/magicenum/meson.build
+++ b/build-config/magicenum/meson.build
@@ -1,6 +1,9 @@
 magic_enum_include = include_directories('include')
 
-magic_enum_args = []
+# SERiFErrorCode includes values up to 1000; magic_enum defaults to 127.
+# Keep this definition on the public dependency so every translation unit
+# instantiates the same reflection tables.
+magic_enum_args = ['-DMAGIC_ENUM_RANGE_MAX=2000']
 
 if get_option('magic_enum_hash')
     magic_enum_args += '-DMAGIC_ENUM_ENABLE_HASH'
@@ -11,3 +14,4 @@ magic_enum_dep = declare_dependency(
     compile_args: magic_enum_args,
 )
 
+meson.override_dependency('magic_enum', magic_enum_dep)
diff --git a/build-config/mfem/meson.build b/build-config/mfem/meson.build
index ce0a82d5..044f550e 100644
--- a/build-config/mfem/meson.build
+++ b/build-config/mfem/meson.build
@@ -406,11 +406,12 @@ if effective_allow_preinstalled
     endforeach
 
     if system_compatible
-      system_compatible = cpp.compiles(
+      system_compatible = cpp.links(
         '#include <mfem.hpp>\n' +
-        'int main(){return 0;}',
+        '#include <sstream>\n' +
+        'int main(){std::istringstream input; mfem::Mesh mesh(input); return mesh.Dimension();}',
         dependencies: system_candidate_dependencies,
-        name: 'system MFEM public headers are consumable',
+        name: 'system MFEM matches the selected C++ ABI',
       )
     endif
 
@@ -597,6 +598,15 @@ else
     bundle_command += ['--c-link-arg=' + platform_link_arg]
     bundle_command += ['--cxx-link-arg=' + platform_link_arg]
   endforeach
+  # Meson does not include the native-file c_ld/cpp_ld selection in *_link_args.
+  # Carry that selection into the external CMake/Make dependency builds too.
+  linker_flags = {'ld.bfd': '-fuse-ld=bfd', 'ld.gold': '-fuse-ld=gold', 'ld.lld': '-fuse-ld=lld'}
+  if linker_flags.has_key(cc.get_linker_id())
+    bundle_command += ['--c-link-arg=' + linker_flags[cc.get_linker_id()]]
+  endif
+  if linker_flags.has_key(cpp.get_linker_id())
+    bundle_command += ['--cxx-link-arg=' + linker_flags[cpp.get_linker_id()]]
+  endif
   foreach i : range(mfem_feature_names.length())
     bundle_command += [
       '--feature',
diff --git a/meson_options.txt b/meson_options.txt
index 395eaa6c..651963d1 100644
--- a/meson_options.txt
+++ b/meson_options.txt
@@ -3,6 +3,7 @@ option('feature_profile', type: 'combo', choices: ['minimal', 'portable', 'full'
 option('build_python', type: 'boolean', value: false, description: 'Build the nanobind Python extension and install native assets inside the Python package.')
 option('build_tests', type: 'boolean', value: true, description: 'Register serial, parallel, and Python smoke tests when applicable.')
 option('build_checks', type: 'boolean', value: true, description: 'Build the compile-time checks for SERiF type invariants.')
+option('sandbox_expected_errors', type: 'array', value: ['finite_sandbox:NON_FINITE_VALUE_ERROR', 'quadrature_sandbox:INVALID_QUADRATURE_ORDER'], description: 'Intentional error examples, as sandbox_name:ERROR_CODE. Other sandboxes must exit successfully.')
 option('install_mfem', type: 'boolean', value: true, description: 'Install the source-built MFEM bundle and its runtime dependencies.')
 option('dependency_prefix', type: 'string', value: '', description: 'Explicit prefix containing optional TPLs. Used even when allow_preinstalled=false.')
 option('jobs', type: 'integer', min: 1, max: 1024, value: 4, description: 'Parallel jobs for the CMake/Make source bundle.')
diff --git a/native/clang.ini b/native/clang.ini
new file mode 100644
index 00000000..21977ae9
--- /dev/null
+++ b/native/clang.ini
@@ -0,0 +1,3 @@
+[binaries]
+c = 'clang'
+cpp = 'clang++'
diff --git a/native/gcc-libcxx.ini b/native/gcc-libcxx.ini
new file mode 100644
index 00000000..a9de4787
--- /dev/null
+++ b/native/gcc-libcxx.ini
@@ -0,0 +1,5 @@
+# Linux GCC 16+ with libc++ development headers in the standard system location.
+# Layer after gcc.ini. Override this include path for a non-system LLVM install.
+[built-in options]
+cpp_args = ['-nostdinc++', '-isystem', '/usr/include/c++/v1']
+cpp_link_args = ['-nostdlib++', '-lc++', '-lc++abi']
diff --git a/native/gcc.ini b/native/gcc.ini
new file mode 100644
index 00000000..3b70f11b
--- /dev/null
+++ b/native/gcc.ini
@@ -0,0 +1,3 @@
+[binaries]
+c = 'gcc'
+cpp = 'g++'
diff --git a/native/ld.ini b/native/ld.ini
new file mode 100644
index 00000000..b0fde818
--- /dev/null
+++ b/native/ld.ini
@@ -0,0 +1,3 @@
+[binaries]
+c_ld = 'bfd'
+cpp_ld = 'bfd'
diff --git a/native/libcxx.ini b/native/libcxx.ini
new file mode 100644
index 00000000..82fe10d6
--- /dev/null
+++ b/native/libcxx.ini
@@ -0,0 +1,4 @@
+# Layer after clang.ini. GCC does not support Clang's -stdlib option.
+[built-in options]
+cpp_args = ['-stdlib=libc++']
+cpp_link_args = ['-stdlib=libc++']
diff --git a/native/libstdcxx.ini b/native/libstdcxx.ini
new file mode 100644
index 00000000..ae18519f
--- /dev/null
+++ b/native/libstdcxx.ini
@@ -0,0 +1,5 @@
+# Layer after clang.ini to explicitly select libstdc++.
+# GCC already uses libstdc++ by default and does not accept -stdlib.
+[built-in options]
+cpp_args = ['-stdlib=libstdc++']
+cpp_link_args = ['-stdlib=libstdc++']
diff --git a/native/lld.ini b/native/lld.ini
new file mode 100644
index 00000000..00a3e0a3
--- /dev/null
+++ b/native/lld.ini
@@ -0,0 +1,3 @@
+[binaries]
+c_ld = 'lld'
+cpp_ld = 'lld'
diff --git a/src/compile_time_checks/serif/discretization/blocks/fields.cpp b/src/compile_time_checks/serif/discretization/blocks/fields.cpp
index 2e76e910..d00804ea 100644
--- a/src/compile_time_checks/serif/discretization/blocks/fields.cpp
+++ b/src/compile_time_checks/serif/discretization/blocks/fields.cpp
@@ -24,5 +24,5 @@ namespace {
     static_assert(DensitySpace::order == 2, "Density should be in a second order space");
     static_assert(DensitySpace::rank::value == 0, "Density should be a scalar field");
     static_assert(std::same_as<blocks::SupportOf<Density>, domain::StellarDomains>, "Density should be supported on stellar domains");
-    static_assert(std::same_as<blocks::QuantityOf<Density>, serif::dimensions::Density>, "Density should have a quantity of type Density");
+    static_assert(std::same_as<blocks::QuantityOf<Density>, serif::dimensions::MassDensity>, "Density should have a quantity of type Density");
 }
diff --git a/src/include/serif/dimensions/dimensions.hpp b/src/include/serif/dimensions/dimensions.hpp
index 22746ee7..8f683e0d 100644
--- a/src/include/serif/dimensions/dimensions.hpp
+++ b/src/include/serif/dimensions/dimensions.hpp
@@ -1,154 +1,346 @@
 #pragma once
+/**
+ * @file dimensions.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief File containing prue dimensions. Note that this is not intended to be (nor does it
+ * have the capability to be) a general dimensional analysis library. There is no
+ * understanding here that, for example, a dimension of length^2 is an area. Rather,
+ * the purpose of this is to force developers to be explicit about what values they
+ * are passing to what functions. Further, most of this system is built to be compile
+ * time constant. This means that the compiler can still apply many optimizations. Note
+ * that the dimensions here do not contain a value, rather these are *just* the concept
+ * of a dimension. To store a value you must use a @ref Quantity and @ref QuantityValue.
+ */
 
 #include <string_view>
 
 #include "serif/dimensions/quantities.hpp"
 
 namespace serif::dimensions {
+    /**
+     * @struct Dimensionless
+     *
+     * @brief Non dimensional type.
+     */
     struct Dimensionless final : PhysicalQuantity {
         static constexpr std::string_view identifier = "dimensionless";
     };
 
+    /**
+     * @struct Mass
+     *
+     * @brief Dimension of mass
+     */
     struct Mass final : PhysicalQuantity {
         static constexpr std::string_view identifier = "mass";
     };
 
+    /**
+     * @struct Length
+     *
+     * @brief Dimension of Length
+     */
     struct Length final : PhysicalQuantity {
         static constexpr std::string_view identifier = "length";
     };
 
+    /**
+     * @struct Time
+     *
+     * @brief Dimension of Time
+     */
     struct Time final : PhysicalQuantity {
         static constexpr std::string_view identifier = "time";
     };
 
+    /**
+     * @struct Area
+     *
+     * @brief Dimension of Area
+     */
     struct Area final : PhysicalQuantity {
         static constexpr std::string_view identifier = "area";
     };
 
+    /**
+     * @struct Volume
+     *
+     * @brief Dimension of Volume
+     */
     struct Volume final : PhysicalQuantity {
         static constexpr std::string_view identifier = "volume";
     };
 
-    struct Density final : ThermodynamicQuantity {
+    /**
+     * @struct MassDensity
+     *
+     * @brief Dimension of Mass Density
+     */
+    struct MassDensity final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "density";
     };
 
-    struct SurfaceDensity final : PhysicalQuantity {
-        static constexpr std::string_view identifier = "surface_density";
-    };
-
+    /**
+     * @struct NumberDensity
+     *
+     * @brief Dimension of Number Density
+     */
     struct NumberDensity final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "number_density";
     };
 
+    /**
+     * @struct Pressure
+     *
+     * @brief Dimension of Pressure
+     */
     struct Pressure final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "pressure";
     };
 
+    /**
+     * @struct Temperature
+     *
+     * @brief Dimension of Temperature
+     */
     struct Temperature final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "temperature";
     };
 
+    /**
+     * @struct Entropy
+     *
+     * @brief Dimension of Entropy
+     */
     struct Entropy final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "entropy";
     };
 
+    /**
+     * @struct SpecificEntropy
+     *
+     * @brief Dimension of Specific Entropy
+     */
     struct SpecificEntropy final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "specific_entropy";
     };
 
+    /**
+     * @struct ChemicalPotential
+     *
+     * @brief Dimension of Chemical Potential
+     */
     struct ChemicalPotential final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "chemical_potential";
     };
 
+    /**
+     * @struct Energy
+     *
+     * @brief Dimension of Energy
+     */
     struct Energy final : PhysicalQuantity {
         static constexpr std::string_view identifier = "energy";
     };
 
+    /**
+     * @struct InternalEnergy
+     *
+     * @brief Dimension of Internal Energy
+     */
     struct InternalEnergy final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "internal_energy";
     };
 
+    /**
+     * @struct SpecificEnergy
+     *
+     * @brief Dimension of Specific Energy
+     */
     struct SpecificEnergy final : PhysicalQuantity {
         static constexpr std::string_view identifier = "specific_energy";
     };
 
+    /**
+     * @struct SpecificInternalEnergy
+     *
+     * @brief Dimension of Specific Internal Energy
+     */
     struct SpecificInternalEnergy final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "specific_internal_energy";
     };
 
+    /**
+     * @struct SpecificEnthalpy
+     *
+     * @brief Dimension of Specific Enthalpy
+     */
     struct SpecificEnthalpy final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "specific_enthalpy";
     };
 
+    /**
+     * @struct EnergyDensity
+     *
+     * @brief Dimension of Energy Density
+     */
     struct EnergyDensity final : ThermodynamicQuantity {
         static constexpr std::string_view identifier = "energy_density";
     };
 
+    /**
+     * @struct GravitationalPotential
+     *
+     * @brief Dimension of Gravitational Potential
+     */
     struct GravitationalPotential final : PhysicalQuantity {
         static constexpr std::string_view identifier = "gravitational_potential";
     };
 
+    /**
+     * @struct Velocity
+     *
+     * @brief Dimension of Velocity
+     */
     struct Velocity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "velocity";
     };
 
+    /**
+     * @struct Acceleration
+     *
+     * @brief Dimension of Acceleration
+     */
     struct Acceleration final : PhysicalQuantity {
         static constexpr std::string_view identifier = "acceleration";
     };
 
+    /**
+     * @struct Frequency
+     *
+     * @brief Dimension of Frequency
+     */
     struct Frequency final : PhysicalQuantity {
         static constexpr std::string_view identifier = "frequency";
     };
 
+    /**
+     * @struct AngularVelocity
+     *
+     * @brief Dimension of Angular Velocity
+     */
     struct AngularVelocity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "angular_velocity";
     };
 
+    /**
+     * @struct Momentum
+     *
+     * @brief Dimension of Momentum
+     */
     struct Momentum final : PhysicalQuantity {
         static constexpr std::string_view identifier = "momentum";
     };
 
+    /**
+     * @struct AngularMomentum
+     *
+     * @brief Dimension of Angular Momentum
+     */
     struct AngularMomentum final : PhysicalQuantity {
         static constexpr std::string_view identifier = "angular_momentum";
     };
 
+    /**
+     * @struct MomentOfInertia
+     *
+     * @brief Dimension of Moment Of Inertia
+     */
     struct MomentOfInertia final : PhysicalQuantity {
         static constexpr std::string_view identifier = "moment_of_inertia";
     };
 
+    /**
+     * @struct Force
+     *
+     * @brief Dimension of Force
+     */
     struct Force final : PhysicalQuantity {
         static constexpr std::string_view identifier = "force";
     };
 
+    /**
+     * @struct Torque
+     *
+     * @brief Dimension of Torque
+     */
     struct Torque final : PhysicalQuantity {
         static constexpr std::string_view identifier = "torque";
     };
 
+    /**
+     * @struct Power
+     *
+     * @brief Dimension of Power
+     */
     struct Power final : PhysicalQuantity {
         static constexpr std::string_view identifier = "power";
     };
 
+    /**
+     * @struct Luminosity
+     *
+     * @brief Dimension of Luminosity
+     */
     struct Luminosity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "luminosity";
     };
 
+    /**
+     * @struct MassFlowRate
+     *
+     * @brief Dimension of Mass Flow Rate
+     */
     struct MassFlowRate final : PhysicalQuantity {
         static constexpr std::string_view identifier = "mass_flow_rate";
     };
 
+    /**
+     * @struct Opacity
+     *
+     * @brief Dimension of Opacity
+     * @note we are not currently using this in the code since we are still using a polytropic
+     * eos. However, when we actually get to using this we will likely want to break this out
+     * into more specific and meaningful opacity dimensions (RMO, per unit wavelength, etc...)
+     */
     struct Opacity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "opacity";
     };
 
+    /**
+     * @struct DynamicViscosity
+     *
+     * @brief Dimension of Dynamic Viscosity
+     */
     struct DynamicViscosity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "dynamic_viscosity";
     };
 
+    /**
+     * @struct KinematicViscosity
+     *
+     * @brief Dimension of Kinematic Viscosity
+     */
     struct KinematicViscosity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "kinematic_viscosity";
     };
 
+    /**
+     * @struct MagneticFluxDensity
+     *
+     * @brief Dimension of Magnetic Flux Density
+     */
     struct MagneticFluxDensity final : PhysicalQuantity {
         static constexpr std::string_view identifier = "magnetic_flux_density";
     };
diff --git a/src/include/serif/dimensions/partial.hpp b/src/include/serif/dimensions/partial.hpp
index 52d14c9a..ffff99b6 100644
--- a/src/include/serif/dimensions/partial.hpp
+++ b/src/include/serif/dimensions/partial.hpp
@@ -1,21 +1,78 @@
 #pragma once
 
+/**
+ * @file partial.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief Compile time tools to represent a partial derivative of one quantity wrt. another quantity.
+ * Note that this system is not responsible for computing partials or in fact derivatives of any kind.
+ * Rather, it is responsible for forcing the developer --- or compiler at least --- to ask
+ * for a specific, explicit, partial and have to spell out what that is. That in turn gives
+ * the compiler enough information to do a lot of downstream things. A motivating example would be automatic verification
+ * of the provided EOS against all partials that might be needed to ensure they are all exported
+ * by the EOS.
+ */
+
 #include "serif/dimensions/quantities.hpp"
 #include "serif/utils/misc/concepts/numeric.hpp"
 
 namespace serif::dimensions {
+    /**
+     * @struct WithRespectTo
+     *
+     * @tparam Quantity A quantity satisfying @ref ThermodynamicQuantityType
+     * @brief A type that is empty at runtime but allows for the developer
+     * to specify to the compiler what quantity the partial is with respect to.
+     */
     template <ThermodynamicQuantityType Quantity>
     struct WithRespectTo final {};
 
+    /**
+     * @class PartialDerivative
+     *
+     * @tparam OutputQuantity
+     * @tparam InputQuantity
+     *
+     * @brief Class to represent a compile time derivative relation between two quantities. Note that
+     * this class does not, in fact cannot, compute any derivatives itself. Rather, this is a way
+     * for the compiler to understand higher level mathematical relations. Unless you are really
+     * in the weeds of the development team you will likely not need to interact with this
+     * directly. Further note that this class does include (unlike the more basic dimension system)
+     * some basic dimensional analysis. Specifically, there are some overloaded operators that
+     * allow for certain specializations of PartialDerivative to operate on one and other
+     * conditioned on what their input and output types are. This does not do any thing like
+     * convert area to length^2. However, it will do more basic things like ensure that you
+     * can only add two partials which agree on their dimensions exactly. Or that while you can
+     * multiple a partial by any scalar, you cannot multiply arbitrary partials.
+     */
     template <ThermodynamicQuantityType OutputQuantity, ThermodynamicQuantityType InputQuantity>
     class PartialDerivative final {
     public:
         explicit constexpr PartialDerivative(const double value) noexcept : m_value(value){}
 
+        /**
+         * @brief Get the underlying value as a compile time constant if possible. This
+         * allows the compiler, where the constructor may be determined at compile time,
+         * to fully optimize PartialDerivative operations as if they were doubles.
+         *
+         * @note calling value gets the underlying numeric value directly thereby bypassing
+         * the limitations on what a valid operation are (e.g. what the operator overloads
+         * enforce).
+         *
+         * @return m_value: This is whatever was passed at constructor time since there is
+         * no way to modify this value after construction without some truly horrific C++ tricks
+         */
         [[nodiscard]] constexpr double value() const noexcept {
             return m_value;
         }
 
+        /**
+         * @brief Overload the + operator between two partials which share the exact same input
+         * and output types
+         * @param lhs One partial
+         * @param rhs Other Partial
+         * @return A new partial with the same input and output types where m_value = lhs.m_value + rhs.m_value
+         */
         friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator+(
             const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
             const PartialDerivative<OutputQuantity, InputQuantity>& rhs
@@ -23,6 +80,13 @@ namespace serif::dimensions {
             return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value + rhs.m_value};
         }
 
+        /**
+         * @brief Overload the - operator between two partials which share the exact same input
+         * and output types
+         * @param lhs One partial
+         * @param rhs Other Partial
+         * @return A new partial with the same input and output types where m_value = lhs.m_value - rhs.m_value
+         */
         friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator-(
             const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
             const PartialDerivative<OutputQuantity, InputQuantity>& rhs
@@ -30,6 +94,14 @@ namespace serif::dimensions {
             return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value - rhs.m_value};
         }
 
+        /**
+         * @brief Overload of the * operator for a left partial and a right scalar
+         * and output types
+         * @tparam Scalar Any value satisfying the @ref IsNumeric concept (integral types and float types)
+         * @param lhs Partial with some input and some output type
+         * @param rhs Some numeric quantity
+         * @return A new partial with the same input and output types where m_value = lhs.m_value * rhs
+         */
         template <utils::misc::concepts::IsNumeric Scalar>
         friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator*(
             const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
@@ -38,6 +110,14 @@ namespace serif::dimensions {
             return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value * rhs};
         }
 
+        /**
+         * @brief Overload of the * operator for a right partial and a left scalar
+         * and output types
+         * @tparam Scalar Any value satisfying the @ref IsNumeric concept (integral types and float types)
+         * @param lhs Some numeric quantity
+         * @param rhs Partial with some input and some output type
+         * @return A new partial with the same input and output types where m_value = lhs * rhs.m_value
+         */
         template <utils::misc::concepts::IsNumeric Scalar>
         friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator*(
             const Scalar& lhs,
@@ -46,6 +126,14 @@ namespace serif::dimensions {
             return PartialDerivative<OutputQuantity, InputQuantity>{lhs * rhs.m_value};
         }
 
+        /**
+         * @brief Overload of the / operator for a left partial and a right scalar
+         * and output types
+         * @tparam Scalar Any value satisfying the @ref IsNumeric concept (integral types and float types)
+         * @param lhs Partial with some input and some output type
+         * @param rhs Some numeric quantity
+         * @return A new partial with the same input and output types where m_value = lhs.m_value / rhs
+         */
         template <utils::misc::concepts::IsNumeric Scalar>
         friend constexpr PartialDerivative<OutputQuantity, InputQuantity> operator/(
             const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
@@ -54,6 +142,13 @@ namespace serif::dimensions {
             return PartialDerivative<OutputQuantity, InputQuantity>{lhs.m_value / rhs};
         }
 
+        /**
+         * @brief Partial ordering overload to compare a left partial and a right scalar
+         * @tparam Scalar Any value satisfying the @ref IsNumeric concept (integral types and float types)
+         * @param lhs Some numeric quantity
+         * @param rhs Some partial with an input and output type
+         * @return partial ordering comparing rhs.m_value to lhs
+         */
         template <utils::misc::concepts::IsNumeric Scalar>
         friend constexpr std::partial_ordering operator<=>(
             const PartialDerivative<OutputQuantity, InputQuantity>& lhs,
@@ -62,6 +157,13 @@ namespace serif::dimensions {
             return lhs.m_value <=> rhs;
         }
 
+        /**
+         * @brief Partial ordering overload to compare a right partial and a left scalar
+         * @tparam Scalar Any value satisfying the @ref IsNumeric concept (integral types and float types)
+         * @param lhs Some numeric quantity
+         * @param rhs Some partial with an input and output type
+         * @return partial ordering comparing rhs.m_value to lhs
+         */
         template <utils::misc::concepts::IsNumeric Scalar>
         friend constexpr std::partial_ordering operator<=>(
             Scalar lhs,
@@ -70,6 +172,12 @@ namespace serif::dimensions {
             return lhs <=> rhs.m_value;
         }
 
+        /**
+         * @brief Partial ordering comparing two partials with exactly the same input and output types
+         * @param lhs Some Partial
+         * @param rhs Other Partial
+         * @return Partial ordering comparing lhs.m_value to rhs.m_value
+         */
         friend constexpr std::partial_ordering operator<=>(
             const PartialDerivative<OutputQuantity, InputQuantity> &lhs,
             const PartialDerivative<OutputQuantity, InputQuantity> &rhs
@@ -78,6 +186,6 @@ namespace serif::dimensions {
         }
 
     private:
-        double m_value;
+        double m_value; ///< Numeric value associated with the partial. If this is known at compile time then the compiler can optimize operations with partial.value() to be equivalent to compile time double operations.
     };
 }
\ No newline at end of file
diff --git a/src/include/serif/dimensions/quantities.hpp b/src/include/serif/dimensions/quantities.hpp
index a186cdb0..b8a8110b 100644
--- a/src/include/serif/dimensions/quantities.hpp
+++ b/src/include/serif/dimensions/quantities.hpp
@@ -1,4 +1,11 @@
 #pragma once
+/**
+ * @file quantities.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief This file contains the code to associate numeric data to a dimension. We call
+ * this a Quantity and a QuantityValue.
+ */
 
 #include <type_traits>
 #include <compare>
@@ -14,19 +21,60 @@ namespace serif::dimensions {
 
     using utils::misc::concepts::IsNumeric;
 
+    /**
+     * @struct PhysicalQuantity
+     *
+     * @breif General type for any physical quantity.
+     */
     struct PhysicalQuantity { };
+
+    /**
+     * @struct ThermodynamicQuantity
+     *
+     * @breif Slightly more specific type for thermodynamic types. Note both of these
+     * are empty at runtime, these are only useful as compile time tools to control
+     * what code the compiler generates.
+     */
     struct ThermodynamicQuantity: PhysicalQuantity { };
 
+    /**
+     * @concept PhysicalQuantityType
+     *
+     * @tparam QuantityCandidate The type that will be validated
+     *
+     * @brief Is compile-time true of the QuantityCandidate can have any const / volatile stripped
+     * and if its derived from PhysicalQuantity. Note all ThermodynamicQuantities will also
+     * satisfy this concept.
+     */
     template <typename QuantityCandidate>
     concept PhysicalQuantityType =
         std::same_as<QuantityCandidate, NeitherConstNorVolatile<QuantityCandidate>> &&
         std::derived_from<QuantityCandidate, PhysicalQuantity>;
 
+    /**
+     * @concept ThermodynamicQuantityType
+     *
+     * @tparam QuantityCandidate
+     *
+     * @brief Is compile-time true of the QuantityCandidate can have any const / volatile stripped
+     * and if its derived from PhysicalQuantity (both enforced through checking @ref PhysicalQuantityType)
+     * and then that is is also derived from ThermodynamicQuantity. Note not all physical quanties will
+     * satisfy this.
+     */
     template <typename QuantityCandidate>
     concept ThermodynamicQuantityType =
         PhysicalQuantityType<QuantityCandidate> &&
         std::derived_from<QuantityCandidate, ThermodynamicQuantity>;
 
+    /**
+     * @class QuantityValue
+     * @tparam Quantity Some valid PhysicalQuantityType. Note all the dimensions in @ref dimensions.hpp
+     * are of type PhysicalQuantity. There should be no other dimensions in the code.
+     *
+     * @brief This is the non derivative equivalent to @ref PartialDerivative. This associates
+     * a dimension with a value and, if that value can be determined at compile time, can be
+     * decomposed to a transparent double through .value()
+     */
     template <PhysicalQuantityType Quantity>
     class QuantityValue final {
     public:
@@ -36,11 +84,25 @@ namespace serif::dimensions {
             return m_value;
         }
 
+
+        /**
+         * @breif Default comparison operator for two QuantityValues. Note that there is
+         * an implicit requirement these be of the same specialization. E.g. each have the same
+         * PhysicalQuantityType Quantity type parameter. Because this is the default and this
+         * class has only one trivially comparable member this just checks if m_value == m_value.
+         * @return true if the values in both are equal.
+         */
         [[nodiscard]] friend constexpr bool operator==(
             const QuantityValue &,
             const QuantityValue &
         ) noexcept = default;
 
+        /**
+         * @breif overload of + operator for two quantities with the same PhysicalQuantityType
+         * @param lhs One Quantity
+         * @param rhs Other Quantity
+         * @return a new quantity with the same PhysicalQuantityType but with a value = lhs.m_value + rhs.m_value
+         */
         friend constexpr QuantityValue operator+(
             const QuantityValue &lhs,
             const QuantityValue &rhs
@@ -48,6 +110,12 @@ namespace serif::dimensions {
             return QuantityValue{lhs.m_value + rhs.m_value};
         }
 
+        /**
+         * @breif overload of - operator for two quantities with the same PhysicalQuantityType
+         * @param lhs One Quantity
+         * @param rhs Other Quantity
+         * @return a new quantity with the same PhysicalQuantityType but with a value = lhs.m_value - rhs.m_value
+         */
         friend constexpr QuantityValue operator-(
             const QuantityValue &lhs,
             const QuantityValue &rhs
@@ -55,6 +123,13 @@ namespace serif::dimensions {
             return QuantityValue{lhs.m_value - rhs.m_value};
         }
 
+        /**
+         * @breif Overload of the * operator for a left quantity and some right numeric, scalar, value.
+         * @tparam Scalar A value satisfying @ref IsNumeric
+         * @param lhs some quantity value
+         * @param rhs some scalar value
+         * @return some quantity value with the same type as lhs but with value = lhs.m_value * scalar
+         */
         template <IsNumeric Scalar>
         friend constexpr QuantityValue operator*(
             const QuantityValue &lhs,
@@ -63,6 +138,13 @@ namespace serif::dimensions {
             return QuantityValue{lhs.m_value * static_cast<double>(rhs)};
         }
 
+        /**
+         * @breif Overload of the * operator for a right quantity and some left numeric, scalar, value.
+         * @tparam Scalar A value satisfying @ref IsNumeric
+         * @param lhs some scalar value
+         * @param rhs some quantity value
+         * @return some quantity value with the same type as lhs but with value = scalar * lhs.m_value
+         */
         template <IsNumeric Scalar>
         friend constexpr QuantityValue operator*(
             const Scalar lhs,
@@ -71,6 +153,13 @@ namespace serif::dimensions {
             return QuantityValue{static_cast<double>(lhs) * rhs.m_value};
         }
 
+        /**
+         * @breif Overload of the / operator for a left quantity and some right numeric, scalar, value.
+         * @tparam Scalar A value satisfying @ref IsNumeric
+         * @param lhs some quantity value
+         * @param rhs some scalar value
+         * @return some quantity value with the same type as lhs but with value = lhs.m_value / scalar
+         */
         template <IsNumeric Scalar>
         friend constexpr QuantityValue operator/(
             const QuantityValue &lhs,
@@ -79,6 +168,13 @@ namespace serif::dimensions {
             return QuantityValue{lhs.m_value / static_cast<double>(rhs)};
         }
 
+        /**
+         * @breif Partial ordering overload for some left quantity and right scalar
+         * @tparam Scalar A value satisfying @ref IsNumeric
+         * @param lhs some quantity value
+         * @param rhs some scalar value
+         * @return The partial ordering comparing lhs.m_value to scalar
+         */
         template <IsNumeric Scalar>
         friend constexpr std::partial_ordering operator<=>(
             const QuantityValue &lhs,
@@ -87,6 +183,13 @@ namespace serif::dimensions {
             return lhs.m_value <=> static_cast<double>(rhs);
         }
 
+        /**
+         * @breif Partial ordering overload for some right quantity and left scalar
+         * @tparam Scalar A value satisfying @ref IsNumeric
+         * @param lhs some scalar value
+         * @param rhs some quantity value
+         * @return The partial ordering comparing scalar to lhs.m_value
+         */
         template <IsNumeric Scalar>
         friend constexpr std::partial_ordering operator<=>(
             const Scalar lhs,
@@ -95,6 +198,12 @@ namespace serif::dimensions {
             return static_cast<double>(lhs) <=> rhs.m_value;
         }
 
+        /**
+         * @breif partial ordering overload to compare two quantity values of the same type
+         * @param lhs Some quantity value
+         * @param rhs Other quantity value
+         * @return the partial ordering of lhs.m_value <=> rhs.m_value
+         */
         friend constexpr std::partial_ordering operator<=>(
             const QuantityValue &lhs,
             const QuantityValue &rhs
@@ -103,19 +212,43 @@ namespace serif::dimensions {
         }
 
     private:
-        double m_value;
+        double m_value; ///< Numeric value associated with the dimension of QuantityValue.
     };
 
+    /**
+     * @struct IsQuantityValue
+     * @tparam Candidate some type that may be a QuantityValue
+     *
+     * @breif Most general specialization, defaults to compile time false.
+     */
     template <typename Candidate>
     struct IsQuantityValue : std::false_type { };
 
+    /**
+     * @struct IsQuantityValue
+     * @tparam Quantity Some quantity that satisfies the @ref PhysicalQuantityType concept
+     *
+     * @breif More specific so if a type satisfies this the compiler will select it first. This evaluates
+     * to compile time true if Quantity satisfies the PhysicalQuantityType concept which then
+     * creats a specialization of IsQuantityValue that is compile time true for any Quantity of that set.
+     */
     template <PhysicalQuantityType Quantity>
     struct IsQuantityValue<QuantityValue<Quantity>> : std::true_type { };
 
+    /**
+     * @concept QuantityValueType
+     * @tparam Candidate Some type that might be a quantity value type
+     *
+     * @brief Turn those compile time structs into a concept that can be used to constrain
+     * a template parameter
+     */
     template <typename Candidate>
     concept QuantityValueType = IsQuantityValue<std::remove_cvref_t<Candidate>>::value;
 
-
+    /**
+     * @struct QuantityOf
+     * @tparam Candidate Some type that may be able to report what its dimension is
+     */
     template <typename Candidate>
     struct QuantityOf;
 
diff --git a/src/include/serif/dimensions/runtime/concepts.hpp b/src/include/serif/dimensions/runtime/concepts.hpp
index 5b7a4420..03d7a6ac 100644
--- a/src/include/serif/dimensions/runtime/concepts.hpp
+++ b/src/include/serif/dimensions/runtime/concepts.hpp
@@ -1,4 +1,13 @@
 #pragma once
+/**
+ * @file concepts.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief Runtime concept for a type that can be projected to a runtime type.
+ * This is useful when converting from the compile time type system to something
+ * that can be inspected at runtime.
+ */
+
 #include <string_view>
 #include <type_traits>
 
@@ -6,12 +15,25 @@
 #include "serif/dimensions/runtime/runtime.hpp"
 
 namespace serif::dimensions::runtime {
+
+    /**
+     * @concept RuntimeIdentifiedThermodynamicQuantity
+     * @breif Concept to limit template arguments such that they are compatible with the expected
+     * ThermodynamicQuantity interface. Specifically, a valid type must contain a member
+     * identifier which can be converted to std::string_view and that member, when converted
+     * to a string_view must be longer than 0 length.
+     */
     template <typename Quantity>
     concept RuntimeIdentifiedThermodynamicQuantity =
         ThermodynamicQuantityType<Quantity> &&
         requires {{Quantity::identifier} -> std::convertible_to<std::string_view>; } &&
         (std::string_view{Quantity::identifier}.size() > 0);
 
+    /**
+     * @brief Extract a compile time identifier from a thermodynamic quantity and project it into a runtime
+     * string_view
+     * @tparam Quantity Type which @ref RuntimeIdentifiedThermodynamicQuantity allows
+     */
     template <RuntimeIdentifiedThermodynamicQuantity Quantity>
     inline constexpr ThermodynamicQuantityID thermodynamicQuantityID{std::string_view{Quantity::identifier}};
 }
\ No newline at end of file
diff --git a/src/include/serif/dimensions/runtime/runtime.hpp b/src/include/serif/dimensions/runtime/runtime.hpp
index 45015b42..263a6bf9 100644
--- a/src/include/serif/dimensions/runtime/runtime.hpp
+++ b/src/include/serif/dimensions/runtime/runtime.hpp
@@ -1,13 +1,30 @@
 #pragma once
+/**
+ * @file runtime.hpp
+ * @author Emily M. Boudreaux
+ *
+ * @brief runtime usable types representing the compile time dimension type system.
+ * Note that generally this is not a dimensional analysis library. Though perhaps
+ * in the future it will grow to something like that. Rather, this system
+ * is intended to ensure that developers are explicit which what numeric values mean.
+ * There is more detail on this in the file comment for @ref dimensions.hpp and @ref quantities.hpp
+ */
+
 #include <string_view>
 #include <span>
 #include <cstdint>
 
 
 namespace serif::dimensions::runtime {
+    /**
+     * @class ThermodynamicQuantityID
+     *
+     * @brief Runtime resolvable ID for a thermodynamic quantity. Used in conjunction with
+     * concepts to resolve compile time thermodynamic types to runtime values.
+     */
     class ThermodynamicQuantityID final {
     public:
-        explicit constexpr ThermodynamicQuantityID(std::string_view name) noexcept : m_name(name) {};
+        explicit constexpr ThermodynamicQuantityID(const std::string_view name) noexcept : m_name(name) {};
 
         [[nodiscard]] constexpr std::string_view name() const noexcept {
             return m_name;
@@ -18,9 +35,17 @@ namespace serif::dimensions::runtime {
         std::string_view m_name; // We can use a string view here since the constructor is constexpr and the string view will be valid for the lifetime of the program.
     };
 
+    /**
+     * @struct RuntimeQuantityValue
+     *
+     * @brief Runtime resolved type for a quantity value. This holds both an ID and a value. This is the runtime
+     * equivalent to a @ref QuantityValue. Note that where these differ is that QuantityValue
+     * is built with template types. This allows the compiler to resolve it to a plain numeric type and still do optimizations
+     * That is generally not the case with RuntimeQuantityValue
+     */
     struct RuntimeQuantityValue final {
         ThermodynamicQuantityID id;
-        double value;
+        double value{};
     };
 
     // TODO: This should be refactored into the eos module rather than the dimensions module
diff --git a/src/include/serif/dimensions/type_alias.hpp b/src/include/serif/dimensions/type_alias.hpp
index 56bb5328..473f294b 100644
--- a/src/include/serif/dimensions/type_alias.hpp
+++ b/src/include/serif/dimensions/type_alias.hpp
@@ -10,8 +10,7 @@ namespace serif::dimensions {
     using TimeValue                                      = QuantityValue<Time>;
     using AreaValue                                      = QuantityValue<Area>;
     using VolumeValue                                    = QuantityValue<Volume>;
-    using DensityValue                                   = QuantityValue<Density>;
-    using SurfaceDensityValue                            = QuantityValue<SurfaceDensity>;
+    using DensityValue                                   = QuantityValue<MassDensity>;
     using NumberDensityValue                             = QuantityValue<NumberDensity>;
     using PressureValue                                  = QuantityValue<Pressure>;
     using TemperatureValue                               = QuantityValue<Temperature>;
diff --git a/src/include/serif/discretization/blocks/fields.hpp b/src/include/serif/discretization/blocks/fields.hpp
index 0059e632..6aefa34e 100644
--- a/src/include/serif/discretization/blocks/fields.hpp
+++ b/src/include/serif/discretization/blocks/fields.hpp
@@ -27,7 +27,7 @@ namespace serif::discretization::blocks {
 
             using FunctionalSpace = functional_space::DiscreteFunctionalSpace<functional_space::L2, 2, 0>; // second order scalar L2 space
             using Support         = domain::StellarDomains;
-            using Quantity        = dimensions::Density;
+            using Quantity        = dimensions::MassDensity;
             static constexpr std::string_view symbol = "ρ";
 
             struct value final : value_block_base {};
diff --git a/src/include/serif/eos/models/polytropic.hpp b/src/include/serif/eos/models/polytropic.hpp
index 3c88f269..5ee5f1fa 100644
--- a/src/include/serif/eos/models/polytropic.hpp
+++ b/src/include/serif/eos/models/polytropic.hpp
@@ -39,7 +39,7 @@ namespace serif::eos::models {
         [[nodiscard]] dimensions::SpecificEnthalpyValue evaluate(relations::SpecificEnthalpyFromPressure, dimensions::PressureValue pressure) const;
         [[nodiscard]] dimensions::DensityValue evaluate(relations::DensityFromSpecificEnthalpy, dimensions::SpecificEnthalpyValue specific_enthalpy) const;
 
-        [[nodiscard]] dimensions::PartialDerivative<dimensions::Density, dimensions::SpecificEnthalpy> partial_derivative(
+        [[nodiscard]] dimensions::PartialDerivative<dimensions::MassDensity, dimensions::SpecificEnthalpy> partial_derivative(
             relations::DensityFromSpecificEnthalpy,
             dimensions::WithRespectTo<dimensions::SpecificEnthalpy>,
             dimensions::SpecificEnthalpyValue specific_enthalpy
@@ -51,9 +51,9 @@ namespace serif::eos::models {
             dimensions::SpecificEnthalpyValue specific_enthalpy
         ) const;
 
-        [[nodiscard]] dimensions::PartialDerivative<dimensions::Pressure, dimensions::Density> partial_derivative(
+        [[nodiscard]] dimensions::PartialDerivative<dimensions::Pressure, dimensions::MassDensity> partial_derivative(
             relations::PressureFromDensity,
-            dimensions::WithRespectTo<dimensions::Density>,
+            dimensions::WithRespectTo<dimensions::MassDensity>,
             dimensions::DensityValue density
         ) const;
 
diff --git a/src/include/serif/eos/relations/relations.hpp b/src/include/serif/eos/relations/relations.hpp
index c0e8bfc6..d6a54ad6 100644
--- a/src/include/serif/eos/relations/relations.hpp
+++ b/src/include/serif/eos/relations/relations.hpp
@@ -6,7 +6,7 @@
 namespace serif::eos::relations {
     using PressureFromDensity = Relation<
         dimensions::Pressure,
-        dimensions::Density>;
+        dimensions::MassDensity>;
 
     using PressureFromSpecificEnthalpy = Relation<
         dimensions::Pressure,
@@ -14,13 +14,13 @@ namespace serif::eos::relations {
 
     using SpecificEnthalpyFromDensity = Relation<
         dimensions::SpecificEnthalpy,
-        dimensions::Density>;
+        dimensions::MassDensity>;
 
     using SpecificEnthalpyFromPressure = Relation<
         dimensions::SpecificEnthalpy,
         dimensions::Pressure>;
 
     using DensityFromSpecificEnthalpy = Relation<
-        dimensions::Density,
+        dimensions::MassDensity,
         dimensions::SpecificEnthalpy>;
 }
\ No newline at end of file
diff --git a/src/include/serif/mapping/compactification/compactification.hpp b/src/include/serif/mapping/compactification/compactification.hpp
new file mode 100644
index 00000000..d960d6bd
--- /dev/null
+++ b/src/include/serif/mapping/compactification/compactification.hpp
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <string_view>
+
+#include "serif/mapping/status.hpp"
+#include "serif/mapping/compactification/input_result.hpp"
+#include "serif/mapping/compactification/variation.hpp"
+
+
+namespace serif::mapping::compactification {
+    class ExteriorDomainMap {
+    public:
+        virtual ~ExteriorDomainMap() = default;
+
+        [[nodiscard]] virtual bool is_affine_in_displacement() const noexcept;
+
+        [[nodiscard]] virtual MappingStatus evaluate(
+            const ExteriorMapInput &input,
+            ExteriorMapResult &result
+        ) const = 0;
+
+        [[nodiscard]] virtual MappingStatus evaluate_variation(
+            const ExteriorMapInput &input,
+            const ExteriorMapResult &result,
+            const ExteriorMapDirection &direction,
+            ExteriorMapVariation &variation
+        ) const = 0;
+
+        [[nodiscard]] virtual std::string_view get_name() const noexcept = 0;
+
+    };
+}
\ No newline at end of file
diff --git a/src/include/serif/mapping/compactification/input_result.hpp b/src/include/serif/mapping/compactification/input_result.hpp
new file mode 100644
index 00000000..792e7fff
--- /dev/null
+++ b/src/include/serif/mapping/compactification/input_result.hpp
@@ -0,0 +1,19 @@
+#pragma once
+
+#include "mfem/linalg/vector.hpp"
+#include "mfem/linalg/densemat.hpp"
+
+namespace serif::mapping::compactification {
+    struct ExteriorMapInput {
+        const mfem::Vector &reference_position;
+        const mfem::Vector &displaced_position;
+        const mfem::DenseMatrix &displacement_jacobian;
+        const mfem::Vector &compactification_coordinate_gradient;
+        double compactification_coordinate;
+    };
+
+    struct ExteriorMapResult {
+        mfem::Vector physical_position;
+        mfem::DenseMatrix mapping_jacobian;
+    };
+}
diff --git a/src/include/serif/mapping/compactification/kelvin.hpp b/src/include/serif/mapping/compactification/kelvin.hpp
new file mode 100644
index 00000000..c8a1689c
--- /dev/null
+++ b/src/include/serif/mapping/compactification/kelvin.hpp
@@ -0,0 +1,62 @@
+#pragma once
+
+#include <string_view>
+
+#include "mfem/linalg/vector.hpp"
+
+#include "serif/mapping/status.hpp"
+#include "serif/mapping/compactification/compactification.hpp"
+#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
+
+namespace serif::mapping::compactification {
+    class KelvinCompactification final : public ExteriorDomainMap {
+    public:
+        explicit KelvinCompactification(const utils::user_options::KelvinCompactificationOptions& options);
+
+        [[nodiscard]] bool is_affine_in_displacement() const noexcept override;
+        [[nodiscard]] MappingStatus evaluate(
+            const ExteriorMapInput& input,
+            ExteriorMapResult& result
+        ) const override;
+
+        [[nodiscard]] MappingStatus evaluate_variation(
+            const ExteriorMapInput& input,
+            const ExteriorMapResult& result,
+            const ExteriorMapDirection& direction,
+            ExteriorMapVariation& variation
+        ) const override;
+
+        [[nodiscard]] std::string_view get_name() const noexcept override;
+
+        [[nodiscard]] double get_reference_stellar_radius() const noexcept;
+        [[nodiscard]] double get_reference_infinity_radius() const noexcept;
+        [[nodiscard]] double get_coordinate_tolerance() const noexcept;
+
+        [[nodiscard]] MappingStatus evaluate_jacobian_variation_factors(
+            double compactification_coordinate,
+            const mfem::Vector &coordinateGradient,
+            double &scale,
+            mfem::Vector &scaleGradient
+        ) const;
+
+    private:
+        struct RadialFactors {
+            double compactification_coordinate;
+            double computational_radius;
+            double scale;
+            double scale_derivative;
+        };
+
+    private:
+        [[nodiscard]] MappingStatus compute_radial_factors(
+            double compactification_coordinate,
+            RadialFactors &radial_factors
+        ) const;
+
+    private:
+        utils::user_options::KelvinCompactificationOptions m_options;
+        const std::string_view m_name = "KelvinCompactification";
+    };
+
+
+}
\ No newline at end of file
diff --git a/src/include/serif/mapping/compactification/variation.hpp b/src/include/serif/mapping/compactification/variation.hpp
new file mode 100644
index 00000000..53e1f1ce
--- /dev/null
+++ b/src/include/serif/mapping/compactification/variation.hpp
@@ -0,0 +1,16 @@
+#pragma once
+
+#include "mfem/linalg/vector.hpp"
+#include "mfem/linalg/densemat.hpp"
+
+namespace serif::mapping::compactification {
+    struct ExteriorMapDirection {
+        const mfem::Vector &displaced_position_variation;
+        const mfem::DenseMatrix &displacement_jacobian_variation;
+    };
+
+    struct ExteriorMapVariation {
+        mfem::Vector physical_position_variation;
+        mfem::DenseMatrix mapping_jacobian_variation;
+    };
+}
diff --git a/src/include/serif/mapping/context.hpp b/src/include/serif/mapping/context.hpp
new file mode 100644
index 00000000..faea57cd
--- /dev/null
+++ b/src/include/serif/mapping/context.hpp
@@ -0,0 +1,42 @@
+#pragma once
+
+#include "mfem/linalg/densemat.hpp"
+#include "mfem/linalg/vector.hpp"
+
+namespace serif::mapping {
+    struct VolumeQuadratureContext {
+        mfem::DenseMatrix J_inv;
+        double detJ;
+        double weight;
+    };
+
+    struct FaceQuadratureContext {
+        mfem::Vector normal;
+        double ds;
+        double v_dot_n_scale;
+    };
+
+    struct MappingPointContext {
+        mfem::Vector reference_position;
+        mfem::Vector displaced_position;
+        mfem::Vector physical_position;
+        mfem::DenseMatrix displacement_jacobian;
+        mfem::DenseMatrix mapping_jacobian;
+        mfem::DenseMatrix inverse_mapping_jacobian;
+        double mapping_determinate{0.0};
+        bool compactified{false};
+    };
+
+    struct VolumeMappingContext {
+        MappingPointContext mapping;
+        VolumeQuadratureContext quadrature;
+    };
+
+    struct FaceMappingContext {
+        MappingPointContext mapping;
+        FaceQuadratureContext quadrature;
+        mfem::Vector reference_normal;
+        double reference_surface_weight{0.0};
+        double physical_surface_weight{0.0};
+    };
+}
\ No newline at end of file
diff --git a/src/include/serif/mapping/status.hpp b/src/include/serif/mapping/status.hpp
new file mode 100644
index 00000000..ea123b53
--- /dev/null
+++ b/src/include/serif/mapping/status.hpp
@@ -0,0 +1,55 @@
+#pragma once
+
+#include <unordered_map>
+#include <cstdint>
+#include <utility>
+
+#include "serif/utils/error/codes.hpp"
+#include "serif/meta/enums.hpp"
+
+namespace serif::mapping
+{
+    enum class MappingStatus : uint16_t {
+        VALID,
+
+        // Start this at 400 so this is parallel to the mapping error codes. There must be a better way to
+        // do this but it is not worth the time to sort out now
+        INVALID_DIMENSION = 400,
+        NON_FINITE_MAPPING_INPUT,
+        INVALID_REFERENCE_RADIUS,
+        AT_COMPACTIFIED_INFINITY,
+        OUTSIDE_REFERENCE_DOMAIN,
+        NON_FINITE_MAPPING_RESULT,
+        NON_POSITIVE_DETERMINATE
+    };
+
+    static auto MappingStatus_enum_string_map = meta::enum_to_enum_string_map<MappingStatus>();
+    static auto MappingStatus_string_enum_map = meta::enum_to_string_enum_map<MappingStatus>();
+
+    inline std::unordered_map<MappingStatus, utils::errors::SERiFErrorCode> make_MappingStatus_SERiFErrorCode_map() {
+        std::unordered_map<MappingStatus, utils::errors::SERiFErrorCode> ret;
+        for (const auto& [MS, MS_name] : MappingStatus_enum_string_map) {
+            if (std::to_underlying(MS) < 400) {
+                continue; // Skip VALID since it is not in the error code set
+            }
+            ret.insert({MS, utils::errors::SERiFErrorCodeNameToTypeMap.at(MS_name)});
+        }
+        return ret;
+    }
+
+    inline std::unordered_map<utils::errors::SERiFErrorCode, MappingStatus> make_SERIFErrorCode_MappingStatus_map() {
+
+        std::unordered_map<utils::errors::SERiFErrorCode, MappingStatus> ret;
+        for (const auto& [MS, MS_name] : MappingStatus_enum_string_map) {
+            if (std::to_underlying(MS) < 400) {
+                continue; // Skip VALID since it is not in the error code set
+            }
+            ret.insert({utils::errors::SERiFErrorCodeNameToTypeMap.at(MS_name), MS});
+        }
+        return ret;
+    }
+
+    static auto MappingStatus_SERiFErrorCode_map = make_MappingStatus_SERiFErrorCode_map();
+    static auto SERiFErrorCode_MappingStatus_map = make_SERIFErrorCode_MappingStatus_map();
+}
+
diff --git a/src/include/serif/mapping/variation.hpp b/src/include/serif/mapping/variation.hpp
new file mode 100644
index 00000000..02fbdcd3
--- /dev/null
+++ b/src/include/serif/mapping/variation.hpp
@@ -0,0 +1,29 @@
+#pragma once
+
+#include "mfem/linalg/densemat.hpp"
+#include "mfem/linalg/vector.hpp"
+
+namespace serif::mapping {
+    struct MappingPointVariation {
+        mfem::Vector displacement_variation;
+        mfem::Vector physical_position_variation;
+        mfem::DenseMatrix displacement_jacobian_variation;
+        mfem::DenseMatrix mapping_jacobian_variation;
+        mfem::DenseMatrix inverse_mapping_jacobian_variation;
+        double mapping_determinate_variation{0.0};
+    };
+
+    struct VolumeMappingVariation {
+        MappingPointVariation mapping;
+        mfem::DenseMatrix inverse_element_jacobian_variation;
+        double weight_variation{0.0};
+    };
+
+    struct FaceMappingVariation {
+        MappingPointVariation mapping;
+        mfem::Vector physical_normal_variation;
+
+        double physical_surface_weight_variation{0.0};
+        double normal_flux_scale_variation{0.0};
+    };
+}
\ No newline at end of file
diff --git a/src/include/serif/utils/checks/finite.hpp b/src/include/serif/utils/checks/finite.hpp
new file mode 100644
index 00000000..a036dde8
--- /dev/null
+++ b/src/include/serif/utils/checks/finite.hpp
@@ -0,0 +1,31 @@
+#pragma once
+
+#include <cmath>
+
+#include "mfem/linalg/vector.hpp"
+#include "mfem/linalg/densemat.hpp"
+
+namespace serif::utils::checks {
+    // If we have a new backend / new vector and matrix type in the future we can add new overloads
+    inline bool is_finite(const mfem::Vector& vector) {
+        for (int i = 0; i < vector.Size(); ++i) {
+            if (!std::isfinite(vector(i))) {
+                return false;
+            }
+        }
+
+        return true;
+    }
+
+    inline bool is_finite(const mfem::DenseMatrix& matrix) {
+        for (int i = 0; i < matrix.Height(); ++i) {
+            for (int j = 0; j < matrix.Width(); ++j) {
+                if (!std::isfinite(matrix(i, j))) {
+                    return false;
+                }
+            }
+        }
+
+        return true;
+    }
+}
\ No newline at end of file
diff --git a/src/include/serif/utils/error/codes.hpp b/src/include/serif/utils/error/codes.hpp
index 510a9549..73fa596e 100644
--- a/src/include/serif/utils/error/codes.hpp
+++ b/src/include/serif/utils/error/codes.hpp
@@ -42,6 +42,18 @@ namespace serif::utils::errors {
         EOS_NONFINITE_RESULT,
         EOS_INVALID_THERMODYNAMIC_INPUT,
 
+        // Mapping Errors
+        INVALID_DIMENSION = 400,
+        NON_FINITE_MAPPING_INPUT,
+        INVALID_REFERENCE_RADIUS,
+        AT_COMPACTIFIED_INFINITY,
+        OUTSIDE_REFERENCE_DOMAIN,
+        NON_FINITE_MAPPING_RESULT,
+        NON_POSITIVE_DETERMINATE,
+
+        IMPROPERLY_ORDERED_REFERENCE_RADII = 450,
+        INVALID_COORDINATE_TOLERANCE,
+
         // COMPILE TIME FAILURES (>1000)
         // Type Invariant Failures
         FAILED_TYPE_INVARIANT = 1000
diff --git a/src/include/serif/utils/types/tensor/backend/mfem/tensor.hpp b/src/include/serif/utils/types/tensor/backend/mfem/tensor.hpp
new file mode 100644
index 00000000..e708031f
--- /dev/null
+++ b/src/include/serif/utils/types/tensor/backend/mfem/tensor.hpp
@@ -0,0 +1,7 @@
+#pragma once
+
+#include "serif/utils/types/tensor/tensor.hpp"
+
+namespace serif::utils::types::backend::mfem {
+
+}
\ No newline at end of file
diff --git a/src/include/serif/utils/types/tensor/tensor.hpp b/src/include/serif/utils/types/tensor/tensor.hpp
new file mode 100644
index 00000000..d39d9e2a
--- /dev/null
+++ b/src/include/serif/utils/types/tensor/tensor.hpp
@@ -0,0 +1,31 @@
+#pragma once
+#include <cstddef>
+
+namespace serif::utils::types {
+    template <typename T>
+    struct VectorBase {
+        using Type = T;
+    };
+
+    template <typename T>
+    struct DenseMatrixBase {
+        using Type = T;
+    };
+
+    template <typename T>
+    struct SparseMatrixBase{
+        using Type = T;
+    };
+
+    template <typename T, std::size_t Rank>
+    struct DenseTensorBase {
+        using Type = T;
+        static constexpr std::size_t rank = Rank;
+    };
+
+    template <typename T, std::size_t Rank>
+    struct SparseTensorBase {
+        using Type = T;
+        static constexpr std::size_t rank = Rank;
+    };
+}
diff --git a/src/include/serif/utils/user_options/mapping/compactification/kelvin.hpp b/src/include/serif/utils/user_options/mapping/compactification/kelvin.hpp
new file mode 100644
index 00000000..eb8f9202
--- /dev/null
+++ b/src/include/serif/utils/user_options/mapping/compactification/kelvin.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+namespace serif::utils::user_options {
+    struct KelvinCompactificationOptions
+    {
+        double r_star_ref{1.0};
+        double r_inf_ref{2.0};
+        double coordinate_tolerance{1.0e-12};
+    };
+
+}
\ No newline at end of file
diff --git a/src/lib/serif/eos/models/polytropic.cpp b/src/lib/serif/eos/models/polytropic.cpp
index cb887967..863a0aa3 100644
--- a/src/lib/serif/eos/models/polytropic.cpp
+++ b/src/lib/serif/eos/models/polytropic.cpp
@@ -104,24 +104,24 @@ namespace serif::eos::models {
         };
     }
 
-    [[nodiscard]] PartialDerivative<Density, SpecificEnthalpy> Polytrope::partial_derivative(DensityFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
+    [[nodiscard]] PartialDerivative<MassDensity, SpecificEnthalpy> Polytrope::partial_derivative(DensityFromSpecificEnthalpy, WithRespectTo<SpecificEnthalpy>, const SpecificEnthalpyValue specific_enthalpy) const {
         validate_finite<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(specific_enthalpy.value());
 
         if (specific_enthalpy.value() < 0.0) {
             // dρ/dh = 0
-            return PartialDerivative<Density, SpecificEnthalpy>{0.0};
+            return PartialDerivative<MassDensity, SpecificEnthalpy>{0.0};
         }
 
         if (specific_enthalpy.value() == 0.0) {
             // dρ/dh = n == 1 -> 1/hs
             // dρ/dh = else   -> 0
-            return PartialDerivative<Density, SpecificEnthalpy>{
+            return PartialDerivative<MassDensity, SpecificEnthalpy>{
                 m_polytropic_index == 1.0 ? 1.0 / m_enthalpy_scale : 0.0
             };
         }
 
         // dρ/dh = (n/hs) * (h/hs)^(n-1)
-        return PartialDerivative<Density, SpecificEnthalpy>{
+        return PartialDerivative<MassDensity, SpecificEnthalpy>{
             m_polytropic_index / m_enthalpy_scale * std::pow(specific_enthalpy.value() / m_enthalpy_scale, m_polytropic_index - 1.0)
         };
     }
@@ -133,16 +133,16 @@ namespace serif::eos::models {
         return PartialDerivative<Pressure, SpecificEnthalpy>{density.value()};
     }
 
-    [[nodiscard]] PartialDerivative<Pressure, Density> Polytrope::partial_derivative(PressureFromDensity, WithRespectTo<Density>, const DensityValue density) const {
+    [[nodiscard]] PartialDerivative<Pressure, MassDensity> Polytrope::partial_derivative(PressureFromDensity, WithRespectTo<MassDensity>, const DensityValue density) const {
         validate_nonnegativity<EOSEvaluationError, SERiFErrorCode::EOS_INVALID_THERMODYNAMIC_INPUT>(density.value());
 
         if (density.value() == 0.0) {
             // dP/dρ = 0
-            return PartialDerivative<Pressure, Density>{0.0};
+            return PartialDerivative<Pressure, MassDensity>{0.0};
         }
 
         // dP/dρ = K * (1 + 1/n) * (ρ)^(1/n)
-        return PartialDerivative<Pressure, Density>{
+        return PartialDerivative<Pressure, MassDensity>{
             m_polytropic_constant * (1.0 + 1.0 / m_polytropic_index) * std::pow(density.value(), 1.0 / m_polytropic_index)
         };
     }
diff --git a/src/lib/serif/mapping/compactification/kelvin.cpp b/src/lib/serif/mapping/compactification/kelvin.cpp
new file mode 100644
index 00000000..b7fc57e3
--- /dev/null
+++ b/src/lib/serif/mapping/compactification/kelvin.cpp
@@ -0,0 +1,269 @@
+#include <cmath>
+
+#include "serif/utils/checks/finite.hpp"
+#include "serif/mapping/compactification/kelvin.hpp"
+#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
+#include "serif/utils/error/errors.hpp"
+
+namespace serif::mapping::compactification {
+    using namespace serif::utils::user_options;
+    using namespace serif::utils::errors;
+    using namespace serif::utils::checks;
+
+    KelvinCompactification::KelvinCompactification(const KelvinCompactificationOptions& options) : m_options(options) {
+        if (!std::isfinite(m_options.r_star_ref)) {
+            serif_error<std::invalid_argument, SERiFErrorCode::NON_FINITE_MAPPING_INPUT>("Kelvin compactification requires a finite reference stellar radius. Received {} ", m_options.r_star_ref);
+        }
+
+        if (!std::isfinite(m_options.r_inf_ref)) {
+            serif_error<std::invalid_argument, SERiFErrorCode::NON_FINITE_MAPPING_INPUT>("Kelvin compactification requires a finite reference infinity radius. Received {} ", m_options.r_star_ref);
+        }
+
+        if (m_options.r_star_ref <= 0 || m_options.r_inf_ref < m_options.r_star_ref)
+        {
+            serif_error<std::invalid_argument, SERiFErrorCode::IMPROPERLY_ORDERED_REFERENCE_RADII>("Kelvin compactification requires 0 < r_star_ref < r_inf_ref. Got r_star_ref={} and r_inf_ref={}", m_options.r_star_ref, m_options.r_inf_ref);
+        }
+
+        if (!std::isfinite(m_options.coordinate_tolerance) || m_options.coordinate_tolerance < 0.0 || m_options.coordinate_tolerance > 1.0)
+        {
+            serif_error<std::invalid_argument, SERiFErrorCode::INVALID_COORDINATE_TOLERANCE>("Kelvin compactification coordinate tolerance must be finite and on [0, 1). Received: {}", m_options.coordinate_tolerance);
+        }
+    }
+
+    bool KelvinCompactification::is_affine_in_displacement() const noexcept {
+        return true; // Always true for kelvin compactification
+    }
+
+    MappingStatus KelvinCompactification::evaluate(const ExteriorMapInput& input, ExteriorMapResult& result) const {
+        const int dimension = input.reference_position.Size();
+
+        if (dimension <= 0 || input.displaced_position.Size() != dimension || input.compactification_coordinate_gradient.Size() != dimension) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        if (input.displacement_jacobian.Height() != dimension || input.displacement_jacobian.Width() != dimension) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        const bool reference_pos_is_finite = is_finite(input.reference_position);
+        const bool displaced_pos_is_finite = is_finite(input.displaced_position);
+        const bool coord_grad_is_finite = is_finite(input.compactification_coordinate_gradient);
+        const bool displacement_jacobian_is_finite = is_finite(input.displacement_jacobian);
+
+        if (!(reference_pos_is_finite && displaced_pos_is_finite && coord_grad_is_finite && displacement_jacobian_is_finite)) {
+            return MappingStatus::NON_FINITE_MAPPING_INPUT;
+        }
+
+        RadialFactors factors{};
+
+        // The key radial factors we use are the scale (which stretches the finite computational domain into the infinite physical domain) and the scale derivative (which is used to compute the mapping jacobian)
+        const MappingStatus factor_status = compute_radial_factors(input.compactification_coordinate, factors);
+        if (factor_status != MappingStatus::VALID) {
+            return factor_status;
+        }
+
+        result.physical_position.SetSize(dimension);
+        result.mapping_jacobian.SetSize(dimension, dimension);
+
+        for (int i = 0; i < dimension; ++i) {
+            // Note how the physical position is just the product of the displaced position and the scale factor.
+            result.physical_position(i) = factors.scale * input.displaced_position(i);
+
+            // The mapping jacobian comes from trivial application of the product rule
+            // recall: r_i = scale * x_i where r is the physical position and x is the displaced position.
+            // then we can differentiate wrt. X_j holding nothing fixed. Note the capital X here, this is the mesh coordinate not the displaced position.
+            // Lets call this jacobian F
+            // F = \frac{\partial r_i}{\partial X_{j}}
+            // F then tells us how the physical position changes as we move along mesh coordinates
+            // Lets then apply this to the function we have for the kelvin compactification
+            // F = scale * \frac{\partial x_i}{\partial X_j} + x_i * \frac{\partial scale}{\partial X_j}
+            // Below you can see the displacement jacobian (\frac{\partial x_i}{\partial X_j}) and the scale derivative (\frac{\partial scale}{\partial X_j}) being applied to compute the mapping jacobian.
+            for (int j = 0; j < dimension; ++j) {
+                const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
+                result.mapping_jacobian(i, j) = factors.scale * input.displacement_jacobian(i, j) + input.displaced_position(i) * scale_gradient;
+            }
+        }
+
+        if (!is_finite(result.physical_position) || !is_finite(result.mapping_jacobian)) {
+            return MappingStatus::NON_FINITE_MAPPING_RESULT;
+        }
+
+        const double mapping_determinant = result.mapping_jacobian.Det();
+
+        if (!std::isfinite(mapping_determinant)) {
+            return MappingStatus::NON_FINITE_MAPPING_RESULT;
+        }
+        if (mapping_determinant <= 0.0) {
+            return MappingStatus::NON_POSITIVE_DETERMINATE;
+        }
+        return MappingStatus::VALID;
+    }
+
+    MappingStatus KelvinCompactification::evaluate_variation(const ExteriorMapInput& input, const ExteriorMapResult& result, const ExteriorMapDirection& direction, ExteriorMapVariation& variation) const {
+        const int dimension = input.reference_position.Size();
+
+        if (dimension <= 0 || input.displaced_position.Size() != dimension || input.compactification_coordinate_gradient.Size() != dimension) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        if (input.displacement_jacobian.Height() != dimension || input.displacement_jacobian.Width() != dimension) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        if (result.physical_position.Size() != dimension || result.mapping_jacobian.Height() != dimension || result.mapping_jacobian.Width() != dimension) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        if (direction.displaced_position_variation.Size() != dimension || direction.displacement_jacobian_variation.Height() != dimension || direction.displacement_jacobian_variation.Width() != dimension) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        const bool reference_position_is_finite = is_finite(input.reference_position);
+        const bool displaced_position_is_finite = is_finite(input.displaced_position);
+        const bool compactification_coord_grad_is_finite = is_finite(input.compactification_coordinate_gradient);
+        const bool displacement_jacobian_is_finite = is_finite(input.displacement_jacobian);
+
+        if (!(reference_position_is_finite && displaced_position_is_finite && compactification_coord_grad_is_finite && displacement_jacobian_is_finite)) {
+            return MappingStatus::NON_FINITE_MAPPING_INPUT;
+        }
+
+        const bool result_physical_position_is_finite = is_finite(result.physical_position);
+        const bool result_mapping_jacobian_is_finite = is_finite(result.mapping_jacobian);
+        const bool direction_displaced_position_is_finite = is_finite(direction.displaced_position_variation);
+        const bool direction_displacement_jacobian_is_finite = is_finite(direction.displacement_jacobian_variation);
+        if (!(result_physical_position_is_finite && result_mapping_jacobian_is_finite && direction_displaced_position_is_finite && direction_displacement_jacobian_is_finite)) {
+            return MappingStatus::NON_FINITE_MAPPING_INPUT; // Note we are checking a parameter called result here; however, this is the result from an evaluate call, it is an input to this function.
+        }
+
+        RadialFactors factors{};
+        const MappingStatus factor_status = compute_radial_factors(input.compactification_coordinate, factors);
+        if (factor_status != MappingStatus::VALID) {
+            return factor_status;
+        }
+
+        variation.physical_position_variation.SetSize(dimension);
+        variation.mapping_jacobian_variation.SetSize(dimension, dimension);
+
+        for (int i = 0; i < dimension; ++i) {
+            variation.physical_position_variation(i) = factors.scale * direction.displaced_position_variation(i);
+
+            for (int j = 0; j < dimension; ++j) {
+                const double scale_gradient = factors.scale_derivative * input.compactification_coordinate_gradient(j);
+                variation.mapping_jacobian_variation(i, j) = factors.scale * direction.displacement_jacobian_variation(i, j) + direction.displaced_position_variation(i) * scale_gradient;
+            }
+        }
+
+        if (!is_finite(variation.physical_position_variation) || !is_finite(variation.mapping_jacobian_variation)) {
+            return MappingStatus::NON_FINITE_MAPPING_RESULT;
+        }
+
+        return MappingStatus::VALID;
+    }
+
+    [[nodiscard]] std::string_view KelvinCompactification::get_name() const noexcept {
+        return m_name;
+    }
+
+    [[nodiscard]] double KelvinCompactification::get_reference_stellar_radius() const noexcept {
+        return m_options.r_star_ref;
+    }
+
+    [[nodiscard]] double KelvinCompactification::get_reference_infinity_radius() const noexcept {
+        return m_options.r_inf_ref;
+    }
+
+    [[nodiscard]] double KelvinCompactification::get_coordinate_tolerance() const noexcept {
+        return m_options.coordinate_tolerance;
+    }
+
+    [[nodiscard]] MappingStatus KelvinCompactification::evaluate_jacobian_variation_factors(double compactification_coordinate, const mfem::Vector& coordinateGradient, double& scale, mfem::Vector& scaleGradient) const {
+        if (coordinateGradient.Size() <= 0) {
+            return MappingStatus::INVALID_DIMENSION;
+        }
+
+        if (!is_finite(coordinateGradient)) {
+            return MappingStatus::NON_FINITE_MAPPING_INPUT;
+        }
+
+        RadialFactors factors{};
+        const MappingStatus status = compute_radial_factors(compactification_coordinate, factors);
+
+        if (status != MappingStatus::VALID) {
+            return status;
+        }
+
+        scale = factors.scale;
+        scaleGradient = coordinateGradient;
+        scaleGradient *= factors.scale_derivative;
+
+        if (!is_finite(scaleGradient)) {
+            return MappingStatus::NON_FINITE_MAPPING_RESULT;
+        }
+
+        return MappingStatus::VALID;
+    }
+
+    MappingStatus KelvinCompactification::compute_radial_factors(double compactification_coordinate, RadialFactors& radial_factors) const {
+        if (!std::isfinite(compactification_coordinate)) {
+            return MappingStatus::NON_FINITE_MAPPING_INPUT;
+        }
+
+        const double tolerance = m_options.coordinate_tolerance;
+
+        if (compactification_coordinate < -tolerance || compactification_coordinate > 1.0 + tolerance) {
+            return MappingStatus::OUTSIDE_REFERENCE_DOMAIN;
+        }
+
+        double coordinate = compactification_coordinate;
+        if (coordinate < 0.0) {
+            coordinate = 0.0; // Not sure this is the best behavior here...
+        }
+
+        if (coordinate >= 1.0 - tolerance) {
+            return MappingStatus::AT_COMPACTIFIED_INFINITY;
+        }
+
+        // Here we need to do some transformations from the options defined on the mesh to useful computational coordinates
+        // r_inf_ref is the computational / reference radius of the infinity surface (the edge of the entire domain) and r_star_ref is the radius of the spherical
+        // stellar model inscribed within. Therefore radial extent is the computational radial distance between the stellar surface and the infinity surface.
+        // Note that this is separate from the parameterize compactification coordinate.
+        const double radial_extent = m_options.r_inf_ref - m_options.r_star_ref;
+
+        // This places us at the correct spot in computational space given the current compactification coordinate. Say you have compactification = 0.5,
+        // an r_star_ref of 2 and a radial extent of 3, this this will place you at 2 + 0.5 * 3 = 3.5 in computational space, which is half way between the stellar surface and the infinity surface.
+        const double computational_radius = m_options.r_star_ref + coordinate * radial_extent;
+
+        if (!std::isfinite(computational_radius) || computational_radius <= 0.0) {
+            return MappingStatus::INVALID_REFERENCE_RADIUS;
+        }
+
+        // Invert the exterior coordinate so it runs from 0 at the star to 1 at compactified infinity
+        const double one_minus_coordinate = 1.0 - coordinate;
+        const double denominator = computational_radius * one_minus_coordinate;
+
+        if (!std::isfinite(denominator) || denominator <= 0.0) {
+            return MappingStatus::NON_FINITE_MAPPING_RESULT;
+        }
+
+        // The scale here is the factor which stretches the finite computational domain into the infinite physical domain. Properties we need this to have
+        // include that it should go to 1 at the stellar surface and go to infinity at the compactified infinity.
+        // Mathematically this is scale = |r|/|x| where r is the physical radius and x is the computational radius.
+        // Put another way, scale is the ratio of the target physical radius for the current compactification coordinate
+        // to the current mesh radius.
+        const double scale = m_options.r_star_ref / denominator;
+        const double scale_derivative = scale * (1.0 / one_minus_coordinate - radial_extent / computational_radius);
+
+        if (!std::isfinite(scale) || !std::isfinite(scale_derivative)) {
+            return MappingStatus::NON_FINITE_MAPPING_RESULT;
+        }
+
+        radial_factors.compactification_coordinate           = coordinate;
+        radial_factors.computational_radius = computational_radius;
+        radial_factors.scale                = scale;
+        radial_factors.scale_derivative     = scale_derivative;
+
+        return MappingStatus::VALID;
+
+
+    }
+}
diff --git a/src/meson.build b/src/meson.build
index 7bad5497..6d51e18e 100644
--- a/src/meson.build
+++ b/src/meson.build
@@ -5,6 +5,7 @@ serif_sources = files(
   'lib/serif/discretization/domain/schema/validation/results.cpp',
   'lib/serif/discretization/quadrature/policy.cpp',
   'lib/serif/eos/models/polytropic.cpp',
+  'lib/serif/mapping/compactification/kelvin.cpp',
 )
 
 serif_lib = library(
@@ -45,6 +46,6 @@ if get_option('build_checks')
     'serif_checks',
     checks_srcs,
     include_directories: [include_dir, nameof_include_dir, magic_enum_include] ,
-    dependencies: [mfem_dep],
+    dependencies: [mfem_dep, magic_enum_dep],
   )
 endif
diff --git a/subprojects/cli11.wrap b/subprojects/cli11.wrap
new file mode 100644
index 00000000..00725907
--- /dev/null
+++ b/subprojects/cli11.wrap
@@ -0,0 +1,10 @@
+[wrap-file]
+directory = CLI11-2.6.1
+source_url = https://github.com/CLIUtils/CLI11/archive/refs/tags/v2.6.1.tar.gz
+source_filename = CLI11-2.6.1.tar.gz
+source_hash = 377691f3fac2b340f12a2f79f523c780564578ba3d6eaf5238e9f35895d5ba95
+source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/cli11_2.6.1-1/CLI11-2.6.1.tar.gz
+wrapdb_version = 2.6.1-1
+
+[provide]
+dependency_names = CLI11
\ No newline at end of file
diff --git a/subprojects/libconfig.wrap b/subprojects/libconfig.wrap
new file mode 100644
index 00000000..86b0e32e
--- /dev/null
+++ b/subprojects/libconfig.wrap
@@ -0,0 +1,4 @@
+[wrap-git]
+url = https://github.com/4D-STAR/libconfig.git
+revision = 11ed6d5a0f57b3516ee2cadf1bce20d380fa07ea
+depth = 1
diff --git a/subprojects/stroid.wrap b/subprojects/stroid.wrap
new file mode 100644
index 00000000..e72f2cc7
--- /dev/null
+++ b/subprojects/stroid.wrap
@@ -0,0 +1,8 @@
+[wrap-git]
+url = https://github.com/4D-STAR/stroid.git
+revision = 10d8b17c63512bcb348100346b8a9bfcba1d148a
+depth = 1
+patch_directory = stroid
+
+[provide]
+stroid = stroid_dep
diff --git a/tests/meson.build b/tests/meson.build
index dfb28571..26f9bd6b 100644
--- a/tests/meson.build
+++ b/tests/meson.build
@@ -4,10 +4,19 @@ if get_option('build_tests')
 #    test_include_dir = include_directories('include')
 #
 #    catch2_dep = dependency('Catch2', required:true)
-    stroid_dep = dependency('stroid', required:true)
+    # Stroid exchanges MFEM and standard-library objects with the sandboxes.
+    # Always build it against the parent's MFEM and C++ ABI, never an unrelated
+    # installed archive (even when pkg-config reports one).
+    stroid_project = subproject('stroid', default_options: [
+      'build_tests=false',
+      'build_tools=false',
+      'build_python=false',
+      'pkg_config=false',
+    ])
+    stroid_dep = stroid_project.get_variable('stroid_dep')
 #
 #    test_sources = discretization_test_sources
 #    executable('serif_tests', test_sources, dependencies: [serif_dep, catch2_dep, stroid_dep], include_directories: test_include_dir)
 
     subdir('sandbox')
-endif
\ No newline at end of file
+endif
diff --git a/tests/sandbox/eos_sandbox.cpp b/tests/sandbox/eos_sandbox.cpp
index 0fe42b86..06b2630a 100644
--- a/tests/sandbox/eos_sandbox.cpp
+++ b/tests/sandbox/eos_sandbox.cpp
@@ -10,17 +10,15 @@
 #include "serif/eos/runtime/concepts.hpp"
 #include "serif/eos/runtime/views.hpp"
 int main() {
-    const serif::eos::models::Polytrope polytrope(1, 0.6);
-    std::println("Density: {}", serif::eos::evaluate<serif::dimensions::Density>(polytrope, serif::dimensions::SpecificEnthalpyValue{1.0}).value());
-    const serif::eos::runtime::EOSView view(polytrope);
-    const auto result = view.try_evaluate<serif::dimensions::Density>(serif::dimensions::SpecificEnthalpyValue{1.0});
-    if (result.has_value()) {
-        std::println("From runtime view: Density: {}", result->value());
-    } else {
-        std::println("From runtime view: EOS view is unavailable.");
-    }
-
-    std::println("{}", serif::discretization::blocks::gravity::name);
-
-
+    const serif::eos::models::Polytrope polytrope(-1, 0.6);
+    // std::println("Density: {}", serif::eos::evaluate<serif::dimensions::Density>(polytrope, serif::dimensions::SpecificEnthalpyValue{1.0}).value());
+    // const serif::eos::runtime::EOSView view(polytrope);
+    // const auto result = view.try_evaluate<serif::dimensions::Density>(serif::dimensions::SpecificEnthalpyValue{1.0});
+    // if (result.has_value()) {
+    //     std::println("From runtime view: Density: {}", result->value());
+    // } else {
+    //     std::println("From runtime view: EOS view is unavailable.");
+    // }
+    //
+    // std::println("{}", serif::discretization::blocks::gravity::name);
 }
\ No newline at end of file
diff --git a/tests/sandbox/mapping_sandbox.cpp b/tests/sandbox/mapping_sandbox.cpp
new file mode 100644
index 00000000..ed16a5f6
--- /dev/null
+++ b/tests/sandbox/mapping_sandbox.cpp
@@ -0,0 +1,12 @@
+#include <print>
+
+#include "serif/mapping/context.hpp"
+#include "serif/mapping/variation.hpp"
+#include "serif/mapping/compactification/compactification.hpp"
+#include "serif/mapping/status.hpp"
+#include "serif/mapping/compactification/kelvin.hpp"
+
+#include "serif/utils/user_options/mapping/compactification/kelvin.hpp"
+int main() {
+    std::println("Hello, World!");
+}
\ No newline at end of file
diff --git a/tests/sandbox/meson.build b/tests/sandbox/meson.build
index c1439a68..88ff87e7 100644
--- a/tests/sandbox/meson.build
+++ b/tests/sandbox/meson.build
@@ -5,13 +5,33 @@ sandboxes = [
     'dimensions_sandbox',
     'quadrature_sandbox',
     'meta_sandbox',
+    'mapping_sandbox'
 ]
 
+expected_errors = {}
+foreach entry : get_option('sandbox_expected_errors')
+  parts = entry.split(':')
+  if parts.length() != 2 or parts[0] not in sandboxes or parts[1] == ''
+    error('Invalid sandbox_expected_errors entry: ' + entry)
+  endif
+  expected_errors += {parts[0]: parts[1]}
+endforeach
+
 foreach sandbox: sandboxes
-  executable(
+  sandbox_exe = executable(
     sandbox,
     sandbox + '.cpp',
     dependencies: [serif_dep, stroid_dep, nameof_dep],
     build_rpath: mfem_runtime_prefix == '' ? '' : mfem_runtime_prefix / 'lib',
   )
+  if expected_errors.has_key(sandbox)
+    test(
+      sandbox, python_build,
+      args: [files('../../tools/check_sandbox_error.py'), sandbox_exe, expected_errors[sandbox]],
+      workdir: meson.project_source_root(),
+      suite: 'sandbox',
+    )
+  else
+    test(sandbox, sandbox_exe, workdir: meson.project_source_root(), suite: 'sandbox')
+  endif
 endforeach
diff --git a/tools/build_mfem_bundle.py b/tools/build_mfem_bundle.py
index 62e49b44..476a3764 100644
--- a/tools/build_mfem_bundle.py
+++ b/tools/build_mfem_bundle.py
@@ -384,7 +384,7 @@ def common_cmake_args(
             )
     c_flags = list(args.c_arg)
     cxx_flags = list(args.cxx_arg)
-    link_flags = list(dict.fromkeys([*args.c_link_arg, *args.cxx_link_arg]))
+    link_flags = []
     if args.host_system == "darwin":
         link_flags.append("-Wl,-headerpad_max_install_names")
     if args.host_system == "emscripten":
@@ -410,10 +410,24 @@ def common_cmake_args(
         result.append(f"-DCMAKE_C_FLAGS={joined_flags(c_flags)}")
     if cxx_flags:
         result.append(f"-DCMAKE_CXX_FLAGS={joined_flags(cxx_flags)}")
-    if link_flags:
-        rendered_link_flags = joined_flags(link_flags)
-        result.append(f"-DCMAKE_EXE_LINKER_FLAGS={rendered_link_flags}")
-        result.append(f"-DCMAKE_SHARED_LINKER_FLAGS={rendered_link_flags}")
+    # CMake's global linker flags apply to both compiler drivers. In particular,
+    # passing Clang's -stdlib=libc++ to a GCC C compiler breaks its ABI check.
+    # Apply Meson's language-specific options only to the matching link language.
+    link_options_file = args.work_dir / "meson-link-options.cmake"
+    link_options_file.parent.mkdir(parents=True, exist_ok=True)
+    link_options = []
+    for language, flags in (("C", args.c_link_arg), ("CXX", args.cxx_link_arg)):
+        for flag in flags:
+            escaped = (flag.replace("\\", "\\\\").replace('"', '\\"')
+                       .replace(";", "\\;").replace("$", "\\$"))
+            # Libraries must follow the object files for --as-needed linkers.
+            command = "link_libraries" if flag.startswith("-l") else "add_link_options"
+            link_options.append(f'{command}("$<$<LINK_LANGUAGE:{language}>:{escaped}>")')
+    link_options_file.write_text("\n".join(link_options) + "\n")
+    result.append(f"-DCMAKE_PROJECT_INCLUDE={link_options_file}")
+    rendered_link_flags = joined_flags(link_flags)
+    for kind in ("EXE", "SHARED", "MODULE"):
+        result.append(f"-DCMAKE_{kind}_LINKER_FLAGS={rendered_link_flags}")
     rpath = runtime_rpath(args.host_system)
     if rpath:
         result.extend(
@@ -586,7 +600,7 @@ def build_mpich(
 
     local_env = dict(env)
     mpich_link_flags = [f"-L{args.prefix / 'lib'}"]
-    mpich_link_flags.extend([*args.c_link_arg, *args.cxx_link_arg])
+    mpich_link_flags.extend(args.c_link_arg)
     if args.host_system == "darwin":
         mpich_link_flags.append("-Wl,-headerpad_max_install_names")
     local_env.update(
@@ -597,7 +611,7 @@ def build_mpich(
                 filter(None, [compiler_flags(build_type), joined_flags(args.c_arg)])
             ),
             "CXXFLAGS": " ".join(
-                filter(None, [compiler_flags(build_type), joined_flags(args.cxx_arg)])
+                filter(None, [compiler_flags(build_type), joined_flags([*args.cxx_arg, *args.cxx_link_arg])])
             ),
             "CPPFLAGS": f"-I{args.prefix / 'include'}",
             "LDFLAGS": " ".join(mpich_link_flags),
@@ -966,7 +980,8 @@ def build_libceed(
             ],
         )
     )
-    link_flags = [*args.c_link_arg, *args.cxx_link_arg]
+    # libCEED's shared library is linked by CC, including its CPU backends.
+    link_flags = list(args.c_link_arg)
     if features["openmp"]:
         # libCEED's default clang flag (-fopenmp) does not work with Apple Clang.
         command.append(f"OMP_FLAG={joined_flags(args.openmp_c_flags)}")
diff --git a/tools/check_sandbox_error.py b/tools/check_sandbox_error.py
new file mode 100644
index 00000000..b737d8e0
--- /dev/null
+++ b/tools/check_sandbox_error.py
@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+"""Verify a sandbox's deliberate uncaught exception, never arbitrary failure."""
+
+import re
+import resource
+import signal
+import subprocess
+import sys
+
+
+def main() -> int:
+    executable, expected_code = sys.argv[1:]
+    # These examples deliberately terminate; avoid writing core dumps.
+    resource.setrlimit(resource.RLIMIT_CORE, (0, 0))
+    result = subprocess.run([executable], capture_output=True, text=True, timeout=25)
+    sys.stdout.write(result.stdout)
+    sys.stderr.write(result.stderr)
+    diagnostic = rf"Error: .*\(code No\.: \d+ \[{re.escape(expected_code)}\]\)"
+    if result.returncode == -signal.SIGABRT and re.search(diagnostic, result.stderr):
+        print(f"Verified expected exception: {expected_code}")
+        return 0
+    print(
+        f"Expected SIGABRT with SERiF error {expected_code}; got exit {result.returncode}",
+        file=sys.stderr,
+    )
+    return 1
+
+
+if __name__ == '__main__':
+    raise SystemExit(main())
diff --git a/tools/validate_matrix.sh b/tools/validate_matrix.sh
index 27fdd205..000b363e 100755
--- a/tools/validate_matrix.sh
+++ b/tools/validate_matrix.sh
@@ -1,66 +1,42 @@
 #!/bin/sh
+# Build and run all six sandboxes with each Linux compiler/library/linker pair.
 set -eu
 
 project_root=$(CDPATH= cd -- "$(dirname -- "$0")/.." && pwd)
 meson_bin=${MESON_BIN:-meson}
-meson_dir=$(CDPATH= cd -- "$(dirname -- "$meson_bin")" && pwd)
-profile=${MFEM_MATRIX_PROFILE:-portable}
+profile=${MFEM_MATRIX_PROFILE:-minimal}
 jobs=${MFEM_JOBS:-4}
+mode=${SERIF_BUILDTYPE:-debug}
 
-compiler_family() {
-    "$1" --version 2>&1 | sed -n '1p'
-}
-
-require_clang() {
-    version=$(compiler_family "$1")
-    case "$version" in
-        *clang*|*Clang*) ;;
-        *) echo "error: $1 is not a Clang compiler: $version" >&2; exit 2 ;;
+# GCC/libc++ requires a recent GCC and Linux headers at /usr/include/c++/v1.
+# Set SERIF_TOOLCHAINS to a space-separated subset to test selected toolchains.
+for toolchain in ${SERIF_TOOLCHAINS:-gcc gcc-libcxx clang-libstdcxx clang-libcxx}; do
+    case "$toolchain" in
+        gcc) compiler=gcc; library= ;;
+        gcc-libcxx) compiler=gcc; library=gcc-libcxx ;;
+        clang-libstdcxx) compiler=clang; library=libstdcxx ;;
+        clang-libcxx) compiler=clang; library=libcxx ;;
+        *) echo "Unknown toolchain: $toolchain" >&2; exit 2 ;;
     esac
-}
-
-require_gnu() {
-    version=$(compiler_family "$1")
-    case "$version" in
-        *gcc*|*GCC*|*"Free Software Foundation"*) ;;
-        *) echo "error: $1 is not a GNU compiler: $version" >&2; exit 2 ;;
-    esac
-}
-
-run_configuration() {
-    name=$1
-    c_compiler=$2
-    cxx_compiler=$3
-    mode=$4
-    build_dir="$project_root/build-$name"
-
-    if test -f "$build_dir/meson-private/coredata.dat"; then
-        wipe_flag=--wipe
-    else
-        wipe_flag=
-    fi
-    env PATH="$meson_dir:$PATH" CC="$c_compiler" CXX="$cxx_compiler" \
-        "$meson_bin" setup $wipe_flag "$build_dir" "$project_root" \
-            --buildtype="$mode" \
-            -Dallow_preinstalled=false \
-            -Dfeature_profile="$profile" \
-            -Dbuild_examples=true \
-            -Dbuild_tests=true \
-            -Dbuild_python=false \
-            -Djobs="$jobs"
-    env PATH="$meson_dir:$PATH" "$meson_bin" compile -C "$build_dir"
-    env PATH="$meson_dir:$PATH" "$meson_bin" test -C "$build_dir" --print-errorlogs
-}
-
-clang_cc=${CLANG_CC:-clang}
-clang_cxx=${CLANG_CXX:-clang++}
-gcc_cc=${GNU_CC:-gcc}
-gcc_cxx=${GNU_CXX:-g++}
-
-require_clang "$clang_cxx"
-require_gnu "$gcc_cxx"
-
-run_configuration clang-debug "$clang_cc" "$clang_cxx" debug
-run_configuration clang-release "$clang_cc" "$clang_cxx" release
-run_configuration gcc-debug "$gcc_cc" "$gcc_cxx" debug
-run_configuration gcc-release "$gcc_cc" "$gcc_cxx" release
+    for linker in ${SERIF_LINKERS:-ld lld}; do
+        case "$linker" in ld|lld) ;; *) echo "Unknown linker: $linker" >&2; exit 2 ;; esac
+        build_dir="$project_root/build/matrix/$toolchain-$linker"
+        set -- "$build_dir" "$project_root" \
+            "--native-file=$project_root/native/$compiler.ini" \
+            "--native-file=$project_root/native/$linker.ini" \
+            "--buildtype=$mode" -Dallow_preinstalled=false \
+            "-Dfeature_profile=$profile" -Dbuild_tests=true "-Djobs=$jobs"
+        if test -n "$library"; then
+            set -- "$@" "--native-file=$project_root/native/$library.ini"
+        fi
+        if test -n "${SERIF_SANDBOX_EXPECTED_ERRORS:-}"; then
+            set -- "$@" "-Dsandbox_expected_errors=$SERIF_SANDBOX_EXPECTED_ERRORS"
+        fi
+        if test -f "$build_dir/meson-private/coredata.dat"; then
+            set -- --reconfigure "$@"
+        fi
+        "$meson_bin" setup "$@"
+        "$meson_bin" compile -C "$build_dir" -j "$jobs"
+        "$meson_bin" test -C "$build_dir" --print-errorlogs
+    done
+done