Merging changes from the trunk into the change branch.
This commit is contained in:
commit
fd245bb194
18 changed files with 762 additions and 461 deletions
|
|
@ -38,14 +38,27 @@ namespace Cantera {
|
|||
|
||||
public:
|
||||
|
||||
//! Type definition for the iterator class that is
|
||||
//! can be used by Array2D types.
|
||||
/*!
|
||||
* this is just equal to vector_fp iterator.
|
||||
*/
|
||||
typedef vector_fp::iterator iterator;
|
||||
|
||||
|
||||
//! Type definition for the const_iterator class that is
|
||||
//! can be used by Array2D types.
|
||||
/*!
|
||||
* this is just equal to vector_fp const_iterator.
|
||||
*/
|
||||
typedef vector_fp::const_iterator const_iterator;
|
||||
|
||||
/**
|
||||
* Default constructor. Create an empty array.
|
||||
*/
|
||||
Array2D() : m_nrows(0), m_ncols(0) { m_data.clear(); }
|
||||
|
||||
Array2D() : m_nrows(0), m_ncols(0) {
|
||||
m_data.clear();
|
||||
}
|
||||
|
||||
//! Constructor.
|
||||
/*!
|
||||
|
|
@ -86,7 +99,7 @@ namespace Cantera {
|
|||
return *this;
|
||||
}
|
||||
|
||||
//! resize the array, and fill the new entries with 'v'
|
||||
//! Resize the array, and fill the new entries with 'v'
|
||||
/*!
|
||||
* @param n This is the number of rows
|
||||
* @param m This is the number of columns in the new matrix
|
||||
|
|
@ -98,7 +111,14 @@ namespace Cantera {
|
|||
m_data.resize(n*m, v);
|
||||
}
|
||||
|
||||
/// append a column
|
||||
//! Append a column to the existing matrix using a std vector
|
||||
/*!
|
||||
* This operation will add a column onto the existing matrix.
|
||||
*
|
||||
* @param c This vector<doublereal> is the entries in the
|
||||
* column to be added. It must have a length
|
||||
* equal to m_nrows or greater.
|
||||
*/
|
||||
void appendColumn(const vector_fp& c) {
|
||||
m_ncols++;
|
||||
m_data.resize(m_nrows*m_ncols);
|
||||
|
|
@ -106,37 +126,66 @@ namespace Cantera {
|
|||
for (m = 0; m < m_nrows; m++) value(m_ncols, m) = c[m];
|
||||
}
|
||||
|
||||
/// append a column
|
||||
void appendColumn(doublereal* c) {
|
||||
//! Append a column to the existing matrix
|
||||
/*!
|
||||
* This operation will add a column onto the existing matrix.
|
||||
*
|
||||
* @param c This vector of doubles is the entries in the
|
||||
* column to be added. It must have a length
|
||||
* equal to m_nrows or greater.
|
||||
*/
|
||||
void appendColumn(const doublereal* const c) {
|
||||
m_ncols++;
|
||||
m_data.resize(m_nrows*m_ncols);
|
||||
int m;
|
||||
for (m = 0; m < m_nrows; m++) value(m_ncols, m) = c[m];
|
||||
}
|
||||
|
||||
/// set the nth row to array rw
|
||||
void setRow(int n, doublereal* rw) {
|
||||
//! Set the nth row to array rw
|
||||
/*!
|
||||
* @param n Index of the row to be changed
|
||||
* @param rw Vector for the row. Must have a length of m_ncols.
|
||||
*/
|
||||
void setRow(int n, const doublereal* const rw) {
|
||||
for (int j = 0; j < m_ncols; j++) {
|
||||
m_data[m_nrows*j + n] = rw[j];
|
||||
}
|
||||
}
|
||||
|
||||
/// get the nth row
|
||||
void getRow(int n, doublereal* rw) {
|
||||
//! Get the nth row and return it in a vector
|
||||
/*!
|
||||
* @param n Index of the row to be returned.
|
||||
* @param rw Return Vector for the operation.
|
||||
* Must have a length of m_ncols.
|
||||
*/
|
||||
void getRow(int n, doublereal* const rw) {
|
||||
for (int j = 0; j < m_ncols; j++) {
|
||||
rw[j] = m_data[m_nrows*j + n];
|
||||
}
|
||||
}
|
||||
|
||||
/// set the values in column m to those in array col
|
||||
void setColumn(int m, doublereal* col) {
|
||||
//! Set the values in column m to those in array col
|
||||
/*!
|
||||
* A(i,m) = col(i)
|
||||
*
|
||||
* @param m Column to set
|
||||
* @param col pointer to a col vector. Vector
|
||||
* must have a length of m_nrows.
|
||||
*/
|
||||
void setColumn(int m, doublereal* const col) {
|
||||
for (int i = 0; i < m_nrows; i++) {
|
||||
m_data[m_nrows*m + i] = col[i];
|
||||
}
|
||||
}
|
||||
|
||||
/// get the values in column m
|
||||
void getColumn(int m, doublereal* col) {
|
||||
//! Get the values in column m
|
||||
/*!
|
||||
* col(i) = A(i,m)
|
||||
*
|
||||
* @param m Column to set
|
||||
* @param col pointer to a col vector that will be returned
|
||||
*/
|
||||
void getColumn(int m, doublereal* const col) {
|
||||
for (int i = 0; i < m_nrows; i++) {
|
||||
col[i] = m_data[m_nrows*m + i];
|
||||
}
|
||||
|
|
@ -147,10 +196,18 @@ namespace Cantera {
|
|||
* heap.
|
||||
*/
|
||||
virtual ~Array2D(){}
|
||||
|
||||
|
||||
/**
|
||||
* Evaluate a*x + y.
|
||||
|
||||
//! Evaluate z = a*x + y.
|
||||
/*!
|
||||
* This function evaluates the AXPY operation, and stores
|
||||
* the result in the object's Array2D object.
|
||||
* It's assumed that all 3 objects have the same dimensions,
|
||||
* but no error checking is done.
|
||||
*
|
||||
* @param a scalar to multiply x with
|
||||
* @param x First Array2D object to be used
|
||||
* @param y Second Array2D object to be used
|
||||
*
|
||||
*/
|
||||
void axpy(doublereal a, const Array2D& x, const Array2D& y) {
|
||||
iterator b = begin();
|
||||
|
|
@ -169,20 +226,43 @@ namespace Cantera {
|
|||
*/
|
||||
doublereal& operator()( int i, int j) { return value(i,j); }
|
||||
|
||||
/**
|
||||
* Allows retrieving elements using the syntax x = A(i,j).
|
||||
|
||||
//! Allows retrieving elements using the syntax x = A(i,j).
|
||||
/*!
|
||||
* @param i Index for the row to be retrieved
|
||||
* @param j Index for the column to be retrieved.
|
||||
*
|
||||
* @return Returns the value of the matrix entry
|
||||
*/
|
||||
doublereal operator() ( int i, int j) const {return value(i,j);}
|
||||
doublereal operator() (int i, int j) const {
|
||||
return value(i,j);
|
||||
}
|
||||
|
||||
//! Returns a changeable reference to position in the matrix
|
||||
/*!
|
||||
* This is a key entry. Returns a reference to the matrixes (i,j)
|
||||
* element. This may be used as an L value.
|
||||
*
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*
|
||||
* @return Returns a changeable reference to the matrix entry
|
||||
*/
|
||||
doublereal& value(int i, int j) {
|
||||
return m_data[m_nrows*j + i];
|
||||
}
|
||||
|
||||
//! Returns the value of a single matrix entry
|
||||
/*!
|
||||
* This is a key entry. Returns the value of the matrix position (i,j)
|
||||
* element.
|
||||
*
|
||||
* @param i The row index
|
||||
* @param j The column index
|
||||
*/
|
||||
doublereal& value( int i, int j) {return m_data[m_nrows*j + i];}
|
||||
doublereal value( int i, int j) const {return m_data[m_nrows*j + i];}
|
||||
doublereal value(int i, int j) const {
|
||||
return m_data[m_nrows*j + i];
|
||||
}
|
||||
|
||||
/// Number of rows
|
||||
size_t nRows() const { return m_nrows; }
|
||||
|
|
@ -208,10 +288,25 @@ namespace Cantera {
|
|||
/// Return a const reference to the data vector
|
||||
const vector_fp& data() const { return m_data; }
|
||||
|
||||
/// Return a pointer to the top of column j, columns are contiguous
|
||||
/// in memory
|
||||
//! Return a pointer to the top of column j, columns are contiguous
|
||||
//! in memory
|
||||
/*!
|
||||
* @param j Value of the column
|
||||
*
|
||||
* @return Returns a pointer to the top of the column
|
||||
*/
|
||||
doublereal * ptrColumn(int j) { return &(m_data[m_nrows*j]); }
|
||||
const doublereal * ptrColumn(int j) const { return &(m_data[m_nrows*j]); }
|
||||
|
||||
//! Return a const pointer to the top of column j, columns are contiguous
|
||||
//! in memory
|
||||
/*!
|
||||
* @param j Value of the column
|
||||
*
|
||||
* @return Returns a const pointer to the top of the column
|
||||
*/
|
||||
const doublereal * ptrColumn(int j) const {
|
||||
return &(m_data[m_nrows*j]);
|
||||
}
|
||||
|
||||
protected:
|
||||
|
||||
|
|
@ -225,7 +320,16 @@ namespace Cantera {
|
|||
int m_ncols;
|
||||
};
|
||||
|
||||
/// output the array
|
||||
//! Output the current contents of the Array2D object
|
||||
/*!
|
||||
* Example of usage:
|
||||
* s << m << endl;
|
||||
*
|
||||
* @param s Reference to the ostream to write to
|
||||
* @param m Object of type Array2D that you are querying
|
||||
*
|
||||
* @return Returns a reference to the ostream.
|
||||
*/
|
||||
inline std::ostream& operator<<(std::ostream& s, const Array2D& m) {
|
||||
int nr = static_cast<int>(m.nRows());
|
||||
int nc = static_cast<int>(m.nColumns());
|
||||
|
|
|
|||
|
|
@ -1519,10 +1519,6 @@ protected:
|
|||
app()->writelog(msg);
|
||||
}
|
||||
|
||||
void writelogAM(const std::string& msg) {
|
||||
app()->writelog(msg);
|
||||
}
|
||||
|
||||
// Write a message to the screen
|
||||
void Application::Messages::writelog(const std::string& msg) {
|
||||
logwriter->write(msg);
|
||||
|
|
@ -1532,9 +1528,6 @@ protected:
|
|||
void writelog(const char* msg) {
|
||||
app()->writelog(msg);
|
||||
}
|
||||
void writelogAM(const char* msg) {
|
||||
app()->writelog(msg);
|
||||
}
|
||||
|
||||
// Write a message to the screen
|
||||
void Application::Messages::writelog(const char* pszmsg) {
|
||||
|
|
|
|||
|
|
@ -39,9 +39,13 @@ namespace Cantera {
|
|||
std::copy(x.begin(), x.begin() + n, y.begin());
|
||||
}
|
||||
|
||||
/**
|
||||
* Divide each element of x by the corresponding element of y.
|
||||
//! Divide each element of x by the corresponding element of y.
|
||||
/*!
|
||||
* This function replaces x[n] by x[n]/y[n], for 0 <= n < x.size()
|
||||
*
|
||||
* @param x Numerator object of the division operation with template type T
|
||||
* At the end of the calculation, it contains the result.
|
||||
* @param y Denominator object of the division template type T
|
||||
*/
|
||||
template<class T>
|
||||
inline void divide_each(T& x, const T& y) {
|
||||
|
|
@ -49,9 +53,15 @@ namespace Cantera {
|
|||
x.begin(), std::divides<TYPENAME_KEYWORD T::value_type>());
|
||||
}
|
||||
|
||||
/**
|
||||
* multiply each element of x by the corresponding element of y.
|
||||
//! Multiply each element of x by the corresponding element of y.
|
||||
/*!
|
||||
* This function replaces x[n] by x[n]*y[n], for 0 <= n < x.size()
|
||||
* This is a templated function with just one template type.
|
||||
*
|
||||
* @param x First object of the multiplication with template type T
|
||||
* At the end of the calculation, it contains the result.
|
||||
* @param y Second object of the multiplication with template type T
|
||||
*
|
||||
*/
|
||||
template<class T>
|
||||
inline void multiply_each(T& x, const T& y) {
|
||||
|
|
@ -59,32 +69,52 @@ namespace Cantera {
|
|||
x.begin(), std::multiplies<TYPENAME_KEYWORD T::value_type>());
|
||||
}
|
||||
|
||||
/**
|
||||
* Multiply each element of x by scale_factor.
|
||||
//! Multiply each element of x by scale_factor.
|
||||
/*!
|
||||
* This function replaces x[n] by x[n]*scale_factor, for 0 <= n < x.size()
|
||||
*
|
||||
* @param x First object of the multiplication with template type T
|
||||
* At the end of the calculation, it contains the result.
|
||||
* @param scale_factor scale factor with template type S
|
||||
*/
|
||||
template<class T, class S>
|
||||
inline void scale(T& x, S scale_factor) {
|
||||
scale(x.begin(), x.end(), x.begin(), scale_factor);
|
||||
}
|
||||
|
||||
/**
|
||||
//! Return the templated dot product of two objects
|
||||
/*!
|
||||
* Returns the sum of x[n]*y[n], for 0 <= n < x.size().
|
||||
*
|
||||
* @param x First object of the dot product with template type T
|
||||
* At the end of the calculation, it contains the result.
|
||||
* @param y Second object of the dot product with template type T
|
||||
*/
|
||||
template<class T>
|
||||
inline doublereal dot_product(const T& x, const T& y) {
|
||||
return std::inner_product(x.begin(), x.end(), y.begin(), 0.0);
|
||||
}
|
||||
|
||||
//! Returns the templated dot ratio of two objects
|
||||
/**
|
||||
* Returns the sum of x[n]/y[n], for 0 <= n < x.size().
|
||||
*
|
||||
* @param x First object of the dot product with template type T
|
||||
* At the end of the calculation, it contains the result.
|
||||
* @param y Second object of the dot product with template type T
|
||||
*/
|
||||
template<class T>
|
||||
inline doublereal dot_ratio(const T& x, const T& y) {
|
||||
return _dot_ratio(x.begin(), x.end(), y.begin(), 0.0);
|
||||
}
|
||||
|
||||
//! Returns a templated addition operation of two objects
|
||||
/**
|
||||
* Replaces x[n] by x[n] + y[n] for 0 <= n < x.size()
|
||||
*
|
||||
* @param x First object of the addition with template type T
|
||||
* At the end of the calculation, it contains the result.
|
||||
* @param y Second object of the addition with template type T
|
||||
*/
|
||||
template<class T>
|
||||
inline void add_each(T& x, const T& y) {
|
||||
|
|
@ -119,9 +149,13 @@ namespace Cantera {
|
|||
return start_value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Finds the entry in a vector with maximum absolute
|
||||
* value, and return this value.
|
||||
|
||||
//! Finds the entry in a vector with maximum absolute
|
||||
//! value, and return this value.
|
||||
/*!
|
||||
* @param v Vector to be queried for maximum value, with template type T
|
||||
*
|
||||
* @return Returns an object of type T that is the maximum value,
|
||||
*/
|
||||
template<class T>
|
||||
inline T absmax(const std::vector<T>& v) {
|
||||
|
|
|
|||
|
|
@ -496,7 +496,8 @@ namespace Cantera {
|
|||
|
||||
//! Adds moles of a certain species to the mixture
|
||||
/*!
|
||||
*
|
||||
* @param indexS Index of the species in the MultiPhase object
|
||||
* @param addedMoles Value of the moles that are added to the species.
|
||||
*/
|
||||
void addSpeciesMoles(const int indexS, const doublereal addedMoles);
|
||||
|
||||
|
|
|
|||
|
|
@ -15,6 +15,8 @@
|
|||
#ifndef _VCS_INTERNAL_H
|
||||
#define _VCS_INTERNAL_H
|
||||
|
||||
#include <cstring>
|
||||
|
||||
#include "vcs_defs.h"
|
||||
#include "vcs_DoubleStarStar.h"
|
||||
#include "vcs_Exception.h"
|
||||
|
|
@ -324,7 +326,6 @@ namespace VCSnonideal {
|
|||
//! available if this ever fails.
|
||||
#define USE_MEMSET
|
||||
#ifdef USE_MEMSET
|
||||
#include <cstring>
|
||||
|
||||
//! Zero a double vector
|
||||
/*!
|
||||
|
|
@ -473,6 +474,15 @@ namespace VCSnonideal {
|
|||
*/
|
||||
void vcs_print_line(const char *str, int num);
|
||||
|
||||
//! Returns a const char string representing the type of the
|
||||
//! species given by the first argument
|
||||
/*!
|
||||
* @param speciesStatus Species status integer representing the type
|
||||
* of the species.
|
||||
* @param length Maximum length of the string to be returned.
|
||||
* Shorter values will yield abbreviated strings.
|
||||
* Defaults to a value of 100.
|
||||
*/
|
||||
const char *vcs_speciesType_string(int speciesStatus, int length = 100);
|
||||
|
||||
//! Print a string within a given space limit
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@ using namespace std;
|
|||
|
||||
#else
|
||||
|
||||
#ifdef SUNDIALS_VERSION_23
|
||||
#if defined(SUNDIALS_VERSION_23) || defined (SUNDIALS_VERSION_24)
|
||||
#include <sundials/sundials_types.h>
|
||||
#include <sundials/sundials_math.h>
|
||||
#include <sundials/sundials_nvector.h>
|
||||
|
|
@ -42,383 +42,471 @@ unsupported sundials version!
|
|||
|
||||
#endif
|
||||
|
||||
#if defined (SUNDIALS_VERSION_24)
|
||||
#define CV_SS 1
|
||||
#define CV_SV 2
|
||||
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
||||
inline static N_Vector nv(void* x) {
|
||||
return reinterpret_cast<N_Vector>(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;
|
||||
};
|
||||
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;
|
||||
}
|
||||
/**
|
||||
* 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);
|
||||
/**
|
||||
* 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;
|
||||
|
||||
|
||||
/// 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;
|
||||
}
|
||||
//delete[] m_iopt;
|
||||
}
|
||||
|
||||
double& CVodesIntegrator::solution(int k){
|
||||
return NV_Ith_S(nv(m_y),k);
|
||||
}
|
||||
double& CVodesIntegrator::solution(int k){
|
||||
return NV_Ith_S(nv(m_y),k);
|
||||
}
|
||||
|
||||
double* CVodesIntegrator::solution(){ return NV_DATA_S(nv(m_y));
|
||||
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));
|
||||
}
|
||||
|
||||
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;
|
||||
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::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::setSensitivityTolerances(double reltol, double abstol) {
|
||||
m_reltolsens = reltol;
|
||||
m_abstolsens = abstol;
|
||||
}
|
||||
|
||||
void CVodesIntegrator::setProblemType(int probtype) {
|
||||
m_type = probtype;
|
||||
}
|
||||
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::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::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::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::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::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();
|
||||
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;
|
||||
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 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);
|
||||
int flag;
|
||||
|
||||
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
|
||||
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));
|
||||
#elif defined(SUNDIALS_VERSION_24)
|
||||
flag = CVodeSensInit(m_cvode_mem, m_np, CV_STAGGERED,
|
||||
CVSensRhsFn (0), m_yS);
|
||||
|
||||
if (flag != CV_SUCCESS) {
|
||||
throw CVodesErr("Error in CVodeSensMalloc");
|
||||
}
|
||||
vector_fp atol(m_np, m_abstolsens);
|
||||
double rtol = m_reltolsens;
|
||||
flag = CVodeSensSStolerances(m_cvode_mem, rtol, DATA_PTR(atol));
|
||||
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
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);
|
||||
|
||||
/*
|
||||
* Specify the method and the iteration type:
|
||||
* Cantera Defaults:
|
||||
* CV_BDF - Use BDF methods
|
||||
* CV_NEWTON - use newton's method
|
||||
*/
|
||||
m_cvode_mem = CVodeCreate(m_method, m_iter);
|
||||
if (!m_cvode_mem) throw CVodesErr("CVodeCreate failed.");
|
||||
|
||||
int flag = 0;
|
||||
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
|
||||
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.");
|
||||
}
|
||||
#elif defined(SUNDIALS_VERSION_24)
|
||||
|
||||
flag = CVodeInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y));
|
||||
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 CVodeInit input argument.");
|
||||
} else {
|
||||
throw CVodesErr("CVodeInit failed.");
|
||||
}
|
||||
}
|
||||
|
||||
if (m_itol == CV_SV) {
|
||||
flag = CVodeSVtolerances(m_cvode_mem, m_reltol, nv(m_abstol));
|
||||
} else {
|
||||
flag = CVodeSStolerances(m_cvode_mem, 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 CVodeInit input argument.");
|
||||
} else {
|
||||
throw CVodesErr("CVodeInit failed.");
|
||||
}
|
||||
}
|
||||
#else
|
||||
printf("unknown sundials verson\n");
|
||||
exit(-1);
|
||||
#endif
|
||||
|
||||
|
||||
|
||||
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;
|
||||
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
|
||||
flag = CVodeSetFdata(m_cvode_mem, (void*)m_fdata);
|
||||
if (flag != CV_SUCCESS) {
|
||||
throw CVodesErr("CVodeSetFdata failed.");
|
||||
}
|
||||
#elif defined(SUNDIALS_VERSION_24)
|
||||
flag = CVodeSetUserData(m_cvode_mem, (void*)m_fdata);
|
||||
if (flag != CV_SUCCESS) {
|
||||
throw CVodesErr("CVodeSetUserData failed.");
|
||||
}
|
||||
#endif
|
||||
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 defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
|
||||
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("CVodeReInit failed. result = "+int2str(result));
|
||||
}
|
||||
#elif defined(SUNDIALS_VERSION_24)
|
||||
result = CVodeReInit(m_cvode_mem, m_t0, nv(m_y));
|
||||
if (result != CV_SUCCESS) {
|
||||
throw CVodesErr("CVodeReInit failed. result = "+int2str(result));
|
||||
}
|
||||
#endif
|
||||
|
||||
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;
|
||||
double tretn;
|
||||
flag = CVode(m_cvode_mem, tout, nv(m_y), &t, CV_NORMAL);
|
||||
if (flag != CV_SUCCESS)
|
||||
throw CVodesErr(" CVodes error encountered.");
|
||||
#if defined(SUNDIALS_VERSION_22) || defined(SUNDIALS_VERSION23)
|
||||
if (m_np > 0) {
|
||||
CVodeGetSens(m_cvode_mem, tout, m_yS);
|
||||
}
|
||||
#elif defined(SUNDIALS_VERSION_24)
|
||||
if (m_np > 0) {
|
||||
CVodeGetSens(m_cvode_mem, &tretn, m_yS);
|
||||
if (fabs(tretn - tout) > 1.0E-5) {
|
||||
throw CVodesErr("Time of Sensitivities different than time of tout");
|
||||
}
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -53,6 +53,16 @@ namespace Cantera {
|
|||
|
||||
};
|
||||
|
||||
//! Prints out the current internal state of the Crystal ThermoPhase object
|
||||
/*!
|
||||
* Example of usage:
|
||||
* s << x << endl;
|
||||
*
|
||||
* @param s Reference to the ostream to write to
|
||||
* @param x Object of type Crystal that you are querying
|
||||
*
|
||||
* @return Returns a reference to the ostream.
|
||||
*/
|
||||
inline std::ostream& operator<<(std::ostream& s, Cantera::Crystal& x) {
|
||||
size_t ip;
|
||||
for (ip = 0; ip < x.nPhases(); ip++) {
|
||||
|
|
|
|||
|
|
@ -551,7 +551,6 @@ namespace Cantera {
|
|||
* It can be shown that the expression
|
||||
*
|
||||
*
|
||||
*
|
||||
* \f[
|
||||
* B^{\phi}_{ca} = \beta^{(0)}_{ca} + \beta^{(1)}_{ca} \exp{(- \alpha^{(1)}_{ca} \sqrt{I})}
|
||||
* + \beta^{(2)}_{ca} \exp{(- \alpha^{(2)}_{ca} \sqrt{I} )}
|
||||
|
|
@ -3201,6 +3200,7 @@ namespace Cantera {
|
|||
|
||||
//! gamma_o value for the cutoff process at the zero solvent point
|
||||
doublereal MC_X_o_min_;
|
||||
|
||||
//! Parameter in the Molality Exp cutoff treatment
|
||||
/*!
|
||||
* This is the slope of the p function at the zero solvent point
|
||||
|
|
@ -3223,10 +3223,16 @@ namespace Cantera {
|
|||
//! Parameter in the Molality Exp cutoff treatment
|
||||
doublereal MC_cpCut_;
|
||||
|
||||
//! Parameter in the Molality Exp cutoff treatment
|
||||
doublereal CROP_ln_gamma_o_min;
|
||||
|
||||
//! Parameter in the Molality Exp cutoff treatment
|
||||
doublereal CROP_ln_gamma_o_max;
|
||||
|
||||
//! Parameter in the Molality Exp cutoff treatment
|
||||
doublereal CROP_ln_gamma_k_min;
|
||||
|
||||
//! Parameter in the Molality Exp cutoff treatment
|
||||
doublereal CROP_ln_gamma_k_max;
|
||||
|
||||
//! This is a boolean-type vector indicating whether
|
||||
|
|
@ -3500,7 +3506,10 @@ namespace Cantera {
|
|||
|
||||
//! Precalculate the IMS Cutoff parameters for typeCutoff = 2
|
||||
void calcIMSCutoffParams_();
|
||||
|
||||
//! Calculate molality cut-off parameters
|
||||
void calcMCCutoffParams_();
|
||||
|
||||
//! Utility function to assign an integer value from a string
|
||||
//! for the ElectrolyteSpeciesType field.
|
||||
/*!
|
||||
|
|
|
|||
|
|
@ -970,7 +970,7 @@ namespace Cantera {
|
|||
numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies();
|
||||
moleFractions_.resize(m_kk);
|
||||
fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
|
||||
fm_invert_ionForNeutral.resize(numNeutralMoleculeSpecies_);
|
||||
fm_invert_ionForNeutral.resize(m_kk);
|
||||
NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_);
|
||||
cationList_.resize(m_kk);
|
||||
anionList_.resize(m_kk);
|
||||
|
|
|
|||
|
|
@ -335,7 +335,12 @@ namespace Cantera {
|
|||
//! True if the number species has been set
|
||||
bool ready() const;
|
||||
|
||||
|
||||
//! Every time the mole fractions have changed, this routine
|
||||
//! will increment the stateMFNumber
|
||||
/*!
|
||||
* @param forceChange If this is true then the stateMFNumber always
|
||||
* changes. This defaults to false.
|
||||
*/
|
||||
void stateMFChangeCalc(bool forceChange = false);
|
||||
|
||||
//! Return the state number
|
||||
|
|
@ -423,7 +428,7 @@ namespace Cantera {
|
|||
|
||||
};
|
||||
|
||||
|
||||
//! Return the State Mole Fraction Number
|
||||
inline int State::stateMFNumber() const {
|
||||
return m_stateNum;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -127,11 +127,6 @@ namespace Cantera {
|
|||
|
||||
|
||||
|
||||
void WaterSSTP::constructPhase() {
|
||||
throw CanteraError("WaterSSTP::constructPhase()", "unimplemented");
|
||||
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* @param infile XML file containing the description of the
|
||||
|
|
|
|||
|
|
@ -400,13 +400,24 @@ namespace Cantera {
|
|||
*/
|
||||
virtual doublereal vaporFraction() const;
|
||||
|
||||
|
||||
//! Set the temperature of the phase
|
||||
/*!
|
||||
* The density and composition of the phase is constant during this
|
||||
* operator.
|
||||
*
|
||||
* @param temp Temperature (Kelvin)
|
||||
*/
|
||||
virtual void setTemperature(const doublereal temp);
|
||||
|
||||
//! Set the density of the phase
|
||||
/*!
|
||||
* The temperature and composition of the phase is constant during this
|
||||
* operator.
|
||||
*
|
||||
* @param dens value of the density in kg m-3
|
||||
*/
|
||||
virtual void setDensity(const doublereal dens);
|
||||
|
||||
void constructPhase();
|
||||
|
||||
|
||||
//! Initialization of a pure water phase using an
|
||||
//! xml file.
|
||||
|
|
|
|||
|
|
@ -47,6 +47,7 @@ typedef int ftnlen; // Fortran hidden string length type
|
|||
#undef HAS_SUNDIALS
|
||||
#undef SUNDIALS_VERSION_22
|
||||
#undef SUNDIALS_VERSION_23
|
||||
#undef SUNDIALS_VERSION_24
|
||||
|
||||
//-------- LAPACK / BLAS ---------
|
||||
|
||||
|
|
|
|||
74
configure
vendored
74
configure
vendored
File diff suppressed because one or more lines are too long
48
configure.in
48
configure.in
|
|
@ -365,29 +365,45 @@ sundials_lib_dep=
|
|||
|
||||
|
||||
if test ${use_sundials} = 1; then
|
||||
AC_DEFINE(HAS_SUNDIALS)
|
||||
echo "using CVODES from SUNDIALS... Sensitivity analysis enabled."
|
||||
|
||||
CVODE_LIBS='-lsundials_cvodes -lsundials_nvecserial'
|
||||
IDA_LIBS='-lsundials_ida -lsundials_nvecserial'
|
||||
|
||||
if test "$SUNDIALS_VERSION" = "2.2"; then
|
||||
AC_DEFINE(SUNDIALS_VERSION_22)
|
||||
sundials_include='-I'${SUNDIALS_HOME}'/include -I'${SUNDIALS_HOME}'/include/sundials -I'${SUNDIALS_HOME}'/include/cvodes -I'${SUNDIALS_HOME}'/include/ida'
|
||||
echo "sundials include directory: " ${sundials_include}
|
||||
echo "sundials library directory: " $SUNDIALS_LIB_DIR
|
||||
sundials_lib_dir=$SUNDIALS_LIB_DIR
|
||||
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
|
||||
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
|
||||
AC_DEFINE(HAS_SUNDIALS)
|
||||
AC_DEFINE(SUNDIALS_VERSION_22)
|
||||
sundials_include='-I'${SUNDIALS_HOME}'/include -I'${SUNDIALS_HOME}'/include/sundials -I'${SUNDIALS_HOME}'/include/cvodes -I'${SUNDIALS_HOME}'/include/ida'
|
||||
echo "sundials include directory: " ${sundials_include}
|
||||
echo "sundials library directory: " $SUNDIALS_LIB_DIR
|
||||
sundials_lib_dir=$SUNDIALS_LIB_DIR
|
||||
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
|
||||
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
|
||||
elif test "$SUNDIALS_VERSION" = "2.3"; then
|
||||
AC_DEFINE(HAS_SUNDIALS)
|
||||
AC_DEFINE(SUNDIALS_VERSION_23)
|
||||
sundials_include='-I'${SUNDIALS_INC_DIR}
|
||||
echo "sundials include directory: " ${sundials_include}
|
||||
echo "sundials library directory: " $SUNDIALS_LIB_DIR
|
||||
sundials_lib_dir=$SUNDIALS_LIB_DIR
|
||||
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
|
||||
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
|
||||
# python tools/src/sundials_version.py $SUNDIALS_HOME
|
||||
elif test "$SUNDIALS_VERSION" = "2.4"; then
|
||||
AC_DEFINE(HAS_SUNDIALS)
|
||||
AC_DEFINE(SUNDIALS_VERSION_24)
|
||||
sundials_include='-I'${SUNDIALS_INC_DIR}
|
||||
echo "sundials include directory: " ${sundials_include}
|
||||
echo "sundials library directory: " $SUNDIALS_LIB_DIR
|
||||
sundials_lib_dir=$SUNDIALS_LIB_DIR
|
||||
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
|
||||
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
|
||||
# python tools/src/sundials_version.py $SUNDIALS_HOME
|
||||
else
|
||||
AC_DEFINE(SUNDIALS_VERSION_23)
|
||||
sundials_include='-I'${SUNDIALS_INC_DIR}
|
||||
echo "sundials include directory: " ${sundials_include}
|
||||
echo "sundials library directory: " $SUNDIALS_LIB_DIR
|
||||
sundials_lib_dir=$SUNDIALS_LIB_DIR
|
||||
sundials_lib="-L$SUNDIALS_LIB_DIR -lsundials_cvodes -lsundials_ida -lsundials_nvecserial"
|
||||
sundials_lib_dep="$SUNDIALS_LIB_DIR/libsundials_cvodes.a $SUNDIALS_LIB_DIR/libsundials_ida.a $SUNDIALS_LIB_DIR/libsundials_nvecserial.a"
|
||||
# python tools/src/sundials_version.py $SUNDIALS_HOME
|
||||
echo "ERROR: unknown or unsupported sundials version #: $SUNDIALS_VERSION"
|
||||
echo " Supported versions are 2.2, 2.3, and 2.4"
|
||||
echo " Please fix or turn off the sundials option by setting USE_SUNDIALS to no"
|
||||
use_sundials=0
|
||||
fi
|
||||
fi
|
||||
|
||||
|
|
|
|||
|
|
@ -315,11 +315,11 @@ USE_SUNDIALS=${USE_SUNDIALS:='default'}
|
|||
|
||||
|
||||
# It is recommended that you install the newest release of sundials
|
||||
# (currently 2.3.0) before building Cantera. But if you want to use an
|
||||
# (currently 2.4.0) before building Cantera. But if you want to use an
|
||||
# older version, set SUNDIALS_VERSION to the version you have.
|
||||
# Acceptable values are '2.2' and '2.3' only; anything else will cause
|
||||
# Cantera to not use sundials.
|
||||
SUNDIALS_VERSION=${SUNDIALS_VERSION:='2.3'}
|
||||
# Acceptable values are '2.2', '2.3', or '2.4' ; anything else will cause
|
||||
# Cantera to
|
||||
SUNDIALS_VERSION=${SUNDIALS_VERSION:='2.4'}
|
||||
|
||||
#-----------------------------------------------------------------
|
||||
# BLAS and LAPACK
|
||||
|
|
|
|||
|
|
@ -21,9 +21,9 @@ Number of reactions = 8
|
|||
3 4.8454e-06 3.0554e-06 6.8484e+04 4.0505e+04 Csoot-*
|
||||
4 2.6364e-05 2.1519e-05 1.1250e+04 5.8320e+03 Csoot-*
|
||||
5 1.3012e-03 1.2749e-03 1.4427e+03 7.3573e+02 Csoot-*
|
||||
6 4.7372e+00 4.7359e+00 5.8161e-01 2.9647e-01 Csoot-*
|
||||
7 6.1771e-08 3.1736e-08
|
||||
FIN 7 6.1771e-08 2.1591e-11 -- success
|
||||
6 4.7372e+00 4.7359e+00 5.8160e-01 2.9647e-01 Csoot-*
|
||||
7 6.4112e-08 3.2378e-08
|
||||
FIN 7 6.4112e-08 3.6822e-12 -- success
|
||||
Gas Temperature = 1.4e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
@ -81,8 +81,8 @@ Sum of coverages = 1
|
|||
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
|
||||
-----------------------------------------------------------------------------------
|
||||
1 5.3218e+03 2.7005e+03
|
||||
2 8.9325e-06 4.0777e-06
|
||||
FIN 2 8.9325e-06 6.2499e-12 -- success
|
||||
2 8.8765e-06 4.0574e-06
|
||||
FIN 2 8.8765e-06 8.4527e-11 -- success
|
||||
Gas Temperature = 1.4e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
@ -140,8 +140,8 @@ Sum of coverages = 1
|
|||
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
|
||||
-----------------------------------------------------------------------------------
|
||||
1 2.1569e+05 9.5571e+04
|
||||
2 1.7622e-04 2.0108e-04
|
||||
FIN 2 1.7622e-04 2.0335e-10 -- success
|
||||
2 7.8671e-04 4.6185e-04
|
||||
FIN 2 7.8671e-04 1.5306e-10 -- success
|
||||
Gas Temperature = 1.5e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
@ -198,8 +198,8 @@ Sum of coverages = 1
|
|||
|
||||
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
|
||||
-----------------------------------------------------------------------------------
|
||||
1 2.8875e-10 1.4792e-10
|
||||
FIN 1 2.8875e-10 1.2324e-10 -- success
|
||||
1 3.0635e-10 1.5694e-10
|
||||
FIN 1 3.0635e-10 1.2314e-10 -- success
|
||||
Gas Temperature = 1.5e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
|
|||
|
|
@ -19,9 +19,9 @@ Number of reactions = 8
|
|||
3 4.8454e-06 3.0554e-06 6.8484e+04 4.0505e+04 Csoot-*
|
||||
4 2.6364e-05 2.1519e-05 1.1250e+04 5.8320e+03 Csoot-*
|
||||
5 1.3012e-03 1.2749e-03 1.4427e+03 7.3573e+02 Csoot-*
|
||||
6 4.7372e+00 4.7359e+00 5.8161e-01 2.9647e-01 Csoot-*
|
||||
7 6.1771e-08 3.1736e-08
|
||||
FIN 7 6.1771e-08 2.1591e-11 -- success
|
||||
6 4.7372e+00 4.7359e+00 5.8160e-01 2.9647e-01 Csoot-*
|
||||
7 6.4112e-08 3.2378e-08
|
||||
FIN 7 6.4112e-08 3.6822e-12 -- success
|
||||
Gas Temperature = 1.4e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
@ -72,8 +72,8 @@ Sum of coverages = 1
|
|||
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
|
||||
-----------------------------------------------------------------------------------
|
||||
1 5.3218e+03 2.7005e+03
|
||||
2 8.9325e-06 4.0777e-06
|
||||
FIN 2 8.9325e-06 6.2499e-12 -- success
|
||||
2 8.8765e-06 4.0574e-06
|
||||
FIN 2 8.8765e-06 8.4527e-11 -- success
|
||||
Gas Temperature = 1.4e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
@ -124,8 +124,8 @@ Sum of coverages = 1
|
|||
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
|
||||
-----------------------------------------------------------------------------------
|
||||
1 2.1569e+05 9.5571e+04
|
||||
2 1.7622e-04 2.0108e-04
|
||||
FIN 2 1.7622e-04 2.0335e-10 -- success
|
||||
2 7.8671e-04 4.6185e-04
|
||||
FIN 2 7.8671e-04 1.5306e-10 -- success
|
||||
Gas Temperature = 1.5e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
@ -175,8 +175,8 @@ Sum of coverages = 1
|
|||
|
||||
Iter Time Del_t Damp DelX Resid Name-Time Name-Damp
|
||||
-----------------------------------------------------------------------------------
|
||||
1 2.8875e-10 1.4792e-10
|
||||
FIN 1 2.8875e-10 1.2324e-10 -- success
|
||||
1 3.0635e-10 1.5694e-10
|
||||
FIN 1 3.0635e-10 1.2314e-10 -- success
|
||||
Gas Temperature = 1.5e+03
|
||||
Gas Pressure = 1.01e+05
|
||||
Gas Phase: gas (0)
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue