research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Eliminate unnecessary use of reinterpret_cast

Browse source
This commit is contained in:
Ray Speth 2013-05-31 15:55:48 +00:00
parent e43a4a8b3c
commit e059b5c60c
4 changed files with 62 additions and 73 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
include/cantera/numerics/IDA_Solver.h
View file

@ -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;

44
src/numerics/CVodeInt.cpp
View file

@ -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<N_Vector>(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<void*>(N_VNew(m_nabs, 0));
m_abstol = N_VNew(m_nabs, 0);
}
for (int i=0; i<m_nabs; i++) {
N_VIth(nv(m_abstol), i) = abstol[i];
N_VIth(m_abstol, i) = abstol[i];
}
m_reltol = reltol;
}
@ -201,14 +195,14 @@ void CVodeInt::initialize(double t0, FuncEval& func)
m_t0 = t0;
if (m_y) {
N_VFree(nv(m_y)); // free solution vector if already allocated
N_VFree(m_y); // free solution vector if already allocated
}
m_y = reinterpret_cast<void*>(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));
}

3
src/numerics/CVodeInt.h
View file

@ -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;

78
src/numerics/IDA_Solver.cpp
View file

@ -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<N_Vector>(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<void*>(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<void*>(N_VNew_Serial(m_neq));
m_ydot = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
m_constraints = reinterpret_cast<void*>(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<m_neq; i++) {
NV_Ith_S(nv(m_y), i) = 0.0;
NV_Ith_S(nv(m_ydot), i) = 0.0;
NV_Ith_S(nv(m_constraints), i) = 0.0;
NV_Ith_S(m_y, i) = 0.0;
NV_Ith_S(m_ydot, i) = 0.0;
NV_Ith_S(m_constraints, i) = 0.0;
}
// get the initial conditions
m_resid.getInitialConditions(m_t0, NV_DATA_S(nv(m_y)), NV_DATA_S(nv(m_ydot)));
m_resid.getInitialConditions(m_t0, NV_DATA_S(m_y), NV_DATA_S(m_ydot));
if (m_ida_mem) {
IDAFree(&m_ida_mem);
@ -360,8 +354,8 @@ void IDA_Solver::init(doublereal t0)
if (m_itol == IDA_SV) {
#if SUNDIALS_VERSION <= 23
// vector atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, nv(m_abstol));
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, m_y, m_ydot,
m_itol, m_reltol, m_abstol);
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
throw IDA_Err("Memory allocation failed.");
@ -373,7 +367,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.");
@ -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) {