| 502 | } |
| 503 | |
| 504 | void ReactorNet::evalJacobian(double t, double* y, double* ydot, double* p, Array2D* j) |
| 505 | { |
| 506 | //evaluate the unperturbed ydot |
| 507 | eval(t, y, ydot, p); |
| 508 | for (size_t n = 0; n < m_nv; n++) { |
| 509 | // perturb x(n) |
| 510 | double ysave = y[n]; |
| 511 | double dy = m_atol[n] + fabs(ysave)*m_rtol; |
| 512 | y[n] = ysave + dy; |
| 513 | dy = y[n] - ysave; |
| 514 | |
| 515 | // calculate perturbed residual |
| 516 | eval(t, y, m_ydot.data(), p); |
| 517 | |
| 518 | // compute nth column of Jacobian |
| 519 | for (size_t m = 0; m < m_nv; m++) { |
| 520 | j->value(m,n) = (m_ydot[m] - ydot[m])/dy; |
| 521 | } |
| 522 | y[n] = ysave; |
| 523 | } |
| 524 | } |
| 525 | |
| 526 | void ReactorNet::updateState(double* y) |
| 527 | { |
no test coverage detected