From 704e0843470cec9f3550e4a0e79731616cfe9386 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 6 Nov 2006 15:20:34 +0000 Subject: [PATCH] added support for modifying NASA polynomial coefficients --- Cantera/src/Constituents.h | 21 ---------- Cantera/src/NasaPoly1.h | 33 ++++++++++++---- Cantera/src/NasaThermo.h | 57 +++++++++++++++++++++++---- Cantera/src/SpeciesThermo.h | 3 ++ Cantera/src/SpeciesThermoFactory.cpp | 10 ++--- Cantera/src/SpeciesThermoInterpType.h | 3 ++ 6 files changed, 85 insertions(+), 42 deletions(-) diff --git a/Cantera/src/Constituents.h b/Cantera/src/Constituents.h index a46e57457..3832c647b 100755 --- a/Cantera/src/Constituents.h +++ b/Cantera/src/Constituents.h @@ -194,27 +194,6 @@ namespace Cantera { */ doublereal size(int k) const { return m_speciesSize[k]; } -#ifdef INCL_DEPRECATED_METHODS - /** - * Return a SpeciesData structure containing species data. - */ - SpeciesData species(int k) const { - if (k < 0 || k >= nSpecies()) - throw SpeciesRangeError("Constituents::charge",k,nSpecies()); - SpeciesData s; - s.name = m_speciesNames[k]; - s.phase = m_speciesPhase[k]; - int offset = m_mm * k; - s.atoms.resize(m_mm); - for (int m = 0; m < m_mm; m++) { - s.atoms[m] = m_speciesComp[offset + m]; - } - s.charge = m_speciesCharge[k]; - s.molecularWeight = m_weight[k]; - return s; - } -#endif - /** * Prohibit addition of more species, and prepare for * calculations with this set of elements and species. diff --git a/Cantera/src/NasaPoly1.h b/Cantera/src/NasaPoly1.h index d639d7269..6613e2e73 100755 --- a/Cantera/src/NasaPoly1.h +++ b/Cantera/src/NasaPoly1.h @@ -94,9 +94,9 @@ namespace Cantera { return (SpeciesThermoInterpType *) np; } - doublereal minTemp() const { return m_lowT;} - doublereal maxTemp() const { return m_highT;} - doublereal refPressure() const { return m_Pref; } + virtual doublereal minTemp() const { return m_lowT;} + virtual doublereal maxTemp() const { return m_highT;} + virtual doublereal refPressure() const { return m_Pref; } virtual int reportType() const { return NASA1; } /** @@ -116,7 +116,7 @@ namespace Cantera { * tt[4] = 1.0/t; * tt[5] = log(t); */ - void updateProperties(const doublereal* tt, + virtual void updateProperties(const doublereal* tt, doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const { doublereal ct0 = m_coeff[2]; // a0 @@ -148,7 +148,7 @@ namespace Cantera { * * (note: this is slow, but it is general) */ - void updatePropertiesTemp(const doublereal temp, + virtual void updatePropertiesTemp(const doublereal temp, doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const { double tPoly[6]; @@ -161,7 +161,7 @@ namespace Cantera { updateProperties(tPoly, cp_R, h_RT, s_R); } - void reportParameters(int &n, int &type, + virtual void reportParameters(int &n, int &type, doublereal &tlow, doublereal &thigh, doublereal &pref, doublereal* const coeffs) const { @@ -170,11 +170,28 @@ namespace Cantera { tlow = m_lowT; thigh = m_highT; pref = m_Pref; - for (int i = 0; i < 7; i++) { - coeffs[i] = m_coeff[i]; + coeffs[5] = m_coeff[0]; + coeffs[6] = m_coeff[1]; + for (int i = 2; i < 7; i++) { + coeffs[i-2] = m_coeff[i]; } +#ifdef WARN_ABOUT_CHANGES_FROM_VERSION_1_6 + cout << "************************************************\n" + cout << "Warning: NasaPoly1::reportParameters now returns \n" + << "the coefficient array in the same order as in\n" + << "the input file. See file NasaPoly1.h" << endl; + cout << "************************************************\n" +#endif } + virtual void modifyParameters(doublereal* coeffs) { + m_coeff[0] = coeffs[5]; + m_coeff[1] = coeffs[6]; + for (int i = 0; i < 5; i++) { + m_coeff[i+2] = coeffs[i]; + } + } + protected: doublereal m_lowT; // lowest valid temperature diff --git a/Cantera/src/NasaThermo.h b/Cantera/src/NasaThermo.h index 87ebdccca..41a79a357 100755 --- a/Cantera/src/NasaThermo.h +++ b/Cantera/src/NasaThermo.h @@ -74,6 +74,7 @@ namespace Cantera { doublereal minTemp, doublereal maxTemp, doublereal refPressure) { + m_name[index] = name; int imid = int(c[0]); // midpoint temp converted to integer int igrp = m_index[imid]; // has this value been seen before? if (igrp == 0) { // if not, prepare new group @@ -97,13 +98,21 @@ namespace Cantera { vector_fp chigh(7); copy(c + 8, c + 15, chigh.begin()); - checkContinuity(name, tmid, clow, &chigh[0]); - m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, pref, &chigh[0])); m_low[igrp-1].push_back(NasaPoly1(index, tlow, tmid, pref, clow)); + vector_fp clu(7), chu(7); + clu[5] = clow[0]; + clu[6] = clow[1]; + copy(clow+2, clow+7, clu.begin()); + chu[5] = chigh[0]; + chu[6] = chigh[1]; + copy(chigh.begin()+2, chigh.begin()+7, chu.begin()); + + checkContinuity(name, tmid, &clu[0], &chu[0]); + if (tlow > m_tlow_max) m_tlow_max = tlow; if (thigh < m_thigh_min) m_thigh_min = thigh; m_tlow.push_back(tlow); @@ -244,6 +253,37 @@ namespace Cantera { } } + + /** + * This utility function modifies the array of coefficients. + * The array is the same as that returned by reportParams, so + * a call can first be made to reportParams to populate the + * array, and then modifyParams can be called to alter + * selected values. For the NASA object, there are 15 + * coefficients. + */ + virtual void modifyParams(int index, doublereal *c) { + int type = reportType(index); + if (type == NASA) { + int grp = m_group_map[index]; + int pos = m_posInGroup_map[index]; + vector &mlg = m_low[grp-1]; + vector &mhg = m_high[grp-1]; + NasaPoly1 *lowPoly = &(mlg[pos]); + NasaPoly1 *highPoly = &(mhg[pos]); + doublereal tmid = lowPoly->maxTemp(); + if (c[0] != tmid) { + throw CanteraError(" ", "Tmid cannot be changed"); + } + lowPoly->modifyParameters(c + 1); + highPoly->modifyParameters(c + 8); + checkContinuity(m_name[index], c[0], c + 1, c + 8); + } else { + throw CanteraError(" ", "confused"); + } + } + + protected: vector > m_high; @@ -271,6 +311,7 @@ namespace Cantera { * temperature polynomials for that species are storred. */ mutable map m_posInGroup_map; + mutable map m_name; private: @@ -280,16 +321,16 @@ namespace Cantera { /// for internal use by checkContinuity doublereal enthalpy_RT(double t, const doublereal* c) { - return c[2] + 0.5*c[3]*t + OneThird*c[4]*t*t - + 0.25*c[5]*t*t*t + 0.2*c[6]*t*t*t*t - + c[0]/t; + return c[0] + 0.5*c[1]*t + OneThird*c[2]*t*t + + 0.25*c[3]*t*t*t + 0.2*c[4]*t*t*t*t + + c[5]/t; } /// for internal use by checkContinuity doublereal entropy_R(double t, const doublereal* c) { - return c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t - + OneThird*c[5]*t*t*t + 0.25*c[6]*t*t*t*t - + c[1]; + return c[0]*log(t) + c[1]*t + 0.5*c[2]*t*t + + OneThird*c[3]*t*t*t + 0.25*c[4]*t*t*t*t + + c[6]; } }; diff --git a/Cantera/src/SpeciesThermo.h b/Cantera/src/SpeciesThermo.h index c19193cd1..b2efc6b8f 100755 --- a/Cantera/src/SpeciesThermo.h +++ b/Cantera/src/SpeciesThermo.h @@ -159,6 +159,9 @@ namespace Cantera { doublereal &minTemp, doublereal &maxTemp, doublereal &refPressure)=0; + + virtual void modifyParams(int index, doublereal *c) {} + }; } diff --git a/Cantera/src/SpeciesThermoFactory.cpp b/Cantera/src/SpeciesThermoFactory.cpp index 95bf505c0..63ef63b8d 100755 --- a/Cantera/src/SpeciesThermoFactory.cpp +++ b/Cantera/src/SpeciesThermoFactory.cpp @@ -160,11 +160,11 @@ namespace Cantera { doublereal* chigh) { // heat capacity - doublereal cplow = poly4(tmid, clow+2); - doublereal cphigh = poly4(tmid, chigh+2); + doublereal cplow = poly4(tmid, clow); + doublereal cphigh = poly4(tmid, chigh); doublereal delta = cplow - cphigh; if (fabs(delta/cplow) > 0.001) { - writelog("\n**** WARNING ****\nFor species "+name+ + writelog("\n\n**** WARNING ****\nFor species "+name+ ", discontinuity in cp/R detected at Tmid = " +fp2str(tmid)+"\n"); writelog("\tValue computed using low-temperature polynomial: " @@ -178,7 +178,7 @@ namespace Cantera { doublereal hrthigh = enthalpy_RT(tmid, chigh); delta = hrtlow - hrthigh; if (fabs(delta/hrtlow) > 0.001) { - writelog("\n**** WARNING ****\nFor species "+name+ + writelog("\n\n**** WARNING ****\nFor species "+name+ ", discontinuity in h/RT detected at Tmid = " +fp2str(tmid)+"\n"); writelog("\tValue computed using low-temperature polynomial: " @@ -192,7 +192,7 @@ namespace Cantera { doublereal srhigh = entropy_R(tmid, chigh); delta = srlow - srhigh; if (fabs(delta/srlow) > 0.001) { - writelog("\n**** WARNING ****\nFor species "+name+ + writelog("\n\n**** WARNING ****\nFor species "+name+ ", discontinuity in s/R detected at Tmid = " +fp2str(tmid)+"\n"); writelog("\tValue computed using low-temperature polynomial: " diff --git a/Cantera/src/SpeciesThermoInterpType.h b/Cantera/src/SpeciesThermoInterpType.h index 9a443e92f..07e62ce3a 100644 --- a/Cantera/src/SpeciesThermoInterpType.h +++ b/Cantera/src/SpeciesThermoInterpType.h @@ -47,6 +47,9 @@ namespace Cantera { doublereal &tlow, doublereal &thigh, doublereal &pref, doublereal* const coeffs) const = 0; + + virtual void modifyParameters(doublereal* coeffs) {} + }; }