oved files to numerics subdirectory
This commit is contained in:
Dave Goodwin
parent
0762b467c2
commit
e0f96378ff
28 changed files with 4762 additions and 0 deletions
153
Cantera/src/numerics/ArrayViewer.h
Executable file
153
Cantera/src/numerics/ArrayViewer.h
Executable file
|
|
@ -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<int>(m.nRows());
|
||||
int nc = static_cast<int>(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
|
||||
|
||||
|
||||
|
||||
181
Cantera/src/numerics/BandMatrix.cpp
Executable file
181
Cantera/src/numerics/BandMatrix.cpp
Executable file
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
153
Cantera/src/numerics/BandMatrix.h
Executable file
153
Cantera/src/numerics/BandMatrix.h
Executable file
|
|
@ -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
|
||||
310
Cantera/src/numerics/CVode.cpp
Executable file
310
Cantera/src/numerics/CVode.cpp
Executable file
|
|
@ -0,0 +1,310 @@
|
|||
/**
|
||||
* @file CVode.cpp
|
||||
*
|
||||
* $Id$
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
|
||||
#include "CVode.h"
|
||||
#include <iostream>
|
||||
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<N_Vector>(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<void*>(N_VNew(n, 0));
|
||||
}
|
||||
for (int i=0; i<n; i++) {
|
||||
N_VIth(nv(m_abstol), i) = abstol[i];
|
||||
}
|
||||
m_reltol = reltol;
|
||||
}
|
||||
|
||||
void CVodeInt::setTolerances(double reltol, double abstol) {
|
||||
m_itol = 0;
|
||||
m_reltol = reltol;
|
||||
m_abstols = abstol;
|
||||
}
|
||||
|
||||
void CVodeInt::setProblemType(int probtype) {
|
||||
m_type = probtype;
|
||||
}
|
||||
|
||||
void CVodeInt::setMethod(MethodType t) {
|
||||
if (t == BDF_Method)
|
||||
m_method = BDF;
|
||||
else if (t == Adams_Method)
|
||||
m_method = ADAMS;
|
||||
else
|
||||
throw CVodeErr("unknown method");
|
||||
}
|
||||
|
||||
void CVodeInt::setMaxStepSize(doublereal hmax) {
|
||||
m_hmax = hmax;
|
||||
m_ropt[HMAX] = hmax;
|
||||
}
|
||||
|
||||
void CVodeInt::setMinStepSize(doublereal hmin) {
|
||||
m_hmin = hmin;
|
||||
m_ropt[HMIN] = hmin;
|
||||
}
|
||||
|
||||
void CVodeInt::setMaxSteps(int nmax) {
|
||||
m_maxsteps = nmax;
|
||||
m_iopt[MXSTEP] = m_maxsteps;
|
||||
}
|
||||
|
||||
void CVodeInt::setIterator(IterType t) {
|
||||
if (t == Newton_Iter)
|
||||
m_iter = NEWTON;
|
||||
else if (t == Functional_Iter)
|
||||
m_iter = FUNCTIONAL;
|
||||
else
|
||||
throw CVodeErr("unknown iterator");
|
||||
}
|
||||
|
||||
void CVodeInt::initialize(double t0, FuncEval& func)
|
||||
{
|
||||
m_neq = func.neq();
|
||||
m_t0 = t0;
|
||||
|
||||
if (m_y) {
|
||||
N_VFree(nv(m_y)); // free solution vector if already allocated
|
||||
}
|
||||
m_y = reinterpret_cast<void*>(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]; }
|
||||
}
|
||||
|
||||
|
||||
|
||||
93
Cantera/src/numerics/CVode.h
Executable file
93
Cantera/src/numerics/CVode.h
Executable file
|
|
@ -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
|
||||
424
Cantera/src/numerics/CVodesIntegrator.cpp
Normal file
424
Cantera/src/numerics/CVodesIntegrator.cpp
Normal file
|
|
@ -0,0 +1,424 @@
|
|||
/**
|
||||
* @file CVodesIntegrator.cpp
|
||||
*
|
||||
*/
|
||||
|
||||
// Copyright 2001 California Institute of Technology
|
||||
#include "config.h"
|
||||
|
||||
#include "CVodesIntegrator.h"
|
||||
#include "stringUtils.h"
|
||||
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#ifdef SUNDIALS_VERSION_22
|
||||
|
||||
#include <sundials_types.h>
|
||||
#include <sundials_math.h>
|
||||
#include <cvodes.h>
|
||||
#include <cvodes_dense.h>
|
||||
#include <cvodes_diag.h>
|
||||
#include <cvodes_spgmr.h>
|
||||
#include <cvodes_band.h>
|
||||
#include <nvector_serial.h>
|
||||
|
||||
#else
|
||||
|
||||
#ifdef SUNDIALS_VERSION_23
|
||||
#include <sundials/sundials_types.h>
|
||||
#include <sundials/sundials_math.h>
|
||||
#include <sundials/sundials_nvector.h>
|
||||
#include <nvector/nvector_serial.h>
|
||||
#include <cvodes/cvodes.h>
|
||||
#include <cvodes/cvodes_dense.h>
|
||||
#include <cvodes/cvodes_diag.h>
|
||||
#include <cvodes/cvodes_spgmr.h>
|
||||
#include <cvodes/cvodes_band.h>
|
||||
|
||||
#else
|
||||
|
||||
unsupported sundials version!
|
||||
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
inline static N_Vector nv(void* x) {
|
||||
return reinterpret_cast<N_Vector>(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<void*>(N_VNew_Serial(n));
|
||||
}
|
||||
for (int i=0; i<n; i++) {
|
||||
NV_Ith_S(nv(m_abstol), i) = abstol[i];
|
||||
}
|
||||
m_reltol = reltol;
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setTolerances(double reltol, double abstol) {
|
||||
m_itol = CV_SS;
|
||||
m_reltol = reltol;
|
||||
m_abstols = abstol;
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) {
|
||||
m_reltolsens = reltol;
|
||||
m_abstolsens = abstol;
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setProblemType(int probtype) {
|
||||
m_type = probtype;
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setMethod(MethodType t) {
|
||||
if (t == BDF_Method)
|
||||
m_method = CV_BDF;
|
||||
else if (t == Adams_Method)
|
||||
m_method = CV_ADAMS;
|
||||
else
|
||||
throw CVodesErr("unknown method");
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setMaxStepSize(doublereal hmax) {
|
||||
m_hmax = hmax;
|
||||
if (m_cvode_mem)
|
||||
CVodeSetMaxStep(m_cvode_mem, hmax);
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setMinStepSize(doublereal hmin) {
|
||||
m_hmin = hmin;
|
||||
if (m_cvode_mem)
|
||||
CVodeSetMinStep(m_cvode_mem, hmin);
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setMaxSteps(int nmax) {
|
||||
m_maxsteps = nmax;
|
||||
if (m_cvode_mem)
|
||||
CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setIterator(IterType t) {
|
||||
if (t == Newton_Iter)
|
||||
m_iter = CV_NEWTON;
|
||||
else if (t == Functional_Iter)
|
||||
m_iter = CV_FUNCTIONAL;
|
||||
else
|
||||
throw CVodesErr("unknown iterator");
|
||||
}
|
||||
|
||||
void CVodesIntegrator::sensInit(double t0, FuncEval& func) {
|
||||
m_np = func.nparams();
|
||||
long int nv = func.neq();
|
||||
|
||||
doublereal* data;
|
||||
int n, j;
|
||||
N_Vector y;
|
||||
y = N_VNew_Serial(nv);
|
||||
m_yS = N_VCloneVectorArray_Serial(m_np, y);
|
||||
for (n = 0; n < m_np; n++) {
|
||||
data = NV_DATA_S(m_yS[n]);
|
||||
for (j = 0; j < nv; j++) {
|
||||
data[j] =0.0;
|
||||
}
|
||||
}
|
||||
|
||||
int flag;
|
||||
flag = CVodeSensMalloc(m_cvode_mem, m_np, CV_STAGGERED, m_yS);
|
||||
if (flag != CV_SUCCESS)
|
||||
throw CVodesErr("Error in CVodeSensMalloc");
|
||||
vector_fp atol(m_np, m_abstolsens);
|
||||
double rtol = m_reltolsens;
|
||||
flag = CVodeSetSensTolerances(m_cvode_mem, CV_SS, rtol, DATA_PTR(atol));
|
||||
}
|
||||
|
||||
void CVodesIntegrator::initialize(double t0, FuncEval& func)
|
||||
{
|
||||
m_neq = func.neq();
|
||||
m_t0 = t0;
|
||||
|
||||
if (m_y) {
|
||||
N_VDestroy_Serial(nv(m_y)); // free solution vector if already allocated
|
||||
}
|
||||
m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq)); // allocate solution vector
|
||||
for (int i=0; i<m_neq; i++) {
|
||||
NV_Ith_S(nv(m_y), i) = 0.0;
|
||||
}
|
||||
// check abs tolerance array size
|
||||
if (m_itol == CV_SV && m_nabs < m_neq)
|
||||
throw CVodesErr("not enough absolute tolerance values specified.");
|
||||
|
||||
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
|
||||
|
||||
if (m_cvode_mem) CVodeFree(&m_cvode_mem);
|
||||
m_cvode_mem = CVodeCreate(m_method, m_iter);
|
||||
if (!m_cvode_mem) throw CVodesErr("CVodeCreate failed.");
|
||||
|
||||
int flag = 0;
|
||||
if (m_itol == CV_SV) {
|
||||
// vector atol
|
||||
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
|
||||
m_reltol, nv(m_abstol));
|
||||
}
|
||||
else {
|
||||
// scalar atol
|
||||
flag = CVodeMalloc(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y), m_itol,
|
||||
m_reltol, &m_abstols);
|
||||
}
|
||||
if (flag != CV_SUCCESS) {
|
||||
if (flag == CV_MEM_FAIL) {
|
||||
throw CVodesErr("Memory allocation failed.");
|
||||
}
|
||||
else if (flag == CV_ILL_INPUT) {
|
||||
throw CVodesErr("Illegal value for CVodeMalloc input argument.");
|
||||
}
|
||||
else
|
||||
throw CVodesErr("CVodeMalloc failed.");
|
||||
}
|
||||
|
||||
|
||||
if (m_type == DENSE + NOJAC) {
|
||||
long int N = m_neq;
|
||||
CVDense(m_cvode_mem, N);
|
||||
}
|
||||
else if (m_type == DIAG) {
|
||||
CVDiag(m_cvode_mem);
|
||||
}
|
||||
else if (m_type == GMRES) {
|
||||
CVSpgmr(m_cvode_mem, PREC_NONE, 0);
|
||||
}
|
||||
else 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 {
|
||||
throw CVodesErr("unsupported option");
|
||||
}
|
||||
|
||||
// pass a pointer to func in m_data
|
||||
m_fdata = new FuncData(&func, func.nparams());
|
||||
|
||||
//m_data = (void*)&func;
|
||||
|
||||
flag = CVodeSetFdata(m_cvode_mem, (void*)m_fdata);
|
||||
if (flag != CV_SUCCESS)
|
||||
throw CVodesErr("CVodeSetFdata failed.");
|
||||
|
||||
if (func.nparams() > 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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
118
Cantera/src/numerics/CVodesIntegrator.h
Normal file
118
Cantera/src/numerics/CVodesIntegrator.h
Normal file
|
|
@ -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 <nvector_serial.h>
|
||||
#else
|
||||
#include <sundials/sundials_nvector.h>
|
||||
#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
|
||||
217
Cantera/src/numerics/DAE_Solver.h
Normal file
217
Cantera/src/numerics/DAE_Solver.h
Normal file
|
|
@ -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 <vector>
|
||||
|
||||
#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
|
||||
30
Cantera/src/numerics/DAE_solvers.cpp
Normal file
30
Cantera/src/numerics/DAE_solvers.cpp
Normal file
|
|
@ -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
|
||||
282
Cantera/src/numerics/DASPK.cpp
Executable file
282
Cantera/src/numerics/DASPK.cpp
Executable file
|
|
@ -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 <iostream>
|
||||
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<void **>(&(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<void *>(&m_resid);
|
||||
void **iddr_ipar = reinterpret_cast<void **>(&(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<int>& 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); }
|
||||
}
|
||||
|
||||
|
||||
|
||||
116
Cantera/src/numerics/DASPK.h
Executable file
116
Cantera/src/numerics/DASPK.h
Executable file
|
|
@ -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 <vector>
|
||||
|
||||
#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<int>& 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
|
||||
|
||||
148
Cantera/src/numerics/DenseMatrix.cpp
Executable file
148
Cantera/src/numerics/DenseMatrix.cpp
Executable file
|
|
@ -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<int>(nRows()),
|
||||
static_cast<int>(nRows()), 1.0, ptrColumn(0), //begin(),
|
||||
static_cast<int>(nRows()), b, 1, 0.0, prod, 1);
|
||||
}
|
||||
|
||||
void DenseMatrix::leftMult(const double* b, double* prod) const {
|
||||
int nc = static_cast<int>(nColumns());
|
||||
int nr = static_cast<int>(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<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()), &A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRF returned INFO = "+int2str(info));
|
||||
ct_dgetrs(ctlapack::NoTranspose,
|
||||
static_cast<int>(A.nRows()), 1, A.ptrColumn(0), //begin(),
|
||||
static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b,
|
||||
static_cast<int>(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<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), A.ptrColumn(0),
|
||||
static_cast<int>(A.nRows()), &A.ipiv()[0], info);
|
||||
if (info != 0)
|
||||
throw CanteraError("DenseMatrix::solve",
|
||||
"DGETRF returned INFO = "+int2str(info));
|
||||
ct_dgetrs(ctlapack::NoTranspose, static_cast<int>(A.nRows()),
|
||||
static_cast<int>(b.nColumns()),
|
||||
A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0], b.ptrColumn(0),
|
||||
static_cast<int>(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<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns())));
|
||||
ct_dgelss(static_cast<int>(A.nRows()),
|
||||
static_cast<int>(A.nColumns()), 1, A.ptrColumn(0),
|
||||
static_cast<int>(A.nRows()), b,
|
||||
static_cast<int>(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<int>(A.nRows()), static_cast<int>(A.nColumns()), 1.0,
|
||||
A.ptrColumn(0), static_cast<int>(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<int>(A.nRows()), static_cast<int>(A.nRows()), 1.0,
|
||||
A.ptrColumn(0), static_cast<int>(A.nRows()), b, 1, 1.0, prod, 1);
|
||||
}
|
||||
|
||||
int invert(DenseMatrix& A, int nn) {
|
||||
integer n = (nn > 0 ? nn : static_cast<int>(A.nRows()));
|
||||
int info=0;
|
||||
ct_dgetrf(n, n, A.ptrColumn(0), static_cast<int>(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<int>(work.size());
|
||||
ct_dgetri(n, A.ptrColumn(0), static_cast<int>(A.nRows()),
|
||||
&A.ipiv()[0],
|
||||
&work[0], lwork, info);
|
||||
if (info != 0)
|
||||
throw CanteraError("invert",
|
||||
"DGETRI returned INFO="+int2str(info));
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
112
Cantera/src/numerics/DenseMatrix.h
Executable file
112
Cantera/src/numerics/DenseMatrix.h
Executable file
|
|
@ -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 <iostream>
|
||||
//#include <vector>
|
||||
//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
|
||||
|
||||
|
||||
|
||||
405
Cantera/src/numerics/Func1.cpp
Normal file
405
Cantera/src/numerics/Func1.cpp
Normal file
|
|
@ -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 "<unknown " + int2str(ID()) + ">("+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));
|
||||
}
|
||||
|
||||
}
|
||||
651
Cantera/src/numerics/Func1.h
Normal file
651
Cantera/src/numerics/Func1.h
Normal file
|
|
@ -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 <iostream>
|
||||
#include <string>
|
||||
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
|
||||
69
Cantera/src/numerics/FuncEval.h
Executable file
69
Cantera/src/numerics/FuncEval.h
Executable file
|
|
@ -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
|
||||
256
Cantera/src/numerics/IDA_Solver.cpp
Normal file
256
Cantera/src/numerics/IDA_Solver.cpp
Normal file
|
|
@ -0,0 +1,256 @@
|
|||
|
||||
/**
|
||||
* @file IDA_Solver.cpp
|
||||
*
|
||||
*/
|
||||
|
||||
// Copyright 2006 California Institute of Technology
|
||||
|
||||
#include "IDA_Solver.h"
|
||||
#include "stringUtils.h"
|
||||
|
||||
#include <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include <sundials_types.h>
|
||||
#include <sundials_math.h>
|
||||
#include <ida.h>
|
||||
#include <ida_dense.h>
|
||||
#include <ida_spgmr.h>
|
||||
#include <ida_band.h>
|
||||
#include <nvector_serial.h>
|
||||
|
||||
inline static N_Vector nv(void* x) {
|
||||
return reinterpret_cast<N_Vector>(x);
|
||||
}
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* A simple class to hold an array of parameter values and a pointer to
|
||||
* an instance of a subclass of ResidEval.
|
||||
*/
|
||||
class ResidData {
|
||||
|
||||
|
||||
|
||||
|
||||
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<void*>(N_VNew_Serial(m_neq));
|
||||
for (int i=0; i < m_neq; i++) {
|
||||
NV_Ith_S(nv(m_abstol), i) = abstol[i];
|
||||
}
|
||||
m_reltol = reltol;
|
||||
}
|
||||
|
||||
void IDA_Solver::setTolerances(double reltol, double abstol) {
|
||||
m_itol = IDA_SS;
|
||||
m_reltol = reltol;
|
||||
m_abstols = abstol;
|
||||
}
|
||||
|
||||
void IDA_Solver::setLinearSolverType(int solverType) {
|
||||
m_type = solverType;
|
||||
}
|
||||
|
||||
void IDA_Solver::init(double t0)
|
||||
{
|
||||
m_t0 = t0;
|
||||
|
||||
if (m_y) N_VDestroy_Serial(nv(m_y));
|
||||
if (m_ydot) N_VDestroy_Serial(nv(m_ydot));
|
||||
if (m_id) N_VDestroy_Serial(nv(m_id));
|
||||
if (m_constraints) N_VDestroy_Serial(nv(m_constraints));
|
||||
|
||||
m_y = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
|
||||
m_ydot = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
|
||||
m_constraints = reinterpret_cast<void*>(N_VNew_Serial(m_neq));
|
||||
|
||||
for (int i=0; i<m_neq; i++) {
|
||||
NV_Ith_S(nv(m_y), i) = 0.0;
|
||||
NV_Ith_S(nv(m_ydot), i) = 0.0;
|
||||
NV_Ith_S(nv(m_constraints), i) = 0.0;
|
||||
}
|
||||
|
||||
// get the initial conditions
|
||||
m_resid.getInitialConditions(m_t0, NV_DATA_S(nv(m_ydot)),
|
||||
NV_DATA_S(nv(m_y)));
|
||||
|
||||
if (m_ida_mem) IDAFree(&m_ida_mem);
|
||||
m_ida_mem = IDACreate();
|
||||
|
||||
int flag = 0;
|
||||
if (m_itol == IDA_SV) {
|
||||
// vector atol
|
||||
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
|
||||
m_itol, m_reltol, nv(m_abstol));
|
||||
}
|
||||
else {
|
||||
// scalar atol
|
||||
flag = IDAMalloc(m_ida_mem, ida_resid, m_t0, nv(m_y), nv(m_ydot),
|
||||
m_itol, m_reltol, &m_abstols);
|
||||
}
|
||||
if (flag != IDA_SUCCESS) {
|
||||
if (flag == IDA_MEM_FAIL) {
|
||||
throw IDA_Err("Memory allocation failed."); }
|
||||
else if (flag == IDA_ILL_INPUT) {
|
||||
throw IDA_Err("Illegal value for IDAMalloc input argument.");
|
||||
}
|
||||
else
|
||||
throw IDA_Err("IDAMalloc failed.");
|
||||
}
|
||||
|
||||
//-----------------------------------
|
||||
// set the linear solver type
|
||||
//-----------------------------------
|
||||
|
||||
if (m_type == 1) {
|
||||
long int N = m_neq;
|
||||
IDADense(m_ida_mem, N);
|
||||
}
|
||||
else if (m_type == 2) {
|
||||
long int N = m_neq;
|
||||
long int nu = m_mupper;
|
||||
long int nl = m_mlower;
|
||||
IDABand(m_ida_mem, N, nu, nl);
|
||||
}
|
||||
else {
|
||||
throw IDA_Err("unsupported linear solver type");
|
||||
}
|
||||
|
||||
|
||||
// pass a pointer to func in m_data
|
||||
m_fdata = new FuncData(&func, func.nparams());
|
||||
|
||||
flag = IDASetRdata(m_ida_mem, (void*)m_fdata);
|
||||
if (flag != IDA_SUCCESS)
|
||||
throw IDA_Err("IDASetRdata failed.");
|
||||
|
||||
// set options
|
||||
//if (m_maxord > 0)
|
||||
// flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
|
||||
//if (m_maxsteps > 0)
|
||||
// flag = CVodeSetMaxNumSteps(m_cvode_mem, m_maxsteps);
|
||||
//if (m_hmax > 0)
|
||||
// flag = CVodeSetMaxStep(m_cvode_mem, m_hmax);
|
||||
}
|
||||
|
||||
void IDA_Solver::solve(double tout)
|
||||
{
|
||||
double t;
|
||||
int flag;
|
||||
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_NORMAL);
|
||||
if (flag != IDA_SUCCESS)
|
||||
throw IDA_Err(" IDA error encountered.");
|
||||
}
|
||||
|
||||
double IDA_Solver::step(double tout)
|
||||
{
|
||||
double t;
|
||||
int flag;
|
||||
flag = IDASolve(m_ida_mem, tout, &t, nv(m_y), nv(m_ydot), IDA_ONE_STEP);
|
||||
if (flag != IDA_SUCCESS)
|
||||
throw IDA_Err(" IDA error encountered.");
|
||||
return t;
|
||||
}
|
||||
|
||||
doublereal IDA_Solver::getOutputParameter(int flag) {
|
||||
switch (flag) {
|
||||
case REAL_WORKSPACE_SIZE:
|
||||
flag = IDAGetWorkSpace(m_ida_mem, &lenrw, &leniw);
|
||||
return doublereal(lenrw);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
137
Cantera/src/numerics/IDA_Solver.h
Normal file
137
Cantera/src/numerics/IDA_Solver.h
Normal file
|
|
@ -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 <vector>
|
||||
|
||||
#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
|
||||
|
||||
188
Cantera/src/numerics/Integrator.h
Executable file
188
Cantera/src/numerics/Integrator.h
Executable file
|
|
@ -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
|
||||
29
Cantera/src/numerics/ODE_integrators.cpp
Normal file
29
Cantera/src/numerics/ODE_integrators.cpp
Normal file
|
|
@ -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);
|
||||
}
|
||||
}
|
||||
}
|
||||
100
Cantera/src/numerics/ResidEval.h
Executable file
100
Cantera/src/numerics/ResidEval.h
Executable file
|
|
@ -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<int, int> m_alg;
|
||||
std::map<int, int> m_constrain;
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
276
Cantera/src/numerics/ctlapack.h
Executable file
276
Cantera/src/numerics/ctlapack.h
Executable file
|
|
@ -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 <vecLib/cblas.h>
|
||||
|
||||
// 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 <Accelerate.h>
|
||||
const CBLAS_ORDER cblasOrder[2] = { CblasRowMajor, CblasColMajor };
|
||||
const CBLAS_TRANSPOSE cblasTrans[2] = { CblasNoTrans, CblasTrans };
|
||||
#endif
|
||||
|
||||
//#ifdef DARWIN
|
||||
//#include <Accelerate.h>
|
||||
//#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
|
||||
79
Cantera/src/numerics/funcs.cpp
Executable file
79
Cantera/src/numerics/funcs.cpp
Executable file
|
|
@ -0,0 +1,79 @@
|
|||
|
||||
// miscellaneous functions
|
||||
|
||||
#ifdef WIN32
|
||||
#pragma warning(disable:4786)
|
||||
#pragma warning(disable:4503)
|
||||
#endif
|
||||
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
11
Cantera/src/numerics/funcs.h
Normal file
11
Cantera/src/numerics/funcs.h
Normal file
|
|
@ -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
|
||||
62
Cantera/src/numerics/lapack.h
Executable file
62
Cantera/src/numerics/lapack.h
Executable file
|
|
@ -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
|
||||
20
Cantera/src/numerics/polyfit.h
Executable file
20
Cantera/src/numerics/polyfit.h
Executable file
|
|
@ -0,0 +1,20 @@
|
|||
|
||||
/* C interface for Fortran DPOLFT subroutine */
|
||||
|
||||
#ifndef CT_POLYFIT_H
|
||||
#define CT_POLYFIT_H
|
||||
|
||||
#include <vector>
|
||||
//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
|
||||
|
||||
|
||||
118
Cantera/src/numerics/sort.cpp
Executable file
118
Cantera/src/numerics/sort.cpp
Executable file
|
|
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
24
Cantera/src/numerics/sort.h
Executable file
24
Cantera/src/numerics/sort.h
Executable file
|
|
@ -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
|
||||
Loading…
Add table
Reference in a new issue