..

ThirdParty/sundials/src/sunmatrix/sparse/sunmatrix_sparse.c

@@ -621,6 +621,8 @@ SUNErrCode SUNMatCopy_Sparse(SUNMatrix A, SUNMatrix B)
   return SUN_SUCCESS;
 }
 
+// AMICI: this is the 5.8.0 version of SUNMatScaleAddI_Sparse
+//  see also https://github.com/LLNL/sundials/issues/581
 SUNErrCode SUNMatScaleAddI_Sparse(sunrealtype c, SUNMatrix A)
 {
     sunindextype j, i, p, nz, newvals, M, N, cend;

include/amici/defines.h

@@ -56,9 +56,6 @@ constexpr double pi = M_PI;
 
 // clang-format off
 
-constexpr int AMICI_ONEOUTPUT=                 5;
-
-// TODO: still needed? check whether all are still in sync
 // Return codes
 //
 // NOTE: When adding / removing / renaming return codes,
@@ -73,6 +70,8 @@ constexpr int AMICI_LSETUP_FAIL=              -6;
 constexpr int AMICI_RHSFUNC_FAIL=             -8;
 constexpr int AMICI_FIRST_RHSFUNC_ERR=        -9;
 constexpr int AMICI_CONSTR_FAIL=             -15;
+constexpr int AMICI_CVODES_CONSTR_FAIL=      -15;
+constexpr int AMICI_IDAS_CONSTR_FAIL=        -11;
 constexpr int AMICI_ILL_INPUT=               -22;
 constexpr int AMICI_ERROR=                   -99;
 constexpr int AMICI_NO_STEADY_STATE=         -81;

include/amici/model_dae.h

@@ -50,7 +50,7 @@ class Model_DAE : public Model {
               model_dimensions, simulation_parameters, o2mode, idlist, z2event,
               state_independent_events
           ) {
-        auto sunctx = derived_state_.x_pos_tmp_.sunctx_;
+        SUNContext sunctx = derived_state_.sunctx_;
         derived_state_.M_ = SUNMatrixWrapper(nx_solver, nx_solver, sunctx);
         auto M_nnz = static_cast<sunindextype>(
             std::reduce(idlist.begin(), idlist.end())

include/amici/model_state.h

@@ -208,8 +208,6 @@ struct ModelStateDerived {
      * We could pass the solver's context to the model during initialize() and
      * only create ModelStateDerived there, but this caused issue in tests with
      * FSA for unclear reasons.
-     * TODO: Do we need to attach an error handler here? Not sure if
-     * NVector and SUNMatrix would even use it.
      */
     sundials::Context sunctx_;
 

include/amici/newton_solver.h

@@ -36,7 +36,7 @@ class NewtonSolver {
      *
      * @param simulationSolver solver with settings
      * @param model pointer to the model instance
-     * @param sunctx SUNDIALS context
+     * @return NewtonSolver according to the specified linsolType
      */
     static std::unique_ptr<NewtonSolver>
     getSolver(Solver const& simulationSolver, Model const& model);

include/amici/serialization.h

@@ -19,8 +19,13 @@
 #include <boost/serialization/map.hpp>
 #include <boost/serialization/vector.hpp>
 
-/** @file serialization.h Helper functions and forward declarations for
- * boost::serialization */
+/**
+ * @file serialization.h Helper functions and forward declarations for
+ * boost::serialization
+ *
+ * FIXME: Many of the serialization functions seem to be outdated and miss
+ * some fields.
+ */
 namespace boost {
 namespace serialization {
 
@@ -284,6 +289,10 @@ void serialize(
 
 /**
  * @brief Serialize AmiVector to a boost archive
+ *
+ * NOTE: The deserialized AmiVector/NVector will not have a valid SUNContext.
+ *  This needs to be set afterwards.
+ *
  * @param ar archive
  * @param v AmiVector
  */
@@ -294,14 +303,9 @@ void serialize(
     if (Archive::is_loading::value) {
         std::vector<amici::realtype> tmp;
         ar & tmp;
-        // TODO: how do we get a new sunctx in here??
-        //  for now, create one for construction of the vector, set it to
-        //  nullptr
         sundials::Context sunctx;
         v = amici::AmiVector(tmp, sunctx);
-        if (v.getNVector()) {
-            v.getNVector()->sunctx = nullptr;
-        }
+        v.set_ctx(sunctx);
     } else {
         auto tmp = v.getVector();
         ar & tmp;

include/amici/steadystateproblem.h

@@ -378,15 +378,6 @@ class SteadystateProblem {
      */
     void getNewtonStep(Model& model);
 
-    /**
-     * @brief SUNDIALS context.
-     *
-     * We use a different context here, because we might
-     * use a different solver than the one used for the simulation, and we
-     * shouldn't use multiple solvers on the same context.
-     */
-    // TODO: attach error handler?
-    sundials::Context sunctx_;
     /** newton step */
     AmiVector delta_;
     /** previous newton step */

include/amici/vector.h

@@ -40,29 +40,31 @@ class AmiVector {
      */
     AmiVector() = default;
 
-    /** Creates an std::vector<realtype> and attaches the
+    /**
+     * @brief Construct empty vector of given size
+     *
+     * Creates an std::vector<realtype> and attaches the
      * data pointer to a newly created N_Vector_Serial.
      * Using N_VMake_Serial ensures that the N_Vector
      * module does not try to deallocate the data vector
      * when calling N_VDestroy_Serial
-     * @brief empty constructor
      * @param length number of elements in vector
      * @param sunctx SUNDIALS context
      */
     explicit AmiVector(long int const length, SUNContext sunctx)
-        : sunctx_(sunctx)
-        , vec_(static_cast<decltype(vec_)::size_type>(length), 0.0)
+        : vec_(static_cast<decltype(vec_)::size_type>(length), 0.0)
         , nvec_(N_VMake_Serial(length, vec_.data(), sunctx)) {}
 
-    /** Moves data from std::vector and constructs an nvec that points to the
+    /**
+     * @brief Constructor from std::vector
+     *
+     * Moves data from std::vector and constructs an nvec that points to the
      * data
-     * @brief constructor from std::vector,
      * @param rvec vector from which the data will be moved
      * @param sunctx SUNDIALS context
      */
     explicit AmiVector(std::vector<realtype> rvec, SUNContext sunctx)
-        : sunctx_(sunctx)
-        , vec_(std::move(rvec))
+        : vec_(std::move(rvec))
         , nvec_(N_VMake_Serial(
               gsl::narrow<long int>(vec_.size()), vec_.data(), sunctx
           )) {}
@@ -80,8 +82,7 @@ class AmiVector {
      * @param vold vector from which the data will be copied
      */
     AmiVector(AmiVector const& vold)
-        : sunctx_(vold.sunctx_)
-        , vec_(vold.vec_) {
+        : vec_(vold.vec_) {
         if (vold.nvec_ == nullptr) {
             nvec_ = nullptr;
             return;
@@ -97,10 +98,9 @@ class AmiVector {
      * @param other vector from which the data will be moved
      */
     AmiVector(AmiVector&& other) noexcept
-        : sunctx_(other.sunctx_)
+        : vec_(std::move(other.vec_))
         , nvec_(nullptr) {
-        vec_ = std::move(other.vec_);
-        synchroniseNVector();
+        synchroniseNVector(other.get_ctx());
     }
 
     /**
@@ -249,8 +249,26 @@ class AmiVector {
         Archive& ar, AmiVector& s, unsigned int version
     );
 
-    /** SUNDIALS context */
-    SUNContext sunctx_{nullptr};
+    /**
+     * @brief Get SUNContext
+     * @return The current SUNContext or nullptr, if this AmiVector is empty
+     */
+    SUNContext get_ctx() const {
+        return nvec_ == nullptr ? nullptr : nvec_->sunctx;
+    }
+
+    /**
+     * @brief Set SUNContext
+     *
+     * If this AmiVector is non-empty, changes the current SUNContext of the
+     * associated N_Vector. If empty, do nothing.
+     *
+     * @param ctx SUNDIALS context to set
+     */
+    void set_ctx(SUNContext ctx) {
+        if (nvec_)
+            nvec_->sunctx = ctx;
+    }
 
   private:
     /** main data storage */
@@ -261,8 +279,9 @@ class AmiVector {
 
     /**
      * @brief reconstructs nvec such that data pointer points to vec data array
+     * @param sunctx SUNDIALS context
      */
-    void synchroniseNVector();
+    void synchroniseNVector(SUNContext sunctx);
 };
 
 /**
@@ -395,10 +414,10 @@ class AmiVectorArray {
      */
     void copy(AmiVectorArray const& other);
 
+  private:
     /** SUNDIALS context */
     SUNContext sunctx_{nullptr};
 
-  private:
     /** main data storage */
     std::vector<AmiVector> vec_array_;
 

src/amici.cpp

@@ -36,6 +36,8 @@ std::map<int, std::string> simulation_status_to_str_map = {
     {AMICI_CONV_FAILURE, "AMICI_CONV_FAILURE"},
     {AMICI_FIRST_RHSFUNC_ERR, "AMICI_FIRST_RHSFUNC_ERR"},
     {AMICI_CONSTR_FAIL, "AMICI_CONSTR_FAIL"},
+    {AMICI_CVODES_CONSTR_FAIL, "AMICI_CVODES_CONSTR_FAIL"},
+    {AMICI_IDAS_CONSTR_FAIL, "AMICI_IDAS_CONSTR_FAIL"},
     {AMICI_RHSFUNC_FAIL, "AMICI_RHSFUNC_FAIL"},
     {AMICI_ILL_INPUT, "AMICI_ILL_INPUT"},
     {AMICI_ERROR, "AMICI_ERROR"},

src/newton_solver.cpp

@@ -26,7 +26,9 @@ NewtonSolver::getSolver(Solver const& simulationSolver, Model const& model) {
 
     /* DIRECT SOLVERS */
     case LinearSolver::dense:
-        solver.reset(new NewtonSolverDense(model, simulationSolver.getSunContext()));
+        solver.reset(
+            new NewtonSolverDense(model, simulationSolver.getSunContext())
+        );
         break;
 
     case LinearSolver::band:
@@ -55,7 +57,9 @@ NewtonSolver::getSolver(Solver const& simulationSolver, Model const& model) {
     case LinearSolver::SuperLUMT:
         throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");
     case LinearSolver::KLU:
-        solver.reset(new NewtonSolverSparse(model, simulationSolver.getSunContext()));
+        solver.reset(
+            new NewtonSolverSparse(model, simulationSolver.getSunContext())
+        );
         break;
     default:
         throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");

src/solver.cpp

@@ -1,5 +1,6 @@
 #include "amici/solver.h"
 
+#include "amici/amici.h"
 #include "amici/exception.h"
 #include "amici/model.h"
 #include "amici/symbolic_functions.h"
@@ -12,78 +13,44 @@
 
 namespace amici {
 
+/* Error handler passed to SUNDIALS. */
 void wrapErrHandlerFn(
     [[maybe_unused]] int line, char const* func, char const* file,
     char const* msg, SUNErrCode err_code, void* err_user_data,
     [[maybe_unused]] SUNContext sunctx
 ) {
     constexpr int BUF_SIZE = 250;
-    char buffer[BUF_SIZE];
-    char buffid[BUF_SIZE];
+
+    char msg_buffer[BUF_SIZE];
+    char id_buffer[BUF_SIZE];
     static_assert(
         std::is_same<SUNErrCode, int>::value, "Must update format string"
     );
+    // for debug builds, include full file path and line numbers
 #ifdef NDEBUG
-    snprintf(buffer, BUF_SIZE, "%s:%d: %s (%d)", file, line, msg, err_code);
+    snprintf(msg_buffer, BUF_SIZE, "%s:%d: %s (%d)", file, line, msg, err_code);
 #else
-    snprintf(buffer, BUF_SIZE, "%s", msg);
+    snprintf(msg_buffer, BUF_SIZE, "%s", msg);
 #endif
     // we need a matlab-compatible message ID
     // i.e. colon separated and only  [A-Za-z0-9_]
+    // see https://mathworks.com/help/matlab/ref/mexception.html
     std::filesystem::path path(file);
     auto file_stem = path.stem().string();
 
-    switch (err_code) {
-    case 99:
-        snprintf(buffid, BUF_SIZE, "%s:%s:WARNING", file_stem.c_str(), func);
-        break;
-
-    case AMICI_TOO_MUCH_WORK:
-        snprintf(
-            buffid, BUF_SIZE, "%s:%s:TOO_MUCH_WORK", file_stem.c_str(), func
-        );
-        break;
-
-    case AMICI_TOO_MUCH_ACC:
-        snprintf(
-            buffid, BUF_SIZE, "%s:%s:TOO_MUCH_ACC", file_stem.c_str(), func
-        );
-        break;
-
-    case AMICI_ERR_FAILURE:
-        snprintf(
-            buffid, BUF_SIZE, "%s:%s:ERR_FAILURE", file_stem.c_str(), func
-        );
-        break;
-
-    case AMICI_CONV_FAILURE:
-        snprintf(
-            buffid, BUF_SIZE, "%s:%s:CONV_FAILURE", file_stem.c_str(), func
-        );
-        break;
-
-    case AMICI_RHSFUNC_FAIL:
-        snprintf(
-            buffid, BUF_SIZE, "%s:%s:RHSFUNC_FAIL", file_stem.c_str(), func
-        );
-        break;
-
-    case AMICI_FIRST_RHSFUNC_ERR:
-        snprintf(
-            buffid, BUF_SIZE, "%s:%s:FIRST_RHSFUNC_ERR", file_stem.c_str(), func
-        );
-        break;
-    default:
-        snprintf(buffid, BUF_SIZE, "%s:%s:OTHER", file_stem.c_str(), func);
-        break;
-    }
+    auto err_code_str = simulation_status_to_str(err_code);
+    snprintf(
+        id_buffer, BUF_SIZE, "%s:%s:%s", file_stem.c_str(), func,
+        err_code_str.c_str()
+    );
 
     if (!err_user_data) {
         throw std::runtime_error("eh_data unset");
     }
+
     auto solver = static_cast<Solver const*>(err_user_data);
     if (solver->logger)
-        solver->logger->log(LogSeverity::debug, buffid, buffer);
+        solver->logger->log(LogSeverity::debug, id_buffer, msg_buffer);
 }
 
 Solver::Solver(Solver const& other)
@@ -127,11 +94,8 @@ Solver::Solver(Solver const& other)
     , max_step_size_(other.max_step_size_)
     , maxstepsB_(other.maxstepsB_)
     , sensi_(other.sensi_) {
-    // AmiVector.setContext()... check for nullptr
-    if (constraints_.data()) {
-        constraints_.sunctx_ = sunctx_;
-        constraints_.getNVector()->sunctx = sunctx_;
-    }
+    // update to our own context
+    constraints_.set_ctx(sunctx_);
 }
 
 SUNContext Solver::getSunContext() const { return sunctx_; }

src/solver_idas.cpp

@@ -65,9 +65,11 @@ static_assert(
     amici::AMICI_LSETUP_FAIL == IDA_LSETUP_FAIL,
     "AMICI_LSETUP_FAIL != IDA_LSETUP_FAIL"
 );
-// FIXME: this does not match CVODE, we need separate return values
-// static_assert(amici::AMICI_CONSTR_FAIL == IDA_CONSTR_FAIL, "AMICI_CONSTR_FAIL
-// != IDA_CONSTR_FAIL");
+// This does not match the CVODE code, we need separate return values
+static_assert(
+    amici::AMICI_IDAS_CONSTR_FAIL == IDA_CONSTR_FAIL,
+    "AMICI_IDAS_CONSTR_FAIL != IDA_CONSTR_FAIL"
+);
 
 /*
  * The following static members are callback function to IDAS.