diff --git a/builds/make.type.disk b/builds/make.type.disk
index c2efaf0f6..f2e3f0ec1 100644
--- a/builds/make.type.disk
+++ b/builds/make.type.disk
@@ -21,7 +21,6 @@ DFLAGS += -DGRAVITY_5_POINTS_GRADIENT
 
 #DFLAGS += -DSTATIC_GRAV
 
-#DFLAGS    += -DOUTPUT_ALWAYS
 DFLAGS    += -DCUDA
 DFLAGS    += -DMPI_CHOLLA
 DFLAGS    += -DPRECISION=2
diff --git a/builds/make.type.mhd b/builds/make.type.mhd
index 1849722d7..4459819e8 100644
--- a/builds/make.type.mhd
+++ b/builds/make.type.mhd
@@ -51,6 +51,3 @@ DFLAGS    += $(MPI_GPU)
 
 # Limit the number of steps to evolve.
 # DFLAGS += -DN_STEPS_LIMIT=1000
-
-# Output on every time step
-# DFLAGS += -DOUTPUT_ALWAYS
diff --git a/cholla-tests-data b/cholla-tests-data
index dcd73ff52..71eb66d63 160000
--- a/cholla-tests-data
+++ b/cholla-tests-data
@@ -1 +1 @@
-Subproject commit dcd73ff52b9027627b247c6d888bcdb56840c85e
+Subproject commit 71eb66d63622ac15c0844ae96ec9034cd5e4f4d3
diff --git a/src/dust/dust_cuda.cu b/src/dust/dust_cuda.cu
index d4c6191ec..28f764900 100644
--- a/src/dust/dust_cuda.cu
+++ b/src/dust/dust_cuda.cu
@@ -25,17 +25,20 @@
   #include "../utils/gpu.hpp"
   #include "../utils/hydro_utilities.h"
 
-void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma)
+void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                 Real grain_radius)
 {
   int n_cells = nx * ny * nz;
   int ngrid   = (n_cells + TPB - 1) / TPB;
   dim3 dim1dGrid(ngrid, 1, 1);
   dim3 dim1dBlock(TPB, 1, 1);
-  hipLaunchKernelGGL(Dust_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, dt, gamma);
+  hipLaunchKernelGGL(Dust_Kernel, dim1dGrid, dim1dBlock, 0, 0, dev_conserved, nx, ny, nz, n_ghost, n_fields, dt, gamma,
+                     grain_radius);
   GPU_Error_Check();
 }
 
-__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma)
+__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                            Real grain_radius)
 {
   // get grid indices
   int n_cells = nx * ny * nz;
@@ -99,8 +102,8 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
     // if dual energy is turned on use temp from total internal energy
     temperature = temperature_init;
 
-    Real tau_sp =
-        Calc_Sputtering_Timescale(number_density, temperature) / TIME_UNIT;  // sputtering timescale, kyr (sim units)
+    Real tau_sp = Calc_Sputtering_Timescale(number_density, temperature, grain_radius) /
+                  TIME_UNIT;  // sputtering timescale, kyr (sim units)
 
     dd_dt = Calc_dd_dt(density_dust, tau_sp);  // rate of change in dust density at current timestep
     dd    = dd_dt * dt;                        // change in dust density at current timestep
@@ -126,17 +129,18 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
 }
 
 // McKinnon et al. (2017) sputtering timescale
-__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature)
+__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature, Real grain_radius)
 {
-  Real grain_radius  = 1;          // dust grain size in units of 0.1 micrometers
-  Real temperature_0 = 2e6;        // temp above which the sputtering rate is ~constant in K
-  Real omega         = 2.5;        // controls the low-temperature scaling of the sputtering rate
-  Real A             = 5.3618e15;  // 0.17 Gyr in s
+  Real a             = grain_radius;  // dust grain size in units of 0.1 micrometers
+  Real temperature_0 = 2e6;           // temp above which the sputtering rate is ~constant in K
+  Real omega         = 2.5;           // controls the low-temperature scaling of the sputtering rate
+  Real A             = 5.3618e15;     // 0.17 Gyr in s
 
   number_density /= (6e-4);  // gas number density in units of 10^-27 g/cm^3
 
   // sputtering timescale, s
-  Real tau_sp = A * (grain_radius / number_density) * (pow(temperature_0 / temperature, omega) + 1);
+  printf("%e\n", grain_radius);
+  Real tau_sp = A * (a / number_density) * (pow(temperature_0 / temperature, omega) + 1);
 
   return tau_sp;
 }
diff --git a/src/dust/dust_cuda.h b/src/dust/dust_cuda.h
index fb72007ac..212901e8a 100644
--- a/src/dust/dust_cuda.h
+++ b/src/dust/dust_cuda.h
@@ -27,7 +27,8 @@
  * \param[in] dt Simulation timestep
  * \param[in] gamma Specific heat ratio
  */
-void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma);
+void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                 Real grain_radius);
 
 /*!
  * \brief Compute the change in dust density for a cell and update its value in dev_conserved.
@@ -42,8 +43,8 @@ void Dust_Update(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n
  * \param[in] dt Simulation timestep
  * \param[in] gamma Specific heat ratio
  */
-__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt,
-                            Real gamma);
+__global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_ghost, int n_fields, Real dt, Real gamma,
+                            Real grain_radius);
 
 /*!
  * \brief Compute the sputtering timescale based on a cell's density and temperature.
@@ -53,7 +54,7 @@ __global__ void Dust_Kernel(Real *dev_conserved, int nx, int ny, int nz, int n_g
  *
  * \return Real Sputtering timescale in seconds (McKinnon et al. 2017)
  */
-__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature);
+__device__ __host__ Real Calc_Sputtering_Timescale(Real number_density, Real temperature, Real grain_radius);
 
 /*!
  * \brief Compute the rate of change in dust density based on the current dust density and sputtering timescale.
diff --git a/src/dust/dust_cuda_tests.cpp b/src/dust/dust_cuda_tests.cpp
index aa26f7a01..5b59b2dc0 100644
--- a/src/dust/dust_cuda_tests.cpp
+++ b/src/dust/dust_cuda_tests.cpp
@@ -27,9 +27,11 @@ TEST(tDUSTTestSputteringTimescale, CorrectInputExpectCorrectOutput)
   Real YR_IN_S                     = 3.154e7;
   Real const k_test_number_density = 1;
   Real const k_test_temperature    = pow(10, 5.0);
+  Real const k_test_grain_radius   = 1;
   Real const k_fiducial_num        = 182565146.96398282;
 
-  Real test_num = Calc_Sputtering_Timescale(k_test_number_density, k_test_temperature) / YR_IN_S;  // yr
+  Real test_num =
+      Calc_Sputtering_Timescale(k_test_number_density, k_test_temperature, k_test_grain_radius) / YR_IN_S;  // yr
 
   double abs_diff;
   int64_t ulps_diff;
diff --git a/src/global/global.cpp b/src/global/global.cpp
index 6579d9955..1a867965c 100644
--- a/src/global/global.cpp
+++ b/src/global/global.cpp
@@ -240,6 +240,10 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
   } else if (strcmp(name, "out_float32_GasEnergy") == 0) {
     parms->out_float32_GasEnergy = atoi(value);
 #endif  // DE
+  } else if (strcmp(name, "output_always") == 0) {
+    int tmp = atoi(value);
+    CHOLLA_ASSERT((tmp == 0) or (tmp == 1), "output_always must be 1 or 0.");
+    parms->output_always = tmp;
 #ifdef MHD
   } else if (strcmp(name, "out_float32_magnetic_x") == 0) {
     parms->out_float32_magnetic_x = atoi(value);
@@ -453,6 +457,12 @@ void Parse_Param(char *name, char *value, struct Parameters *parms)
   } else if (strcmp(name, "skewersdir") == 0) {
     strncpy(parms->skewersdir, value, MAXLEN);
   #endif
+#endif
+#ifdef SCALAR
+  #ifdef DUST
+  } else if (strcmp(name, "grain_radius") == 0) {
+    parms->grain_radius = atoi(value);
+  #endif
 #endif
   } else if (!Is_Param_Valid(name)) {
     chprintf("WARNING: %s/%s: Unknown parameter/value pair!\n", name, value);
diff --git a/src/global/global.h b/src/global/global.h
index 6c9524001..3cf58312a 100644
--- a/src/global/global.h
+++ b/src/global/global.h
@@ -204,6 +204,7 @@ struct Parameters {
 #ifdef DE
   int out_float32_GasEnergy = 0;
 #endif
+  bool output_always = false;
 #ifdef STATIC_GRAV
   int custom_grav = 0;  // flag to set specific static gravity field
 #endif
@@ -341,6 +342,11 @@ struct Parameters {
   char skewersdir[MAXLEN];
   #endif
 #endif
+#ifdef SCALAR
+  #ifdef DUST
+  Real grain_radius;
+  #endif
+#endif
 };
 
 /*! \fn void parse_params(char *param_file, struct Parameters * parms);
diff --git a/src/grid/grid3D.cpp b/src/grid/grid3D.cpp
index 93499fea0..f36d74c1d 100644
--- a/src/grid/grid3D.cpp
+++ b/src/grid/grid3D.cpp
@@ -280,6 +280,12 @@ void Grid3D::Initialize(struct Parameters *P)
   H.OUTPUT_SCALE_FACOR = not(P->scale_outputs_file[0] == '\0');
 #endif
 
+#ifdef SCALAR
+  #ifdef DUST
+  H.grain_radius = P->grain_radius;
+  #endif
+#endif
+
   H.Output_Initial = true;
 }
 
@@ -524,7 +530,7 @@ Real Grid3D::Update_Hydro_Grid()
 
   #ifdef DUST
   // ==Apply dust from dust/dust_cuda.h==
-  Dust_Update(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama);
+  Dust_Update(C.device, H.nx, H.ny, H.nz, H.n_ghost, H.n_fields, H.dt, gama, H.grain_radius);
   #endif  // DUST
 
 #endif  // CUDA
diff --git a/src/grid/grid3D.h b/src/grid/grid3D.h
index 5354acf0a..e248f6490 100644
--- a/src/grid/grid3D.h
+++ b/src/grid/grid3D.h
@@ -268,6 +268,12 @@ struct Header {
    *  \brief Flag set to true when all the data will  be written to file
    * (Restart File ) */
   bool Output_Complete_Data;
+
+#ifdef SCALAR
+  #ifdef DUST
+  Real grain_radius;
+  #endif
+#endif
 };
 
 /*! \class Grid3D
diff --git a/src/main.cpp b/src/main.cpp
index cb6d1a72b..016df3f61 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -325,9 +325,7 @@ int main(int argc, char *argv[])
         "%9.3f ms   total time = %9.4f s\n\n",
         G.H.n_step, G.H.t, G.H.dt, (stop_step - start_step) * 1000, G.H.t_wall);
 
-#ifdef OUTPUT_ALWAYS
-    G.H.Output_Now = true;
-#endif
+    if (P.output_always) G.H.Output_Now = true;
 
 #ifdef ANALYSIS
     if (G.Analysis.Output_Now) {
diff --git a/src/reconstruction/plmc_cuda.cu b/src/reconstruction/plmc_cuda.cu
index 41a5ae505..bb31e9904 100644
--- a/src/reconstruction/plmc_cuda.cu
+++ b/src/reconstruction/plmc_cuda.cu
@@ -124,6 +124,9 @@ __global__ __launch_bounds__(TPB) void PLMC_cuda(Real *dev_conserved, Real *dev_
   reconstruction::Primitive interface_L_iph = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, 1.0);
   reconstruction::Primitive interface_R_imh = reconstruction::Calc_Interface_Linear(cell_i, del_m_i, -1.0);
 
+  // Limit the interfaces
+  reconstruction::Plm_Limit_Interfaces(interface_L_iph, interface_R_imh, cell_imo, cell_i, cell_ipo);
+
 #ifndef VL
 
   Real const dtodx = dt / dx;
diff --git a/src/reconstruction/plmc_cuda_tests.cu b/src/reconstruction/plmc_cuda_tests.cu
index 0c567b6d0..68f11b396 100644
--- a/src/reconstruction/plmc_cuda_tests.cu
+++ b/src/reconstruction/plmc_cuda_tests.cu
@@ -204,17 +204,17 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput)
                                                                            {{21, 0.73640639402573249},
                                                                             {85, 3.3462413154443715},
                                                                             {149, 2.1945584994458125},
-                                                                            {213, 0.67418839414138987},
-                                                                            {277, 16.909618487528142},
-                                                                            {341, 2.1533768050263267},
-                                                                            {405, 1.6994195863331925}},
+                                                                            {213, 1.1837630990406585},
+                                                                            {277, 17.570011907061254},
+                                                                            {341, 2.1583975283044725},
+                                                                            {405, 1.7033818819502551}},
                                                                            {{21, 0.25340904981266843},
                                                                             {85, 2.0441984720128734},
                                                                             {149, 1.9959059157695584},
                                                                             {213, 0.45377591914009824},
-                                                                            {277, 23.677832869261188},
-                                                                            {341, 1.5437923271692418},
-                                                                            {405, 1.8141353672443383}}};
+                                                                            {277, 24.018953780483471},
+                                                                            {341, 1.7033818819502551},
+                                                                            {405, 1.8587936590169301}}};
   std::vector<std::unordered_map<int, double>> fiducial_interface_right = {{{20, 0.59023012197434721},
                                                                             {84, 3.0043379408547275},
                                                                             {148, 2.6320759184913625},
@@ -227,18 +227,18 @@ TEST(tMHDPlmcReconstructor, CorrectInputExpectCorrectOutput)
                                                                                {81, 2.5027813113931279},
                                                                                {145, 2.6371119205792346},
                                                                                {209, 1.0210845222961809},
-                                                                               {273, 21.360010722689488},
+                                                                               {273, 21.353253570231175},
                                                                                {337, 2.1634182515826184},
-                                                                               {401, 1.7073441775673177},
+                                                                               {401, 1.7033818819502551},
                                                                            },
                                                                            {
                                                                                {5, 0.92705119413602599},
                                                                                {69, 1.9592598982258778},
                                                                                {133, 0.96653490574340428},
                                                                                {197, 1.3203867992383289},
-                                                                               {261, 8.0057564947791793},
+                                                                               {261, 7.9217487636977353},
                                                                                {325, 1.8629714367312684},
-                                                                               {389, 1.9034519507895218},
+                                                                               {389, 1.8587936590169301},
                                                                            }};
 
   // Loop over different directions
diff --git a/src/reconstruction/reconstruction.h b/src/reconstruction/reconstruction.h
index 07aae21a6..23442a776 100644
--- a/src/reconstruction/reconstruction.h
+++ b/src/reconstruction/reconstruction.h
@@ -673,6 +673,56 @@ Primitive __device__ __host__ __inline__ Calc_Interface_Linear(Primitive const &
 }
 // =====================================================================================================================
 
+// =====================================================================================================================
+/*!
+ * \brief Apply limiting the the primitive interfaces in PLM reconstructions
+ *
+ * \param[in,out] interface_L_iph The unlimited left plus 1/2 interface
+ * \param[in,out] interface_R_imh The unlimited right minus 1/2 interface
+ * \param[in] cell_imo The cell centered values at i-1
+ * \param[in] cell_i The cell centered values at i
+ * \param[in] cell_ipo The cell centered values at i+1
+ */
+void __device__ __host__ __inline__ Plm_Limit_Interfaces(Primitive &interface_L_iph, Primitive &interface_R_imh,
+                                                         Primitive const &cell_imo, Primitive const &cell_i,
+                                                         Primitive const &cell_ipo)
+{
+  auto limiter = [](Real &l_iph, Real &r_imh, Real const &val_imo, Real const &val_i, Real const &val_ipo) {
+    r_imh = fmax(fmin(val_i, val_imo), r_imh);
+    r_imh = fmin(fmax(val_i, val_imo), r_imh);
+    l_iph = fmax(fmin(val_i, val_ipo), l_iph);
+    l_iph = fmin(fmax(val_i, val_ipo), l_iph);
+  };
+
+  limiter(interface_L_iph.density, interface_R_imh.density, cell_imo.density, cell_i.density, cell_ipo.density);
+  limiter(interface_L_iph.velocity_x, interface_R_imh.velocity_x, cell_imo.velocity_x, cell_i.velocity_x,
+          cell_ipo.velocity_x);
+  limiter(interface_L_iph.velocity_y, interface_R_imh.velocity_y, cell_imo.velocity_y, cell_i.velocity_y,
+          cell_ipo.velocity_y);
+  limiter(interface_L_iph.velocity_z, interface_R_imh.velocity_z, cell_imo.velocity_z, cell_i.velocity_z,
+          cell_ipo.velocity_z);
+  limiter(interface_L_iph.pressure, interface_R_imh.pressure, cell_imo.pressure, cell_i.pressure, cell_ipo.pressure);
+
+#ifdef MHD
+  limiter(interface_L_iph.magnetic_y, interface_R_imh.magnetic_y, cell_imo.magnetic_y, cell_i.magnetic_y,
+          cell_ipo.magnetic_y);
+  limiter(interface_L_iph.magnetic_z, interface_R_imh.magnetic_z, cell_imo.magnetic_z, cell_i.magnetic_z,
+          cell_ipo.magnetic_z);
+#endif  // MHD
+
+#ifdef DE
+  limiter(interface_L_iph.gas_energy, interface_R_imh.gas_energy, cell_imo.gas_energy, cell_i.gas_energy,
+          cell_ipo.gas_energy);
+#endif  // DE
+#ifdef SCALAR
+  for (int i = 0; i < NSCALARS; i++) {
+    limiter(interface_L_iph.scalar[i], interface_R_imh.scalar[i], cell_imo.scalar[i], cell_i.scalar[i],
+            cell_ipo.scalar[i]);
+  }
+#endif  // SCALAR
+}
+// =====================================================================================================================
+
 // =====================================================================================================================
 /*!
  * \brief Compute the interface state for the CTU version fo the reconstructor from the slope and cell centered state
diff --git a/src/reconstruction/reconstruction_tests.cu b/src/reconstruction/reconstruction_tests.cu
index 74c0e6896..69e1a7ef6 100644
--- a/src/reconstruction/reconstruction_tests.cu
+++ b/src/reconstruction/reconstruction_tests.cu
@@ -238,7 +238,7 @@ TEST(tALLReconstructionLoadData, CorrectInputExpectCorrectOutput)
   reconstruction::Primitive fiducial_data{13, 3.0769230769230771, 5.1538461538461542, 7.2307692307692308,
                                           39950.641025641031};
   #ifdef DE
-  fiducial_data.pressure = 34274.282506448195;
+  fiducial_data.pressure = 39950.641025641031;
   #endif  // DE
   testing_utilities::Check_Results(fiducial_data.density, test_data.density, "density");
   testing_utilities::Check_Results(fiducial_data.velocity_x, test_data.velocity_x, "velocity_x");
@@ -613,3 +613,74 @@ TEST(tALLReconstructionWriteData, CorrectInputExpectCorrectOutput)
     testing_utilities::Check_Results(fiducial_val, test_val, "Interface at i=" + std::to_string(i));
   }
 }
+
+TEST(tHYDROReconstructionPlmLimitInterfaces, CorrectInputExpectCorrectOutput)
+{
+  // Set up values to test
+  reconstruction::Primitive interface_l_iph, interface_r_imh;
+  reconstruction::Primitive cell_im1, cell_i, cell_ip1;
+  interface_r_imh.density    = -1.94432878387898625e+14;
+  interface_r_imh.velocity_x = 1.42049955114756404e-04;
+  interface_r_imh.velocity_y = -2.61311412306644180e-06;
+  interface_r_imh.velocity_z = -1.99429361865204601e-07;
+  interface_r_imh.pressure   = -2.01130121665840250e-14;
+  interface_l_iph.density    = 1.94433200621991188e+14;
+  interface_l_iph.velocity_x = 1.42025407335853601e-04;
+  interface_l_iph.velocity_y = -2.61311412306644180e-06;
+  interface_l_iph.velocity_z = -6.01154878659959398e-06;
+  interface_l_iph.pressure   = 2.01130321665840277e-14;
+
+  cell_im1.density    = 1.61101072114153951e+08;
+  cell_i.density      = 1.61117046279133737e+08;
+  cell_ip1.density    = 1.61011252191243321e+08;
+  cell_im1.velocity_x = 1.42067642369120116e-04;
+  cell_i.velocity_x   = 1.42037681225305003e-04;
+  cell_ip1.velocity_x = 1.41901817571928041e-04;
+  cell_im1.velocity_y = -2.61228250783092252e-06;
+  cell_i.velocity_y   = -2.61311412306644180e-06;
+  cell_ip1.velocity_y = -2.61155204131260820e-06;
+  cell_im1.velocity_z = 2.71420653365757378e-06;
+  cell_i.velocity_z   = -3.10548907423239929e-06;
+  cell_ip1.velocity_z = -8.91005201578514336e-06;
+  cell_im1.pressure   = 9.99999999999999945e-21;
+  cell_i.pressure     = 9.99999999999999945e-21;
+  cell_ip1.pressure   = 4.70262856027679407e-03;
+
+  // Set fiducial values
+  reconstruction::Primitive interface_r_imh_fiducial, interface_l_iph_fiducial;
+  interface_r_imh_fiducial.density    = 161101072.11415395;
+  interface_r_imh_fiducial.velocity_x = 1.42049955114756404e-04;
+  interface_r_imh_fiducial.velocity_y = -2.61311412306644180e-06;
+  interface_r_imh_fiducial.velocity_z = -1.99429361865204601e-07;
+  interface_r_imh_fiducial.pressure   = 9.99999999999999945e-21;
+  interface_l_iph_fiducial.density    = 1.61117046279133737e+08;
+  interface_l_iph_fiducial.velocity_x = 1.42025407335853601e-04;
+  interface_l_iph_fiducial.velocity_y = -2.61311412306644180e-06;
+  interface_l_iph_fiducial.velocity_z = -6.01154878659959398e-06;
+  interface_l_iph_fiducial.pressure   = 2.0113032166584028e-14;
+
+  // Run function
+  reconstruction::Plm_Limit_Interfaces(interface_l_iph, interface_r_imh, cell_im1, cell_i, cell_ip1);
+
+  // Check values
+  testing_utilities::Check_Results(interface_l_iph_fiducial.density, interface_l_iph.density,
+                                   "Mismatch in l_iph density");
+  testing_utilities::Check_Results(interface_l_iph_fiducial.velocity_x, interface_l_iph.velocity_x,
+                                   "Mismatch in l_iph velocity_x");
+  testing_utilities::Check_Results(interface_l_iph_fiducial.velocity_y, interface_l_iph.velocity_y,
+                                   "Mismatch in l_iph velocity_y");
+  testing_utilities::Check_Results(interface_l_iph_fiducial.velocity_z, interface_l_iph.velocity_z,
+                                   "Mismatch in l_iph velocity_z");
+  testing_utilities::Check_Results(interface_l_iph_fiducial.pressure, interface_l_iph.pressure,
+                                   "Mismatch in l_iph pressure");
+  testing_utilities::Check_Results(interface_r_imh_fiducial.density, interface_r_imh.density,
+                                   "Mismatch in r_imh density");
+  testing_utilities::Check_Results(interface_r_imh_fiducial.velocity_x, interface_r_imh.velocity_x,
+                                   "Mismatch in r_imh velocity_x");
+  testing_utilities::Check_Results(interface_r_imh_fiducial.velocity_y, interface_r_imh.velocity_y,
+                                   "Mismatch in r_imh velocity_y");
+  testing_utilities::Check_Results(interface_r_imh_fiducial.velocity_z, interface_r_imh.velocity_z,
+                                   "Mismatch in r_imh velocity_z");
+  testing_utilities::Check_Results(interface_r_imh_fiducial.pressure, interface_r_imh.pressure,
+                                   "Mismatch in r_imh pressure");
+}
diff --git a/src/system_tests/mhd_system_tests.cpp b/src/system_tests/mhd_system_tests.cpp
index c7a21aaae..4261797b2 100644
--- a/src/system_tests/mhd_system_tests.cpp
+++ b/src/system_tests/mhd_system_tests.cpp
@@ -765,7 +765,7 @@ TEST_P(tMHDSYSTEMParameterizedMpi, AdvectingFieldLoopCorrectInputExpectCorrectOu
   test_runner.numMpiRanks = GetParam();
 
   // Only do the L2 Norm test. The regular cell-to-cell comparison is brittle for this test across systems
-  test_runner.runTest(true, 3.9E-8, 1.6E-6);
+  test_runner.runTest(true, 3.9E-8, 2.25E-6);
 }
 
 /// Test the MHD Blast Wave
diff --git a/src/utils/hydro_utilities.h b/src/utils/hydro_utilities.h
index c0f783e1c..24caff9f7 100644
--- a/src/utils/hydro_utilities.h
+++ b/src/utils/hydro_utilities.h
@@ -35,7 +35,7 @@ inline __host__ __device__ Real Calc_Pressure_Primitive(Real const &E, Real cons
                                                         Real const &vz, Real const &gamma, Real const &magnetic_x = 0.0,
                                                         Real const &magnetic_y = 0.0, Real const &magnetic_z = 0.0)
 {
-  Real pressure = (E - 0.5 * d * (vx * vx + ((vy * vy) + (vz * vz))));
+  Real pressure = E - 0.5 * d * math_utils::SquareMagnitude(vx, vy, vz);
 
 #ifdef MHD
   pressure -= mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z);
@@ -48,7 +48,7 @@ inline __host__ __device__ Real Calc_Pressure_Conserved(Real const &E, Real cons
                                                         Real const &mz, Real const &gamma, Real const &magnetic_x = 0.0,
                                                         Real const &magnetic_y = 0.0, Real const &magnetic_z = 0.0)
 {
-  Real pressure = (E - 0.5 * (mx * mx + my * my + mz * mz) / d);
+  Real pressure = E - 0.5 * math_utils::SquareMagnitude(mx, my, mz) / d;
 
 #ifdef MHD
   pressure -= mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z);
@@ -76,7 +76,7 @@ inline __host__ __device__ Real Calc_Energy_Primitive(Real const &P, Real const
                                                       Real const &magnetic_y = 0.0, Real const &magnetic_z = 0.0)
 {
   // Compute and return energy
-  Real energy = (fmax(P, TINY_NUMBER) / (gamma - 1.)) + 0.5 * d * (vx * vx + vy * vy + vz * vz);
+  Real energy = (fmax(P, TINY_NUMBER) / (gamma - 1.)) + 0.5 * d * math_utils::SquareMagnitude(vx, vy, vz);
 
 #ifdef MHD
   energy += mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z);
@@ -92,7 +92,7 @@ inline __host__ __device__ Real Calc_Energy_Conserved(Real const &P, Real const
 {
   // Compute and return energy
   Real energy = (fmax(P, TINY_NUMBER) / (gamma - 1.)) +
-                (0.5 / d) * (momentum_x * momentum_x + momentum_y * momentum_y + momentum_z * momentum_z);
+                (0.5 / d) * math_utils::SquareMagnitude(momentum_x, momentum_y, momentum_z);
 
 #ifdef MHD
   energy += mhd::utils::computeMagneticEnergy(magnetic_x, magnetic_y, magnetic_z);
@@ -130,7 +130,7 @@ inline __host__ __device__ Real Get_Pressure_From_DE(Real const &E, Real const &
 inline __host__ __device__ Real Calc_Kinetic_Energy_From_Velocity(Real const &d, Real const &vx, Real const &vy,
                                                                   Real const &vz)
 {
-  return 0.5 * d * (vx * vx + vy * vy * vz * vz);
+  return 0.5 * d * math_utils::SquareMagnitude(vx, vy, vz);
 }
 
 /*!
@@ -145,7 +145,7 @@ inline __host__ __device__ Real Calc_Kinetic_Energy_From_Velocity(Real const &d,
 inline __host__ __device__ Real Calc_Kinetic_Energy_From_Momentum(Real const &d, Real const &mx, Real const &my,
                                                                   Real const &mz)
 {
-  return (0.5 / d) * (mx * mx + my * my * mz * mz);
+  return (0.5 / d) * math_utils::SquareMagnitude(mx, my, mz);
 }
 
 /*!
diff --git a/src/utils/hydro_utilities_tests.cpp b/src/utils/hydro_utilities_tests.cpp
index eda204c76..b200ddd8c 100644
--- a/src/utils/hydro_utilities_tests.cpp
+++ b/src/utils/hydro_utilities_tests.cpp
@@ -238,7 +238,7 @@ TEST(tHYDROHydroUtilsGetPressureFromDE, CorrectInputExpectCorrectOutput)
 TEST(tHYDROtMHDCalcKineticEnergyFromVelocity, CorrectInputExpectCorrectOutput)
 {
   TestParams parameters;
-  std::vector<double> fiducialEnergies{0.0, 6.307524975350106e-145, 7.3762470327090601e+249};
+  std::vector<double> fiducialEnergies{0.0, 6.307524975350106e-145, 1.9018677140549924e+150};
   double const coef = 1E-50;
 
   for (size_t i = 0; i < parameters.names.size(); i++) {
@@ -252,7 +252,7 @@ TEST(tHYDROtMHDCalcKineticEnergyFromVelocity, CorrectInputExpectCorrectOutput)
 TEST(tHYDROtMHDCalcKineticEnergyFromMomentum, CorrectInputExpectCorrectOutput)
 {
   TestParams parameters;
-  std::vector<double> fiducialEnergies{0.0, 0.0, 7.2568536478335773e+147};
+  std::vector<double> fiducialEnergies{0.0, 0.0, 3.0042157852278499e+49};
   double const coef = 1E-50;
 
   for (size_t i = 0; i < parameters.names.size(); i++) {
diff --git a/src/utils/math_utilities.h b/src/utils/math_utilities.h
index 1480f852c..68d13f19d 100644
--- a/src/utils/math_utilities.h
+++ b/src/utils/math_utilities.h
@@ -82,4 +82,20 @@ inline __device__ __host__ Real dotProduct(Real const &a1, Real const &a2, Real
 };
 // =========================================================================
 
+// =========================================================================
+/*!
+ * \brief Compute the magnitude of a vector
+ *
+ * \param[in] v1 The first element of the vector
+ * \param[in] v2 The second element of the vector
+ * \param[in] v3 The third element of the vector
+ *
+ * \return Real The dot product of a and b
+ */
+inline __device__ __host__ Real SquareMagnitude(Real const &v1, Real const &v2, Real const &v3)
+{
+  return dotProduct(v1, v2, v3, v1, v2, v3);
+};
+// =========================================================================
+
 }  // namespace math_utils
diff --git a/src/utils/math_utilities_tests.cpp b/src/utils/math_utilities_tests.cpp
index 83fd1d232..a49cd8a41 100644
--- a/src/utils/math_utilities_tests.cpp
+++ b/src/utils/math_utilities_tests.cpp
@@ -56,4 +56,22 @@ TEST(tALLDotProduct, CorrectInputExpectCorrectOutput)
   // Now check results
   testing_utilities::Check_Results(fiducialDotProduct, testDotProduct, "dot product");
 }
+// =========================================================================
+
+// =========================================================================
+/*!
+ * \brief Test the math_utils::dotProduct function
+ *
+ */
+TEST(tALLSquareMagnitude, CorrectInputExpectCorrectOutput)
+{
+  std::vector<double> a = {11.503067766457753, 98.316634031589935, 41.12177317622657};
+
+  double const fiducial_square_magnitude = 11489.481324498336;
+
+  double test_square_magnitude = math_utils::SquareMagnitude(a.at(0), a.at(1), a.at(2));
+
+  // Now check results
+  testing_utilities::Check_Results(fiducial_square_magnitude, test_square_magnitude, "dot product");
+}
 // =========================================================================
\ No newline at end of file
diff --git a/src/utils/mhd_utilities.h b/src/utils/mhd_utilities.h
index 55ecc6f75..f409fd4b0 100644
--- a/src/utils/mhd_utilities.h
+++ b/src/utils/mhd_utilities.h
@@ -18,6 +18,7 @@
 #include "../grid/grid3D.h"
 #include "../utils/cuda_utilities.h"
 #include "../utils/gpu.hpp"
+#include "../utils/math_utilities.h"
 
 namespace mhd::utils
 {
@@ -74,7 +75,7 @@ inline __host__ __device__ Real _magnetosonicSpeed(Real const &density, Real con
 inline __host__ __device__ Real computeMagneticEnergy(Real const &magneticX, Real const &magneticY,
                                                       Real const &magneticZ)
 {
-  return 0.5 * (magneticX * magneticX + ((magneticY * magneticY) + (magneticZ * magneticZ)));
+  return 0.5 * math_utils::SquareMagnitude(magneticX, magneticY, magneticZ);
 }
 // =========================================================================
 
@@ -98,9 +99,7 @@ inline __host__ __device__ Real computeThermalEnergy(Real const &energyTot, Real
                                                      Real const &magneticX, Real const &magneticY,
                                                      Real const &magneticZ, Real const &gamma)
 {
-  return energyTot -
-         0.5 * (momentumX * momentumX + ((momentumY * momentumY) + (momentumZ * momentumZ))) /
-             fmax(density, TINY_NUMBER) -
+  return energyTot - 0.5 * math_utils::SquareMagnitude(momentumX, momentumY, momentumZ) / fmax(density, TINY_NUMBER) -
          computeMagneticEnergy(magneticX, magneticY, magneticZ);
 }
 // =========================================================================