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

Remove try/catch from ReactorNet::evalJacobian

Browse source 
This function isn't actually called by CVODE (which uses its own difference
quotient method for calculating the Jacobian), so the note here isn't correct.
This commit is contained in:
Ray Speth 2014-03-28 23:12:17 +00:00
parent b2c050e545
commit 4587883d5e
1 changed files with 14 additions and 21 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

35
src/zeroD/ReactorNet.cpp
View file

@ -207,31 +207,24 @@ void ReactorNet::evalJacobian(doublereal t, doublereal* y,
doublereal ysave, dy;
Array2D& jac = *j;
// use a try... catch block, since exceptions are not passed
// through CVODE, since it is C code
try {
//evaluate the unperturbed ydot
eval(t, y, ydot, p);
for (size_t n = 0; n < m_nv; n++) {
//evaluate the unperturbed ydot
eval(t, y, ydot, p);
for (size_t n = 0; n < m_nv; n++) {
// perturb x(n)
ysave = y[n];
dy = m_atol[n] + fabs(ysave)*m_rtol;
y[n] = ysave + dy;
dy = y[n] - ysave;
// perturb x(n)
ysave = y[n];
dy = m_atol[n] + fabs(ysave)*m_rtol;
y[n] = ysave + dy;
dy = y[n] - ysave;
// calculate perturbed residual
eval(t, y, DATA_PTR(m_ydot), p);
// calculate perturbed residual
eval(t, y, DATA_PTR(m_ydot), p);
// compute nth column of Jacobian
for (size_t m = 0; m < m_nv; m++) {
jac(m,n) = (m_ydot[m] - ydot[m])/dy;
}
y[n] = ysave;
// compute nth column of Jacobian
for (size_t m = 0; m < m_nv; m++) {
jac(m,n) = (m_ydot[m] - ydot[m])/dy;
}
} catch (CanteraError& err) {
std::cerr << err.what() << std::endl;
error("Terminating execution.");
y[n] = ysave;
}
}