From a81572e92096b8b12d180df9ac009170c6331945 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 24 Oct 2005 16:41:24 +0000 Subject: [PATCH] initial import --- Cantera/src/CVodesIntegrator.cpp | 319 +++++++++++++++++++++++++++++++ Cantera/src/CVodesIntegrator.h | 89 +++++++++ Cantera/src/ODE_integrators.cpp | 25 +++ 3 files changed, 433 insertions(+) create mode 100644 Cantera/src/CVodesIntegrator.cpp create mode 100644 Cantera/src/CVodesIntegrator.h create mode 100644 Cantera/src/ODE_integrators.cpp diff --git a/Cantera/src/CVodesIntegrator.cpp b/Cantera/src/CVodesIntegrator.cpp new file mode 100644 index 000000000..56f66e482 --- /dev/null +++ b/Cantera/src/CVodesIntegrator.cpp @@ -0,0 +1,319 @@ +/** + * @file CVodeInt.cpp + * + */ + +// Copyright 2001 California Institute of Technology + +#include "CVodesIntegrator.h" +#include +using namespace std; + + +// sundials includes +#include +#include +#include +#include +#include +#include +#include +#include + +inline static N_Vector nv(void* x) { + return reinterpret_cast(x); +} + +extern "C" { + + /** + * Function called by cvodes to evaluate ydot given y. The cvode + * integrator allows passing in a void* pointer to access + * external data. This pointer is cast to a pointer to a instance + * of class FuncEval. The equations to be integrated should be + * specified by deriving a class from FuncEval that evaluates the + * desired equations. + * @ingroup odeGroup + */ + static void cvodes_rhs(realtype t, N_Vector y, N_Vector ydot, + void *f_data) { + double* ydata = NV_DATA_S(y); //N_VDATA(y); + double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot); + Cantera::FuncEval* f = (Cantera::FuncEval*)f_data; + f->eval(t, ydata, ydotdata, NULL); + } + +} + +namespace Cantera { + + + /** + * Constructor. Default settings: dense jacobian, no user-supplied + * Jacobian function, Newton iteration. + */ + CVodesIntegrator::CVodesIntegrator() : m_neq(0), + m_cvode_mem(0), + m_t0(0.0), + m_y(0), + m_abstol(0), + m_type(DENSE+NOJAC), + m_itol(CV_SS), + m_method(CV_BDF), + m_iter(CV_NEWTON), + m_maxord(0), + m_reltol(1.e-9), + m_abstols(1.e-15), + m_nabs(0), + m_hmax(0.0), + m_maxsteps(20000) + { + //m_ropt.resize(OPT_SIZE,0.0); + //m_iopt = new long[OPT_SIZE]; + //fill(m_iopt, m_iopt+OPT_SIZE,0); + } + + + /// Destructor. + CVodesIntegrator::~CVodesIntegrator() + { + if (m_cvode_mem) CVodeFree(m_cvode_mem); + if (m_y) N_VDestroy_Serial(nv(m_y)); //N_VFree(nv(m_y)); + if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); //N_VFree(nv(m_abstol)); + //delete[] m_iopt; + } + + double& CVodesIntegrator::solution(int k){ + return NV_Ith_S(nv(m_y),k); + } + + double* CVodesIntegrator::solution(){ return NV_DATA_S(nv(m_y)); + } + + void CVodesIntegrator::setTolerances(double reltol, int n, double* abstol) { + m_itol = CV_SV; + m_nabs = n; + if (n != m_neq) { + if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); + m_abstol = reinterpret_cast(N_VNew_Serial(n)); + } + for (int i=0; i(N_VNew_Serial(m_neq)); // allocate solution vector + 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 CVodesIntegrator::reinitialize(double t0, FuncEval& func) + { + m_t0 = t0; + func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y))); + + // set options + // m_iopt[MXSTEP] = m_maxsteps; + //m_iopt[MAXORD] = m_maxord; + //m_ropt[HMAX] = m_hmax; + + //if (m_cvode_mem) CVodeFree(m_cvode_mem); + + int result; + if (m_itol == CV_SV) { + result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), + m_itol, m_reltol, + nv(m_abstol)); + } + else { + result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), + m_itol, m_reltol, + &m_abstols); + } + + if (result != 0) throw CVodesErr("CVReInit failed."); + + if (m_type == DENSE + NOJAC) { + long int N = m_neq; + CVDense(m_cvode_mem, N); + } + else if (m_type == DIAG) { + CVDiag(m_cvode_mem); + } + else if (m_type == GMRES) { + CVSpgmr(m_cvode_mem, PREC_NONE, 0); + } + else { + throw CVodesErr("unsupported option"); + } + + // pass a pointer to func in m_data + m_data = (void*)&func; + long int flag = CVodeSetFdata(m_cvode_mem, m_data); + if (flag != CV_SUCCESS) + throw CVodesErr("CVodeSetFdata 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 CVodesIntegrator::integrate(double tout) + { + double t; + int flag; + flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL); + if (flag != CV_SUCCESS) + throw CVodesErr(" CVodes error encountered."); + } + + double CVodesIntegrator::step(double tout) + { + double t; + int flag; + flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_ONE_STEP); + if (flag != CV_SUCCESS) + throw CVodesErr(" CVodes error encountered."); + return t; + } + + int CVodesIntegrator::nEvals() const { + long int ne; + return CVodeGetNumRhsEvals(m_cvode_mem, &ne); + return ne; + //return m_iopt[NFE]; + } +} + + diff --git a/Cantera/src/CVodesIntegrator.h b/Cantera/src/CVodesIntegrator.h new file mode 100644 index 000000000..cfb45d145 --- /dev/null +++ b/Cantera/src/CVodesIntegrator.h @@ -0,0 +1,89 @@ +/** + * @file CVodesWrapper.h + */ + +/* $Author$ + * $Date$ + * $Revision$ + */ + +// Copyright 2005 California Institute of Technology + + +#ifndef CT_CVODESWRAPPER_H +#define CT_CVODESWRAPPER_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "Integrator.h" +#include "FuncEval.h" +#include "ctexceptions.h" +#include "ct_defs.h" + +namespace Cantera { + + /** + * Exception thrown when a CVODES error is encountered. + */ + class CVodesErr : public CanteraError { + public: + CVodesErr(string msg) : CanteraError("CVodesIntegrator", msg){} + }; + + + /** + * Wrapper class for 'cvodes' integrator from LLNL. + * + * @see FuncEval.h. Classes that use CVodeInt: + * ImplicitChem, ImplicitSurfChem, Reactor + * + */ + class CVodesIntegrator : public Integrator { + + public: + + CVodesIntegrator(); + virtual ~CVodesIntegrator(); + virtual void setTolerances(double reltol, int n, double* abstol); + virtual void setTolerances(double reltol, double abstol); + virtual void setProblemType(int probtype); + virtual void initialize(double t0, FuncEval& func); + virtual void reinitialize(double t0, FuncEval& func); + virtual void integrate(double tout); + virtual doublereal step(double tout); + virtual double& solution(int k); + virtual double* solution(); + virtual int nEquations() const { return m_neq;} + virtual int nEvals() const; + virtual void setMaxOrder(int n) { m_maxord = n; } + virtual void setMethod(MethodType t); + virtual void setIterator(IterType t); + virtual void setMaxStepSize(double hmax); + virtual void setMinStepSize(double hmin); + virtual void setMaxSteps(int nmax); + + private: + + int m_neq; + void* m_cvode_mem; + double m_t0; + void *m_y, *m_abstol; + int m_type; + int m_itol; + int m_method; + int m_iter; + int m_maxord; + double m_reltol; + double m_abstols; + int m_nabs; + double m_hmax, m_hmin; + int m_maxsteps; + void* m_data; + }; + +} // namespace + +#endif diff --git a/Cantera/src/ODE_integrators.cpp b/Cantera/src/ODE_integrators.cpp new file mode 100644 index 000000000..b3d34d5b3 --- /dev/null +++ b/Cantera/src/ODE_integrators.cpp @@ -0,0 +1,25 @@ +#include "ct_defs.h" +#include "Integrator.h" + +#ifdef HAS_SUNDIALS +#include "CVodesIntegrator.cpp" +#else +#include "CVode.cpp" +#endif + +// namespace Cantera { + +// Integrator* newIntegrator(string itype) { +// if (itype == "CVODE") { +// #ifdef HAS_SUNDIALS +// return new CVodesIntegrator(); +// #else +// return new CVodeInt(); +// #endif +// } +// else { +// throw CanteraError("newIntegrator", +// "unknown ODE integrator: "+itype); +// } +// } +// }