Catch NaNs in static solver (rice-solar-physics#94)

* catch NaNs in static solver

* fix bug with integer

* add test to check for NaNs with adaptive solver

* merge tests into one function

---------

Co-authored-by: Jeffrey Reep <[email protected]>

source/main.cpp

@@ -106,6 +106,10 @@ int main(int argc, char *argv[])
   {
     // Constant timestep integration
     num_steps = boost::numeric::odeint::integrate_const( controlled_stepper, loop->CalculateDerivs, state, loop->parameters.tau, loop->parameters.total_time, loop->parameters.tau, obs->Observe);
+    if(obs->CheckNan(state))
+    {
+        throw std::runtime_error("NaNs were detected in the output.  Check the input configuration.");
+    }
   }
 
   //Print results to file

source/observer.cpp

@@ -60,3 +60,18 @@ int Observer::CheckNan(state_type &state, double &time, double &tau, double old_
   // Pass otherwise
   return 0;
 }
+
+int Observer::CheckNan(state_type &state)
+{
+  // Check for NaNs in the state
+  for(int j=0; j<state.size(); j++)
+  {
+    if(std::isnan(state[j]))
+    {
+      return 1;
+    }
+  }
+
+  // Pass otherwise
+  return 0;
+}

source/observer.h

@@ -55,7 +55,10 @@ class Observer {
   // Boost integrator does not check for NaNs so this is done manually. If a
   // NaN is found anywhere in the state vector, the state and time are set 
   // back to the previous step and the timestep is reduced.
+  // The overloaded function, for use with the static time step solver, only checks for NaNs
+  // and does not reset the state or time.  
   int CheckNan(state_type &state, double &time, double &tau, double old_time, double old_tau);
+  int CheckNan(state_type &state);
 };
 // Pointer to the <Observer> class
 typedef Observer* OBSERVER;

tests/test_solver.py

@@ -82,4 +82,16 @@ def test_insufficient_heating(base_config, value):
     config['use_adaptive_solver'] = False
     config['heating']['background'] = value
     with pytest.raises(EbtelPlusPlusError):
-        run_ebtelplusplus(config)
+        run_ebtelplusplus(config)
+
+@pytest.mark.parametrize('use_adaptive_solver', [True, False])
+def test_NaNs_in_solver(base_config, use_adaptive_solver):
+    config = base_config.copy()
+    config['use_adaptive_solver'] = use_adaptive_solver
+    config['heating']['events'] = [
+                {'event': {'rise_start': 0.0, 'rise_end': 100.0, 'decay_start': 100.0,
+                           'decay_end': 200.0, 'magnitude': -10.0}}
+            ]
+    with pytest.raises(EbtelPlusPlusError):
+        run_ebtelplusplus(config)
+