diff --git a/Cantera/src/base/ctml.cpp b/Cantera/src/base/ctml.cpp index d11ad58a6..b705927c0 100755 --- a/Cantera/src/base/ctml.cpp +++ b/Cantera/src/base/ctml.cpp @@ -883,17 +883,27 @@ namespace ctml { * units converter is used. * @param nodeName XML Name of the XML node to read. * The default value for the node name is floatArray + * + * @return Returns the number of floats read */ - void getFloatArray(const Cantera::XML_Node& node, vector_fp& v, const bool convert, + int getFloatArray(const Cantera::XML_Node& node, vector_fp& v, const bool convert, const std::string unitsString, const std::string nodeName) { string::size_type icom; string numstr; doublereal dtmp; string nn = node.name(); - if (nn != nodeName) - throw CanteraError("getFloatArray", - "wrong xml element type/name: was expecting " + const Cantera::XML_Node *readNode = &node; + if (nn != nodeName) { + vector ll; + node.getChildren(nodeName, ll); + if (ll.size() == 0) { + throw CanteraError("getFloatArray", + "wrong xml element type/name: was expecting " + nodeName + "but accessed " + node.name()); + } else { + readNode = ll[0]; + } + } v.clear(); doublereal vmin = Undef, vmax = Undef; @@ -902,7 +912,7 @@ namespace ctml { /* * Get the attributes field, units, from the XML node */ - std::string units = node["units"]; + std::string units = (*readNode)["units"]; if (units != "" && convert) { if (unitsString == "actEnergy" && units != "") { funit = actEnergyToSI(units); @@ -911,13 +921,13 @@ namespace ctml { } } - if (node["min"] != "") - vmin = atofCheck(node["min"].c_str()); - if (node["max"] != "") - vmax = atofCheck(node["max"].c_str()); + if ((*readNode)["min"] != "") + vmin = atofCheck((*readNode)["min"].c_str()); + if ((*readNode)["max"] != "") + vmax = atofCheck((*readNode)["max"].c_str()); doublereal vv; - std::string val = node.value(); + std::string val = readNode->value(); while (1 > 0) { icom = val.find(','); if (icom != string::npos) { @@ -956,6 +966,7 @@ namespace ctml { for (int n = 0; n < nv; n++) { v[n] *= funit; } + return v.size(); } // This routine is used to interpret the value portions of XML diff --git a/Cantera/src/base/ctml.h b/Cantera/src/base/ctml.h index 4b3d4ebd9..25fbe3a9c 100755 --- a/Cantera/src/base/ctml.h +++ b/Cantera/src/base/ctml.h @@ -363,8 +363,9 @@ namespace ctml { * units converter is used. * @param nodeName XML Name of the XML node to read. * The default value for the node name is floatArray + * @return Returns the number of floats read into v. */ - void getFloatArray(const Cantera::XML_Node& node, Cantera::vector_fp& v, + int getFloatArray(const Cantera::XML_Node& node, Cantera::vector_fp& v, const bool convert=true, const std::string typeString="", const std::string nodeName = "floatArray"); diff --git a/Cantera/src/oneD/Domain1D.h b/Cantera/src/oneD/Domain1D.h index b201fe304..7d57583ec 100644 --- a/Cantera/src/oneD/Domain1D.h +++ b/Cantera/src/oneD/Domain1D.h @@ -58,7 +58,7 @@ namespace Cantera { m_jstart(0), m_left(0), m_right(0), - m_id("-"), m_desc("-"), + m_id(""), m_desc(""), m_refiner(0), m_bw(-1) { resize(nv, points); } @@ -412,9 +412,11 @@ namespace Cantera { */ void setID(const std::string& s) {m_id = s;} - std::string id() { - if (m_id != "") return m_id; - else return std::string("domain ") + int2str(m_index); + std::string id() const { + if (m_id != "") + return m_id; + else + return std::string("domain ") + int2str(m_index); } /** @@ -519,7 +521,10 @@ namespace Cantera { int m_jstart; Domain1D *m_left, *m_right; - std::string m_id, m_desc; + + //! Identity tag for the domain + std::string m_id; + std::string m_desc; Refiner* m_refiner; vector_int m_td; std::vector m_name; diff --git a/Cantera/src/oneD/StFlow.cpp b/Cantera/src/oneD/StFlow.cpp index 826cfd3f7..29fc73328 100644 --- a/Cantera/src/oneD/StFlow.cpp +++ b/Cantera/src/oneD/StFlow.cpp @@ -32,1170 +32,1170 @@ using namespace std; namespace Cantera { - //------------------- importSolution ------------------------ + //------------------- importSolution ------------------------ - /** - * Import a previous solution to use as an initial estimate. The - * previous solution may have been computed using a different - * reaction mechanism. Species in the old and new mechanisms are - * matched by name, and any species in the new mechanism that were - * not in the old one are set to zero. The new solution is created - * with the same number of grid points as in the old solution. - */ - void importSolution(int points, - doublereal* oldSoln, igthermo_t& oldmech, - int size_new, doublereal* newSoln, igthermo_t& newmech) { + /** + * Import a previous solution to use as an initial estimate. The + * previous solution may have been computed using a different + * reaction mechanism. Species in the old and new mechanisms are + * matched by name, and any species in the new mechanism that were + * not in the old one are set to zero. The new solution is created + * with the same number of grid points as in the old solution. + */ + void importSolution(int points, + doublereal* oldSoln, igthermo_t& oldmech, + int size_new, doublereal* newSoln, igthermo_t& newmech) { - // Number of components in old and new solutions - int nv_old = oldmech.nSpecies() + 4; - int nv_new = newmech.nSpecies() + 4; + // Number of components in old and new solutions + int nv_old = oldmech.nSpecies() + 4; + int nv_new = newmech.nSpecies() + 4; - if (size_new < nv_new*points) { - throw CanteraError("importSolution", - "new solution array must have length "+ - int2str(nv_new*points)); - } + if (size_new < nv_new*points) { + throw CanteraError("importSolution", + "new solution array must have length "+ + int2str(nv_new*points)); + } - int n, j, knew; - string nm; + int n, j, knew; + string nm; - // copy u,V,T,lambda - for (j = 0; j < points; j++) - for (n = 0; n < 4; n++) - newSoln[nv_new*j + n] = oldSoln[nv_old*j + n]; + // copy u,V,T,lambda + for (j = 0; j < points; j++) + for (n = 0; n < 4; n++) + newSoln[nv_new*j + n] = oldSoln[nv_old*j + n]; - // copy mass fractions - int nsp0 = oldmech.nSpecies(); - //int nsp1 = newmech.nSpecies(); + // copy mass fractions + int nsp0 = oldmech.nSpecies(); + //int nsp1 = newmech.nSpecies(); - // loop over the species in the old mechanism - for (int k = 0; k < nsp0; k++) { - nm = oldmech.speciesName(k); // name + // loop over the species in the old mechanism + for (int k = 0; k < nsp0; k++) { + nm = oldmech.speciesName(k); // name - // location of this species in the new mechanism. - // If < 0, then the species is not in the new mechanism. - knew = newmech.speciesIndex(nm); + // location of this species in the new mechanism. + // If < 0, then the species is not in the new mechanism. + knew = newmech.speciesIndex(nm); - // copy this species from the old to the new solution vectors - if (knew >= 0) { - for (j = 0; j < points; j++) { - newSoln[nv_new*j + 4 + knew] = oldSoln[nv_old*j + 4 + k]; - } - } - } - - - // normalize mass fractions - for (j = 0; j < points; j++) { - newmech.setMassFractions(&newSoln[nv_new*j + 4]); - newmech.getMassFractions(&newSoln[nv_new*j + 4]); - } + // copy this species from the old to the new solution vectors + if (knew >= 0) { + for (j = 0; j < points; j++) { + newSoln[nv_new*j + 4 + knew] = oldSoln[nv_old*j + 4 + k]; + } + } } - static void st_drawline() { - writelog("\n-------------------------------------" - "------------------------------------------"); + // normalize mass fractions + for (j = 0; j < points; j++) { + newmech.setMassFractions(&newSoln[nv_new*j + 4]); + newmech.getMassFractions(&newSoln[nv_new*j + 4]); + } + } + + + static void st_drawline() { + writelog("\n-------------------------------------" + "------------------------------------------"); + } + + StFlow::StFlow(igthermo_t* ph, int nsp, int points) : + Domain1D(nsp+4, points), + m_inlet_u(0.0), + m_inlet_V(0.0), + m_inlet_T(-1.0), + m_surface_T(-1.0), + m_press(-1.0), + m_nsp(nsp), + m_thermo(0), + m_kin(0), + m_trans(0), + m_jac(0), + m_ok(false), + m_do_soret(false), + m_transport_option(-1), + m_efctr(0.0) + { + m_type = cFlowType; + + m_points = points; + m_thermo = ph; + + if (ph == 0) return; // used to create a dummy object + + int nsp2 = m_thermo->nSpecies(); + if (nsp2 != m_nsp) { + m_nsp = nsp2; + Domain1D::resize(m_nsp+4, points); } - StFlow::StFlow(igthermo_t* ph, int nsp, int points) : - Domain1D(nsp+4, points), - m_inlet_u(0.0), - m_inlet_V(0.0), - m_inlet_T(-1.0), - m_surface_T(-1.0), - m_press(-1.0), - m_nsp(nsp), - m_thermo(0), - m_kin(0), - m_trans(0), - m_jac(0), - m_ok(false), - m_do_soret(false), - m_transport_option(-1), - m_efctr(0.0) - { - m_type = cFlowType; - m_points = points; - m_thermo = ph; + // make a local copy of the species molecular weight vector + m_wt = m_thermo->molecularWeights(); - if (ph == 0) return; // used to create a dummy object + // the species mass fractions are the last components in the solution + // vector, so the total number of components is the number of species + // plus the offset of the first mass fraction. + m_nv = c_offset_Y + m_nsp; - int nsp2 = m_thermo->nSpecies(); - if (nsp2 != m_nsp) { - m_nsp = nsp2; - Domain1D::resize(m_nsp+4, points); - } + // enable all species equations by default + m_do_species.resize(m_nsp, true); + + // but turn off the energy equation at all points + m_do_energy.resize(m_points,false); + + m_diff.resize(m_nsp*m_points); + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_surfdot.resize(m_nsp, 0.0); + m_ybar.resize(m_nsp); - // make a local copy of the species molecular weight vector - m_wt = m_thermo->molecularWeights(); + //-------------- default solution bounds -------------------- - // the species mass fractions are the last components in the solution - // vector, so the total number of components is the number of species - // plus the offset of the first mass fraction. - m_nv = c_offset_Y + m_nsp; - - // enable all species equations by default - m_do_species.resize(m_nsp, true); - - // but turn off the energy equation at all points - m_do_energy.resize(m_points,false); - - m_diff.resize(m_nsp*m_points); - m_multidiff.resize(m_nsp*m_nsp*m_points); - m_flux.resize(m_nsp,m_points); - m_wdot.resize(m_nsp,m_points, 0.0); - m_surfdot.resize(m_nsp, 0.0); - m_ybar.resize(m_nsp); - - - //-------------- default solution bounds -------------------- - - vector_fp vmin(m_nv), vmax(m_nv); + vector_fp vmin(m_nv), vmax(m_nv); - // no bounds on u - vmin[0] = -1.e20; - vmax[0] = 1.e20; + // no bounds on u + vmin[0] = -1.e20; + vmax[0] = 1.e20; - // V - vmin[1] = -1.e20; - vmax[1] = 1.e20; + // V + vmin[1] = -1.e20; + vmax[1] = 1.e20; - // temperature bounds - vmin[2] = 200.0; - vmax[2]= 1.e9; + // temperature bounds + vmin[2] = 200.0; + vmax[2]= 1.e9; - // lamda should be negative - vmin[3] = -1.e20; - vmax[3] = 1.e20; + // lamda should be negative + vmin[3] = -1.e20; + vmax[3] = 1.e20; - // mass fraction bounds - int k; - for (k = 0; k < m_nsp; k++) { - vmin[4+k] = -1.0e-5; - vmax[4+k] = 1.0e5; - } - setBounds(vmin.size(), DATA_PTR(vmin), vmax.size(), DATA_PTR(vmax)); - - - //-------------------- default error tolerances ---------------- - vector_fp rtol(m_nv, 1.0e-8); - vector_fp atol(m_nv, 1.0e-15); - setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),false); - setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),true); - - //-------------------- grid refinement ------------------------- - m_refiner->setActive(0, false); - m_refiner->setActive(1, false); - m_refiner->setActive(2, false); - m_refiner->setActive(3, false); - - vector_fp gr; - for (int ng = 0; ng < m_points; ng++) gr.push_back(1.0*ng/m_points); - setupGrid(m_points, DATA_PTR(gr)); - setID("stagnation flow"); + // mass fraction bounds + int k; + for (k = 0; k < m_nsp; k++) { + vmin[4+k] = -1.0e-5; + vmax[4+k] = 1.0e5; } + setBounds(vmin.size(), DATA_PTR(vmin), vmax.size(), DATA_PTR(vmax)); - /** - * Change the grid size. Called after grid refinement. - */ - void StFlow::resize(int ncomponents, int points) { - Domain1D::resize(ncomponents, points); - m_rho.resize(m_points, 0.0); - m_wtm.resize(m_points, 0.0); - m_cp.resize(m_points, 0.0); - m_enth.resize(m_points, 0.0); - m_visc.resize(m_points, 0.0); - m_tcon.resize(m_points, 0.0); + //-------------------- default error tolerances ---------------- + vector_fp rtol(m_nv, 1.0e-8); + vector_fp atol(m_nv, 1.0e-15); + setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),false); + setTolerances(rtol.size(), DATA_PTR(rtol), atol.size(), DATA_PTR(atol),true); - if (m_transport_option == c_Mixav_Transport) { - m_diff.resize(m_nsp*m_points); - } - else { - m_multidiff.resize(m_nsp*m_nsp*m_points); - m_diff.resize(m_nsp*m_points); - } - m_flux.resize(m_nsp,m_points); - m_wdot.resize(m_nsp,m_points, 0.0); - m_do_energy.resize(m_points,false); + //-------------------- grid refinement ------------------------- + m_refiner->setActive(0, false); + m_refiner->setActive(1, false); + m_refiner->setActive(2, false); + m_refiner->setActive(3, false); - m_fixedy.resize(m_nsp, m_points); - m_fixedtemp.resize(m_points); + vector_fp gr; + for (int ng = 0; ng < m_points; ng++) gr.push_back(1.0*ng/m_points); + setupGrid(m_points, DATA_PTR(gr)); + setID("stagnation flow"); + } - m_dz.resize(m_points-1); - m_z.resize(m_points); - } + + /** + * Change the grid size. Called after grid refinement. + */ + void StFlow::resize(int ncomponents, int points) { + Domain1D::resize(ncomponents, points); + m_rho.resize(m_points, 0.0); + m_wtm.resize(m_points, 0.0); + m_cp.resize(m_points, 0.0); + m_enth.resize(m_points, 0.0); + m_visc.resize(m_points, 0.0); + m_tcon.resize(m_points, 0.0); + + if (m_transport_option == c_Mixav_Transport) { + m_diff.resize(m_nsp*m_points); + } + else { + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_diff.resize(m_nsp*m_points); + } + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_do_energy.resize(m_points,false); + + m_fixedy.resize(m_nsp, m_points); + m_fixedtemp.resize(m_points); + + m_dz.resize(m_points-1); + m_z.resize(m_points); + } - void StFlow::setupGrid(int n, const doublereal* z) { - resize(m_nv, n); - int j; + void StFlow::setupGrid(int n, const doublereal* z) { + resize(m_nv, n); + int j; - m_z[0] = z[0]; - for (j = 1; j < m_points; j++) { - m_z[j] = z[j]; - m_dz[j-1] = m_z[j] - m_z[j-1]; - } + m_z[0] = z[0]; + for (j = 1; j < m_points; j++) { + m_z[j] = z[j]; + m_dz[j-1] = m_z[j] - m_z[j-1]; } + } - /** - * Install a transport manager. - */ - void StFlow::setTransport(Transport& trans, bool withSoret) { - m_trans = &trans; - m_do_soret = withSoret; + /** + * Install a transport manager. + */ + void StFlow::setTransport(Transport& trans, bool withSoret) { + m_trans = &trans; + m_do_soret = withSoret; - if (m_trans->model() == cMulticomponent) { - m_transport_option = c_Multi_Transport; - m_multidiff.resize(m_nsp*m_nsp*m_points); - m_diff.resize(m_nsp*m_points); - m_dthermal.resize(m_nsp, m_points, 0.0); - } - else if (m_trans->model() == cMixtureAveraged) { - m_transport_option = c_Mixav_Transport; - m_diff.resize(m_nsp*m_points); - if (withSoret) - throw CanteraError("setTransport", - "Thermal diffusion (the Soret effect) " - "requires using a multicomponent transport model."); - } - else - throw CanteraError("setTransport","unknown transport model."); + if (m_trans->model() == cMulticomponent) { + m_transport_option = c_Multi_Transport; + m_multidiff.resize(m_nsp*m_nsp*m_points); + m_diff.resize(m_nsp*m_points); + m_dthermal.resize(m_nsp, m_points, 0.0); } - - void StFlow::enableSoret(bool withSoret) { - if (m_transport_option == c_Multi_Transport) - m_do_soret = withSoret; - else { - throw CanteraError("setTransport", - "Thermal diffusion (the Soret effect) " - "requires using a multicomponent transport model."); - } + else if (m_trans->model() == cMixtureAveraged) { + m_transport_option = c_Mixav_Transport; + m_diff.resize(m_nsp*m_points); + if (withSoret) + throw CanteraError("setTransport", + "Thermal diffusion (the Soret effect) " + "requires using a multicomponent transport model."); } + else + throw CanteraError("setTransport","unknown transport model."); + } - - /** - * Set the gas object state to be consistent with the solution at - * point j. - */ - void StFlow::setGas(const doublereal* x,int j) { - m_thermo->setTemperature(T(x,j)); - const doublereal* yy = x + m_nv*j + c_offset_Y; - m_thermo->setMassFractions_NoNorm(yy); - m_thermo->setPressure(m_press); + void StFlow::enableSoret(bool withSoret) { + if (m_transport_option == c_Multi_Transport) + m_do_soret = withSoret; + else { + throw CanteraError("setTransport", + "Thermal diffusion (the Soret effect) " + "requires using a multicomponent transport model."); } + } - /** - * Set the gas state to be consistent with the solution at the - * midpoint between j and j + 1. - */ - void StFlow::setGasAtMidpoint(const doublereal* x,int j) { - m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); - const doublereal* yyj = x + m_nv*j + c_offset_Y; - const doublereal* yyjp = x + m_nv*(j+1) + c_offset_Y; - for (int k = 0; k < m_nsp; k++) - m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); - m_thermo->setMassFractions_NoNorm(DATA_PTR(m_ybar)); - m_thermo->setPressure(m_press); + /** + * Set the gas object state to be consistent with the solution at + * point j. + */ + void StFlow::setGas(const doublereal* x,int j) { + m_thermo->setTemperature(T(x,j)); + const doublereal* yy = x + m_nv*j + c_offset_Y; + m_thermo->setMassFractions_NoNorm(yy); + m_thermo->setPressure(m_press); + } + + + /** + * Set the gas state to be consistent with the solution at the + * midpoint between j and j + 1. + */ + void StFlow::setGasAtMidpoint(const doublereal* x,int j) { + m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); + const doublereal* yyj = x + m_nv*j + c_offset_Y; + const doublereal* yyjp = x + m_nv*(j+1) + c_offset_Y; + for (int k = 0; k < m_nsp; k++) + m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); + m_thermo->setMassFractions_NoNorm(DATA_PTR(m_ybar)); + m_thermo->setPressure(m_press); + } + + + 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(); + } - 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 + * flow. If jpt is less than zero, the residual function is + * evaluated at all grid points. If jpt >= 0, then the residual + * function is only evaluated at grid points jpt-1, jpt, and + * jpt+1. This option is used to efficiently evaluate the + * Jacobian numerically. + * + */ - //------------------------------------------------------ + void AxiStagnFlow::eval(int jg, doublereal* xg, + doublereal* rg, integer* diagg, doublereal rdt) { - /** - * Evaluate the residual function for axisymmetric stagnation - * flow. If jpt is less than zero, the residual function is - * evaluated at all grid points. If jpt >= 0, then the residual - * function is only evaluated at grid points jpt-1, jpt, and - * jpt+1. This option is used to efficiently evaluate the - * Jacobian numerically. - * - */ + // if evaluating a Jacobian, and the global point is outside + // the domain of influence for this domain, then skip + // evaluating the residual + if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; - void AxiStagnFlow::eval(int jg, doublereal* xg, - doublereal* rg, integer* diagg, doublereal rdt) { + // if evaluating a Jacobian, compute the steady-state residual + if (jg >= 0) rdt = 0.0; - // if evaluating a Jacobian, and the global point is outside - // the domain of influence for this domain, then skip - // evaluating the residual - if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; - - // if evaluating a Jacobian, compute the steady-state residual - if (jg >= 0) rdt = 0.0; - - // 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(); - if (jg < 0) { // evaluate all points - jmin = 0; - jmax = m_points - 1; - } - else { // evaluate points for Jacobian - jmin = max(jpt-1, 0); - jmax = min(jpt+1,m_points-1); - } + if (jg < 0) { // evaluate all points + jmin = 0; + jmax = m_points - 1; + } + else { // evaluate points for Jacobian + jmin = max(jpt-1, 0); + jmax = min(jpt+1,m_points-1); + } - // properties are computed for grid points from j0 to j1 - int j0 = max(jmin-1,0); - int j1 = min(jmax+1,m_points-1); + // properties are computed for grid points from j0 to j1 + int j0 = max(jmin-1,0); + int j1 = min(jmax+1,m_points-1); - int j, k; + int j, k; - //----------------------------------------------------- - // update properties - //----------------------------------------------------- + //----------------------------------------------------- + // update properties + //----------------------------------------------------- - // update thermodynamic properties only if a Jacobian is not - // being evaluated - if (jpt < 0) { //if (jpt < 0 || (m_transport_option == c_Multi_Transport)) { - updateThermo(x, j0, j1); + // update thermodynamic properties only if a Jacobian is not + // being evaluated + if (jpt < 0) { //if (jpt < 0 || (m_transport_option == c_Multi_Transport)) { + updateThermo(x, j0, j1); - // update transport properties only if a Jacobian is not being - // evaluated - updateTransport(x, j0, j1); - } + // update transport properties only if a Jacobian is not being + // evaluated + updateTransport(x, j0, j1); + } - // update the species diffusive mass fluxes whether or not a - // Jacobian is being evaluated - updateDiffFluxes(x, j0, j1); + // update the species diffusive mass fluxes whether or not a + // Jacobian is being evaluated + 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, dtdzj; + doublereal sum, sum2, dtdzj; - for (j = jmin; j <= jmax; j++) { + for (j = jmin; j <= jmax; j++) { - //---------------------------------------------- - // left boundary - //---------------------------------------------- + //---------------------------------------------- + // left boundary + //---------------------------------------------- - if (j == 0) { + if (j == 0) { - // these may be modified by a boundary object + // these may be modified by a boundary object - // Continuity. This propagates information right-to-left, - // since rho_u at point 0 is dependent on rho_u at point 1, - // but not on mdot from the inlet. - rsd[index(c_offset_U,0)] = - -(rho_u(x,1) - rho_u(x,0))/m_dz[0] - -(density(1)*V(x,1) + density(0)*V(x,0)); + // Continuity. This propagates information right-to-left, + // since rho_u at point 0 is dependent on rho_u at point 1, + // but not on mdot from the inlet. + rsd[index(c_offset_U,0)] = + -(rho_u(x,1) - rho_u(x,0))/m_dz[0] + -(density(1)*V(x,1) + density(0)*V(x,0)); - // the inlet (or other) object connected to this one - // will modify these equations by subtracting its values - // for V, T, and mdot. As a result, these residual equations - // will force the solution variables to the values for - // the boundary object - rsd[index(c_offset_V,0)] = V(x,0); - rsd[index(c_offset_T,0)] = T(x,0); - rsd[index(c_offset_L,0)] = -rho_u(x,0); + // the inlet (or other) object connected to this one + // will modify these equations by subtracting its values + // for V, T, and mdot. As a result, these residual equations + // will force the solution variables to the values for + // the boundary object + rsd[index(c_offset_V,0)] = V(x,0); + rsd[index(c_offset_T,0)] = T(x,0); + rsd[index(c_offset_L,0)] = -rho_u(x,0); - // The default boundary condition for species is zero - // flux. However, the boundary object may modify - // this. - sum = 0.0; - for (k = 0; k < m_nsp; k++) { - sum += Y(x,k,0); - rsd[index(c_offset_Y + k, 0)] = - -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); - } - rsd[index(c_offset_Y, 0)] = 1.0 - sum; - } + // The default boundary condition for species is zero + // flux. However, the boundary object may modify + // this. + sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,0); + rsd[index(c_offset_Y + k, 0)] = + -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); + } + rsd[index(c_offset_Y, 0)] = 1.0 - sum; + } - //---------------------------------------------- - // - // right boundary - // - //---------------------------------------------- + //---------------------------------------------- + // + // right boundary + // + //---------------------------------------------- - else if (j == m_points - 1) { + else if (j == m_points - 1) { - // the boundary object connected to the right of this - // one may modify or replace these equations. The - // default boundary conditions are zero u, V, and T, - // and zero diffusive flux for all species. + // the boundary object connected to the right of this + // one may modify or replace these equations. The + // default boundary conditions are zero u, V, and T, + // and zero diffusive flux for all species. - rsd[index(0,j)] = rho_u(x,j); - rsd[index(1,j)] = V(x,j); - rsd[index(2,j)] = T(x,j); - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - diag[index(c_offset_L, j)] = 0; - doublereal sum = 0.0; - for (k = 0; k < m_nsp; k++) { - sum += Y(x,k,j); - rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); - } - rsd[index(4,j)] = 1.0 - sum; - diag[index(4,j)] = 0; + rsd[index(0,j)] = rho_u(x,j); + rsd[index(1,j)] = V(x,j); + rsd[index(2,j)] = T(x,j); + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + doublereal sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); + } + rsd[index(4,j)] = 1.0 - sum; + diag[index(4,j)] = 0; - } + } - //------------------------------------------ - // interior points - //------------------------------------------ + //------------------------------------------ + // interior points + //------------------------------------------ - else { + else { - //---------------------------------------------- - // Continuity equation - // - // Note that this propagates the mass flow rate - // information to the left (j+1 -> j) from the - // value specified at the right boundary. The - // lambda information propagates in the opposite - // direction. - // - // d(\rho u)/dz + 2\rho V = 0 - // - //------------------------------------------------ + //---------------------------------------------- + // Continuity equation + // + // Note that this propagates the mass flow rate + // information to the left (j+1 -> j) from the + // value specified at the right boundary. The + // lambda information propagates in the opposite + // direction. + // + // d(\rho u)/dz + 2\rho V = 0 + // + //------------------------------------------------ - rsd[index(c_offset_U,j)] = - -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] - -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); + rsd[index(c_offset_U,j)] = + -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); - //algebraic constraint - diag[index(c_offset_U, j)] = 0; + //algebraic constraint + diag[index(c_offset_U, j)] = 0; - //------------------------------------------------ - // Radial momentum equation - // - // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda - // - //------------------------------------------------- - rsd[index(c_offset_V,j)] - = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) - - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] - - rdt*(V(x,j) - V_prev(j)); - diag[index(c_offset_V, j)] = 1; + //------------------------------------------------ + // Radial momentum equation + // + // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda + // + //------------------------------------------------- + rsd[index(c_offset_V,j)] + = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) + - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] + - rdt*(V(x,j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; - //------------------------------------------------- - // Species equations - // - // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k - // - //------------------------------------------------- - getWdot(x,j); + //------------------------------------------------- + // Species equations + // + // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k + // + //------------------------------------------------- + getWdot(x,j); - doublereal convec, diffus; - for (k = 0; k < m_nsp; 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; - } + doublereal convec, diffus; + for (k = 0; k < m_nsp; 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; + } - //----------------------------------------------- - // energy equation - //----------------------------------------------- + //----------------------------------------------- + // energy equation + //----------------------------------------------- - if (m_do_energy[j]) { + if (m_do_energy[j]) { - setGas(x,j); + setGas(x,j); - // heat release term - const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); - const vector_fp& cp_R = m_thermo->cp_R_ref(); + // heat release term + const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); + const vector_fp& cp_R = m_thermo->cp_R_ref(); - sum = 0.0; - sum2 = 0.0; - 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 = dTdz(x,j); - sum2 *= GasConstant * dtdzj; + sum = 0.0; + sum2 = 0.0; + 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 = 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)] = + - 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)] = - rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j)); + rsd[index(c_offset_T, j)] = + rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j)); - rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); - diag[index(c_offset_T, j)] = 1; - } + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + diag[index(c_offset_T, j)] = 1; + } - // residual equations if the energy equation is disabled + // residual equations if the energy equation is disabled - if (!m_do_energy[j]) { - rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); - diag[index(c_offset_T, 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; + } - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - diag[index(c_offset_L, j)] = 0; - } - } + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + } } + } - /** - * Update the transport properties at grid points in the range - * from j0 to j1, based on solution x. - */ - void StFlow::updateTransport(doublereal* x,int j0, int j1) { - int j,k,m; + /** + * Update the transport properties at grid points in the range + * from j0 to j1, based on solution x. + */ + void StFlow::updateTransport(doublereal* x,int j0, int j1) { + int j,k,m; - if (m_transport_option == c_Mixav_Transport) { - for (j = j0; j < j1; j++) { - setGasAtMidpoint(x,j); - m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); - m_trans->getMixDiffCoeffs(DATA_PTR(m_diff) + j*m_nsp); - m_tcon[j] = m_trans->thermalConductivity(); - } - } - else if (m_transport_option == c_Multi_Transport) { - doublereal sum, sumx, wtm, dz; - doublereal eps = 1.0e-12; - for (m = j0; m < j1; m++) { - setGasAtMidpoint(x,m); - dz = m_z[m+1] - m_z[m]; - wtm = m_thermo->meanMolecularWeight(); - - m_visc[m] = (m_dovisc ? m_trans->viscosity() : 0.0); - - m_trans->getMultiDiffCoeffs(m_nsp, - DATA_PTR(m_multidiff) + mindex(0,0,m)); - - for (k = 0; k < m_nsp; k++) { - sum = 0.0; - sumx = 0.0; - for (j = 0; j < m_nsp; j++) { - if (j != k) { - sum += m_wt[j]*m_multidiff[mindex(k,j,m)]* - ((X(x,j,m+1) - X(x,j,m))/dz + eps); - sumx += (X(x,j,m+1) - X(x,j,m))/dz; - } - } - m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps)); - } - - m_tcon[m] = m_trans->thermalConductivity(); - if (m_do_soret) { - m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp); - } - } - } + if (m_transport_option == c_Mixav_Transport) { + for (j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_visc[j] = (m_dovisc ? m_trans->viscosity() : 0.0); + m_trans->getMixDiffCoeffs(DATA_PTR(m_diff) + j*m_nsp); + m_tcon[j] = m_trans->thermalConductivity(); + } } + else if (m_transport_option == c_Multi_Transport) { + doublereal sum, sumx, wtm, dz; + doublereal eps = 1.0e-12; + for (m = j0; m < j1; m++) { + setGasAtMidpoint(x,m); + dz = m_z[m+1] - m_z[m]; + wtm = m_thermo->meanMolecularWeight(); + + m_visc[m] = (m_dovisc ? m_trans->viscosity() : 0.0); + + m_trans->getMultiDiffCoeffs(m_nsp, + DATA_PTR(m_multidiff) + mindex(0,0,m)); + + for (k = 0; k < m_nsp; k++) { + sum = 0.0; + sumx = 0.0; + for (j = 0; j < m_nsp; j++) { + if (j != k) { + sum += m_wt[j]*m_multidiff[mindex(k,j,m)]* + ((X(x,j,m+1) - X(x,j,m))/dz + eps); + sumx += (X(x,j,m+1) - X(x,j,m))/dz; + } + } + m_diff[k + m*m_nsp] = sum/(wtm*(sumx+eps)); + } + + m_tcon[m] = m_trans->thermalConductivity(); + if (m_do_soret) { + m_trans->getThermalDiffCoeffs(m_dthermal.ptrColumn(0) + m*m_nsp); + } + } + } + } - //------------------------------------------------------ + //------------------------------------------------------ - /** - * Evaluate the residual function for axisymmetric stagnation - * flow. If jpt is less than zero, the residual function is - * evaluated at all grid points. If jpt >= 0, then the residual - * function is only evaluated at grid points jpt-1, jpt, and - * jpt+1. This option is used to efficiently evaluate the - * Jacobian numerically. - * - */ + /** + * Evaluate the residual function for axisymmetric stagnation + * flow. If jpt is less than zero, the residual function is + * evaluated at all grid points. If jpt >= 0, then the residual + * function is only evaluated at grid points jpt-1, jpt, and + * jpt+1. This option is used to efficiently evaluate the + * Jacobian numerically. + * + */ - void FreeFlame::eval(int jg, doublereal* xg, - doublereal* rg, integer* diagg, doublereal rdt) { + void FreeFlame::eval(int jg, doublereal* xg, + doublereal* rg, integer* diagg, doublereal rdt) { - // if evaluating a Jacobian, and the global point is outside - // the domain of influence for this domain, then skip - // evaluating the residual - if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; + // if evaluating a Jacobian, and the global point is outside + // the domain of influence for this domain, then skip + // evaluating the residual + if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; - // if evaluating a Jacobian, compute the steady-state residual - if (jg >= 0) rdt = 0.0; + // if evaluating a Jacobian, compute the steady-state residual + if (jg >= 0) rdt = 0.0; - // 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(); - if (jg < 0) { // evaluate all points - jmin = 0; - jmax = m_points - 1; - } - else { // evaluate points for Jacobian - jmin = max(jpt-1, 0); - jmax = min(jpt+1,m_points-1); - } + if (jg < 0) { // evaluate all points + jmin = 0; + jmax = m_points - 1; + } + else { // evaluate points for Jacobian + jmin = max(jpt-1, 0); + jmax = min(jpt+1,m_points-1); + } - // properties are computed for grid points from j0 to j1 - int j0 = max(jmin-1,0); - int j1 = min(jmax+1,m_points-1); + // properties are computed for grid points from j0 to j1 + int j0 = max(jmin-1,0); + int j1 = min(jmax+1,m_points-1); - int j, k; + int j, k; - //----------------------------------------------------- - // update properties - //----------------------------------------------------- + //----------------------------------------------------- + // update properties + //----------------------------------------------------- - // update thermodynamic properties only if a Jacobian is not - // being evaluated - if (jpt < 0) { - updateThermo(x, j0, j1); - updateTransport(x, j0, j1); - } + // update thermodynamic properties only if a Jacobian is not + // being evaluated + if (jpt < 0) { + updateThermo(x, j0, j1); + updateTransport(x, j0, j1); + } - // update the species diffusive mass fluxes whether or not a - // Jacobian is being evaluated - updateDiffFluxes(x, j0, j1); + // update the species diffusive mass fluxes whether or not a + // Jacobian is being evaluated + 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, dtdzj; + doublereal sum, sum2, dtdzj; - for (j = jmin; j <= jmax; j++) { + for (j = jmin; j <= jmax; j++) { - //---------------------------------------------- - // left boundary - //---------------------------------------------- + //---------------------------------------------- + // left boundary + //---------------------------------------------- - if (j == 0) { + if (j == 0) { - // these may be modified by a boundary object + // these may be modified by a boundary object - // Continuity. This propagates information right-to-left, - // since rho_u at point 0 is dependent on rho_u at point 1, - // but not on mdot from the inlet. - rsd[index(c_offset_U,0)] = - -(rho_u(x,1) - rho_u(x,0))/m_dz[0] - -(density(1)*V(x,1) + density(0)*V(x,0)); + // Continuity. This propagates information right-to-left, + // since rho_u at point 0 is dependent on rho_u at point 1, + // but not on mdot from the inlet. + rsd[index(c_offset_U,0)] = + -(rho_u(x,1) - rho_u(x,0))/m_dz[0] + -(density(1)*V(x,1) + density(0)*V(x,0)); - // the inlet (or other) object connected to this one - // will modify these equations by subtracting its values - // for V, T, and mdot. As a result, these residual equations - // will force the solution variables to the values for - // the boundary object - rsd[index(c_offset_V,0)] = V(x,0); - rsd[index(c_offset_T,0)] = T(x,0); - rsd[index(c_offset_L,0)] = -rho_u(x,0); + // the inlet (or other) object connected to this one + // will modify these equations by subtracting its values + // for V, T, and mdot. As a result, these residual equations + // will force the solution variables to the values for + // the boundary object + rsd[index(c_offset_V,0)] = V(x,0); + rsd[index(c_offset_T,0)] = T(x,0); + rsd[index(c_offset_L,0)] = -rho_u(x,0); - // The default boundary condition for species is zero - // flux - sum = 0.0; - for (k = 0; k < m_nsp; k++) { - sum += Y(x,k,0); - rsd[index(c_offset_Y + k, 0)] = - -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); - } - rsd[index(c_offset_Y, 0)] = 1.0 - sum; - } + // The default boundary condition for species is zero + // flux + sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,0); + rsd[index(c_offset_Y + k, 0)] = + -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); + } + rsd[index(c_offset_Y, 0)] = 1.0 - sum; + } - //---------------------------------------------- - // - // right boundary - // - //---------------------------------------------- + //---------------------------------------------- + // + // right boundary + // + //---------------------------------------------- - else if (j == m_points - 1) { + else if (j == m_points - 1) { - // the boundary object connected to the right of this - // one may modify or replace these equations. The - // default boundary conditions are zero u, V, and T, - // and zero diffusive flux for all species. + // the boundary object connected to the right of this + // one may modify or replace these equations. The + // default boundary conditions are zero u, V, and T, + // and zero diffusive flux for all species. - // zero gradient - rsd[index(0,j)] = rho_u(x,j) - rho_u(x,j-1); - rsd[index(1,j)] = V(x,j); - rsd[index(2,j)] = T(x,j) - T(x,j-1); - doublereal sum = 0.0; - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - diag[index(c_offset_L, j)] = 0; - for (k = 0; k < m_nsp; k++) { - sum += Y(x,k,j); - rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); - } - rsd[index(4,j)] = 1.0 - sum; - diag[index(4,j)] = 0; - } + // zero gradient + rsd[index(0,j)] = rho_u(x,j) - rho_u(x,j-1); + rsd[index(1,j)] = V(x,j); + rsd[index(2,j)] = T(x,j) - T(x,j-1); + doublereal sum = 0.0; + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); + } + rsd[index(4,j)] = 1.0 - sum; + diag[index(4,j)] = 0; + } - //------------------------------------------ - // interior points - //------------------------------------------ + //------------------------------------------ + // interior points + //------------------------------------------ - else { + else { - //---------------------------------------------- - // Continuity equation - //---------------------------------------------- + //---------------------------------------------- + // Continuity equation + //---------------------------------------------- - if (grid(j) > m_zfixed){ - rsd[index(c_offset_U,j)] = - - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1] - - (density(j-1)*V(x,j-1) + density(j)*V(x,j)); - } + if (grid(j) > m_zfixed){ + rsd[index(c_offset_U,j)] = + - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1] + - (density(j-1)*V(x,j-1) + density(j)*V(x,j)); + } - else if (grid(j) == m_zfixed){ - if (m_do_energy[j]) { - rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); - } - else { - rsd[index(c_offset_U,j)] = (rho_u(x,j) - - m_rho[0]*0.3); - } - } - else if(grid(j) < m_zfixed){ - rsd[index(c_offset_U,j)] = - - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j] - - (density(j+1)*V(x,j+1) + density(j)*V(x,j)); - } - //algebraic constraint - diag[index(c_offset_U, j)] = 0; + else if (grid(j) == m_zfixed){ + if (m_do_energy[j]) { + rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + } + else { + rsd[index(c_offset_U,j)] = (rho_u(x,j) + - m_rho[0]*0.3); + } + } + else if(grid(j) < m_zfixed){ + rsd[index(c_offset_U,j)] = + - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + - (density(j+1)*V(x,j+1) + density(j)*V(x,j)); + } + //algebraic constraint + diag[index(c_offset_U, j)] = 0; - //------------------------------------------------ - // Radial momentum equation - // - // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda - // - //------------------------------------------------- - rsd[index(c_offset_V,j)] - = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) - - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] - - rdt*(V(x,j) - V_prev(j)); - diag[index(c_offset_V, j)] = 1; + //------------------------------------------------ + // Radial momentum equation + // + // \rho u dV/dz + \rho V^2 = d(\mu dV/dz)/dz - lambda + // + //------------------------------------------------- + rsd[index(c_offset_V,j)] + = (shear(x,j) - lambda(x,j) - rho_u(x,j)*dVdz(x,j) + - m_rho[j]*V(x,j)*V(x,j))/m_rho[j] + - rdt*(V(x,j) - V_prev(j)); + diag[index(c_offset_V, j)] = 1; - //------------------------------------------------- - // Species equations - // - // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k - // - //------------------------------------------------- - getWdot(x,j); + //------------------------------------------------- + // Species equations + // + // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k + // + //------------------------------------------------- + getWdot(x,j); - doublereal convec, diffus; - for (k = 0; k < m_nsp; 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; - } + doublereal convec, diffus; + for (k = 0; k < m_nsp; 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; + } - //----------------------------------------------- - // energy equation - //----------------------------------------------- + //----------------------------------------------- + // energy equation + //----------------------------------------------- - if (m_do_energy[j]) { + if (m_do_energy[j]) { - setGas(x,j); + setGas(x,j); - // heat release term - const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); - const vector_fp& cp_R = m_thermo->cp_R_ref(); + // heat release term + const vector_fp& h_RT = m_thermo->enthalpy_RT_ref(); + const vector_fp& cp_R = m_thermo->cp_R_ref(); - sum = 0.0; - sum2 = 0.0; - 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 = dTdz(x,j); - sum2 *= GasConstant * dtdzj; + sum = 0.0; + sum2 = 0.0; + 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 = 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)] = + - 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)] = - rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j)); + rsd[index(c_offset_T, j)] = + rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j)); - rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); - diag[index(c_offset_T, j)] = 1; - } - // residual equations if the energy equation is disabled - else { - rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); - diag[index(c_offset_T, j)] = 0; - } + rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j)); + diag[index(c_offset_T, j)] = 1; + } + // residual equations if the energy equation is disabled + else { + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); + diag[index(c_offset_T, j)] = 0; + } - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - diag[index(c_offset_L, j)] = 0; - } - } + rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); + diag[index(c_offset_L, j)] = 0; + } } + } - /** - * Print the solution. - */ - void StFlow::showSolution(const doublereal* x) { - int nn = m_nv/5; - int i, j, n; - //char* buf = new char[100]; - char buf[100]; + /** + * Print the solution. + */ + void StFlow::showSolution(const doublereal* x) { + int nn = m_nv/5; + int i, j, n; + //char* buf = new char[100]; + char buf[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); - sprintf(buf, " Pressure: %10.4g Pa \n", m_press); - writelog(buf); - for (i = 0; i < nn; i++) { - st_drawline(); - sprintf(buf, "\n z "); - writelog(buf); - for (n = 0; n < 5; n++) { - sprintf(buf, " %10s ",componentName(i*5 + n).c_str()); - writelog(buf); - } - st_drawline(); - for (j = 0; j < m_points; j++) { - sprintf(buf, "\n %10.4g ",m_z[j]); - writelog(buf); - for (n = 0; n < 5; n++) { - sprintf(buf, " %10.4g ",component(x, i*5+n,j)); - writelog(buf); - } - } - writelog("\n"); - } - int nrem = m_nv - 5*nn; - st_drawline(); - sprintf(buf, "\n z "); - writelog(buf); - for (n = 0; n < nrem; n++) { - sprintf(buf, " %10s ", componentName(nn*5 + n).c_str()); - writelog(buf); - } - st_drawline(); - for (j = 0; j < m_points; j++) { - sprintf(buf, "\n %10.4g ",m_z[j]); - writelog(buf); - for (n = 0; n < nrem; n++) { - sprintf(buf, " %10.4g ",component(x, nn*5+n,j)); - writelog(buf); - } - } - writelog("\n"); + sprintf(buf, " Pressure: %10.4g Pa \n", m_press); + writelog(buf); + for (i = 0; i < nn; i++) { + st_drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10s ",componentName(i*5 + n).c_str()); + writelog(buf); + } + st_drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10.4g ",component(x, i*5+n,j)); + writelog(buf); + } + } + writelog("\n"); } + int nrem = m_nv - 5*nn; + st_drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10s ", componentName(nn*5 + n).c_str()); + writelog(buf); + } + st_drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10.4g ",component(x, nn*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } - /** - * Update the diffusive mass fluxes. - */ - void StFlow::updateDiffFluxes(const doublereal* x, int j0, int j1) { - int j, k, m; - doublereal sum, wtm, rho, dz, gradlogT; + /** + * Update the diffusive mass fluxes. + */ + void StFlow::updateDiffFluxes(const doublereal* x, int j0, int j1) { + int j, k, m; + doublereal sum, wtm, rho, dz, gradlogT; - switch (m_transport_option) { + switch (m_transport_option) { - case c_Mixav_Transport: - case c_Multi_Transport: - for (j = j0; j < j1; j++) { - sum = 0.0; - wtm = m_wtm[j]; - rho = density(j); - dz = z(j+1) - z(j); + case c_Mixav_Transport: + case c_Multi_Transport: + for (j = j0; j < j1; j++) { + sum = 0.0; + wtm = m_wtm[j]; + rho = density(j); + dz = z(j+1) - z(j); - for (k = 0; k < m_nsp; k++) { - m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); - m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; - sum -= m_flux(k,j); - } - // correction flux to insure that \sum_k Y_k V_k = 0. - for (k = 0; k < m_nsp; k++) m_flux(k,j) += sum*Y(x,k,j); - } - break; + for (k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); + m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; + sum -= m_flux(k,j); + } + // correction flux to insure that \sum_k Y_k V_k = 0. + for (k = 0; k < m_nsp; k++) m_flux(k,j) += sum*Y(x,k,j); + } + break; - default: - throw CanteraError("updateDiffFluxes","unknown transport model"); - } - - if (m_do_soret) { - for (m = j0; m < j1; m++) { - gradlogT = 2.0*(T(x,m+1) - T(x,m))/(T(x,m+1) + T(x,m)); - for (k = 0; k < m_nsp; k++) { - m_flux(k,m) -= m_dthermal(k,m)*gradlogT; - } - } - } + default: + throw CanteraError("updateDiffFluxes","unknown transport model"); } - - string StFlow::componentName(int n) const { - switch(n) { - case 0: return "u"; - case 1: return "V"; - case 2: return "T"; - case 3: return "lambda"; - default: - if (n >= (int) c_offset_Y && n < (int) (c_offset_Y + m_nsp)) { - return m_thermo->speciesName(n - c_offset_Y); - } - else - return ""; - } + if (m_do_soret) { + for (m = j0; m < j1; m++) { + gradlogT = 2.0*(T(x,m+1) - T(x,m))/(T(x,m+1) + T(x,m)); + for (k = 0; k < m_nsp; k++) { + m_flux(k,m) -= m_dthermal(k,m)*gradlogT; + } + } } + } - //added by Karl Meredith - int StFlow::componentIndex(string name) const { + string StFlow::componentName(int n) const { + switch(n) { + case 0: return "u"; + case 1: return "V"; + case 2: return "T"; + case 3: return "lambda"; + default: + if (n >= (int) c_offset_Y && n < (int) (c_offset_Y + m_nsp)) { + return m_thermo->speciesName(n - c_offset_Y); + } + else + return ""; + } + } + + + //added by Karl Meredith + int StFlow::componentIndex(string name) const { - if(name=="u") {return 0;} - else if (name=="V") {return 1;} - else if (name=="T") {return 2;} - else if (name=="lambda") {return 3;} - else { - for (int n=4;n ignored; - int nsp = m_thermo->nSpecies(); - vector_int did_species(nsp, 0); - - vector str; - dom.getChildren("string",str); - int nstr = static_cast(str.size()); - for (int istr = 0; istr < nstr; istr++) { - const XML_Node& nd = *str[istr]; - writelog(nd["title"]+": "+nd.value()+"\n"); - } - - //map params; - double pp = -1.0; - pp = getFloat(dom, "pressure", "pressure"); - setPressure(pp); + return -1; + } - vector d; - dom.child("grid_data").getChildren("floatArray",d); - int nd = static_cast(d.size()); + void StFlow::restore(const XML_Node& dom, doublereal* soln) { - vector_fp x; - int n, np = 0, j, ks, k; - 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(); - writelog("Grid contains "+int2str(np)+ - " points.\n"); - readgrid = true; - setupGrid(np, DATA_PTR(x)); - } - } - if (!readgrid) { - throw CanteraError("StFlow::restore", - "domain contains no grid points."); - } + vector ignored; + int nsp = m_thermo->nSpecies(); + vector_int did_species(nsp, 0); - 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") { - ; // already read grid - } - 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]; - - // 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()); - setFixedTempProfile(zz, x); - } - 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 = static_cast(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)+" "); - } - } - - return; - error: - throw CanteraError("StFlow::restore","Data size error"); + vector str; + dom.getChildren("string",str); + int nstr = static_cast(str.size()); + for (int istr = 0; istr < nstr; istr++) { + const XML_Node& nd = *str[istr]; + writelog(nd["title"]+": "+nd.value()+"\n"); } + //map params; + double pp = -1.0; + pp = getFloat(dom, "pressure", "pressure"); + setPressure(pp); - void StFlow::save(XML_Node& o, const doublereal * const sol) { - int k; + vector d; + dom.child("grid_data").getChildren("floatArray",d); + int nd = static_cast(d.size()); - ArrayViewer soln(m_nv, m_points, const_cast(sol) + loc()); - - XML_Node& flow = (XML_Node&)o.addChild("domain"); - flow.addAttribute("type",flowType()); - flow.addAttribute("id",m_id); - flow.addAttribute("points",m_points); - flow.addAttribute("components",m_nv); - - if (m_desc != "") addString(flow,"description",m_desc); - XML_Node& gv = flow.addChild("grid_data"); - addFloat(flow, "pressure", m_press, "Pa", "pressure"); - addFloatArray(gv,"z",m_z.size(),DATA_PTR(m_z), - "m","length"); - vector_fp x(static_cast(soln.nColumns())); - - soln.getRow(0,DATA_PTR(x)); - addFloatArray(gv,"u",x.size(),DATA_PTR(x),"m/s","velocity"); - - soln.getRow(1,DATA_PTR(x)); - addFloatArray(gv,"V", - x.size(),DATA_PTR(x),"1/s","rate"); - - soln.getRow(2,DATA_PTR(x)); - addFloatArray(gv,"T",x.size(),DATA_PTR(x),"K","temperature",0.0); - - soln.getRow(3,DATA_PTR(x)); - addFloatArray(gv,"L",x.size(),DATA_PTR(x),"N/m^4"); - - for (k = 0; k < m_nsp; k++) { - soln.getRow(4+k,DATA_PTR(x)); - addFloatArray(gv,m_thermo->speciesName(k), - x.size(),DATA_PTR(x),"","massFraction",0.0,1.0); - } + vector_fp x; + int n, np = 0, j, ks, k; + 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(); + writelog("Grid contains "+int2str(np)+ + " points.\n"); + readgrid = true; + setupGrid(np, DATA_PTR(x)); + } + } + if (!readgrid) { + throw CanteraError("StFlow::restore", + "domain 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") { + ; // already read grid + } + 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]; - void StFlow::setJac(MultiJac* jac) { - m_jac = jac; + // 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()); + setFixedTempProfile(zz, x); + } + 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 = static_cast(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)+" "); + } + } + + return; + error: + throw CanteraError("StFlow::restore","Data size error"); + } + + + + void StFlow::save(XML_Node& o, const doublereal * const sol) { + int k; + + ArrayViewer soln(m_nv, m_points, const_cast(sol) + loc()); + + XML_Node& flow = (XML_Node&)o.addChild("domain"); + flow.addAttribute("type",flowType()); + flow.addAttribute("id",m_id); + flow.addAttribute("points",m_points); + flow.addAttribute("components",m_nv); + + if (m_desc != "") addString(flow,"description",m_desc); + XML_Node& gv = flow.addChild("grid_data"); + addFloat(flow, "pressure", m_press, "Pa", "pressure"); + addFloatArray(gv,"z",m_z.size(),DATA_PTR(m_z), + "m","length"); + vector_fp x(static_cast(soln.nColumns())); + + soln.getRow(0,DATA_PTR(x)); + addFloatArray(gv,"u",x.size(),DATA_PTR(x),"m/s","velocity"); + + soln.getRow(1,DATA_PTR(x)); + addFloatArray(gv,"V", + x.size(),DATA_PTR(x),"1/s","rate"); + + soln.getRow(2,DATA_PTR(x)); + addFloatArray(gv,"T",x.size(),DATA_PTR(x),"K","temperature",0.0); + + soln.getRow(3,DATA_PTR(x)); + addFloatArray(gv,"L",x.size(),DATA_PTR(x),"N/m^4"); + + for (k = 0; k < m_nsp; k++) { + soln.getRow(4+k,DATA_PTR(x)); + addFloatArray(gv,m_thermo->speciesName(k), + x.size(),DATA_PTR(x),"","massFraction",0.0,1.0); + } + } + + + void StFlow::setJac(MultiJac* jac) { + m_jac = jac; + } + } // namespace diff --git a/Cantera/src/thermo/PDSS_SSVol.cpp b/Cantera/src/thermo/PDSS_SSVol.cpp index 78a4260dd..21ec884a6 100644 --- a/Cantera/src/thermo/PDSS_SSVol.cpp +++ b/Cantera/src/thermo/PDSS_SSVol.cpp @@ -133,12 +133,20 @@ namespace Cantera { m_constMolarVolume = getFloat(*ss, "molarVolume", "toSI"); } else if (model == "temperature_polynomial") { volumeModel_ = cSSVOLUME_TPOLY; - getFloatArray(*ss, TCoeff_, true, "", "floatArray"); + int num = getFloatArray(*ss, TCoeff_, true, "", "volumeTemperaturePolynomial"); + if (num != 4) { + throw CanteraError("PDSS_SSVol::constructPDSSXML", + " Didn't get 3 density polynomial numbers for species " + speciesNode.name()); + } } else if (model == "density_temperature_polynomial") { volumeModel_ = cSSVOLUME_DENSITY_TPOLY; - getFloatArray(*ss, TCoeff_, true, "", "floatArray"); + int num = getFloatArray(*ss, TCoeff_, true, "", "densityTemperaturePolynomial"); + if (num != 4) { + throw CanteraError("PDSS_SSVol::constructPDSSXML", + " Didn't get 3 density polynomial numbers for species " + speciesNode.name()); + } } else { - throw CanteraError("PDSS_SSVol::initThermoXML", + throw CanteraError("PDSS_SSVol::constructPDSSXML", "standardState model for species isn't constant_incompressible: " + speciesNode.name()); } std::string id = ""; @@ -327,15 +335,15 @@ namespace Cantera { if (volumeModel_ == cSSVOLUME_CONSTANT ) { m_Vss_ptr[m_spindex] = m_constMolarVolume; } else if (volumeModel_ == cSSVOLUME_TPOLY) { - m_Vss_ptr[m_spindex] = TCoeff_[0] + m_temp * (TCoeff_[1] + m_temp * TCoeff_[2]); - dVdT_ = TCoeff_[1] + 2.0 * m_temp * TCoeff_[2]; - d2VdT2_ = 2.0 * TCoeff_[2]; + m_Vss_ptr[m_spindex] = TCoeff_[0] + m_temp * (TCoeff_[1] + m_temp * (TCoeff_[2] + m_temp * TCoeff_[3])); + dVdT_ = TCoeff_[1] + 2.0 * m_temp * TCoeff_[2] + 3.0 * m_temp * m_temp * TCoeff_[3]; + d2VdT2_ = 2.0 * TCoeff_[2] + 6.0 * m_temp * TCoeff_[3]; } else if (volumeModel_ == cSSVOLUME_DENSITY_TPOLY) { - doublereal dens = (TCoeff_[0] + m_temp * (TCoeff_[1] + m_temp * TCoeff_[2])); + doublereal dens = TCoeff_[0] + m_temp * (TCoeff_[1] + m_temp * (TCoeff_[2] + m_temp * TCoeff_[3])); m_Vss_ptr[m_spindex] = m_mw / dens; doublereal dens2 = dens * dens; - doublereal ddensdT = TCoeff_[1] + 2.0 * m_temp * TCoeff_[2]; - doublereal d2densdT2 = 2.0 * TCoeff_[2]; + doublereal ddensdT = TCoeff_[1] + 2.0 * m_temp * TCoeff_[2] + 3.0 * m_temp * m_temp * TCoeff_[3]; + doublereal d2densdT2 = 2.0 * TCoeff_[2] + 6.0 * m_temp * TCoeff_[3]; dVdT_ = - m_mw / (dens2) * (ddensdT); d2VdT2_ = 2.0 * m_mw / (dens2 * dens) * ddensdT * ddensdT - m_mw / dens2 * d2densdT2; } else { diff --git a/Cantera/src/thermo/VPSSMgrFactory.cpp b/Cantera/src/thermo/VPSSMgrFactory.cpp index 0993635bb..186bf5576 100644 --- a/Cantera/src/thermo/VPSSMgrFactory.cpp +++ b/Cantera/src/thermo/VPSSMgrFactory.cpp @@ -114,6 +114,10 @@ namespace Cantera { } else if (ssModel == "constant_incompressible" || ssModel == "constantVolume") { has_nasa_constVol++; + } else if (ssModel == "temperature_polynomial" || + ssModel == "density_temperature_polynomial" || + ssModel == "constant") { + has_other++; } else { throw UnknownVPSSMgrModel("getVPSSMgrTypes:", spNode->attrib("name")); @@ -126,6 +130,10 @@ namespace Cantera { } else if (ssModel == "constant_incompressible" || ssModel == "constantVolume") { has_shomate_constVol++; + } else if (ssModel == "temperature_polynomial" || + ssModel == "density_temperature_polynomial" || + ssModel == "constant") { + has_other++; } else { throw UnknownVPSSMgrModel("getVPSSMgrTypes:", spNode->attrib("name")); @@ -138,6 +146,10 @@ namespace Cantera { } else if (ssModel == "constant_incompressible" || ssModel == "constantVolume") { has_simple_constVol++; + } else if (ssModel == "temperature_polynomial" || + ssModel == "density_temperature_polynomial" || + ssModel == "constant") { + has_other++; } else { throw UnknownVPSSMgrModel("getVPSSMgrTypes:", spNode->attrib("name"));