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

*** empty log message ***

Browse source
This commit is contained in:
Dave Goodwin 2003-09-13 16:24:09 +00:00
parent 88da08f3bb
commit b663e860a9
6 changed files with 487 additions and 908 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

310
Cantera/src/Surf1D.h
View file

@ -1,310 +0,0 @@
dep
#ifndef CT_SURF1D_H
#define CT_SURF1D_H
#include "Resid1D.h"
#include "SurfPhase.h"
#include "InterfaceKinetics.h"
#include "StFlow.h"
#include "OneDim.h"
#include "ctml.h"
namespace Cantera {
// A class for surface domains in one-dimensional simulations, The
// surface is zero-dimensional, and defined by a set of surface
// species coverages.
class Surf1D : public Resid1D {
public:
Surf1D(InterfaceKinetics* skin = 0) : Resid1D(1, 1, 0.0) {
m_type = cSurfType;
m_flow_left = 0;
m_flow_right = 0;
m_kin = 0;
m_sphase = 0;
if (skin) setKinetics(skin);
}
virtual ~Surf1D(){}
// Set the kinetics manager for the surface.
void setKinetics(InterfaceKinetics* kin) {
m_kin = kin;
int np = kin->nPhases();
m_sphase = 0;
for (int n = 0; n < np; n++) {
if (kin->phase(n).eosType() == cSurf) {
m_sphase = (SurfPhase*)&m_kin->phase(n);
m_nsurf = n;
}
else {
m_bulk.push_back(&kin->phase(n));
m_nbulk.push_back(n);
}
}
if (!m_sphase)
throw CanteraError("setKinetics","no surface phase defined");
m_nsp = m_sphase->nSpecies();
resize(m_nsp,1);
if (m_bulk.size() == 1) {
m_bulk.push_back(0);
}
}
void fixSpecies(int k, doublereal c) {
if (c >= 0.0) m_fixed_cov[k] = c;
m_do_surf_species[k] = false;
needJacUpdate();
}
void solveSpecies(int k) {
m_do_surf_species[k] = true;
needJacUpdate();
}
/// Set the surface temperature
void setTemperature(doublereal t) {
m_sphase->setTemperature(t);
needJacUpdate();
}
/// Temperature [K].
doublereal temperature() {
return m_sphase->temperature();
}
void setCoverages(doublereal* c) {
m_sphase->setCoverages(c);
copy(c, c + m_nsp, m_fixed_cov.begin());
}
void setMultiplier(int k, doublereal f) {
m_mult[k] = f;
needJacUpdate();
}
doublereal multiplier(int k) { return m_mult[k]; }
virtual string componentName(int n) const {
return m_sphase->speciesName(n);
}
virtual void init() {
if (m_index < 0) {
throw CanteraError("Surf1D",
"install in container before calling init.");
}
m_nsp = m_sphase->nSpecies();
resize(m_nsp,1);
m_mult.resize(m_nsp, 1.0);
m_do_surf_species.resize(m_nsp, true);
m_fixed_cov.resize(m_nsp, 1.0/m_nsp);
// set bounds
vector_fp lower(m_nsp, -1.e-3);
vector_fp upper(m_nsp, 1.0);
setBounds(m_nsp, lower.begin(), m_nsp, upper.begin());
// set tolerances
vector_fp rtol(m_nsp, 1e-4);
vector_fp atol(m_nsp, 1.e-10);
setTolerances(m_nsp, rtol.begin(), m_nsp, atol.begin());
m_left_nsp = 0;
m_right_nsp = 0;
// check for left and right flow objects
if (m_index > 0) {
Resid1D& r = container().domain(m_index-1);
if (r.domainType() == cFlowType) {
m_flow_left = (StFlow*)&r;
m_left_nv = m_flow_left->nComponents();
m_left_points = m_flow_left->nPoints();
m_left_loc = container().start(m_index-1);
m_left_nsp = m_left_nv - 4;
m_phase_left = &m_flow_left->phase();
m_molwt_left = m_phase_left->molecularWeights().begin();
if (m_phase_left == m_bulk[0])
m_start_left = m_kin->start(m_nbulk[0]);
else if (m_phase_left == m_bulk[1])
m_start_left = m_kin->start(m_nbulk[1]);
else
throw CanteraError("Surf1D::init",
"left gas does not match one in surface mechanism");
}
else
throw CanteraError("Surf1D::init",
"Surface domains can only be "
"connected to flow domains.");
}
if (m_index < container().nDomains() - 1) {
Resid1D& r = container().domain(m_index+1);
if (r.domainType() == cFlowType) {
m_flow_right = (StFlow*)&r;
m_right_nv = m_flow_right->nComponents();
m_right_loc = container().start(m_index+1);
m_right_nsp = m_right_nv - 4;
m_phase_right = &m_flow_right->phase();
m_molwt_right = m_phase_right->molecularWeights().begin();
if (m_phase_right == m_bulk[0])
m_start_right = m_kin->start(m_nbulk[0]);
else if (m_phase_right == m_bulk[1])
m_start_right = m_kin->start(m_nbulk[1]);
else
throw CanteraError("Surf1D::init",
"right gas does not match one in surface mechanism");
}
else
throw CanteraError("Surf1D::init",
"Surface domains can only be "
"connected to flow domains.");
}
m_work.resize(m_kin->nSpecies());
}
virtual void eval(int jg, doublereal* xg, doublereal* rg,
integer* diagg, doublereal rdt) {
int k;
if (jg >= 0 && (jg < firstPoint() - 2
|| jg > lastPoint() + 2)) return;
// start of local part of global arrays
doublereal* x = xg + loc();
doublereal* r = rg + loc();
integer* diag = diagg + loc();
// set the coverages
doublereal sum = 0.0;
for (k = 0; k < m_nsp; k++) {
m_work[k] = x[k];
sum += x[k];
}
m_sphase->setCoverages(m_work.begin());
// set the left gas state to the adjacent point
int leftloc = 0, rightloc = 0;
int pnt = 0;
if (m_flow_left) {
leftloc = m_flow_left->loc();
pnt = m_flow_left->nPoints() - 1;
m_flow_left->setGas(xg + leftloc, pnt);
}
if (m_flow_right) {
rightloc = m_flow_right->loc();
m_flow_right->setGas(xg + rightloc, 0);
}
m_kin->getNetProductionRates(m_work.begin());
doublereal rs0 = 1.0/m_sphase->siteDensity();
scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]);
bool enabled = true;
int ioffset = m_kin->start(m_nsurf); // m_left_nsp + m_right_nsp;
doublereal maxx = -1.0;
int imx = -1;
for (k = 0; k < m_nsp; k++) {
r[k] = m_work[k + ioffset] * m_sphase->size(k) * rs0;
r[k] -= rdt*(x[k] - prevSoln(k,0));
diag[k] = 1;
if (x[k] > maxx) {
maxx = x[k];
imx = k;
}
if (!m_do_surf_species[k]) {
r[k] = x[k] - m_fixed_cov[k];
diag[k] = 0;
enabled = false;
}
}
if (enabled) {
r[imx] = 1.0 - sum;
diag[imx] = 0;
}
// gas-phase residuals
doublereal rho;
if (m_flow_left) {
rho = m_phase_left->density();
doublereal rdz = 2.0/
(m_flow_left->z(m_left_points-1) -
m_flow_left->z(m_left_points - 2));
for (k = 0; k < m_left_nsp; k++)
m_work[k + m_start_left] *= m_molwt_left[k];
int ileft = loc() - m_left_nv;
// if the energy equation is enabled at this point,
// set the gas temperature to the surface temperature
if (m_flow_left->doEnergy(pnt)) {
rg[ileft + 2] = xg[ileft + 2] - m_sphase->temperature();
}
for (k = 1; k < m_left_nsp; k++) {
if (enabled && m_flow_left->doSpecies(k)) {
rg[ileft + 4 + k] += m_work[k + m_start_left];
//+= rdz*m_work[k + m_sp_left]/rho;
}
}
}
if (m_flow_right) {
for (k = 0; k < m_right_nsp; k++)
m_work[k + m_start_right] *= m_molwt_right[k];
int iright = loc() + m_nsp;
rg[iright + 2] -= m_sphase->temperature();
//r[iright + 3] = x[iright];
for (k = 0; k < m_right_nsp; k++) {
rg[iright + 4 + k] -= m_work[k + m_start_right];
}
}
}
virtual void save(XML_Node& o, doublereal* soln) {
doublereal* s = soln + loc();
XML_Node& surf = o.addChild("surface");
for (int k = 0; k < m_nsp; k++) {
ctml::addFloat(surf, componentName(k), s[k], "", "coverage",
0.0, 1.0);
}
}
protected:
InterfaceKinetics* m_kin;
SurfPhase* m_sphase;
StFlow *m_flow_left, *m_flow_right;
int m_left_nv, m_right_nv;
int m_left_loc, m_right_loc;
int m_left_points;
int m_nsp, m_left_nsp, m_right_nsp;
vector_fp m_work;
const doublereal *m_molwt_right, *m_molwt_left;
int m_sp_left, m_sp_right;
int m_start_left, m_start_right, m_start_surf;
ThermoPhase *m_phase_left, *m_phase_right;
vector<ThermoPhase*> m_bulk;
vector<int> m_nbulk;
int m_nsurf;
vector_fp m_mult;
vector<bool> m_do_surf_species;
vector_fp m_fixed_cov;
};
}
#endif

846
Cantera/src/oneD/StFlow.cpp
View file

@ -301,6 +301,28 @@ namespace Cantera {
}
void StFlow::_finalize(const doublereal* x) {
int k, j;
doublereal zz, tt;
int nz = m_zfix.size();
bool e = m_do_energy[0];
for (j = 0; j < m_points; j++) {
if (e || nz == 0)
setTemperature(j, T(x, j));
else {
zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
tt = linearInterp(zz, m_zfix, m_tfix);
setTemperature(j, tt);
}
for (k = 0; k < m_nsp; k++) {
setMassFraction(j, k, Y(x, k, j));
}
}
if (e) solveEnergyEqn();
}
//------------------------------------------------------
/**
* Evaluate the residual function for axisymmetric stagnation
@ -587,249 +609,240 @@ namespace Cantera {
}
void OneDFlow::eval(int jg, doublereal* xg, doublereal* rg, integer* diagg,
doublereal rdt) {
// void OneDFlow::eval(int jg, doublereal* xg, doublereal* rg, integer* diagg,
// doublereal rdt) {
static double elapsed;
// doublereal rtau = 1.e5;
// static double elapsed;
clock_t t0 = clock();
// clock_t t0 = clock();
// doublereal rdt_save = rdt;
if (jg >= 0) rdt = 0.0;
// if (jg >= 0) rdt = 0.0;
if (jg >= 0 && (jg < firstPoint() || jg > lastPoint())) return;
// if (jg >= 0 && (jg < firstPoint() || jg > lastPoint())) return;
// start of local part of global arrays
doublereal* x = xg + loc();
doublereal* rsd = rg + loc();
integer* diag = diagg + loc();
// // start of local part of global arrays
// doublereal* x = xg + loc();
// doublereal* rsd = rg + loc();
// integer* diag = diagg + loc();
int jmin, jmax, jpt;
jpt = jg - firstPoint();
// int jmin, jmax, jpt;
// jpt = jg - firstPoint();
for (int jj = 0; jj < m_points*m_nv; jj++) {
if (x[jj] < -1.e20 || x[jj] > 1.e20) {
showSolution(cout, x);
throw CanteraError("tlt","tlt");
}
}
// // the residual function is evaluated for jmin <= j <= jmax, and
// // properties and evaluated for j0 <= j <= j1.
// the residual function is evaluated for jmin <= j <= jmax, and
// properties and evaluated for j0 <= j <= j1.
// if (jg < 0) {
// jmin = 0;
// jmax = m_points - 1;
// }
// else {
// jmin = max(jpt-1,0);
// jmax = min(jpt+1,m_points-1);
// }
// int j0 = max(jmin-1,0);
// int j1 = min(jmax+1,m_points-1);
if (jg < 0) {
jmin = 0;
jmax = m_points - 1;
}
else {
jmin = max(jpt-1,0);
jmax = min(jpt+1,m_points-1);
}
int j0 = max(jmin-1,0);
int j1 = min(jmax+1,m_points-1);
// int j, k;
int j, k;
// //-----------------------------------------------------
// // compute properties needed in the residual equations
// //-----------------------------------------------------
//-----------------------------------------------------
// compute properties needed in the residual equations
//-----------------------------------------------------
// // for each point, synchronize the state of the fluid object
// // with the current solution values, and then use this object
// // to compute the density, mean molecular weight, and mean
// // specific heat at constant pressure.
// if (jpt < 0) updateThermo(x, j0, j1);
// for each point, synchronize the state of the fluid object
// with the current solution values, and then use this object
// to compute the density, mean molecular weight, and mean
// specific heat at constant pressure.
if (jpt < 0) updateThermo(x, j0, j1);
// // skip updating transport properties if a Jacobian is
// // being evaluated
// if (jpt < 0) updateTransport(x, j0, j1);
// skip updating transport properties if a Jacobian is
// being evaluated
if (jpt < 0) updateTransport(x, j0, j1);
// update the species diffusive mass fluxes
updateDiffFluxes(x, j0, j1);
// // update the species diffusive mass fluxes
// updateDiffFluxes(x, j0, j1);
//----------------------------------------------------
// evaluate the residual equations at all required
// grid points
//----------------------------------------------------
// //----------------------------------------------------
// // evaluate the residual equations at all required
// // grid points
// //----------------------------------------------------
doublereal sum, sum2, deltaz, dtdzj;
// doublereal sum, sum2, deltaz, dtdzj;
for (j = jmin; j <= jmax; j++) {
// for (j = jmin; j <= jmax; j++) {
//----------------------------------------------
// boundaries
//----------------------------------------------
// //----------------------------------------------
// // boundaries
// //----------------------------------------------
if (j == 0) {
setGas(x,0);
m_boundary[0]->eval(x, m_rho[0], m_flux.begin(),
rsd);
}
// if (j == 0) {
// setGas(x,0);
// m_boundary[0]->eval(x, m_rho[0], m_flux.begin(),
// rsd);
// }
else if (j == m_points - 1) {
m_boundary[1]->eval(x + index(0, j), m_rho[j],
m_flux.begin() + m_nsp*(j-1),
rsd + index(0, j));
}
// else if (j == m_points - 1) {
// m_boundary[1]->eval(x + index(0, j), m_rho[j],
// m_flux.begin() + m_nsp*(j-1),
// rsd + index(0, j));
// }
//------------------------------------------
// interior points
//------------------------------------------
// //------------------------------------------
// // interior points
// //------------------------------------------
else {
// else {
// continuity
rsd[index(c_offset_U,j)] = (rho_u(x,j-1) - rho_u(x,j));
// // continuity
// rsd[index(c_offset_U,j)] = (rho_u(x,j-1) - rho_u(x,j));
// radial velocity = 0
rsd[index(c_offset_V,j)] = V(x,j);
// // radial velocity = 0
// rsd[index(c_offset_V,j)] = V(x,j);
// species equations
getWdot(x,j);
// // species equations
// getWdot(x,j);
doublereal convec, diffus;
for (k = 0; k < m_nsp; k++) {
if (m_do_species[k]) {
// doublereal convec, diffus;
// for (k = 0; k < m_nsp; k++) {
// if (m_do_species[k]) {
convec = rho_u(x,j) * dYdz(x,k,j);
diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))/(z(j+1) - z(j-1));
rsd[index(c_offset_Y + k, j)] =
(m_wt[k]*wdot(k,j) - convec - diffus)/m_rho[j]
- rdt*(Y(x,k,j) - Y_prev(k,j));
diag[index(c_offset_Y + k, j)] = 1;
}
else
rsd[index(c_offset_Y+k,j)] = (Y(x,k,j) - Y_fixed(k,j));
}
// convec = rho_u(x,j) * dYdz(x,k,j);
// diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))/(z(j+1) - z(j-1));
// rsd[index(c_offset_Y + k, j)] =
// (m_wt[k]*wdot(k,j) - convec - diffus)/m_rho[j]
// - rdt*(Y(x,k,j) - Y_prev(k,j));
// diag[index(c_offset_Y + k, j)] = 1;
// }
// else
// rsd[index(c_offset_Y+k,j)] = (Y(x,k,j) - Y_fixed(k,j));
// }
// energy equation
// // energy equation
if (m_do_energy[j]) {
setGas(x,j);
// if (m_do_energy[j]) {
// setGas(x,j);
// heat release term
const vector_fp& h_RT = m_thermo->enthalpy_RT();
const vector_fp& cp_R = m_thermo->cp_R();
sum = 0.0;
sum2 = 0.0;
deltaz = (z(j+1) - z(j-1));
doublereal flxk;
for (k = 0; k < m_nsp; k++) {
flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
// // heat release term
// const vector_fp& h_RT = m_thermo->enthalpy_RT();
// const vector_fp& cp_R = m_thermo->cp_R();
// sum = 0.0;
// sum2 = 0.0;
// deltaz = (z(j+1) - z(j-1));
// doublereal flxk;
// for (k = 0; k < m_nsp; k++) {
// flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
sum += wdot(k,j)*h_RT[k];
sum2 += flxk*cp_R[k]/m_wt[k];
}
sum *= GasConstant * T(x,j);
dtdzj = (T(x,j+1) - T(x,j-1))/deltaz; // dTdz(x,j) + (m_dz[j-1]/deltaz)*(dTdz(x,j+1) - dTdz(x,j));
sum2 *= GasConstant * dtdzj;
rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj
- divHeatFlux(x,j) - sum - sum2;
rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
diag[index(c_offset_T, j)] = 1;
}
// sum += wdot(k,j)*h_RT[k];
// sum2 += flxk*cp_R[k]/m_wt[k];
// }
// sum *= GasConstant * T(x,j);
// dtdzj = (T(x,j+1) - T(x,j-1))/deltaz; // dTdz(x,j) + (m_dz[j-1]/deltaz)*(dTdz(x,j+1) - dTdz(x,j));
// sum2 *= GasConstant * dtdzj;
// rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj
// - divHeatFlux(x,j) - sum - sum2;
// rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
// rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
// diag[index(c_offset_T, j)] = 1;
// }
// lambda = 0
rsd[index(c_offset_L, j)] = lambda(x,j);
// // lambda = 0
// rsd[index(c_offset_L, j)] = lambda(x,j);
}
for (k = 0; k < m_nsp; k++) {
if (!m_do_species[k]) {
rsd[index(c_offset_Y+k,j)] =
(Y(x,k,j) - Y_fixed(k,j));
diag[index(c_offset_Y+k, j)] = 0;
}
}
if (!m_do_energy[j]) {
rsd[index(c_offset_T, j)] = (T(x,j) - T_fixed(j));
diag[index(c_offset_T, j)] = 0;
}
// }
// for (k = 0; k < m_nsp; k++) {
// if (!m_do_species[k]) {
// rsd[index(c_offset_Y+k,j)] =
// (Y(x,k,j) - Y_fixed(k,j));
// diag[index(c_offset_Y+k, j)] = 0;
// }
// }
// if (!m_do_energy[j]) {
// rsd[index(c_offset_T, j)] = (T(x,j) - T_fixed(j));
// diag[index(c_offset_T, j)] = 0;
// }
}
clock_t t1 = clock();
elapsed += double(t1 - t0)/CLOCKS_PER_SEC;
}
// }
// clock_t t1 = clock();
// elapsed += double(t1 - t0)/CLOCKS_PER_SEC;
// }
/**
* Update the transport properties at grid points in the range
* from j0 to j1, based on solution x.
*/
void OneDFlow::updateTransport(doublereal* x,int j0, int j1) {
int j;
if (m_transport_option == c_Mixav_Transport) {
for (j = j0; j < j1; j++) {
setGasAtMidpoint(x,j);
m_trans->getMixDiffCoeffs(m_diff.begin() + j*m_nsp);
m_tcon[j] = m_trans->thermalConductivity();
}
}
else if (m_transport_option == c_Multi_Transport) {
for (j = j0; j < j1; j++) {
setGasAtMidpoint(x,j);
m_trans->getMultiDiffCoeffs(m_nsp, m_diff.begin() + mindex(0,0,j));
m_tcon[j] = m_trans->thermalConductivity();
}
}
}
// /**
// * Update the transport properties at grid points in the range
// * from j0 to j1, based on solution x.
// */
// void OneDFlow::updateTransport(doublereal* x,int j0, int j1) {
// int j;
// if (m_transport_option == c_Mixav_Transport) {
// for (j = j0; j < j1; j++) {
// setGasAtMidpoint(x,j);
// m_trans->getMixDiffCoeffs(m_diff.begin() + j*m_nsp);
// m_tcon[j] = m_trans->thermalConductivity();
// }
// }
// else if (m_transport_option == c_Multi_Transport) {
// for (j = j0; j < j1; j++) {
// setGasAtMidpoint(x,j);
// m_trans->getMultiDiffCoeffs(m_nsp, m_diff.begin() + mindex(0,0,j));
// m_tcon[j] = m_trans->thermalConductivity();
// }
// }
// }
/**
* Print the solution.
*/
void StFlow::showSolution(ostream& s, const doublereal* x) {
int nn = m_nv/5;
int i, j, n;
char* buf = new char[100];
// /**
// * Print the solution.
// */
// void StFlow::showSolution(ostream& s, const doublereal* x) {
// int nn = m_nv/5;
// int i, j, n;
// char* buf = new char[100];
// The mean molecular weight is needed to convert
updateThermo(x, 0, m_points-1);
// // The mean molecular weight is needed to convert
// updateThermo(x, 0, m_points-1);
for (i = 0; i < nn; i++) {
drawline(s);
sprintf(buf, "\n z ");
s << buf;
for (n = 0; n < 5; n++) {
sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
s << buf;
}
drawline(s);
for (j = 0; j < m_points; j++) {
sprintf(buf, "\n %10.4g ",m_z[j]);
s << buf;
for (n = 0; n < 5; n++) {
sprintf(buf, " %10.4g ",component(x, i*5+n,j));
s << buf;
}
}
s << endl;
}
int nrem = m_nv - 5*nn;
drawline(s);
sprintf(buf, "\n z ");
s << buf;
for (n = 0; n < nrem; n++) {
sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
s << buf;
}
drawline(s);
for (j = 0; j < m_points; j++) {
sprintf(buf, "\n %10.4g ",m_z[j]);
s << buf;
for (n = 0; n < nrem; n++) {
sprintf(buf, " %10.4g ",component(x, nn*5+n,j));
s << buf;
}
}
s << endl;
}
// for (i = 0; i < nn; i++) {
// drawline(s);
// sprintf(buf, "\n z ");
// s << buf;
// for (n = 0; n < 5; n++) {
// sprintf(buf, " %10s ",componentName(i*5 + n).c_str());
// s << buf;
// }
// drawline(s);
// for (j = 0; j < m_points; j++) {
// sprintf(buf, "\n %10.4g ",m_z[j]);
// s << buf;
// for (n = 0; n < 5; n++) {
// sprintf(buf, " %10.4g ",component(x, i*5+n,j));
// s << buf;
// }
// }
// s << endl;
// }
// int nrem = m_nv - 5*nn;
// drawline(s);
// sprintf(buf, "\n z ");
// s << buf;
// for (n = 0; n < nrem; n++) {
// sprintf(buf, " %10s ", componentName(nn*5 + n).c_str());
// s << buf;
// }
// drawline(s);
// for (j = 0; j < m_points; j++) {
// sprintf(buf, "\n %10.4g ",m_z[j]);
// s << buf;
// for (n = 0; n < nrem; n++) {
// sprintf(buf, " %10.4g ",component(x, nn*5+n,j));
// s << buf;
// }
// }
// s << endl;
// }
/**
@ -945,32 +958,32 @@ namespace Cantera {
void StFlow::outputTEC(ostream &s, const doublereal* x,
string title, int zone) {
int j,k;
s << "TITLE = \"" + title + "\"" << endl;
s << "VARIABLES = \"Z (m)\"" << endl;
s << "\"u (m/s)\"" << endl;
s << "\"V (1/s)\"" << endl;
s << "\"T (K)\"" << endl;
s << "\"lambda\"" << endl;
// void StFlow::outputTEC(ostream &s, const doublereal* x,
// string title, int zone) {
// int j,k;
// s << "TITLE = \"" + title + "\"" << endl;
// s << "VARIABLES = \"Z (m)\"" << endl;
// s << "\"u (m/s)\"" << endl;
// s << "\"V (1/s)\"" << endl;
// s << "\"T (K)\"" << endl;
// s << "\"lambda\"" << endl;
for (k = 0; k < m_nsp; k++) {
s << "\"" << m_thermo->speciesName(k) << "\"" << endl;
}
s << "ZONE T=\"c" << zone << "\"" << endl;
s << " I=" << m_points << ",J=1,K=1,F=POINT" << endl;
s << "DT=(SINGLE SINGLE SINGLE SINGLE";
for (k = 0; k < m_nsp; k++) s << " SINGLE";
s << " )" << endl;
for (j = 0; j < m_points; j++) {
s << z(j) << " ";
for (k = 0; k < m_nv; k++) {
s << component(x, k, j) << " ";
}
s << endl;
}
}
// for (k = 0; k < m_nsp; k++) {
// s << "\"" << m_thermo->speciesName(k) << "\"" << endl;
// }
// s << "ZONE T=\"c" << zone << "\"" << endl;
// s << " I=" << m_points << ",J=1,K=1,F=POINT" << endl;
// s << "DT=(SINGLE SINGLE SINGLE SINGLE";
// for (k = 0; k < m_nsp; k++) s << " SINGLE";
// s << " )" << endl;
// for (j = 0; j < m_points; j++) {
// s << z(j) << " ";
// for (k = 0; k < m_nv; k++) {
// s << component(x, k, j) << " ";
// }
// s << endl;
// }
// }
string StFlow::componentName(int n) const {
@ -983,207 +996,202 @@ namespace Cantera {
if (n >= (int) c_offset_Y && n < (int) (c_offset_Y + m_nsp)) {
return m_thermo->speciesName(n - c_offset_Y);
}
// if (m_do_species[n - c_offset_Y])
// return m_thermo->speciesName(n - c_offset_Y)+" ";
// else
// return m_thermo->speciesName(n - c_offset_Y)+" *";
//}
else
return "<unknown>";
}
}
/**
* Returns true if all necessary parameters have been set; otherwise it
* throws an exception.
*/
bool StFlow::ready() {
if (m_press < 0.0) {
throw CanteraError("StFlow::ready",
"pressure not specified - call setPressure");
return false;
}
if (m_points == 0) {
throw CanteraError("StFlow::ready",
"grid not specified - call setupGrid");
return false;
}
if (m_nsp < 0) {
throw CanteraError("StFlow::ready",
"fluid not specified - call specifyFluid");
return false;
}
if (m_boundary[0] == 0 || m_boundary[1] == 0) {
throw CanteraError("StFlow::ready",
"boundaries not specified - call setBoundary");
return false;
}
m_ok = true;
return m_ok;
}
// /**
// * Returns true if all necessary parameters have been set; otherwise it
// * throws an exception.
// */
// bool StFlow::ready() {
// if (m_press < 0.0) {
// throw CanteraError("StFlow::ready",
// "pressure not specified - call setPressure");
// return false;
// }
// if (m_points == 0) {
// throw CanteraError("StFlow::ready",
// "grid not specified - call setupGrid");
// return false;
// }
// if (m_nsp < 0) {
// throw CanteraError("StFlow::ready",
// "fluid not specified - call specifyFluid");
// return false;
// }
// if (m_boundary[0] == 0 || m_boundary[1] == 0) {
// throw CanteraError("StFlow::ready",
// "boundaries not specified - call setBoundary");
// return false;
// }
// m_ok = true;
// return m_ok;
// }
void StFlow::restore(int job,
string fname, string id, int& size_z, doublereal* z,
int& size_soln, doublereal* soln) {
// void StFlow::restore(int job,
// string fname, string id, int& size_z, doublereal* z,
// int& size_soln, doublereal* soln) {
vector<string> ignored;
int nsp = m_thermo->nSpecies();
vector_int did_species(nsp, 0);
// vector<string> ignored;
// int nsp = m_thermo->nSpecies();
// vector_int did_species(nsp, 0);
ifstream s(fname.c_str());
if (!s)
throw CanteraError("StFlow::restore",
"could not open input file "+fname);
// ifstream s(fname.c_str());
// if (!s)
// throw CanteraError("StFlow::restore",
// "could not open input file "+fname);
XML_Node root;
root.build(s);
s.close();
int k;
// XML_Node root;
// root.build(s);
// s.close();
// int k;
const XML_Node* f = root.findID(id);
if (!f) {
throw CanteraError("StFlow::restore","No solution with id = "+id);
}
// const XML_Node* f = root.findID(id);
// if (!f) {
// throw CanteraError("StFlow::restore","No solution with id = "+id);
// }
const XML_Node& flow = f->child("domain");
f = &flow;
// const XML_Node& flow = f->child("domain");
// f = &flow;
//if (f->name() != "flowfield") {
// throw CanteraError("StFlow::restore","The element with id "
// +id+" does not contain flowfield data.");
//}
// //if (f->name() != "flowfield") {
// // throw CanteraError("StFlow::restore","The element with id "
// // +id+" does not contain flowfield data.");
// //}
vector<XML_Node*> str;
f->getChildren("string",str);
int nstr = str.size();
for (int istr = 0; istr < nstr; istr++) {
const XML_Node& nd = *str[istr];
writelog(nd["title"]+": "+nd.value()+"\n");
}
// vector<XML_Node*> str;
// f->getChildren("string",str);
// int nstr = str.size();
// for (int istr = 0; istr < nstr; istr++) {
// const XML_Node& nd = *str[istr];
// writelog(nd["title"]+": "+nd.value()+"\n");
// }
vector<XML_Node*> d;
f->child("grid_data").getChildren("floatArray",d);
int nd = d.size();
// vector<XML_Node*> d;
// f->child("grid_data").getChildren("floatArray",d);
// int nd = d.size();
vector_fp x;
int n, np = 0, j, ks;
string nm;
bool readgrid = false, wrote_header = false;
for (n = 0; n < nd; n++) {
const XML_Node& fa = *d[n];
nm = fa["title"];
if (nm == "z") {
getFloatArray(fa,x,false);
np = x.size();
if (job == -1) {
size_z = np;
//size_soln = (nd - 1)*np;
size_soln = (m_nsp + 4)*np;
return;
}
writelog("Grid contains "+int2str(np)+
" points.\n");
if (size_z < np) {
throw CanteraError("restore",
"grid array must be have length at least "
+int2str(np));
}
// if (size_soln < (m_nsp + 4)*np) {
// vector_fp x;
// int n, np = 0, j, ks;
// string nm;
// bool readgrid = false, wrote_header = false;
// for (n = 0; n < nd; n++) {
// const XML_Node& fa = *d[n];
// nm = fa["title"];
// if (nm == "z") {
// getFloatArray(fa,x,false);
// np = x.size();
// if (job == -1) {
// size_z = np;
// //size_soln = (nd - 1)*np;
// size_soln = (m_nsp + 4)*np;
// return;
// }
// writelog("Grid contains "+int2str(np)+
// " points.\n");
// if (size_z < np) {
// throw CanteraError("restore",
// "solution array must have length at least "
// +int2str((m_nsp + 4)*np));
// "grid array must be have length at least "
// +int2str(np));
// }
copy(x.begin(), x.end(), z);
readgrid = true;
}
}
if (!readgrid) {
throw CanteraError("StFlow::restore",
"solution contains no grid points.");
}
// // if (size_soln < (m_nsp + 4)*np) {
// // throw CanteraError("restore",
// // "solution array must have length at least "
// // +int2str((m_nsp + 4)*np));
// // }
// copy(x.begin(), x.end(), z);
// readgrid = true;
// }
// }
// if (!readgrid) {
// throw CanteraError("StFlow::restore",
// "solution contains no grid points.");
// }
writelog("Importing datasets:\n");
for (n = 0; n < nd; n++) {
const XML_Node& fa = *d[n];
nm = fa["title"];
getFloatArray(fa,x,false);
if (nm == "u") {
writelog("axial velocity ");
if ((int) x.size() == np) {
for (j = 0; j < np; j++) {
soln[index(0,j)] = x[j];
}
}
else {
goto error;
}
}
else if (nm == "z") {
;
}
else if (nm == "V") {
writelog("radial velocity ");
if ((int) x.size() == np) {
for (j = 0; j < np; j++)
soln[index(1,j)] = x[j];
}
else goto error;
}
else if (nm == "T") {
writelog("temperature ");
if ((int) x.size() == np) {
for (j = 0; j < np; j++)
soln[index(2,j)] = x[j];
}
else goto error;
}
else if (nm == "L") {
writelog("lambda ");
if ((int) x.size() == np) {
for (j = 0; j < np; j++)
soln[index(3,j)] = x[j];
}
else goto error;
}
else if (m_thermo->speciesIndex(nm) >= 0) {
writelog(nm+" ");
if ((int) x.size() == np) {
k = m_thermo->speciesIndex(nm);
did_species[k] = 1;
for (j = 0; j < np; j++)
soln[index(k+4,j)] = x[j];
}
}
else
ignored.push_back(nm);
}
// writelog("Importing datasets:\n");
// for (n = 0; n < nd; n++) {
// const XML_Node& fa = *d[n];
// nm = fa["title"];
// getFloatArray(fa,x,false);
// if (nm == "u") {
// writelog("axial velocity ");
// if ((int) x.size() == np) {
// for (j = 0; j < np; j++) {
// soln[index(0,j)] = x[j];
// }
// }
// else {
// goto error;
// }
// }
// else if (nm == "z") {
// ;
// }
// else if (nm == "V") {
// writelog("radial velocity ");
// if ((int) x.size() == np) {
// for (j = 0; j < np; j++)
// soln[index(1,j)] = x[j];
// }
// else goto error;
// }
// else if (nm == "T") {
// writelog("temperature ");
// if ((int) x.size() == np) {
// for (j = 0; j < np; j++)
// soln[index(2,j)] = x[j];
// }
// else goto error;
// }
// else if (nm == "L") {
// writelog("lambda ");
// if ((int) x.size() == np) {
// for (j = 0; j < np; j++)
// soln[index(3,j)] = x[j];
// }
// else goto error;
// }
// else if (m_thermo->speciesIndex(nm) >= 0) {
// writelog(nm+" ");
// if ((int) x.size() == np) {
// k = m_thermo->speciesIndex(nm);
// did_species[k] = 1;
// for (j = 0; j < np; j++)
// soln[index(k+4,j)] = x[j];
// }
// }
// else
// ignored.push_back(nm);
// }
if (ignored.size() != 0) {
writelog("\n\n");
writelog("Ignoring datasets:\n");
int nn = ignored.size();
for (int n = 0; n < nn; n++) {
writelog(ignored[n]+" ");
}
}
// if (ignored.size() != 0) {
// writelog("\n\n");
// writelog("Ignoring datasets:\n");
// int nn = ignored.size();
// for (int n = 0; n < nn; n++) {
// writelog(ignored[n]+" ");
// }
// }
for (ks = 0; ks < nsp; ks++) {
if (did_species[ks] == 0) {
if (!wrote_header) {
writelog("Missing data for species:\n");
wrote_header = true;
}
writelog(m_thermo->speciesName(ks)+" ");
}
}
// for (ks = 0; ks < nsp; ks++) {
// if (did_species[ks] == 0) {
// if (!wrote_header) {
// writelog("Missing data for species:\n");
// wrote_header = true;
// }
// writelog(m_thermo->speciesName(ks)+" ");
// }
// }
writelog("\n\nFinished importing solution.\n\n");
return;
error:
throw CanteraError("StFlow::restore","Data size error");
}
// writelog("\n\nFinished importing solution.\n\n");
// return;
// error:
// throw CanteraError("StFlow::restore","Data size error");
// }
@ -1267,10 +1275,14 @@ namespace Cantera {
for (j = 0; j < np; j++)
soln[index(2,j)] = x[j];
// For fixed-temperature simulations, use the imported temperature profile by default.
// If this is not desired, call setFixedTempProfile *after* restoring the solution.
// For fixed-temperature simulations, use the
// imported temperature profile by default. If
// this is not desired, call setFixedTempProfile
// *after* restoring the solution.
vector_fp zz(np);
for (int jj = 0; jj < np; jj++) zz[jj] = (grid(jj) - zmin())/(zmax() - zmin());
for (int jj = 0; jj < np; jj++)
zz[jj] = (grid(jj) - zmin())/(zmax() - zmin());
setFixedTempProfile(zz, x);
}
else goto error;
@ -1365,20 +1377,16 @@ namespace Cantera {
m_jac = jac;
}
//void StFlow::requestJacUpdate() {
// if (m_jac) m_jac->setAge(10000);
//}
void StFlow::setEnergyFactor(doublereal efctr) {
doublereal de = efctr - m_efctr;
m_efctr = efctr;
int strt = loc();
int jg;
for (int j = 1; j < m_points - 1; j++) {
jg = strt + index(c_offset_T, j);
m_jac->incrementDiagonal(jg, -de);
}
}
// void StFlow::setEnergyFactor(doublereal efctr) {
// doublereal de = efctr - m_efctr;
// m_efctr = efctr;
// int strt = loc();
// int jg;
// for (int j = 1; j < m_points - 1; j++) {
// jg = strt + index(c_offset_T, j);
// m_jac->incrementDiagonal(jg, -de);
// }
// }
}

82
Cantera/src/oneD/StFlow.h
View file

@ -17,7 +17,7 @@
#include "../transport/TransportBase.h"
#include "Domain1D.h"
#include "../Array.h"
#include "../sort.h"
//#include "../sort.h"
#include "../IdealGasPhase.h"
#include "../Kinetics.h"
#include "../funcs.h"
@ -56,10 +56,8 @@ namespace Cantera {
/**
* A class for one-dimensional reacting stagnation-point
* flows. This class implements the one-dimensional similarity
* solution for a chemically-reacting, axisymmetric,
* stagnation-point flow.
* This class implements the one-dimensional similarity solution
* for a chemically-reacting, axisymmetric, flow.
*/
class StFlow : public Domain1D {
@ -82,7 +80,6 @@ namespace Cantera {
virtual void setupGrid(int n, const doublereal* z);
//thermo_t& phase() { return *m_phase; }
thermo_t& phase() { return *m_thermo; }
kinetics_t& kinetics() { return *m_kin; }
@ -90,9 +87,7 @@ namespace Cantera {
* Set the thermo manager. Note that the flow equations assume
* the ideal gas equation.
*/
void setThermo(igthermo_t& th) {
m_thermo = &th;
}
void setThermo(igthermo_t& th) { m_thermo = &th; }
/// set the kinetics manager
void setKinetics(kinetics_t& kin) { m_kin = &kin; }
@ -104,7 +99,7 @@ namespace Cantera {
void setPressure(doublereal p) { m_press = p; }
/// Check that all required parameters have been set.
bool ready();
//bool ready();
virtual void setState(int point, const doublereal* state) {
setTemperature(point, state[2]);
@ -125,26 +120,7 @@ namespace Cantera {
}
}
virtual void _finalize(const doublereal* x) {
int k, j;
doublereal zz, tt;
int nz = m_zfix.size();
bool e = m_do_energy[0];
for (j = 0; j < m_points; j++) {
if (e || nz == 0)
setTemperature(j, T(x, j));
else {
zz = (z(j) - z(0))/(z(m_points - 1) - z(0));
tt = linearInterp(zz, m_zfix, m_tfix);
setTemperature(j, tt);
}
for (k = 0; k < m_nsp; k++) {
setMassFraction(j, k, Y(x, k, j));
}
}
if (e) solveEnergyEqn();
}
virtual void _finalize(const doublereal* x);
void setFixedTempProfile(vector_fp& zfixed, vector_fp& tfixed) {
m_zfix = zfixed;
@ -165,6 +141,7 @@ namespace Cantera {
* Set the mass fraction fixed point for species k at grid
* point j, and disable the species equation so that the
* solution will be held to this value.
* note: in practice, the species are hardly ever held fixed.
*/
void setMassFraction(int j, int k, doublereal y) {
m_fixedy(k,j) = y;
@ -181,23 +158,23 @@ namespace Cantera {
*/
doublereal Y_fixed(int k, int j) const {return m_fixedy(k,j);}
virtual string componentName(int n) const;
/**
* Write a Tecplot zone corresponding to the current solution.
* May be called multiple times to generate animation.
*/
void outputTEC(ostream &s, const doublereal* x,
string title, int zone);
// /**
// * Write a Tecplot zone corresponding to the current solution.
// * May be called multiple times to generate animation.
// */
// void outputTEC(ostream &s, const doublereal* x,
// string title, int zone);
virtual void showSolution(ostream& s, const doublereal* x);
// virtual void showSolution(ostream& s, const doublereal* x);
virtual void showSolution(const doublereal* x);
//void save(string fname, string id, string desc, doublereal* soln);
virtual void save(XML_Node& o, doublereal* sol);
void restore(int job, string fname, string id, int& size_z,
doublereal* z, int& size_soln, doublereal* soln);
// void restore(int job, string fname, string id, int& size_z,
// doublereal* z, int& size_soln, doublereal* soln);
virtual void restore(const XML_Node& dom, doublereal* soln);
@ -240,7 +217,7 @@ namespace Cantera {
needJacUpdate();
}
void setEnergyFactor(doublereal efctr);
// void setEnergyFactor(doublereal efctr);
void fixSpecies(int k=-1) {
if (k == -1) {
@ -259,22 +236,23 @@ namespace Cantera {
virtual void setFixedPoint(int j0, doublereal t0){}
virtual void setBoundaries(FlowBdry::Boundary* left,
FlowBdry::Boundary* right) {
if (left) {
m_boundary[0] = left;
left->faceRight();
}
if (right) {
m_boundary[1] = right;
right->faceLeft();
}
}
// virtual void setBoundaries(FlowBdry::Boundary* left,
// FlowBdry::Boundary* right) {
// if (left) {
// m_boundary[0] = left;
// left->faceRight();
// }
// if (right) {
// m_boundary[1] = right;
// right->faceLeft();
// }
// }
void setJac(MultiJac* jac);
void setGas(const doublereal* x,int j);
void setGasAtMidpoint(const doublereal* x,int j);
protected:
// used to write mass fractions to plot files.

28
Cantera/src/units.h
View file

@ -1,3 +1,9 @@
/**
* @file units.h
*
* Unit conversions.
*/
#ifndef CT_UNITS_H
#define CT_UNITS_H
@ -20,9 +26,13 @@ namespace Cantera {
}
}
virtual ~Unit() {
}
virtual ~Unit() {}
/**
* Return the multiplier required to convert an activation
* energy to SI units.
* @param units activation energy units
*/
doublereal actEnergyToSI(string units) {
if (m_act_u.find(units) != m_act_u.end()) {
return m_act_u[units];
@ -32,15 +42,25 @@ namespace Cantera {
}
}
/**
* Return the multiplier required to convert a dimensional quantity
* with units specified by string 'units' to SI units.
*/
doublereal toSI(string units) {
// if dimensionless, return 1.0
if (units == "") return 1.0;
doublereal f = 1.0, fctr;
int tsize;
string u = units, tok, tsub;
int k;
char action = '-';
//if (units[0] == '/') action = '/';
while (1 > 0) {
// get token consisting of all characters up to the next
// dash, slash, or the end of the string
k = u.find_first_of("/-");
if (k >= 0)
tok = u.substr(0,k);
@ -79,6 +99,8 @@ namespace Cantera {
fctr = m_u[tok];
}
// tok is not one of the entries in map m_u, then
// m_u[tok] returns 0.0. Check for this.
if (fctr == 0)
throw CanteraError("toSI","unknown unit: "+tsub);
if (action == '-') f *= fctr;

113
Cantera/src/updaters.h
View file

@ -1,113 +0,0 @@
#ifndef CT_UPDATERS_H
#define CT_UPDATERS_H
#include "PropertyUpdater.h"
namespace Cantera {
//--------------------------------------------------------
// Property Updaters
//--------------------------------------------------------
/**
* Invokes method 'update_T' of the object it is initialized with.
*/
template<class S>
struct T_Updater : public Updater {
T_Updater(S& s) : m_s(s) {}
void update() { m_s.update_T(); }
S& m_s;
};
/**
* Invokes method 'updateMoleFractions' of the object it is
* initialized with.
*/
template<class S>
struct UpdateMoleFractions : public Updater {
UpdateMoleFractions(S& s) : Updater(), m_s(s) {}
void update() { m_s.updateMoleFractions(); }
S& m_s;
};
/**
* Invokes method 'updateMW' of the object it is
* initialized with.
*/
template<class S>
struct UpdateMolWt : public Updater {
UpdateMolWt(S& s) : Updater(), m_s(s) {}
void update() { m_s.updateMW(); }
S& m_s;
};
/**
* Updater responsible for updating the species standard-state
* thermodynamic properties.
* @ingroup updategroup
*/
template<class S>
struct UpdateThermo : public Updater {
UpdateThermo(S& s) : Updater(), m_s(s) {}
void update() { m_s._updateThermo(); }
S& m_s;
};
/**
* Updater responsible for updating the temperature-dependent
* parts of the transport properties.
* @ingroup updategroup
*/
template<class S>
struct UpdateTransport_T : public Updater {
UpdateTransport_T(S& s) : Updater(), m_s(s) {}
void update() { m_s._update_transport_T(); }
S& m_s;
};
/**
* Updater responsible for updating the concentration-dependent
* parts of the transport properties.
* @ingroup updategroup
*/
template<class S>
struct UpdateTransport_C : public Updater {
UpdateTransport_C(S& s) : Updater(), m_s(s) {}
void update() { m_s._update_transport_C(); }
S& m_s;
};
/**
* Template for an updater subclass that calls method _updateRates_T()
* of the object it is initialized with.
* @ingroup updategroup
*/
template<class S>
struct UpdateRates_T : public Updater {
UpdateRates_T(S& s) : Updater(), m_s(s) {}
void update() { m_s._update_rates_T(); }
S& m_s;
};
/**
* Template for an updater subclass that calls method _updateRates_C()
* of the object it is initialized with.
* @ingroup updategroup
*/
template<class S>
struct UpdateRates_C : public Updater {
UpdateRates_C(S& s) : Updater(), m_s(s) {}
void update() { m_s._update_rates_C(); }
S& m_s;
};
}
#endif

16
include/Interface.h
View file

@ -15,19 +15,13 @@ namespace Cantera {
public:
Interface(string infile, string id, vector<ThermoPhase*> phases)
: m_ok(false), m_r(0) {
string path = findInputFile(infile);
ifstream fin(path.c_str());
if (!fin) {
throw CanteraError("Interface","could not open "
+path+" for reading.");
}
m_r = new XML_Node("-");
m_r->build(fin);
XML_Node* x = find_XML("", m_r, id, "", "");
m_r = get_XML_File(infile);
if (id == "-") id = "";
XML_Node* x = get_XML_Node("#"+id, m_r);
if (!x)
throw CanteraError("Interface","error in find_XML");
throw CanteraError("Interface","error in get_XML_Node");
importPhase(*x, this);
phases.push_back(this);