diff --git a/CMakeLists.txt b/CMakeLists.txt
index 51edffe..d505d9d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -34,6 +34,7 @@ if (${CMAKE_BUILD_TYPE} STREQUAL "Debug")
      "Bounds checking on Kokkos views")
   set(Kokkos_ENABLE_DEBUG_DUALVIEW_MODIFY_CHECK ON CACHE BOOL
       "Sanity checks on Kokkos DualView")
+    add_compile_definitions(ATHELAS_DEBUG)
 endif()
 
 set(Kokkos_ENABLE_AGGRESSIVE_VECTORIZATION ON CACHE BOOL
diff --git a/src/basis/polynomial_basis.hpp b/src/basis/polynomial_basis.hpp
index 05b73a0..a7326c3 100644
--- a/src/basis/polynomial_basis.hpp
+++ b/src/basis/polynomial_basis.hpp
@@ -1,5 +1,5 @@
-#ifndef _POLYNOMIAL_BASIS_HPP_
-#define _POLYNOMIAL_BASIS_HPP_
+#ifndef POLYNOMIAL_BASIS_HPP_
+#define POLYNOMIAL_BASIS_HPP_
 
 /**
  * File     :  polynomial_basis.hpp
@@ -72,4 +72,4 @@ class ModalBasis {
   Real ( *dfunc )( int n, Real x, Real x_c );
 };
 
-#endif // _POLYNOMIAL_BASIS_HPP_
+#endif // POLYNOMIAL_BASIS_HPP_
diff --git a/src/bc/boundary_conditions.cpp b/src/bc/boundary_conditions.cpp
index 176843f..ecf0223 100644
--- a/src/bc/boundary_conditions.cpp
+++ b/src/bc/boundary_conditions.cpp
@@ -13,8 +13,16 @@
 
 #include "boundary_conditions.hpp"
 #include "grid.hpp"
+#include "utilities.hpp"
 
-// Apply Boundary Conditions to fluid fields
+/**
+ * Apply Boundary Conditions to fluid fields
+ * Supported Options:
+ *  reflecting
+ *  periodic
+ *  shockless_noh
+ *  homogenous (Default)
+ **/
 void ApplyBC( View3D<Real> U, GridStructure *Grid, const int order,
               const std::string BC ) {
 
@@ -26,7 +34,7 @@ void ApplyBC( View3D<Real> U, GridStructure *Grid, const int order,
   const int nG = Grid->Get_Guard( );
 
   // ! ? How to correctly implement reflecting BC ? !
-  if ( BC == "Reflecting" ) {
+  if ( utilities::to_lower( BC ) == "reflecting" ) {
     for ( int iCF = 0; iCF < nvars; iCF++ ) {
       // Inner Boudnary
       for ( int iX = 0; iX < ilo; iX++ )
@@ -50,14 +58,15 @@ void ApplyBC( View3D<Real> U, GridStructure *Grid, const int order,
           }
         }
     }
-  } else if ( BC == "Periodic" ) {
+  } else if ( utilities::to_lower( BC ) == "periodic" ) {
     for ( int iCF = 0; iCF < nvars; iCF++ )
       for ( int iX = 0; iX < ilo; iX++ )
         for ( int k = 0; k < order; k++ ) {
           U( iCF, ilo - 1 - iX, k ) = U( iCF, ihi - iX, k );
           U( iCF, ihi + 1 + iX, k ) = U( iCF, ilo + iX, k );
         }
-  } else if ( BC == "ShocklessNoh" ) /* Special case for ShocklessNoh test */
+  } else if ( utilities::to_lower( BC ) ==
+              "shockless_noh" ) /* Special case for ShocklessNoh test */
   {
     // for ( int iCF = 0; iCF < 3; iCF++ )
     for ( int iX = 0; iX < ilo; iX++ )
@@ -121,12 +130,14 @@ void ApplyBC( View3D<Real> U, GridStructure *Grid, const int order,
           U( 2, ihi + 1 + iX, k ) = U( 2, ihi - iX, k );
         }
       }
-  } else {
+  } else if ( utilities::to_lower( BC ) == "homogenous" ) {
     for ( int iCF = 0; iCF < nvars; iCF++ )
       for ( int iX = 0; iX < ilo; iX++ )
         for ( int k = 0; k < order; k++ ) {
           U( iCF, ilo - 1 - iX, k ) = U( iCF, ilo + iX, k );
           U( iCF, ihi + 1 + iX, k ) = U( iCF, ihi - iX, k );
         }
+  } else {
+    throw Error( " ! Error: BC not supported!" );
   }
 }
diff --git a/src/driver.cpp b/src/driver.cpp
index 7cb0bfd..81ca907 100644
--- a/src/driver.cpp
+++ b/src/driver.cpp
@@ -139,6 +139,10 @@ int main( int argc, char *argv[] ) {
                            &Basis, &eos, &S_Limiter, opts );
       }
 
+#ifdef ATHELAS_DEBUG
+      check_state( &state, Grid.Get_ihi( ) );
+#endif
+
       t += dt;
 
       // Write state
diff --git a/src/eos/eos.hpp b/src/eos/eos.hpp
index 2562a53..e8d71ee 100644
--- a/src/eos/eos.hpp
+++ b/src/eos/eos.hpp
@@ -9,11 +9,14 @@
 
 #include "abstractions.hpp"
 #include "eos_base.hpp"
+#include "error.hpp"
 
 class IdealGas : public EosBase<IdealGas> {
  public:
   IdealGas( ) = default;
-  IdealGas( double gm ) : gamma( gm ) {}
+  IdealGas( double gm ) : gamma( gm ) {
+    if ( gamma <= 0.0 ) throw Error( " ! Adiabatic gamma <= 0.0!" );
+  }
 
   Real PressureFromConserved( const Real Tau, const Real V, const Real EmT,
                               Real *lambda ) const;
diff --git a/src/fluid/fluid_discretization.cpp b/src/fluid/fluid_discretization.cpp
index cc63730..85ea481 100644
--- a/src/fluid/fluid_discretization.cpp
+++ b/src/fluid/fluid_discretization.cpp
@@ -31,13 +31,11 @@ void ComputeIncrement_Fluid_Divergence(
   const auto &ihi    = Grid.Get_ihi( );
   const int nvars    = U.extent( 0 );
 
-  // Real rho_L, rho_R, P_L, P_R, Cs_L, Cs_R, lam_L, lam_R, P;
-
   // --- Interpolate Conserved Variable to Interfaces ---
 
   // Left/Right face states
   Kokkos::parallel_for(
-      "Interface States",
+      "Fluid :: Interface States",
       Kokkos::MDRangePolicy<Kokkos::Rank<2>>( { ilo, 0 }, { ihi + 2, nvars } ),
       KOKKOS_LAMBDA( const int iX, const int iCF ) {
         uCF_F_L( iCF, iX ) = Basis->BasisEval( U, iX - 1, iCF, nNodes + 1 );
@@ -46,7 +44,7 @@ void ComputeIncrement_Fluid_Divergence(
 
   // --- Calc numerical flux at all faces
   Kokkos::parallel_for(
-      "Numerical Fluxes", Kokkos::RangePolicy<>( ilo, ihi + 2 ),
+      "Fluid :: Numerical Fluxes", Kokkos::RangePolicy<>( ilo, ihi + 2 ),
       KOKKOS_LAMBDA( int iX ) {
         auto uCF_L = Kokkos::subview( uCF_F_L, Kokkos::ALL, iX );
         auto uCF_R = Kokkos::subview( uCF_F_R, Kokkos::ALL, iX );
@@ -54,6 +52,16 @@ void ComputeIncrement_Fluid_Divergence(
         const Real rho_L = 1.0 / uCF_L( 0 );
         const Real rho_R = 1.0 / uCF_R( 0 );
 
+        // Debug mode assertions.
+        assert( rho_L > 0.0 && !std::isnan( rho_L ) &&
+                "fluid_discretization :: Numerical Fluxes bad rho." );
+        assert( rho_R > 0.0 && !std::isnan( rho_R ) &&
+                "fluid_discretization :: Numerical Fluxes bad rho." );
+        assert( uCF_L( 2 ) > 0.0 && !std::isnan( uCF_L( 2 ) ) &&
+                "fluid_discretization :: Numerical Fluxes bad energy." );
+        assert( uCF_R( 2 ) > 0.0 && !std::isnan( uCF_R( 2 ) ) &&
+                "fluid_discretization :: Numerical Fluxes bad energy." );
+
         auto lambda      = nullptr;
         const Real P_L   = eos->PressureFromConserved( uCF_L( 0 ), uCF_L( 1 ),
                                                        uCF_L( 2 ), lambda );
@@ -85,7 +93,7 @@ void ComputeIncrement_Fluid_Divergence(
 
   // --- Surface Term ---
   Kokkos::parallel_for(
-      "Surface Term",
+      "Fluid :: Surface Term",
       Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                               { order, ihi + 1, nvars } ),
       KOKKOS_LAMBDA( const int k, const int iX, const int iCF ) {
@@ -103,7 +111,7 @@ void ComputeIncrement_Fluid_Divergence(
   if ( order > 1 ) {
     // --- Compute Flux_q everywhere for the Volume term ---
     Kokkos::parallel_for(
-        "Flux_q",
+        "Fluid :: Flux_q",
         Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                                 { nNodes, ihi + 1, nvars } ),
         KOKKOS_LAMBDA( const int iN, const int iX, const int iCF ) {
@@ -119,7 +127,7 @@ void ComputeIncrement_Fluid_Divergence(
     // --- Volume Term ---
     // TODO: Make Flux_q a function?
     Kokkos::parallel_for(
-        "Volume Term",
+        "Fluid :: Volume Term",
         Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                                 { order, ihi + 1, nvars } ),
         KOKKOS_LAMBDA( const int k, const int iX, const int iCF ) {
@@ -148,7 +156,7 @@ void ComputeIncrement_Fluid_Geometry( const View3D<Real> U, GridStructure &Grid,
   const int ihi    = Grid.Get_ihi( );
 
   Kokkos::parallel_for(
-      "Geometry Term; Fluid",
+      "Fluid :: Geometry Term",
       Kokkos::MDRangePolicy<Kokkos::Rank<2>>( { 0, ilo }, { order, ihi + 1 } ),
       KOKKOS_LAMBDA( const int k, const int iX ) {
         Real local_sum = 0.0;
@@ -183,7 +191,7 @@ void ComputeIncrement_Fluid_Rad( const View3D<Real> uCF, const View3D<Real> uCR,
   const int ihi    = Grid.Get_ihi( );
 
   Kokkos::parallel_for(
-      "Fluid Source Term; Rad",
+      "Fluid :: Source Term; Rad",
       Kokkos::MDRangePolicy<Kokkos::Rank<2>>( { 0, ilo }, { order, ihi + 1 } ),
       KOKKOS_LAMBDA( const int k, const int iX ) {
         Real local_sum1 = 0.0;
diff --git a/src/fluid/fluid_utilities.cpp b/src/fluid/fluid_utilities.cpp
index 1609855..7ced840 100644
--- a/src/fluid/fluid_utilities.cpp
+++ b/src/fluid/fluid_utilities.cpp
@@ -55,6 +55,8 @@ void ComputePrimitiveFromConserved( View3D<Real> uCF, View3D<Real> uPF,
  * TODO: Flux_Fluid needs streamlining
  **/
 Real Flux_Fluid( const Real V, const Real P, const int iCF ) {
+  assert( iCF == 0 || iCF == 1 || iCF == 2 );
+  assert( P > 0.0 && "Flux_Flux :: negative pressure" );
   if ( iCF == 0 ) {
     return -V;
   } else if ( iCF == 1 ) {
@@ -89,6 +91,8 @@ Real Source_Fluid_Rad( Real D, Real V, Real T, Real X, Real kappa, Real E,
 void NumericalFlux_Gudonov( const Real vL, const Real vR, const Real pL,
                             const Real pR, const Real zL, const Real zR,
                             Real &Flux_U, Real &Flux_P ) {
+  assert( pL > 0.0 && pR > 0.0 &&
+          "NumericalFlux_Gudonov :: negative pressure" );
   Flux_U = ( pL - pR + zR * vR + zL * vL ) / ( zR + zL );
   Flux_P = ( zR * pL + zL * pR + zL * zR * ( vL - vR ) ) / ( zR + zL );
 }
@@ -128,6 +132,9 @@ Real ComputeTimestep_Fluid( const View3D<Real> U, const GridStructure *Grid,
         Real vel_x  = U( 1, iX, 0 );
         Real eint_x = U( 2, iX, 0 );
 
+        assert( tau_x > 0.0 && "Compute Timestep :: bad specific volume" );
+        assert( eint_x > 0.0 && "Compute Timestep :: bad specific energy" );
+
         Real dr = Grid->Get_Widths( iX );
 
         auto lambda = nullptr;
@@ -144,8 +151,8 @@ Real ComputeTimestep_Fluid( const View3D<Real> U, const GridStructure *Grid,
   dt = std::max( CFL * dt, MIN_DT );
   dt = std::min( dt, MAX_DT );
 
-  // Triggers on NaN
-  if ( dt != dt ) {
+  // could be assertion?
+  if ( std::isnan( dt ) ) {
     throw Error( " ! nan encountered in ComputeTimestep.\n" );
   }
 
diff --git a/src/geometry/geometry.hpp b/src/geometry/geometry.hpp
index ee56c6b..303d903 100644
--- a/src/geometry/geometry.hpp
+++ b/src/geometry/geometry.hpp
@@ -1,5 +1,5 @@
-#ifndef _GEOMETRY_HPP_
-#define _GEOMETRY_HPP_
+#ifndef GEOMETRY_HPP_
+#define GEOMETRY_HPP_
 
 /**
  * File     :  geometry.hpp
@@ -13,4 +13,4 @@
 namespace geometry {
 enum Geometry { Planar, Spherical };
 }
-#endif // _GEOMETRY_HPP_
+#endif // GEOMETRY_HPP_
diff --git a/src/geometry/grid.hpp b/src/geometry/grid.hpp
index c04e895..93585bd 100644
--- a/src/geometry/grid.hpp
+++ b/src/geometry/grid.hpp
@@ -1,5 +1,5 @@
-#ifndef _GRID_HPP_
-#define _GRID_HPP_
+#ifndef GRID_HPP_
+#define GRID_HPP_
 
 /**
  * File     :  grid.hpp
@@ -81,4 +81,4 @@ class GridStructure {
   Kokkos::View<Real **> Grid;
 };
 
-#endif // _GRID_HPP_
+#endif // GRID_HPP_
diff --git a/src/pgen/problem_in.cpp b/src/pgen/problem_in.cpp
index 9f04a7d..bb09d0a 100644
--- a/src/pgen/problem_in.cpp
+++ b/src/pgen/problem_in.cpp
@@ -24,6 +24,8 @@ ProblemIn::ProblemIn( const std::string fn ) {
     throw Error( " ! Issue reading input deck!" );
   }
 
+  std::printf( " ~ Loading Input Deck ...\n" );
+
   /* Grab as std::optional<type> */
 
   // problem
@@ -65,19 +67,34 @@ ProblemIn::ProblemIn( const std::string fn ) {
   if ( pn ) {
     ProblemName = pn.value( );
   } else {
-    throw Error( " ! Error: problem not supplied in input deck." );
+    throw Error(
+        " ! Initialization Error: problem not supplied in input deck." );
   }
+  // Validity of ProblemName checked in initialization.
+
   if ( bc ) {
-    BC = bc.value( );
+    BC = utilities::to_lower( bc.value( ) );
   } else {
-    throw Error( " ! Error: boundary condition not supplied in input deck." );
+    throw Error( " ! Initialization Error: boundary condition not supplied in "
+                 "input deck." );
+  }
+  if ( BC != "homogenous" && BC != "reflecting" && BC != "shockless_noh" &&
+       BC != "periodic" ) {
+    throw Error(
+        " ! Initialization Error: Bad boundary condition choice. Choose: \n"
+        " - homogenous \n"
+        " - reflecting \n"
+        " - periodic \n"
+        " - shockless_noh" );
   }
+
   if ( geom ) {
     Geometry = ( utilities::to_lower( geom.value( ) ) == "spherical" )
                    ? geometry::Spherical
                    : geometry::Planar;
   } else {
-    Geometry = geometry::Planar;
+    std::printf( "   - Defaulting to planar geometry!\n" );
+    Geometry = geometry::Planar; // default
   }
   if ( basis ) {
     Basis = ( utilities::to_lower( basis.value( ) ) == "legendre" )
@@ -85,30 +102,36 @@ ProblemIn::ProblemIn( const std::string fn ) {
                 : PolyBasis::Taylor;
   } else {
     Basis = PolyBasis::Legendre;
+    std::printf( "   - Defaulting to Legendre polynomial basis!\n" );
   }
 
   if ( x1 ) {
     xL = x1.value( );
   } else {
-    throw Error( " ! Error: xL not supplied in input deck." );
+    throw Error( " ! Initialization Error: xL not supplied in input deck." );
   }
   if ( x2 ) {
     xR = x2.value( );
   } else {
-    throw Error( " ! Error: xR not supplied in input deck." );
+    throw Error( " ! Initialization Error: xR not supplied in input deck." );
   }
+  if ( x1 >= x2 ) throw Error( " ! Initialization Error: x1 >= xz2" );
+
   if ( tf ) {
     t_end = tf.value( );
   } else {
-    throw Error( " ! Error: t_end not supplied in input deck." );
+    throw Error( " ! Initialization Error: t_end not supplied in input deck." );
   }
-  CFL = cfl.value_or( 0.5 );
+  if ( tf <= 0.0 ) throw Error( " ! Initialization Error: tf <= 0.0" );
+
   if ( nX ) {
     nElements = nX.value( );
   } else {
-    throw Error( " ! Error: nX not supplied in innput deck." );
+    throw Error( " ! Initialization Error: nX not supplied in innput deck." );
   }
 
+  // many defaults not mentioned below...
+  CFL     = cfl.value_or( 0.5 );
   Restart = rest.value_or( false );
   do_rad  = rad.value_or( false );
 
@@ -124,6 +147,17 @@ ProblemIn::ProblemIn( const std::string fn ) {
   gamma_l        = gamma1.value_or( 0.005 );
   gamma_i        = gamma2.value_or( 0.990 );
   gamma_r        = gamma3.value_or( 0.005 );
+
+  // varous checks
+  if ( CFL <= 0.0 ) throw Error( " ! Initialization : CFL <= 0.0!" );
+  if ( nGhost <= 0 ) throw Error( " ! Initialization : nGhost <= 0!" );
+  if ( pOrder <= 0 ) throw Error( " ! Initialization : pOrder <= 0!" );
+  if ( nNodes <= 0 ) throw Error( " ! Initialization : nNodes <= 0!" );
+  if ( tOrder <= 0 ) throw Error( " ! Initialization : tOrder <= 0!" );
+  if ( nStages <= 0 ) throw Error( " ! Initialization : nStages <= 0!" );
+  if ( TCI_Threshold <= 0.0 )
+    throw Error( " ! Initialization : TCI_Threshold <= 0.0!" );
+
   if ( ( gamma2 && !gamma1 ) || ( gamma2 && !gamma3 ) ) {
     gamma_i = gamma2.value( );
     gamma_l = ( 1.0 - gamma_i ) / 2.0;
@@ -133,7 +167,9 @@ ProblemIn::ProblemIn( const std::string fn ) {
   if ( std::fabs( sum_g - 1.0 ) > 1.0e-10 ) {
     std::fprintf( stderr, "{gamma}, sum gamma = { %.10f %.10f %.10f }, %.18e\n",
                   gamma_l, gamma_i, gamma_r, 1.0 - sum_g );
-    throw Error( " ! Linear weights must sum to unity." );
+    throw Error( " ! Initialization Error: Linear weights must sum to unity." );
   }
   weno_r = wenor.value_or( 2.0 );
+
+  std::printf( " ~ Configuration ... Complete\n\n" );
 }
diff --git a/src/radiation/rad_discretization.cpp b/src/radiation/rad_discretization.cpp
index 6a84550..dbbe9d8 100644
--- a/src/radiation/rad_discretization.cpp
+++ b/src/radiation/rad_discretization.cpp
@@ -32,14 +32,12 @@ void ComputeIncrement_Rad_Divergence(
   const auto &ihi    = Grid.Get_ihi( );
   const int nvars    = uCR.extent( 0 );
 
-  // Real rho_L, rho_R, P_L, P_R, Cs_L, Cs_R, lam_L, lam_R, P;
-
   // --- Interpolate Conserved Variable to Interfaces ---
 
   // Left/Right face states
   // TODO: Can this just be moved into the below kernel with a iCF loop?
   Kokkos::parallel_for(
-      "Interface States; Rad",
+      "Radiation :: Interface States",
       Kokkos::MDRangePolicy<Kokkos::Rank<2>>( { ilo, 0 }, { ihi + 2, nvars } ),
       KOKKOS_LAMBDA( const int iX, const int iCR ) {
         uCR_F_L( iCR, iX ) = Basis->BasisEval( uCR, iX - 1, iCR, nNodes + 1 );
@@ -48,7 +46,7 @@ void ComputeIncrement_Rad_Divergence(
 
   // --- Calc numerical flux at all faces
   Kokkos::parallel_for(
-      "Numerical Fluxes; Rad", Kokkos::RangePolicy<>( ilo, ihi + 2 ),
+      "Radiation :: Numerical Fluxes", Kokkos::RangePolicy<>( ilo, ihi + 2 ),
       KOKKOS_LAMBDA( int iX ) {
         auto uCR_L = Kokkos::subview( uCR_F_L, Kokkos::ALL, iX );
         auto uCR_R = Kokkos::subview( uCR_F_R, Kokkos::ALL, iX );
@@ -56,6 +54,20 @@ void ComputeIncrement_Rad_Divergence(
         const Real tauR = Basis->BasisEval( uCF, iX, 0, 0 );
         const Real tauL = Basis->BasisEval( uCF, iX - 1, 0, nNodes + 1 );
 
+        // Debug mode assertions.
+        assert( tauL > 0.0 && !std::isnan( tauL ) &&
+                "rad_discretization :: Numerical Fluxes bad specific volume." );
+        assert( tauR > 0.0 && !std::isnan( tauR ) &&
+                "rad_discretization :: Numerical Fluxes bad specific volume." );
+        assert( uCR_L( 0 ) > 0.0 && !std::isnan( uCR_L( 0 ) ) &&
+                "rad_Discretization :: Numerical Fluxes bad energy." );
+        assert( uCR_R( 0 ) > 0.0 && !std::isnan( uCR_R( 0 ) ) &&
+                "rad_Discretization :: Numerical Fluxes bad energy." );
+        assert( !std::isnan( uCR_L( 1 ) ) &&
+                "rad_Discretization :: Numerical Fluxes bad flux." );
+        assert( !std::isnan( uCR_R( 1 ) ) &&
+                "rad_Discretization :: Numerical Fluxes bad flux." );
+
         const Real Em_L = uCR_L( 0 ) / tauL;
         const Real Fm_L = uCR_L( 1 ) / tauL;
         const Real Em_R = uCR_R( 0 ) / tauR;
@@ -87,7 +99,7 @@ void ComputeIncrement_Rad_Divergence(
 
   // --- Surface Term ---
   Kokkos::parallel_for(
-      "Surface Term; Rad",
+      "Radiation :: Surface Term",
       Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                               { order, ihi + 1, nvars } ),
       KOKKOS_LAMBDA( const int k, const int iX, const int iCR ) {
@@ -105,7 +117,7 @@ void ComputeIncrement_Rad_Divergence(
   if ( order > 1 ) {
     // --- Compute Flux_q everywhere for the Volume term ---
     Kokkos::parallel_for(
-        "Flux_q; Rad",
+        "Radiation :: Flux_q",
         Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                                 { nNodes, ihi + 1, nvars } ),
         KOKKOS_LAMBDA( const int iN, const int iX, const int iCR ) {
@@ -122,7 +134,7 @@ void ComputeIncrement_Rad_Divergence(
     // --- Volume Term ---
     // TODO: Make Flux_q a function?
     Kokkos::parallel_for(
-        "Volume Term; Rad",
+        "Radiation :: Volume Term",
         Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                                 { order, ihi + 1, nvars } ),
         KOKKOS_LAMBDA( const int k, const int iX, const int iCR ) {
@@ -153,7 +165,7 @@ void ComputeIncrement_Rad_Source( const View3D<Real> uCR,
   const int nvars  = uCR.extent( 0 );
 
   Kokkos::parallel_for(
-      "Rad::Source",
+      "Rad :: Source",
       Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                               { order, ihi + 1, nvars } ),
       KOKKOS_LAMBDA( const int k, const int iX, const int iCR ) {
@@ -221,7 +233,7 @@ void Compute_Increment_Explicit_Rad(
 
   // --- First: Zero out dU  ---
   Kokkos::parallel_for(
-      "Rad::Zero dU",
+      "Rad :: Zero dU",
       Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, 0, 0 },
                                               { order, ihi + 1, nvars } ),
       KOKKOS_LAMBDA( const int k, const int iX, const int iCR ) {
@@ -238,7 +250,7 @@ void Compute_Increment_Explicit_Rad(
 
   // --- Divide update by mass mastrix ---
   Kokkos::parallel_for(
-      "Rad::Divide Update / Mass Matrix",
+      "Rad :: Divide Update / Mass Matrix",
       Kokkos::MDRangePolicy<Kokkos::Rank<3>>( { 0, ilo, 0 },
                                               { order, ihi + 1, nvars } ),
       KOKKOS_LAMBDA( const int k, const int iX, const int iCR ) {
diff --git a/src/radiation/rad_utilities.cpp b/src/radiation/rad_utilities.cpp
index 18861c7..6d3e129 100644
--- a/src/radiation/rad_utilities.cpp
+++ b/src/radiation/rad_utilities.cpp
@@ -22,6 +22,8 @@
  * radiation flux factor
  **/
 Real FluxFactor( const Real E, const Real F ) {
+  assert( E > 0.0 &&
+          "Radiation :: FluxFactor :: non positive definite energy density." );
   const Real c = constants::c_cgs;
   return std::fabs( F ) / ( c * E );
 }
@@ -31,7 +33,9 @@ Real FluxFactor( const Real E, const Real F ) {
  * Here E and F are per unit volume
  **/
 Real Flux_Rad( Real E, Real F, Real P, Real V, int iCR ) {
-  assert( iCR == 0 || iCR == 1 );
+  assert( iCR == 0 || iCR == 1 && "Radiation :: FluxFactor :: bad iCR." );
+  assert( E > 0.0 &&
+          "Radiation :: FluxFactor :: non positive definite energy density." );
 
   const Real c = constants::c_cgs;
   return ( iCR == 0 ) ? c * F - E * V : 1.0 * c * P - F * V;
@@ -39,10 +43,27 @@ Real Flux_Rad( Real E, Real F, Real P, Real V, int iCR ) {
 
 /**
  * Radiation 4 force for rad-matter interactions
+ * D : Density
+ * V : Velocity
+ * T : Temperature
+ * kappa : kappa
+ * E : radiation energy density
+ * F : radiation momentum density
+ * Pr : radiation momentum closure
  * TODO: total opacity X?
  **/
 void RadiationFourForce( Real D, Real V, Real T, Real kappa, Real E, Real F,
                          Real Pr, Real &G0, Real &G ) {
+  assert(
+      D >= 0.0 &&
+      "Radiation :: RadiationFourFource :: Non positive definite density." );
+  assert( T > 0.0 &&
+          "Radiation :: RadiationFourFource :: Non positive temperature." );
+  assert( kappa > 0.0 &&
+          "Radiation :: RadiationFourFource :: Non positive opacity." );
+  assert( E > 0.0 && "Radiation :: RadiationFourFource :: Non positive "
+                     "definite radiation energy density." );
+
   const Real a = constants::a;
   const Real c = constants::c_cgs;
 
@@ -66,7 +87,14 @@ void RadiationFourForce( Real D, Real V, Real T, Real kappa, Real E, Real F,
  **/
 Real Source_Rad( Real D, Real V, Real T, Real X, Real kappa, Real E, Real F,
                  Real Pr, int iCR ) {
-  assert( iCR == 0 || iCR == 1 );
+  assert( iCR == 0 || iCR == 1 && "Radiation :: FluxFactor :: bad iCR." );
+  assert(
+      D >= 0.0 &&
+      "Radiation :: RadiationFourFource :: Non positive definite density." );
+  assert( T > 0.0 &&
+          "Radiation :: RadiationFourFource :: Non positive temperature." );
+  assert( E > 0.0 && "Radiation :: RadiationFourFource :: Non positive "
+                     "definite radiation energy density." );
 
   const Real c = constants::c_cgs;
 
@@ -95,6 +123,8 @@ Real ComputeOpacity( const Real D, const Real V, const Real Em ) {
 /* pressure tensor closure */
 // TODO: check Closure
 Real ComputeClosure( const Real E, const Real F ) {
+  assert( E > 0.0 && "Radiation :: ComputeClosure :: Non positive definite "
+                     "radiation energy density." );
   if ( E == 0.0 ) return 0.0; // This is a hack
   const Real f = FluxFactor( E, F );
   const Real chi =
diff --git a/src/utils/error.hpp b/src/utils/error.hpp
index a24db94..24419c5 100644
--- a/src/utils/error.hpp
+++ b/src/utils/error.hpp
@@ -3,19 +3,67 @@
  * --------------
  *
  * Author  : Brandon L. Barker
- * Purpose : Print error messages ...
+ * Purpose : Error throwing class, state checking
  **/
 
-#ifndef _ERROR_HPP_
-#define _ERROR_HPP_
+#ifndef ERROR_HPP_
+#define ERROR_HPP_
 
 #include <assert.h> /* assert */
+#include <exception>
 #include <stdexcept>
 
+#include "constants.hpp"
+#include "state.hpp"
+
 class Error : public std::runtime_error {
 
  public:
   Error( const std::string &message ) : std::runtime_error( message ) {}
 };
 
-#endif // _ERROR_HPP_
+template <typename T>
+void check_state( T state, const int ihi ) {
+
+  auto uCR = state->Get_uCR( );
+  auto uCF = state->Get_uCF( );
+
+  const Real c = constants::c_cgs;
+
+  Kokkos::parallel_for(
+      "check_state", ihi + 2, KOKKOS_LAMBDA( const int iX ) {
+        const Real tau   = uCF( 0, iX, 0 ); // cell averages checked
+        const Real vel   = uCF( 1, iX, 0 );
+        const Real e_m   = uCF( 2, iX, 0 );
+        const Real e_rad = uCR( 0, iX, 0 );
+        const Real f_rad = uCR( 1, iX, 0 );
+
+        assert( tau > 0.0 && " ! Check state :: Negative or zero density!" );
+        assert( !std::isnan( tau ) &&
+                " ! Check state :: Specific volume NaN!" );
+
+        assert(
+            std::fabs( vel ) < c &&
+            " ! Check state :: Velocity reached or exceeded speed of light!" );
+        assert( !std::isnan( vel ) && " ! Check state :: Velocity NaN!" );
+
+        assert( e_m > 0.0 &&
+                " Check state :: ! Negative or zero specific total energy!" );
+        assert( !std::isnan( e_m ) &&
+                " ! Check state :: Specific energy NaN!" );
+
+        assert(
+            e_rad > 0.0 &&
+            " ! Check state :: Negative or zero radiation energy density!" );
+        assert( !std::isnan( e_rad ) &&
+                " ! Check state :: Radiation energy NaN!" );
+
+        // TODO: radiation flux bound
+        // assert ( f_rad >= 0.0 && " Check state :: ! Negative or zero
+        // radiation energy density!" );
+        assert( !std::isnan( f_rad ) &&
+                " ! Check state :: Radiation flux NaN!" );
+      } );
+}
+
+#endif // ERROR_HPP_
diff --git a/src/utils/utilities.hpp b/src/utils/utilities.hpp
index a0e9682..6e50c4b 100644
--- a/src/utils/utilities.hpp
+++ b/src/utils/utilities.hpp
@@ -1,5 +1,5 @@
-#ifndef _UTILITIES_HPP_
-#define _UTILITIES_HPP_
+#ifndef UTILITIES_HPP_
+#define UTILITIES_HPP_
 
 #include <algorithm>
 #include <cctype>
@@ -46,4 +46,4 @@ T to_lower( T data ) {
 
 } // namespace utilities
 
-#endif // _UTILITIES_HPP_
+#endif // UTILITIES_HPP_