diff --git a/include/amici/defines.h b/include/amici/defines.h index 747553da9b..411e58e7c0 100644 --- a/include/amici/defines.h +++ b/include/amici/defines.h @@ -73,6 +73,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; diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index 08d8bb2a5c..2e2e7f2bcc 100644 --- a/include/amici/model_dae.h +++ b/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( std::reduce(idlist.begin(), idlist.end()) diff --git a/include/amici/vector.h b/include/amici/vector.h index 23642a2429..e3336867e9 100644 --- a/include/amici/vector.h +++ b/include/amici/vector.h @@ -40,29 +40,31 @@ class AmiVector { */ AmiVector() = default; - /** Creates an std::vector and attaches the + /** + * @brief Construct empty vector of given size + * + * Creates an std::vector 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(length), 0.0) + : vec_(static_cast(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 rvec, SUNContext sunctx) - : sunctx_(sunctx) - , vec_(std::move(rvec)) + : vec_(std::move(rvec)) , nvec_(N_VMake_Serial( gsl::narrow(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,23 @@ 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 +276,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 +411,11 @@ class AmiVectorArray { */ void copy(AmiVectorArray const& other); + + private: /** SUNDIALS context */ SUNContext sunctx_{nullptr}; - private: /** main data storage */ std::vector vec_array_; diff --git a/src/amici.cpp b/src/amici.cpp index 20012d5dd3..702f9a4df5 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -36,6 +36,8 @@ std::map 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"}, diff --git a/src/solver.cpp b/src/solver.cpp index 368c760db0..961d285f8e 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -3,6 +3,7 @@ #include "amici/exception.h" #include "amici/model.h" #include "amici/symbolic_functions.h" +#include "amici/amici.h" #include @@ -12,78 +13,41 @@ 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::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); #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(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 +91,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_; } diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp index 3fc3fa22fb..5fa2fa42bb 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -65,9 +65,8 @@ 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. diff --git a/src/sundials_linsol_wrapper.cpp b/src/sundials_linsol_wrapper.cpp index 3ebc97d1ca..835ce3454d 100644 --- a/src/sundials_linsol_wrapper.cpp +++ b/src/sundials_linsol_wrapper.cpp @@ -160,9 +160,9 @@ SUNLinSolBand::SUNLinSolBand(N_Vector x, SUNMatrix A) } SUNLinSolBand::SUNLinSolBand(AmiVector const& x, int ubw, int lbw) - : A_(SUNMatrixWrapper(x.getLength(), ubw, lbw, x.sunctx_)) { + : A_(SUNMatrixWrapper(x.getLength(), ubw, lbw, x.get_ctx())) { solver_ - = SUNLinSol_Band(const_cast(x.getNVector()), A_, x.sunctx_); + = SUNLinSol_Band(const_cast(x.getNVector()), A_, x.get_ctx()); if (!solver_) throw AmiException("Failed to create solver."); } @@ -170,9 +170,9 @@ SUNLinSolBand::SUNLinSolBand(AmiVector const& x, int ubw, int lbw) SUNMatrix SUNLinSolBand::getMatrix() const { return A_.get(); } SUNLinSolDense::SUNLinSolDense(AmiVector const& x) - : A_(SUNMatrixWrapper(x.getLength(), x.getLength(), x.sunctx_)) { + : A_(SUNMatrixWrapper(x.getLength(), x.getLength(), x.get_ctx())) { solver_ - = SUNLinSol_Dense(const_cast(x.getNVector()), A_, x.sunctx_); + = SUNLinSol_Dense(const_cast(x.getNVector()), A_, x.get_ctx()); if (!solver_) throw AmiException("Failed to create solver."); } @@ -189,7 +189,7 @@ SUNLinSolKLU::SUNLinSolKLU( AmiVector const& x, int nnz, int sparsetype, StateOrdering ordering ) : A_(SUNMatrixWrapper( - x.getLength(), x.getLength(), nnz, sparsetype, x.sunctx_ + x.getLength(), x.getLength(), nnz, sparsetype, x.get_ctx() )) { solver_ = SUNLinSol_KLU( const_cast(x.getNVector()), A_, A_.get()->sunctx @@ -250,7 +250,7 @@ SUNLinSolSPBCGS::SUNLinSolSPBCGS(N_Vector x, int pretype, int maxl) SUNLinSolSPBCGS::SUNLinSolSPBCGS(AmiVector const& x, int pretype, int maxl) { solver_ = SUNLinSol_SPBCGS( - const_cast(x.getNVector()), pretype, maxl, x.sunctx_ + const_cast(x.getNVector()), pretype, maxl, x.get_ctx() ); if (!solver_) throw AmiException("Failed to create solver."); @@ -284,7 +284,7 @@ N_Vector SUNLinSolSPBCGS::getResid() const { SUNLinSolSPFGMR::SUNLinSolSPFGMR(AmiVector const& x, int pretype, int maxl) : SUNLinSolWrapper(SUNLinSol_SPFGMR( - const_cast(x.getNVector()), pretype, maxl, x.sunctx_ + const_cast(x.getNVector()), pretype, maxl, x.get_ctx() )) { if (!solver_) throw AmiException("Failed to create solver."); @@ -318,7 +318,7 @@ N_Vector SUNLinSolSPFGMR::getResid() const { SUNLinSolSPGMR::SUNLinSolSPGMR(AmiVector const& x, int pretype, int maxl) : SUNLinSolWrapper(SUNLinSol_SPGMR( - const_cast(x.getNVector()), pretype, maxl, x.sunctx_ + const_cast(x.getNVector()), pretype, maxl, x.get_ctx() )) { if (!solver_) throw AmiException("Failed to create solver."); @@ -358,7 +358,7 @@ SUNLinSolSPTFQMR::SUNLinSolSPTFQMR(N_Vector x, int pretype, int maxl) SUNLinSolSPTFQMR::SUNLinSolSPTFQMR(AmiVector const& x, int pretype, int maxl) { solver_ = SUNLinSol_SPTFQMR( - const_cast(x.getNVector()), pretype, maxl, x.sunctx_ + const_cast(x.getNVector()), pretype, maxl, x.get_ctx() ); if (!solver_) throw AmiException("Failed to create solver."); diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp index 1228f72801..a5847d982f 100644 --- a/src/sundials_matrix_wrapper.cpp +++ b/src/sundials_matrix_wrapper.cpp @@ -149,8 +149,7 @@ void SUNMatrixWrapper::reallocate(sunindextype NNZ) { if (int ret = SUNSparseMatrix_Reallocate(matrix_, NNZ) != SUN_SUCCESS) throw std::runtime_error( - "SUNSparseMatrix_Reallocate failed with " - "error code " + "SUNSparseMatrix_Reallocate failed with error code " + std::to_string(ret) + "." ); @@ -165,8 +164,7 @@ void SUNMatrixWrapper::realloc() { "CSR_MAT."); if (int ret = SUNSparseMatrix_Realloc(matrix_) != SUN_SUCCESS) throw std::runtime_error( - "SUNSparseMatrix_Realloc failed with " - "error code " + "SUNSparseMatrix_Realloc failed with error code " + std::to_string(ret) + "." ); @@ -510,7 +508,7 @@ void SUNMatrixWrapper::sparse_sum(std::vector const& mats) { for (acol = 0; acol < columns(); acol++) { set_indexptr(acol, nnz); /* column j of A starts here */ - for (auto& mat : mats) + for (auto const& mat : mats) nnz = mat.scatter( acol, 1.0, w.data(), gsl::make_span(x), acol + 1, this, nnz ); diff --git a/src/vector.cpp b/src/vector.cpp index a731fb3029..3e4e2852a7 100644 --- a/src/vector.cpp +++ b/src/vector.cpp @@ -6,9 +6,8 @@ namespace amici { AmiVector& AmiVector::operator=(AmiVector const& other) { - sunctx_ = other.sunctx_; vec_ = other.vec_; - synchroniseNVector(); + synchroniseNVector(other.get_ctx()); return *this; } @@ -56,13 +55,13 @@ void AmiVector::copy(AmiVector const& other) { std::copy(other.vec_.begin(), other.vec_.end(), vec_.begin()); } -void AmiVector::synchroniseNVector() { +void AmiVector::synchroniseNVector(SUNContext sunctx) { if (nvec_) N_VDestroy_Serial(nvec_); nvec_ = vec_.empty() ? nullptr : N_VMake_Serial( - gsl::narrow(vec_.size()), vec_.data(), sunctx_ + gsl::narrow(vec_.size()), vec_.data(), sunctx ); }