From bc7f8dd5817c022bb53a0f6c71a2f885657a1aa3 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 17 Nov 2003 13:51:16 +0000 Subject: [PATCH] cleanup --- Cantera/matlab/Makefile.in | 3 +- Cantera/python/Cantera/Transport.py | 21 +- Cantera/python/examples/flame2.py | 2 +- Cantera/src/oneD/StFlow.cpp | 475 +--------------------------- Cantera/src/oneD/StFlow.h | 90 +++--- 5 files changed, 51 insertions(+), 540 deletions(-) diff --git a/Cantera/matlab/Makefile.in b/Cantera/matlab/Makefile.in index 67dc84aef..34500356e 100644 --- a/Cantera/matlab/Makefile.in +++ b/Cantera/matlab/Makefile.in @@ -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 diff --git a/Cantera/python/Cantera/Transport.py b/Cantera/python/Cantera/Transport.py index def692d14..9e46b1943 100755 --- a/Cantera/python/Cantera/Transport.py +++ b/Cantera/python/Cantera/Transport.py @@ -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, diff --git a/Cantera/python/examples/flame2.py b/Cantera/python/examples/flame2.py index 8917f8383..7b45893b3 100755 --- a/Cantera/python/examples/flame2.py +++ b/Cantera/python/examples/flame2.py @@ -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 diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index b82338e2e..c130e7c11 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -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 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 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 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 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 diff --git a/Cantera/src/oneD/StFlow.h b/Cantera/src/oneD/StFlow.h index 9f544c2b2..a1c06ce5e 100644 --- a/Cantera/src/oneD/StFlow.h +++ b/Cantera/src/oneD/StFlow.h @@ -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;