diff --git a/Cantera/src/Surf1D.h b/Cantera/src/Surf1D.h deleted file mode 100755 index a2feb1400..000000000 --- a/Cantera/src/Surf1D.h +++ /dev/null @@ -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 m_bulk; - vector m_nbulk; - int m_nsurf; - vector_fp m_mult; - vector m_do_surf_species; - vector_fp m_fixed_cov; - }; - -} - -#endif diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index 78268e89a..1b0f3b770 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -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 ""; } } - /** - * 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 ignored; - int nsp = m_thermo->nSpecies(); - vector_int did_species(nsp, 0); +// 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); +// 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 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 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 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); +// } +// } } diff --git a/Cantera/src/oneD/StFlow.h b/Cantera/src/oneD/StFlow.h index 86461ac77..e5637aca5 100644 --- a/Cantera/src/oneD/StFlow.h +++ b/Cantera/src/oneD/StFlow.h @@ -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. diff --git a/Cantera/src/units.h b/Cantera/src/units.h index f0ea6a54f..1ddd437a1 100644 --- a/Cantera/src/units.h +++ b/Cantera/src/units.h @@ -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; diff --git a/Cantera/src/updaters.h b/Cantera/src/updaters.h deleted file mode 100755 index 7985645e7..000000000 --- a/Cantera/src/updaters.h +++ /dev/null @@ -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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - struct UpdateRates_C : public Updater { - UpdateRates_C(S& s) : Updater(), m_s(s) {} - void update() { m_s._update_rates_C(); } - S& m_s; - }; - - -} - -#endif - - diff --git a/include/Interface.h b/include/Interface.h index a79dd647b..2f0414ec7 100644 --- a/include/Interface.h +++ b/include/Interface.h @@ -15,19 +15,13 @@ namespace Cantera { public: Interface(string infile, string id, vector 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);