research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

support for sensitivity analysis

Browse source
This commit is contained in:
Dave Goodwin 2005-11-10 15:06:33 +00:00
parent ed770b98c8
commit 6547372987
27 changed files with 356 additions and 136 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
Cantera/src/CVode.cpp
View file

@ -40,7 +40,7 @@ extern "C" {
double* ydata = N_VDATA(y);
double* ydotdata = N_VDATA(ydot);
Cantera::FuncEval* f = (Cantera::FuncEval*)f_data;
f->eval(t, ydata, ydotdata);
f->eval(t, ydata, ydotdata, NULL);
}
@ -71,7 +71,7 @@ extern "C" {
dy = 1.0/ewtdata[j];
ydata[j] = ysave + dy;
dy = ydata[j] - ysave;
func->eval(t, ydata, ydot);
func->eval(t, ydata, ydot, NULL);
for (i=0; i < N; i++) {
col_j[i] = (ydot[i] - fydata[i])/dy;
}
@ -212,6 +212,8 @@ namespace Cantera {
m_iter, m_itol, &m_reltol,
&m_abstols, m_data, NULL, TRUE, m_iopt,
m_ropt.begin(), NULL);
cout << "m_reltol = " << m_reltol << endl;
cout << "m_abstols = " << m_abstols << endl;
}
if (!m_cvode_mem) throw CVodeErr("CVodeMalloc failed.");

154
Cantera/src/CVodesIntegrator.cpp
View file

@ -1,11 +1,13 @@
/**
* @file CVodeInt.cpp
* @file CVodesIntegrator.cpp
*
*/
// Copyright 2001 California Institute of Technology
#include "CVodesIntegrator.h"
#include "stringUtils.h"
#include <iostream>
using namespace std;
@ -24,6 +26,21 @@ 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" {
/**
@ -39,8 +56,18 @@ extern "C" {
void *f_data) {
double* ydata = NV_DATA_S(y); //N_VDATA(y);
double* ydotdata = NV_DATA_S(ydot); //N_VDATA(ydot);
Cantera::FuncEval* f = (Cantera::FuncEval*)f_data;
f->eval(t, ydata, ydotdata, NULL);
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, d->m_pars.begin());
//}
//catch (...) {
//Cantera::showErrors();
//Cantera::error("Teminating execution");
//}
}
}
@ -64,9 +91,11 @@ namespace Cantera {
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_maxsteps(20000), m_np(0)
{
//m_ropt.resize(OPT_SIZE,0.0);
//m_iopt = new long[OPT_SIZE];
@ -77,9 +106,15 @@ namespace Cantera {
/// Destructor.
CVodesIntegrator::~CVodesIntegrator()
{
if (m_cvode_mem) CVodeFree(m_cvode_mem);
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)); //N_VFree(nv(m_y));
if (m_abstol) N_VDestroy_Serial(nv(m_abstol)); //N_VFree(nv(m_abstol));
delete m_fdata;
//delete[] m_iopt;
}
@ -109,6 +144,11 @@ namespace Cantera {
m_abstols = abstol;
}
void CVodesIntegrator::setSensitivityTolerances(double reltol, double abstol) {
m_reltolsens = reltol;
m_abstolsens = abstol;
}
void CVodesIntegrator::setProblemType(int probtype) {
m_type = probtype;
}
@ -151,6 +191,29 @@ namespace Cantera {
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;
m_yS = N_VNewVectorArray_Serial(m_np, nv);
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;
cout << "atol = " << atol[0] << " " << atol[m_np-1] << endl;
flag = CVodeSetSensTolerances(m_cvode_mem, CV_SS, rtol, atol.begin());
}
void CVodesIntegrator::initialize(double t0, FuncEval& func)
{
m_neq = func.neq();
@ -166,11 +229,13 @@ namespace Cantera {
// 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)));
//m_iopt[MXSTEP] = m_maxsteps;
//m_iopt[MAXORD] = m_maxord;
//m_ropt[HMAX] = m_hmax;
//try {
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
//}
//catch (CanteraError) {
//showErrors();
//error("Teminating execution");
// }
if (m_cvode_mem) CVodeFree(m_cvode_mem);
@ -196,12 +261,17 @@ namespace Cantera {
// &m_abstols, m_data, NULL, TRUE, m_iopt,
// m_ropt.begin(), NULL);
}
if (flag == CV_MEM_FAIL) {
throw CVodesErr("Memory allocation failed.");
}
else if (flag == CV_ILL_INPUT) {
throw CVodesErr("Illegal value for CVodeMalloc input argument.");
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.");
}
cout << "returned from CVodeMalloc. m_cvode_mem = " << m_cvode_mem << endl;
if (m_type == DENSE + NOJAC) {
@ -219,11 +289,20 @@ namespace Cantera {
}
// pass a pointer to func in m_data
m_data = (void*)&func;
flag = CVodeSetFdata(m_cvode_mem, 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, m_fdata->m_pars.begin(),
NULL, NULL);
}
// set options
if (m_maxord > 0)
flag = CVodeSetMaxOrd(m_cvode_mem, m_maxord);
@ -237,16 +316,15 @@ namespace Cantera {
void CVodesIntegrator::reinitialize(double t0, FuncEval& func)
{
m_t0 = t0;
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
//try {
func.getInitialConditions(m_t0, m_neq, NV_DATA_S(nv(m_y)));
//}
//catch (CanteraError) {
//showErrors();
//error("Teminating execution");
//}
// set options
// m_iopt[MXSTEP] = m_maxsteps;
//m_iopt[MAXORD] = m_maxord;
//m_ropt[HMAX] = m_hmax;
//if (m_cvode_mem) CVodeFree(m_cvode_mem);
int result;
int result, flag;
if (m_itol == CV_SV) {
result = CVodeReInit(m_cvode_mem, cvodes_rhs, m_t0, nv(m_y),
m_itol, m_reltol,
@ -257,8 +335,8 @@ namespace Cantera {
m_itol, m_reltol,
&m_abstols);
}
if (result != 0) throw CVodesErr("CVReInit failed.");
cout << "returned from CVodeReInit. m_cvode_mem = " << m_cvode_mem << endl;
if (result != CV_SUCCESS) throw CVodesErr("CVReInit failed. result = "+int2str(result));
if (m_type == DENSE + NOJAC) {
long int N = m_neq;
@ -274,11 +352,6 @@ namespace Cantera {
throw CVodesErr("unsupported option");
}
// pass a pointer to func in m_data
m_data = (void*)&func;
long int flag = CVodeSetFdata(m_cvode_mem, m_data);
if (flag != CV_SUCCESS)
throw CVodesErr("CVodeSetFdata failed.");
// set options
if (m_maxord > 0)
@ -296,7 +369,10 @@ namespace Cantera {
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)
{
@ -310,10 +386,18 @@ namespace Cantera {
int CVodesIntegrator::nEvals() const {
long int ne;
return CVodeGetNumRhsEvals(m_cvode_mem, &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);
}
}

17
Cantera/src/CVodesIntegrator.h
View file

@ -23,8 +23,13 @@
#include "ctexceptions.h"
#include "ct_defs.h"
#include <nvector.h>
#include <nvector_serial.h>
namespace Cantera {
class FuncData;
/**
* Exception thrown when a CVODES error is encountered.
*/
@ -49,6 +54,7 @@ namespace Cantera {
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);
@ -65,8 +71,13 @@ namespace Cantera {
virtual void setMinStepSize(double hmin);
virtual void setMaxSteps(int nmax);
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;
@ -78,10 +89,14 @@ namespace Cantera {
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;
void* m_data;
FuncData* m_fdata;
N_Vector* m_yS;
int m_np;
};
} // namespace

4
Cantera/src/ChemEquil.cpp
View file

@ -352,6 +352,7 @@ namespace Cantera {
m_p2 = new DensityCalculator<thermo_t>;
break;
default:
endLogGroup();
throw CanteraError("equilibrate","illegal property pair.");
}
@ -364,6 +365,7 @@ namespace Cantera {
if (tempFixed) {
double tfixed = s.temperature();
if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) {
endLogGroup();
throw CanteraError("ChemEquil","Specified temperature ("
+fp2str(m_thermo->temperature())+" K) outside "
"valid range of "+fp2str(m_thermo->minTemp())+" K to "
@ -573,6 +575,8 @@ namespace Cantera {
fail++;
if (fail > 3) {
addLogEntry("dampStep","Failed 3 times. Giving up.");
endLogGroup(); // iteration
endLogGroup(); // equilibrate
s.restoreState(state);
throw CanteraError("equilibrate",
"Cannot find an acceptable Newton damping coefficient.");

5
Cantera/src/EdgeKinetics.cpp
View file

@ -263,7 +263,10 @@ namespace Cantera {
//cout << "i, beta = " << i << " " << m_beta[i] << endl;
if (eamod != 0.0 && m_E[i] != 0.0) {
ea = GasConstant * m_E[i];
if (eamod + ea < 0.0) eamod = -ea;
if (eamod + ea < 0.0) {
writelog("Warning: act energy mod too large");
eamod = -ea;
}
kf[irxn] *= exp(-eamod*rrt);
}
}

6
Cantera/src/FalloffMgr.h
View file

@ -86,8 +86,10 @@ namespace Cantera {
*/
void updateTemp(doublereal t, workPtr work) {
int i;
for (i = 0; i < m_n; i++) m_falloff[i]->updateTemp(t,
work + m_offset[i]);
for (i = 0; i < m_n; i++) {
m_falloff[i]->updateTemp(t,
work + m_offset[i]);
}
}
/**

43
Cantera/src/Func1.h
View file

@ -32,20 +32,30 @@ namespace Cantera {
/**
* Base class for 'functor' classes that evaluate a function of
* one variable.
* one variable.
*/
class Func1 {
public:
Func1() {}
virtual ~Func1() {}
/// Calls method eval to evaluate the function
doublereal operator()(doublereal t) { return eval(t); }
/// Evaluate the function.
virtual doublereal eval(doublereal t) { return 0.0; }
protected:
int m_n;
private:
};
/**
* 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) {
@ -86,23 +96,6 @@ namespace Cantera {
return r;
}
// virtual string show(doublereal t) {
// int n;
// string s = "";
// doublereal r = m_c[m_n-1];
// for (n = m_n-1; n >= 0; n--) {
// s += fp2str(m_c[n]);
// if (n > 0) s += "*x";
// if (n > 1) s += "^"+int2str(n);
// if (n > 0) {
// if (m_c[n] < 0.0) s += " - ";
// else
// r *= t;
// r += m_c[m_n - n - 1];
// }
// return r;
// }
protected:
int m_n;
vector_fp m_c;
@ -110,7 +103,12 @@ namespace Cantera {
/**
* Fourier cosine/sin series.
* 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:
@ -146,6 +144,9 @@ namespace Cantera {
/**
* 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:

6
Cantera/src/FuncEval.h
View file

@ -37,8 +37,9 @@ namespace Cantera {
* @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)=0;
virtual void eval(double t, double* y, double* ydot, double* p)=0;
/**
* Fill the solution vector with the initial conditions
@ -49,6 +50,9 @@ namespace Cantera {
/** Number of equations. */
virtual int neq()=0;
/// Number of parameters.
virtual int nparams() { return 0; }
protected:
private:

8
Cantera/src/GRI_30_Kinetics.cpp
View file

@ -25,10 +25,6 @@
using namespace std;
#ifdef HAVE_INTEL_MKL
#include "mkl_vml.h"
#endif
namespace Cantera {
/**
@ -44,11 +40,9 @@ namespace Cantera {
doublereal logT = log(T);
m_kdata->m_logc_ref = m_kdata->m_logp_ref - logT;
update_rates(T, logT, m_kdata->m_rfn.begin());
m_falloff_low_rates.update(T, logT, m_kdata->m_rfn_low.begin());
m_falloff_low_rates.update(T, logT, m_kdata->m_rfn_low.begin());
m_falloff_high_rates.update(T, logT, m_kdata->m_rfn_high.begin());
m_falloffn.updateTemp(T, m_kdata->falloff_work.begin());
m_kdata->m_temp = T;
gri30_updateKc();
m_kdata->m_ROP_ok = false;

7
Cantera/src/ImplicitChem.cpp
View file

@ -16,7 +16,7 @@
#endif
#include "ImplicitChem.h"
#include "CVode.h"
#include "Integrator.h"
namespace Cantera {
@ -24,7 +24,7 @@ namespace Cantera {
: FuncEval(), m_kin(&kin), m_thermo(&therm), m_integ(0),
m_atol(1.e-15), m_rtol(1.e-7), m_maxstep(0.0), m_energy(false)
{
m_integ = new CVodeInt;
m_integ = newIntegrator("CVODE"); //CVodeInt;
//m_mix = &kin.phase();
m_wt = m_thermo->molecularWeights();
@ -74,7 +74,8 @@ namespace Cantera {
/**
* Called by the integrator to evaluate ydot given y at time 'time'.
*/
void ImplicitChem::eval(doublereal time, doublereal* y, doublereal* ydot)
void ImplicitChem::eval(doublereal time, doublereal* y,
doublereal* ydot, doublereal* p)
{
updateState(y); // synchronize the mixture state with y
m_thermo->setPressure(m_press);

5
Cantera/src/ImplicitChem.h
View file

@ -17,7 +17,7 @@
#endif
#include "FuncEval.h"
#include "CVode.h"
#include "Integrator.h"
#include "Kinetics.h"
#include "ThermoPhase.h"
@ -83,7 +83,8 @@ namespace Cantera {
// overloaded methods of class FuncEval
virtual int neq() { return m_nsp; }
virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
virtual void eval(doublereal t, doublereal* y, doublereal* ydot,
doublereal* p);
virtual void getInitialConditions(doublereal t0, size_t leny,
doublereal* y);

10
Cantera/src/ImplicitSurfChem.cpp
View file

@ -17,7 +17,8 @@
#endif
#include "ImplicitSurfChem.h"
#include "CVode.h"
#include "Integrator.h"
namespace Cantera {
@ -41,7 +42,7 @@ namespace Cantera {
nt = k[n]->nTotalSpecies();
if (nt > ntmax) ntmax = nt;
}
m_integ = new CVodeInt;
m_integ = newIntegrator("CVODE");// CVodeInt;
// use backward differencing, with a full Jacobian computed
// numerically, and use a Newton linear iterator
@ -88,12 +89,13 @@ namespace Cantera {
* Called by the integrator to evaluate ydot given y at time 'time'.
*/
void ImplicitSurfChem::eval(doublereal time, doublereal* y,
doublereal* ydot)
doublereal* ydot, doublereal* p)
{
int n;
updateState(y); // synchronize the surface state(s) with y
doublereal rs0, sum;
int loc, k, kstart;
for (int n = 0; n < m_nsurf; n++) {
for (n = 0; n < m_nsurf; n++) {
rs0 = 1.0/m_surf[n]->siteDensity();
m_kin[n]->getNetProductionRates(m_work.begin());
kstart = m_kin[n]->kineticsSpeciesIndex(0,m_surfindex[n]);

5
Cantera/src/ImplicitSurfChem.h
View file

@ -20,7 +20,7 @@
#endif
#include "FuncEval.h"
#include "CVode.h"
#include "Integrator.h"
#include "InterfaceKinetics.h"
#include "SurfPhase.h"
@ -82,7 +82,8 @@ namespace Cantera {
// overloaded methods of class FuncEval
virtual int neq() { return m_nv; }
virtual void eval(doublereal t, doublereal* y, doublereal* ydot);
virtual void eval(doublereal t, doublereal* y, doublereal* ydot,
doublereal* p);
virtual void getInitialConditions(doublereal t0,
size_t leny, doublereal* y);

13
Cantera/src/Integrator.h
View file

@ -85,6 +85,9 @@ namespace Cantera {
virtual void setTolerances(doublereal reltol, doublereal abstol)
{ warn("setTolerances"); }
virtual void setSensitivityTolerances(doublereal reltol, doublereal abstol)
{ warn("setSensitivityTolerances"); }
/**
* Set problem type.
*/
@ -157,6 +160,13 @@ namespace Cantera {
virtual void setMaxSteps(int nmax)
{ warn("setMaxStep"); }
virtual int nSensParams()
{ warn("nSensParams()"); return 0; }
virtual double sensitivity(int k, int p) {
warn("sensitivity"); return 0.0;
}
private:
doublereal m_dummy;
@ -167,6 +177,9 @@ namespace Cantera {
};
// defined in ODE_integrators.cpp
Integrator* newIntegrator(string itype);
} // namespace
#endif

9
Cantera/src/Makefile.in
View file

@ -18,6 +18,8 @@ CANTERA_LIB = @buildlib@/libcantera.a
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) $(LOCAL_DEFNS)
EXT = ../../ext
do_ranlib = @DO_RANLIB@
USE_SUNDIALS = @use_sundials@
SUNDIALS_INC = @sundials_include@
#----------------------
# kernel components
@ -87,8 +89,8 @@ RPATH = $(RPATH_OBJ)
# solvers
SOLVERS_OBJ = CVode.o BandMatrix.o
SOLVERS_H = CVode.h BandMatrix.h Integrator.h
SOLVERS_OBJ = ODE_integrators.o BandMatrix.o
SOLVERS_H = BandMatrix.h Integrator.h
SOLVERS = $(SOLVERS_OBJ)
# 1D flow models
@ -183,6 +185,9 @@ ALL_H = $(BASE_H) $(THERMO_H) $(KINETICS_H) $(HETEROKIN_H) \
.cpp.o:
@CXX@ -c $< $(CXX_INCLUDES) $(CXX_FLAGS)
ODE_integrators.o:
@CXX@ -c ODE_integrators.cpp $(CXX_INCLUDES) $(SUNDIALS_INC) $(CXX_FLAGS)
lib: $(OBJ_LIB)
$(RM) $(CANTERA_LIB)
@ARCHIVE@ $(CANTERA_LIB) *.o > /dev/null

2
Cantera/src/MultiPhase.cpp
View file

@ -534,7 +534,6 @@ namespace Cantera {
"No convergence for T");
}
else if (XY == SP) {
writelog("SP\n");
s0 = entropy();
start = true;
Tlow = 1.0; // m_Tmin; // lower bound on T
@ -610,7 +609,6 @@ namespace Cantera {
"No convergence for T");
}
// else if (XY == SP) {
// if (loglevel > 0) {
// addLogEntry("problem type","fixed S,P");

3
Cantera/src/NasaPoly1.h
View file

@ -13,6 +13,7 @@
#ifndef CT_NASAPOLY1_H
#define CT_NASAPOLY1_H
#include "global.h"
#include "SpeciesThermoInterpType.h"
namespace Cantera {
@ -136,6 +137,8 @@ namespace Cantera {
cp_R[m_index] = cp;
h_RT[m_index] = h;
s_R[m_index] = s;
//writelog("NASA1: for species "+int2str(m_index)+", h_RT = "+
// fp2str(h)+"\n");
}
/**

3
Cantera/src/NasaThermo.h
View file

@ -17,6 +17,7 @@
#include "NasaPoly1.h"
#include "speciesThermoTypes.h"
#include "polyfit.h"
#include "global.h"
namespace Cantera {
@ -70,6 +71,8 @@ namespace Cantera {
doublereal minTemp, doublereal maxTemp,
doublereal refPressure) {
//writelog("installing NASA thermo for species "+name+"\n");
//writelog(" index = "+int2str(index)+"\n");
int imid = int(c[0]); // midpoint temp converted to integer
int igrp = m_index[imid]; // has this value been seen before?
if (igrp == 0) { // if not, prepare new group

30
Cantera/src/ODE_integrators.cpp
View file

@ -7,19 +7,19 @@
#include "CVode.cpp"
#endif
// namespace Cantera {
namespace Cantera {
// Integrator* newIntegrator(string itype) {
// if (itype == "CVODE") {
// #ifdef HAS_SUNDIALS
// return new CVodesIntegrator();
// #else
// return new CVodeInt();
// #endif
// }
// else {
// throw CanteraError("newIntegrator",
// "unknown ODE integrator: "+itype);
// }
// }
// }
Integrator* newIntegrator(string itype) {
if (itype == "CVODE") {
#ifdef HAS_SUNDIALS
return new CVodesIntegrator();
#else
return new CVodeInt();
#endif
}
else {
throw CanteraError("newIntegrator",
"unknown ODE integrator: "+itype);
}
}
}

102
Cantera/src/ReactionPath.cpp
View file

@ -450,7 +450,7 @@ namespace Cantera {
}
}
s << " label = " << "\"" << "Scale = "
<< flmax << "\";" << endl; //\\l\\l created with Cantera (www.cantera.org)\\l\";"
<< flmax << "\\l " << title << "\";" << endl; //created with Cantera (www.cantera.org)\\l\";"
s << " fontname = \""+m_font+"\";" << endl << "}" << endl;
}
@ -658,21 +658,72 @@ namespace Cantera {
g.fmt(out, m_elementSymbols);
}
void ReactionPathBuilder::findElements(Kinetics& kin) {
string ename;
m_enamemap.clear();
m_nel = 0;
int i, np = kin.nPhases();
ThermoPhase* p;
map<string, int> enamemap;
for (i = 0; i < np; i++) {
p = &kin.thermo(i);
// iterate over the elements in this phase
int m, nel = p->nElements();
for (m = 0; m < nel; m++) {
ename = p->elementName(m);
// if no entry is found for this element name, then
// it is a new element. In this case, add the name
// to the list of names, increment the element count,
// and add an entry to the name->(index+1) map.
if (m_enamemap.find(ename) == m_enamemap.end()) {
m_enamemap[ename] = m_nel + 1;
m_elementSymbols.push_back(ename);
m_nel++;
}
}
}
m_atoms.resize(kin.nTotalSpecies(), m_nel, 0.0);
string sym;
int k, ip, nsp, mlocal, kp, m;
// iterate over the elements
for (m = 0; m < m_nel; m++) {
sym = m_elementSymbols[m];
k = 0;
// iterate over the phases
for (ip = 0; ip < np; ip++) {
phase_t* p = &kin.thermo(ip);
nsp = p->nSpecies();
mlocal = p->elementIndex(sym);
for (kp = 0; kp < nsp; kp++) {
if (mlocal >= 0) {
m_atoms(k, m) = p->nAtoms(kp, mlocal);
}
k++;
}
}
}
}
int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) {
//m_warn.clear();
m_transfer.clear();
const Kinetics::thermo_t& ph = kin.thermo();
m_nel = ph.nElements();
m_ns = ph.nSpecies();
m_nr = kin.nReactions();
//const Kinetics::thermo_t& ph = kin.thermo();
m_elementSymbols.clear();
findElements(kin);
//m_nel = ph.nElements();
m_ns = kin.nTotalSpecies(); //ph.nSpecies();
m_nr = kin.nReactions();
int m, i;
for (m = 0; m < m_nel; m++) {
m_elementSymbols.push_back(ph.elementName(m));
}
//for (m = 0; m < m_nel; m++) {
// m_elementSymbols.push_back(ph.elementName(m));
//}
// all reactants / products, even ones appearing on both sides
// of the reaction
@ -736,7 +787,7 @@ namespace Cantera {
for (n = 0; n < nrnet; n++) {
k = m_reac[i][n];
for (int m = 0; m < m_nel; m++) {
m_elatoms(m,i) += ph.nAtoms(k,m);
m_elatoms(m,i) += m_atoms(k,m); //ph.nAtoms(k,m);
}
}
}
@ -746,7 +797,7 @@ namespace Cantera {
m_sgroup.resize(m_ns);
int j;
for (j = 0; j < m_ns; j++) {
for (int m = 0; m < m_nel; m++) comp[m] = int(ph.nAtoms(j,m));
for (int m = 0; m < m_nel; m++) comp[m] = int(m_atoms(j,m)); //ph.nAtoms(j,m));
m_sgroup[j] = Group(comp);
}
@ -768,10 +819,11 @@ namespace Cantera {
nar = 0;
nap = 0;
for (j = 0; j < nr; j++) {
if (ph.nAtoms(m_reac[i][j],m) > 0) nar++;
// if (ph.nAtoms(m_reac[i][j],m) > 0) nar++;
if (m_atoms(m_reac[i][j],m) > 0) nar++;
}
for (j = 0; j < np; j++) {
if (ph.nAtoms(m_prod[i][j],m) > 0) nap++;
if (m_atoms(m_prod[i][j],m) > 0) nap++;
}
if (nar > 1 && nap > 1) {
m_determinate[i] = false; break;
@ -812,19 +864,19 @@ namespace Cantera {
doublereal threshold = 0.0;
bool fwd_incl, rev_incl, force_incl;
const Kinetics::thermo_t& ph = s.thermo();
int m = ph.elementIndex(element);
// const Kinetics::thermo_t& ph = s.thermo();
int m = m_enamemap[element]-1; //ph.elementIndex(element);
r.element = element;
if (m < 0) return -1;
//int k;
int kk = ph.nSpecies();
int kk = s.nTotalSpecies();
s.getFwdRatesOfProgress(m_ropf.begin());
s.getRevRatesOfProgress(m_ropr.begin());
ph.getMoleFractions(m_x.begin());
//ph.getMoleFractions(m_x.begin());
//doublereal sum = 0.0;
//for (k = 0; k < kk; k++) {
@ -869,7 +921,7 @@ namespace Cantera {
revlabel = "";
for (l = 0; l < np; l++) {
if (l != kp)
revlabel += " + "+ ph.speciesName(m_prod[i][l]);
revlabel += " + "+ s.kineticsSpeciesName(m_prod[i][l]);
}
if (s.reactionType(i) == THREE_BODY_RXN)
revlabel += " + M ";
@ -882,8 +934,8 @@ namespace Cantera {
// element m, and both are allowed to appear in
// the diagram
if ((kkr != kkp) && (ph.nAtoms(kkr,m) > 0
&& ph.nAtoms(kkp,m) > 0)
if ((kkr != kkp) && (m_atoms(kkr,m) > 0
&& m_atoms(kkp,m) > 0)
&& status[kkr] >= 0 && status[kkp] >= 0)
{
@ -894,8 +946,8 @@ namespace Cantera {
// reactant species was the source of a
// given m-atom in the product
if ( (ph.nAtoms(kkp,m) < m_elatoms(m, i)) &&
(ph.nAtoms(kkr,m) < m_elatoms(m, i)) )
if ( (m_atoms(kkp,m) < m_elatoms(m, i)) &&
(m_atoms(kkr,m) < m_elatoms(m, i)) )
{
map<int, map<int, Group> >& g = m_transfer[i];
if (g.empty()) {
@ -926,7 +978,7 @@ namespace Cantera {
// the same expression.
else {
f = ph.nAtoms(kkp,m) * ph.nAtoms(kkr,m) / m_elatoms(m, i);
f = m_atoms(kkp,m) * m_atoms(kkr,m) / m_elatoms(m, i);
}
fwd = ropf*f;
@ -940,10 +992,10 @@ namespace Cantera {
if (fwd_incl || rev_incl)
{
if (!r.hasNode(kkr)) {
r.addNode(kkr, ph.speciesName(kkr), m_x[kkr]);
r.addNode(kkr, s.kineticsSpeciesName(kkr), m_x[kkr]);
}
if (!r.hasNode(kkp)) {
r.addNode(kkp, ph.speciesName(kkp), m_x[kkp]);
r.addNode(kkp, s.kineticsSpeciesName(kkp), m_x[kkp]);
}
}
if (fwd_incl) {

3
Cantera/src/ReactionPath.h
View file

@ -252,6 +252,7 @@ namespace Cantera {
void writeGroup(ostream& out, const Group& g);
protected:
void findElements(Kinetics& kin);
int m_nr;
int m_ns;
@ -268,6 +269,8 @@ namespace Cantera {
// map<int, int> m_warn;
map<int, map<int, map<int, Group> > > m_transfer;
vector<bool> m_determinate;
Array2D m_atoms;
map<string,int> m_enamemap;
};
}

3
Cantera/src/SpeciesThermoFactory.cpp
View file

@ -72,6 +72,7 @@ namespace Cantera {
iother = 1;
}
if (iother) {
writelog("returning new GeneralSpeciesThermo");
return new GeneralSpeciesThermo();
}
return newSpeciesThermo(NASA*inasa
@ -90,6 +91,7 @@ namespace Cantera {
}
}
if (iother) {
writelog("returning new GeneralSpeciesThermo");
return new GeneralSpeciesThermo();
}
return newSpeciesThermo(NASA*inasa
@ -110,6 +112,7 @@ namespace Cantera {
}
}
if (iother) {
writelog("returning new GeneralSpeciesThermo");
return new GeneralSpeciesThermo();
}
return newSpeciesThermo(NASA*inasa

3
Cantera/src/ThermoPhase.cpp
View file

@ -123,13 +123,14 @@ namespace Cantera {
dt = (u - intEnergy_mass())/cv_mass();
if (dt > 100.0) dt = 100.0;
else if (dt < -100.0) dt = -100.0;
setTemperature(temperature() + dt);
setTemperature(temperature() + 0.5*dt);
if (fabs(dt) < tol) {
return;
}
}
throw CanteraError("setState_UV",
"no convergence. dt = " + fp2str(dt)+"\n"
+"tol = "+fp2str(tol)+"\n"
+"u = "+fp2str(u)+" v = "+fp2str(v)+"\n");
}

12
Cantera/src/ctvector.h
View file

@ -22,6 +22,10 @@
#include <iostream>
#include "config.h"
#ifdef DEBUG_MODE
#include "ctexceptions.h"
#endif
namespace ct {
/**
@ -44,7 +48,13 @@ namespace ct {
ctvector_fp operator=(const ctvector_fp& x);
virtual ~ctvector_fp();
value_type operator[](size_t n) const { return _data[n]; }
value_type operator[](size_t n) const {
#ifdef DEBUG_MODE
if (n < 0 || n >= _size)
throw CanteraError("ctvector_fp","index out of range");
#endif
return _data[n];
}
value_type& operator[](size_t n) { return _data[n]; }
void resize(size_t n);

12
Cantera/src/misc.cpp
View file

@ -772,10 +772,11 @@ namespace Cantera {
if (title != "" && title != __app->loggroups.back()) {
writelog("Logfile error."
"\n beginLogGroup: "+ __app->loggroups.back()+
"\n endLogGroup; "+title+"\n");
"\n endLogGroup: "+title+"\n");
cout << "calling write_logfile..." << endl;
write_logfile("logerror");
__app->loggroups.clear();
__app->loglevels.clear();
//__app->loggroups.clear();
//__app->loglevels.clear();
}
else if (__app->loggroups.size() == 1) {
write_logfile(__app->loggroups.back()+"_log");
@ -796,7 +797,9 @@ namespace Cantera {
/// file will be overwritten. will be appended to the name.
/// @ingroup HTML_logs
void write_logfile(string file) {
if (!__app->xmllog) return;
if (!__app->xmllog) {
return;
}
string::size_type idot = file.rfind('.');
string ext = "";
string nm = file;
@ -827,6 +830,7 @@ namespace Cantera {
// Now we have a file name that does not correspond to any
// existing file. Open it as an output stream, and dump the
// XML (HTML) tree to it.
if (__app->xmllog) {
ofstream f(fname.c_str());
// go to the top of the tree, and write it all.

2
Cantera/src/transport/SolidTransport.cpp
View file

@ -62,7 +62,7 @@ namespace Cantera {
/**
* Compute the mobilities of the species from the diffusion coefficients,
* usind the Einstein relation.
* using the Einstein relation.
*/
void SolidTransport::getMobilities(doublereal* mobil) {
int k;

19
Cantera/src/transport/TransportBase.h
View file

@ -34,6 +34,7 @@ namespace Cantera {
const int cUserTransport = 500;
const int cFtnTransport = 600;
// forward reference
class XML_Writer;
@ -143,16 +144,26 @@ namespace Cantera {
/**
* Get the molar fluxes [kmol/m^2/s], given the thermodynamic
* state at two nearby points. @param state1 Array of
* temperature, density, and mass fractions for state 1.
* state at two nearby points.
* @param state1 Array of temperature, density, and mass
* fractions for state 1.
* @param state2 Array of temperature, density, and mass
* fractions for state 2. @param delta Distance from state 1
* to state 2 (m).
* fractions for state 2.
* @param delta Distance from state 1 to state 2 (m).
*/
virtual void getMolarFluxes(const doublereal* state1,
const doublereal* state2, doublereal delta,
doublereal* fluxes) { err("getMolarFluxes"); }
/**
* Get the mass fluxes [kg/m^2/s], given the thermodynamic
* state at two nearby points.
* @param state1 Array of temperature, density, and mass
* fractions for state 1.
* @param state2 Array of temperature, density, and mass
* fractions for state 2.
* @param delta Distance from state 1 to state 2 (m).
*/
virtual void getMassFluxes(const doublereal* state1,
const doublereal* state2, doublereal delta,
doublereal* fluxes) { err("getMassFluxes"); }