自学教程：C++ thermo函数代码示例

/** * Update the equilibrium constants in molar units. */void AqueousKinetics::updateKc(){    doublereal rt = GasConstant * m_temp;    thermo().getStandardChemPotentials(&m_grt[0]);    fill(m_rkcn.begin(), m_rkcn.end(), 0.0);    for (size_t k = 0; k < thermo().nSpecies(); k++) {        doublereal logStandConc_k = thermo().logStandardConc(k);        m_grt[k] -= rt * logStandConc_k;    }    // compute Delta G^0 for all reversible reactions    m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);    //doublereal logStandConc = m_kdata->m_logStandConc;    doublereal rrt = 1.0/(GasConstant * thermo().temperature());    for (size_t i = 0; i < m_nrev; i++) {        size_t irxn = m_revindex[i];        m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);    }    for (size_t i = 0; i != m_nirrev; ++i) {        m_rkcn[ m_irrev[i] ] = 0.0;    }}

void GasKinetics::_update_rates_C(){    thermo().getActivityConcentrations(&m_conc[0]);    doublereal ctot = thermo().molarDensity();    // 3-body reactions    if (!concm_3b_values.empty()) {        m_3b_concm.update(m_conc, ctot, &concm_3b_values[0]);    }    // Falloff reactions    if (!concm_falloff_values.empty()) {        m_falloff_concm.update(m_conc, ctot, &concm_falloff_values[0]);    }    // P-log reactions    if (m_plog_rates.nReactions()) {        double logP = log(thermo().pressure());        m_plog_rates.update_C(&logP);    }    // Chebyshev reactions    if (m_cheb_rates.nReactions()) {        double log10P = log10(thermo().pressure());        m_cheb_rates.update_C(&log10P);    }    m_ROP_ok = false;}

    /**     * Get the equilibrium constants of all reactions, whether     * reversible or not.     */    void InterfaceKinetics::getEquilibriumConstants(doublereal* kc) {        int i;        int n, nsp, k, ik=0;        doublereal rt = GasConstant*thermo(0).temperature();        doublereal rrt = 1.0/rt;        int np = nPhases();        for (n = 0; n < np; n++) {            thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);            nsp = thermo(n).nSpecies();            for (k = 0; k < nsp; k++) {                m_mu0[ik] -= rt*thermo(n).logStandardConc(k);                m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);                ik++;            }        }        fill(kc, kc + m_ii, 0.0);        m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), kc);        for (i = 0; i < m_ii; i++) {            kc[i] = exp(-kc[i]*rrt);        }    }

  void InterfaceKinetics::getExchangeCurrentQuantities() {    /*     * Combine vectors of the standard Gibbs free energies of the     * species and the standard concentrations to get change in      * standard chemical potential for reaction.     * Outputs:     *   - m_deltaG0  -- stores standard state chemical potentials     *   - m_ProdStanConcReac -- stores products of standard concentrations     */    int ik = 0;      int np = nPhases();    for (int n = 0; n < np; n++) {      thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);      int nsp = thermo(n).nSpecies();      for (int k = 0; k < nsp; k++) {        m_StandardConc[ik] = thermo(n).standardConcentration(k);        ik++;      }    }    // m_deltaG0 will be used to pass electrochemical potentials    m_rxnstoich.getReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_deltaG0));        for (int i = 0; i < m_ii; i++) {       m_ProdStanConcReac[i] = 1.0;    }    m_rxnstoich.multiplyReactants(DATA_PTR(m_StandardConc), DATA_PTR(m_ProdStanConcReac));  }

void InterpKinetics::update_rates_T(){    if (!m_ntemps) {        rebuildInterpData(200, 250, 3000);    }    double Tnow = thermo().temperature();    if (Tnow >= m_tmax) {        rebuildInterpData(m_ntemps + 10, m_tmin, Tnow + 100.0);    } else if (Tnow <= m_tmin) {        rebuildInterpData(m_ntemps + 2, std::max(Tnow - 20.0, 10.0), m_tmax);    }    size_t n = static_cast<size_t>(floor((Tnow-m_tmin)/(m_tmax-m_tmin)*(m_ntemps-1)));    double dT = Tnow - n * (m_tmax - m_tmin) / (m_ntemps - 1) - m_tmin;    assert(dT >= 0 && dT <= (m_tmax - m_tmin) / (m_ntemps - 1));    assert(n < m_ntemps);    size_t nFall = m_falloff_low_rates.nReactions();    vec_map(&m_rfn[0], nReactions()) = m_rfn_const.col(n) + m_rfn_slope.col(n) * dT;    vec_map(&m_rkcn[0], nReactions()) = m_rkcn_const.col(n) + m_rkcn_slope.col(n) * dT;    vec_map(&m_rfn_low[0], nFall) = m_rfn_low_const.col(n) + m_rfn_low_slope.col(n) * dT;    vec_map(&m_rfn_high[0], nFall) = m_rfn_high_const.col(n) + m_rfn_high_slope.col(n) * dT;    vec_map(&falloff_work[0], nFall) = m_falloff_work_const.col(n) + m_falloff_work_slope.col(n) * dT;    m_logStandConc = log(thermo().standardConcentration());    m_temp = Tnow;    m_ROP_ok = false;}

//====================================================================================================================void GasKinetics::_update_rates_T(){    doublereal T = thermo().temperature();    m_logStandConc = log(thermo().standardConcentration());    doublereal logT = log(T);    if (!m_rfn.empty()) {        m_rates.update(T, logT, &m_rfn[0]);    }    if (!m_rfn_low.empty()) {        m_falloff_low_rates.update(T, logT, &m_rfn_low[0]);        m_falloff_high_rates.update(T, logT, &m_rfn_high[0]);    }    if (!falloff_work.empty()) {        m_falloffn.updateTemp(T, &falloff_work[0]);    }    if (m_plog_rates.nReactions()) {        m_plog_rates.update(T, logT, &m_rfn[0]);    }    if (m_cheb_rates.nReactions()) {        m_cheb_rates.update(T, logT, &m_rfn[0]);    }    m_temp = T;    updateKc();    m_ROP_ok = false;};

//====================================================================================================================void GasKinetics::init(){    m_kk = thermo().nSpecies();    m_rrxn.resize(m_kk);    m_prxn.resize(m_kk);    m_conc.resize(m_kk);    m_grt.resize(m_kk);    m_logp_ref = log(thermo().refPressure()) - log(GasConstant);}

void BulkKinetics::getDeltaSSEnthalpy(doublereal* deltaH){    // Get the standard state enthalpies of the species.    thermo().getEnthalpy_RT(&m_grt[0]);    for (size_t k = 0; k < m_kk; k++) {        m_grt[k] *= thermo().RT();    }    // Use the stoichiometric manager to find deltaH for each reaction.    getReactionDelta(&m_grt[0], deltaH);}

 /**  * This function looks up the string name of a species and  * returns a reference to the ThermoPhase object of the  * phase where the species resides.  * Will throw an error if the species string doesn't match.  */ thermo_t& Kinetics::speciesPhase(std::string nm) {   int np = static_cast<int>(m_thermo.size());   int k;   string id;   for (int n = 0; n < np; n++) {     k = thermo(n).speciesIndex(nm);     if (k >= 0) return thermo(n);   }   throw CanteraError("speciesPhase", "unknown species "+nm);   return thermo(0); }

/** * This routine will look up a species number based on the input * std::string nm. The lookup of species will occur for all phases * listed in the kinetics object. * *  return *   - If a match is found, the position in the species list is returned. *   - If no match is found, the value -1 is returned. * * @param nm   Input string name of the species */size_t Kinetics::kineticsSpeciesIndex(const std::string& nm) const{    for (size_t n = 0; n < m_thermo.size(); n++) {        string id = thermo(n).id();        // Check the ThermoPhase object for a match        size_t k = thermo(n).speciesIndex(nm);        if (k != npos) {            return k + m_start[n];        }    }    return npos;}

void Init_Tisch1(){  int t;  for (t=0; t <= 255;t++) set_rgb_color(t,pal[t].r,pal[t].g,pal[t].b );    Kurven = 0;  Lichter1(250) = 0;  Lichter1(251) = 0;  Lichter1(252) = 0;  Lichter2(247) = 0;  Lichter2(248) = 0;  Lichter2(249) = 0;  Lichter3(244) = 0;  Lichter3(245) = 0;  Lichter3(246) = 0;  Licht4   = 0;  PushUp = TRUE;  Bonus    = 0;  MAXfarbe = 234; // {235-255}  temp = 3;  thermo(temp);  PCSspe[1] = 0;  PCSspe[2] = 0;  PCSspe[3] = 0;  special = 0;  BumpCount = 0;  // Common variables with different values depending on table  tnr = '1';  // {tablenr}  ArmXLinks = 79;  ArmYLinks = 400+135;  ArmXRechts = 159;  ArmYRechts = 400+135;}

/** * Update properties that depend on concentrations. Currently only * the enhanced collision partner concentrations are updated here. */void AqueousKinetics::_update_rates_C(){    thermo().getActivityConcentrations(&m_conc[0]);    m_ROP_ok = false;}

void AqueousKinetics::_update_rates_T(){    doublereal T = thermo().temperature();    m_rates.update(T, log(T), m_rfn.data());    m_temp = T;    updateKc();    m_ROP_ok = false;}

  /**   * Update the equilibrium constants in molar units for all   * reversible reactions. Irreversible reactions have their    * equilibrium constant set to zero.   * For reactions involving charged species the equilibrium    * constant is adjusted according to the electrostatis potential.   */  void InterfaceKinetics::updateKc() {    int i, irxn;    vector_fp& m_rkc = m_kdata->m_rkcn;    fill(m_rkc.begin(), m_rkc.end(), 0.0);    //static vector_fp mu(nTotalSpecies());    if (m_nrev > 0) {      /*       * Get the vector of electrochemical potentials and store it in m_mu0       */      int n, nsp, k, ik = 0;      doublereal rt = GasConstant*thermo(0).temperature();      doublereal rrt = 1.0 / rt;      int np = nPhases();      for (n = 0; n < np; n++) {	thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);	nsp = thermo(n).nSpecies();	for (k = 0; k < nsp; k++) {	  m_mu0[ik] -= rt * thermo(n).logStandardConc(k);	  m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);	  ik++;	}      }      // compute Delta mu^0 for all reversible reactions      m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), 				      DATA_PTR(m_rkc));      for (i = 0; i < m_nrev; i++) {	irxn = m_revindex[i];	if (irxn < 0 || irxn >= nReactions()) {	  throw CanteraError("InterfaceKinetics",			     "illegal value: irxn = "+int2str(irxn));	}	m_rkc[irxn] = exp(m_rkc[irxn]*rrt);      }      for (i = 0; i != m_nirrev; ++i) {	m_rkc[ m_irrev[i] ] = 0.0;      }    }  }

void InterpKinetics::rebuildInterpData(size_t nTemps, double Tmin, double Tmax){    m_ntemps = nTemps;    m_tmin = Tmin;    m_tmax = Tmax;    m_rfn_const = dmatrix::Zero(nReactions(), m_ntemps);    m_rfn_slope = dmatrix::Zero(nReactions(), m_ntemps);    m_rkcn_const = dmatrix::Zero(nReactions(), m_ntemps);    m_rkcn_slope = dmatrix::Zero(nReactions(), m_ntemps);    size_t nFall = m_falloff_low_rates.nReactions();    m_rfn_low_const = dmatrix::Zero(nFall, m_ntemps);    m_rfn_low_slope = dmatrix::Zero(nFall, m_ntemps);    m_rfn_high_const = dmatrix::Zero(nFall, m_ntemps);    m_rfn_high_slope = dmatrix::Zero(nFall, m_ntemps);    m_falloff_work_const = dmatrix::Zero(nFall, m_ntemps);    m_falloff_work_slope = dmatrix::Zero(nFall, m_ntemps);    double T_save = thermo().temperature();    double rho_save = thermo().density();    for (size_t n = 0; n < m_ntemps; n++) {        thermo().setState_TR(Tmin + ((double) n)/(m_ntemps - 1) * (Tmax - Tmin),                             rho_save);        GasKinetics::update_rates_T();        m_rfn_const.col(n) = vec_map(&m_rfn[0], nReactions());        m_rkcn_const.col(n) = vec_map(&m_rkcn[0], nReactions());        m_rfn_low_const.col(n) = vec_map(&m_rfn_low[0], nFall);        m_rfn_high_const.col(n) = vec_map(&m_rfn_high[0], nFall);        m_falloff_work_const.col(n) = vec_map(&falloff_work[0], nFall);    }    double dT = (Tmax - Tmin) / (m_ntemps - 1);    for (size_t n = 0; n < m_ntemps - 1; n++) {        m_rfn_slope.col(n) = (m_rfn_const.col(n+1) - m_rfn_const.col(n)) / dT;        m_rkcn_slope.col(n) = (m_rkcn_const.col(n+1) - m_rkcn_const.col(n)) / dT;        m_rfn_low_slope.col(n) = (m_rfn_low_const.col(n+1) - m_rfn_low_const.col(n)) / dT;        m_rfn_high_slope.col(n) = (m_rfn_high_const.col(n+1) - m_rfn_high_const.col(n)) / dT;        m_falloff_work_slope.col(n) = (m_falloff_work_const.col(n+1) - m_falloff_work_const.col(n)) / dT;    }    thermo().setState_TR(T_save, rho_save);}

/** * This routine will look up a species number based on the input * std::string nm. The lookup of species will occur in the specified * phase of the object, or all phases if ph is "<any>". * *  return *   - If a match is found, the position in the species list is returned. *   - If no match is found, the value npos (-1) is returned. * * @param nm   Input string name of the species * @param ph   Input string name of the phase. */size_t Kinetics::kineticsSpeciesIndex(const std::string& nm,                                      const std::string& ph) const{    if (ph == "<any>") {        return kineticsSpeciesIndex(nm);    }    for (size_t n = 0; n < m_thermo.size(); n++) {        string id = thermo(n).id();        if (ph == id) {            size_t k = thermo(n).speciesIndex(nm);            if (k == npos) {                return npos;            }            return k + m_start[n];        }    }    return npos;}

  void AntiochMixture<KineticsThermoCurveFit>::build_stat_mech_ref_correction()  {    Antioch::StatMechThermodynamics<libMesh::Real> thermo( *(this->_antioch_gas.get()) );    _h_stat_mech_ref_correction.resize(this->n_species());    for( unsigned int s = 0; s < this->n_species(); s++ )      {        _h_stat_mech_ref_correction[s] = -thermo.h_tot( s, 298.15 ) + thermo.e_0(s);      }  }