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

Fixed sensitivity calculation when using CVodesIntegrator::step

Browse source 
Previously, the sensitivity coefficient matrix m_yS was being updated only when
the `integrate` method was used to advance the system. Now, the sensitivities
are updated when a coefficient is requested after taking a timestep, regardless
of the method used to advance the system.
This commit is contained in:
Ray Speth 2012-10-24 15:41:59 +00:00
parent a6dc994439
commit 7b76e762d2
2 changed files with 27 additions and 21 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

43
src/numerics/CVodesIntegrator.cpp
View file

@ -135,7 +135,8 @@ CVodesIntegrator::CVodesIntegrator() :
m_maxsteps(20000),
m_fdata(0),
m_np(0),
m_mupper(0), m_mlower(0)
m_mupper(0), m_mlower(0),
m_sens_ok(false)
{
//m_ropt.resize(OPT_SIZE,0.0);
//m_iopt = new long[OPT_SIZE];
@ -257,6 +258,7 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func)
{
m_np = func.nparams();
size_t nv = func.neq();
m_sens_ok = false;
doublereal* data;
N_Vector y;
@ -297,6 +299,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
{
m_neq = func.neq();
m_t0 = t0;
m_time = t0;
if (m_y) {
N_VDestroy_Serial(nv(m_y)); // free solution vector if already allocated
@ -430,6 +433,7 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func)
void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
{
m_t0 = t0;
m_time = t0;
//try {
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
//}
@ -491,39 +495,24 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
void CVodesIntegrator::integrate(double tout)
{
double t;
int flag;
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL);
int flag = CVode(m_cvode_mem, tout, nv(m_y), &m_time, CV_NORMAL);
if (flag != CV_SUCCESS) {
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) {
CVodeGetSens(m_cvode_mem, tout, m_yS);
}
#elif SUNDIALS_VERSION >= 24
double tretn;
if (m_np > 0) {
CVodeGetSens(m_cvode_mem, &tretn, m_yS);
if (fabs(tretn - tout) > 1.0E-5) {
throw CVodesErr("Time of Sensitivities different than time of tout");
}
}
#endif
m_sens_ok = false;
}
double CVodesIntegrator::step(double tout)
{
double t;
int flag;
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP);
int flag = CVode(m_cvode_mem, tout, nv(m_y), &m_time, CV_ONE_STEP);
if (flag != CV_SUCCESS) {
throw CVodesErr(" CVodes error encountered. Error code: " + int2str(flag) +
"\nComponents with largest weighted error estimates:\n" + getErrorInfo(10));
}
return t;
m_sens_ok = false;
return m_time;
}
int CVodesIntegrator::nEvals() const
@ -536,6 +525,18 @@ int CVodesIntegrator::nEvals() const
double CVodesIntegrator::sensitivity(size_t k, size_t p)
{
if (!m_sens_ok && m_np) {
#if SUNDIALS_VERSION <= 23
int flag = CVodeGetSens(m_cvode_mem, m_time, m_yS);
#elif SUNDIALS_VERSION >= 24
int flag = CVodeGetSens(m_cvode_mem, &m_time, m_yS);
#endif
if (flag != CV_SUCCESS) {
throw CVodesErr("CVodeGetSens failed. Error code: " + int2str(flag));
}
m_sens_ok = true;
}
if (k >= m_neq) {
throw CVodesErr("sensitivity: k out of range ("+int2str(p)+")");
}

5
src/numerics/CVodesIntegrator.h
View file

@ -92,6 +92,7 @@ private:
size_t m_neq;
void* m_cvode_mem;
double m_t0;
double m_time; //!< The current integrator time
void* m_y, *m_abstol;
int m_type;
int m_itol;
@ -108,6 +109,10 @@ private:
N_Vector* m_yS;
size_t m_np;
int m_mupper, m_mlower;
//! Indicates whether the sensitivities stored in m_yS have been updated
//! for at the current integrator time.
bool m_sens_ok;
};
} // namespace