From 2082b576994b36e7c6211d693b7a8bdc70290ced Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 30 Jul 2011 03:08:10 +0000 Subject: [PATCH] Worked on getting the IDA_Solver interface operable. Previously, it hadn't been. It now works on one test problem. --- Cantera/src/numerics/DAE_Solver.h | 361 +++++++------- Cantera/src/numerics/DAE_solvers.cpp | 9 +- Cantera/src/numerics/FuncEval.h | 2 +- Cantera/src/numerics/IDA_Solver.cpp | 689 ++++++++++++++++++-------- Cantera/src/numerics/IDA_Solver.h | 361 ++++++++++---- Cantera/src/numerics/Integrator.h | 45 +- Cantera/src/numerics/Makefile.in | 4 +- Cantera/src/numerics/ResidEval.h | 9 +- Cantera/src/numerics/ResidJacEval.cpp | 7 + Cantera/src/numerics/ResidJacEval.h | 14 + 10 files changed, 1009 insertions(+), 492 deletions(-) diff --git a/Cantera/src/numerics/DAE_Solver.h b/Cantera/src/numerics/DAE_Solver.h index f57d28bbe..cbc4d4b03 100644 --- a/Cantera/src/numerics/DAE_Solver.h +++ b/Cantera/src/numerics/DAE_Solver.h @@ -12,207 +12,236 @@ * Copyright 2006 California Institute of Technology * */ - -#undef DAE_DEVEL - #ifndef CT_DAE_Solver_H #define CT_DAE_Solver_H #include #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 diff --git a/Cantera/src/numerics/DAE_solvers.cpp b/Cantera/src/numerics/DAE_solvers.cpp index 4e0c2d7c7..db3b8b339 100644 --- a/Cantera/src/numerics/DAE_solvers.cpp +++ b/Cantera/src/numerics/DAE_solvers.cpp @@ -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); } } } diff --git a/Cantera/src/numerics/FuncEval.h b/Cantera/src/numerics/FuncEval.h index d88181533..e835eb9ab 100644 --- a/Cantera/src/numerics/FuncEval.h +++ b/Cantera/src/numerics/FuncEval.h @@ -55,7 +55,7 @@ namespace Cantera { */ virtual int neq()=0; - /// Number of parameters. + //! Number of parameters. virtual int nparams() { return 0; } protected: diff --git a/Cantera/src/numerics/IDA_Solver.cpp b/Cantera/src/numerics/IDA_Solver.cpp index eec7f9226..6f2575138 100644 --- a/Cantera/src/numerics/IDA_Solver.cpp +++ b/Cantera/src/numerics/IDA_Solver.cpp @@ -1,4 +1,3 @@ - /** * @file IDA_Solver.cpp * @@ -10,8 +9,8 @@ #include "stringUtils.h" #include -using namespace std; +#ifdef SUNDIALS_VERSION_22 #include #include #include @@ -19,238 +18,506 @@ using namespace std; #include #include #include +#else +#include +#include +#include +#include +#include +#include +#include +#endif + +using namespace std; inline static N_Vector nv(void* x) { - return reinterpret_cast(x); + return reinterpret_cast(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(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(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(N_VNew_Serial(m_neq)); + m_ydot = reinterpret_cast(N_VNew_Serial(m_neq)); + m_constraints = reinterpret_cast(N_VNew_Serial(m_neq)); + + for (int i=0; i(N_VNew_Serial(m_neq)); - m_ydot = reinterpret_cast(N_VNew_Serial(m_neq)); - m_constraints = reinterpret_cast(N_VNew_Serial(m_neq)); - - for (int i=0; i 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; + } + //==================================================================================================================== + } - - diff --git a/Cantera/src/numerics/IDA_Solver.h b/Cantera/src/numerics/IDA_Solver.h index 59546d991..03d60c904 100644 --- a/Cantera/src/numerics/IDA_Solver.h +++ b/Cantera/src/numerics/IDA_Solver.h @@ -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 #include "DAE_Solver.h" #include "ctexceptions.h" + +#ifdef SUNDIALS_VERSION_22 +#include +#else +#include + +// 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; + }; } diff --git a/Cantera/src/numerics/Integrator.h b/Cantera/src/numerics/Integrator.h index 55e8185c3..be04bc2d4 100644 --- a/Cantera/src/numerics/Integrator.h +++ b/Cantera/src/numerics/Integrator.h @@ -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"); } diff --git a/Cantera/src/numerics/Makefile.in b/Cantera/src/numerics/Makefile.in index bfd885d08..8787bb236 100644 --- a/Cantera/src/numerics/Makefile.in +++ b/Cantera/src/numerics/Makefile.in @@ -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 diff --git a/Cantera/src/numerics/ResidEval.h b/Cantera/src/numerics/ResidEval.h index 3e25f8129..4f1eb5e61 100644 --- a/Cantera/src/numerics/ResidEval.h +++ b/Cantera/src/numerics/ResidEval.h @@ -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: diff --git a/Cantera/src/numerics/ResidJacEval.cpp b/Cantera/src/numerics/ResidJacEval.cpp index cb6d08988..ea5ad953d 100644 --- a/Cantera/src/numerics/ResidJacEval.cpp +++ b/Cantera/src/numerics/ResidJacEval.cpp @@ -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 diff --git a/Cantera/src/numerics/ResidJacEval.h b/Cantera/src/numerics/ResidJacEval.h index a1799253e..5277388fb 100644 --- a/Cantera/src/numerics/ResidJacEval.h +++ b/Cantera/src/numerics/ResidJacEval.h @@ -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.