[Reactor] Suppress recoverable errors in RHS function by default

Accumulate errors and show them if the integrator actually fails. In "verbose"
mode, the errors are shown as they occur.

include/cantera/numerics/CVodesIntegrator.h

@@ -94,6 +94,7 @@ private:
 
     size_t m_neq;
     void* m_cvode_mem;
+    FuncEval* m_func;
     double m_t0;
     double m_time; //!< The current integrator time
     N_Vector m_y, m_abstol;

include/cantera/numerics/FuncEval.h

@@ -26,7 +26,7 @@ namespace Cantera
 class FuncEval
 {
 public:
-    FuncEval() {}
+    FuncEval();
     virtual ~FuncEval() {}
 
     /**
@@ -38,6 +38,17 @@ public:
      */
     virtual void eval(double t, double* y, double* ydot, double* p)=0;
 
+    //! Evaluate the right-hand side using return code to indicate status.
+    /*!
+     * Errors are indicated using the return value, rather than by throwing
+     *  exceptions. This method is used when calling from a C-based integrator
+     *  such as CVODES. Exceptions may either be stored or printed, based on the
+     *  setting of suppressErrors().
+     *  @returns 0 for a successful evaluation; 1 after a potentially-
+     *      recoverable error; -1 after an unrecoverable error.
+     */
+    int eval_nothrow(double t, double* y, double* ydot);
+
     /**
      * Fill the solution vector with the initial conditions
      * at initial time t0.
@@ -62,6 +73,24 @@ public:
         return m_sens_params.size();
     }
 
+    //! Enable or disable suppression of errors when calling eval()
+    void suppressErrors(bool suppress) {
+        m_suppress_errors = suppress;
+    }
+
+    //! Get current state of error suppression
+    bool suppressErrors() const {
+        return m_suppress_errors;
+    };
+
+    //! Return a string containing the text of any suppressed errors
+    std::string getErrors() const;
+
+    //! Clear any previously-stored suppressed errors
+    void clearErrors() {
+        m_errors.clear();
+    };
+
     //! Values for the problem parameters for which sensitivities are computed
     //! This is the array which is perturbed and passed back as the fourth
     //! argument to eval().
@@ -69,6 +98,13 @@ public:
 
     //! Scaling factors for each sensitivity parameter
     vector_fp m_paramScales;
+
+protected:
+    // If true, errors are accumulated in m_errors. Otherwise, they are printed
+    bool m_suppress_errors;
+
+    //! Errors occuring during function evaluations
+    std::vector<std::string> m_errors;
 };
 
 }

include/cantera/zeroD/ReactorNet.h

@@ -105,6 +105,7 @@ public:
     //! reactor network.
     void setVerbose(bool v = true) {
         m_verbose = v;
+        suppressErrors(!m_verbose);
     }
 
     //! Return a reference to the integrator.

src/numerics/CVodesIntegrator.cpp

@@ -46,21 +46,8 @@ extern "C" {
      */
     static int cvodes_rhs(realtype t, N_Vector y, N_Vector ydot, void* f_data)
     {
-        try {
+        FuncEval* f = (FuncEval*) f_data;
-            FuncEval* f = (FuncEval*) f_data;
+        return f->eval_nothrow(t, NV_DATA_S(y), NV_DATA_S(ydot));
-            f->eval(t, NV_DATA_S(y), NV_DATA_S(ydot), f->m_sens_params.data());
-        } catch (CanteraError& err) {
-            std::cerr << err.what() << std::endl;
-            return 1; // possibly recoverable error
-        } catch (std::exception& err) {
-            std::cerr << "cvodes_rhs: unhandled exception:" << std::endl;
-            std::cerr << err.what() << std::endl;
-            return -1; // unrecoverable error
-        } catch (...) {
-            std::cerr << "cvodes_rhs: unhandled exception of uknown type" << std::endl;
-            return -1; // unrecoverable error
-        }
-        return 0; // successful evaluation
     }
 
     //! Function called by CVodes when an error is encountered instead of
@@ -78,6 +65,7 @@ extern "C" {
 CVodesIntegrator::CVodesIntegrator() :
     m_neq(0),
     m_cvode_mem(0),
+    m_func(0),
     m_t0(0.0),
     m_y(0),
     m_abstol(0),
@@ -251,6 +239,8 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
     m_neq = func.neq();
     m_t0 = t0;
     m_time = t0;
+    m_func = &func;
+    func.clearErrors();
 
     if (m_y) {
         N_VDestroy_Serial(m_y); // free solution vector if already allocated
@@ -333,6 +323,8 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
     m_t0 = t0;
     m_time = t0;
     func.getState(NV_DATA_S(m_y));
+    m_func = &func;
+    func.clearErrors();
 
     int result = CVodeReInit(m_cvode_mem, m_t0, m_y);
     if (result != CV_SUCCESS) {
@@ -393,10 +385,15 @@ void CVodesIntegrator::integrate(double tout)
     }
     int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_NORMAL);
     if (flag != CV_SUCCESS) {
+        string f_errs = m_func->getErrors();
+        if (!f_errs.empty()) {
+            f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
+        }
         throw CanteraError("CVodesIntegrator::integrate",
             "CVodes error encountered. Error code: {}\n{}\n"
+            "{}"
             "Components with largest weighted error estimates:\n{}",
-            flag, m_error_message, getErrorInfo(10));
+            flag, m_error_message, f_errs, getErrorInfo(10));
     }
     m_sens_ok = false;
 }
@@ -405,10 +402,15 @@ double CVodesIntegrator::step(double tout)
 {
     int flag = CVode(m_cvode_mem, tout, m_y, &m_time, CV_ONE_STEP);
     if (flag != CV_SUCCESS) {
+        string f_errs = m_func->getErrors();
+        if (!f_errs.empty()) {
+            f_errs = "Exceptions caught during RHS evaluation:\n" + f_errs;
+        }
         throw CanteraError("CVodesIntegrator::step",
             "CVodes error encountered. Error code: {}\n{}\n"
+            "{}"
             "Components with largest weighted error estimates:\n{}",
-            flag, m_error_message, getErrorInfo(10));
+            flag, f_errs, m_error_message, getErrorInfo(10));
 
     }
     m_sens_ok = false;

src/numerics/FuncEval.cpp

@@ -0,0 +1,54 @@
+#include "cantera/numerics/FuncEval.h"
+#include <sstream>
+
+namespace Cantera
+{
+
+FuncEval::FuncEval()
+    : m_suppress_errors(false)
+{
+}
+
+int FuncEval::eval_nothrow(double t, double* y, double* ydot)
+{
+    try {
+        eval(t, y, ydot, m_sens_params.data());
+    } catch (CanteraError& err) {
+        if (suppressErrors()) {
+            m_errors.push_back(err.getMessage());
+        } else {
+            writelog(err.what());
+        }
+        return 1; // possibly recoverable error
+    } catch (std::exception& err) {
+        if (suppressErrors()) {
+            m_errors.push_back(err.what());
+        } else {
+            writelog("FuncEval::eval_nothrow: unhandled exception:\n");
+            writelog(err.what());
+            writelogendl();
+        }
+        return -1; // unrecoverable error
+    } catch (...) {
+        std::string msg = "FuncEval::eval_nothrow: unhandled exception"
+            " of unknown type\n";
+        if (suppressErrors()) {
+            m_errors.push_back(msg);
+        } else {
+            writelog(msg);
+        }
+        return -1; // unrecoverable error
+    }
+    return 0; // successful evaluation
+}
+
+std::string FuncEval::getErrors() const {
+    std::stringstream errs;
+    for (const auto& err : m_errors) {
+        errs << err;
+        errs << "\n";
+    }
+    return errs.str();
+}
+
+}

src/zeroD/ReactorNet.cpp

@@ -22,6 +22,7 @@ ReactorNet::ReactorNet() :
     m_verbose(false)
 {
     m_integ = newIntegrator("CVODE");
+    suppressErrors(true);
 
     // use backward differencing, with a full Jacobian computed
     // numerically, and use a Newton linear iterator