<commit_before>/**
 * @file MultiPhaseEquil.cpp
 */
#include "MultiPhaseEquil.h"
#include "MultiPhase.h"
#ifdef WITH_ELECTROLYTES
#include "MolalityVPSSTP.h"
#endif
#include "sort.h"
#include "global.h"

#include <math.h>
#include <iostream>
using namespace std;


namespace Cantera {

    const doublereal TINY = 1.0e-20;

#if defined(WITH_HTML_LOGS)
    /// Used to print reaction equations. Given a stoichiometric
    /// coefficient 'nu' and a chemical symbol 'sym', return a string
    /// for this species in the reaction. 
    /// @param first  if this is false, then a " + " string will be
    /// added to the beginning of the string.
    /// @param nu  Stoichiometric coefficient. May be positive or negative. The
    /// absolute value will be used in the string. 
    /// @param sym Species chemical symbol.
    /// 
    static string coeffString(bool first, doublereal nu, string sym) {
        if (nu == 0.0) return "";
        string strt = " + ";
        if (first) strt = "";
        if (nu == 1.0 || nu == -1.0) 
            return strt + sym;
        string s = fp2str(fabs(nu));
        return strt + s + " " + sym;
    }
#endif

    /// Constructor. Construct a multiphase equilibrium manager for a
    /// multiphase mixture.
    /// @param mix Pointer to a multiphase mixture object.
    /// @param start If true, the initial composition will be
    /// determined by a linear Gibbs minimization, otherwise the
    /// initial mixture composition will be used.
    MultiPhaseEquil::MultiPhaseEquil(mix_t* mix, bool start, int loglevel) : m_mix(mix)
    {
        // the multi-phase mixture
        //        m_mix = mix;

        // store some mixture parameters locally
        m_nel_mix = mix->nElements();
        m_nsp_mix = mix->nSpecies();
        m_np = mix->nPhases();
        m_press = mix->pressure();
        m_temp = mix->temperature();

        index_t m, k;
        m_nel = 0;
        m_nsp = 0;
        m_eloc = 1000;
        m_incl_species.resize(m_nsp_mix,1);
        m_incl_element.resize(m_nel_mix,1);
        for (m = 0; m < m_nel_mix; m++) {
            string enm = mix->elementName(m);
            // element 'E' or 'e' represents an electron; this
            // requires special handling, so save its index
            // for later use
            if (enm == "E" || enm == "e") {
                m_eloc = m;
            }
            // if an element other than electrons is not present in
            // the mixture, then exclude it and all species containing
            // it from the calculation. Electrons are a special case,
            // since a species can have a negative number of 'atoms'
            // of electrons (positive ions).
            if (m_mix->elementMoles(m) <= 0.0) {
                if (m != m_eloc) {
                    m_incl_element[m] = 0;
                    for (k = 0; k < m_nsp_mix; k++) {
                        if (m_mix->nAtoms(k,m) != 0.0) {
                            m_incl_species[k] = 0;
                        }
                    }
                }
            }
        }
        
        // Now build the list of elements to be included, starting with
        // electrons, if they are present.
        if (m_eloc < m_nel_mix) {
            m_element.push_back(m_eloc);
            m_nel++;
        }
        // add the included elements other than electrons
        for (m = 0; m < m_nel_mix; m++) {
            if (m_incl_element[m] == 1 && m != m_eloc) {
                m_nel++;
                m_element.push_back(m);
            }
        }

        // include pure single-constituent phases only if their thermo
        // data are valid for this temperature. This is necessary,
        // since some thermo polynomial fits are done only for a
        // limited temperature range. For example, using the NASA
        // polynomial fits for solid ice and liquid water, if this
        // were not done the calculation would predict solid ice to be
        // present far above its melting point, since the thermo
        // polynomial fits only extend to 273.15 K, and give
        // unphysical results above this temperature, leading
        // (incorrectly) to Gibbs free energies at high temperature
        // lower than for liquid water.
        index_t ip;
        for (k = 0; k < m_nsp_mix; k++) {
            ip = m_mix->speciesPhaseIndex(k);
            if (!m_mix->solutionSpecies(k) && 
                !m_mix->tempOK(ip)) {
                m_incl_species[k] = 0;
                if (m_mix->speciesMoles(k) > 0.0) {
                    throw CanteraError("MultiPhaseEquil",
                        "condensed-phase species"+ m_mix->speciesName(k) 
                        + " is excluded since its thermo properties are \n"
                        "not valid at this temperature, but it has "
                        "non-zero moles in the initial state.");
                }
            }
        }

        // Now build the list of all species to be included in the
        // calculation.
        for (k = 0; k < m_nsp_mix; k++) {
            if (m_incl_species[k] ==1) {
                m_nsp++;
                m_species.push_back(k);
            }
        }

        // some work arrays for internal use
        m_work.resize(m_nsp);
        m_work2.resize(m_nsp);
        m_work3.resize(m_nsp_mix);
        m_mu.resize(m_nsp_mix);

        // number of moles of each species
        m_moles.resize(m_nsp);
        m_lastmoles.resize(m_nsp);
        m_dxi.resize(m_nsp - m_nel);

        // initialize the mole numbers to the mixture composition
        index_t ik;
        for (ik = 0; ik < m_nsp; ik++) {
            m_moles[ik] = m_mix->speciesMoles(m_species[ik]);
        }

        // Delta G / RT for each reaction
        m_deltaG_RT.resize(m_nsp - m_nel, 0.0);

        m_majorsp.resize(m_nsp);
        m_sortindex.resize(m_nsp,0);
        m_lastsort.resize(m_nel);
        m_solnrxn.resize(m_nsp - m_nel);
        m_A.resize(m_nel, m_nsp, 0.0);
        m_N.resize(m_nsp, m_nsp - m_nel);
        m_order.resize(m_nsp, 0);
	for (k = 0; k < m_nsp; k++) {
	  m_order[k] = k;
	}

        // if the 'start' flag is set, estimate the initial mole
        // numbers by doing a linear Gibbs minimization. In this case,
        // only the elemental composition of the initial mixture state
        // matters.
        if (start) {
            setInitialMoles(loglevel-1);
        }
        computeN();

        // Take a very small step in composition space, so that no
        // species has precisely zero moles.
        vector_fp dxi(m_nsp - m_nel, 1.0e-20);
        multiply(m_N, DATA_PTR(dxi), DATA_PTR(m_work));
        unsort(m_work);

        for (k = 0; k < m_nsp; k++) {
            m_moles[k] += m_work[k];<commit_msg>Fixed illegal access to empty vector in MultiPhaseEquil constructor<commit_after>/**
 * @file MultiPhaseEquil.cpp
 */
#include "MultiPhaseEquil.h"
#include "MultiPhase.h"
#ifdef WITH_ELECTROLYTES
#include "MolalityVPSSTP.h"
#endif
#include "sort.h"
#include "global.h"

#include <math.h>
#include <iostream>
using namespace std;


namespace Cantera {

    const doublereal TINY = 1.0e-20;

#if defined(WITH_HTML_LOGS)
    /// Used to print reaction equations. Given a stoichiometric
    /// coefficient 'nu' and a chemical symbol 'sym', return a string
    /// for this species in the reaction. 
    /// @param first  if this is false, then a " + " string will be
    /// added to the beginning of the string.
    /// @param nu  Stoichiometric coefficient. May be positive or negative. The
    /// absolute value will be used in the string. 
    /// @param sym Species chemical symbol.
    /// 
    static string coeffString(bool first, doublereal nu, string sym) {
        if (nu == 0.0) return "";
        string strt = " + ";
        if (first) strt = "";
        if (nu == 1.0 || nu == -1.0) 
            return strt + sym;
        string s = fp2str(fabs(nu));
        return strt + s + " " + sym;
    }
#endif

    /// Constructor. Construct a multiphase equilibrium manager for a
    /// multiphase mixture.
    /// @param mix Pointer to a multiphase mixture object.
    /// @param start If true, the initial composition will be
    /// determined by a linear Gibbs minimization, otherwise the
    /// initial mixture composition will be used.
    MultiPhaseEquil::MultiPhaseEquil(mix_t* mix, bool start, int loglevel) : m_mix(mix)
    {
        // the multi-phase mixture
        //        m_mix = mix;

        // store some mixture parameters locally
        m_nel_mix = mix->nElements();
        m_nsp_mix = mix->nSpecies();
        m_np = mix->nPhases();
        m_press = mix->pressure();
        m_temp = mix->temperature();

        index_t m, k;
        m_nel = 0;
        m_nsp = 0;
        m_eloc = 1000;
        m_incl_species.resize(m_nsp_mix,1);
        m_incl_element.resize(m_nel_mix,1);
        for (m = 0; m < m_nel_mix; m++) {
            string enm = mix->elementName(m);
            // element 'E' or 'e' represents an electron; this
            // requires special handling, so save its index
            // for later use
            if (enm == "E" || enm == "e") {
                m_eloc = m;
            }
            // if an element other than electrons is not present in
            // the mixture, then exclude it and all species containing
            // it from the calculation. Electrons are a special case,
            // since a species can have a negative number of 'atoms'
            // of electrons (positive ions).
            if (m_mix->elementMoles(m) <= 0.0) {
                if (m != m_eloc) {
                    m_incl_element[m] = 0;
                    for (k = 0; k < m_nsp_mix; k++) {
                        if (m_mix->nAtoms(k,m) != 0.0) {
                            m_incl_species[k] = 0;
                        }
                    }
                }
            }
        }
        
        // Now build the list of elements to be included, starting with
        // electrons, if they are present.
        if (m_eloc < m_nel_mix) {
            m_element.push_back(m_eloc);
            m_nel++;
        }
        // add the included elements other than electrons
        for (m = 0; m < m_nel_mix; m++) {
            if (m_incl_element[m] == 1 && m != m_eloc) {
                m_nel++;
                m_element.push_back(m);
            }
        }

        // include pure single-constituent phases only if their thermo
        // data are valid for this temperature. This is necessary,
        // since some thermo polynomial fits are done only for a
        // limited temperature range. For example, using the NASA
        // polynomial fits for solid ice and liquid water, if this
        // were not done the calculation would predict solid ice to be
        // present far above its melting point, since the thermo
        // polynomial fits only extend to 273.15 K, and give
        // unphysical results above this temperature, leading
        // (incorrectly) to Gibbs free energies at high temperature
        // lower than for liquid water.
        index_t ip;
        for (k = 0; k < m_nsp_mix; k++) {
            ip = m_mix->speciesPhaseIndex(k);
            if (!m_mix->solutionSpecies(k) && 
                !m_mix->tempOK(ip)) {
                m_incl_species[k] = 0;
                if (m_mix->speciesMoles(k) > 0.0) {
                    throw CanteraError("MultiPhaseEquil",
                        "condensed-phase species"+ m_mix->speciesName(k) 
                        + " is excluded since its thermo properties are \n"
                        "not valid at this temperature, but it has "
                        "non-zero moles in the initial state.");
                }
            }
        }

        // Now build the list of all species to be included in the
        // calculation.
        for (k = 0; k < m_nsp_mix; k++) {
            if (m_incl_species[k] ==1) {
                m_nsp++;
                m_species.push_back(k);
            }
        }

        // some work arrays for internal use
        m_work.resize(m_nsp);
        m_work2.resize(m_nsp);
        m_work3.resize(m_nsp_mix);
        m_mu.resize(m_nsp_mix);

        // number of moles of each species
        m_moles.resize(m_nsp);
        m_lastmoles.resize(m_nsp);
        m_dxi.resize(m_nsp - m_nel);

        // initialize the mole numbers to the mixture composition
        index_t ik;
        for (ik = 0; ik < m_nsp; ik++) {
            m_moles[ik] = m_mix->speciesMoles(m_species[ik]);
        }

        // Delta G / RT for each reaction
        m_deltaG_RT.resize(m_nsp - m_nel, 0.0);

        m_majorsp.resize(m_nsp);
        m_sortindex.resize(m_nsp,0);
        m_lastsort.resize(m_nel);
        m_solnrxn.resize(m_nsp - m_nel);
        m_A.resize(m_nel, m_nsp, 0.0);
        m_N.resize(m_nsp, m_nsp - m_nel);
        m_order.resize(m_nsp, 0);
	for (k = 0; k < m_nsp; k++) {
	  m_order[k] = k;
	}

        // if the 'start' flag is set, estimate the initial mole
        // numbers by doing a linear Gibbs minimization. In this case,
        // only the elemental composition of the initial mixture state
        // matters.
        if (start) {
            setInitialMoles(loglevel-1);
        }
        computeN();

        // Take a very small step in composition space, so that no
        // species has precisely zero moles.
        vector_fp dxi(m_nsp - m_nel, 1.0e-20);
        if (!dxi.empty()) {
          multiply(m_N, DATA_PTR(dxi), DATA_PTR(m_work));
          unsort(m_work);
        }

        for (k = 0; k < m_nsp; k++) {
            m_moles[k] += m_work[k];
{'fraction_non_alphanumeric': 0.10609845701689934, 'fraction_numerical': 0.00822924320352682, 'mean_word_length': 1.6563231850117095, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 37, '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': '5984819', 'n_tokens_mistral': 4303, 'n_tokens_neox': 4102, 'n_words': 1289}