<commit_before>/**
 *  @file GasKinetics.cpp
 *
 * Homogeneous kinetics in ideal gases
 *
 */

// Copyright 2001  California Institute of Technology

#include "cantera/kinetics/GasKinetics.h"

#include "cantera/kinetics/ReactionData.h"
#include "cantera/kinetics/Enhanced3BConc.h"
#include "cantera/kinetics/ThirdBodyMgr.h"
#include "cantera/kinetics/RateCoeffMgr.h"

#include <iostream>
using namespace std;

namespace Cantera
{

//====================================================================================================================
/*
 * Construct an empty reaction mechanism.
 */
GasKinetics::
GasKinetics(thermo_t* thermo) :
    Kinetics(),
    m_nfall(0),
    m_nirrev(0),
    m_nrev(0),
    m_logp_ref(0.0),
    m_logc_ref(0.0),
    m_logStandConc(0.0),
    m_ROP_ok(false),
    m_temp(0.0),
    m_finalized(false)
{
    if (thermo != 0) {
        addPhase(*thermo);
    }
    m_temp = 0.0;
}

//====================================================================================================================
GasKinetics::GasKinetics(const GasKinetics& right) :
    Kinetics(),
    m_nfall(0),
    m_nirrev(0),
    m_nrev(0),
    m_logp_ref(0.0),
    m_logc_ref(0.0),
    m_logStandConc(0.0),
    m_ROP_ok(false),
    m_temp(0.0),
    m_finalized(false)
{
    m_temp = 0.0;
    *this = right;
}
//====================================================================================================================
GasKinetics::~GasKinetics()
{
}
//====================================================================================================================
GasKinetics& GasKinetics::operator=(const GasKinetics& right)
{
    if (this == &right) {
        return *this;
    }

    Kinetics::operator=(right);

    m_nfall = right.m_nfall;
    m_fallindx = right.m_fallindx;
    m_falloff_low_rates = right.m_falloff_low_rates;
    m_falloff_high_rates = right.m_falloff_high_rates;
    m_rates = right.m_rates;
    m_index = right.m_index;
    m_falloffn = right.m_falloffn;
    m_3b_concm = right.m_3b_concm;
    m_falloff_concm = right.m_falloff_concm;
    m_irrev = right.m_irrev;
    m_plog_rates = right.m_plog_rates;
    m_cheb_rates = right.m_cheb_rates;

    m_rxnstoich = right.m_rxnstoich;

    m_fwdOrder = right.m_fwdOrder;
    m_nirrev = right.m_nirrev;
    m_nrev = right.m_nrev;
    m_rgroups = right.m_rgroups;
    m_pgroups = right.m_pgroups;
    m_rxntype = right.m_rxntype;
    m_rrxn = right.m_rrxn;
    m_prxn = right.m_prxn;
    m_dn = right.m_dn;
    m_revindex = right.m_revindex;
    m_rxneqn = right.m_rxneqn;

    m_logp_ref = right.m_logp_ref;
    m_logc_ref  = right.m_logc_ref;
    m_logStandConc = right.m_logStandConc;
    m_ropf = right.m_ropf;
    m_ropr = right.m_ropr;
    m_ropnet = right.m_ropnet;
    m_rfn_low = right.m_rfn_low;
    m_rfn_high = right.m_rfn_high;
    m_ROP_ok  = right.m_ROP_ok;
    m_temp = right.m_temp;
    m_rfn  = right.m_rfn;
    falloff_work = right.falloff_work;
    concm_3b_values = right.concm_3b_values;
    concm_falloff_values = right.concm_falloff_values;
    m_rkcn = right.m_rkcn;

    m_conc = right.m_conc;
    m_grt = right.m_grt;
    m_finalized = right.m_finalized;

    throw CanteraError("GasKinetics::operator=()",
                       "Unfinished implementation");

    return *this;
}
//====================================================================================================================
// Duplication routine for objects which inherit from Kinetics
/*
 *  This virtual routine can be used to duplicate %Kinetics objects
 *  inherited from %Kinetics even if the application only has
 *  a pointer to %Kinetics to work with.
 *
 *  These routines are basically wrappers around the derived copy
 *  constructor.
 *
 * @param  tpVector Vector of shallow pointers to ThermoPhase objects. this is the
 *                  m_thermo vector within this object
 */
Kinetics* GasKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
    GasKinetics* gK = new GasKinetics(*this);
    gK->assignShallowPointers(tpVector);
    return gK;
}
//====================================================================================================================
/**
 * Update temperature-dependent portions of reaction rates and
 * falloff functions.
 */
void GasKinetics::update_T()
{
}
//====================================================================================================================
void GasKinetics::
update_C() {}
//====================================================================================================================
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::
_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;
}
//====================================================================================================================
/**
 * Update the equilibrium constants in molar units.
 */
void GasKinetics::updateKc()
{
    thermo().getStandardChemPotentials(&m_grt[0]);
    fill(m_rkcn.begin(), m_rkcn.end(), 0.0);

    // compute Delta G^0 for all reversible reactions
    m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);

    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 - m_dn[irxn]*m_logStandConc);
    }

    for (size_t i = 0; i != m_nirrev; ++i) {
        m_rkcn[ m_irrev[i] ] = 0.0;
    }
}
//====================================================================================================================
/**
 * Get the equilibrium constants of all reactions, whether
 * reversible or not.
 */
void GasKinetics::getEquilibriumConstants(doublereal* kc)
{
    _update_rates_T();
    thermo().getStandardChemPotentials(&m_grt[0]);
    fill(m_rkcn.begin(), m_rkcn.end(), 0.0);

    // compute Delta G^0 for all reactions
    m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);

    doublereal rrt = 1.0/(GasConstant * thermo().temperature());
    for (size_t i = 0; i < m_ii; i++) {
        kc[i] = exp(-m_rkcn[i]*rrt + m_dn[i]*m_logStandConc);
    }
<commit_msg>[Kinetics] Ensure that equilibrium constants are not zero<commit_after>/**
 *  @file GasKinetics.cpp
 *
 * Homogeneous kinetics in ideal gases
 *
 */

// Copyright 2001  California Institute of Technology

#include "cantera/kinetics/GasKinetics.h"

#include "cantera/kinetics/ReactionData.h"
#include "cantera/kinetics/Enhanced3BConc.h"
#include "cantera/kinetics/ThirdBodyMgr.h"
#include "cantera/kinetics/RateCoeffMgr.h"

#include <iostream>
using namespace std;

namespace Cantera
{

//====================================================================================================================
/*
 * Construct an empty reaction mechanism.
 */
GasKinetics::
GasKinetics(thermo_t* thermo) :
    Kinetics(),
    m_nfall(0),
    m_nirrev(0),
    m_nrev(0),
    m_logp_ref(0.0),
    m_logc_ref(0.0),
    m_logStandConc(0.0),
    m_ROP_ok(false),
    m_temp(0.0),
    m_finalized(false)
{
    if (thermo != 0) {
        addPhase(*thermo);
    }
    m_temp = 0.0;
}

//====================================================================================================================
GasKinetics::GasKinetics(const GasKinetics& right) :
    Kinetics(),
    m_nfall(0),
    m_nirrev(0),
    m_nrev(0),
    m_logp_ref(0.0),
    m_logc_ref(0.0),
    m_logStandConc(0.0),
    m_ROP_ok(false),
    m_temp(0.0),
    m_finalized(false)
{
    m_temp = 0.0;
    *this = right;
}
//====================================================================================================================
GasKinetics::~GasKinetics()
{
}
//====================================================================================================================
GasKinetics& GasKinetics::operator=(const GasKinetics& right)
{
    if (this == &right) {
        return *this;
    }

    Kinetics::operator=(right);

    m_nfall = right.m_nfall;
    m_fallindx = right.m_fallindx;
    m_falloff_low_rates = right.m_falloff_low_rates;
    m_falloff_high_rates = right.m_falloff_high_rates;
    m_rates = right.m_rates;
    m_index = right.m_index;
    m_falloffn = right.m_falloffn;
    m_3b_concm = right.m_3b_concm;
    m_falloff_concm = right.m_falloff_concm;
    m_irrev = right.m_irrev;
    m_plog_rates = right.m_plog_rates;
    m_cheb_rates = right.m_cheb_rates;

    m_rxnstoich = right.m_rxnstoich;

    m_fwdOrder = right.m_fwdOrder;
    m_nirrev = right.m_nirrev;
    m_nrev = right.m_nrev;
    m_rgroups = right.m_rgroups;
    m_pgroups = right.m_pgroups;
    m_rxntype = right.m_rxntype;
    m_rrxn = right.m_rrxn;
    m_prxn = right.m_prxn;
    m_dn = right.m_dn;
    m_revindex = right.m_revindex;
    m_rxneqn = right.m_rxneqn;

    m_logp_ref = right.m_logp_ref;
    m_logc_ref  = right.m_logc_ref;
    m_logStandConc = right.m_logStandConc;
    m_ropf = right.m_ropf;
    m_ropr = right.m_ropr;
    m_ropnet = right.m_ropnet;
    m_rfn_low = right.m_rfn_low;
    m_rfn_high = right.m_rfn_high;
    m_ROP_ok  = right.m_ROP_ok;
    m_temp = right.m_temp;
    m_rfn  = right.m_rfn;
    falloff_work = right.falloff_work;
    concm_3b_values = right.concm_3b_values;
    concm_falloff_values = right.concm_falloff_values;
    m_rkcn = right.m_rkcn;

    m_conc = right.m_conc;
    m_grt = right.m_grt;
    m_finalized = right.m_finalized;

    throw CanteraError("GasKinetics::operator=()",
                       "Unfinished implementation");

    return *this;
}
//====================================================================================================================
// Duplication routine for objects which inherit from Kinetics
/*
 *  This virtual routine can be used to duplicate %Kinetics objects
 *  inherited from %Kinetics even if the application only has
 *  a pointer to %Kinetics to work with.
 *
 *  These routines are basically wrappers around the derived copy
 *  constructor.
 *
 * @param  tpVector Vector of shallow pointers to ThermoPhase objects. this is the
 *                  m_thermo vector within this object
 */
Kinetics* GasKinetics::duplMyselfAsKinetics(const std::vector<thermo_t*> & tpVector) const
{
    GasKinetics* gK = new GasKinetics(*this);
    gK->assignShallowPointers(tpVector);
    return gK;
}
//====================================================================================================================
/**
 * Update temperature-dependent portions of reaction rates and
 * falloff functions.
 */
void GasKinetics::update_T()
{
}
//====================================================================================================================
void GasKinetics::
update_C() {}
//====================================================================================================================
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::
_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;
}
//====================================================================================================================
/**
 * Update the equilibrium constants in molar units.
 */
void GasKinetics::updateKc()
{
    thermo().getStandardChemPotentials(&m_grt[0]);
    fill(m_rkcn.begin(), m_rkcn.end(), 0.0);

    // compute Delta G^0 for all reversible reactions
    m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);

    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] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc),
                                BigNumber);
    }

    for (size_t i = 0; i != m_nirrev; ++i) {
        m_rkcn[ m_irrev[i] ] = 0.0;
    }
}
//====================================================================================================================
/**
 * Get the equilibrium constants of all reactions, whether
 * reversible or not.
 */
void GasKinetics::getEquilibriumConstants(doublereal* kc)
{
    _update_rates_T();
    thermo().getStandardChemPotentials(&m_grt[0]);
    fill(m_rkcn.begin(), m_rkcn.end(), 0.0);

    // compute Delta G^0 for all reactions
    m_rxnstoich.getReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]);

    doublereal rrt = 1.0/(GasConstant * thermo().temperature());
    for (size_t i = 0; i < m_ii; i++) {
        kc[i] = exp(-m_rkcn[i]*rrt + m_dn[i]*m_logStandConc);
    }

{'fraction_non_alphanumeric': 0.2770661291371351, 'fraction_numerical': 0.00988360751674361, 'mean_word_length': 3.9789575914535447, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5189606', 'n_tokens_mistral': 5595, 'n_tokens_neox': 5186, 'n_words': 876}