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

// Copyright 2001  California Institute of Technology


// turn off warnings under Windows
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif

#include "GasKinetics.h"

#include "ReactionData.h"
#include "Enhanced3BConc.h"
#include "ThirdBodyMgr.h"
#include "RateCoeffMgr.h"

//#include "../user/grirxnstoich.h"

#include <iostream>
using namespace std;


namespace Cantera {


    /**
     * Construct an empty reaction mechanism.
     */    
    GasKinetics::
    GasKinetics(thermo_t* thermo) :
        Kinetics(),
        m_kk(0), 
        m_nfall(0), 
        m_dt_threshold(0.0), // 1.e-6),
        m_nirrev(0), 
        m_nrev(0),
        m_finalized(false)
    {
        if (thermo != 0) addPhase(*thermo);
        m_kdata = new GasKineticsData;
        m_kdata->m_temp = 0.0;
        m_rxnstoich = new ReactionStoichMgr;
    }

  GasKinetics::
  ~GasKinetics() {delete m_kdata; delete m_rxnstoich;}

    /**
     * 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_kdata->m_logStandConc = log(thermo().standardConcentration()); 
        if (fabs(T - m_kdata->m_temp) > 0.0) { // m_dt_threshold) {
            doublereal logT = log(T);
            //m_kdata->m_logp0 - logT;
            m_rates.update(T, logT, m_kdata->m_rfn.begin());
            m_falloff_low_rates.update(T, logT, m_kdata->m_rfn_low.begin()); 
            m_falloff_high_rates.update(T, logT, m_kdata->m_rfn_high.begin());
            m_falloffn.updateTemp(T, m_kdata->falloff_work.begin());
            m_kdata->m_temp = T;
            updateKc();
            m_kdata->m_ROP_ok = false;
        }
    };


    /**
     * Update properties that depend on concentrations. Currently only
     * the enhanced collision partner concentrations are updated here.
     */         
    void GasKinetics::
    _update_rates_C() {
        thermo().getActivityConcentrations(m_conc.begin());
        doublereal ctot = thermo().molarDensity();
        m_3b_concm.update(m_conc, ctot, m_kdata->concm_3b_values.begin());
        m_falloff_concm.update(m_conc, ctot, 
            m_kdata->concm_falloff_values.begin());
        m_kdata->m_ROP_ok = false;
    }

    /**
     * Update the equilibrium constants in molar units.
     */
    void GasKinetics::updateKc() {
        int i, irxn;
        vector_fp& m_rkc = m_kdata->m_rkcn;
        
        thermo().getStandardChemPotentials(m_grt.begin());
        fill(m_rkc.begin(), m_rkc.end(), 0.0);

        // compute Delta G^0 for all reversible reactions
        m_rxnstoich->getRevReactionDelta(m_ii, m_grt.begin(), m_rkc.begin());
 
        doublereal logStandConc = m_kdata->m_logStandConc;
        doublereal rrt = 1.0/(GasConstant * thermo().temperature());
        for (i = 0; i < m_nrev; i++) {
            irxn = m_revindex[i];
            m_rkc[irxn] = exp(m_rkc[irxn]*rrt - m_dn[irxn]*logStandConc);
        }

        for(i = 0; i != m_nirrev; ++i) {
            m_rkc[ m_irrev[i] ] = 0.0;
        }
    }

    /**
     * Get the equilibrium constants of all reactions, whether
     * reversible or not.
     */
    void GasKinetics::getEquilibriumConstants(doublereal* kc) {
        int i;
        _update_rates_T();
        vector_fp& rkc = m_kdata->m_rkcn;
        //thermo().getGibbs_RT(m_grt.begin());
        thermo().getStandardChemPotentials(m_grt.begin());
        fill(rkc.begin(), rkc.end(), 0.0);
        
        // compute Delta G^0 for all reactions
        m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), rkc.begin());
 
        doublereal logStandConc = m_kdata->m_logStandConc;
        doublereal rrt = 1.0/(GasConstant * thermo().temperature());
        for (i = 0; i < m_ii; i++) {
            kc[i] = exp(-rkc[i]*rrt + m_dn[i]*logStandConc);
        }
    }

    /**
     *
     * getDeltaGibbs():
     *
     * Return the vector of values for the reaction gibbs free energy
     * change
     * These values depend upon the concentration
     * of the ideal gas.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaGibbs(doublereal* deltaG) {
	/*
	 * Get the chemical potentials of the species in the 
	 * ideal gas solution.
	 */
	thermo().getChemPotentials(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaG);
    }
    
    /**
     *
     * getDeltaEnthalpy():
     * 
     * Return the vector of values for the reactions change in
     * enthalpy.
     * These values depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaEnthalpy(doublereal* deltaH) {
	/*
	 * Get the partial molar enthalpy of all species in the 
	 * ideal gas.
	 */
	thermo().getPartialMolarEnthalpies(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaH);
    }

    /************************************************************************
     *
     * getDeltaEntropy():
     *
     * Return the vector of values for the reactions change in
     * entropy.
     * These values depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1 Kelvin-1
     */
    void GasKinetics::getDeltaEntropy( doublereal* deltaS) {
	/*
	 * Get the partial molar entropy of all species in the
	 * solid solution.
	 */
	thermo().getPartialMolarEntropies(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaS for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaS);
    }

    /**
     *
     * getDeltaSSGibbs():
     *
     * Return the vector of values for the reaction 
     * standard state gibbs free energy change.
     * These values don't depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaSSGibbs(doublereal* deltaG) {
	/*
	 *  Get the standard state chemical potentials of the species.
         *  This is the array of chemical potentials at unit activity 
	 *  We define these here as the chemical potentials of the pure
	 *  species at the temperature and pressure of the solution.
	 */
        thermo().getStandardChemPotentials(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaG);
    }

    /**
     *
     * getDeltaSSEnthalpy():
     *
     * Return the vector of values for the change in the
     * standard state enthalpies of reaction.
     * These values don't depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH) {
	/*
	 *  Get the standard state enthalpies of the species.
         *  This is the array of chemical potentials at unit activity 
	 *  We define these here as the enthalpies of the pure
	 *  species at the temperature and pressure of the solution.
	 */
	thermo().getEnthalpy_RT(m_grt.begin());
	doublereal RT = thermo().temperature() * GasConstant;
	for (int k = 0; k < m_kk; k++) {
	  m_grt[k] *= RT;
	}
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaH);
    }

    /*********************************************************************
     *
     * getDeltaSSEntropy():
     *
     * Return the vector of values for the change in the
     * standard state entropies for each reaction.
     * These values don't depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1 Kelvin-1
     */
    void GasKinetics::getDeltaSSEntropy(doublereal* deltaS) {
	/*
	 *  Get the standard state entropy of the species.
	 *  We define these here as the entropies of the pure
	 *  species at the temperature and pressure of the solution.
	 */
	thermo().getEntropy_R(m_grt.begin());
	doublereal R = GasConstant;
	for (int k = 0; k < m_kk; k++) {
	  m_grt[k] *= R;
	}
	/*
	 * Use the stoichiometric manager to find deltaS for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaS);
    }

    void GasKinetics::processFalloffReactions() {

        int i;
        const vector_fp& fc = m_kdata->concm_falloff_values;
        const array_fp& m_rf_low = m_kdata->m_rfn_low;
        const array_fp& m_rf_high = m_kdata->m_rfn_high;

        // use m_ropr for temporary storage of reduced pressure
        array_fp& pr = m_kdata->m_ropr;

        array_fp& ropf = m_kdata->m_ropf;

        for (i = 0; i < m_nfall; i++) {
            pr[i] = fc[i] * m_rf_low[i] / m_rf_high[i];
        }

        m_falloffn.pr_to_falloff( pr.begin(), m_kdata->falloff_work.begin() );
        
        for (i = 0; i < m_nfall; i++) {
            pr[i] *= m_rf_high[i]; 
        }

        scatter_copy(pr.begin(), pr.begin() + m_nfall, 
            ropf.begin(), m_fallindx.begin());
    }


    void GasKinetics::updateROP() {

        _update_rates_T();
        _update_rates_C();

        if (m_kdata->m_ROP_ok) return;

        const vector_fp& rf = m_kdata->m_rfn;
        const vector_fp& m_rkc = m_kdata->m_rkcn;
        array_fp& ropf = m_kdata->m_ropf;
        array_fp& ropr = m_kdata->m_ropr;
        array_fp& ropnet = m_kdata->m_ropnet;

        // copy rate coefficients into ropf
        copy(rf.begin(), rf.end(), ropf.begin());

        // multiply ropf by enhanced 3b conc for all 3b rxns
        m_3b_concm.multiply( ropf.begin(), m_kdata->concm_3b_values.begin() );

        processFalloffReactions();

        // multiply by perturbation factor
        multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
           
        // copy the forward rates to the reverse rates                
        copy(ropf.begin(), ropf.end(), ropr.begin());
        
        // for reverse rates computed from thermochemistry, multiply
        // the forward rates copied into m_ropr by the reciprocals of
        // the equilibrium constants
        multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());

        // multiply ropf by concentration products
        m_rxnstoich->multiplyReactants(m_conc.begin(), ropf.begin()); 
        //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); 

        // for reversible reactions, multiply ropr by concentration
        // products
        m_rxnstoich->multiplyRevProducts(m_conc.begin(), ropr.begin()); 
        //m_revProductStoich.multiply(m_conc.begin(), ropr.begin());

        for (int j = 0; j != m_ii; ++j) {
            ropnet[j] = ropf[j] - ropr[j];
        }

        m_kdata->m_ROP_ok = true;
    }

    /**
     *
     * getFwdRateConstants():
     *
     * Update the rate of progress for the reactions.
     * This key routine makes sure that the rate of progress vectors
     * located in the solid kinetics data class are up to date.
     */
    void GasKinetics::
    getFwdRateConstants(doublereal *kfwd) {
        _update_rates_T();
	_update_rates_C();

	// copy rate coefficients into ropf
	const vector_fp& rf = m_kdata->m_rfn;
	array_fp& ropf = m_kdata->m_ropf;
        copy(rf.begin(), rf.end(), ropf.begin());

        // multiply ropf by enhanced 3b conc for all 3b rxns
        m_3b_concm.multiply(ropf.begin(), m_kdata->concm_3b_values.begin() );

	/*
	 * This routine is hardcoded to replace some of the values
	 * of the ropf vector.
	 */
        processFalloffReactions();

        // multiply by perturbation factor
        multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
       
	for (int i = 0; i < m_ii; i++) {
	  kfwd[i] = ropf[i];
	}
    }

    /**
     *
     * getRevRateConstants():
     *
     * Return a vector of the reverse reaction rate constants
     *
     * Length is the number of reactions. units depends
     * on many issues. Note, this routine will return rate constants
     * for irreversible reactions if the default for
     * doIrreversible is overridden.
     */
    void GasKinetics::
    getRevRateConstants(doublereal *krev, bool doIrreversible) {
	/*
	 * go get the forward rate constants. -> note, we don't
	 * really care about speed or redundancy in these
	 * informational routines.
	 */
        getFwdRateConstants(krev);

	if (doIrreversible) {
	  doublereal *tmpKc = m_kdata->m_ropnet.begin();
	  getEquilibriumConstants(tmpKc);
	  for (int i = 0; i < m_ii; i++) {
	    krev[i] /=  tmpKc[i];
	  }
	} else {
	  /*
	   * m_rkc[] is zero for irreversibly reactions
	   */
	  const vector_fp& m_rkc = m_kdata->m_rkcn;
	  for (int i = 0; i < m_ii; i++) {
	    krev[i] *= m_rkc[i];
	  }
	}
    }
<commit_msg>added statements to getEquilibriumConstants to force a T update on the next call. Array m_rkcn was being left containing incorrect data. Problem is now fixed.<commit_after>/**
 *  @file GasKinetics.cpp 
 *
 * Homogeneous kinetics in ideal gases
 * 
 */

// Copyright 2001  California Institute of Technology


// turn off warnings under Windows
#ifdef WIN32
#pragma warning(disable:4786)
#pragma warning(disable:4503)
#endif

#include "GasKinetics.h"

#include "ReactionData.h"
#include "Enhanced3BConc.h"
#include "ThirdBodyMgr.h"
#include "RateCoeffMgr.h"

//#include "../user/grirxnstoich.h"

#include <iostream>
using namespace std;


namespace Cantera {


    /**
     * Construct an empty reaction mechanism.
     */    
    GasKinetics::
    GasKinetics(thermo_t* thermo) :
        Kinetics(),
        m_kk(0), 
        m_nfall(0), 
        m_dt_threshold(0.0), // 1.e-6),
        m_nirrev(0), 
        m_nrev(0),
        m_finalized(false)
    {
        if (thermo != 0) addPhase(*thermo);
        m_kdata = new GasKineticsData;
        m_kdata->m_temp = 0.0;
        m_rxnstoich = new ReactionStoichMgr;
    }

  GasKinetics::
  ~GasKinetics() {delete m_kdata; delete m_rxnstoich;}

    /**
     * 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_kdata->m_logStandConc = log(thermo().standardConcentration()); 
        if (fabs(T - m_kdata->m_temp) > 0.0) { // m_dt_threshold) {
            doublereal logT = log(T);
            //m_kdata->m_logp0 - logT;
            m_rates.update(T, logT, m_kdata->m_rfn.begin());
            m_falloff_low_rates.update(T, logT, m_kdata->m_rfn_low.begin()); 
            m_falloff_high_rates.update(T, logT, m_kdata->m_rfn_high.begin());
            m_falloffn.updateTemp(T, m_kdata->falloff_work.begin());
            m_kdata->m_temp = T;
            updateKc();
            m_kdata->m_ROP_ok = false;
        }
    };


    /**
     * Update properties that depend on concentrations. Currently only
     * the enhanced collision partner concentrations are updated here.
     */         
    void GasKinetics::
    _update_rates_C() {
        thermo().getActivityConcentrations(m_conc.begin());
        doublereal ctot = thermo().molarDensity();
        m_3b_concm.update(m_conc, ctot, m_kdata->concm_3b_values.begin());
        m_falloff_concm.update(m_conc, ctot, 
            m_kdata->concm_falloff_values.begin());
        m_kdata->m_ROP_ok = false;
    }

    /**
     * Update the equilibrium constants in molar units.
     */
    void GasKinetics::updateKc() {
        int i, irxn;
        vector_fp& m_rkc = m_kdata->m_rkcn;
        
        thermo().getStandardChemPotentials(m_grt.begin());
        fill(m_rkc.begin(), m_rkc.end(), 0.0);

        // compute Delta G^0 for all reversible reactions
        m_rxnstoich->getRevReactionDelta(m_ii, m_grt.begin(), m_rkc.begin());
 
        doublereal logStandConc = m_kdata->m_logStandConc;
        doublereal rrt = 1.0/(GasConstant * thermo().temperature());
        for (i = 0; i < m_nrev; i++) {
            irxn = m_revindex[i];
            m_rkc[irxn] = exp(m_rkc[irxn]*rrt - m_dn[irxn]*logStandConc);
        }

        for(i = 0; i != m_nirrev; ++i) {
            m_rkc[ m_irrev[i] ] = 0.0;
        }
    }

    /**
     * Get the equilibrium constants of all reactions, whether
     * reversible or not.
     */
    void GasKinetics::getEquilibriumConstants(doublereal* kc) {
        int i;
        _update_rates_T();
        vector_fp& rkc = m_kdata->m_rkcn;
        //thermo().getGibbs_RT(m_grt.begin());
        thermo().getStandardChemPotentials(m_grt.begin());
        fill(rkc.begin(), rkc.end(), 0.0);
        
        // compute Delta G^0 for all reactions
        m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), rkc.begin());
 
        doublereal logStandConc = m_kdata->m_logStandConc;
        doublereal rrt = 1.0/(GasConstant * thermo().temperature());
        for (i = 0; i < m_ii; i++) {
            kc[i] = exp(-rkc[i]*rrt + m_dn[i]*logStandConc);
        }

        // force an update of T-dependent properties, so that m_rkcn will
        // be updated before it is used next.
        m_kdata->m_temp = 0.0;
    }

    /**
     *
     * getDeltaGibbs():
     *
     * Return the vector of values for the reaction gibbs free energy
     * change
     * These values depend upon the concentration
     * of the ideal gas.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaGibbs(doublereal* deltaG) {
	/*
	 * Get the chemical potentials of the species in the 
	 * ideal gas solution.
	 */
	thermo().getChemPotentials(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaG);
    }
    
    /**
     *
     * getDeltaEnthalpy():
     * 
     * Return the vector of values for the reactions change in
     * enthalpy.
     * These values depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaEnthalpy(doublereal* deltaH) {
	/*
	 * Get the partial molar enthalpy of all species in the 
	 * ideal gas.
	 */
	thermo().getPartialMolarEnthalpies(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaH);
    }

    /************************************************************************
     *
     * getDeltaEntropy():
     *
     * Return the vector of values for the reactions change in
     * entropy.
     * These values depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1 Kelvin-1
     */
    void GasKinetics::getDeltaEntropy( doublereal* deltaS) {
	/*
	 * Get the partial molar entropy of all species in the
	 * solid solution.
	 */
	thermo().getPartialMolarEntropies(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaS for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaS);
    }

    /**
     *
     * getDeltaSSGibbs():
     *
     * Return the vector of values for the reaction 
     * standard state gibbs free energy change.
     * These values don't depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaSSGibbs(doublereal* deltaG) {
	/*
	 *  Get the standard state chemical potentials of the species.
         *  This is the array of chemical potentials at unit activity 
	 *  We define these here as the chemical potentials of the pure
	 *  species at the temperature and pressure of the solution.
	 */
        thermo().getStandardChemPotentials(m_grt.begin());
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaG);
    }

    /**
     *
     * getDeltaSSEnthalpy():
     *
     * Return the vector of values for the change in the
     * standard state enthalpies of reaction.
     * These values don't depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1
     */
    void GasKinetics::getDeltaSSEnthalpy(doublereal* deltaH) {
	/*
	 *  Get the standard state enthalpies of the species.
         *  This is the array of chemical potentials at unit activity 
	 *  We define these here as the enthalpies of the pure
	 *  species at the temperature and pressure of the solution.
	 */
	thermo().getEnthalpy_RT(m_grt.begin());
	doublereal RT = thermo().temperature() * GasConstant;
	for (int k = 0; k < m_kk; k++) {
	  m_grt[k] *= RT;
	}
	/*
	 * Use the stoichiometric manager to find deltaG for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaH);
    }

    /*********************************************************************
     *
     * getDeltaSSEntropy():
     *
     * Return the vector of values for the change in the
     * standard state entropies for each reaction.
     * These values don't depend upon the concentration
     * of the solution.
     *
     *  units = J kmol-1 Kelvin-1
     */
    void GasKinetics::getDeltaSSEntropy(doublereal* deltaS) {
	/*
	 *  Get the standard state entropy of the species.
	 *  We define these here as the entropies of the pure
	 *  species at the temperature and pressure of the solution.
	 */
	thermo().getEntropy_R(m_grt.begin());
	doublereal R = GasConstant;
	for (int k = 0; k < m_kk; k++) {
	  m_grt[k] *= R;
	}
	/*
	 * Use the stoichiometric manager to find deltaS for each
	 * reaction.
	 */
	m_rxnstoich->getReactionDelta(m_ii, m_grt.begin(), deltaS);
    }

    void GasKinetics::processFalloffReactions() {

        int i;
        const vector_fp& fc = m_kdata->concm_falloff_values;
        const array_fp& m_rf_low = m_kdata->m_rfn_low;
        const array_fp& m_rf_high = m_kdata->m_rfn_high;

        // use m_ropr for temporary storage of reduced pressure
        array_fp& pr = m_kdata->m_ropr;

        array_fp& ropf = m_kdata->m_ropf;

        for (i = 0; i < m_nfall; i++) {
            pr[i] = fc[i] * m_rf_low[i] / m_rf_high[i];
        }

        m_falloffn.pr_to_falloff( pr.begin(), m_kdata->falloff_work.begin() );
        
        for (i = 0; i < m_nfall; i++) {
            pr[i] *= m_rf_high[i]; 
        }

        scatter_copy(pr.begin(), pr.begin() + m_nfall, 
            ropf.begin(), m_fallindx.begin());
    }


    void GasKinetics::updateROP() {

        _update_rates_T();
        _update_rates_C();

        if (m_kdata->m_ROP_ok) return;

        const vector_fp& rf = m_kdata->m_rfn;
        const vector_fp& m_rkc = m_kdata->m_rkcn;
        array_fp& ropf = m_kdata->m_ropf;
        array_fp& ropr = m_kdata->m_ropr;
        array_fp& ropnet = m_kdata->m_ropnet;

        // copy rate coefficients into ropf
        copy(rf.begin(), rf.end(), ropf.begin());

        // multiply ropf by enhanced 3b conc for all 3b rxns
        m_3b_concm.multiply( ropf.begin(), m_kdata->concm_3b_values.begin() );

        processFalloffReactions();

        // multiply by perturbation factor
        multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
           
        // copy the forward rates to the reverse rates                
        copy(ropf.begin(), ropf.end(), ropr.begin());
        
        // for reverse rates computed from thermochemistry, multiply
        // the forward rates copied into m_ropr by the reciprocals of
        // the equilibrium constants
        multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());

        // multiply ropf by concentration products
        m_rxnstoich->multiplyReactants(m_conc.begin(), ropf.begin()); 
        //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); 

        // for reversible reactions, multiply ropr by concentration
        // products
        m_rxnstoich->multiplyRevProducts(m_conc.begin(), ropr.begin()); 
        //m_revProductStoich.multiply(m_conc.begin(), ropr.begin());

        for (int j = 0; j != m_ii; ++j) {
            ropnet[j] = ropf[j] - ropr[j];
        }

        m_kdata->m_ROP_ok = true;
    }

    /**
     *
     * getFwdRateConstants():
     *
     * Update the rate of progress for the reactions.
     * This key routine makes sure that the rate of progress vectors
     * located in the solid kinetics data class are up to date.
     */
    void GasKinetics::
    getFwdRateConstants(doublereal *kfwd) {
        _update_rates_T();
	_update_rates_C();

	// copy rate coefficients into ropf
	const vector_fp& rf = m_kdata->m_rfn;
	array_fp& ropf = m_kdata->m_ropf;
        copy(rf.begin(), rf.end(), ropf.begin());

        // multiply ropf by enhanced 3b conc for all 3b rxns
        m_3b_concm.multiply(ropf.begin(), m_kdata->concm_3b_values.begin() );

	/*
	 * This routine is hardcoded to replace some of the values
	 * of the ropf vector.
	 */
        processFalloffReactions();

        // multiply by perturbation factor
        multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
       
	for (int i = 0; i < m_ii; i++) {
	  kfwd[i] = ropf[i];
	}
    }

    /**
     *
     * getRevRateConstants():
     *
     * Return a vector of the reverse reaction rate constants
     *
     * Length is the number of reactions. units depends
     * on many issues. Note, this routine will return rate constants
     * for irreversible reactions if the default for
     * doIrreversible is overridden.
     */
    void GasKinetics::
    getRevRateConstants(doublereal *krev, bool doIrreversible) {
	/*
	 * go get the forward rate constants. -> note, we don't
	 * really care about speed or redundancy in these
	 * informational routines.
	 */
        getFwdRateConstants(krev);

	if (doIrreversible) {
	  doublereal *tmpKc = m_kdata->m_ropnet.begin();
	  getEquilibriumConstants(tmpKc);
	  for (int i = 0; i < m_ii; i++) {
	    krev[i] /=  tmpKc[i];
	  }
	} else {
	  /*
	   * m_rkc[] is zero for irreversibly reactions
	   */
	  const vector_fp& m_rkc = m_kdata->m_rkcn;
	  for (int i = 0; i < m_ii; i++) {
              krev[i] *= m_rkc[i];
	  }
	}
    }

{'fraction_non_alphanumeric': 0.12554493307839387, 'fraction_numerical': 0.005430210325047801, 'mean_word_length': 2.409517601043025, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 102, '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': '6857576', 'n_tokens_mistral': 9269, 'n_tokens_neox': 8611, 'n_words': 2171}