diff --git a/Cantera/src/numerics/ArrayViewer.h b/Cantera/src/numerics/ArrayViewer.h new file mode 100755 index 000000000..112648530 --- /dev/null +++ b/Cantera/src/numerics/ArrayViewer.h @@ -0,0 +1,153 @@ +/** + * @file ArrayViewer.h + * + * Header file for class ArrayViewer + */ + +/* $Author$ + * $Revision$ + * $Date$ + * + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_ARRAYVIEWER_H +#define CT_ARRAYVIEWER_H + +#include "ct_defs.h" +#include "ctexceptions.h" +#include "stringUtils.h" +#include "utilities.h" + +namespace Cantera { + + + /** + * An interface for 2D arrays stored in column-major + * (Fortran-compatible) form. This class is designed for + * situations when you have a Fortran-compatible 2D array that + * you want to be able to easily get/set individual elements or + * entire rows or columns. Instances of ArrayViewer store only a + * pointer to the first element in the array; no copy is made of + * the array. @todo This class should probably be renamed, since + * it not only views, but also can modify, array elements. This + * class is hardly used; candidate for removal. + */ + class ArrayViewer { + + public: + + typedef doublereal* iterator; + typedef const doublereal* const_iterator; + + + /** + * Default constructor. Create an empty array viewer. + */ + ArrayViewer() : m_nrows(0), m_ncols(0) { data = 0; } + + + /** + * Constructor. Create an \c m by \c n array viewer for array v. + */ + ArrayViewer(int m, int n, doublereal* v) + : m_nrows(m), m_ncols(n) { + data = v; + } + + /// resize the array viewer + void resize(int n, int m) { + m_nrows = n; + m_ncols = m; + } + + /// set the nth row to array rw + void setRow(int n, doublereal* rw) { + for (int j = 0; j < m_ncols; j++) { + data[m_nrows*j + n] = rw[j]; + } + } + + /// get the nth row + void getRow(int n, doublereal* rw) { + for (int j = 0; j < m_ncols; j++) { + rw[j] = data[m_nrows*j + n]; + } + } + + /// set the values in column m to those in array col + void setColumn(int m, doublereal* col) { + for (int i = 0; i < m_nrows; i++) { + data[m_nrows*m + i] = col[i]; + } + } + + /// get the values in column m + void getColumn(int m, doublereal* col) { + for (int i = 0; i < m_nrows; i++) { + col[i] = data[m_nrows*m + i]; + } + } + + /// Destructor. Does nothing. + virtual ~ArrayViewer(){} + + doublereal& operator()( int i, int j) {return value(i,j);} + doublereal operator() ( int i, int j) const {return value(i,j);} + + /// Return a reference to the (i,j) array element. + doublereal& value( int i, int j) {return data[m_nrows*j + i];} + + /// Return the value of the (i,j) array element. + doublereal value( int i, int j) const {return data[m_nrows*j + i];} + + /// Number of rows + size_t nRows() const { return m_nrows; } + + /// Number of columns + size_t nColumns() const { return m_ncols; } + + iterator begin() { return data; } + iterator end() { return data + m_nrows*m_ncols; } + const_iterator begin() const { return data; } + const_iterator end() const { return data + m_nrows*m_ncols; } + + doublereal* data; + + protected: + + int m_nrows, m_ncols; + }; + + /// output the array + inline std::ostream& operator<<(std::ostream& s, const ArrayViewer& m) { + int nr = static_cast(m.nRows()); + int nc = static_cast(m.nColumns()); + int i,j; + for (i = 0; i < nr; i++) { + for (j = 0; j < nc; j++) { + s << m(i,j) << ", "; + } + s << std::endl; + } + return s; + } + + /// Multiply the array by a constant. + inline void operator*=(ArrayViewer& m, doublereal a) { + scale(m.begin(), m.end(), m.begin(), a); + } + + /// Increment the entire array by a constant. + inline void operator+=(ArrayViewer& x, const ArrayViewer& y) { + sum_each(x.begin(), x.end(), y.begin()); + } + +} + +#endif + + + diff --git a/Cantera/src/numerics/BandMatrix.cpp b/Cantera/src/numerics/BandMatrix.cpp new file mode 100755 index 000000000..b643fe1b2 --- /dev/null +++ b/Cantera/src/numerics/BandMatrix.cpp @@ -0,0 +1,181 @@ +/** + * @file BandMatrix.cpp + * + * Banded matrices. + */ + +// Copyright 2001 California Institute of Technology + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "BandMatrix.h" +#include "ctlapack.h" +#include "utilities.h" +#include "ctexceptions.h" +#include "stringUtils.h" +#include "global.h" + +using namespace std; + +namespace Cantera { + + + /// Default constructor. + BandMatrix::BandMatrix() : m_factored(false), m_n(0), + m_kl(0), m_ku(0), m_zero(0.0) { + data.clear(); ludata.clear(); + } + + + /** + * Constructor. Create an n by n banded matrix. + * @param n number of rows and columns + * @param kl number of subdiagonals + * @param ku number of superdiagonals + * @param v initial value (default = 0.0) + */ + BandMatrix::BandMatrix(int n, int kl, int ku, doublereal v) + : m_factored(false), m_n(n), m_kl(kl), m_ku(ku) { + data.resize(n*(2*kl + ku + 1)); + ludata.resize(n*(2*kl + ku + 1)); + fill(data.begin(), data.end(), v); + fill(ludata.begin(), ludata.end(), 0.0); + m_ipiv.resize(m_n); + } + + /// copy constructor + BandMatrix::BandMatrix(const BandMatrix& y) { + m_n = y.m_n; + m_kl = y.m_kl; + m_ku = y.m_ku; + data = y.data; + ludata = y.ludata; + m_factored = y.m_factored; + m_ipiv = y.m_ipiv; + } + + BandMatrix& BandMatrix::operator=(const BandMatrix& y) { + if (&y == this) return *this; + m_n = y.m_n; + m_kl = y.m_kl; + m_ku = y.m_ku; + m_ipiv = y.m_ipiv; + data = y.data; + ludata = y.ludata; + m_factored = y.m_factored; + return *this; + } + + void BandMatrix::resize(int n, int kl, int ku, doublereal v) { + m_n = n; + m_kl = kl; + m_ku = ku; + data.resize(n*(2*kl + ku + 1)); + ludata.resize(n*(2*kl + ku + 1)); + m_ipiv.resize(m_n); + fill(data.begin(), data.end(), v); + fill(data.begin(), data.end(), 0.0); + m_factored = false; + } + + + /** + * Multiply A*b and write result to \c prod. + */ + void BandMatrix::mult(const double* b, double* prod) const { + int nr = rows(); + int m, j; + double sum = 0.0; + for (m = 0; m < nr; m++) { + sum = 0.0; + for (j = m - m_kl; j <= m + m_ku; j++) { + if (j >= 0 && j < m_n) + sum += _value(m,j)*b[j]; + } + prod[m] = sum; + } + } + + + /** + * Multiply b*A and write result to \c prod. + */ + void BandMatrix::leftMult(const double* b, double* prod) const { + int nc = columns(); + int n, i; + double sum = 0.0; + for (n = 0; n < nc; n++) { + sum = 0.0; + for (i = n - m_ku; i <= n + m_kl; i++) { + if (i >= 0 && i < m_n) + sum += _value(i,n)*b[i]; + } + prod[n] = sum; + } + } + + + /** + * Perform an LU decomposition. LAPACK routine DGBTRF is used. + * The factorization is saved in ludata. + */ + int BandMatrix::factor() { + int info=0; + copy(data.begin(), data.end(), ludata.begin()); + ct_dgbtrf(rows(), columns(), nSubDiagonals(), nSuperDiagonals(), + DATA_PTR(ludata), ldim(), DATA_PTR(ipiv()), info); + + // if info = 0, LU decomp succeeded. + if (info == 0) { + m_factored = true; + } + else { + m_factored = false; + ofstream fout("bandmatrix.csv"); + fout << *this << endl; + fout.close(); + } + return info; + } + + + + int BandMatrix::solve(int n, const doublereal* b, doublereal* x) { + copy(b, b+n, x); + return solve(n, x); + } + + int BandMatrix::solve(int n, doublereal* b) { + int info = 0; + if (!m_factored) info = factor(); + if (info == 0) + ct_dgbtrs(ctlapack::NoTranspose, columns(), nSubDiagonals(), + nSuperDiagonals(), 1, DATA_PTR(ludata), ldim(), + DATA_PTR(ipiv()), b, columns(), info); + + // error handling + if (info != 0) { + ofstream fout("bandmatrix.csv"); + fout << *this << endl; + fout.close(); + } + return info; + } + + ostream& operator<<(ostream& s, const BandMatrix& m) { + int nr = m.rows(); + int nc = m.columns(); + int i,j; + for (i = 0; i < nr; i++) { + for (j = 0; j < nc; j++) { + s << m(i,j) << ", "; + } + s << endl; + } + return s; + } +} + diff --git a/Cantera/src/numerics/BandMatrix.h b/Cantera/src/numerics/BandMatrix.h new file mode 100755 index 000000000..54d5211a6 --- /dev/null +++ b/Cantera/src/numerics/BandMatrix.h @@ -0,0 +1,153 @@ +/** + * @file BandMatrix.h + * + * Banded matrices. + */ + +/* + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_BANDMATRIX_H +#define CT_BANDMATRIX_H + +#include "ct_defs.h" +#include "ctlapack.h" +#include "utilities.h" +#include "ctexceptions.h" + +namespace Cantera { + + /** + * A class for banded matrices. + */ + class BandMatrix { + + public: + + BandMatrix(); + BandMatrix(int n, int kl, int ku, doublereal v = 0.0); + + /// copy constructor + BandMatrix(const BandMatrix& y); + + /// Destructor. Does nothing. + virtual ~BandMatrix(){} + + /// assignment. + BandMatrix& operator=(const BandMatrix& y); + + void resize(int n, int kl, int ku, doublereal v = 0.0); + + void bfill(doublereal v) { + std::fill(data.begin(), data.end(), v); + m_factored = false; + } + + doublereal& operator()( int i, int j) { + return value(i,j); + } + + doublereal operator() ( int i, int j) const { + return value(i,j); + } + + /// Return a reference to element (i,j). Since this method may + /// alter the element value, it may need to be refactored, so + /// the flag m_factored is set to false. + doublereal& value( int i, int j) { + m_factored = false; + if (i < j - m_ku || i > j + m_kl) { + m_zero = 0.0; + return m_zero; + } + return data[index(i,j)]; + } + + /// Return the value of element (i,j). This method does not + /// alter the array. + doublereal value( int i, int j) const { + if (i < j - m_ku || i > j + m_kl) return 0.0; + return data[index(i,j)]; + } + + /// Return the location in the internal 1D array corresponding to + /// the (i,j) element in the banded array. + int index(int i, int j) const { + int rw = m_kl + m_ku + i - j; + return (2*m_kl + m_ku + 1)*j + rw; + } + + /// Return the value of the (i,j) element for (i,j) within the + /// bandwidth. For efficiency, this method does not check that + /// (i,j) are within the bandwidth; it is up to the calling + /// program to insure that this is true. + doublereal _value(int i, int j) const { + return data[index(i,j)]; + } + + /// Number of rows + int nRows() const { return m_n; } + /// @deprecated Redundant. + int rows() const { return m_n; } + + /// Number of columns + int nColumns() const { return m_n; } + /// @deprecated Redundant. + int columns() const { return m_n; } + + /// Number of subdiagonals + int nSubDiagonals() const { return m_kl; } + + /// Number of superdiagonals + int nSuperDiagonals() const { return m_ku; } + + int ldim() const { return 2*m_kl + m_ku + 1; } + vector_int& ipiv() { return m_ipiv; } + + /// Multiply A*b and write result to prod. + void mult(const double* b, double* prod) const; + + /// Multiply b*A and write result to prod. + void leftMult(const double* b, double* prod) const; + + int factor(); + + //void solve(const vector_fp& b, vector_fp& x); + + int solve(int n, const doublereal* b, doublereal* x); + int solve(int n, doublereal* b); + + vector_fp::iterator begin() { + m_factored = false; + return data.begin(); + } + vector_fp::iterator end() { + m_factored = false; + return data.end(); + } + vector_fp::const_iterator begin() const { return data.begin(); } + vector_fp::const_iterator end() const { return data.end(); } + + vector_fp data; + vector_fp ludata; + bool m_factored; + + protected: + + int m_n, m_kl, m_ku; + doublereal m_zero; + vector_int m_ipiv; + + }; + + std::ostream& operator<<(std::ostream& s, const BandMatrix& m); + +} + +#endif diff --git a/Cantera/src/numerics/CVode.cpp b/Cantera/src/numerics/CVode.cpp new file mode 100755 index 000000000..13f7c1551 --- /dev/null +++ b/Cantera/src/numerics/CVode.cpp @@ -0,0 +1,310 @@ +/** + * @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 new file mode 100755 index 000000000..8e7e62723 --- /dev/null +++ b/Cantera/src/numerics/CVode.h @@ -0,0 +1,93 @@ +/** + * @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/CVodesIntegrator.cpp b/Cantera/src/numerics/CVodesIntegrator.cpp new file mode 100644 index 000000000..b08806560 --- /dev/null +++ b/Cantera/src/numerics/CVodesIntegrator.cpp @@ -0,0 +1,424 @@ +/** + * @file CVodesIntegrator.cpp + * + */ + +// Copyright 2001 California Institute of Technology +#include "config.h" + +#include "CVodesIntegrator.h" +#include "stringUtils.h" + +#include +using namespace std; + +#ifdef SUNDIALS_VERSION_22 + +#include +#include +#include +#include +#include +#include +#include +#include + +#else + +#ifdef SUNDIALS_VERSION_23 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#else + +unsupported sundials version! + +#endif + +#endif + +inline static N_Vector nv(void* x) { + return reinterpret_cast(x); +} + +namespace Cantera { + + class FuncData { + public: + FuncData(FuncEval* f, int npar = 0) { + m_pars.resize(npar, 1.0); + m_func = f; + } + virtual ~FuncData() {} + vector_fp m_pars; + FuncEval* m_func; + }; +} + + +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 int 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::FuncData* d = (Cantera::FuncData*)f_data; + Cantera::FuncEval* f = d->m_func; + //try { + if (d->m_pars.size() == 0) + f->eval(t, ydata, ydotdata, NULL); + else + f->eval(t, ydata, ydotdata, DATA_PTR(d->m_pars)); + //} + //catch (...) { + //Cantera::showErrors(); + //Cantera::error("Teminating execution"); + //} + return 0; + } + +} + +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_reltolsens(1.0e-5), + m_abstolsens(1.0e-4), + m_nabs(0), + m_hmax(0.0), + m_maxsteps(20000), m_np(0), + m_mupper(0), m_mlower(0) + { + //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) { + if (m_np > 0) + CVodeSensFree(m_cvode_mem); + CVodeFree(&m_cvode_mem); + } + if (m_y) N_VDestroy_Serial(nv(m_y)); + if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); + delete m_fdata; + + //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) { + sensInit(t0, func); + flag = CVodeSetSensParams(m_cvode_mem, DATA_PTR(m_fdata->m_pars), + NULL, NULL); + } + + // 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::reinitialize(double t0, FuncEval& func) + { + m_t0 = t0; + //try { + func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y))); + //} + //catch (CanteraError) { + //showErrors(); + //error("Teminating execution"); + //} + + int result, flag; + 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 != CV_SUCCESS) + throw CVodesErr("CVReInit failed. result = "+int2str(result)); + + 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 == BAND + NOJAC) { + long int N = m_neq; + long int nu = m_mupper; + long int nl = m_mlower; + CVBand(m_cvode_mem, N, nu, nl); + } + else if (m_type == GMRES) { + CVSpgmr(m_cvode_mem, PREC_NONE, 0); + } + else { + throw CVodesErr("unsupported option"); + } + + + // 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."); + if (m_np > 0) { + CVodeGetSens(m_cvode_mem, tout, m_yS); + } + } + + 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; + CVodeGetNumRhsEvals(m_cvode_mem, &ne); + return ne; + //return m_iopt[NFE]; + } + + double CVodesIntegrator::sensitivity(int k, int p) { + if (k < 0 || k >= m_neq) + throw CVodesErr("sensitivity: k out of range ("+int2str(p)+")"); + if (p < 0 || p >= m_np) + throw CVodesErr("sensitivity: p out of range ("+int2str(p)+")"); + return NV_Ith_S(m_yS[p],k); + } +} + + diff --git a/Cantera/src/numerics/CVodesIntegrator.h b/Cantera/src/numerics/CVodesIntegrator.h new file mode 100644 index 000000000..73a8516bb --- /dev/null +++ b/Cantera/src/numerics/CVodesIntegrator.h @@ -0,0 +1,118 @@ +/** + * @file CVodesWrapper.h + */ + +/* $Author$ + * $Date$ + * $Revision$ + */ + +// Copyright 2005 California Institute of Technology + + +#ifndef CT_CVODESWRAPPER_H +#define CT_CVODESWRAPPER_H + +#ifdef HAS_SUNDIALS + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "Integrator.h" +#include "FuncEval.h" +#include "ctexceptions.h" +#include "ct_defs.h" + +#ifdef SUNDIALS_VERSION_22 +#include +#else +#include +#endif + +namespace Cantera { + + class FuncData; + + /** + * Exception thrown when a CVODES error is encountered. + */ + class CVodesErr : public CanteraError { + public: + CVodesErr(std::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 setSensitivityTolerances(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); + virtual void setBandwidth(int N_Upper, int N_Lower) { + m_mupper = N_Upper; + m_mlower = N_Lower; + } + virtual int nSensParams() { return m_np; } + virtual double sensitivity(int k, int p); + + private: + + void sensInit(double t0, FuncEval& func); + + 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; + double m_reltolsens, m_abstolsens; + int m_nabs; + double m_hmax, m_hmin; + int m_maxsteps; + FuncData* m_fdata; + N_Vector* m_yS; + int m_np; + int m_mupper, m_mlower; + }; + +} // namespace + +#else + +No sundials! + +#endif + +#endif diff --git a/Cantera/src/numerics/DAE_Solver.h b/Cantera/src/numerics/DAE_Solver.h new file mode 100644 index 000000000..cb37eee20 --- /dev/null +++ b/Cantera/src/numerics/DAE_Solver.h @@ -0,0 +1,217 @@ +/** + * + * @file DAE_Solver.h + * + * Header file for class DAE_Solver + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * 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 "global.h" + +namespace Cantera { + + 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; + }; + + const int cDirect = 0; + const int cKrylov = 1; + + + + /** + * Wrapper for DAE solvers + */ + class DAE_Solver { + public: + + DAE_Solver(ResidEval& 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 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 void correctInitial_Y_given_Yp() { + warn("correctInitial_Y_given_Yp"); + } + + virtual void correctInitial_YaYp_given_Yd() { + warn("correctInitial_YaYp_given_Yd"); + } + + /** + * Solve the system of equations up to time tout. + */ + virtual int solve(doublereal tout) { + warn("solve"); return 0; + } + + /** + * Take one internal step. + */ + virtual int step(doublereal tout) { + warn("step"); return 0; + } + + /// Number of equations. + int nEquations() const { return m_resid.nEquations(); } + + /** + * 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; + + + private: + void warn(std::string msg) const { + writelog(">>>> Warning: method "+msg+" of base class " + +"DAE_Solver called. Nothing done.\n"); + } + }; + +} + +#endif diff --git a/Cantera/src/numerics/DAE_solvers.cpp b/Cantera/src/numerics/DAE_solvers.cpp new file mode 100644 index 000000000..80d8a2a1e --- /dev/null +++ b/Cantera/src/numerics/DAE_solvers.cpp @@ -0,0 +1,30 @@ + +#include "ct_defs.h" +#include "DAE_Solver.h" + +#ifdef DAE_DEVEL + +#ifdef HAS_SUNDIALS +#include "IDA_Solver.cpp" +#endif + +namespace Cantera { + + DAE_Solver* newDAE_Solver(string itype) { + if (itype == "IDA") { +#ifdef HAS_SUNDIALS + return new IDA_Solver(); +#else + raise CanteraError("newDAE_Solver","IDA solver requires sundials" + " package, but Cantera was not built with sundials."); +#endif + } + else { + throw CanteraError("newDAE_Solver", + "unknown DAE solver: "+itype); + } + } +} + +# +#endif diff --git a/Cantera/src/numerics/DASPK.cpp b/Cantera/src/numerics/DASPK.cpp new file mode 100755 index 000000000..5fedb370d --- /dev/null +++ b/Cantera/src/numerics/DASPK.cpp @@ -0,0 +1,282 @@ +/** + * @file DASPK.cpp + * + */ + +// Copyright 2001 California Institute of Technology + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "DASPK.h" +#include "ctexceptions.h" +#include "stringUtils.h" + +#include +using namespace std; + +extern "C" { + + typedef void (*ResidFunc)(const doublereal* t, + const doublereal* y, const doublereal* yprime, + const doublereal* cj, doublereal* delta, + integer* ires, doublereal* rpar, integer* ipar); + + typedef void (*JacFunc)(); + typedef void (*PsolFunc)(); + + extern void ddaspk_(ResidFunc res, integer* neq, doublereal* t, + doublereal* y, doublereal* yprime, doublereal* tout, integer* info, + doublereal* rtol, doublereal* atol, integer* idid, doublereal* rwork, + integer* lrw, integer* iwork, integer* liw, doublereal* rpar, + integer* ipar, JacFunc jac, PsolFunc psol); + + + /** + * Function called by DASPK to evaluate the residual. + */ + static void ddaspk_res(const doublereal* t, + const doublereal* y, const doublereal* yprime, + const doublereal* cj, doublereal* delta, + integer* ires, doublereal* rpar, integer* ipar) { + void **iddres_res = reinterpret_cast(&(ipar[0])); + void *hndl = *iddres_res; + Cantera::ResidEval* f = (Cantera::ResidEval*)hndl; + double delta_t = 0.0; + f->evalResid(*t, delta_t, y, yprime, delta); + } + + static void ddaspk_jac() {} + static void ddaspk_psol() {} + +} + + +namespace Cantera { + + class DASPKErr : public CanteraError { + public: + DASPKErr(string proc, string msg) + : CanteraError("DASPK::"+proc,msg) {} + virtual ~DASPKErr(){} + }; + + /** + * Constructor. Default settings: dense jacobian, no user-supplied + * Jacobian function, Newton iteration. + */ + DASPK::DASPK(ResidEval& f) : + m_resid(f), + m_idid(0), + m_lrw(0), + m_liw(0), + m_ml(0), + m_mu(0), + m_lenwp(0), + m_ok(false), + m_init(false), + m_time(0.0) + { + m_info.resize(20); + m_neq = f.neq(); + m_rwork.resize(20); // will be reset later + m_iwork.resize(20); // " + m_ipar.resize(2); + m_rpar.resize(2); + void *iddr = static_cast(&m_resid); + void **iddr_ipar = reinterpret_cast(&(m_ipar[0])); + *iddr_ipar = iddr; + setTolerances(1.e-7, 1.e-15); + } + + + /// Destructor. + DASPK::~DASPK(){} + + void DASPK::setTolerances(int nr, double* reltol, int na, double* abstol) { + // scalar tolerances + if (nr == 1 && na == 1) { + setInfo(2,0); + m_rtol.resize(1); + m_rtol[0] = reltol[0]; + m_atol.resize(1); + m_atol[0] = abstol[0]; + } + // vector tolerances + else { + setInfo(2,1); + m_rtol.resize(neq()); + m_atol.resize(neq()); + copy(reltol, reltol + nr, m_rtol.begin()); + copy(abstol, abstol + na, m_atol.begin()); + } + } + + void DASPK::setTolerances(double reltol, double abstol) { + doublereal rtol = reltol; + doublereal atol = abstol; + setTolerances(1, &rtol, 1, &atol); + } + + void DASPK::setJacobian(Jacobian& jac) { + + // No Jacobian evaluation function is supplied, so let DASPK + // compute the Jacobian by numerical finite-difference + if (!jac.supplied()) setInfo(5,0); + else { + setInfo(5,1); + } + + if (jac.isBanded()) { + setInfo(6,1); + setIwork(1, jac.lowerBandWidth()); + setIwork(2, jac.upperBandWidth()); + } + else setInfo(6,0); + } + + void DASPK::setMethod(int methodType) { + if (methodType == cDirect) + setInfo(12,0); + else if (methodType == cKrylov) + setInfo(12,1); + else + throw DASPKErr("setMethod", + "method must be either cDirect " + "or cKrylov"); + } + + void DASPK::setMaxTime(doublereal tmax) { + setInfo(4,1); + setRwork(1,tmax); + } + + void DASPK::setMaxStepSize(doublereal dtmax) { + setInfo(7,1); + setRwork(2,dtmax); + } + + void DASPK::setInitialIntStepSize(doublereal h0) { + setInfo(8,1); + setRwork(3,h0); + } + + void DASPK::setMaxOrder(int n) { + setInfo(9,1); + setIwork(3,n); + } + + void DASPK::estimateInitial_Y_given_Yp() { + setInfo(11,2); + } + + void DASPK::estimateInitial_YaYp_given_Yd( + const vector& vartypes) { + setInfo(11,2); + int m, n = neq(); + int lid = ((info(10) == 0 || info(10) == 2) ? + 41 : 41 + neq()); + if (int(m_iwork.size()) < lid + neq()) + m_iwork.resize(lid + neq()); + for (m = 0; m < n; m++) { + setIwork(lid + m, vartypes[m]); + } + } + + void DASPK::sizeRwork() { + int base; + if (info(12) == 0) { + base = 50 + 9*neq(); + if (info(6) == 0) + base += neq()*neq(); + else { + base += (2*m_ml + m_mu + 1)*neq(); + if (info(5) == 0) + base += 2*(neq()/(m_ml + m_mu + 1) + 1); + } + } + else { + base = 91 + 18*neq() + m_lenwp; + } + if (info(16) == 1) base += neq(); + + /// @todo fix this! + base = 2000000; // tmp + + m_rwork.resize(base, 0.0); + m_lrw = base; + } + + void DASPK::sizeIwork() { + int base; + if (info(12) == 0) { + base = 40 + neq(); + } + else { + base = 40 + m_lenwp; + } + if (info(10) == 1 || info(10) == 3) base += neq(); + if (info(11) == 1 || info(16) == 1) base += neq(); + m_iwork.resize(base); + m_liw = base; + } + + + void DASPK::init(doublereal t0) + { + m_init = true; + m_time = t0; + setInfo(1,0); // tells DASPK to initialize + sizeRwork(); + sizeIwork(); + //m_resid.init(t0); + } + + int DASPK::integrate(doublereal tout) { + if (!m_init) init(0.0); + + doublereal tfinal = tout; + setInfo(3,0); // don't want intermediate output + + ddaspk_(ddaspk_res, &m_neq, &m_time, m_resid.solution(), + m_resid.solution_dot(), &tfinal, m_info.begin(), + m_rtol.begin(), m_atol.begin(), &m_idid, + m_rwork.begin(), &m_lrw, m_iwork.begin(), &m_liw, + m_rpar.begin(), m_ipar.begin(), ddaspk_jac, ddaspk_psol); + + return m_idid; + } + + void DASPK::step(double tout) + { + setInfo(3,1); // do want intermediate output + doublereal tfinal = tout; + // setInfo(3,0); // don't want intermediate output + + ddaspk_(ddaspk_res, &m_neq, &m_time, m_resid.solution(), + m_resid.solution_dot(), &tfinal, m_info.begin(), + m_rtol.begin(), m_atol.begin(), &m_idid, + m_rwork.begin(), &m_lrw, m_iwork.begin(), &m_liw, + m_rpar.begin(), m_ipar.begin(), ddaspk_jac, ddaspk_psol); + if (m_idid < 0) { + throw DASPKErr("step", + "DASPK returned IDID = "+int2str(m_idid)); + m_ok = false; + } + else if (m_idid == 1 || m_idid == 2 || m_idid == 3) { + m_ok = true; + } + else { + m_ok = false; + } + return; + } + + int DASPK::nEvals() const { return iwork(12); } +} + + + diff --git a/Cantera/src/numerics/DASPK.h b/Cantera/src/numerics/DASPK.h new file mode 100755 index 000000000..d119420b1 --- /dev/null +++ b/Cantera/src/numerics/DASPK.h @@ -0,0 +1,116 @@ +/** + * + * @file DASPK.h + * + * Header file for class DASPK + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2001 California Institute of Technology + * + */ + +#ifndef CT_DASPK_H +#define CT_DASPK_H + +#include + +#include "ct_defs.h" +#include "ResidEval.h" + +namespace Cantera { + + 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 BandedJac : public Jacobian { + public: + BandedJac(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 ResidEval; + + const int cDirect = 0; + const int cKrylov = 1; + + /** + * Wrapper for DASPK 2.0 DAE solver of Petzold et al. + */ + class DASPK { + public: + + DASPK(ResidEval& f); + virtual ~DASPK(); + + integer iwork(int n) const { return m_iwork[n-1];} + doublereal rwork(int n) const { return m_rwork[n-1];} + void setIwork(int n, integer m) { m_iwork[n-1] = m; } + void setRwork(int n, doublereal v) { m_rwork[n-1] = v; } + integer info(int n) { return m_info[n-1]; } + void setInfo(int n, integer m) { m_info[n-1] = m; } + + void setTolerances(int nr, double* reltol, int na, double* abstol); + void setTolerances(double reltol, double abstol); + void setJacobian(Jacobian& jac); + void setMethod(int methodType); + void setMaxTime(doublereal tmax); + void setMaxStepSize(doublereal dtmax); + void setMaxOrder(int n); + void setInitialIntStepSize(doublereal h0); + void estimateInitial_Y_given_Yp(); + void estimateInitial_YaYp_given_Yd(const vector& vartypes); + void sizeRwork(); + void sizeIwork(); + int integrate(doublereal tout); + void step(doublereal tout); + int nEvals() const; + int neq() { return m_resid.neq(); } + void init(doublereal t0); + + protected: + + ResidEval& m_resid; + vector_int m_info; + vector_int m_iwork; + vector_int m_ipar; + + vector_fp m_rwork; + vector_fp m_atol; + vector_fp m_rtol; + vector_fp m_rpar; + + integer m_idid; + integer m_neq; + integer m_lrw; + integer m_liw; + integer m_ml; + integer m_mu; + integer m_lenwp; + + bool m_ok; + bool m_init; + doublereal m_time; + }; + +} + +#endif + diff --git a/Cantera/src/numerics/DenseMatrix.cpp b/Cantera/src/numerics/DenseMatrix.cpp new file mode 100755 index 000000000..2e83a27c8 --- /dev/null +++ b/Cantera/src/numerics/DenseMatrix.cpp @@ -0,0 +1,148 @@ +/** + * @file DenseMatrix.cpp + */ + +// Copyright 2001 California Institute of Technology + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ct_defs.h" +#include "ctlapack.h" +#include "utilities.h" +#include "DenseMatrix.h" + +namespace Cantera { + + /// assignment. + DenseMatrix& DenseMatrix::operator=(const DenseMatrix& y) { + if (&y == this) return *this; + Array2D::operator=(y); + m_ipiv = y.ipiv(); + return *this; + } + + void DenseMatrix::resize(int n, int m, doublereal v) { + Array2D::resize(n,m,v); + m_ipiv.resize( max(n,m) ); + } + + void DenseMatrix::mult(const double* b, double* prod) const { + ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, + static_cast(nRows()), + static_cast(nRows()), 1.0, ptrColumn(0), //begin(), + static_cast(nRows()), b, 1, 0.0, prod, 1); + } + + void DenseMatrix::leftMult(const double* b, double* prod) const { + int nc = static_cast(nColumns()); + int nr = static_cast(nRows()); + int n, i; + double sum = 0.0; + for (n = 0; n < nc; n++) { + sum = 0.0; + for (i = 0; i < nr; i++) { + sum += value(i,n)*b[i]; + } + prod[n] = sum; + } + } + + int solve(DenseMatrix& A, double* b) { + int info=0; + ct_dgetrf(static_cast(A.nRows()), + static_cast(A.nColumns()), A.ptrColumn(0), //begin(), + static_cast(A.nRows()), &A.ipiv()[0], info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRF returned INFO = "+int2str(info)); + ct_dgetrs(ctlapack::NoTranspose, + static_cast(A.nRows()), 1, A.ptrColumn(0), //begin(), + static_cast(A.nRows()), + &A.ipiv()[0], b, + static_cast(A.nColumns()), info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRS returned INFO = "+int2str(info)); + return 0; + } + + int solve(DenseMatrix& A, DenseMatrix& b) { + int info=0; + ct_dgetrf(static_cast(A.nRows()), + static_cast(A.nColumns()), A.ptrColumn(0), + static_cast(A.nRows()), &A.ipiv()[0], info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRF returned INFO = "+int2str(info)); + ct_dgetrs(ctlapack::NoTranspose, static_cast(A.nRows()), + static_cast(b.nColumns()), + A.ptrColumn(0), static_cast(A.nRows()), + &A.ipiv()[0], b.ptrColumn(0), + static_cast(b.nRows()), info); + if (info != 0) + throw CanteraError("DenseMatrix::solve", + "DGETRS returned INFO = "+int2str(info)); + return 0; + } + + +#ifdef INCL_LEAST_SQUARES + /** @todo fix lwork */ + int leastSquares(DenseMatrix& A, double* b) { + int info = 0; + int rank = 0; + double rcond = -1.0; + // fix this! + int lwork = 6000; // 2*(3*min(m,n) + max(2*min(m,n), max(m,n))); + vector_fp work(lwork); + vector_fp s(min(static_cast(A.nRows()), + static_cast(A.nColumns()))); + ct_dgelss(static_cast(A.nRows()), + static_cast(A.nColumns()), 1, A.ptrColumn(0), + static_cast(A.nRows()), b, + static_cast(A.nColumns()), &s[0], //.begin(), + rcond, rank, &work[0], work.size(), info); + if (info != 0) + throw CanteraError("DenseMatrix::leaseSquares", + "DGELSS returned INFO = "+int2str(info)); + return 0; + } +#endif + + void multiply(const DenseMatrix& A, const double* b, double* prod) { + ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, + static_cast(A.nRows()), static_cast(A.nColumns()), 1.0, + A.ptrColumn(0), static_cast(A.nRows()), b, 1, 0.0, prod, 1); + } + + void increment(const DenseMatrix& A, + const double* b, double* prod) { + ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, + static_cast(A.nRows()), static_cast(A.nRows()), 1.0, + A.ptrColumn(0), static_cast(A.nRows()), b, 1, 1.0, prod, 1); + } + + int invert(DenseMatrix& A, int nn) { + integer n = (nn > 0 ? nn : static_cast(A.nRows())); + int info=0; + ct_dgetrf(n, n, A.ptrColumn(0), static_cast(A.nRows()), + &A.ipiv()[0], info); + if (info != 0) + throw CanteraError("invert", + "DGETRF returned INFO="+int2str(info)); + + vector_fp work(n); + integer lwork = static_cast(work.size()); + ct_dgetri(n, A.ptrColumn(0), static_cast(A.nRows()), + &A.ipiv()[0], + &work[0], lwork, info); + if (info != 0) + throw CanteraError("invert", + "DGETRI returned INFO="+int2str(info)); + return 0; + } +} + diff --git a/Cantera/src/numerics/DenseMatrix.h b/Cantera/src/numerics/DenseMatrix.h new file mode 100755 index 000000000..bb335ab9e --- /dev/null +++ b/Cantera/src/numerics/DenseMatrix.h @@ -0,0 +1,112 @@ +/** + * @file DenseMatrix.h + * + * Dense (not sparse) matrices. + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + * + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_DENSEMATRIX_H +#define CT_DENSEMATRIX_H + +//#include +//#include +//using namespace std; + +#include "ct_defs.h" +#include "Array.h" + +namespace Cantera { + + /** + * A class for full (non-sparse) matrices with Fortran-compatible + * data storage. Adds matrix operations to class Array2D. + */ + class DenseMatrix : public Array2D { + + public: + + DenseMatrix(){} + + /** + * Constructor. Create an \c n by \c m matrix, and initialize + * all elements to \c v. + */ + DenseMatrix(int n, int m, doublereal v = 0.0) : Array2D(n,m,v) { + m_ipiv.resize( max(n, m) ); + } + + /// copy constructor + DenseMatrix(const DenseMatrix& y) : Array2D(y) { + m_ipiv = y.ipiv(); + } + + /// assignment. + DenseMatrix& operator=(const DenseMatrix& y); + + void resize(int n, int m, doublereal v = 0.0); + + /// Destructor. Does nothing. + virtual ~DenseMatrix(){} + + + /** + * Multiply A*b and write result to \c prod. + */ + virtual void mult(const double* b, double* prod) const; + + /** + * Left-multiply the matrix by transpose(b), and write the + * result to prod. + */ + virtual void leftMult(const double* b, double* prod) const; + + vector_int& ipiv() { return m_ipiv; } + const vector_int& ipiv() const { return m_ipiv; } + + protected: + + vector_int m_ipiv; + }; + + + /** + * Solve Ax = b. Array b is overwritten on exit with x. + */ + int solve(DenseMatrix& A, double* b); + + /** Solve Ax = b for multiple right-hand-side vectors. */ + int solve(DenseMatrix& A, DenseMatrix& b); + +#ifdef INCL_LEAST_SQUARES + /** @todo fix lwork */ + int leastSquares(DenseMatrix& A, double* b); +#endif + /** + * Multiply \c A*b and return the result in \c prod. Uses BLAS + * routine DGEMV. + */ + void multiply(const DenseMatrix& A, const double* b, double* prod); + + void increment(const DenseMatrix& A, + const double* b, double* prod); + + /** + * invert A. A is overwritten with A^-1. + */ + int invert(DenseMatrix& A, int nn=-1); + +} + +#endif + + + diff --git a/Cantera/src/numerics/Func1.cpp b/Cantera/src/numerics/Func1.cpp new file mode 100644 index 000000000..6fcc8d53f --- /dev/null +++ b/Cantera/src/numerics/Func1.cpp @@ -0,0 +1,405 @@ +#include "Func1.h" +#include "stringUtils.h" +#include "global.h" + +using namespace std; + +namespace Cantera { + + static Func1* checkDupl(Func1& f) { + if (f.parent() != 0) + return &f.duplicate(); + else + return &f; + } + + Func1& Sin1::derivative() const { + Func1* c = new Cos1(m_c); + Func1* r = &newTimesConstFunction(*c, m_c); +#ifdef DEBUG_FUNC + cout << "Sin1::derivative: \n"; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + r->write("x") + "\'\n"; +#endif + return *r; + } + + Func1& Cos1::derivative() const { + Func1* s = new Sin1(m_c); + Func1* r = &newTimesConstFunction(*s, -m_c); +#ifdef DEBUG_FUNC + cout << "Cos1::derivative: \n"; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + r->write("x") + "\'\n"; +#endif + return *r; + } + + Func1& Exp1::derivative() const { + Func1* f = new Exp1(m_c); + if (m_c != 1.0) + return newTimesConstFunction(*f, m_c); + else + return *f; + } + + Func1& Pow1::derivative() const { + Func1* r; + if (m_c == 0.0) { + r = new Const1(0.0); + } + else if (m_c == 1.0) { + r = new Const1(1.0); + } + else { + Func1* f = new Pow1(m_c - 1.0); + r = &newTimesConstFunction(*f, m_c); + } +#ifdef DEBUG_FUNC + cout << "Pow1::derivative: \n"; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + r->write("x") + "\'\n"; +#endif + return *r; + } + + string Func1::write(std::string arg) const { + return "("+arg+")"; + } + + string Sin1::write(string arg) const { + string c = ""; + if (m_c != 1.0) c = fp2str(m_c); + return "\\sin("+c+arg+")"; + } + + string Cos1::write(string arg) const { + string c = ""; + if (m_c != 1.0) c = fp2str(m_c); + return "\\cos("+c+arg+")"; + } + + string Pow1::write(string arg) const { + //cout << "Pow1" << endl; + string c = ""; + if (m_c == 0.5) { + return "\\sqrt{" + arg + "}"; + } + if (m_c == -0.5) { + return "\\frac{1}{\\sqrt{" + arg + "}}"; + } + if (m_c != 1.0) { + c = fp2str(m_c); + return "\\left("+arg+"\\right)^{"+c+"}"; + } + else { + return arg; + } + } + + string Exp1::write(string arg) const { + string c = ""; + if (m_c != 1.0) c = fp2str(m_c); + return "\\exp("+c+arg+")"; + } + + string Const1::write(string arg) const { + //cout << "Const1" << endl; + string c = ""; + c = fp2str(m_c); + return c; + } + + string Ratio1::write(string arg) const { + //cout << "Ratio1" << endl; + return "\\frac{" + m_f1->write(arg) + "}{" + + m_f2->write(arg) + "}"; + } + + string Product1::write(string arg) const { + //cout << "Product1" << endl; + string s = m_f1->write(arg); + if (m_f1->order() < order()) s = "\\left(" + s + "\\right)"; + string s2 = m_f2->write(arg); + if (m_f2->order() < order()) s2 = "\\left(" + s2 + "\\right)"; + return s + " " + s2; + } + + string Sum1::write(string arg) const { + //cout << "Sum1" << endl; + string s1 = m_f1->write(arg); + string s2 = m_f2->write(arg); + if (s2[0] == '-') return s1 + " - " + s2.substr(1,s2.size()); + else return s1 + " + " + s2; + } + + string Diff1::write(string arg) const { + //cout << "Diff1" << endl; + string s1 = m_f1->write(arg); + string s2 = m_f2->write(arg); + if (s2[0] == '-') return s1 + " + " + s2.substr(1,s2.size()); + else return s1 + " - " + s2; + } + + string Composite1::write(string arg) const { + //cout << "Composite1" << endl; + string g = m_f2->write(arg); + return m_f1->write(g); + } + + string TimesConstant1::write(string arg) const { + //cout << "TimesConstant1" << endl; + string s = m_f1->write(arg); + if (m_f1->order() < order()) s = "\\left(" + s + "\\right)"; + if (m_c == 1.0) return s; + if (m_c == -1.0) return "-"+s; + char n = s[0]; + if (n >= '0' && n <= '9') + s = "\\left(" + s + "\\right)"; + return fp2str(m_c) + s; + } + + string PlusConstant1::write(string arg) const { + //cout << "PlusConstant1" << endl; + if (m_c == 0.0) return m_f1->write(arg); + return m_f1->write(arg) + " + " + fp2str(m_c); + } + + doublereal Func1::isProportional(TimesConstant1& other) { + if (isIdentical(other.func1())) return other.c(); + return 0.0; + } + doublereal Func1::isProportional(Func1& other) { + if (isIdentical(other)) return 1.0; + else return 0.0; + } + + static bool isConstant(Func1& f) { + if (f.ID() == ConstFuncType) + return true; + else + return false; + } + + static bool isZero(Func1& f) { + if (f.ID() == ConstFuncType && f.c() == 0.0) + return true; + else + return false; + } + + static bool isOne(Func1& f) { + if (f.ID() == ConstFuncType && f.c() == 1.0) + return true; + else + return false; + } + + static bool isTimesConst(Func1& f) { + if (f.ID() == TimesConstantFuncType) + return true; + else + return false; + } + + static bool isExp(Func1& f) { + if (f.ID() == ExpFuncType) + return true; + else + return false; + } + + static bool isPow(Func1& f) { + if (f.ID() == PowFuncType) + return true; + else + return false; + } + + Func1& newSumFunction(Func1& f1, Func1& f2) { + if (f1.isIdentical(f2)) + return newTimesConstFunction(f1, 2.0); + if (isZero(f1)) { + delete &f1; + return f2; + } + if (isZero(f2)) { + delete &f2; + return f1; + } + doublereal c = f1.isProportional(f2); + if (c != 0) { + if (c == -1.0) + return *(new Const1(0.0)); + else { + return newTimesConstFunction(f1, c + 1.0); + } + } + return *(new Sum1(f1, f2)); + } + + Func1& newDiffFunction(Func1& f1, Func1& f2) { + if (isZero(f2)) { + delete &f2; return f1; + } + if (f1.isIdentical(f2)) { + delete &f1; delete &f2; + return *(new Const1(0.0)); + } + doublereal c = f1.isProportional(f2); + if (c != 0.0) { + if (c == 1.0) return *(new Const1(0.0)); + else return newTimesConstFunction(f1, 1.0 - c); + } + return *(new Diff1(f1, f2)); + } + + Func1& newProdFunction(Func1& f1, Func1& f2) { + if (isOne(f1)) { + delete &f1; return f2; + } + if (isOne(f2)) { + delete &f2; return f1; + } + if (isZero(f1) || isZero(f2)) { + delete &f1; delete &f2; + return *(new Const1(0.0)); + } + if (isConstant(f1) && isConstant(f2)) { + doublereal c1c2 = f1.c() * f2.c(); + delete &f1; delete &f2; + return *(new Const1(c1c2)); + } + if (isConstant(f1)) { + doublereal c = f1.c(); + delete &f1; + return newTimesConstFunction(f2, c); + } + if (isConstant(f2)) { + doublereal c = f2.c(); + delete &f2; + return newTimesConstFunction(f1, c); + } + + if (isPow(f1) && isPow(f2)) { + Func1& p = *(new Pow1(f1.c() + f2.c())); + delete &f1; delete &f2; + return p; + } + + if (isExp(f1) && isExp(f2)) { + Func1& p = *(new Exp1(f1.c() + f2.c())); + delete &f1; delete &f2; + return p; + } + + bool tc1 = isTimesConst(f1); + bool tc2 = isTimesConst(f2); + + if (tc1 || tc2) { + doublereal c1 = 1.0, c2 = 1.0; + Func1 *ff1 = 0, *ff2 = 0; + if (tc1) { + c1 = f1.c(); + ff1 = &f1.func1_dup(); + delete &f1; + } + else ff1 = &f1; + if (tc2) { + c2 = f2.c(); + ff2 = &f2.func1_dup(); + delete &f2; + } + else ff2 = &f2; + Func1& p = newProdFunction(*ff1, *ff2); + + if (c1*c2 != 1.0) { + return newTimesConstFunction(p, c1*c2); + } + else + return p; + } + else + return *(new Product1(f1, f2)); + } + + Func1& newRatioFunction(Func1& f1, Func1& f2) { + if (isOne(f2)) return f1; + if (isZero(f1)) return *(new Const1(0.0)); + if (f1.isIdentical(f2)) { + delete &f1; delete &f2; + return *(new Const1(1.0)); + } + if (f1.ID() == PowFuncType && f2.ID() == PowFuncType) { + return *(new Pow1(f1.c() - f2.c())); + } + if (f1.ID() == ExpFuncType && f2.ID() == ExpFuncType) { + return *(new Exp1(f1.c() - f2.c())); + } + return *(new Ratio1(f1, f2)); + } + + Func1& newCompositeFunction(Func1& f1, Func1& f2) { + //#ifdef DEBUG_FUNC + //cout << "creating new composite function." << endl; + //cout << "f1 = " << f1.write("x") << " " << f1.ID() << endl; + //cout << "f2 = " << f2.write("x") << " " << f2.ID() << endl; + //#endif + if (isZero(f1)) { + delete &f1; delete &f2; + return *(new Const1(0.0)); + } + if (isConstant(f1)) { + delete &f2; + return f1; + } + if (isPow(f1) && f1.c() == 1.0) { + delete &f1; + return f2; + } + if (isPow(f1) && f1.c() == 0.0) { + delete &f1; + delete &f2; + return *(new Const1(1.0)); + } + if (isPow(f1) && isPow(f2)) { + doublereal c1c2 = f1.c() * f2.c(); + delete &f1; + delete &f2; + return *(new Pow1(c1c2)); + } + return *(new Composite1(f1, f2)); + } + + Func1& newTimesConstFunction(Func1& f, doublereal c) { + if (c == 0.0) { + delete &f; + return *(new Const1(0.0)); + } + if (c == 1.0) { + return f; + } + if (f.ID() == TimesConstantFuncType) { + f.setC(f.c() * c); + return f; + } + return *(new TimesConstant1(f, c)); + } + + Func1& newPlusConstFunction(Func1& f, doublereal c) { + if (c == 0.0) { + return f; + } + if (isConstant(f)) { + doublereal cc = f.c() + c; + delete &f; + return *(new Const1(cc)); + } + if (f.ID() == PlusConstantFuncType) { + f.setC(f.c() + c); + return f; + } + return *(new PlusConstant1(f, c)); + } + +} diff --git a/Cantera/src/numerics/Func1.h b/Cantera/src/numerics/Func1.h new file mode 100644 index 000000000..246845d3e --- /dev/null +++ b/Cantera/src/numerics/Func1.h @@ -0,0 +1,651 @@ +/** + * @file Func1.h + * + * $Author$ + * $Date$ + * $Revision$ + */ + +// Copyright 2001 California Institute of Technology + +#ifndef CT_FUNC1_H +#define CT_FUNC1_H + +#undef DEBUG_FUNC + + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ct_defs.h" + +#include +#include +using namespace std; + +namespace Cantera { + + const int FourierFuncType = 1; + const int PolyFuncType = 2; + const int ArrheniusFuncType = 3; + const int GaussianFuncType = 4; + const int SumFuncType = 20; + const int DiffFuncType = 25; + const int ProdFuncType = 30; + const int RatioFuncType = 40; + const int PeriodicFuncType = 50; + const int CompositeFuncType = 60; + const int TimesConstantFuncType = 70; + const int PlusConstantFuncType = 80; + const int SinFuncType = 100; + const int CosFuncType = 102; + const int ExpFuncType = 104; + const int PowFuncType = 106; + const int ConstFuncType = 110; + + class Sin1; + class Cos1; + class Exp1; + class Pow1; + class TimesConstant1; + + /** + * Base class for 'functor' classes that evaluate a function of + * one variable. + */ + class Func1 { + public: + Func1() : m_c(0.0), m_f1(0), m_f2(0), m_parent(0) {} + virtual ~Func1() {} + virtual int ID() const { return 0; } + + virtual Func1& duplicate() { cout << "DUPL ERR: ID = " << ID() << endl; + return *(new Func1);} + + /// Calls method eval to evaluate the function + doublereal operator()(doublereal t) const { return eval(t); } + + /// Evaluate the function. + virtual doublereal eval(doublereal t) const { return 0.0; } + + virtual Func1& derivative() const { + cout << "derivative error... ERR: ID = " << ID() << endl; + cout << write("x") << endl; + return *(new Func1); + } + + bool isIdentical(Func1& other) const { + if ((ID() != other.ID()) || (m_c != other.m_c)) + return false; + if (m_f1) { + if (!other.m_f1) return false; + if (!m_f1->isIdentical(*other.m_f1)) return false; + } + if (m_f2) { + if (!other.m_f2) return false; + if (!m_f2->isIdentical(*other.m_f2)) return false; + } + return true; + } + + virtual doublereal isProportional(TimesConstant1& other); + virtual doublereal isProportional(Func1& other); + + virtual std::string write(std::string arg) const; + + doublereal c() const { return m_c; } + void setC(doublereal c) { m_c = c; } + Func1& func1() { return *m_f1; } + Func1& func2() { return *m_f2; } + virtual int order() const { return 3; } + Func1& func1_dup() const { return m_f1->duplicate(); } + Func1& func2_dup() const { return m_f2->duplicate(); } + Func1* parent() { return m_parent; } + void setParent(Func1* p) { m_parent = p; } + + protected: + doublereal m_c; + Func1 *m_f1, *m_f2; + Func1* m_parent; + + private: + }; + + + Func1& newSumFunction(Func1& f1, Func1& f2); + Func1& newDiffFunction(Func1& f1, Func1& f2); + Func1& newProdFunction(Func1& f1, Func1& f2); + Func1& newRatioFunction(Func1& f1, Func1& f2); + Func1& newCompositeFunction(Func1& f1, Func1& f2); + Func1& newTimesConstFunction(Func1& f1, doublereal c); + Func1& newPlusConstFunction(Func1& f1, doublereal c); + + /// sin + class Sin1 : public Func1 { + public: + Sin1(doublereal omega = 1.0) { + m_c = omega; + } + virtual ~Sin1() {} + virtual std::string write(std::string arg) const; + virtual Func1& duplicate() { return *(new Sin1(m_c)); } + virtual int ID() const { return SinFuncType; } + virtual doublereal eval(doublereal t) const { + return sin(m_c*t); + } + virtual Func1& derivative() const; + + protected: + + }; + + /// cos + class Cos1 : public Func1 { + public: + Cos1(doublereal omega = 1.0) { + m_c = omega; + } + virtual ~Cos1() {} + virtual std::string write(std::string arg) const; + virtual Func1& duplicate() { return *(new Cos1(m_c)); } + virtual int ID() const { return CosFuncType; } + virtual doublereal eval(doublereal t) const { + return cos(m_c * t); + } + virtual Func1& derivative() const; + + protected: + }; + + /// exp + class Exp1 : public Func1 { + public: + Exp1(doublereal A = 1.0) {m_c = A;} + virtual ~Exp1() {} + virtual std::string write(std::string arg) const; + virtual int ID() const { return ExpFuncType; } + virtual Func1& duplicate() { return *(new Exp1(m_c)); } + virtual doublereal eval(doublereal t) const { + return exp(m_c*t); + } + + virtual Func1& derivative() const; + + protected: + + }; + + /// pow + class Pow1 : public Func1 { + public: + Pow1(doublereal n) {m_c = n;} + virtual ~Pow1() {} + virtual std::string write(std::string arg) const; + virtual int ID() const { return PowFuncType; } + virtual Func1& duplicate() { return *(new Pow1(m_c)); } + virtual doublereal eval(doublereal t) const { + return pow(t, m_c); + } + virtual Func1& derivative() const; + + protected: + + }; + + /** + * Constant. + */ + class Const1 : public Func1 { + public: + Const1(doublereal A) { + m_c = A; + } + virtual ~Const1() {} + virtual std::string write(std::string arg) const; + virtual int ID() const { return ConstFuncType; } + virtual doublereal eval(doublereal t) const { + return m_c; + } + virtual Func1& duplicate() { return *(new Const1(m_c)); } + virtual Func1& derivative() const { + Func1* z = new Const1(0.0); + return *z; + } + + protected: + }; + + + + /** + * Sum of two functions. + */ + class Sum1 : public Func1 { + public: + Sum1(Func1& f1, Func1& f2) { + m_f1 = &f1; + m_f2 = &f2; + m_f1->setParent(this); + m_f2->setParent(this); + } + virtual ~Sum1() { + delete m_f1; + delete m_f2; + } + virtual int ID() const { return SumFuncType; } + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) + m_f2->eval(t); + } + virtual Func1& duplicate() { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newSumFunction(f1d, f2d); + return dup; + } + virtual Func1& derivative() const { + Func1& d1 = m_f1->derivative(); + Func1& d2 = m_f2->derivative(); + Func1& d = newSumFunction(d1, d2); + return d; + } + virtual int order() const { return 0; } + virtual std::string write(std::string arg) const; + protected: + }; + + + /** + * Difference of two functions. + */ + class Diff1 : public Func1 { + public: + Diff1(Func1& f1, Func1& f2) { + m_f1 = &f1; + m_f2 = &f2; + } + virtual ~Diff1() { + delete m_f1; + delete m_f2; + } + virtual int ID() const { return DiffFuncType; } + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) - m_f2->eval(t); + } + virtual Func1& duplicate() { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newDiffFunction(f1d, f2d); + return dup; + } + virtual Func1& derivative() const { + Func1& d = newDiffFunction(m_f1->derivative(), m_f2->derivative()); + return d; + } + virtual int order() const { return 0; } + virtual std::string write(std::string arg) const; + + protected: + }; + + + /** + * Product of two functions. + */ + class Product1 : public Func1 { + public: + Product1(Func1& f1, Func1& f2) { + m_f1 = &f1; + m_f2 = &f2; + } + + virtual ~Product1() { + delete m_f1; + delete m_f2; + } + virtual int ID() const { return ProdFuncType; } + virtual Func1& duplicate() { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newProdFunction(f1d, f2d); + return dup; + } + virtual std::string write(std::string arg) const; + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) * m_f2->eval(t); + } + virtual Func1& derivative() const { + Func1& a1 = newProdFunction(m_f1->duplicate(), m_f2->derivative()); + Func1& a2 = newProdFunction(m_f2->duplicate(), m_f1->derivative()); + Func1& s = newSumFunction(a1, a2); + return s; + } + virtual int order() const { return 1; } + + protected: + }; + + /** + * Product of two functions. + */ + class TimesConstant1 : public Func1 { + public: + TimesConstant1(Func1& f1, doublereal A) { + m_f1 = &f1; + m_c = A; + } + + virtual ~TimesConstant1() { + delete m_f1; + } + virtual int ID() const { return TimesConstantFuncType; } + virtual Func1& duplicate() { + Func1& f1 = m_f1->duplicate(); + Func1* dup = new TimesConstant1(f1, m_c); + return *dup; + } + virtual doublereal isProportional(TimesConstant1& other) { + if (func1().isIdentical(other.func1())) + return (other.c()/c()); + else + return 0.0; + } + virtual doublereal isProportional(Func1& other) { + if (func1().isIdentical(other)) return 1.0/c(); + else return 0.0; + } + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) * m_c; + } + virtual Func1& derivative() const { + Func1& f1d = m_f1->derivative(); + Func1* d = &newTimesConstFunction(f1d, m_c); + return *d; + } + virtual std::string write(std::string arg) const; + virtual int order() const { return 0; } + protected: + }; + + /** + * A function plus a constant. + */ + class PlusConstant1 : public Func1 { + public: + PlusConstant1(Func1& f1, doublereal A) { + m_f1 = &f1; + m_c = A; + } + + virtual ~PlusConstant1() { + delete m_f1; + } + virtual int ID() const { return PlusConstantFuncType; } + virtual Func1& duplicate() { + Func1& f1 = m_f1->duplicate(); + Func1* dup = new PlusConstant1(f1, m_c); + return *dup; + } + + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) + m_c; + } + virtual Func1& derivative() const { + Func1& f1d = m_f1->derivative(); + return f1d; + } + virtual std::string write(std::string arg) const; + virtual int order() const { return 0; } + + protected: + }; + + + /** + * Ratio of two functions. + */ + class Ratio1 : public Func1 { + public: + Ratio1(Func1& f1, Func1& f2) { + m_f1 = &f1; + m_f2 = &f2; + } + virtual ~Ratio1() { + delete m_f1; + delete m_f2; + } + virtual int ID() const { return RatioFuncType; } + virtual doublereal eval(doublereal t) const { + return m_f1->eval(t) / m_f2->eval(t); + } + virtual Func1& duplicate() { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newRatioFunction(f1d, f2d); + return dup; + } + virtual Func1& derivative() const { + Func1& a1 = newProdFunction(m_f1->derivative(), m_f2->duplicate()); + Func1& a2 = newProdFunction(m_f1->duplicate(), m_f2->derivative()); + Func1& s = newDiffFunction(a1, a2); + Func1& p = newProdFunction(m_f2->duplicate(), m_f2->duplicate()); + Func1& r = newRatioFunction(s, p); + return r; + } + virtual std::string write(std::string arg) const; + virtual int order() const { return 1; } + + protected: + }; + + /** + * Composite function. + */ + class Composite1 : public Func1 { + public: + Composite1(Func1& f1, Func1& f2) { + m_f1 = &f1; + m_f2 = &f2; + } + virtual ~Composite1() { + delete m_f1; + delete m_f2; + } + virtual int ID() const { return CompositeFuncType; } + virtual doublereal eval(doublereal t) const { + return m_f1->eval( m_f2->eval(t) ); + } + virtual Func1& duplicate() { + Func1& f1d = m_f1->duplicate(); + Func1& f2d = m_f2->duplicate(); + Func1& dup = newCompositeFunction(f1d, f2d); + return dup; + } + virtual Func1& derivative() const { + Func1* d1 = &m_f1->derivative(); + + Func1* d3 = &newCompositeFunction(*d1, m_f2->duplicate()); + Func1* d2 = &m_f2->derivative(); + Func1* p = &newProdFunction(*d3, *d2); +#ifdef DEBUG_FUNC + cout << "Composite1::derivative: \n"; + cout << "f1 = " << m_f1->write("x") << endl; + cout << "f2 = " << m_f2->write("x") << endl; + cout << "d1 = " << d1 << " " << d1->write("x") << endl; + cout << "d3 = " << d3->write("x") << endl; + cout << "d2 = " << d2->write("x") << endl; + cout << "function = \'" + write("x") + "\'\n"; + cout << "derivative = \'" + p->write("x") + "\'\n"; +#endif + return *p; + } + virtual std::string write(std::string arg) const; + virtual int order() const { return 2; } + protected: + }; + + // + // The functors below are the old-style ones. They still work, + // but can't do derivatives. + // + + /** + * A Gaussian. + * \f[ + * f(t) = A e^{-[(t - t_0)/\tau]^2} + * \f] + * where \f[ \tau = \frac{fwhm}{2\sqrt{\ln 2}} \f] + * @param A peak value + * @param t0 offset + * @param fwhm full width at half max + */ + class Gaussian : public Func1 { + public: + Gaussian(double A, double t0, double fwhm) { + m_A = A; + m_t0 = t0; + m_tau = fwhm/(2.0*std::sqrt(std::log(2.0))); + } + virtual ~Gaussian() {} + virtual doublereal eval(doublereal t) const { + doublereal x = (t - m_t0)/m_tau; + return m_A*std::exp(-x*x); + } + protected: + doublereal m_A, m_t0, m_tau; + private: + }; + + + /** + * Polynomial of degree n. + */ + class Poly1 : public Func1 { + public: + Poly1(int n, doublereal* c) { + m_n = n+1; + m_c.resize(n+1); + std::copy(c, c+m_n, m_c.begin()); + } + virtual ~Poly1() {} + + virtual doublereal eval(doublereal t) const { + int n; + doublereal r = m_c[m_n-1]; + for (n = 1; n < m_n; n++) { + r *= t; + r += m_c[m_n - n - 1]; + } + return r; + } + + protected: + int m_n; + vector_fp m_c; + }; + + + /** + * Fourier cosine/sine series. + * + * \f[ + * f(t) = \frac{A_0}{2} + + * \sum_{n=1}^N A_n \cos (n \omega t) + B_n \sin (n \omega t) + * \f] + */ + class Fourier1 : public Func1 { + public: + Fourier1(int n, doublereal omega, doublereal a0, + doublereal* a, doublereal* b) { + m_n = n; + m_omega = omega; + m_a0_2 = 0.5*a0; + m_ccos.resize(n); + m_csin.resize(n); + std::copy(a, a+n, m_ccos.begin()); + std::copy(b, b+n, m_csin.begin()); + } + virtual ~Fourier1() {} + + virtual doublereal eval(doublereal t) const { + int n, nn; + doublereal sum = m_a0_2; + for (n = 0; n < m_n; n++) { + nn = n + 1; + sum += m_ccos[n]*std::cos(m_omega*nn*t) + + m_csin[n]*std::sin(m_omega*nn*t); + } + return sum; + } + + protected: + int m_n; + doublereal m_omega, m_a0_2; + vector_fp m_ccos, m_csin; + }; + + + /** + * Sum of Arrhenius terms. + * \f[ + * f(T) = \sum_{n=1}^N A_n T^b_n \exp(-E_n/T) + * \f] + */ + class Arrhenius1 : public Func1 { + public: + Arrhenius1(int n, doublereal* c) { + m_n = n; + m_A.resize(n); + m_b.resize(n); + m_E.resize(n); + int loc; + for (int i = 0; i < n; i++) { + loc = 3*i; + m_A[i] = c[loc]; + m_b[i] = c[loc+1]; + m_E[i] = c[loc+2]; + } + } + virtual ~Arrhenius1() {} + + virtual doublereal eval(doublereal t) const { + int n; + doublereal sum = 0.0; + for (n = 0; n < m_n; n++) { + sum += m_A[n]*std::pow(t,m_b[n])*std::exp(-m_E[n]/t); + } + return sum; + } + + protected: + int m_n; + vector_fp m_A, m_b, m_E; + }; + + /** + * Periodic function. Takes any function and makes it + * periodic with period T. + */ + class Periodic1 : public Func1 { + public: + Periodic1(Func1& f, doublereal T) { + m_func = &f; + m_c = T; + } + virtual ~Periodic1() { delete m_func; } + virtual doublereal eval(doublereal t) const { + int np = int(t/m_c); + doublereal time = t - np*m_c; + return m_func->eval(time); + } + protected: + Func1* m_func; + + private: + }; + +} + + +#endif diff --git a/Cantera/src/numerics/FuncEval.h b/Cantera/src/numerics/FuncEval.h new file mode 100755 index 000000000..d88181533 --- /dev/null +++ b/Cantera/src/numerics/FuncEval.h @@ -0,0 +1,69 @@ +/** + * @file FuncEval.h + * + */ + +// Copyright 2001 California Institute of Technology + +#ifndef CT_FUNCEVAL_H +#define CT_FUNCEVAL_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ct_defs.h" + +namespace Cantera { + + + /** + * Virtual base class for ODE right-hand-side function evaluators. + * Classes derived from FuncEval evaluate the right-hand-side function + * \f$ \vec{F}(t,\vec{y})\f$ in + * \f[ + * \dot{\vec{y}} = \vec{F}(t,\vec{y}). + * \f] + * @ingroup odeGroup + */ + class FuncEval { + + public: + + FuncEval() {} + virtual ~FuncEval() {} + + /** + * Evaluate the right-hand-side function. Called by the + * integrator. + * @param t time. (input) + * @param y solution vector. (input) + * @param ydot rate of change of solution vector. (output) + * @param p parameter vector + */ + virtual void eval(double t, double* y, double* ydot, double* p)=0; + + /** + * Fill the solution vector with the initial conditions + * at initial time t0. + */ + virtual void getInitialConditions(double t0, size_t leny, double* y)=0; + + /** + * Number of equations. + */ + virtual int neq()=0; + + /// Number of parameters. + virtual int nparams() { return 0; } + + protected: + + private: + + }; + +} + +#endif diff --git a/Cantera/src/numerics/IDA_Solver.cpp b/Cantera/src/numerics/IDA_Solver.cpp new file mode 100644 index 000000000..0458477b3 --- /dev/null +++ b/Cantera/src/numerics/IDA_Solver.cpp @@ -0,0 +1,256 @@ + +/** + * @file IDA_Solver.cpp + * + */ + +// Copyright 2006 California Institute of Technology + +#include "IDA_Solver.h" +#include "stringUtils.h" + +#include +using namespace std; + +#include +#include +#include +#include +#include +#include +#include + +inline static N_Vector nv(void* 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 { + + + + + public: + ResidData(ResidEval* f, int npar = 0) { + m_func = f; + } + virtual ~ResidData() {} + ResidEval* 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; + } + +} + +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; + } + + 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) 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; + } + + 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(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; + } + + doublereal IDA_Solver::getOutputParameter(int flag) { + switch (flag) { + case REAL_WORKSPACE_SIZE: + flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw); + return doublereal(lenrw); + } + +} + + diff --git a/Cantera/src/numerics/IDA_Solver.h b/Cantera/src/numerics/IDA_Solver.h new file mode 100644 index 000000000..2820b723e --- /dev/null +++ b/Cantera/src/numerics/IDA_Solver.h @@ -0,0 +1,137 @@ +/** + * + * @file IDA_Solver.h + * + * Header file for class IDA_Solver + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2006 California Institute of Technology + * + */ + +#ifndef CT_IDA_Solver_H +#define CT_IDA_Solver_H + +#include + +#include "DAE_Solver.h" +#include "ctexceptions.h" + +namespace Cantera { + + /** + * Exception thrown when a IDA error is encountered. + */ + class IDA_Err : public CanteraError { + public: + IDA_Err(string msg) : CanteraError("IDA_Solver", msg){} + }; + + + class ResidData; // forward reference + + class IDA_Solver : public DAE_Solver { + public: + + 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 setMaxStepSize(doublereal dtmax); + + virtual void setMaxOrder(int n); + + virtual void setMaxNumSteps(int n); + virtual void setInitialStepSize(doublereal h0); + virtual void setMaxStepSize(doublereal hmax); + 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); + + + /** + * 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); + + /** + * 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); + + + virtual int solve(doublereal tout); + + virtual int step(doublereal tout); + + virtual void init(doublereal t0); + + /// the current value of solution component k. + 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; + + protected: + + int m_neq; + void* m_ida_mem; + double m_t0; + void *m_y, *m_ydot, *m_id, *m_constraints, *m_abstol; + int m_type; + int m_itol; + int m_iter; + double m_reltol; + double m_abstols; + int m_nabs; + double m_hmax, m_hmin; + int m_maxsteps, m_maxord; + ResidData* m_fdata; + int m_mupper, m_mlower; + }; + +} + +#endif + diff --git a/Cantera/src/numerics/Integrator.h b/Cantera/src/numerics/Integrator.h new file mode 100755 index 000000000..c2c0f1800 --- /dev/null +++ b/Cantera/src/numerics/Integrator.h @@ -0,0 +1,188 @@ +/** + * @file Integrator.h + * + * $Author$ + * $Date$ + * $Revision$ + * + */ + +/** + * @defgroup odeGroup ODE Integrators + */ + +// Copyright 2001 California Institute of Technology + + +#ifndef CT_INTEGRATOR_H +#define CT_INTEGRATOR_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "FuncEval.h" + +#include "ct_defs.h" +#include "global.h" + +#define DIAG 1 +#define DENSE 2 +#define NOJAC 4 +#define JAC 8 +#define GMRES 16 +#define BAND 32 + +namespace Cantera { + + /** + * Specifies the method used to integrate the system of equations. + * Not all methods are supported by all integrators. + */ + enum MethodType { + BDF_Method, /**< Backward Differentiation */ + Adams_Method /**< Adams */ + }; + + /** + * Specifies the method used for iteration. + * Not all methods are supported by all integrators. + */ + enum IterType { + Newton_Iter, /**< Newton iteration */ + Functional_Iter /**< Functional iteration */ + }; + + + /** + * Abstract base class for ODE system integrators. + * @ingroup odeGroup + */ + class Integrator { + + public: + + virtual ~Integrator() {} + + /** Set or reset the number of equations. */ + //virtual void resize(int n)=0; + + /** + * Set error tolerances. + * @param reltol scalar relative tolerance + * @param number of equations + * @param abstol array of N absolute tolerance values + */ + virtual void setTolerances(doublereal reltol, int n, + doublereal* abstol) { warn("setTolerances"); } + + /** + * Set error tolerances. + * @param reltol scalar relative tolerance + * @param abstol scalar absolute tolerance + */ + virtual void setTolerances(doublereal reltol, doublereal abstol) + { warn("setTolerances"); } + + virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol) + {}// { warn("setSensitivityTolerances"); } + + /** + * Set problem type. + */ + virtual void setProblemType(int probtype) { warn("setProblemType"); } + + /** + * Initialize the integrator for a new problem. Call after + * all options have been set. + * @param t0 initial time + * @param func RHS evaluator object for system of equations. + */ + virtual void initialize(doublereal t0, FuncEval& func) + { warn("initialize"); } + + 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. + */ + virtual void integrate(doublereal tout) + { warn("integrate"); } + + /** + * Integrate the system of equations. + * @param tout integrate to this time. Note that this is the + * absolute time value, not a time interval. + */ + virtual doublereal step(doublereal tout) + { warn("step"); return 0.0; } + + /** The current value of the solution of equation k. */ + virtual doublereal& solution(int k) + { warn("solution"); return m_dummy; } + + /** The current value of the solution of the system of equations. */ + virtual doublereal* solution() + { warn("solution"); return 0; } + + /** The number of equations. */ + virtual int nEquations() const + { warn("nEquations"); return 0; } + + /** The number of function evaluations. */ + virtual int nEvals() const + { warn("nEvals"); return 0; } + + /** Set the maximum integration order that will be used. **/ + virtual void setMaxOrder(int n) + { warn("setMaxorder"); } + + /** Set the solution method */ + virtual void setMethod(MethodType t) + { warn("setMethodType"); } + + /** Set the linear iterator. */ + virtual void setIterator(IterType t) + { warn("setInterator"); } + + /** Set the maximum step size */ + virtual void setMaxStepSize(double hmax) + { warn("setMaxStepSize"); } + + /** Set the minimum step size */ + virtual void setMinStepSize(double hmin) + { warn("setMinStepSize"); } + + virtual void setMaxSteps(int nmax) + { warn("setMaxStep"); } + + virtual void setBandwidth(int N_Upper, int N_Lower) + { warn("setBandwidth"); } + + virtual int nSensParams() + { warn("nSensParams()"); return 0; } + + virtual double sensitivity(int k, int p) { + warn("sensitivity"); return 0.0; + } + + private: + + doublereal m_dummy; + void warn(std::string msg) const { + writelog(">>>> Warning: method "+msg+" of base class " + +"Integrator called. Nothing done.\n"); + } + + }; + + // defined in ODE_integrators.cpp + Integrator* newIntegrator(std::string itype); + +} // namespace + +#endif diff --git a/Cantera/src/numerics/ODE_integrators.cpp b/Cantera/src/numerics/ODE_integrators.cpp new file mode 100644 index 000000000..8aaf4ad8b --- /dev/null +++ b/Cantera/src/numerics/ODE_integrators.cpp @@ -0,0 +1,29 @@ +#include "ct_defs.h" +#include "Integrator.h" + +#ifdef NO_SUNDIALS +#undef HAS_SUNDIALS +#endif + +#ifdef HAS_SUNDIALS +#include "CVodesIntegrator.cpp" +#else +#include "CVode.cpp" +#endif + +namespace Cantera { + + Integrator* newIntegrator(std::string itype) { + if (itype == "CVODE") { +#ifdef HAS_SUNDIALS + return new CVodesIntegrator(); +#else + return new CVodeInt(); +#endif + } + else { + throw CanteraError("newIntegrator", + "unknown ODE integrator: "+itype); + } + } +} diff --git a/Cantera/src/numerics/ResidEval.h b/Cantera/src/numerics/ResidEval.h new file mode 100755 index 000000000..d2c3be01f --- /dev/null +++ b/Cantera/src/numerics/ResidEval.h @@ -0,0 +1,100 @@ +/** + * @file ResidEval.h + * + */ + +// Copyright 2006 California Institute of Technology + +#ifndef CT_RESIDEVAL_H +#define CT_RESIDEVAL_H + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ct_defs.h" + +namespace Cantera { + + const int c_NONE = 0; + const int c_GE_ZERO = 1; + const int c_GT_ZERO = 2; + const int c_LE_ZERO = -1; + const int c_LT_ZERO = -2; + + /** + * Virtual base class for DAE residual function evaluators. + * Classes derived from ResidEval evaluate the residual function + * \f[ + \vec{F}(t,\vec{y}, \vec{y^\prime}) + * \f] + * The DAE solver attempts to find a solution y(t) such that F = 0. + * @ingroup DAE_Group + */ + class ResidEval { + + public: + + ResidEval() {} + virtual ~ResidEval() {} + + /** + * Constrain solution component k. Possible values for + * 'flag' are: + * - c_NONE no constraint + * - c_GE_ZERO >= 0 + * - c_GT_ZERO > 0 + * - c_LE_ZERO <= 0 + * - c_LT_ZERO < 0 + */ + virtual void constrain(int k, int flag) { m_constrain[k] = flag; } + int constraint(int k) { return m_constrain[k]; } + + /** + * Specify that solution component k is purely algebraic - + * that is, the derivative of this component does not appear + * in the residual function. + */ + virtual void setAlgebraic(int k) { m_alg[k] = 1; } + virtual bool isAlgebraic(int k) {return (m_alg[k] == 1); } + + + /** + * 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(double t, const double* y, + const double* ydot, double* r)=0; + + /** + * Fill the solution and derivative vectors with the initial + * conditions at initial time t0. If these do not satisfy the + * residual equation, call one of the "corrrectInitial_xxx" + * methods before calling solve. + */ + virtual void getInitialConditions(double t0, double* y, + doublereal* ydot)=0; + + /** + * Number of equations. + */ + virtual int nEquations()=0; + + + protected: + + std::map m_alg; + std::map m_constrain; + + private: + + }; + +} + +#endif diff --git a/Cantera/src/numerics/ctlapack.h b/Cantera/src/numerics/ctlapack.h new file mode 100755 index 000000000..436fdb9c0 --- /dev/null +++ b/Cantera/src/numerics/ctlapack.h @@ -0,0 +1,276 @@ +/** + * @file ctlapack.h + */ + +/* $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology. + +#ifndef CT_CTLAPACK_H +#define CT_CTLAPACK_H + +#ifdef DARWIN +#undef USE_CBLAS +#undef NO_FTN_STRING_LEN_AT_END +#endif + +#include "ct_defs.h" + +//#include + +// map BLAS names to names with or without a trailing underscore. +#ifndef LAPACK_FTN_TRAILING_UNDERSCORE + +#define _DGEMV_ dgemv +#define _DGETRF_ dgetrf +#define _DGETRS_ dgetrs +#define _DGETRI_ dgetri +#define _DGELSS_ dgelss +#define _DGBSV_ dgbsv +#define _DGBTRF_ dgbtrf +#define _DGBTRS_ dgbtrs + +#define _DSCAL_ dscal + +#else + +#define _DGEMV_ dgemv_ +#define _DGETRF_ dgetrf_ +#define _DGETRS_ dgetrs_ +#define _DGETRI_ dgetri_ +#define _DGELSS_ dgelss_ +#define _DGBSV_ dgbsv_ +#define _DGBTRF_ dgbtrf_ +#define _DGBTRS_ dgbtrs_ + +#define _DSCAL_ dscal_ + +#endif + + +namespace ctlapack { + typedef enum {Transpose = 1, NoTranspose = 0} transpose_t; + typedef enum {ColMajor = 1, RowMajor = 0} storage_t; +} +const char no_yes[2] = {'N', 'T'}; + +#ifdef USE_CBLAS +#include +const CBLAS_ORDER cblasOrder[2] = { CblasRowMajor, CblasColMajor }; +const CBLAS_TRANSPOSE cblasTrans[2] = { CblasNoTrans, CblasTrans }; +#endif + +//#ifdef DARWIN +//#include +//#else + +// C interfaces for Fortran Lapack routines +extern "C" { + +#ifdef LAPACK_FTN_STRING_LEN_AT_END + int _DGEMV_(const char* transpose, + const integer* m, const integer* n, const doublereal* alpha, + const doublereal* a, const integer* lda, const doublereal* x, + const integer* incX, const doublereal* beta, doublereal* y, + const integer* incY, ftnlen trsize); +#else + + int _DGEMV_(const char* transpose, ftnlen trsize, + const integer* m, const integer* n, const doublereal* alpha, + const doublereal* a, const integer* lda, const doublereal* x, + const integer* incX, const doublereal* beta, doublereal* y, + const integer* incY); +#endif + + int _DGETRF_(const integer* m, const integer* n, + doublereal* a, integer* lda, integer* ipiv, + integer* info); + +#ifdef LAPACK_FTN_STRING_LEN_AT_END + + int _DGETRS_(const char* transpose, const integer* n, + const integer* nrhs, doublereal* a, const integer* lda, + integer* ipiv, doublereal* b, const integer* ldb, + integer* info, ftnlen trsize); + +#else + + int _DGETRS_(const char* transpose, ftnlen trsize, const integer* n, + const integer* nrhs, const doublereal* a, const integer* lda, + integer* ipiv, doublereal* b, const integer* ldb, integer* info); + +#endif + + int _DGETRI_(const integer* n, doublereal* a, const integer* lda, + integer* ipiv, doublereal* work, integer* lwork, integer* info); + + int _DGELSS_(integer *m, integer *n, integer *nrhs, + doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * + s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, + integer *info); + + + int _DGBSV_(integer *n, integer *kl, integer *ku, integer *nrhs, + doublereal *a, integer *lda, integer *ipiv, doublereal *b, + integer *ldb, integer *info); + + int _DGBTRF_(integer* m, integer *n, integer *kl, integer *ku, + doublereal *a, integer *lda, integer *ipiv, integer *info); + +#ifdef LAPACK_FTN_STRING_LEN_AT_END + int _DGBTRS_(const char* trans, integer *n, integer *kl, integer *ku, + integer *nrhs, doublereal *a, integer *lda, integer *ipiv, + doublereal *b, integer *ldb, integer *info, ftnlen trsize); +#else + int _DGBTRS_(const char* trans, ftnlen trsize, + integer *n, integer *kl, integer *ku, + integer *nrhs, doublereal *a, integer *lda, integer *ipiv, + doublereal *b, integer *ldb, integer *info); +#endif + + int _DSCAL_(integer *n, doublereal *da, doublereal *dx, integer *incx); +void cblas_dscal(const int N, const double alpha, double *X, const int incX); + +} +//#endif + +namespace Cantera { + + inline void ct_dgemv(ctlapack::storage_t storage, + ctlapack::transpose_t trans, + int m, int n, doublereal alpha, const doublereal* a, int lda, + const doublereal* x, int incX, doublereal beta, + doublereal* y, int incY) + { +#ifdef USE_CBLAS + cblas_dgemv(cblasOrder[storage], cblasTrans[trans], m, n, alpha, + a, lda, x, incX, beta, y, incY); +#else + integer f_m = m, f_n = n, f_lda = lda, f_incX = incX, f_incY = incY; + doublereal f_alpha = alpha, f_beta = beta; + ftnlen trsize = 1; +#ifdef NO_FTN_STRING_LEN_AT_END + _DGEMV_(&no_yes[trans], &f_m, &f_n, &f_alpha, a, + &f_lda, x, &f_incX, &f_beta, y, &f_incY); +#else +#ifdef LAPACK_FTN_STRING_LEN_AT_END + _DGEMV_(&no_yes[trans], &f_m, &f_n, &f_alpha, a, + &f_lda, x, &f_incX, &f_beta, y, &f_incY, trsize); +#else + _DGEMV_(&no_yes[trans], trsize, &f_m, &f_n, &f_alpha, a, + &f_lda, x, &f_incX, &f_beta, y, &f_incY); +#endif +#endif +#endif + + } + + + inline void ct_dgbsv(int n, int kl, int ku, int nrhs, + doublereal* a, int lda, integer* ipiv, doublereal* b, int ldb, + int& info) { + integer f_n = n, f_kl = kl, f_ku = ku, f_nrhs = nrhs, f_lda = lda, + f_ldb = ldb, f_info = info; + _DGBSV_(&f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv, + b, &f_ldb, &f_info); + info = f_info; + } + + inline void ct_dgbtrf(int m, int n, int kl, int ku, + doublereal* a, int lda, integer* ipiv, int& info) { + integer f_m = m, f_n = n, f_kl = kl, f_ku = ku, + f_lda = lda, f_info = info; + _DGBTRF_(&f_m, &f_n, &f_kl, &f_ku, a, &f_lda, ipiv, &f_info); + info = f_info; + } + + inline void ct_dgbtrs(ctlapack::transpose_t trans, int n, + int kl, int ku, int nrhs, doublereal* a, int lda, + integer* ipiv, doublereal* b, int ldb, int& info) { + integer f_n = n, f_kl = kl, f_ku = ku, f_nrhs = nrhs, f_lda = lda, + f_ldb = ldb, f_info = info; + char tr = no_yes[trans]; +#ifdef NO_FTN_STRING_LEN_AT_END + _DGBTRS_(&tr, &f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv, + b, &f_ldb, &f_info); +#else + ftnlen trsize = 1; +#ifdef LAPACK_FTN_STRING_LEN_AT_END + _DGBTRS_(&tr, &f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv, + b, &f_ldb, &f_info, trsize); +#else + _DGBTRS_(&tr, trsize, &f_n, &f_kl, &f_ku, &f_nrhs, a, &f_lda, ipiv, + b, &f_ldb, &f_info); +#endif +#endif + info = f_info; + } + + inline void ct_dgetrf(int m, int n, + doublereal* a, int lda, integer* ipiv, int& info) { + integer mm = m; + integer nn = n; + integer ldaa = lda; + integer infoo = info; + _DGETRF_(&mm, &nn, a, &ldaa, ipiv, &infoo); + info = infoo; + } + + inline void ct_dgetrs(ctlapack::transpose_t trans, int n, + int nrhs, doublereal* a, int lda, + integer* ipiv, doublereal* b, int ldb, int& info) + { + integer f_n = n, f_lda = lda, f_nrhs = nrhs, f_ldb = ldb, + f_info = info; + char tr = no_yes[trans]; + +#ifdef NO_FTN_STRING_LEN_AT_END + _DGETRS_(&tr, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb, + &f_info); +#else + ftnlen trsize = 1; +#ifdef LAPACK_FTN_STRING_LEN_AT_END + _DGETRS_(&tr, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb, + &f_info, trsize); +#else + _DGETRS_(&tr, trsize, &f_n, &f_nrhs, a, &f_lda, ipiv, b, &f_ldb, + &f_info); +#endif +#endif + info = f_info; + } + + inline void ct_dgetri(int n, doublereal* a, int lda, integer* ipiv, + doublereal* work, int lwork, int& info) { + integer f_n = n, f_lda = lda, f_lwork = lwork, f_info = info; + _DGETRI_(&f_n, a, &f_lda, ipiv, work, &f_lwork, &f_info); + } + + inline void ct_dgelss(int m, int n, int nrhs, doublereal* a, + int lda, doublereal* b, int ldb, doublereal* s, + doublereal rcond, int& rank, doublereal* work, int lwork, + int& info) { + doublereal f_rcond = rcond; + integer f_m = m, f_n = n, f_nrhs = nrhs, f_lda = lda, f_ldb = ldb, + f_rank = rank, f_info = info, f_lwork = lwork; + //f_lwork = 2*(3*min(m,n) + max(2*min(m,n), max(m,n))); + _DGELSS_(&f_m, &f_n, &f_nrhs, a, &f_lda, b, &f_ldb, s, &f_rcond, + &f_rank, work, &f_lwork, &f_info); + info = f_info; + rank = f_rank; + } + + inline void ct_dscal(int n, doublereal da, doublereal* dx, int incx) { + //integer f_n = n, f_incx = incx; + //_DSCAL_(&f_n, &da, dx, &f_incx); + cblas_dscal(n, da, dx, incx); + + } + + +} + +#endif diff --git a/Cantera/src/numerics/funcs.cpp b/Cantera/src/numerics/funcs.cpp new file mode 100755 index 000000000..0aaa08160 --- /dev/null +++ b/Cantera/src/numerics/funcs.cpp @@ -0,0 +1,79 @@ + +// miscellaneous functions + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include +#include + +using namespace std; + +#include "ct_defs.h" +#include "ctexceptions.h" +#include "stringUtils.h" + +extern "C" { + + int dpolft_(integer* n, doublereal* x, doublereal* y, doublereal* w, + integer* maxdeg, integer* ndeg, doublereal* eps, doublereal* r, + integer* ierr, doublereal* a); + + int dpcoef_(integer* l, doublereal* c, doublereal* tc, doublereal* a); +} + +namespace Cantera { + + /** + * Linearly interpolate a function defined on a discrete grid. + * vector xpts contains a monotonic sequence of grid points, and + * vector fpts contains function values defined at these points. + * The value returned is the linear interpolate at point x. + * If x is outside the range of xpts, the value of fpts at the + * nearest end is returned. + */ + + doublereal linearInterp(doublereal x, const vector_fp& xpts, + const vector_fp& fpts) { + if (x <= xpts[0]) + return fpts[0]; + if (x >= xpts.back()) + return fpts.back(); + vector_fp::const_iterator loc = + lower_bound(xpts.begin(), xpts.end(), x); + int iloc = int(loc - xpts.begin()) - 1; + doublereal ff = fpts[iloc] + + (x - xpts[iloc])*(fpts[iloc + 1] + - fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]); + return ff; + } + + + + doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w, + int maxdeg, int& ndeg, doublereal eps, doublereal* r) { + integer nn = n; + integer mdeg = maxdeg; + integer ndg = ndeg; + doublereal epss = eps; + integer ierr; + int worksize = 3*n + 3*maxdeg + 3; + vector_fp awork(worksize,0.0); + vector_fp coeffs(n+1, 0.0); + doublereal zer = 0.0; + + dpolft_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0], + &ierr, &awork[0]); + if (ierr != 1) throw CanteraError("polyfit", + "DPOLFT returned error code IERR = " + int2str(ierr) + + "while attempting to fit " + int2str(n) + " data points " + + "to a polynomial of degree " + int2str(maxdeg)); + ndeg = ndg; + dpcoef_(&ndg, &zer, r, &awork[0]); + return epss; + } + + +} diff --git a/Cantera/src/numerics/funcs.h b/Cantera/src/numerics/funcs.h new file mode 100644 index 000000000..d7b98e161 --- /dev/null +++ b/Cantera/src/numerics/funcs.h @@ -0,0 +1,11 @@ +#ifndef CT_FUNCS_H +#define CT_FUNCS_H + +#include "ct_defs.h" + +namespace Cantera { + doublereal linearInterp(doublereal x, const vector_fp& xpts, + const vector_fp& fpts); +} + +#endif diff --git a/Cantera/src/numerics/lapack.h b/Cantera/src/numerics/lapack.h new file mode 100755 index 000000000..f4ad15a4c --- /dev/null +++ b/Cantera/src/numerics/lapack.h @@ -0,0 +1,62 @@ + +// Copyright 2001 California Institute of Technology. +// All rights reserved. + +#include "ct_defs.h" + +#ifndef LAPACK_H +#define LAPACK_H + +#if defined(NEEDS_F77_TRANSLATION) + +#if defined(F77EXTERNS_UPPERCASE_NOTRAILINGBAR) +#define dgelss_ DGELSS +#define dgetrs_ DGETRS +#define dgetrf_ DGETRF +#define dgetri_ DGETRI +#define dvode_ DVODE +#define ddassl_ DDASSL +#define simplx_ SIMPLX +#define splin2_ SPLIN2 +#define splie2_ SPLIE2 +#define dgelss_ DGELSS +#endif + +#endif + +extern "C" +{ + + // /* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, + // doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * + // s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, + // integer *info); + + ///* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs, +// doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer * +// ldb, integer *info, int len); + +///* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer * +// lda, integer *ipiv, integer *info); + +///* Subroutine */ int dgetri_(integer *n, doublereal *a, integer * +// lda, integer *ipiv, doublereal *work, integer *lwork, integer *info); + +// /* Subroutine */ int dgemv_(char *trans, integer *m, integer *n, doublereal * +// alpha, doublereal *a, integer *lda, doublereal *x, integer *incx, +// doublereal *beta, doublereal *y, integer *incy, int len); + +void dcopy_(const integer *n, const doublereal *dx, const integer *incx, doublereal *dy, const integer *incy); + //doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy); + +doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy); + +void daxpy_(integer* n, doublereal* a, doublereal* x, integer* incx, + doublereal* y, integer* incy); +void dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx); +integer idamax_(integer* n, doublereal* a, integer* incx); + +} + + +#endif diff --git a/Cantera/src/numerics/polyfit.h b/Cantera/src/numerics/polyfit.h new file mode 100755 index 000000000..916f87929 --- /dev/null +++ b/Cantera/src/numerics/polyfit.h @@ -0,0 +1,20 @@ + +/* C interface for Fortran DPOLFT subroutine */ + +#ifndef CT_POLYFIT_H +#define CT_POLYFIT_H + +#include +//using namespace std; + +#include "ct_defs.h" + +namespace Cantera { + + doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w, + int maxdeg, int& ndeg, doublereal eps, doublereal* r); + +} +#endif + + diff --git a/Cantera/src/numerics/sort.cpp b/Cantera/src/numerics/sort.cpp new file mode 100755 index 000000000..ed6fc21c4 --- /dev/null +++ b/Cantera/src/numerics/sort.cpp @@ -0,0 +1,118 @@ +/** + * @file sort.cpp + * + * $Id$ + */ + +#ifdef WIN32 +#pragma warning(disable:4786) +#endif + +#include "sort.h" + +namespace Cantera { + + // sort (x,y) pairs by x + + void heapsort(vector_fp& x, vector_int& y) { + int n = x.size(); + if (n < 2) return; + doublereal rra; + integer rrb; + int ll = n/2; + int iret = n-1; + + while (1 > 0) { + if (ll > 0) { + ll--; + rra = x[ll]; + rrb = y[ll]; + } + else { + rra = x[iret]; + rrb = y[iret]; + x[iret] = x[0]; + y[iret] = y[0]; + iret--; + if (iret == 0) { + x[0] = rra; + y[0] = rrb; + return; + } + } + + int i = ll; + int j = ll + ll + 1; + + while (j <= iret) { + if (j < iret) { + if (x[j] < x[j+1]) + j++; + } + if (rra < x[j]) { + x[i] = x[j]; + y[i] = y[j]; + i = j; + j = j + j + 1; + } + else { + j = iret + 1; + } + } + x[i] = rra; + y[i] = rrb; + } + } + + void heapsort(vector_fp& x, vector_fp& y) { + int n = x.size(); + if (n < 2) return; + doublereal rra; + doublereal rrb; + int ll = n/2; + int iret = n-1; + + while (1 > 0) { + if (ll > 0) { + ll--; + rra = x[ll]; + rrb = y[ll]; + } + else { + rra = x[iret]; + rrb = y[iret]; + x[iret] = x[0]; + y[iret] = y[0]; + iret--; + if (iret == 0) { + x[0] = rra; + y[0] = rrb; + return; + } + } + + int i = ll; + int j = ll + ll + 1; + + while (j <= iret) { + if (j < iret) { + if (x[j] < x[j+1]) + j++; + } + if (rra < x[j]) { + x[i] = x[j]; + y[i] = y[j]; + i = j; + j = j + j + 1; + } + else { + j = iret + 1; + } + } + x[i] = rra; + y[i] = rrb; + } + } + +} + diff --git a/Cantera/src/numerics/sort.h b/Cantera/src/numerics/sort.h new file mode 100755 index 000000000..e83ab53b0 --- /dev/null +++ b/Cantera/src/numerics/sort.h @@ -0,0 +1,24 @@ +/** + * @file sort.h + * + * $Id$ + */ + +#ifndef CT_SORT_H +#define CT_SORT_H + +#include "ct_defs.h" + +namespace Cantera { + + /// Given two arrays x and y, sort the (x,y) pairs by the x + /// values. This version is for floating-point x, and integer y. + void heapsort(vector_fp& x, vector_int& y); + + /// Given two arrays x and y, sort the (x,y) pairs by the x + /// values. This version is for floating-point x, and + /// floating-point y. + void heapsort(vector_fp& x, vector_fp& y); +} + +#endif