From 87062d115d02550016ea03b145f582bd47f7ef48 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sat, 11 Aug 2012 23:58:37 +0000
Subject: [PATCH] Integration failures with CVODES now generate a more
 informative exception

The exception message now lists the solution components and the corresponding
weighted error estimate which should indicate which variables are responsible
for the integration failure.
---
 src/numerics/CVodesIntegrator.cpp | 34 +++++++++++++++++++++++++++++--
 src/numerics/CVodesIntegrator.h   |  6 ++++++
 2 files changed, 38 insertions(+), 2 deletions(-)

diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp
index 05b533e01..224dd6aa2 100644
--- a/src/numerics/CVodesIntegrator.cpp
+++ b/src/numerics/CVodesIntegrator.cpp
@@ -47,6 +47,9 @@ using namespace std;
 
 #endif
 
+#include <sstream>
+#include <algorithm>
+
 inline static N_Vector nv(void* x)
 {
     return reinterpret_cast<N_Vector>(x);
@@ -492,7 +495,8 @@ void CVodesIntegrator::integrate(double tout)
     int flag;
     flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL);
     if (flag != CV_SUCCESS) {
-        throw CVodesErr(" CVodes error encountered. Error code: " + int2str(flag));
+        throw CVodesErr(" CVodes error encountered. Error code: " + int2str(flag) +
+                        "\nComponents with largest weighted error estimates:\n" + getErrorInfo(10));
     }
 #if SUNDIALS_VERSION <= 23
     if (m_np > 0) {
@@ -515,7 +519,9 @@ double CVodesIntegrator::step(double tout)
     int flag;
     flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP);
     if (flag != CV_SUCCESS) {
-        throw CVodesErr(" CVodes error encountered. Error code: " + int2str(flag));
+        throw CVodesErr(" CVodes error encountered. Error code: " + int2str(flag) +
+                        "\nComponents with largest weighted error estimates:\n" + getErrorInfo(10));
+
     }
     return t;
 }
@@ -538,6 +544,30 @@ double CVodesIntegrator::sensitivity(size_t k, size_t p)
     }
     return NV_Ith_S(m_yS[p],k);
 }
+
+string CVodesIntegrator::getErrorInfo(int N) {
+    N_Vector errs = N_VNew_Serial(m_neq);
+    N_Vector errw = N_VNew_Serial(m_neq);
+    CVodeGetErrWeights(m_cvode_mem, errw);
+    CVodeGetEstLocalErrors(m_cvode_mem, errs);
+
+    vector<pair<pair<double, double>, size_t> > weightedErrors;
+    for (size_t i=0; i<m_neq; i++) {
+        double err = NV_Ith_S(errs, i) * NV_Ith_S(errw, i);
+        weightedErrors.push_back(make_pair(make_pair(-abs(err), err), i));
+    }
+    N_VDestroy(errs);
+    N_VDestroy(errw);
+
+    sort(weightedErrors.begin(), weightedErrors.end());
+    stringstream s;
+    for (int i=0; i<N; i++) {
+        s << weightedErrors[i].second << ": "
+          << weightedErrors[i].first.second << endl;
+    }
+    return s.str();
+}
+
 }
 
 
diff --git a/src/numerics/CVodesIntegrator.h b/src/numerics/CVodesIntegrator.h
index cb6f69029..dee382e49 100644
--- a/src/numerics/CVodesIntegrator.h
+++ b/src/numerics/CVodesIntegrator.h
@@ -79,6 +79,12 @@ public:
     }
     virtual double sensitivity(size_t k, size_t p);
 
+    //! Returns a string listing the weighted error estimates associated
+    //! with each solution component.
+    //! This information can be used to identify which variables are
+    //! responsible for integrator failures or unexpected small timesteps.
+    virtual std::string getErrorInfo(int N);
+
 private:
 
     void sensInit(double t0, FuncEval& func);