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

More integration of IDA into Cantera.

Browse source 
added cj parameter to some existing routines.
  Added correct delta_t to the parameter list of  ResidJacEval object
This commit is contained in:
Harry Moffat 2011-08-01 15:26:20 +00:00
parent 2f97180753
commit 78eb7b4127
6 changed files with 40 additions and 15 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/numerics/BEulerInt.cpp
View file

@ -656,7 +656,7 @@ namespace Cantera {
/********************************************************************
* Call the function to get a jacobian.
*/
m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f);
m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f);
#ifdef DEBUG_HKM
//double dddd = J(89, 89);
//checkFinite(dddd);

28
Cantera/src/numerics/IDA_Solver.cpp
View file

@ -66,6 +66,13 @@ extern "C" {
* the sundials data types N_Vector, etc.), we define this function as the single function that IDA always calls. The
* real evaluation of the residual is done by an instance of a subclass of ResidEval, passed in to this
* function as a pointer in the parameters.
*
* FROM IDA WRITEUP -> What the IDA solver expects as a return flag from its residual routines ------
* A IDAResFn res should return a value of 0 if successful, a positive
* value if a recoverable error occured (e.g. yy has an illegal value),
* or a negative value if a nonrecoverable error occured. In the latter
* case, the program halts. If a recoverable error occured, the integrator
* will attempt to correct and retry.
*/
static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void *f_data) {
double* ydata = NV_DATA_S(y);
@ -74,8 +81,15 @@ extern "C" {
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
Cantera::ResidJacEval* f = d->m_func;
Cantera::IDA_Solver *s = d->m_solver;
f->eval(t, ydata, ydotdata, rdata);
return 0;
double delta_t = s->getCurrentStepFromIDA();
// TODO evaluate evalType. Assumed to be Base_ResidEval
int retn = 0;
int flag = f->evalResidNJ(t, delta_t, ydata, ydotdata, rdata);
if (flag < 0) {
// This signals to IDA that a nonrecoverable error has occurred.
retn = flag;
}
return retn;
}
//! Function called by by IDA to evaluate the Jacobian, given y and ydot.
@ -102,7 +116,9 @@ extern "C" {
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
Cantera::ResidJacEval* f = d->m_func;
doublereal * const * colPts = Jac->cols;
doublereal delta_t = 0.0;
Cantera::IDA_Solver *s = d->m_solver;
double delta_t = s->getCurrentStepFromIDA();
// printf(" delta_t = %g 1/cj = %g\n", delta_t, 1.0/c_j);
f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata);
return 0;
}
@ -241,6 +257,12 @@ namespace Cantera {
m_tstop = tstop;
}
//====================================================================================================================
doublereal IDA_Solver::getCurrentStepFromIDA() {
doublereal hcur;
IDAGetCurrentStep(m_ida_mem, &hcur);
return hcur;
}
//====================================================================================================================
void IDA_Solver::setJacobianType(int formJac) {
m_formJac = formJac;
if (m_ida_mem) {

6
Cantera/src/numerics/IDA_Solver.h
View file

@ -109,6 +109,12 @@ namespace Cantera {
*/
virtual void setStopTime(doublereal tstop);
//! Get the current step size from IDA via a call
/*!
* @return Returns the current step size.
*/
virtual double getCurrentStepFromIDA();
//! Set the form of the jacobian
/*!

3
Cantera/src/numerics/NonlinearSolver.cpp
View file

@ -2979,6 +2979,7 @@ namespace Cantera {
int info;
doublereal ysave, ydotsave, dy;
int retn = 1;
/*
* Clear the factor flag
*/
@ -2987,7 +2988,7 @@ namespace Cantera {
/********************************************************************
* Call the function to get a jacobian.
*/
info = m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f);
info = m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f);
m_nJacEval++;
m_nfe++;
if (info != 1) {

6
Cantera/src/numerics/ResidJacEval.cpp
View file

@ -328,17 +328,13 @@ namespace Cantera {
* @param resid Value of the residual that is computed (output)
*/
int ResidJacEval::
evalJacobian(const doublereal t, const doublereal delta_t,
evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
const doublereal * const y,
const doublereal * const ydot,
SquareMatrix &J,
doublereal * const resid)
{
doublereal * const * jac_colPts = J.colPts();
doublereal cj = 0.0;
if (delta_t > 0.0) {
cj = 1.0/delta_t;
}
return evalJacobianDP(t, delta_t, cj, y, ydot, jac_colPts, resid);
}
//====================================================================================================================

10
Cantera/src/numerics/ResidJacEval.h
View file

@ -318,6 +318,7 @@ namespace Cantera {
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param cj Coefficient of yprime used in the evalulation of the jacobian
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param J Reference to the SquareMatrix object to be calculated (output)
@ -327,11 +328,9 @@ namespace Cantera {
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual int evalJacobian(const doublereal t, const doublereal delta_t,
const doublereal* const y,
const doublereal* const ydot,
SquareMatrix &J,
doublereal * const resid);
virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
const doublereal* const y, const doublereal* const ydot,
SquareMatrix &J, doublereal * const resid);
//! Calculate an analytical jacobian and the residual at the current time and values.
/*!
@ -339,6 +338,7 @@ namespace Cantera {
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param cj Coefficient of yprime used in the evalulation of the jacobian
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param J Reference to the SquareMatrix object to be calculated (output)