This commit is contained in:
Dave Goodwin
2003-11-17 13:51:16 +00:00
parent 7225d42000
commit bc7f8dd581
5 changed files with 51 additions and 540 deletions

View File

@@ -43,7 +43,8 @@ else
@PYTHON_CMD@ setup_winmatlab.py
endif
(@MATLAB_CMD@ -nodesktop -nojvm -nosplash -r setup)
rm -f setup.m
rm -f setup.m
install:
@INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera

View File

@@ -4,10 +4,7 @@ import exceptions
class Transport:
"""Transport property manager. In most cases, this class is used
as a base class to provide transport properties, and is not
instantiated directly.
"""Transport property manager.
A transport property manager is responsible for computing transport
properties.
"""
@@ -16,9 +13,10 @@ class Transport:
"""Create a transport property manager.
xml_phase --- XML phase element
phase --- ThermoPhase instance representing phase transport
properties are for
model --- transport model
phase --- ThermoPhase instance representing the phase that the
transport properties are for
model --- string specifying transport model. If omitted,
model will be taken from the input file.
loglevel --- controls amount of diagnostic output
"""
@@ -38,10 +36,11 @@ class Transport:
self._models[self.model] = self.__tr_id
def __del__(self):
try:
_cantera.tran_delete(self.__tr_id)
except:
pass
for m in self._models.keys():
try:
_cantera.tran_delete(self._models[m])
except:
pass
def addTransportModel(self, model, loglevel=1):
new_id = _cantera.Transport(model,

View File

@@ -26,7 +26,7 @@ tol_ss = [1.0e-5, 1.0e-9] # [rtol atol] for steady-state
# problem
tol_ts = [1.0e-5, 1.0e-4] # [rtol atol] for time stepping
loglevel = 5 # amount of diagnostic output (0
loglevel = 1 # amount of diagnostic output (0
# to 5)
refine_grid = 1 # 1 to enable refinement, 0 to

View File

@@ -595,7 +595,7 @@ namespace Cantera {
*/
void AxiStagnFlow::updateTransport(doublereal* x,int j0, int j1) {
int j,k,m;
// writelog("\nentered updateTransport\n");
if (m_transport_option == c_Mixav_Transport) {
for (j = j0; j < j1; j++) {
setGasAtMidpoint(x,j);
@@ -610,7 +610,6 @@ namespace Cantera {
for (m = j0; m < j1; m++) {
setGasAtMidpoint(x,m);
dz = m_z[m+1] - m_z[m];
//dz = m_z[j+1] - m_z[j];
wtm = m_thermo->meanMolecularWeight();
m_visc[m] = m_trans->viscosity();
@@ -637,246 +636,9 @@ namespace Cantera {
}
}
}
//writelog("leaving updateTransport\n");
}
// void OneDFlow::eval(int jg, doublereal* xg, doublereal* rg, integer* diagg,
// doublereal rdt) {
// static double elapsed;
// clock_t t0 = clock();
// if (jg >= 0) rdt = 0.0;
// 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();
// int jmin, jmax, jpt;
// jpt = jg - firstPoint();
// // 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);
// int j, k;
// //-----------------------------------------------------
// // 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);
// // 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);
// //----------------------------------------------------
// // evaluate the residual equations at all required
// // grid points
// //----------------------------------------------------
// doublereal sum, sum2, deltaz, dtdzj;
// for (j = jmin; j <= jmax; j++) {
// //----------------------------------------------
// // boundaries
// //----------------------------------------------
// 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));
// }
// //------------------------------------------
// // interior points
// //------------------------------------------
// else {
// // 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);
// // species equations
// getWdot(x,j);
// 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));
// }
// // energy equation
// 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));
// 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);
// }
// 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;
// }
// /**
// * 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];
// // 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;
// }
/**
* Print the solution.
*/
@@ -993,35 +755,6 @@ 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;
// 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 {
switch(n) {
case 0: return "u";
@@ -1038,199 +771,6 @@ namespace Cantera {
}
// /**
// * 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) {
// 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);
// 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& flow = f->child("domain");
// f = &flow;
// //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*> 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) {
// // 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);
// }
// 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)+" ");
// }
// }
// writelog("\n\nFinished importing solution.\n\n");
// return;
// error:
// throw CanteraError("StFlow::restore","Data size error");
// }
void StFlow::restore(const XML_Node& dom, doublereal* soln) {
vector<string> ignored;
@@ -1413,16 +953,5 @@ namespace Cantera {
m_jac = jac;
}
// 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);
// }
// }
}
} // namespace

View File

@@ -17,7 +17,6 @@
#include "../transport/TransportBase.h"
#include "Domain1D.h"
#include "../Array.h"
//#include "../sort.h"
#include "../IdealGasPhase.h"
#include "../Kinetics.h"
#include "../funcs.h"
@@ -25,18 +24,17 @@
namespace Cantera {
typedef IdealGasPhase igthermo_t;
class MultiJac;
//------------------------------------------
// constants
//------------------------------------------
/**
* Offsets of solution components in the solution array.
*/
// Offsets of solution components in the solution array.
const unsigned int c_offset_U = 0; // axial velocity
const unsigned int c_offset_V = 1; // strain rate
const unsigned int c_offset_T = 2; // temperature
@@ -48,16 +46,17 @@ namespace Cantera {
const int c_Multi_Transport = 1;
const int c_Soret = 2;
//-----------------------------------------------------------
// Class StFlow
//-----------------------------------------------------------
/**
* This class implements the one-dimensional similarity solution
* for a chemically-reacting, axisymmetric, flow.
* This class represents 1D flow domains that satisfy the
* one-dimensional similarity solution for chemically-reacting,
* axisymmetric, flows.
*/
class StFlow : public Domain1D {
@@ -67,7 +66,11 @@ namespace Cantera {
// construction and destruction
//--------------------------------
// Constructor.
/// Constructor. Create a new flow domain.
/// @param gas Object representing the gas phase. This object
/// will be used to evaluate all thermodynamic, kinetic, and transport
/// properties.
/// @param nsp Number of species.
StFlow(igthermo_t* ph = 0, int nsp = 1, int points = 1);
/// Destructor.
@@ -89,18 +92,19 @@ namespace Cantera {
*/
void setThermo(igthermo_t& th) { m_thermo = &th; }
/// set the kinetics manager
/// Set the kinetics manager. The kinetics manager must
void setKinetics(kinetics_t& kin) { m_kin = &kin; }
/// set the transport manager
void setTransport(Transport& trans, bool withSoret = false);
/// set the pressure
/// Set the pressure. Since the flow equations are for the limit of
/// small Mach number, the pressure is very nearly constant
/// throughout the flow.
void setPressure(doublereal p) { m_press = p; }
/// Check that all required parameters have been set.
//bool ready();
/// @todo remove? may be unused
virtual void setState(int point, const doublereal* state) {
setTemperature(point, state[2]);
int k;
@@ -109,7 +113,9 @@ namespace Cantera {
}
}
/// Write the initial solution estimate into
/// array x.
virtual void _getInitialSoln(doublereal* x) {
int k, j;
for (j = 0; j < m_points; j++) {
@@ -121,7 +127,12 @@ namespace Cantera {
}
virtual void _finalize(const doublereal* x);
/// Sometimes it is desired to carry out the simulation
/// using a specified temperature profile, rather than
/// computing it by solving the energy equation. This
/// method specifies this profile.
void setFixedTempProfile(vector_fp& zfixed, vector_fp& tfixed) {
m_zfix = zfixed;
m_tfix = tfixed;
@@ -147,35 +158,23 @@ namespace Cantera {
m_fixedy(k,j) = y;
m_do_species[k] = true; // false;
}
/**
* The fixed temperature value at point j.
*/
/// The fixed temperature value at point j.
doublereal T_fixed(int j) const {return m_fixedtemp[j];}
/**
* The fixed mass fraction value of species k at point j.
*/
/// The fixed mass fraction value of species k at point j.
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);
// virtual void showSolution(ostream& s, const doublereal* x);
virtual void showSolution(const doublereal* x);
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);
virtual void restore(const XML_Node& dom, doublereal* soln);
// overloaded in subclasses
@@ -217,8 +216,6 @@ namespace Cantera {
needJacUpdate();
}
// void setEnergyFactor(doublereal efctr);
void fixSpecies(int k=-1) {
if (k == -1) {
for (int i = 0; i < m_nsp; i++)
@@ -226,7 +223,6 @@ namespace Cantera {
}
else m_do_species[k] = false;
needJacUpdate();
//m_jac->setAge(10000);
}
void integrateChem(doublereal* x,doublereal dt);
@@ -236,18 +232,6 @@ 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();
// }
// }
void setJac(MultiJac* jac);
void setGas(const doublereal* x,int j);
void setGasAtMidpoint(const doublereal* x,int j);
@@ -255,8 +239,6 @@ namespace Cantera {
protected:
// used to write mass fractions to plot files.
doublereal component(const doublereal* x, int i, int j) const {
doublereal xx = x[index(i,j)];
return xx;