Add more optional error checking (#24)

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

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_

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!" );
   }
 }

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

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;

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,14 +44,24 @@ 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 );
 
         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;

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" );
   }
 

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_

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_

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,50 +67,71 @@ 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" )
                 ? PolyBasis::Legendre
                 : 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" );
 }