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

Worked on getting the IDA_Solver interface operable. Previously,

Browse source 
it hadn't been. It now works on one test problem.
This commit is contained in:
Harry Moffat 2011-07-30 03:08:10 +00:00
parent 9d46965e7d
commit 2082b57699
10 changed files with 1009 additions and 492 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

361
Cantera/src/numerics/DAE_Solver.h
View file

@ -12,207 +12,236 @@
* Copyright 2006 California Institute of Technology
*
*/
#undef DAE_DEVEL
#ifndef CT_DAE_Solver_H
#define CT_DAE_Solver_H
#include <vector>
#include "ct_defs.h"
#include "ResidEval.h"
#include "ResidJacEval.h"
#include "global.h"
namespace Cantera {
#define DAE_DEVEL
#ifdef DAE_DEVEL
class Jacobian {
public:
Jacobian(){}
virtual ~Jacobian(){}
virtual bool supplied() { return false; }
virtual bool isBanded() { return false; }
virtual int lowerBandWidth() { return 0; }
virtual int upperBandWidth() { return 0; }
};
class Jacobian {
public:
Jacobian(){}
virtual ~Jacobian(){}
virtual bool supplied() { return false; }
virtual bool isBanded() { return false; }
virtual int lowerBandWidth() { return 0; }
virtual int upperBandWidth() { return 0; }
};
class BandedJacobian : public Jacobian {
public:
BandedJacobian(int ml, int mu) {
m_ml = ml; m_mu = mu;
}
virtual bool supplied() { return false; }
virtual bool isBanded() { return true; }
virtual int lowerBandWidth() { return m_ml; }
virtual int upperBandWidth() { return m_mu; }
protected:
int m_ml, m_mu;
};
class BandedJacobian : public Jacobian {
public:
BandedJacobian(int ml, int mu) {
m_ml = ml; m_mu = mu;
}
virtual bool supplied() { return false; }
virtual bool isBanded() { return true; }
virtual int lowerBandWidth() { return m_ml; }
virtual int upperBandWidth() { return m_mu; }
protected:
int m_ml, m_mu;
};
const int cDirect = 0;
const int cKrylov = 1;
const int cDirect = 0;
const int cKrylov = 1;
/**
* Wrapper for DAE solvers
*/
class DAE_Solver {
public:
DAE_Solver(ResidJacEval& f) :
m_resid(f),
m_neq(f.nEquations()),
m_time(0.0)
{
}
virtual ~DAE_Solver(){}
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a vector absolute tolerance.
*/
virtual void setTolerances(doublereal reltol,
doublereal* abstol) {
warn("setTolerances");
}
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a scalar absolute tolerance.
*/
virtual void setTolerances(doublereal reltol, doublereal abstol) {
warn("setTolerances");
}
/**
* Specify a Jacobian evaluator. If this method is not called,
* the Jacobian will be computed by finite difference.
*/
void setJacobian(Jacobian& jac) {
warn("setJacobian");
}
virtual void setLinearSolverType(int solverType) {
warn("setLinearSolverType");
}
virtual void setDenseLinearSolver() {
warn("setDenseLinearSolver");
}
virtual void setBandedLinearSolver(int m_upper, int m_lower) {
warn("setBandedLinearSolver");
}
virtual void setMaxStepSize(doublereal dtmax) {
warn("setMaxStepSize");
}
virtual void setMaxOrder(int n) {
warn("setMaxOrder");
}
virtual void setMaxNumSteps(int n) {
warn("setMaxNumSteps");
}
virtual void setInitialStepSize(doublereal h0) {
warn("setInitialStepSize");
}
virtual void setStopTime(doublereal tstop) {
warn("setStopTime");
}
virtual void setMaxErrTestFailures(int n) {
warn("setMaxErrTestFailures");
}
virtual void setMaxNonlinIterations(int n) {
warn("setMaxNonlinIterations");
}
virtual void setMaxNonlinConvFailures(int n) {
warn("setMaxNonlinConvFailures");
}
virtual void inclAlgebraicInErrorTest(bool yesno) {
warn("inclAlgebraicInErrorTest");
}
/**
* Wrapper for DAE solvers
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the derivatives
* of all variables, this method computes the initial y
* values.
*/
class DAE_Solver {
public:
virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp,
doublereal tout) {
warn("correctInitial_Y_given_Yp");
}
/**
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the initial
* values of all differential variables, it computes the
* initial values of all algebraic variables and the initial
* derivatives of all differential variables.
*/
virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp,
doublereal tout)
{
warn("correctInitial_YaYp_given_Yd");
}
DAE_Solver(ResidEval& f) : m_resid(f),
m_neq(f.nEquations()),
m_time(0.0) {}
/**
* Solve the system of equations up to time tout.
*/
virtual int solve(doublereal tout) {
warn("solve"); return 0;
}
virtual ~DAE_Solver(){}
/**
* Take one internal step.
*/
virtual doublereal step(doublereal tout) {
warn("step"); return 0;
}
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a vector absolute tolerance.
*/
virtual void setTolerances(doublereal reltol,
doublereal* abstol) {
warn("setTolerances");
}
/// Number of equations.
int nEquations() const {
return m_resid.nEquations();
}
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a scalar absolute tolerance.
*/
virtual void setTolerances(doublereal reltol, doublereal abstol) {
warn("setTolerances");
}
/**
* initialize. Base class method does nothing.
*/
virtual void init(doublereal t0) {}
/**
* Specify a Jacobian evaluator. If this method is not called,
* the Jacobian will be computed by finite difference.
*/
void setJacobian(Jacobian& jac) {
warn("setJacobian");
}
/**
* Set a solver-specific input parameter.
*/
virtual void setInputParameter(int flag, doublereal value) {
warn("setInputParameter");
}
virtual void setLinearSolverType(int solverType) {
warn("setLinearSolverType");
}
/**
* Get the value of a solver-specific output parameter.
*/
virtual doublereal getOutputParameter(int flag) const {
warn("getOutputParameter"); return 0.0;
}
virtual void setDenseLinearSolver() {
warn("setDenseLinearSolver");
}
/// the current value of solution component k.
virtual doublereal solution(int k) const {
warn("solution"); return 0.0;
}
virtual void setBandedLinearSolver(int m_upper, int m_lower) {
warn("setBandedLinearSolver");
}
virtual void setMaxTime(doublereal tmax) {
warn("setMaxTime");
}
virtual void setMaxStepSize(doublereal dtmax) {
warn("setMaxStepSize");
}
virtual void setMaxOrder(int n) {
warn("setMaxOrder");
}
virtual void setMaxNumSteps(int n) {
warn("setMaxNumSteps");
}
virtual void setInitialStepSize(doublereal h0) {
warn("setInitialStepSize");
}
virtual void setStopTime(doublereal tstop) {
warn("setStopTime");
}
virtual void setMaxErrTestFailures(int n) {
warn("setMaxErrTestFailures");
}
virtual void setMaxNonlinIterations(int n) {
warn("setMaxNonlinIterations");
}
virtual void setMaxNonlinConvFailures(int n) {
warn("setMaxNonlinConvFailures");
}
virtual void inclAlgebraicInErrorTest(bool yesno) {
warn("inclAlgebraicInErrorTest");
}
virtual const doublereal* solutionVector() const {
warn("solutionVector"); return &m_dummy;
}
virtual void correctInitial_Y_given_Yp() {
warn("correctInitial_Y_given_Yp");
}
/// the current value of the derivative of solution component k.
virtual doublereal derivative(int k) const {
warn("derivative"); return 0.0;
}
virtual void correctInitial_YaYp_given_Yd() {
warn("correctInitial_YaYp_given_Yd");
}
virtual const doublereal* derivativeVector() const {
warn("derivativeVector"); return &m_dummy;
}
/**
* Solve the system of equations up to time tout.
*/
virtual int solve(doublereal tout) {
warn("solve"); return 0;
}
protected:
/**
* Take one internal step.
*/
virtual doublereal step(doublereal tout) {
warn("step"); return 0;
}
doublereal m_dummy;
/// Number of equations.
int nEquations() const { return m_resid.nEquations(); }
ResidJacEval& m_resid;
/**
* initialize. Base class method does nothing.
*/
virtual void init(doublereal t0) {}
/**
* Set a solver-specific input parameter.
*/
virtual void setInputParameter(int flag, doublereal value) {
warn("setInputParameter");
}
/**
* Get the value of a solver-specific output parameter.
*/
virtual doublereal getOutputParameter(int flag) const {
warn("getOutputParameter"); return 0.0;
}
/// the current value of solution component k.
virtual doublereal solution(int k) const {
warn("solution"); return 0.0;
}
virtual const doublereal* solutionVector() const {
warn("solutionVector"); return &m_dummy;
}
/// the current value of the derivative of solution component k.
virtual doublereal derivative(int k) const {
warn("derivative"); return 0.0;
}
virtual const doublereal* derivativeVector() const {
warn("derivativeVector"); return &m_dummy;
}
protected:
doublereal m_dummy;
ResidEval& m_resid;
integer m_neq;
doublereal m_time;
//! Number of total equations in the system
integer m_neq;
doublereal m_time;
private:
void warn(std::string msg) const {
writelog(">>>> Warning: method "+msg+" of base class "
+"DAE_Solver called. Nothing done.\n");
}
};
private:
void warn(std::string msg) const {
writelog(">>>> Warning: method "+msg+" of base class "
+"DAE_Solver called. Nothing done.\n");
}
};
//! Factor method for choosing a DAE solver
/*!
*
* @param itype String identifying the type
* (IDA is the only option)
* @param f Residual function to be solved by the DAE algorithm
*
* @return Returns a point to the instantiated DAE_Solver object
*/
DAE_Solver* newDAE_Solver(std::string itype, ResidJacEval& f);
#endif

9
Cantera/src/numerics/DAE_solvers.cpp
View file

@ -3,6 +3,7 @@
#include "DAE_Solver.h"
// DAE_DEVEL is turned off at the current time
#define DAE_DEVEL
#ifdef DAE_DEVEL
#ifdef HAS_SUNDIALS
@ -11,18 +12,18 @@
namespace Cantera {
DAE_Solver* newDAE_Solver(string itype) {
DAE_Solver* newDAE_Solver(string itype, ResidJacEval& f) {
if (itype == "IDA") {
#ifdef HAS_SUNDIALS
return new IDA_Solver();
return new IDA_Solver(f);
#else
raise CanteraError("newDAE_Solver","IDA solver requires sundials"
" package, but Cantera was not built with sundials.");
" package, but Cantera was not built with sundials.");
#endif
}
else {
throw CanteraError("newDAE_Solver",
"unknown DAE solver: "+itype);
"unknown DAE solver: "+itype);
}
}
}

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

@ -55,7 +55,7 @@ namespace Cantera {
*/
virtual int neq()=0;
/// Number of parameters.
//! Number of parameters.
virtual int nparams() { return 0; }
protected:

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

@ -1,4 +1,3 @@
/**
* @file IDA_Solver.cpp
*
@ -10,8 +9,8 @@
#include "stringUtils.h"
#include <iostream>
using namespace std;
#ifdef SUNDIALS_VERSION_22
#include <sundials_types.h>
#include <sundials_math.h>
#include <ida.h>
@ -19,238 +18,506 @@ using namespace std;
#include <ida_spgmr.h>
#include <ida_band.h>
#include <nvector_serial.h>
#else
#include <sundials/sundials_types.h>
#include <sundials/sundials_math.h>
#include <ida/ida.h>
#include <ida/ida_dense.h>
#include <ida/ida_spgmr.h>
#include <ida/ida_band.h>
#include <nvector/nvector_serial.h>
#endif
using namespace std;
inline static N_Vector nv(void* x) {
return reinterpret_cast<N_Vector>(x);
return reinterpret_cast<N_Vector>(x);
}
namespace Cantera {
/**
* A simple class to hold an array of parameter values and a pointer to
* an instance of a subclass of ResidEval.
*/
class ResidData {
/**
* A simple class to hold an array of parameter values and a pointer to
* an instance of a subclass of ResidEval.
*/
class ResidData {
public:
ResidData(ResidJacEval* f, int npar = 0) {
m_func = f;
}
virtual ~ResidData() {
}
public:
ResidData(ResidEval* f, int npar = 0) {
m_func = f;
}
virtual ~ResidData() {}
ResidEval* m_func;
};
ResidJacEval* m_func;
};
}
//======================================================================================================================
extern "C" {
/**
* Function called by IDA to evaluate the residual, given y and
* ydot. IDA allows passing in a void* pointer to access
* external data. Instead of requiring the user to provide a
* residual function directly to IDA (which would require using
* 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.
*/
static int ida_resid(realtype t, N_Vector y, N_Vector ydot,
N_Vector r, void *f_data) {
double* ydata = NV_DATA_S(y);
double* ydotdata = NV_DATA_S(ydot);
double* rdata = NV_DATA_S(r);
Cantera::ResidData* d = (Cantera::ResidData*)f_data;
Cantera::ResidEval* f = d->m_func;
f->eval(t, ydata, ydotdata, rdata);
return 0;
}
//! Function called by IDA to evaluate the residual, given y and ydot.
/*!
* IDA allows passing in a void* pointer to access external data. Instead of requiring the user to provide a
* residual function directly to IDA (which would require using
* 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.
*/
static int ida_resid(realtype t, N_Vector y, N_Vector ydot, N_Vector r, void *f_data) {
double* ydata = NV_DATA_S(y);
double* ydotdata = NV_DATA_S(ydot);
double* rdata = NV_DATA_S(r);
Cantera::ResidData* d = (Cantera::ResidData*) f_data;
Cantera::ResidJacEval* f = d->m_func;
f->eval(t, ydata, ydotdata, rdata);
return 0;
}
}
namespace Cantera {
/**
* Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration.
*/
IDA_Solver::IDA_Solver(ResidEval& f) : DAE_Solver(f),
m_neq(0),
m_ida_mem(0),
m_t0(0.0),
m_y(0),
m_ydot(0),
m_abstol(0),
m_type(0),
m_itol(IDA_SS),
m_iter(0),
m_maxord(0),
m_reltol(1.e-9),
m_abstols(1.e-15),
m_nabs(0),
m_hmax(0.0),
m_maxsteps(20000),
m_mupper(0),
m_mlower(0) {}
/// Destructor.
IDA_Solver::~IDA_Solver()
{
if (m_ida_mem) {
IDAFree(&m_ida_mem);
}
if (m_y) N_VDestroy_Serial(nv(m_y));
if (m_ydot) N_VDestroy_Serial(nv(m_ydot));
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
delete m_fdata;
//====================================================================================================================
/*
* Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration.
*/
IDA_Solver::IDA_Solver(ResidJacEval& f) :
DAE_Solver(f),
m_ida_mem(0),
m_t0(0.0),
m_y(0),
m_ydot(0),
m_id(0),
m_constraints(0),
m_abstol(0),
m_type(0),
m_itol(IDA_SS),
m_iter(0),
m_reltol(1.e-9),
m_abstols(1.e-15),
m_nabs(0),
m_hmax(0.0),
m_hmin(0.0),
m_h0(0.0),
m_maxsteps(20000),
m_maxord(0),
m_tstop(0.0),
m_maxErrTestFails(-1),
m_maxNonlinIters(0),
m_maxNonlinConvFails(-1),
m_setSuppressAlg(0),
m_fdata(0),
m_mupper(0),
m_mlower(0)
{
}
//====================================================================================================================
IDA_Solver::~IDA_Solver()
{
if (m_ida_mem) {
IDAFree(&m_ida_mem);
}
doublereal IDA_Solver::solution(int k) const {
return NV_Ith_S(nv(m_y),k);
if (m_y) N_VDestroy_Serial(nv(m_y));
if (m_ydot) N_VDestroy_Serial(nv(m_ydot));
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
if (m_constraints) N_VDestroy_Serial(nv(m_constraints));
delete m_fdata;
}
//====================================================================================================================
doublereal IDA_Solver::solution(int k) const {
return NV_Ith_S(nv(m_y),k);
}
//====================================================================================================================
const doublereal* IDA_Solver::solutionVector() const {
return NV_DATA_S(nv(m_y));
}
//====================================================================================================================
doublereal IDA_Solver::derivative(int k) const {
return NV_Ith_S(nv(m_ydot),k);
}
//====================================================================================================================
const doublereal* IDA_Solver::derivativeVector() const {
return NV_DATA_S(nv(m_ydot));
}
//====================================================================================================================
void IDA_Solver::setTolerances(double reltol, double* abstol) {
m_itol = IDA_SV;
if (!m_abstol) {
m_abstol = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
}
const doublereal* IDA_Solver::solutionVector() const { return NV_DATA_S(nv(m_y));}
doublereal IDA_Solver::derivative(int k) const {
return NV_Ith_S(nv(m_ydot),k);
for (int i = 0; i < m_neq; i++) {
NV_Ith_S(nv(m_abstol), i) = abstol[i];
}
const doublereal* IDA_Solver::derivativeVector() const { return NV_DATA_S(nv(m_ydot));}
m_reltol = reltol;
int flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol));
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
}
}
//====================================================================================================================
void IDA_Solver::setTolerances(doublereal reltol, doublereal abstol) {
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.");
}
}
//====================================================================================================================
void IDA_Solver::setLinearSolverType(int solverType) {
m_type = solverType;
}
//====================================================================================================================
void IDA_Solver::setDenseLinearSolver() {
setLinearSolverType(0);
}
//====================================================================================================================
void IDA_Solver::setBandedLinearSolver(int m_upper, int m_lower) {
m_type = 2;
m_upper = m_mupper;
m_mlower = m_lower;
}
//====================================================================================================================
void IDA_Solver::setMaxOrder(int n) {
m_maxord = n;
}
//====================================================================================================================
void IDA_Solver::setMaxNumSteps(int n) {
m_maxsteps = n;
}
//====================================================================================================================
void IDA_Solver::setInitialStepSize(doublereal h0) {
m_h0 = h0;
}
//====================================================================================================================
void IDA_Solver::setStopTime(doublereal tstop) {
m_tstop = tstop;
}
//====================================================================================================================
void IDA_Solver::setMaxErrTestFailures(int maxErrTestFails) {
m_maxErrTestFails = maxErrTestFails;
}
//====================================================================================================================
void IDA_Solver::setMaxNonlinIterations(int n) {
m_maxNonlinIters = n;
}
//====================================================================================================================
void IDA_Solver::setMaxNonlinConvFailures(int n) {
m_maxNonlinConvFails = n;
}
//====================================================================================================================
void IDA_Solver::inclAlgebraicInErrorTest(bool yesno) {
if (yesno) {
m_setSuppressAlg = 0;
} else {
m_setSuppressAlg = 1;
}
}
//====================================================================================================================
void IDA_Solver::init(doublereal t0) {
void IDA_Solver::setTolerances(double reltol, double* abstol) {
m_itol = IDA_SV;
if (m_abstol) N_VDestroy_Serial(nv(m_abstol));
m_abstol = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
for (int i=0; i < m_neq; i++) {
NV_Ith_S(nv(m_abstol), i) = abstol[i];
}
m_reltol = reltol;
m_t0 = t0;
if (m_y) {
N_VDestroy_Serial(nv(m_y));
}
if (m_ydot) N_VDestroy_Serial(nv(m_ydot));
if (m_id) N_VDestroy_Serial(nv(m_id));
if (m_constraints) N_VDestroy_Serial(nv(m_constraints));
m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
m_ydot = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
m_constraints = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
for (int i=0; i<m_neq; i++) {
NV_Ith_S(nv(m_y), i) = 0.0;
NV_Ith_S(nv(m_ydot), i) = 0.0;
NV_Ith_S(nv(m_constraints), i) = 0.0;
}
// get the initial conditions
m_resid.getInitialConditions(m_t0, NV_DATA_S(nv(m_y)), NV_DATA_S(nv(m_ydot)));
if (m_ida_mem) {
IDAFree(&m_ida_mem);
}
void IDA_Solver::setTolerances(double reltol, double abstol) {
m_itol = IDA_SS;
m_reltol = reltol;
m_abstols = abstol;
}
void IDA_Solver::setLinearSolverType(int solverType) {
m_type = solverType;
}
void IDA_Solver::init(double t0)
{
m_t0 = t0;
if (m_y) N_VDestroy_Serial(nv(m_y));
if (m_ydot) N_VDestroy_Serial(nv(m_ydot));
if (m_id) N_VDestroy_Serial(nv(m_id));
if (m_constraints) N_VDestroy_Serial(nv(m_constraints));
m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
m_ydot = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
m_constraints = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
for (int i=0; i<m_neq; i++) {
NV_Ith_S(nv(m_y), i) = 0.0;
NV_Ith_S(nv(m_ydot), i) = 0.0;
NV_Ith_S(nv(m_constraints), i) = 0.0;
}
// get the initial conditions
m_resid.getInitialConditions(m_t0, NV_DATA_S(nv(m_ydot)),
NV_DATA_S(nv(m_y)));
if (m_ida_mem) IDAFree(&m_ida_mem);
m_ida_mem = IDACreate();
int flag = 0;
if (m_itol == IDA_SV) {
// vector atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, nv(m_abstol));
}
else {
// scalar atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, &m_abstols);
}
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
throw IDA_Err("Memory allocation failed."); }
else if (flag == IDA_ILL_INPUT) {
throw IDA_Err("Illegal value for IDAMalloc input argument.");
}
else
throw IDA_Err("IDAMalloc failed.");
}
//-----------------------------------
// set the linear solver type
//-----------------------------------
if (m_type == 1) {
long int N = m_neq;
IDADense(m_ida_mem, N);
}
else if (m_type == 2) {
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
IDABand(m_ida_mem, N, nu, nl);
}
else {
throw IDA_Err("unsupported linear solver type");
}
// pass a pointer to func in m_data
m_fdata = new ResidData(&func, func.nparams());
flag = IDASetRdata(m_ida_mem, (void*)m_fdata);
if (flag != IDA_SUCCESS)
throw IDA_Err("IDASetRdata failed.");
// set options
//if (m_maxord > 0)
// flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
//if (m_maxsteps > 0)
// flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
//if (m_hmax > 0)
// flag = CVodeSetMaxStep(m_cvode_mem, m_hmax);
}
void IDA_Solver::solve(double tout)
{
double t;
int flag;
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_NORMAL);
if (flag != IDA_SUCCESS)
throw IDA_Err(" IDA error encountered.");
}
double IDA_Solver::step(double tout)
{
double t;
int flag;
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP);
if (flag != IDA_SUCCESS)
throw IDA_Err(" IDA error encountered.");
return t;
/* Call IDACreate */
m_ida_mem = IDACreate();
int flag = 0;
if (m_itol == IDA_SV) {
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
// vector atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, nv(m_abstol));
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
throw IDA_Err("Memory allocation failed.");
} else if (flag == IDA_ILL_INPUT) {
throw IDA_Err("Illegal value for IDAMalloc input argument.");
} else
throw IDA_Err("IDAMalloc failed.");
}
doublereal IDA_Solver::getOutputParameter(int flag) {
switch (flag) {
case REAL_WORKSPACE_SIZE:
flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw);
return doublereal(lenrw);
}
#elif defined(SUNDIALS_VERSION_24)
flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
throw IDA_Err("Memory allocation failed.");
} else if (flag == IDA_ILL_INPUT) {
throw IDA_Err("Illegal value for IDAMalloc input argument.");
}
else
throw IDA_Err("IDAMalloc failed.");
}
flag = IDASVtolerances(m_ida_mem, m_reltol, nv(m_abstol));
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
}
#endif
}
else {
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
// scalar atol
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
m_itol, m_reltol, &m_abstols);
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
throw IDA_Err("Memory allocation failed."); }
else if (flag == IDA_ILL_INPUT) {
throw IDA_Err("Illegal value for IDAMalloc input argument.");
}
else
throw IDA_Err("IDAMalloc failed.");
}
#elif defined(SUNDIALS_VERSION_24)
flag = IDAInit(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
if (flag == IDA_MEM_FAIL) {
throw IDA_Err("Memory allocation failed."); }
else if (flag == IDA_ILL_INPUT) {
throw IDA_Err("Illegal value for IDAMalloc input argument.");
}
else
throw IDA_Err("IDAMalloc failed.");
}
flag = IDASStolerances(m_ida_mem, m_reltol, m_abstols);
if (flag != IDA_SUCCESS) {
throw IDA_Err("Memory allocation failed.");
}
#endif
}
//-----------------------------------
// set the linear solver type
//-----------------------------------
if (m_type == 1 || m_type == 0) {
long int N = m_neq;
flag = IDADense(m_ida_mem, N);
if (flag) {
throw IDA_Err("IDADense failed");
}
}
else if (m_type == 2) {
long int N = m_neq;
long int nu = m_mupper;
long int nl = m_mlower;
IDABand(m_ida_mem, N, nu, nl);
}
else {
throw IDA_Err("unsupported linear solver type");
}
// pass a pointer to func in m_data
m_fdata = new ResidData(&m_resid, m_resid.nparams());
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION_23)
flag = IDASetRdata(m_ida_mem, (void*)m_fdata);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetRdata failed.");
}
#elif defined(SUNDIALS_VERSION_24)
flag = IDASetUserData(m_ida_mem, (void*)m_fdata);
if (flag != IDA_SUCCESS)
throw IDA_Err("IDASetUserData failed.");
#endif
// set options
if (m_maxord > 0) {
flag = IDASetMaxOrd(m_ida_mem, m_maxord);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetMaxOrd failed.");
}
}
if (m_maxsteps > 0) {
flag = IDASetMaxNumSteps(m_ida_mem, m_maxsteps);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetMaxNumSteps failed.");
}
}
if (m_h0 > 0.0) {
flag = IDASetInitStep(m_ida_mem, m_h0);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetInitStep failed.");
}
}
if (m_tstop > 0.0) {
flag = IDASetStopTime(m_ida_mem, m_tstop);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetStopTime failed.");
}
}
if (m_maxErrTestFails >= 0) {
flag = IDASetMaxErrTestFails(m_ida_mem, m_maxErrTestFails);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetMaxErrTestFails failed.");
}
}
if (m_maxNonlinIters >= 0) {
flag = IDASetMaxNonlinIters(m_ida_mem, m_maxNonlinIters);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetmaxNonlinIters failed.");
}
}
if (m_maxNonlinConvFails >= 0) {
flag = IDASetMaxConvFails(m_ida_mem, m_maxNonlinConvFails);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetMaxConvFails failed.");
}
}
if (m_setSuppressAlg != 0) {
flag = IDASetSuppressAlg(m_ida_mem, m_setSuppressAlg);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDASetSuppressAlg failed.");
}
}
}
//====================================================================================================================
// Calculate consistent value of the starting solution given the starting solution derivatives
/*
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the derivatives
* of all variables, this method computes the initial y
* values.
*/
void IDA_Solver::correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, doublereal tout) {
int icopt = IDA_Y_INIT;
doublereal tout1 = tout;
if (tout == 0.0) {
double h0 = 1.0E-5;
if (m_h0 > 0.0) {
h0 = m_h0;
}
tout1 = m_t0 + h0;
}
int flag = IDACalcIC(m_ida_mem, icopt, tout1);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDACalcIC failed: error = " + int2str(flag));
}
flag = IDAGetSolution(m_ida_mem, tout1, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag));
}
doublereal *yy = NV_DATA_S(nv(m_y));
doublereal *yyp = NV_DATA_S(nv(m_ydot));
for (int i = 0; i < m_neq; i++) {
y[i] = yy[i];
yp[i] = yyp[i];
}
}
//====================================================================================================================
/*
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the initial
* values of all differential variables, it computes the
* initial values of all algebraic variables and the initial
* derivatives of all differential variables.
*
* @param y Calculated value of the solution vector after the procedure ends
* @param yp Calculated value of the solution derivative after the procedure
* @param The first value of t at which a soluton will be
* requested (from IDASolve). (This is needed here to
* determine the direction of integration and rough scale
* in the independent variable t.
*/
void IDA_Solver::correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, doublereal tout) {
int icopt = IDA_YA_YDP_INIT;
doublereal tout1 = tout;
if (tout == 0.0) {
double h0 = 1.0E-5;
if (m_h0 > 0.0) {
h0 = m_h0;
}
tout1 = m_t0 + h0;
}
int flag = IDACalcIC(m_ida_mem, icopt, tout1);
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDACalcIC failed: error = " + int2str(flag));
}
flag = IDAGetSolution(m_ida_mem, tout1, nv(m_y), nv(m_ydot));
if (flag != IDA_SUCCESS) {
throw IDA_Err("IDAGetSolution failed: error = " + int2str(flag));
}
doublereal *yy = NV_DATA_S(nv(m_y));
doublereal *yyp = NV_DATA_S(nv(m_ydot));
for (int i = 0; i < m_neq; i++) {
y[i] = yy[i];
yp[i] = yyp[i];
}
}
//====================================================================================================================
int IDA_Solver::solve(double tout)
{
double t;
int flag;
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_NORMAL);
if (flag != IDA_SUCCESS)
throw IDA_Err(" IDA error encountered.");
return flag;
}
//====================================================================================================================
double IDA_Solver::step(double tout)
{
double t;
int flag;
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP);
if (flag != IDA_SUCCESS)
throw IDA_Err(" IDA error encountered.");
return t;
}
//====================================================================================================================
doublereal IDA_Solver::getOutputParameter(int flag) const {
long int lenrw, leniw;
switch (flag) {
case REAL_WORKSPACE_SIZE:
flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw);
return doublereal(lenrw);
break;
}
return 0.0;
}
//====================================================================================================================
}

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

@ -13,121 +13,304 @@
*
*/
#ifndef CT_IDA_Solver_H
#define CT_IDA_Solver_H
#ifndef CT_IDA_SOLVER_H
#define CT_IDA_SOLVER_H
#include <vector>
#include "DAE_Solver.h"
#include "ctexceptions.h"
#ifdef SUNDIALS_VERSION_22
#include <nvector_serial.h>
#else
#include <sundials/sundials_nvector.h>
// These constants are defined internally in the ida package, ida.c
#define IDA_NN 0
#define IDA_SS 1
#define IDA_SV 2
#define IDA_WF 3
#endif
#if defined(SUNDIALS_VERSION_24)
#define REAL_WORKSPACE_SIZE 0
#endif
namespace Cantera {
/**
* Exception thrown when a IDA error is encountered.
/**
* Exception thrown when a IDA error is encountered.
*/
class IDA_Err : public CanteraError {
public:
IDA_Err(std::string msg) : CanteraError("IDA_Solver", msg){}
};
class ResidData; // forward reference
class IDA_Solver : public DAE_Solver {
public:
//! Constructor.
/*!
* Default settings: dense jacobian, no user-supplied Jacobian function, Newton iteration.
*
* @param f Function that will supply the time dependent residual to be solved
*/
class IDA_Err : public CanteraError {
public:
IDA_Err(std::string msg) : CanteraError("IDA_Solver", msg){}
};
IDA_Solver(ResidJacEval& f);
virtual ~IDA_Solver();
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a vector absolute tolerance.
*/
virtual void setTolerances(doublereal reltol,
doublereal* abstol);
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a scalar absolute tolerance.
*/
virtual void setTolerances(doublereal reltol, doublereal abstol);
virtual void setLinearSolverType(int solverType);
//! Set up the problem to use a dense linear direct solver
virtual void setDenseLinearSolver();
//! Set up the problem to use a band solver
/*!
* @param m_upper upper band width of the matrix
* @param m_lower lower band width of the matrix
*/
virtual void setBandedLinearSolver(int m_upper, int m_lower);
virtual void setMaxOrder(int n);
//! Set the maximum number of time steps
/*!
* @param n input of maximum number of time steps
*/
virtual void setMaxNumSteps(int n);
//! Sset the initial step size
/*!
* @param h0 initial step size value
*/
virtual void setInitialStepSize(doublereal h0);
//! Set the stop time
/*!
* @param tstop the independent variable value past which the solution is not to proceed.
*/
virtual void setStopTime(doublereal tstop);
class ResidData; // forward reference
virtual void setMaxErrTestFailures(int n);
class IDA_Solver : public DAE_Solver {
public:
//! Set the maximum number of nonlinear iterations on a timestep
/*!
* @param n Set the max iterations. The default is 4, which seems awefully low to me.
*/
virtual void setMaxNonlinIterations(int n);
IDA_Solver(ResidEval& f);
virtual ~IDA_Solver();
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a vector absolute tolerance.
*/
virtual void setTolerances(doublereal reltol,
doublereal* abstol);
/**
* Set error tolerances. This version specifies a scalar
* relative tolerance, and a scalar absolute tolerance.
*/
virtual void setTolerances(doublereal reltol, doublereal abstol);
virtual void setLinearSolverType(int solverType);
virtual void setDenseLinearSolver();
virtual void setBandedLinearSolver(int m_upper, int m_lower);
virtual void setMaxTime(doublereal tmax);
virtual void setMaxOrder(int n);
virtual void setMaxNumSteps(int n);
virtual void setInitialStepSize(doublereal h0);
virtual void setStopTime(doublereal tstop);
virtual void setMaxErrTestFailures(int n);
virtual void setMaxNonlinIterations(int n);
virtual void setMaxNonlinConvFailures(int n);
virtual void inclAlgebraicInErrorTest(bool yesno);
virtual void setInputParameter(int flag, doublereal value);
virtual doublereal getOutputParameter(int flag);
//! Set the maximum number of nonlinear solver convergence failures
/*!
* @param n Value of nonlin failures. If value is exceeded, the calculation terminates.
*/
virtual void setMaxNonlinConvFailures(int n);
/**
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the derivatives
* of all variables, this method computes the initial y
* values.
*/
virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp,
doublereal tout);
virtual void inclAlgebraicInErrorTest(bool yesno);
/**
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the initial
* values of all differential variables, it computes the
* initial values of all algebraic variables and the initial
* derivatives of all differential variables.
*/
virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp,
doublereal tout);
/**
* Get the value of a solver-specific output parameter.
*/
virtual doublereal getOutputParameter(int flag) const;
//! Calculate consistent value of the starting solution given the starting solution derivatives
/*!
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the derivatives
* of all variables, this method computes the initial y
* values.
*/
virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp,
doublereal tout);
//! Calculate consistent value of the algebraic constraints and derivatives at the start of the problem
/*!
* This method may be called if the initial conditions do not
* satisfy the residual equation F = 0. Given the initial
* values of all differential variables, it computes the
* initial values of all algebraic variables and the initial
* derivatives of all differential variables.
*/
virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, doublereal tout);
//! Step the system to a final value of the time
/*!
* @param tout Final value of the time
*
* @return Returns the IDASolve() return flag
*
* The return values for IDASolve are described below.
* (The numerical return values are defined above in this file.)
* All unsuccessful returns give a negative return value.
*
* IDA_SUCCESS
* IDASolve succeeded and no roots were found.
*
* IDA_ROOT_RETURN: IDASolve succeeded, and found one or more roots.
* If nrtfn > 1, call IDAGetRootInfo to see which g_i were found
* to have a root at (*tret).
*
* IDA_TSTOP_RETURN:
* IDASolve returns computed results for the independent variable
* value tstop. That is, tstop was reached.
*
* IDA_MEM_NULL:
* The IDA_mem argument was NULL.
*
* IDA_ILL_INPUT:
* One of the inputs to IDASolve is illegal. This includes the
* situation when a component of the error weight vectors
* becomes < 0 during internal stepping. It also includes the
* situation where a root of one of the root functions was found
* both at t0 and very near t0. The ILL_INPUT flag
* will also be returned if the linear solver function IDA---
* (called by the user after calling IDACreate) failed to set one
* of the linear solver-related fields in ida_mem or if the linear
* solver's init routine failed. In any case, the user should see
* the printed error message for more details.
*
*
* IDA_TOO_MUCH_WORK:
* The solver took mxstep internal steps but could not reach tout.
* The default value for mxstep is MXSTEP_DEFAULT = 500.
*
* IDA_TOO_MUCH_ACC:
* The solver could not satisfy the accuracy demanded by the user
* for some internal step.
*
* IDA_ERR_FAIL:
* Error test failures occurred too many times (=MXETF = 10) during
* one internal step.
*
* IDA_CONV_FAIL:
* Convergence test failures occurred too many times (= MXNCF = 10)
* during one internal step.
*
* IDA_LSETUP_FAIL:
* The linear solver's setup routine failed
* in an unrecoverable manner.
*
* IDA_LSOLVE_FAIL:
* The linear solver's solve routine failed
* in an unrecoverable manner.
*
* IDA_CONSTR_FAIL:
* The inequality constraints were violated,
* and the solver was unable to recover.
*
* IDA_REP_RES_ERR:
* The user's residual function repeatedly returned a recoverable
* error flag, but the solver was unable to recover.
*
* IDA_RES_FAIL:
* The user's residual function returned a nonrecoverable error
* flag.
*
*/
virtual int solve(doublereal tout);
virtual doublereal step(doublereal tout);
virtual void init(doublereal t0);
//! the current value of solution component k.
/*!
* @param k index of the solution
*/
virtual doublereal solution(int k) const;
virtual const doublereal* solutionVector() const;
//! the current value of the derivative of solution component k.
virtual doublereal derivative(int k) const;
virtual const doublereal* derivativeVector() const;
void *IDAMemory() {
return m_ida_mem;
}
protected:
//! Pointer to the IDA memory for the problem
void* m_ida_mem;
//! Initial value of the time
doublereal m_t0;
//! Current value of the solution vector
void *m_y;
//! Current value of the derivative of the solution vector
void *m_ydot;
void *m_id;
void *m_constraints;
void *m_abstol;
int m_type;
virtual int solve(doublereal tout);
int m_itol;
int m_iter;
doublereal m_reltol;
doublereal m_abstols;
int m_nabs;
virtual doublereal step(doublereal tout);
//! Maximum value of the timestep allowed
doublereal m_hmax;
virtual void init(doublereal t0);
//! Minimum value of the timestep allowd
doublereal m_hmin;
/// the current value of solution component k.
virtual doublereal solution(int k) const;
//! Value of the initial time step
doublereal m_h0;
virtual const doublereal* solutionVector() const;
//! Maximum number of time steps allowed
int m_maxsteps;
/// the current value of the derivative of solution component k.
virtual doublereal derivative(int k) const;
//! maximum time step order of the method
int m_maxord;
virtual const doublereal* derivativeVector() const;
//! maximum time
doublereal m_tstop;
protected:
//! maximum number of error test failures
int m_maxErrTestFails;
int m_neq;
void* m_ida_mem;
doublereal m_t0;
void *m_y, *m_ydot, *m_id, *m_constraints, *m_abstol;
int m_type;
int m_itol;
int m_iter;
doublereal m_reltol;
doublereal m_abstols;
int m_nabs;
doublereal m_hmax, m_hmin;
int m_maxsteps, m_maxord;
ResidData* m_fdata;
int m_mupper, m_mlower;
};
//! Maximum number of nonlinear solver iterations at one solution
/*!
* If zero, this is the default of 4.
*/
int m_maxNonlinIters;
//! Maximum number of nonlinear convergence failures
int m_maxNonlinConvFails;
//! If true, the algebraic variables don't contribute to error tolerances
int m_setSuppressAlg;
ResidData* m_fdata;
int m_mupper;
int m_mlower;
};
}

45
Cantera/src/numerics/Integrator.h
View file

@ -47,18 +47,20 @@ namespace Cantera {
Adams_Method /**< Adams */
};
/**
* Specifies the method used for iteration.
//! Specifies the method used for iteration.
/*!
* Not all methods are supported by all integrators.
*/
enum IterType {
Newton_Iter, /**< Newton iteration */
Functional_Iter /**< Functional iteration */
//! Newton Iteration
Newton_Iter,
//! Functional Iteration
Functional_Iter
};
/**
* Abstract base class for ODE system integrators.
//! Abstract base class for ODE system integrators.
/*!
* @ingroup odeGroup
*/
class Integrator {
@ -77,8 +79,8 @@ namespace Cantera {
/** Set or reset the number of equations. */
//virtual void resize(int n)=0;
/**
* Set error tolerances.
//! Set error tolerances.
/*!
* @param reltol scalar relative tolerance
* @param number of equations
* @param abstol array of N absolute tolerance values
@ -96,13 +98,20 @@ namespace Cantera {
virtual void setTolerances(doublereal reltol, doublereal abstol)
{ warn("setTolerances"); }
virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol)
{}// { warn("setSensitivityTolerances"); }
/**
* Set problem type.
//! Set the sensitvity error tolerances
/*!
* @param reltol scalar relative tolerance
* @param abstol scalar absolute tolerance
*/
virtual void setProblemType(int probtype) { warn("setProblemType"); }
virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol)
{ }
//! Set the problem type.
/*!
* @param probtype Type of the problem
*/
virtual void setProblemType(int probtype)
{ warn("setProblemType"); }
/**
* Initialize the integrator for a new problem. Call after
@ -116,10 +125,10 @@ namespace Cantera {
virtual void reinitialize(doublereal t0, FuncEval& func)
{ warn("reinitialize"); }
/**
* Integrate the system of equations.
* @param tout integrate to this time. Note that this is the
* absolute time value, not a time interval.
//! Integrate the system of equations.
/*!
* @param tout Integrate to this time. Note that this is the
* absolute time value, not a time interval.
*/
virtual void integrate(doublereal tout)
{ warn("integrate"); }

4
Cantera/src/numerics/Makefile.in
View file

@ -37,14 +37,14 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
NUMERICS_OBJ = DenseMatrix.o funcs.o Func1.o \
ODE_integrators.o BandMatrix.o DAE_solvers.o \
funcs.o sort.o SquareMatrix.o ResidJacEval.o NonlinearSolver.o \
solveProb.o BEulerInt.o RootFind.o
solveProb.o BEulerInt.o RootFind.o IDA_Solver.o
NUMERICS_H = ArrayViewer.h DenseMatrix.h \
funcs.h ctlapack.h Func1.h FuncEval.h \
polyfit.h\
BandMatrix.h Integrator.h DAE_Solver.h ResidEval.h sort.h \
SquareMatrix.h ResidJacEval.h NonlinearSolver.h \
solveProb.h BEulerInt.h RootFind.h
solveProb.h BEulerInt.h RootFind.h IDA_Solver.h
ifeq ($(use_sundials), 1)
ODEPACKAGE_H = CVodesIntegrator.h

9
Cantera/src/numerics/ResidEval.h
View file

@ -155,7 +155,14 @@ namespace Cantera {
}
}
//! Return the number of parameters in the calculation
/*!
* This is the number of parameters in the sensitivity calculation. We have
* set this to zero and have included it for later expansion
*/
int nparams () const {
return 0;
}
protected:

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

@ -309,6 +309,13 @@ namespace Cantera {
return 1;
}
//====================================================================================================================
int ResidJacEval::eval(const doublereal t, const doublereal * const y, const doublereal * const ydot,
doublereal * const r) {
double deltaT = -1.0;
int flag = evalResidNJ(t, deltaT, y, ydot, r);
return flag;
}
//====================================================================================================================
// Calculate an analytical jacobian and the residual at the current time and values.
/*
* Only called if the jacFormation method is set to analytical

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

@ -125,6 +125,20 @@ namespace Cantera {
const doublereal delta_x = 0.0);
/**
* Evaluate the residual function. Called by the
* integrator.
* @param t time. (input)
* @param y solution vector. (input)
* @param ydot rate of change of solution vector. (input)
* @param r residual vector (output)
*/
virtual int eval(const doublereal t, const doublereal * const y,
const doublereal * const ydot,
doublereal * const r);
//! Fill in the initial conditions
/*!
* Values for both the solution and the value of ydot may be provided.