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

Fixed an error in the IDA_Solver class.

Browse source 
Added an analytical jacobian capability to IDA_Solver.
Added column pointers as a member to DenseMatrix
This commit is contained in:
Harry Moffat 2011-07-31 23:44:45 +00:00
parent 82e194a2b0
commit 2f97180753
6 changed files with 200 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/numerics/DenseMatrix.h
View file

@ -113,7 +113,7 @@ namespace Cantera {
* @param m New number of columns
* @param v Default fill value. defaults to zero.
*/
void resize(int n, int m, doublereal v);
void resize(int n, int m, doublereal v = 0.0);
//! Return a vector of const pointers to the columns
/*!

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

@ -44,14 +44,16 @@ namespace Cantera {
public:
ResidData(ResidJacEval* f, int npar = 0) {
ResidData(ResidJacEval* f, IDA_Solver *s, int npar = 0) {
m_func = f;
m_solver = s;
}
virtual ~ResidData() {
}
ResidJacEval* m_func;
IDA_Solver * m_solver;
};
}
@ -71,9 +73,41 @@ extern "C" {
double* rdata = NV_DATA_S(r);
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;
}
//! Function called by by IDA to evaluate the Jacobian, given y and ydot.
/*!
*
*
* typedef int (*IDADlsDenseJacFn)(int N, realtype t, realtype c_j,
* N_Vector y, N_Vector yp, N_Vector r,
* DlsMat Jac, void *user_data,
* N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
*
* A IDADlsDenseJacFn should return
* 0 if successful,
* a positive int if a recoverable error occurred, or
* a negative int if a nonrecoverable error occurred.
* In the case of a recoverable error return, the integrator will
* attempt to recover by reducing the stepsize (which changes cj).
*/
static int ida_jacobian(int nrows, realtype t, realtype c_j, N_Vector y, N_Vector ydot, N_Vector r,
DlsMat Jac, void *f_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) {
doublereal * ydata = NV_DATA_S(y);
doublereal * ydotdata = NV_DATA_S(ydot);
doublereal * rdata = NV_DATA_S(r);
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
Cantera::ResidJacEval* f = d->m_func;
doublereal * const * colPts = Jac->cols;
doublereal delta_t = 0.0;
f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata);
return 0;
}
}
namespace Cantera {
@ -103,7 +137,12 @@ namespace Cantera {
m_h0(0.0),
m_maxsteps(20000),
m_maxord(0),
m_formJac(0),
m_tstop(0.0),
m_told_old(0.0),
m_told(0.0),
m_tcurrent(0.0),
m_deltat(0.0),
m_maxErrTestFails(-1),
m_maxNonlinIters(0),
m_maxNonlinConvFails(-1),
@ -151,10 +190,12 @@ namespace Cantera {
for (int i = 0; i < m_neq; i++) {
NV_Ith_S(nv(m_abstol), i) = abstol[i];
}
m_reltol = reltol;
int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol));
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
m_reltol = reltol;
if (m_ida_mem) {
int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol));
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
}
}
}
//====================================================================================================================
@ -162,9 +203,11 @@ namespace Cantera {
m_itol = IDA_SS;
m_reltol = reltol;
m_abstols = abstol;
int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
if (m_ida_mem) {
int flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
}
}
}
//====================================================================================================================
@ -198,6 +241,18 @@ namespace Cantera {
m_tstop = tstop;
}
//====================================================================================================================
void IDA_Solver::setJacobianType(int formJac) {
m_formJac = formJac;
if (m_ida_mem) {
if (m_formJac == 1) {
int flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDADlsSetDenseJacFn failed.");
}
}
}
}
//====================================================================================================================
void IDA_Solver::setMaxErrTestFailures(int maxErrTestFails) {
m_maxErrTestFails = maxErrTestFails;
}
@ -222,6 +277,9 @@ namespace Cantera {
void IDA_Solver::init(doublereal t0) {
m_t0 = t0;
m_told = t0;
m_told_old = t0;
m_tcurrent = t0;
if (m_y) {
N_VDestroy_Serial(nv(m_y));
}
@ -338,9 +396,15 @@ namespace Cantera {
throw IDA_Err("unsupported linear solver type");
}
if (m_formJac == 1) {
flag = IDADlsSetDenseJacFn(m_ida_mem, ida_jacobian);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDADlsSetDenseJacFn failed.");
}
}
// pass a pointer to func in m_data
m_fdata = new ResidData(&m_resid, m_resid.nparams());
m_fdata = new ResidData(&m_resid, this, m_resid.nparams());
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
flag = IDASetRdata(m_ida_mem, (void*)m_fdata);
if (flag != IDA_SUCCESS) {
@ -383,7 +447,7 @@ namespace Cantera {
throw IDA_Err("IDASetMaxErrTestFails failed.");
}
}
if (m_maxNonlinIters >= 0) {
if (m_maxNonlinIters > 0) {
flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetmaxNonlinIters failed.");
@ -490,11 +554,35 @@ namespace Cantera {
//====================================================================================================================
int IDA_Solver::solve(double tout)
{
double t;
double tretn;
int flag;
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_NORMAL);
if (flag != IDA_SUCCESS)
flag = IDASetStopTime(m_ida_mem, tout);
if (flag != IDA_SUCCESS) {
throw IDA_Err(" IDA error encountered.");
}
do {
if (tout <= m_tcurrent) {
throw IDA_Err(" tout <= tcurrent");
}
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);
if (flag < 0) {
throw IDA_Err(" IDA error encountered.");
} else if (flag == IDA_TSTOP_RETURN) {
// we've reached our goal, and have actually integrated past it
} else if (flag == IDA_ROOT_RETURN) {
// not sure what to do with this yet
} else if (flag == IDA_WARNING) {
throw IDA_Err(" IDA Warning encountered.");
}
m_tcurrent = tretn;
m_deltat = m_tcurrent - m_told;
} while (tretn < tout);
if (flag != IDA_SUCCESS && flag != IDA_TSTOP_RETURN) {
throw IDA_Err(" IDA error encountered.");
}
return flag;
}
//====================================================================================================================
@ -502,9 +590,23 @@ namespace Cantera {
{
double t;
int flag;
if (tout <= m_tcurrent) {
throw IDA_Err(" tout <= tcurrent");
}
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);
if (flag != IDA_SUCCESS)
if (flag < 0) {
throw IDA_Err(" IDA error encountered.");
} else if (flag == IDA_TSTOP_RETURN) {
// we've reached our goal, and have actually integrated past it
} else if (flag == IDA_ROOT_RETURN) {
// not sure what to do with this yet
} else if (flag == IDA_WARNING) {
throw IDA_Err(" IDA Warning encountered.");
}
m_tcurrent = t;
m_deltat = m_tcurrent - m_told;
return t;
}
//====================================================================================================================

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

@ -110,6 +110,17 @@ namespace Cantera {
virtual void setStopTime(doublereal tstop);
//! Set the form of the jacobian
/*!
*
* @param formJac Form of the jacobian
*
* 0 numerical jacobian
* 1 analytical jacobian given by the evalJacobianDP() function
*/
virtual void setJacobianType(int formJac);
virtual void setMaxErrTestFailures(int n);
//! Set the maximum number of nonlinear iterations on a timestep
@ -289,9 +300,30 @@ namespace Cantera {
//! maximum time step order of the method
int m_maxord;
//! Form of the jacobian
/*!
* 0 numerical jacobian created by ida
* 1 analytical jacobian. Must have populated the evalJacobianDP()
* function in the ResidJacEval class.
* 2 numerical jacobian formed by the ResidJacEval class (unimplemented)
*/
int m_formJac;
//! maximum time
doublereal m_tstop;
//! Value of the previous, previous time
doublereal m_told_old;
//! Value of the previous time
doublereal m_told;
//! Value of the current time
doublereal m_tcurrent;
//! Value of deltaT for the current step
doublereal m_deltat;
//! maximum number of error test failures
int m_maxErrTestFails;

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

@ -334,7 +334,35 @@ namespace Cantera {
SquareMatrix &J,
doublereal * const resid)
{
throw CanteraError("ResidJacEval::evalJacobian()", "Not implemented\n");
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);
}
//====================================================================================================================
// Calculate an analytical jacobian and the residual at the current time and values.
/*
* Only called if the jacFormation method is set to analytical
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param c_j The current value of the coefficient of the time derivative
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param jac_colPts Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
int ResidJacEval::
evalJacobianDP(const doublereal t, const doublereal delta_t,
const doublereal c_j,
const doublereal * const y,
const doublereal * const ydot,
doublereal * const * jac_colPts,
doublereal * const resid)
{
throw CanteraError("ResidJacEval::evalJacobianDP()", "Not implemented\n");
return 1;
}
//====================================================================================================================

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

@ -333,6 +333,27 @@ namespace Cantera {
SquareMatrix &J,
doublereal * const resid);
//! Calculate an analytical jacobian and the residual at the current time and values.
/*!
* Only called if the jacFormation method is set to analytical
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @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)
* @param resid Value of the residual that is computed (output)
*
* @return Returns a flag to indicate that operation is successful.
* 1 Means a successful operation
* -0 or neg value Means an unsuccessful operation
*/
virtual int evalJacobianDP(const doublereal t, const doublereal delta_t, doublereal cj,
const doublereal* const y,
const doublereal* const ydot,
doublereal * const *jacobianColPts,
doublereal * const resid);
protected:
//! constant value of atol

2
Cantera/src/thermo/MixtureFugacityTP.h
View file

@ -55,7 +55,7 @@ namespace Cantera {
* steps for efficiently handling mixture of gases that whose standard states
* are defined as ideal gases, but which describe also non-ideal solutions.
* In addition a multicomponent liquid phase below the critical temperature of the
* mixture is also allowed. The main subclass will be a mixture Redlich-kwong class.
* mixture is also allowed. The main subclass is currently a mixture Redlich-Kwong class.
*
* Several concepts are introduced. The first concept is there are temporary
* variables for holding the species standard state values