diff --git a/include/cantera/numerics/IDA_Solver.h b/include/cantera/numerics/IDA_Solver.h index 28f8fdad7..0c7fcbc13 100644 --- a/include/cantera/numerics/IDA_Solver.h +++ b/include/cantera/numerics/IDA_Solver.h @@ -240,13 +240,13 @@ protected: doublereal m_t0; //! Current value of the solution vector - void* m_y; + N_Vector m_y; //! Current value of the derivative of the solution vector - void* m_ydot; - void* m_id; - void* m_constraints; - void* m_abstol; + N_Vector m_ydot; + N_Vector m_id; + N_Vector m_constraints; + N_Vector m_abstol; int m_type; int m_itol; diff --git a/src/numerics/CVodeInt.cpp b/src/numerics/CVodeInt.cpp index 0b72565d0..9ca07140f 100644 --- a/src/numerics/CVodeInt.cpp +++ b/src/numerics/CVodeInt.cpp @@ -15,16 +15,10 @@ using namespace std; #include "../../ext/cvode/include/cvdense.h" #include "../../ext/cvode/include/cvdiag.h" #include "../../ext/cvode/include/cvspgmr.h" -#include "../../ext/cvode/include/nvector.h" #include "../../ext/cvode/include/cvode.h" #include "cantera/base/stringUtils.h" -inline static N_Vector nv(void* x) -{ - return reinterpret_cast(x); -} - extern "C" { /** @@ -110,21 +104,21 @@ CVodeInt::~CVodeInt() CVodeFree(m_cvode_mem); } if (m_y) { - N_VFree(nv(m_y)); + N_VFree(m_y); } if (m_abstol) { - N_VFree(nv(m_abstol)); + N_VFree(m_abstol); } delete[] m_iopt; } double& CVodeInt::solution(size_t k) { - return N_VIth(nv(m_y), int(k)); + return N_VIth(m_y, int(k)); } double* CVodeInt::solution() { - return N_VDATA(nv(m_y)); + return N_VDATA(m_y); } void CVodeInt::setTolerances(double reltol, size_t n, double* abstol) @@ -133,12 +127,12 @@ void CVodeInt::setTolerances(double reltol, size_t n, double* abstol) m_nabs = int(n); if (m_nabs != m_neq) { if (m_abstol) { - N_VFree(nv(m_abstol)); + N_VFree(m_abstol); } - m_abstol = reinterpret_cast(N_VNew(m_nabs, 0)); + m_abstol = N_VNew(m_nabs, 0); } for (int i=0; i(N_VNew(m_neq, 0)); // allocate solution vector + m_y = N_VNew(m_neq, 0); // allocate solution vector // check abs tolerance array size if (m_itol == 1 && m_nabs < m_neq) { throw CVodeErr("not enough absolute tolerance values specified."); } - func.getInitialConditions(m_t0, m_neq, N_VDATA(nv(m_y))); + func.getInitialConditions(m_t0, m_neq, N_VDATA(m_y)); // set options m_iopt[MXSTEP] = m_maxsteps; @@ -223,12 +217,12 @@ void CVodeInt::initialize(double t0, FuncEval& func) m_data = (void*)&func; if (m_itol) { - m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, nv(m_y), m_method, + m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, m_y, m_method, m_iter, m_itol, &m_reltol, - nv(m_abstol), m_data, NULL, 1, m_iopt, + m_abstol, m_data, NULL, 1, m_iopt, DATA_PTR(m_ropt), NULL); } else { - m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, nv(m_y), m_method, + m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, m_y, m_method, m_iter, m_itol, &m_reltol, &m_abstols, m_data, NULL, 1, m_iopt, DATA_PTR(m_ropt), NULL); @@ -255,7 +249,7 @@ void CVodeInt::initialize(double t0, FuncEval& func) void CVodeInt::reinitialize(double t0, FuncEval& func) { m_t0 = t0; - func.getInitialConditions(m_t0, m_neq, N_VDATA(nv(m_y))); + func.getInitialConditions(m_t0, m_neq, N_VDATA(m_y)); // set options m_iopt[MXSTEP] = m_maxsteps; @@ -268,12 +262,12 @@ void CVodeInt::reinitialize(double t0, FuncEval& func) m_data = (void*)&func; int result; if (m_itol) { - result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, nv(m_y), m_method, + result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, m_y, m_method, m_iter, m_itol, &m_reltol, - nv(m_abstol), m_data, NULL, 1, m_iopt, + m_abstol, m_data, NULL, 1, m_iopt, DATA_PTR(m_ropt), NULL); } else { - result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, nv(m_y), m_method, + result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, m_y, m_method, m_iter, m_itol, &m_reltol, &m_abstols, m_data, NULL, 1, m_iopt, DATA_PTR(m_ropt), NULL); @@ -301,7 +295,7 @@ void CVodeInt::integrate(double tout) { double t; int flag; - flag = CVode(m_cvode_mem, tout, nv(m_y), &t, NORMAL); + flag = CVode(m_cvode_mem, tout, m_y, &t, NORMAL); if (flag != SUCCESS) { throw CVodeErr(" CVode error encountered. Error code: " + int2str(flag)); } @@ -311,7 +305,7 @@ double CVodeInt::step(double tout) { double t; int flag; - flag = CVode(m_cvode_mem, tout, nv(m_y), &t, ONE_STEP); + flag = CVode(m_cvode_mem, tout, m_y, &t, ONE_STEP); if (flag != SUCCESS) { throw CVodeErr(" CVode error encountered. Error code: " + int2str(flag)); } diff --git a/src/numerics/CVodeInt.h b/src/numerics/CVodeInt.h index c6f234c60..b06918b01 100644 --- a/src/numerics/CVodeInt.h +++ b/src/numerics/CVodeInt.h @@ -10,6 +10,7 @@ #include "cantera/numerics/FuncEval.h" #include "cantera/base/ctexceptions.h" #include "cantera/base/ct_defs.h" +#include "../../ext/cvode/include/nvector.h" namespace Cantera { @@ -66,7 +67,7 @@ private: int m_neq; void* m_cvode_mem; double m_t0; - void* m_y, *m_abstol; + N_Vector m_y, m_abstol; int m_type; int m_itol; int m_method; diff --git a/src/numerics/IDA_Solver.cpp b/src/numerics/IDA_Solver.cpp index 2ea1408eb..c4cb7b198 100644 --- a/src/numerics/IDA_Solver.cpp +++ b/src/numerics/IDA_Solver.cpp @@ -26,12 +26,6 @@ typedef int sd_size_t; typedef long int sd_size_t; #endif - -inline static N_Vector nv(void* x) -{ - return reinterpret_cast(x); -} - namespace Cantera { @@ -168,52 +162,52 @@ IDA_Solver::~IDA_Solver() IDAFree(&m_ida_mem); } if (m_y) { - N_VDestroy_Serial(nv(m_y)); + N_VDestroy_Serial(m_y); } if (m_ydot) { - N_VDestroy_Serial(nv(m_ydot)); + N_VDestroy_Serial(m_ydot); } if (m_abstol) { - N_VDestroy_Serial(nv(m_abstol)); + N_VDestroy_Serial(m_abstol); } if (m_constraints) { - N_VDestroy_Serial(nv(m_constraints)); + N_VDestroy_Serial(m_constraints); } delete m_fdata; } doublereal IDA_Solver::solution(int k) const { - return NV_Ith_S(nv(m_y),k); + return NV_Ith_S(m_y,k); } const doublereal* IDA_Solver::solutionVector() const { - return NV_DATA_S(nv(m_y)); + return NV_DATA_S(m_y); } doublereal IDA_Solver::derivative(int k) const { - return NV_Ith_S(nv(m_ydot),k); + return NV_Ith_S(m_ydot,k); } const doublereal* IDA_Solver::derivativeVector() const { - return NV_DATA_S(nv(m_ydot)); + return NV_DATA_S(m_ydot); } void IDA_Solver::setTolerances(double reltol, double* abstol) { m_itol = IDA_SV; if (!m_abstol) { - m_abstol = reinterpret_cast(N_VNew_Serial(m_neq)); + m_abstol = N_VNew_Serial(m_neq); } for (int i = 0; i < m_neq; i++) { - NV_Ith_S(nv(m_abstol), i) = abstol[i]; + NV_Ith_S(m_abstol, i) = abstol[i]; } m_reltol = reltol; if (m_ida_mem) { - int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol)); + int flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol); if (flag != IDA_SUCCESS) { throw IDA_Err("Memory allocation failed."); } @@ -321,30 +315,30 @@ void IDA_Solver::init(doublereal t0) m_told_old = t0; m_tcurrent = t0; if (m_y) { - N_VDestroy_Serial(nv(m_y)); + N_VDestroy_Serial(m_y); } if (m_ydot) { - N_VDestroy_Serial(nv(m_ydot)); + N_VDestroy_Serial(m_ydot); } if (m_id) { - N_VDestroy_Serial(nv(m_id)); + N_VDestroy_Serial(m_id); } if (m_constraints) { - N_VDestroy_Serial(nv(m_constraints)); + N_VDestroy_Serial(m_constraints); } - m_y = reinterpret_cast(N_VNew_Serial(m_neq)); - m_ydot = reinterpret_cast(N_VNew_Serial(m_neq)); - m_constraints = reinterpret_cast(N_VNew_Serial(m_neq)); + m_y = N_VNew_Serial(m_neq); + m_ydot = N_VNew_Serial(m_neq); + m_constraints = N_VNew_Serial(m_neq); for (int i=0; i= 24 - flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot)); + flag = IDAInit(m_ida_mem, ida_resid, m_t0, m_y, m_ydot); if (flag != IDA_SUCCESS) { if (flag == IDA_MEM_FAIL) { throw IDA_Err("Memory allocation failed."); @@ -383,7 +377,7 @@ void IDA_Solver::init(doublereal t0) throw IDA_Err("IDAMalloc failed."); } } - flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol)); + flag = IDASVtolerances(m_ida_mem, m_reltol, m_abstol); if (flag != IDA_SUCCESS) { throw IDA_Err("Memory allocation failed."); } @@ -391,7 +385,7 @@ void IDA_Solver::init(doublereal t0) } else { #if SUNDIALS_VERSION <= 23 // scalar atol - flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot), + flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, m_y, m_ydot, m_itol, m_reltol, &m_abstols); if (flag != IDA_SUCCESS) { if (flag == IDA_MEM_FAIL) { @@ -404,7 +398,7 @@ void IDA_Solver::init(doublereal t0) } #elif SUNDIALS_VERSION >= 24 - flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot)); + flag = IDAInit(m_ida_mem, ida_resid, m_t0, m_y, m_ydot); if (flag != IDA_SUCCESS) { if (flag == IDA_MEM_FAIL) { throw IDA_Err("Memory allocation failed."); @@ -530,12 +524,12 @@ void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doubl } - flag = IDAGetConsistentIC(m_ida_mem, nv(m_y), nv(m_ydot)); + flag = IDAGetConsistentIC(m_ida_mem, m_y, m_ydot); if (flag != IDA_SUCCESS) { throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag)); } - doublereal* yy = NV_DATA_S(nv(m_y)); - doublereal* yyp = NV_DATA_S(nv(m_ydot)); + doublereal* yy = NV_DATA_S(m_y); + doublereal* yyp = NV_DATA_S(m_ydot); for (int i = 0; i < m_neq; i++) { y[i] = yy[i]; @@ -561,12 +555,12 @@ void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, dou } - flag = IDAGetConsistentIC(m_ida_mem, nv(m_y), nv(m_ydot)); + flag = IDAGetConsistentIC(m_ida_mem, m_y, m_ydot); if (flag != IDA_SUCCESS) { throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag)); } - doublereal* yy = NV_DATA_S(nv(m_y)); - doublereal* yyp = NV_DATA_S(nv(m_ydot)); + doublereal* yy = NV_DATA_S(m_y); + doublereal* yyp = NV_DATA_S(m_ydot); for (int i = 0; i < m_neq; i++) { y[i] = yy[i]; @@ -588,7 +582,7 @@ int IDA_Solver::solve(double tout) } m_told_old = m_told; m_told = m_tcurrent; - flag = IDASolve(m_ida_mem, tout, &tretn, nv(m_y), nv(m_ydot), IDA_ONE_STEP); + flag = IDASolve(m_ida_mem, tout, &tretn, m_y, m_ydot, IDA_ONE_STEP); if (flag < 0) { throw IDA_Err(" IDA error encountered."); } else if (flag == IDA_TSTOP_RETURN) { @@ -617,7 +611,7 @@ double IDA_Solver::step(double tout) } m_told_old = m_told; m_told = m_tcurrent; - flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP); + flag = IDASolve(m_ida_mem, tout, &t, m_y, m_ydot, IDA_ONE_STEP); if (flag < 0) { throw IDA_Err(" IDA error encountered."); } else if (flag == IDA_TSTOP_RETURN) {