diff --git a/Cantera/src/numerics/CMakeLists.txt b/Cantera/src/numerics/CMakeLists.txt index 1da3015c9..d7c452775 100644 --- a/Cantera/src/numerics/CMakeLists.txt +++ b/Cantera/src/numerics/CMakeLists.txt @@ -9,7 +9,7 @@ ADD_LIBRARY(numerics ${NUMERICS_SRCS}) TARGET_LINK_LIBRARIES (numerics ctbase cvode) -SET(NUMERICS_H ArrayViewer.h CVode.h CVodesIntegrator.h DenseMatrix.h +SET(NUMERICS_H ArrayViewer.h CVodeInt.h CVodesIntegrator.h DenseMatrix.h funcs.h ctlapack.h Func1.h FuncEval.h polyfit.h BandMatrix.h Integrator.h DAE_Solver.h ResidEval.h sort.h) diff --git a/Cantera/src/numerics/CVode.cpp b/Cantera/src/numerics/CVode.cpp deleted file mode 100755 index f27b7e11c..000000000 --- a/Cantera/src/numerics/CVode.cpp +++ /dev/null @@ -1,310 +0,0 @@ -/** - * @file CVode.cpp - * - * $Id$ - */ - -// Copyright 2001 California Institute of Technology - -#include "CVode.h" -#include -using namespace std; - - -// cvode includes -#include "../../../ext/cvode/include/llnltyps.h" -#include "../../../ext/cvode/include/llnlmath.h" -#include "../../../ext/cvode/include/cvode.h" -#include "../../../ext/cvode/include/cvdense.h" -#include "../../../ext/cvode/include/cvdiag.h" -#include "../../../ext/cvode/include/cvspgmr.h" -#include "../../../ext/cvode/include/nvector.h" -#include "../../../ext/cvode/include/cvode.h" - -inline static N_Vector nv(void* x) { - return reinterpret_cast(x); -} - -extern "C" { - - /** - * Function called by cvode 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 cvode_rhs(integer N, real t, N_Vector y, N_Vector ydot, - void *f_data) { - double* ydata = N_VDATA(y); - double* ydotdata = N_VDATA(ydot); - Cantera::FuncEval* f = (Cantera::FuncEval*)f_data; - f->eval(t, ydata, ydotdata, NULL); - } - - - /** - * Function called by cvode to evaluate the Jacobian matrix. - * (temporary) - * @ingroup odeGroup - */ - static void cvode_jac(integer N, DenseMat J, RhsFn f, void *f_data, - real t, N_Vector y, N_Vector fy, N_Vector ewt, real h, real uround, - void *jac_data, long int *nfePtr, N_Vector vtemp1, N_Vector vtemp2, - N_Vector vtemp3) - { - // get pointers to start of data - double* ydata = N_VDATA(y); - double* fydata = N_VDATA(fy); - double* ewtdata = N_VDATA(ewt); - double* ydot = N_VDATA(vtemp1); - - Cantera::FuncEval* func = (Cantera::FuncEval*)f_data; - - int i,j; - double* col_j; - double ysave, dy; - for (j=0; j < N; j++) { - col_j = (J->data)[j]; - ysave = ydata[j]; - dy = 1.0/ewtdata[j]; - ydata[j] = ysave + dy; - dy = ydata[j] - ysave; - func->eval(t, ydata, ydot, NULL); - for (i=0; i < N; i++) { - col_j[i] = (ydot[i] - fydata[i])/dy; - } - ydata[j] = ysave; - } - } -} - -namespace Cantera { - - - /** - * Constructor. Default settings: dense jacobian, no user-supplied - * Jacobian function, Newton iteration. - */ - CVodeInt::CVodeInt() : m_neq(0), - m_cvode_mem(0), - m_t0(0.0), - m_y(0), - m_abstol(0), - m_type(DENSE+NOJAC), - m_itol(0), - m_method(BDF), - m_iter(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. - CVodeInt::~CVodeInt() - { - if (m_cvode_mem) CVodeFree(m_cvode_mem); - if (m_y) N_VFree(nv(m_y)); - if (m_abstol) N_VFree(nv(m_abstol)); - delete[] m_iopt; - } - - double& CVodeInt::solution(int k){ return N_VIth(nv(m_y),k); } - double* CVodeInt::solution(){ return N_VDATA(nv(m_y)); } - - void CVodeInt::setTolerances(double reltol, int n, double* abstol) { - m_itol = 1; - m_nabs = n; - if (n != m_neq) { - if (m_abstol) N_VFree(nv(m_abstol)); - m_abstol = reinterpret_cast(N_VNew(n, 0)); - } - for (int i=0; i(N_VNew(m_neq, 0)); // allocate solution vector - // check abs tolerance array size - if (m_itol == 1 && m_nabs < m_neq) - throw CVodeErr("not enough absolute tolerance values specified."); - func.getInitialConditions(m_t0, m_neq, N_VDATA(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); - - // pass a pointer to func in m_data - m_data = (void*)&func; - - if (m_itol) { - m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, nv(m_y), m_method, - m_iter, m_itol, &m_reltol, - nv(m_abstol), m_data, NULL, TRUE, m_iopt, - DATA_PTR(m_ropt), NULL); - } - else { - m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, nv(m_y), m_method, - m_iter, m_itol, &m_reltol, - &m_abstols, m_data, NULL, TRUE, m_iopt, - DATA_PTR(m_ropt), NULL); - } - - if (!m_cvode_mem) throw CVodeErr("CVodeMalloc failed."); - - if (m_type == DENSE + NOJAC) { - CVDense(m_cvode_mem, NULL, NULL); - } - else if (m_type == DENSE + JAC) { - CVDense(m_cvode_mem, cvode_jac, NULL); - } - else if (m_type == DIAG) { - CVDiag(m_cvode_mem); - } - else if (m_type == GMRES) { - CVSpgmr(m_cvode_mem, NONE, MODIFIED_GS, 0, 0.0, - NULL, NULL, NULL); - } - else { - throw CVodeErr("unsupported option"); - } - } - - - void CVodeInt::reinitialize(double t0, FuncEval& func) - { - m_t0 = t0; - func.getInitialConditions(m_t0, m_neq, N_VDATA(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); - - // pass a pointer to func in m_data - m_data = (void*)&func; - int result; - if (m_itol) { - result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, nv(m_y), m_method, - m_iter, m_itol, &m_reltol, - nv(m_abstol), m_data, NULL, TRUE, m_iopt, - DATA_PTR(m_ropt), NULL); - } - else { - result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, nv(m_y), m_method, - m_iter, m_itol, &m_reltol, - &m_abstols, m_data, NULL, TRUE, m_iopt, - DATA_PTR(m_ropt), NULL); - } - - if (result != 0) throw CVodeErr("CVReInit failed."); - - if (m_type == DENSE + NOJAC) { - CVDense(m_cvode_mem, NULL, NULL); - } - else if (m_type == DENSE + JAC) { - CVDense(m_cvode_mem, cvode_jac, NULL); - } - else if (m_type == DIAG) { - CVDiag(m_cvode_mem); - } - else if (m_type == GMRES) { - CVSpgmr(m_cvode_mem, NONE, MODIFIED_GS, 0, 0.0, - NULL, NULL, NULL); - } - else { - throw CVodeErr("unsupported option"); - } - } - - void CVodeInt::integrate(double tout) - { - double t; - int flag; - flag = CVode(m_cvode_mem, tout, nv(m_y), &t, NORMAL); - if (flag != SUCCESS) - throw CVodeErr(" CVode error encountered."); - } - - double CVodeInt::step(double tout) - { - double t; - int flag; - flag = CVode(m_cvode_mem, tout, nv(m_y), &t, ONE_STEP); - if (flag != SUCCESS) - throw CVodeErr(" CVode error encountered."); - return t; - } - - int CVodeInt::nEvals() const { return m_iopt[NFE]; } -} - - - diff --git a/Cantera/src/numerics/CVode.h b/Cantera/src/numerics/CVode.h deleted file mode 100755 index 8e7e62723..000000000 --- a/Cantera/src/numerics/CVode.h +++ /dev/null @@ -1,93 +0,0 @@ -/** - * @file CVode.h - */ - -/* $Author$ - * $Date$ - * $Revision$ - */ - -// Copyright 2001 California Institute of Technology - - -#ifndef CT_CVODE_H -#define CT_CVODE_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 CVODE error is encountered. - */ - class CVodeErr : public CanteraError { - public: - CVodeErr(std::string msg) : CanteraError("CVodeInt", msg){} - }; - - - /** - * Wrapper class for 'cvode' integrator from LLNL. - * The unmodified cvode code is in directory ext/cvode. - * - * @see FuncEval.h. Classes that use CVodeInt: - * ImplicitChem, ImplicitSurfChem, Reactor - * - */ - class CVodeInt : public Integrator { - - public: - - CVodeInt(); - virtual ~CVodeInt(); - 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; - - vector_fp m_ropt; - long int* m_iopt; - void* m_data; - }; - -} // namespace - -#endif // CT_CVODE diff --git a/Cantera/src/numerics/CVodeInt.cpp b/Cantera/src/numerics/CVodeInt.cpp new file mode 100644 index 000000000..658094a5e --- /dev/null +++ b/Cantera/src/numerics/CVodeInt.cpp @@ -0,0 +1,310 @@ +/** + * @file CVodeInt.cpp + * + * $Id$ + */ + +// Copyright 2001 California Institute of Technology + +#include "CVodeInt.h" +#include +using namespace std; + + +// cvode includes +#include "../../../ext/cvode/include/llnltyps.h" +#include "../../../ext/cvode/include/llnlmath.h" +#include "../../../ext/cvode/include/cvode.h" +#include "../../../ext/cvode/include/cvdense.h" +#include "../../../ext/cvode/include/cvdiag.h" +#include "../../../ext/cvode/include/cvspgmr.h" +#include "../../../ext/cvode/include/nvector.h" +#include "../../../ext/cvode/include/cvode.h" + +inline static N_Vector nv(void* x) { + return reinterpret_cast(x); +} + +extern "C" { + + /** + * Function called by cvode 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 cvode_rhs(integer N, real t, N_Vector y, N_Vector ydot, + void *f_data) { + double* ydata = N_VDATA(y); + double* ydotdata = N_VDATA(ydot); + Cantera::FuncEval* f = (Cantera::FuncEval*)f_data; + f->eval(t, ydata, ydotdata, NULL); + } + + + /** + * Function called by cvode to evaluate the Jacobian matrix. + * (temporary) + * @ingroup odeGroup + */ + static void cvode_jac(integer N, DenseMat J, RhsFn f, void *f_data, + real t, N_Vector y, N_Vector fy, N_Vector ewt, real h, real uround, + void *jac_data, long int *nfePtr, N_Vector vtemp1, N_Vector vtemp2, + N_Vector vtemp3) + { + // get pointers to start of data + double* ydata = N_VDATA(y); + double* fydata = N_VDATA(fy); + double* ewtdata = N_VDATA(ewt); + double* ydot = N_VDATA(vtemp1); + + Cantera::FuncEval* func = (Cantera::FuncEval*)f_data; + + int i,j; + double* col_j; + double ysave, dy; + for (j=0; j < N; j++) { + col_j = (J->data)[j]; + ysave = ydata[j]; + dy = 1.0/ewtdata[j]; + ydata[j] = ysave + dy; + dy = ydata[j] - ysave; + func->eval(t, ydata, ydot, NULL); + for (i=0; i < N; i++) { + col_j[i] = (ydot[i] - fydata[i])/dy; + } + ydata[j] = ysave; + } + } +} + +namespace Cantera { + + + /** + * Constructor. Default settings: dense jacobian, no user-supplied + * Jacobian function, Newton iteration. + */ + CVodeInt::CVodeInt() : m_neq(0), + m_cvode_mem(0), + m_t0(0.0), + m_y(0), + m_abstol(0), + m_type(DENSE+NOJAC), + m_itol(0), + m_method(BDF), + m_iter(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. + CVodeInt::~CVodeInt() + { + if (m_cvode_mem) CVodeFree(m_cvode_mem); + if (m_y) N_VFree(nv(m_y)); + if (m_abstol) N_VFree(nv(m_abstol)); + delete[] m_iopt; + } + + double& CVodeInt::solution(int k){ return N_VIth(nv(m_y),k); } + double* CVodeInt::solution(){ return N_VDATA(nv(m_y)); } + + void CVodeInt::setTolerances(double reltol, int n, double* abstol) { + m_itol = 1; + m_nabs = n; + if (n != m_neq) { + if (m_abstol) N_VFree(nv(m_abstol)); + m_abstol = reinterpret_cast(N_VNew(n, 0)); + } + for (int i=0; i(N_VNew(m_neq, 0)); // allocate solution vector + // check abs tolerance array size + if (m_itol == 1 && m_nabs < m_neq) + throw CVodeErr("not enough absolute tolerance values specified."); + func.getInitialConditions(m_t0, m_neq, N_VDATA(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); + + // pass a pointer to func in m_data + m_data = (void*)&func; + + if (m_itol) { + m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, nv(m_y), m_method, + m_iter, m_itol, &m_reltol, + nv(m_abstol), m_data, NULL, TRUE, m_iopt, + DATA_PTR(m_ropt), NULL); + } + else { + m_cvode_mem = CVodeMalloc(m_neq, cvode_rhs, m_t0, nv(m_y), m_method, + m_iter, m_itol, &m_reltol, + &m_abstols, m_data, NULL, TRUE, m_iopt, + DATA_PTR(m_ropt), NULL); + } + + if (!m_cvode_mem) throw CVodeErr("CVodeMalloc failed."); + + if (m_type == DENSE + NOJAC) { + CVDense(m_cvode_mem, NULL, NULL); + } + else if (m_type == DENSE + JAC) { + CVDense(m_cvode_mem, cvode_jac, NULL); + } + else if (m_type == DIAG) { + CVDiag(m_cvode_mem); + } + else if (m_type == GMRES) { + CVSpgmr(m_cvode_mem, NONE, MODIFIED_GS, 0, 0.0, + NULL, NULL, NULL); + } + else { + throw CVodeErr("unsupported option"); + } + } + + + void CVodeInt::reinitialize(double t0, FuncEval& func) + { + m_t0 = t0; + func.getInitialConditions(m_t0, m_neq, N_VDATA(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); + + // pass a pointer to func in m_data + m_data = (void*)&func; + int result; + if (m_itol) { + result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, nv(m_y), m_method, + m_iter, m_itol, &m_reltol, + nv(m_abstol), m_data, NULL, TRUE, m_iopt, + DATA_PTR(m_ropt), NULL); + } + else { + result = CVReInit(m_cvode_mem, cvode_rhs, m_t0, nv(m_y), m_method, + m_iter, m_itol, &m_reltol, + &m_abstols, m_data, NULL, TRUE, m_iopt, + DATA_PTR(m_ropt), NULL); + } + + if (result != 0) throw CVodeErr("CVReInit failed."); + + if (m_type == DENSE + NOJAC) { + CVDense(m_cvode_mem, NULL, NULL); + } + else if (m_type == DENSE + JAC) { + CVDense(m_cvode_mem, cvode_jac, NULL); + } + else if (m_type == DIAG) { + CVDiag(m_cvode_mem); + } + else if (m_type == GMRES) { + CVSpgmr(m_cvode_mem, NONE, MODIFIED_GS, 0, 0.0, + NULL, NULL, NULL); + } + else { + throw CVodeErr("unsupported option"); + } + } + + void CVodeInt::integrate(double tout) + { + double t; + int flag; + flag = CVode(m_cvode_mem, tout, nv(m_y), &t, NORMAL); + if (flag != SUCCESS) + throw CVodeErr(" CVode error encountered."); + } + + double CVodeInt::step(double tout) + { + double t; + int flag; + flag = CVode(m_cvode_mem, tout, nv(m_y), &t, ONE_STEP); + if (flag != SUCCESS) + throw CVodeErr(" CVode error encountered."); + return t; + } + + int CVodeInt::nEvals() const { return m_iopt[NFE]; } +} + + + diff --git a/Cantera/src/numerics/CVodeInt.h b/Cantera/src/numerics/CVodeInt.h new file mode 100644 index 000000000..ac1836b1a --- /dev/null +++ b/Cantera/src/numerics/CVodeInt.h @@ -0,0 +1,93 @@ +/** + * @file CVodeInt.h + */ + +/* $Author$ + * $Date$ + * $Revision$ + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_CVODEINT_H +#define CT_CVODEINT_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 CVODE error is encountered. + */ + class CVodeErr : public CanteraError { + public: + CVodeErr(std::string msg) : CanteraError("CVodeInt", msg){} + }; + + + /** + * Wrapper class for 'cvode' integrator from LLNL. + * The unmodified cvode code is in directory ext/cvode. + * + * @see FuncEval.h. Classes that use CVodeInt: + * ImplicitChem, ImplicitSurfChem, Reactor + * + */ + class CVodeInt : public Integrator { + + public: + + CVodeInt(); + virtual ~CVodeInt(); + 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; + + vector_fp m_ropt; + long int* m_iopt; + void* m_data; + }; + +} // namespace + +#endif // CT_CVODE diff --git a/Cantera/src/numerics/Makefile.in b/Cantera/src/numerics/Makefile.in index 690fb0fa4..bffd2a159 100644 --- a/Cantera/src/numerics/Makefile.in +++ b/Cantera/src/numerics/Makefile.in @@ -46,8 +46,8 @@ ifeq ($(use_sundials), 1) ODEPACKAGE_H = CVodesIntegrator.h ODEPACKAGE_OBJ = CVodesIntegrator.o else - ODEPACKAGE_H = CVode.h - ODEPACKAGE_OBJ = CVode.o + ODEPACKAGE_H = CVodeInt.h + ODEPACKAGE_OBJ = CVodeInt.o endif TOTAL_H = $(NUMERICS_H) $(ODEPACKAGE_H) diff --git a/Cantera/src/numerics/ODE_integrators.cpp b/Cantera/src/numerics/ODE_integrators.cpp index 2846b5ace..e0c613b09 100644 --- a/Cantera/src/numerics/ODE_integrators.cpp +++ b/Cantera/src/numerics/ODE_integrators.cpp @@ -5,7 +5,7 @@ #ifdef HAS_SUNDIALS #include "CVodesIntegrator.h" #else -#include "CVode.h" +#include "CVodeInt.h" #endif namespace Cantera {