diff --git a/include/cantera/Cantera.h b/include/cantera/Cantera.h new file mode 100644 index 000000000..e01b711aa --- /dev/null +++ b/include/cantera/Cantera.h @@ -0,0 +1,36 @@ +#ifndef CANTERA_H_INCL +#define CANTERA_H_INCL + +// Current 'Cantera.h' header + + +// definitions +#ifndef CANTERA_APP +#define CANTERA_APP +#endif + +namespace Cantera_CXX{ } + +using namespace Cantera_CXX; + +#include "base/ct_defs.h" + +// some useful functions +#include "base/global.h" + +// the CanteraError exception class +#include "base/ctexceptions.h" + +// The Cantera logger class +#include "base/logger.h" + +// Include the timer +#include "base/clockWC.h" + +// Include routines for reading and writing XML files +#include "base/xml.h" + +// Include string utility routines +#include "base/stringUtils.h" + +#endif diff --git a/include/cantera/equil/vcs_DoubleStarStar.h b/include/cantera/equil/vcs_DoubleStarStar.h new file mode 100644 index 000000000..c72376e8b --- /dev/null +++ b/include/cantera/equil/vcs_DoubleStarStar.h @@ -0,0 +1,129 @@ +/** + * @file vcs_DoubleStarStar.h + * + * Header file for class DoubleStarStar + */ +#ifndef VCS_DOUBLESTARSTAR_H +#define VCS_DOUBLESTARSTAR_H + +#include + +namespace VCSnonideal +{ + +using std::size_t; + +//! A class for 2D double arrays stored in column-major +//! (Fortran-compatible) form. +/*! + * In this form, the data entry for an n row, m col + * matrix is + * index = i + (n-1) * j + * where + * Matrix[j][i] + * i = row + * j = column + * The way this is instantiated is via the constructor: + * DoubleStarStar Dmatrix(mcol, mrow); + * + * The way this is referenced is via the notation: + * Dmatrix[icol][irow] + */ +class DoubleStarStar +{ + +public: + + //! Default constructor. Create an empty array. + DoubleStarStar(); + + //! Constructor. + /*! + * Create an \c nrow by \c mcol double array, and initialize + * all elements to \c v. + * + * @param mcol Number of columns + * @param nrow Number of rows + */ + DoubleStarStar(size_t mcol, size_t nrow, double v = 0.0); + + //! copy constructor + /*! + * @param y object to be copied + */ + DoubleStarStar(const DoubleStarStar& y); + + /// assignment operator + /*! + * @param y object to be copied + */ + DoubleStarStar& operator=(const DoubleStarStar& y); + + //! Resize the array, and fill the new entries with 'v' + /*! + * @param mrow This is the number of columns in the new matrix + * @param ncol This is the number of rows + * @param v Default fill value -> defaults to zero. + */ + void resize(size_t mcol, size_t nrow, double v = 0.0); + + //! Pointer to the top of the column + /*! + * @param jcol This is the jth column + * + * @return returns the pointer to the top of the jth column + */ + double* operator[](size_t jcol); + + //! Returns a const Pointer to the top of the jth column + /*! + * @param jcol This is the jth column + * + * @return returns the pointer to the top of the jth column + */ + const double* operator[](size_t jcol) const; + + //! Returns a double ** pointer to the base address + /*! + * This is the second way to get to the data + * This returns a double ** which can later be used in + * Dmatrix[icol][irow] notation to get to the data + */ + double* const* baseDataAddr(); + + //! Returns a const double ** pointer to the base address + /*! + * This is the second way to get to the data + * This returns a double ** which can later be used in + * Dmatrix[icol][irow] notation to get to the data + */ + double const* const* constBaseDataAddr() const; + + //! Number of rows + size_t nRows() const; + + //! Number of columns + size_t nColumns() const; + +private: + //! Storage area + std::vector m_data; + + //! Vector of addresses for the top of the columns + /*! + * Length = mcol + */ + std::vector m_colAddr; + + //! number of rows + size_t m_nrows; + + //! number of columns + size_t m_ncols; +}; + +} + +#endif + + diff --git a/include/cantera/equil/vcs_IntStarStar.h b/include/cantera/equil/vcs_IntStarStar.h new file mode 100644 index 000000000..5e9474f11 --- /dev/null +++ b/include/cantera/equil/vcs_IntStarStar.h @@ -0,0 +1,113 @@ +/** + * @file IntStarStar.h + * + * Header file for class IntStarStar + */ +#ifndef VCS_INTSTARSTAR_H +#define VCS_INTSTARSTAR_H + +#include + +namespace VCSnonideal +{ +using std::size_t; + +//! A class for 2D int arrays stored in column-major +//! (Fortran-compatible) form. +/*! + * In this form, the data entry for an n row, m col + * matrix is + * index = i + (n-1) * j + * where + * Matrix[j][i] + * i = row + * j = column + */ +class IntStarStar +{ + +public: + + //! Default constructor. Create an empty array. + IntStarStar(); + + //! Constructor. + /*! + * Create an \c nrow by \c mcol int array, and initialize + * all elements to \c v. + * + * @param mcol Number of columns + * @param nrow Number of rows + */ + IntStarStar(size_t mcol, size_t nrow, int v = 0); + + //! Copy constructor + /*! + * @param y Object to be copied + */ + IntStarStar(const IntStarStar& y); + + //! Assignment operator + /*! + * @param y Object to be copied + */ + IntStarStar& operator=(const IntStarStar& y); + + //! Resize the array, and fill the new entries with 'v' + /*! + * @param mcol This is the number of columns in the new matrix + * @param nrow This is the number of rows + * @param v Default fill value -> defaults to zero. + */ + void resize(size_t mcol, size_t nrow, int v = 0); + + //! Pointer to the top of the column + /*! + * @param jcol Pointer to the top of the jth column + */ + int* operator[](size_t jcol); + + //! Pointer to the top of the column + /*! + * @param j Pointer to the top of the jth column + */ + const int* operator[](size_t jcol) const; + + //! Returns a int ** pointer to the base address + /*! + * This is the second way to get to the data + * This returns a int ** which can later be used in + * Imatrix[icol][irow] notation to get to the data + */ + int* const* baseDataAddr(); + + //! Number of rows + size_t nRows() const; + + //! Number of columns + size_t nColumns() const; + +private: + //! Storage area for the matrix, layed out in Fortran style, row-inner, column outer format + /*! + * Length = m_nrows * m_ncols + */ + std::vector m_data; + + //! Vector of column addresses + /*! + * Length = number of columns = m_ncols + */ + std::vector m_colAddr; + + //! number of rows + size_t m_nrows; + + //! number of columns + size_t m_ncols; +}; + +} + +#endif + diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h new file mode 100644 index 000000000..972ed1da3 --- /dev/null +++ b/include/cantera/equil/vcs_VolPhase.h @@ -0,0 +1,1038 @@ +/** + * @file vcs_VolPhase.h + * Header for the object representing each phase within vcs + */ +/* + * Copyright (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef VCS_VOLPHASE_H +#define VCS_VOLPHASE_H + +#include "cantera/equil/vcs_DoubleStarStar.h" + +#include +#include + +/* + * Forward references + */ +// Forward reference for ThermoPhase object within the Cantera namespace +namespace Cantera +{ +class ThermoPhase; +} + +namespace VCSnonideal +{ +/* + * Models for the species activity coefficients + * + */ +#define VCS_AC_CONSTANT 0 +//#define VCS_AC_DEBYE_HUCKEL 23 +//#define VCS_AC_REGULAR_SOLN 25 +//#define VCS_AC_MARGULES 300 +#define VCS_AC_UNK_CANTERA -1 +#define VCS_AC_UNK -2 +/* + * + * Models for the standard state volume of each species + */ +#define VCS_SSVOL_IDEALGAS 0 +#define VCS_SSVOL_CONSTANT 1 + +/* + * DEFINITIONS FOR THE vcs_VolPhase structure + * + * + * Equation of State Types + * - Permissible values for the EqnState variable in CPC_PHASE structure + */ +#define VCS_EOS_CONSTANT 0 +#define VCS_EOS_IDEAL_GAS 1 +#define VCS_EOS_STOICH_SUB 5 +#define VCS_EOS_IDEAL_SOLN 22 +#define VCS_EOS_DEBEYE_HUCKEL 23 +#define VCS_EOS_REDLICK_KWONG 24 +#define VCS_EOS_REGULAR_SOLN 25 +#define VCS_EOS_UNK_CANTERA -1 + + +struct VCS_SPECIES; +class vcs_SpeciesProperties; +class VCS_SOLVE; + + +//! Phase information and Phase calculations for vcs. +/*! + * Each phase in a vcs calculation has a vcs_VolPhase object associated + * with it. This object helps to coordinate property evaluations for + * species within the phase. Usually these evaluations must be carried + * out on a per phase basis. However, vcs frequently needs per species + * quantitites. Therefore, we need an interface layer between vcs + * and Cantera's ThermoPhase. + * + * The species stay in the same ordering within this structure. + * The vcs algorithm will change the ordering of species in + * the global species list. However, the indexing of species in this + * list stays the same. This structure contains structures that + * point to the species belonging to this phase in the global + * vcs species list. + * + * This object is considered not to own the underlying Cantera ThermoPhase + * object for the phase. + * + * This object contains an idea of the temperature and pressure. + * It checks to see if if the temperature and pressure has changed before calling + * underlying property evalulation routines. + * + * The object contains values for the electric potential of a phase. + * It coordinates the evalulation of properties wrt when the electric + * potential of a phase has changed. + * + * The object knows about the mole fractions of the phase. It controls + * the values of mole fractions, and coordinates the property evalulation + * wrt to changes in the mole fractions. It also will keep track of the + * likely values of mole fractions in multicomponent phases even when + * the phase doesn't actually exist within the thermo program. + * + * The object knows about the total moles of a phase. It checkes to + * see if the phase currently exists or not, and modifies its behavior + * accordingly. + * + * + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + */ +class vcs_VolPhase +{ +public: + + /************************************************************************* + * FUNCTIONS * + ************************************************************************/ + + //! Base constructor for the class + vcs_VolPhase(VCS_SOLVE* owningSolverObject = 0); + + //! Copy constructor + /*! + * @param b object to be copied + */ + vcs_VolPhase(const vcs_VolPhase& b); + + //! Assignment operator + /*! + * @param b object to be copied + */ + vcs_VolPhase& operator=(const vcs_VolPhase& b); + + //! Destructor + ~vcs_VolPhase(); + + + //! The resize() function fills in all of the initial information if it + //! is not given in the constructor. + /*! + * @param phaseNum index of the phase in the vcs problem + * @param numSpecies Number of species in the phase + * @param phaseName String name for the phase + * @param molesInert kmoles of inert in the phase (defaults to zero) + */ + void resize(const size_t phaseNum, const size_t numSpecies, + const size_t numElem, const char* const phaseName, + const double molesInert = 0.0); + + void elemResize(const size_t numElemConstraints); + + //! Evaluate activity coefficients and return the kspec coefficient + /*! + * We carry out a calculation whenever UpTODate_AC is false. Specifically + * whenever a phase goes zero, we do not carry out calculations on it. + * + * @param kspec species number + */ + double AC_calc_one(size_t kspec) const; + + + //! Set the moles and/or mole fractions within the phase + /*! + * Sets the mole fraction and total moles within the phase + * + * @param molNum total moles in the phase + * @param moleFracVec Vector of input mole fractions + * @param vcsStateStatus Status flag for this update + */ + void setMoleFractionsState(const double molNum, const double* const moleFracVec, + const int vcsStateStatus); + + //! Set the moles within the phase + /*! + * This function takes as input the mole numbers in vcs format, and + * then updates this object with their values. This is essentially + * a gather routine. + * + * @param molesSpeciesVCS Array of mole numbers. Note, the indices + * for species in + * this array may not be contiguous. IndSpecies[] is needed + * to gather the species into the local contiguous vector + * format. + */ + void setMolesFromVCS(const int stateCalc, + const double* molesSpeciesVCS = 0); + + //! Set the moles within the phase + /*! + * This function takes as input the mole numbers in vcs format, and + * then updates this object with their values. This is essentially + * a gather routine. + * Additionally it checks to see that the total moles value in + * TPhMoles[iplace] is equal to the internally computed value. + * If this isn't the case, an error exit is carried out. + * + * + * @param molesSpeciesVCS array of mole numbers. Note, the indices + * for species in + * this array may not be contiguous. IndSpecies[] is needed + * to gather the species into the local contiguous vector + * format. + * @param TPhMoles VCS's array containing the number of moles + * in each phase. + * @param iphase index of the current phase. + * + */ + void setMolesFromVCSCheck(const int stateCalc, + const double* molesSpeciesVCS, + const double* const TPhMoles); + + //! Update the moles within the phase, if necessary + /*! + * This function takes as input the stateCalc value, which + * determines where within VCS_SOLVE to fetch the mole numbers. + * It then updates this object with their values. This is essentially + * a gather routine. + * + * @param stateCalc State calc value either VCS_STATECALC_OLD + * or VCS_STATECALC_NEW. With any other value + * nothing is done. + * + */ + void updateFromVCS_MoleNumbers(const int stateCalc); + + //! Fill in an activity coefficients vector within a VCS_SOLVE object + /*! + * This routine will calculate the activity coefficients for the + * current phase, and fill in the corresponding entries in the + * VCS activity coefficients vector. + * + * @param AC vector of activity coefficients for all of the species + * in all of the phases in a VCS problem. Only the + * entries for the current phase are filled in. + */ + void sendToVCS_ActCoeff(const int stateCalc, double* const AC); + + //! set the electric potential of the phase + /*! + * @param phi electric potential (volts) + */ + void setElectricPotential(const double phi); + + //! Returns the electric field of the phase + /*! + * Units are potential + */ + double electricPotential() const; + + //! Gibbs free energy calculation for standard state of one species + /*! + * Calculate the Gibbs free energies for the standard state + * of the kth species. + * The results are held internally within the object. + * The kth species standard state G is returned + * + * @param kspec Species number (within the phase) + * @param TKelvin Current temperature + * @param pres Current pressure + * + * @return Gstar[kspec] returns the gibbs free energy for the + * standard state of the kth species. + */ + double GStar_calc_one(size_t kspec) const; + + //! Gibbs free energy calculation at a temperature for the reference state + //! of a species, return a value for one species + /*! + * @param kspec species index + * @param TKelvin temperature + * + * @return return value of the gibbs free energy + */ + double G0_calc_one(size_t kspec) const; + + //! Molar volume calculation for standard state of one species + /*! + * Calculate the molar volume for the standard states + * The results are held internally within the object. + * Return the molar volume for one species + * + * @param kspec Species number (within the phase) + * @param TKelvin Current temperature + * @param pres Current pressure + * + * @return molar volume of the kspec species's standard + * state (m**3/kmol) + */ + double VolStar_calc_one(size_t kglob) const; + + //! Fill in the partial molar volume vector for VCS + /*! + * This routine will calculate the partial molar volumes for the + * current phase (if needed), and fill in the corresponding entries in the + * VCS partial molar volumes vector. + * + * @param VolPM vector of partial molar volumes for all of the species + * in all of the phases in a VCS problem. Only the + * entries for the current phase are filled in. + */ + double sendToVCS_VolPM(double* const VolPM) const; + + //! Fill in the partial molar volume vector for VCS + /*! + * This routine will calculate the partial molar volumes for the + * current phase (if needed), and fill in the corresponding entries in the + * VCS partial molar volumes vector. + * + * @param VolPM vector of partial molar volumes for all of the species + * in all of the phases in a VCS problem. Only the + * entries for the current phase are filled in. + */ + void sendToVCS_GStar(double* const gstar) const; + + //! Sets the temperature and pressure in this object and + //! underlying objects + /*! + * Sets the temperature and pressure in this object and + * underlying objects. The underlying objects refers to the + * Cantera's ThermoPhase object for this phase. + * + * @param temperature_Kelvin (Kelvin) + * @param pressure_PA Pressure (MKS units - Pascal) + */ + void setState_TP(const double temperature_Kelvin, const double pressure_PA); + + //! Sets the temperature in this object and + //! underlying objects + /*! + * Sets the temperature and pressure in this object and + * underlying objects. The underlying objects refers to the + * Cantera's ThermoPhase object for this phase. + * + * @param temperature_Kelvin (Kelvin) + */ + void setState_T(const double temperature_Kelvin); + + // Downloads the ln ActCoeff jacobian into the VCS version of the + // ln ActCoeff jacobian. + /* + * + * This is essentially a scatter operation. + * + * @param LnAcJac_VCS jacobian parameter + * The Jacobians are actually d( lnActCoeff) / d (MolNumber); + * dLnActCoeffdMolNumber[j][k] + * + * j = id of the species mole number + * k = id of the species activity coefficient + */ + void sendToVCS_LnActCoeffJac(double* const* const LnACJac_VCS); + + //! Set the pointer for Cantera's ThermoPhase parameter + /*! + * When we first initialize the ThermoPhase object, we read the + * state of the ThermoPhase into vcs_VolPhase object. + * + * @param tp_ptr Pointer to the ThermoPhase object corresponding + * to this phase. + */ + void setPtrThermoPhase(Cantera::ThermoPhase* tp_ptr); + + //! Return a const ThermoPhase pointer corresponding to this phase + /*! + * @return pointer to the ThermoPhase. + */ + const Cantera::ThermoPhase* ptrThermoPhase() const; + + //! Return the total moles in the phase + /*! + * + * Units -> depends on VCS_UnitsFormat variable + * Cantera -> J/kmol + */ + double totalMoles() const; + + //! Returns the mole fraction of the kspec species + /*! + * @param kspec Index of the species in the phase + * + * @return Value of the mole fraction + */ + double molefraction(size_t kspec) const; + + //! Sets the total moles in the phase + /*! + * We don't have to flag the internal state as changing here + * because we have just changed the total moles. + * + * @param totalMols Total moles in the phase (kmol) + */ + void setTotalMoles(const double totalMols); + + //! Sets the mole flag within the object to out of date + /*! + * This will trigger the object to go get the current mole numbers + * when it needs it. + */ + void setMolesOutOfDate(int stateCalc = -1); + + //! Sets the mole flag within the object to be current + /*! + * + */ + void setMolesCurrent(int vcsStateStatus); + +private: + //! Set the mole fractions from a conventional mole fraction vector + /*! + * + * @param xmol Value of the mole fractions for the species + * in the phase. These are contiguous. + */ + void setMoleFractions(const double* const xmol); + +public: + + //! Return a const reference to the mole fractions stored in the + //! object. + const std::vector & moleFractions() const; + + double moleFraction(size_t klocal) const; + + //! Sets the creationMoleNum's within the phase object + /*! + * @param F_k Pointer to a vector of n_k's + */ + void setCreationMoleNumbers(const double* const n_k, const std::vector &creationGlobalRxnNumbers); + + //! Return a const reference to the creationMoleNumbers stored in the object. + /*! + * @return Returns a const reference to the vector of creationMoleNumbers + */ + const std::vector & creationMoleNumbers(std::vector &creationGlobalRxnNumbers) const; + + //! Returns whether the phase is an ideal solution phase + bool isIdealSoln() const; + + //! Returns whether the object is using cantera calls. + bool usingCanteraCalls() const; + + //! Return the index of the species that represents the + //! the voltage of the phase + size_t phiVarIndex() const; + + void setPhiVarIndex(size_t phiVarIndex); + + //! Retrieve the kth Species structure for the species belonging to this phase + /*! + * The index into this vector is the species index within the phase. + * + * @param kindex kth species index. + */ + vcs_SpeciesProperties* speciesProperty(const size_t kindex); + + //! int indicating whether the phase exists or not + /*! + * returns the m_existence int for the phase + * + * - VCS_PHASE_EXIST_ZEROEDPHASE = -6: Set to not exist by fiat from a + * higher level. + * This is used in phase stability boundary calculations + * - VCS_PHASE_EXIST_NO = 0: Doesn't exist currently + * - VCS_PHASE_EXIST_MINORCONC = 1: Exists, but the concentration is + * so low that an alternate + * method is used to calculate the total phase concentrations. + * - VCS_PHASE_EXIST_YES = 2 : Does exist currently + * - VCS_PHASE_EXIST_ALWAYS = 3: Always exists because it contains + * inerts which can't exist in any other phase. Or, + * the phase exists always because it consists of a single + * species, which is identified with the voltage, i.e., + * it's an electron metal phase. + */ + int exists() const; + + //! Set the existence flag in the object + /*! + * Note the total moles of the phase must have been set appropriately + * before calling this routine. + * + * @param existence Phase existence flag + * + * @note try to eliminate this routine + */ + void setExistence(const int existence); + + //! Return the Global VCS index of the kth species in the phase + /*! + * @param spIndex local species index (0 to the number of species + * in the phase) + * + * @return Returns the VCS_SOLVE species index of the species. + * This changes as rearrangements are carried out. + */ + size_t spGlobalIndexVCS(const size_t spIndex) const; + + + //! set the Global VCS index of the kth species in the phase + /*! + * @param spIndex local species index (0 to the number of species + * in the phase) + * + * @return Returns the VCS_SOLVE species index of the that species + * This changes as rearrangements are carried out. + */ + void setSpGlobalIndexVCS(const size_t spIndex, const size_t spGlobalIndex); + + //! Sets the total moles of inert in the phase + /*! + * @param tMolesInert Value of the total kmols of inert species in the + * phase. + */ + void setTotalMolesInert(const double tMolesInert); + + //! returns the value of the total kmol of inert in the phase + /*! + * @return Returns the total value of the kmol of inert in the phase + */ + double totalMolesInert() const; + + //! Returns the global index of the local element index for the phase + size_t elemGlobalIndex(const size_t e) const; + + //! sets a local phase element to a global index value + /*! + * @param eLocal Local phase element index + * @param eGlobal Global phase element index + */ + void setElemGlobalIndex(const size_t eLocal, const size_t eGlobal); + + //! Returns the number of element constraints + size_t nElemConstraints() const; + + //! Name of the element constraint with index \c e. + /*! + * @param e Element index. + */ + std::string elementName(const size_t e) const; + + //! Type of the element constraint with index \c e. + /*! + * @param e Element index. + */ + int elementType(const size_t e) const; + + //! Set the element Type of the element constraint with index \c e. + /*! + * @param e Element index + * @param eType type of the element. + */ + void setElementType(const size_t e, const int eType); + + //! Transfer all of the element information from the + //! ThermoPhase object to the vcs_VolPhase object. + /*! + * Also decide whether we need a new charge neutrality + * element in the phase to enforce a charge neutrality + * constraint. + * + * @param tPhase Pointer to the thermophase object + */ + size_t transferElementsFM(const Cantera::ThermoPhase* const tPhase); + + //! Get a constant form of the Species Formula Matrix + /*! + * Returns a double ** pointer such that + * + * fm[e][f] is the formula matrix entry for element e for species k + */ + double const* const* getFormulaMatrix() const; + + //! Returns the type of the species unknown + /*! + * @param k species index + * + * returns the SpeciesUnknownType[k] = type of species + * Normal -> VCS_SPECIES_TYPE_MOLUNK + * ( unknown is the mole number in the phase) + * metal electron -> VCS_SPECIES_INTERFACIALVOLTAGE + * ( unknown is the interfacial voltage (volts) + */ + int speciesUnknownType(const size_t k) const; + + + int elementActive(const size_t e) const; + + + //! Return the number of species in the phase + size_t nSpecies() const; + +private: + + //! Evaluate the activity coefficients at the current conditions + /*! + * We carry out a calculation whenever UpTODate_AC is false. Specifically + * whenever a phase goes zero, we do not carry out calculations on it. + */ + void _updateActCoeff() const; + + //! Gibbs free energy calculation for standard states + /*! + * Calculate the Gibbs free energies for the standard states + * The results are held internally within the object. + * + * @param TKelvin Current temperature + * @param pres Current pressure + */ + void _updateGStar() const; + + //! Gibbs free energy calculation at a temperature for the reference state + //! of each species + /*! + * + */ + void _updateG0() const; + + //! Molar volume calculation for standard states + /*! + * Calculate the molar volume for the standard states + * The results are held internally within the object. + * + * @param TKelvin Current temperature + * @param pres Current pressure + * + * Units are in m**3/kmol + */ + void _updateVolStar() const; + + //! Calculate the partial molar volumes of all species and return the + //! total volume + /*! + * Calculates these quantitites internally + * + * @return total volume + */ + double _updateVolPM() const; + + //! Evaluation of Activity Coefficient Jacobians + /*! + * This is the derivative of the ln of the activity coefficient + * with respect to mole number of jth species. + * (temp, pressure, and other mole numbers held constant) + * + * We employ a finite difference derivative approach here. + * Because we have to change the mole numbers, this is not + * a const function, even though the paradigm would say that + * it should be. + * + * @param moleNumbers Mole numbers are input. + */ + void _updateLnActCoeffJac(); + + //! Updates the mole fraction depenpencies + /*! + * Whenever the mole fractions change, this routine + * should be called. + */ + void _updateMoleFractionDependencies(); + + + /************************************************************************* + * MEMBER DATA * + ************************************************************************/ + +private: + //! Backtrack value of VCS_SOLVE * + /*! + * Note the default for this is 0. That's a valid value too, since + * VCS_PROB also uses vcs_VolPhase objects. + */ + VCS_SOLVE* m_owningSolverObject; + +public: + //! Original ID of the phase in the problem. + /*! + * If a non-ideal phase splits into two due to a + * miscibility gap, these numbers will stay the + * same after the split. + */ + size_t VP_ID_; + + //! ID of the surface or volume domain in which the + //! this phase exists + /*! + * This ventures into the idea of installing a physical location + * into a thermodynamics program. This unknown is currently not + * being used. + * @deprecated + */ + int Domain_ID; + + //! If true, this phase consists of a single species + bool m_singleSpecies; + + //! If true, this phase is a gas-phase like phase + /*! + * A RTlog(p/1atm) term is added onto the chemical potential for inert + * species if this is true. + */ + bool m_gasPhase; + + //! Type of the equation of state + /*! + * The known types are listed at the top of this file. + */ + int m_eqnState; + + //! This is the element number for the charge neutrality + //! condition of the phase + /*! + * If it has one. If it does not have a charge neutrality + * constraint, then this value is equal to -1 + */ + size_t ChargeNeutralityElement; + + //! Units for the chemical potential data, pressure data, volume, + //! and species amounts + /*! + * All internally stored quantities will have these units. Also, printed + * quantitities will display in these units. Input quantities are expected + * in these units. + * + * Chem_Pot Pres vol moles + * ---------------------------------------------------------------------- + * -1 VCS_UNITS_KCALMOL = kcal/gmol Pa m**3 kmol + * 0 VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol + * 1 VCS_UNITS_KJMOL = kJ / gmol Pa m**3 kmol + * 2 VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol + * 3 VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol + * ---------------------------------------------------------------------- + * + * see vcs_defs.h for more information. + * + * Currently, this value should be the same as the owning VCS_PROB or + * VCS_SOLVE object. There is no code for handling anything else atm. + * + * (This variable is needed for the vcsc code, where it is not equal + * to VCS_UNITS_MKS). + */ + int p_VCS_UnitsFormat; + + //! Convention for the activity formulation + /*! + * 0 = molar based activities (default) + * 1 = Molality based activities + * mu = mu_0 + ln a_molality + * standard state is based on unity molality + */ + int p_activityConvention; + +private: + //! Number of element constraints within the problem + /*! + * This is usually equal to the number of elements. + */ + size_t m_numElemConstraints; + + //! vector of strings containing the element constraint names + /*! + * Length = nElemConstraints + */ + std::vector m_elementNames; + + //! boolean indicating whether an element constraint is active + //! for the current problem + std::vector m_elementActive; + + //! Type of the element constraint + /*! + * m_elType[j] = type of the element + * 0 VCS_ELEM_TYPE_ABSPOS Normal element that is positive + * or zero in all species. + * 1 VCS_ELEM_TPYE_ELECTRONCHARGE element dof that corresponds + * to the charge DOF. + * 2 VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may + * mean that a species has neg 0 or pos value + * of that constraint (other than charge) + */ + std::vector m_elementType; + + //! Formula Matrix for the phase + /*! + * FormulaMatrix[j][kspec] + * = Formula Matrix for the species + * Number of elements, j, + * in the kspec species + */ + DoubleStarStar m_formulaMatrix; + + //! Type of the species unknown + /*! + * SpeciesUnknownType[k] = type of species + * Normal -> VCS_SPECIES_TYPE_MOLUNK + * ( unknown is the mole number in the phase) + * metal electron -> VCS_SPECIES_INTERFACIALVOLTAGE + * ( unknown is the interfacial voltage (volts) + */ + std::vector m_speciesUnknownType; + + //! Index of the element number in the global list of elements + //! stored in VCS_PROB or VCS_SOLVE + std::vector m_elemGlobalIndex; + + //! Number of species in the phase + size_t m_numSpecies; + +public: + //! String name for the phase + std::string PhaseName; + +private: + //! Total moles of inert in the phase + double m_totalMolesInert; + + //! Boolean indicating whether the phase is an ideal solution + //! and therefore its molar-based activity coefficients are + //! uniformly equal to one. + bool m_isIdealSoln; + + //! Current state of existence: + /*! + * VCS_PHASE_EXIST_ZEROEDPHASE = -6: Set to not exist by fiat from a + * higher level. + * This is used in phase stability boundary calculations + * VCS_PHASE_EXIST_NO = 0: Doesn't exist currently + * VCS_PHASE_EXIST_MINORCONC = 1: Exists, but the concentration is + * so low that an alternate + * method is used to calculate the total phase concentrations. + * VCS_PHASE_EXIST_YES = 2 : Does exist currently + * VCS_PHASE_EXIST_ALWAYS = 3: Always exists because it contains + * inerts which can't exist in any other phase. Or, + * the phase exists always because it consists of a single + * species, which is identified with the voltage, i.e., + * its an electron metal phase. + */ + int m_existence; + + // Index of the first MF species in the list of unknowns for this phase + /*! + * This is always equal to zero. + * Am anticipating the case where the phase potential is species # 0, + * for multiphase phases. Right now we have the phase potential equal + * to 0 for single species phases, where we set by hand the mole fraction + * of species 0 to one. + */ + int m_MFStartIndex; + + //! Index into the species vectors + /*! + * Maps the phase species number into the global species number. + * Note, as part of the vcs algorithm, the order of the species + * vector is changed during the algorithm + */ + std::vector IndSpecies; + + //! Vector of Species structures for the species belonging to this phase + /*! + * The index into this vector is the species index within the phase. + */ + std::vector ListSpeciesPtr; + + //! If this is true, then calculations are actually performed within + //! Cantera + bool m_useCanteraCalls; + /** + * If we are using Cantera, this is the + * pointer to the ThermoPhase object. If not, this is null. + */ + Cantera::ThermoPhase* TP_ptr; + + //! Total mols in the phase + /*! + * units are kmol + */ + double v_totalMoles; + + //! Vector of the current mole fractions for species + //! in the phase + std::vector Xmol_; + + //! Vector of current creationMoleNumbers_ + /*! + * These are the actual unknowns in the phase stability problem + */ + std::vector creationMoleNumbers_; + + //! Vector of creation global reaction numbers for the phase stability problem + /*! + * The phase stability problem requires a global reaction number for each + * species in the phase. Usually this is the krxn = kglob - M for species + * in the phase that are not components. For component species, the + * choice of the reaction is one which maximimes the chance that the phase + * pops into (or remains in) existence. + * The index here is the local phase species index. + * the value of the variable is the global vcs reaction number. Note, + * that the global reaction number will go out of order when the species positions + * are swapped. So, this number has to be recalculated. + */ + std::vector creationGlobalRxnNumbers_; + + //! If the potential is a solution variable in VCS, it acts as a species. + //! This is the species index in the phase for the potential + size_t m_phiVarIndex; + + //! Total Volume of the phase + /*! + * units are m**3 + */ + mutable double m_totalVol; + + //! Vector of calculated SS0 chemical potentials for the + //! current Temperature. + /*! + * Note, This is the chemical potential derived strictly from the polynomial + * in temperature. Pressure effects have to be added in to + * get to the standard state. + * + * Units -> depends on VCS_UnitsFormat variable + * Cantera -> J/kmol + */ + mutable std::vector SS0ChemicalPotential; + + //! Vector of calculated Star chemical potentials for the + //! current Temperature and pressure. + /*! + * Note, This is the chemical potential at unit activity. Thus, we can call + * it the standard state chemical potential as well. + * + * Units -> depends on VCS_UnitsFormat variable + * Cantera -> J/kmol + */ + mutable std::vector StarChemicalPotential; + + //! Vector of the Star molar Volumes of the species. + /*! + * units m3 / kmol + */ + mutable std::vector StarMolarVol; + + //! Vector of the Partial molar Volumes of the species. + /*! + * units m3 / kmol + */ + mutable std::vector PartialMolarVol; + + //! Vector of calculated activity coefficients for the current state + /*! + * Whether or not this vector is current is determined by + * the bool m_UpToDate_AC. + */ + mutable std::vector ActCoeff; + + //! Vector of the derivatives of the ln activity coefficient wrt to the + //! current mole number + /*! + * dLnActCoeffdMolNumber[j][k]; + * j = id of the species mole number + * k = id of the species activity coefficient + */ + mutable DoubleStarStar dLnActCoeffdMolNumber; + + //! Status + /*! + * valid values are + * VCS_STATECALC_OLD + * VCS_STATECALC_NEW + * VCS_STATECALC_TMP + */ + int m_vcsStateStatus; + + + //! Value of the potential for the phase (Volts) + double m_phi; + + //! Boolean indicating whether the object has an uptodate mole number vector + //! and potential with respect to the current vcs state calc status + bool m_UpToDate; + + //! Boolean indicating whether activity coefficients are uptodate. + /*! + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + */ + mutable bool m_UpToDate_AC; + + //! Boolean indicating whether Star volumes are uptodate. + /*! + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + * Star volumes are sensitive to temperature and pressure + */ + mutable bool m_UpToDate_VolStar; + + //! Boolean indicating whether partial molar volumes are uptodate. + /*! + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + * partial molar volumes are sensitive to everything + */ + mutable bool m_UpToDate_VolPM; + + //! Boolean indicating whether GStar is uptodate. + /*! + * GStar is sensitive to the temperature and the pressure, only + */ + mutable bool m_UpToDate_GStar; + + + //! Boolean indicating whether G0 is uptodate. + /*! + * G0 is sensitive to the temperature and the pressure, only + */ + mutable bool m_UpToDate_G0; + + //! Current value of the temperature for this object, and underlying objects + double Temp_; + + //! Current value of the pressure for this object, and underlying objects + double Pres_; + + + +}; + +//! Return a string representing the equation of state +/*! + * @param EOSType : integer value of the equation of state + * + * @return returns a string representing the EOS + */ +std::string string16_EOSType(int EOSType); + +} + +#endif diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h new file mode 100644 index 000000000..c4df5bd61 --- /dev/null +++ b/include/cantera/equil/vcs_internal.h @@ -0,0 +1,530 @@ +/** + * @file vcs_internal.h + * Internal declarations for the VCSnonideal package + */ +/* + * Copyright (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef _VCS_INTERNAL_H +#define _VCS_INTERNAL_H + +#include + +#include "cantera/equil/vcs_defs.h" + +#include "cantera/base/global.h" + +#ifndef ALTLINPROG +#define ALTLINPROG +#endif + +namespace VCSnonideal +{ +using Cantera::npos; + +//! Points to the data in a std::vector<> object +#define VCS_DATA_PTR(vvv) (&(vvv[0])) + +//! define this Cantera function to replace printf +/*! + * We can replace this with printf easily + */ +#define plogf Cantera::writelogf + +//! define this Cantera function to replace cout << endl; +/*! + * We use this to place an endl in the log file, and + * ensure that the IO buffers are flushed. + */ +#define plogendl() Cantera::writelogendl() + +//! Global hook for turning on and off time printing. +/*! + * Default is to allow printing. But, you can assign this to zero + * globally to turn off all time printing. + * This is helpful for test suite purposes where you are interested + * in differences in text files. + */ +extern int vcs_timing_print_lvl; + +/* + * Forward references + */ +class VCS_SPECIES_THERMO; +class VCS_PROB; + +//! Amount of extra printing that is done while in debug mode. +/*! + * 0 -> none + * 1 -> some + * 2 -> alot (default) + * 3 -> everything + */ + +//! Class to keep track of time and iterations +/*! + * class keeps all of the counters together. + */ +class VCS_COUNTERS +{ +public: + //! Total number of iterations in the main loop + //! of vcs_TP() to solve for thermo equilibrium + int T_Its; + + //! Current number of iterations in the main loop + //! of vcs_TP() to solve for thermo equilibrium + int Its; + + //! Total number of optimizations of the + //! components basis set done + int T_Basis_Opts; + + //! number of optimizations of the components basis set done + int Basis_Opts; + + //! Current number of times the initial thermo + //! equilibrium estimator has been called + int T_Calls_Inest; + + //! Current number of calls to vcs_TP + int T_Calls_vcs_TP; + + //! Current time spent in vcs_TP + double T_Time_vcs_TP; + + //! Current time spent in vcs_TP + double Time_vcs_TP; + + //! Total Time spent in basopt + double T_Time_basopt; + + //! Current Time spent in basopt + double Time_basopt; + + //! Time spent in initial estimator + double T_Time_inest; + + //! Time spent in the vcs suite of programs + double T_Time_vcs; +}; + +//! Returns the value of the gas constant in the units specified by parameter +/*! + * @param mu_units Specifies the units. + * - VCS_UNITS_KCALMOL: kcal gmol-1 K-1 + * - VCS_UNITS_UNITLESS: 1.0 K-1 + * - VCS_UNITS_KJMOL: kJ gmol-1 K-1 + * - VCS_UNITS_KELVIN: 1.0 K-1 + * - VCS_UNITS_MKS: joules kmol-1 K-1 = kg m2 s-2 kmol-1 K-1 + */ +double vcsUtil_gasConstant(int mu_units); + +//! Invert an n x n matrix and solve m rhs's +/*! + * Solve a square matrix with multiple right hand sides + * + * \f[ + * C X + B = 0; + * \f] + * + * This routine uses Gauss elimination and is optimized for the solution + * of lots of rhs's. A crude form of row pivoting is used here. + * The matrix C is destroyed during the solve. + * + * @return The solution x[] is returned in the matrix B. + * Routine returns an integer representing success: + * - 1 : Matrix is singluar + * - 0 : solution is OK + * + * + * @param c Matrix to be inverted. c is in fortran format, i.e., rows + * are the inner loop. Row numbers equal to idem. + * c[i+j*idem] = c_i_j = Matrix to be inverted: + * - i = row number + * - j = column number + * + * @param idem number of row dimensions in c + * @param n Number of rows and columns in c + * @param b Multiple RHS. Note, b is actually the negative of + * most formulations. Row numbers equal to idem. + * b[i+j*idem] = b_i_j = vectors of rhs's: + * - i = row number + * - j = column number + * (each column is a new rhs) + * @param m number of rhs's + */ +int vcsUtil_mlequ(double* c, size_t idem, size_t n, double* b, size_t m); + +//! Invert an n x n matrix and solve m rhs's +/*! + * Solve a square matrix with multiple right hand sides + * + * \f[ + * C X + B = 0; + * \f] + * + * This routine uses Gauss-Jordan elimination and is optimized for the solution + * of lots of rhs's. Full row and column pivoting is used here. It's been + * shown to be necessary in at least one case. + * The matrix C is destroyed during the solve. + * + * @return The solution x[] is returned in the matrix B. + * Routine returns an integer representing success: + * - 1 : Matrix is singluar + * - 0 : solution is OK + * + * @param c Matrix to be inverted. c is in fortran format, i.e., rows + * are the inner loop. Row numbers equal to idem. + * c[i+j*idem] = c_i_j = Matrix to be inverted: + * - i = row number + * - j = column number + * + * @param idem number of row dimensions in c + * @param n Number of rows and columns in c + * @param b Multiple RHS. Note, b is actually the negative of + * most formulations. Row numbers equal to idem. + * b[i+j*idem] = b_i_j = vectors of rhs's: + * - i = row number + * - j = column number + * (each column is a new rhs) + * @param m number of rhs's + */ +int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m); + + +//! Definition of the function pointer for the root finder +/*! + * see vcsUtil_root1d for a definition of how to use this. + */ +typedef double(*VCS_FUNC_PTR)(double xval, double Vtarget, + int varID, void* fptrPassthrough, + int* err); + +//! One dimensional root finder +/*! + * + * This root finder will find the root of a one dimensional + * equation + * + * \f[ + * f(x) = 0 + * \f] + * where x is a bounded quantity: \f$ x_{min} < x < x_max \f$ + * + * The functional to be minimized must have the following call + * structure: + * + * @verbatim + typedef double (*VCS_FUNC_PTR)(double xval, double Vtarget, + int varID, void *fptrPassthrough, + int *err); @endverbatim + * + * xval is the current value of the x variable. Vtarget is the + * requested value of f(x), usually 0. varID is an integer + * that is passed through. fptrPassthrough is a void pointer + * that is passed through. err is a return error indicator. + * err = 0 is the norm. anything else is considered a fatal + * error. + * The return value of the function is the current value of + * f(xval). + * + * @param xmin Minimum permissible value of the x variable + * @param xmax Maximum permissible value of the x paramerer + * @param itmax Maximum number of iterations + * @param func function pointer, pointing to the function to be + * minimized + * @param fptrPassthrough Pointer to void that gets passed through + * the rootfinder, unchanged, to the func. + * @param FuncTargVal Target value of the function. This is usually set + * to zero. + * @param varID Variable ID. This is usually set to zero. + * @param xbest Pointer to the initial value of x on input. On output + * This contains the root value. + * @param printLvl Print level of the routine. + * + * + * Following is a nontrial example for vcs_root1d() in which the position of a + * cylinder floating on the water is calculated. + * + * @verbatim + #include + #include + + #include "equil/vcs_internal.h" + + const double g_cgs = 980.; + const double mass_cyl = 0.066; + const double diam_cyl = 0.048; + const double rad_cyl = diam_cyl / 2.0; + const double len_cyl = 5.46; + const double vol_cyl = Pi * diam_cyl * diam_cyl / 4 * len_cyl; + const double rho_cyl = mass_cyl / vol_cyl; + const double rho_gas = 0.0; + const double rho_liq = 1.0; + const double sigma = 72.88; + // Contact angle in radians + const double alpha1 = 40.0 / 180. * Pi; + + double func_vert(double theta1, double h_2, double rho_c) { + double f_grav = - Pi * rad_cyl * rad_cyl * rho_c * g_cgs; + double tmp = rad_cyl * rad_cyl * g_cgs; + double tmp1 = theta1 + sin(theta1) * cos(theta1) - 2.0 * h_2 / rad_cyl * sin(theta1); + double f_buoy = tmp * (Pi * rho_gas + (rho_liq - rho_gas) * tmp1); + double f_sten = 2 * sigma * sin(theta1 + alpha1 - Pi); + double f_net = f_grav + f_buoy + f_sten; + return f_net; + } + double calc_h2_farfield(double theta1) { + double rhs = sigma * (1.0 + cos(alpha1 + theta1)); + rhs *= 2.0; + rhs = rhs / (rho_liq - rho_gas) / g_cgs; + double sign = -1.0; + if (alpha1 + theta1 < Pi) sign = 1.0; + double res = sign * sqrt(rhs); + double h2 = res + rad_cyl * cos(theta1); + return h2; + } + double funcZero(double xval, double Vtarget, int varID, void *fptrPassthrough, int *err) { + double theta = xval; + double h2 = calc_h2_farfield(theta); + double fv = func_vert(theta, h2, rho_cyl); + return fv; + } + int main () { + double thetamax = Pi; + double thetamin = 0.0; + int maxit = 1000; + int iconv; + double thetaR = Pi/2.0; + int printLvl = 4; + + iconv = VCSnonideal::vcsUtil_root1d(thetamin, thetamax, maxit, + funcZero, + (void *) 0, 0.0, 0, + &thetaR, printLvl); + printf("theta = %g\n", thetaR); + double h2Final = calc_h2_farfield(thetaR); + printf("h2Final = %g\n", h2Final); + return 0; + } @endverbatim + * + */ +int vcsUtil_root1d(double xmin, double xmax, size_t itmax, VCS_FUNC_PTR func, + void* fptrPassthrough, + double FuncTargVal, int varID, double* xbest, + int printLvl = 0); + +//! Returns the system wall clock time in seconds +/*! + * @return time in seconds. + */ +double vcs_second(); + +//! This define turns on using memset and memcpy. I have not run into +//! any systems where this is a problem. It's the fastest way to do +//! low lvl operations where applicable. There are alternative routines +//! available if this ever fails. +#define USE_MEMSET +#ifdef USE_MEMSET + +//! Zero a double vector +/*! + * @param vec_to vector of doubles + * @param length length of the vector to zero. + */ +inline void vcs_dzero(double* const vec_to, const size_t length) +{ + (void) memset((void*) vec_to, 0, length * sizeof(double)); +} + +//! Zero an int vector +/*! + * @param vec_to vector of ints + * @param length length of the vector to zero. + */ +inline void vcs_izero(int* const vec_to, const size_t length) +{ + (void) memset((void*) vec_to, 0, length * sizeof(int)); +} + +//! Copy a double vector +/*! + * @param vec_to Vector to copy into. This vector must be dimensioned + * at least as large as the vec_from vector. + * @param vec_from Vector to copy from + * @param length Number of doubles to copy. + */ +inline void vcs_dcopy(double* const vec_to, + const double* const vec_from, const size_t length) +{ + (void) memcpy((void*) vec_to, (const void*) vec_from, + (length) * sizeof(double)); +} + + +//! Copy an int vector +/*! + * @param vec_to Vector to copy into. This vector must be dimensioned + * at least as large as the vec_from vector. + * @param vec_from Vector to copy from + * @param length Number of int to copy. + */ +inline void vcs_icopy(int* const vec_to, + const int* const vec_from, const size_t length) +{ + (void) memcpy((void*) vec_to, (const void*) vec_from, + (length) * sizeof(int)); +} + +//! Zero a std double vector +/*! + * @param vec_to vector of doubles + * @param length length of the vector to zero. + */ +inline void vcs_vdzero(std::vector &vec_to, const size_t length) +{ + (void) memset((void*)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(double)); +} + +//! Zero a std int vector +/*! + * @param vec_to vector of ints + * @param length length of the vector to zero. + */ +inline void vcs_vizero(std::vector &vec_to, const size_t length) +{ + (void) memset((void*)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(int)); +} + +//! Copy one std double vector into another +/*! + * This is an inlined function that uses memcpy. memcpy is probably + * the fastest way to do this. This routine requires the vectors to be + * previously dimensioned appropriately. No error checking is done. + * + * @param vec_to Vector to copy into. This vector must be dimensioned + * at least as large as the vec_from vector. + * @param vec_from Vector to copy from + * @param length Number of doubles to copy. + */ +inline void vcs_vdcopy(std::vector & vec_to, + const std::vector & vec_from, size_t length) +{ + (void) memcpy((void*)&(vec_to[0]), (const void*) &(vec_from[0]), + (length) * sizeof(double)); +} + +//! Copy one std integer vector into another +/*! + * This is an inlined function that uses memcpy. memcpy is probably + * the fastest way to do this. This routine requires the + * + * @param vec_to Vector to copy into. This vector must be dimensioned + * at least as large as the vec_from vector. + * @param vec_from Vector to copy from + * @param length Number of integers to copy. + */ +inline void vcs_vicopy(std::vector & vec_to, + const std::vector & vec_from, const int length) +{ + (void) memcpy((void*)&(vec_to[0]), (const void*) &(vec_from[0]), + (length) * sizeof(int)); +} +#else +extern void vcs_dzero(double* const, const int); +extern void vcs_izero(int* const , const int); +extern void vcs_dcopy(double* const, const double* const, const int); +extern void vcs_icopy(int* const, const int* const, const int); +extern void vcs_vdzero(std::vector &vvv, const int len = -1); +extern void vcs_vizero(std::vector &vvv, const int len = -1); +void vcs_vdcopy(std::vector &vec_to, + const std::vector vec_from, const int len = -1); +void vcs_vicopy(std::vector &vec_to, + const std::vector vec_from, const int len = -1); +#endif + +//! determine the l2 norm of a vector of doubles +/*! + * @param vec vector of doubles + * + * @return Returns the l2 norm of the vector + */ +double vcs_l2norm(const std::vector vec); + +//! Finds the location of the maximum component in a double vector +/*! + * @param x pointer to a vector of doubles + * @param xSize pointer to a vector of doubles used as a multiplier + * to x[] + * @param j lowest index to search from + * @param n highest index to search from + * @return Return index of the greatest value on X(i) searched + * j <= i < n + */ +size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n); + +//! Returns the maximum integer in a list +/*! + * @param vector pointer to a vector of ints + * @param length length of the integer vector + * + * @return returns the max integer value in the list + */ +int vcs_max_int(const int* vector, int length); + +//! Prints a line consisting of mutliple occurances of the same string +/*! + * This prints a string num times, and then terminate with a + * end of line character + * + * @param str C string that is null terminated + * @param num number of times the string is to be printed + */ +void vcs_print_line(const char* str, int num); + +//! Returns a const char string representing the type of the +//! species given by the first argument +/*! + * @param speciesStatus Species status integer representing the type + * of the species. + * @param length Maximum length of the string to be returned. + * Shorter values will yield abbreviated strings. + * Defaults to a value of 100. + */ +const char* vcs_speciesType_string(int speciesStatus, int length = 100); + +//! Print a string within a given space limit +/*! + * This routine limits the amount of the string that will be printed to a + * maximum of "space" characters. Printing is done to + * to Cantera's writelog() function. + * + * @param str String, which must be null terminated. + * @param space space limit for the printing. + * @param alignment Alignment of string within the space: + * - 0 centered + * - 1 right aligned + * - 2 left aligned + */ +void vcs_print_stringTrunc(const char* str, size_t space, int alignment); + +//! Simple routine to check whether two doubles are equal up to +//! roundoff error +/*! + * Currently it's set to check for 10 digits of + * relative accuracy. + * + * @param d1 first double + * @param d2 second double + * + * @return returns true if the doubles are "equal" and false otherwise + */ +bool vcs_doubleEqual(double d1, double d2); + +} + +#endif diff --git a/include/cantera/equil/vcs_prob.h b/include/cantera/equil/vcs_prob.h new file mode 100644 index 000000000..a8ef0325b --- /dev/null +++ b/include/cantera/equil/vcs_prob.h @@ -0,0 +1,372 @@ +/** + * @file vcs_prob.h + * Header for the Interface class for the vcs thermo equilibrium solver package, + */ +/* + * Copyright (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + +#ifndef _VCS_PROB_H +#define _VCS_PROB_H + +#include "vcs_DoubleStarStar.h" +#include "vcs_IntStarStar.h" +#include "cantera/equil/vcs_defs.h" +#include +#include + +namespace VCSnonideal +{ + +class vcs_VolPhase; +class VCS_SPECIES_THERMO; + +//! Interface class for the vcs thermo equilibrium solver package, +//! which generally describes the problem to be solved. +/*! + * HKM add: + * HaveEstimate -> 0 no estimate, or estimate that doesn' satisfy elem + * abundances + * 1 have an estimate that satisfies elem_abund. + * 2 Have an estimate that minimizes a subproblem + * and satisfies elem abund. + * solnFound -> True, soln to current problem found and included here + * False, soln has not been found. + */ +class VCS_PROB +{ +public: + + //! Problem type. I.e., the identity of what is held constant. + /*! + * Currently, T and P are held constant, and this input + * is ignored + */ + int prob_type; + + //! Total number of species in the problems + size_t nspecies; + + //! Species number used to malloc data structures + size_t NSPECIES0; + + //! Number of element constraints in the equilibrium problem + size_t ne; + + //! Number of element constraints used to malloc data structures + //! involving elements + size_t NE0; + + //! Number of phases in the problem + size_t NPhase; + + //! Number of phases used to malloc data structures + size_t NPHASE0; + + //! Vector of chemical potentials of the species + /*! + * This is a calculated output quantity + * length = number of species + * units = m_VCS_UnitsFormat; + */ + std::vector m_gibbsSpecies; + + //! Total number of moles of the kth species. + /*! + * This is both an input and an output variable. + * On input, this is an estimate of the mole numbers. + * The actual element abundance vector contains the problem specification. + * + * On output, this contains the solution for the total number of moles + * of the kth species. + * + * units = m_VCS_UnitsFormat + */ + std::vector w; + + //! Mole fraction vector + /*! + * This is a calculated vector, calculated from w[] + * length number of species. + * -> Take out? -> No, useful for storage of a quantity often needed + */ + std::vector mf; + + //! Element abundances for jth element + /*! + * This is input from the input file and is considered a constant from + * thereon within the vcs_solve_TP(). + * units = m_VCS_UnitsFormat + */ + std::vector gai; + + //! Formula Matrix for the problem + /*! + * FormulaMatrix[j][kspec] = Number of elements, j, in the kspec + * species + */ + DoubleStarStar FormulaMatrix; + + //! Specifies the species unknown type + /*! + * There are two types. One is the straightforward + * species, with the mole number w[k], as the + * unknown. The second is the an interfacial + * voltage where w[k] refers to the interfacial + * voltage in volts. + * These species types correspond to metalic + * electrons corresponding to electrodes. + * The voltage and other interfacial conditions + * sets up an interfacial current, which is + * set to zero in this initial treatment. + * Later we may have non-zero interfacial currents. + */ + std::vector SpeciesUnknownType; + + //! Temperature (Kelvin) + /*! + * Specification of the temperature for the equilibrium problem + */ + double T; + + //! Pressure + /*! + * units given by m_VCS_UnitsFormat + * -> are now PA + */ + double PresPA; + + //! Volume of the entire system + /*! + * units given by m_VCS_UnitsFormat + * Note, this is an output variable atm + */ + double Vol; + + //! Partial Molar Volumes of species + /*! + * This is a calculated vector, calculated from w[] + * length number of species. + * -> Take out? -> No, useful for storage of a quantity often needed + */ + std::vector VolPM; + + //! Units for the chemical potential data, pressure data, volume, + //! and species amounts + /*! + * All internally stored quantities will have these units. Also, printed + * quantitities will display in these units. + * + * Chem_Pot Pres vol moles + * ---------------------------------------------------------------------- + * -1 VCS_UNITS_KCALMOL = kcal/mol atm cm**3 gmol + * 0 VCS_UNITS_UNITLESS = MU / RT -> no units atm cm**3 gmol + * 1 VCS_UNITS_KJMOL = kJ / mol atm cm**3 gmol + * 2 VCS_UNITS_KELVIN = KELVIN -> MU / R atm cm**3 gmol + * 3 VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol + * ---------------------------------------------------------------------- + * + * see vcs_defs.h for more information + */ + int m_VCS_UnitsFormat; + + //! Specification of the initial estimate method + /*! + * iest = Initial estimate: 0 user estimate + * 1 user estimate if satisifies elements + * -1 machine estimate + */ + int iest; + + //! Tolerance requirement for major species + double tolmaj; + + //! Tolerance requirement for minor species + double tolmin; + + //! Mapping between the species and the phases + std::vector PhaseID; + + //! Vector of strings containing the species names + std::vector SpName; + + //! vector of strings containing the element names + std::vector ElName; + + //! vector of Element types + std::vector m_elType; + + //! Specifies whether an element constraint is active + /*! + * The default is true + * Length = nelements + */ + std::vector ElActive; + + //! Molecular weight of species + /*! + * WtSpecies[k] = molecular weight of species in gm/mol + */ + std::vector WtSpecies; + + //! Charge of each species + std::vector Charge; + + //! Array of phase structures + std::vector VPhaseList; + + // String containing the title of the run + std::string Title; + + //! Vector of pointers to thermo structures which identify the model + //! and parameters for evaluating the thermodynamic + //! functions for that particular species + std::vector SpeciesThermo; + + //! Number of iterations + /*! + * This is an output variable + */ + int m_Iterations; + + //! Number of basis optimizations used + /*! + * This is an output variable + */ + int m_NumBasisOptimizations; + + //! Print level for print routines + int m_printLvl; + + //! Debug print lvl + int vcs_debug_print_lvl; + + //! Constructor + /*! + * This constructor initializes the sizes within the object + * to parameter values. + * + * @param nsp number of species + * @param nel number of elements + * @param nph number of phases + */ + VCS_PROB(size_t nsp, size_t nel, size_t nph); + + //! Destructor + ~VCS_PROB(); + + //! Resizes all of the phase lists within the structure + /*! + * Note, this doesn't change the number of phases in the problem. + * It will change NPHASE0 if nsp is greater than NPHASE0. + * + * @param nPhase size to dimension all the phase lists to + * @param force If true, this will dimension the size to be equal to nPhase + * even if nPhase is less than the current value of NPHASE0 + */ + void resizePhase(size_t nPhase, int force); + + //! Resizes all of the species lists within the structure + /*! + * Note, this doesn't change the number of species in the problem. + * It will change NSPECIES0 if nsp is greater than NSPECIES0. + * + * @param nsp size to dimension all the species lists to + * @param force If true, this will dimension the size to be equal to nsp + * even if nsp is less than the current value of NSPECIES0 + */ + void resizeSpecies(size_t nsp, int force); + + //! Resizes all of the element lists within the structure + /*! + * Note, this doesn't change the number of element constraints in the problem. + * It will change NE0 if nel is greater than NE0. + * + * @param nel size to dimension all the elements lists + * @param force If true, this will dimension the size to be equal to nel + * even if nel is less than the current value of NEL0 + */ + void resizeElements(size_t nel, int force); + + + //! Calculate the element abundance vector + /*! + * Calculates the element abundance vectors from the mole + * numbers + */ + void set_gai(); + + //! Print out the problem specification in all generality + //! as it currently exists in the VCS_PROB object + /*! + * @param print_lvl Parameter lvl for printing + * 0 - no printing + * 1 - all printing + */ + void prob_report(int print_lvl); + + //! Add elements to the local element list + /*! + * This routine sorts through the elements defined in the + * vcs_VolPhase object. It then adds the new elements to + * the VCS_PROB object, and creates a global map, which is + * stored in the vcs_VolPhase object. + * Id and matching of elements is done strictly via the element name, + * with case not mattering. + * + * The routine also fills in the position of the element + * in the vcs_VolPhase object's ElGlobalIndex field. + * + * @param volPhase Object containing the phase to be added. + * The elements in this phase are parsed for + * addition to the global element list + */ + void addPhaseElements(vcs_VolPhase* volPhase); + + + //! This routine resizes the number of elements in the VCS_PROB object by + //! adding a new element to the end of the element list + /*! + * The element name is added. Formula vector entries ang element + * abundances for the new element are set to zero. + * + * Returns the index number of the new element. + * + * @param elNameNew New name of the element + * @param elType Type of the element + * @param elactive boolean indicating whether the element is active + * + * @return returns the index number of the new element + */ + size_t addElement(const char* elNameNew, int elType, int elactive); + + + //! This routines adds entries for the formula matrix for one species + /*! + * This routines adds entries for the formula matrix for this object + * for one species + * + * This object also fills in the index filed, IndSpecies, within + * the volPhase object. + * + * @param volPhase object containing the species + * @param k Species number within the volPhase k + * @param kT global Species number within this object + * + */ + size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT); + + void reportCSV(const std::string& reportFile); + + //! Set the debug level + /*! + * @param vcs_debug_print_lvl input debug level + */ + void setDebugPrintLvl(int vcs_debug_print_lvl); +}; + +} + +#endif diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h new file mode 100644 index 000000000..af91a8b2b --- /dev/null +++ b/include/cantera/equil/vcs_solve.h @@ -0,0 +1,2084 @@ +/** + * @file vcs_solve.h + * Header file for the internal object that holds the vcs equilibrium problem + * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). + */ +/* + * Copyright (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ + + +#ifndef _VCS_SOLVE_H +#define _VCS_SOLVE_H + +/* +* Index of Symbols +* ------------------- +* irxn -> refers to the species or rxn between the species and +* the components in the problem +* k -> refers to the species +* j -> refers to the element or component +* +* ### -> to be eliminated +*/ +#include +#include + +#include "cantera/base/ct_defs.h" +#include "cantera/equil/vcs_defs.h" +#include "cantera/equil/vcs_DoubleStarStar.h" +#include "cantera/equil/vcs_IntStarStar.h" +#include "cantera/equil/vcs_internal.h" + +namespace VCSnonideal +{ +/* + * Forward references + */ +class vcs_VolPhase; +class VCS_SPECIES_THERMO; +class VCS_PROB; +class VCS_COUNTERS; + + +//! This is the main structure used to hold the internal data +//! used in vcs_solve_TP(), and to solve TP systems. +/*! + * The indices of information in this + * structure may change when the species basis changes or when + * phases pop in and out of existence. Both of these operations + * change the species ordering. + * + */ +class VCS_SOLVE +{ +public: + //! Constructor for the VCS_SOLVE class + VCS_SOLVE(); + + //! Destructor + ~VCS_SOLVE(); + + + //! Initialize the sizes within the VCS_SOLVE object + /*! + * This resizes all of the internal arrays within the object. This routine + * operates in two modes. If all of the parameters are the same as it + * currently exists in the object, nothing is done by this routine; a quick + * exit is carried out and all of the data in the object persists. + * + * IF any of the parameters are different than currently exists in the + * object, then all of the data in the object must be redone. It may not + * be zeroed, but it must be redone. + * + * @param nspecies0 Number of species within the object + * @param nelements Number of element constraints within the problem + * @param nphase0 Number of phases defined within the problem. + * + */ + void vcs_initSizes(const size_t nspecies0, const size_t nelements, const size_t nphase0); + + //! Solve an equilibrium problem + /*! + * This is the main interface routine to the equilibrium solver + * + * Input: + * @param vprob Object containing the equilibrium Problem statement + * + * @param ifunc Determines the operation to be done: Valid values: + * 0 -> Solve a new problem by initializing structures + * first. An initial estimate may or may not have + * been already determined. This is indicated in the + * VCS_PROB structure. + * 1 -> The problem has already been initialized and + * set up. We call this routine to resolve it + * using the problem statement and + * solution estimate contained in + * the VCS_PROB structure. + * 2 -> Don't solve a problem. Destroy all the private + * structures. + * + * @param ipr Printing of results + * ipr = 1 -> Print problem statement and final results to + * standard output + * 0 -> don't report on anything + * @param ip1 Printing of intermediate results + * IP1 = 1 -> Print intermediate results. + * + * @param maxit Maximum number of iterations for the algorithm + * + * Output: + * + * @return + * nonzero value: failure to solve the problem at hand. + * zero : success + */ + int vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit); + + //! Main routine that solves for equilibrium at constant T and P + //! using a variant of the VCS method + /*! + * This is the main routine taht solves for equilibrium at constant T and P + * using a variant of the VCS method. Nonideal phases can be accommodated + * as well. + * + * Any number of single-species phases and multi-species phases + * can be handled by the present version. + * + * Input + * ------------ + * @param print_lvl 1 -> Print results to standard output + * 0 -> don't report on anything + * + * @param printDetails 1 -> Print intermediate results. + * + * @param maxit Maximum number of iterations for the algorithm + * + * @return 0 = Equilibrium Achieved + * 1 = Range space error encountered. The element abundance criteria are + * only partially satisfied. Specifically, the first NC= (number of + * components) conditions are satisfied. However, the full NE + * (number of elements) conditions are not satisfied. The equilibrirum + * condition is returned. + * -1 = Maximum number of iterations is exceeded. Convergence was not + * found. + */ + int vcs_solve_TP(int print_lvl, int printDetails, int maxit); + + + int vcs_PS(VCS_PROB* vprob, int iph, int printLvl, double& feStable); + + void vcs_reinsert_deleted(size_t kspec); + + //! Choose the optimum species basis for the calculations + /*! + * Choose the optimum component species basis for the calculations. + * This is done by choosing the species with the largest mole fraction + * not currently a linear combination of the previous components. + * Then, calculate the stoichiometric coefficient matrix for that + * basis. + * + * Rearranges the solution data to put the component data at the + * front of the species list. + * + * Then, calculates m_stoichCoeffRxnMatrix[irxn][jcomp] the formation reactions + * for all noncomponent species in the mechanism. + * Also calculates DNG(I) and DNL(I), the net mole change for each + * formation reaction. + * Also, initializes IR(I) to the default state. + * + * Input + * --------- + * @param doJustCompoents If true, the m_stoichCoeffRxnMatrix[][] and + * m_deltaMolNumPhase[] are not calculated. + * + * @param aw Vector of mole fractions which will be used to construct an + * optimal basis from. + * + * @param sa Gramm-Schmidt orthog work space (nc in length) sa[j] + * @param ss Gramm-Schmidt orthog work space (nc in length) ss[j] + * @param sm QR matrix work space (nc*ne in length) sm[i+j*ne] + * @param test This is a small negative number dependent upon whether + * an estimate is supplied or not. + * + * Output + * --------- + * @param usedZeroedSpecies = If true, then a species with a zero concentration + * was used as a component. The problem may be + * converged. Or, the problem may have a range space + * error and may not have a proper solution. + * + * Internal Variables calculated by this routine: + * ----------------------------------------------- + * + * m_numComponents + * Number of component species + * + * component species + * This routine calculates the m_numComponent species. It switches + * their positions in the species vector so that they occupy + * the first m_numComponent spots in the species vector. + * + * m_stoichCoeffRxnMatrix[irxn][jcomp] + * Stoichiometric coefficient matrix for the reaction mechanism + * expressed in Reduced Canonical Form. + * jcomp refers to the component number, and irxn + * refers to the irxn_th non-component species. + * + * m_deltaMolNumPhase[irxn] + * Change in the number of total number of moles of species in all phases + * due to the noncomponent formation reaction, irxn. + * + * m_deltaMolNumPhase[irxn][iphase] + * Change in the number of moles in phase, iphase, due to the + * noncomponent formation reaction, irxn. + * + * m_phaseParticipation[irxn] + * This is 1 if the phase, iphase, participates in the + * formation reaction, irxn, and zero otherwise. + * + * @return Returns VCS_SUCCESS if everything went ok. Returns something else if + * there is a problem. + */ + int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], + double ss[], double test, bool* const usedZeroedSpecies); + + //! Choose a species to test for the next component + /*! + * We make the choice based on testing (molNum[i] * spSize[i]) for its maximum value. + * Preference for single species phases is also made. + * + * @param molNum Mole number vector + * @param j index into molNum[] that indicates where the search will start from + * Previous successful components are swapped into the fronto of + * molNum[]. + * @param n Length of molNum[] + */ + size_t vcs_basisOptMax(const double* const molNum, const size_t j, const size_t n); + + //! Evaluate the species category for the indicated species + /*! + * All evaluations are done using the "old" version of the solution. + * + * @param kspec Species to be evaluated + * + * @return Returns the calculated species type + */ + int vcs_species_type(const size_t kspec) const; + + bool vcs_evaluate_speciesType(); + + //! We calculate the dimensionless chemical potentials of all species + //! in a single phase. + /*! + * We calculate the dimensionless chemical potentials of all species + * in a single phase. + * + * Note, for multispecies phases which are currently zeroed out, + * the chemical potential is filled out with the standard chemical + * potential. + * + * For species in multispecies phases whose concentration is zero, + * we need to set the mole fraction to a very low value. + * Its chemical potential + * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration + * to keep numbers positive. + * + * Formula: + * --------------- + * + * Ideal Mixtures: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * + * ( This is equivalent to the adding the log of the + * mole fraction onto the standard chemical + * potential. ) + * + * Non-Ideal Mixtures: + * ActivityConvention = 0: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * ( This is equivalent to the adding the log of the + * mole fraction multiplied by the activity coefficient + * onto the standard chemical potential. ) + * + * ActivityConvention = 1: -> molality activity formulation + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) + * - ln(Mnaught * m_units) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * note: m_SSfeSpecies(I) is the molality based standard state. + * However, ActCoeff[I] is the molar based activity coefficient + * We have used the formulas; + * + * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] + * where Xmol[N] is the mole fraction of the solvent + * ActCoeff_M[I] is the molality based act coeff. + * + * note: This is equivalent to the "normal" molality formulation: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff_M[I] * m(I)) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase] + * where m[I] is the molality of the ith solute + * + * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) + * + * + * note: z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction + * of i in the phase. + * + * + * NOTE: + * As per the discussion in vcs_dfe(), for small species where the mole + * fraction is small: + * + * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF + * + * The chemical potential is calculated as: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) + * + * Input + * -------- + * iph : Phase to be calculated + * molNum(i) : Number of moles of species i + * (VCS species order) + * ff : standard state chemical potentials. These are the + * chemical potentials of the standard states at + * the same T and P as the solution. + * (VCS species order) + * Output + * ------- + * ac[] : Activity coefficients for species in phase + * (VCS species order) + * mu_i[] : Dimensionless chemical potentials for phase species + * (VCS species order) + * + */ + void vcs_chemPotPhase(const int stateCalc, const size_t iph, const double* const molNum, + double* const ac, double* const mu_i, + const bool do_deleted = false); + + //! Calculate the dimensionless chemical potentials of all species or + //! of certain groups of species, at a fixed temperature and pressure. + /*! + * We calculate the dimensionless chemical potentials of all species + * or certain groups of species here, at a fixed temperature and pressure, + * for the input mole vector z[] in the parameter list. + * Nondimensionalization is achieved by division by RT. + * + * Note, for multispecies phases which are currently zeroed out, + * the chemical potential is filled out with the standard chemical + * potential. + * + * For species in multispecies phases whose concentration is zero, + * we need to set the mole fraction to a very low value. + * Its chemical potential + * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration + * to keep numbers positive. + * + * + * Formula: + * --------------- + * + * Ideal Mixtures: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * + * ( This is equivalent to the adding the log of the + * mole fraction onto the standard chemical + * potential. ) + * + * Non-Ideal Mixtures: -> molar activity formulation + * ActivityConvention = 0: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * + * ( This is equivalent to the adding the log of the + * mole fraction multiplied by the activity coefficient + * onto the standard chemical potential. ) + * + * note: z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction + * of i in the phase. + * + * ActivityConvention = 1: -> molality activity formulation + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) + * - ln(Mnaught * m_units) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * + * note: m_SSfeSpecies(I) is the molality based standard state. + * However, ActCoeff[I] is the molar based activity coefficient + * We have used the formulas; + * + * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] + * where Xmol[N] is the mole fraction of the solvent + * ActCoeff_M[I] is the molality based act coeff. + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff_M[I] * m(I)) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * where m[I] is the molality of the ith solute + * + * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) + * + * + * Handling of Small Species: + * ------------------------------ + * As per the discussion above, for small species where the mole + * fraction + * + * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF + * + * The chemical potential is calculated as: + * + * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) + * + * Species in the following categories are treated as "small species" + * + * - VCS_SPECIES_DELETED + * - VCS_SPECIES_ACTIVEBUTZERO + * . + * + * Handling of Small Species: + * ------------------------------ + * For species in multispecies phases which are currently not active, the + * treatment is different. These species are in the following species categories: + * + * - VCS_SPECIES_ZEROEDMS + * - VCS_SPECIES_ZEROEDPHASE + * . + * + * For these species, the ln( ActCoeff[I] X[I]) term is + * dropped altogether. The following equation is used. + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * + * + * Handling of "Species" Representing Interfacial Voltages + * --------------------------------------------------------- + * + * These species have species types of VCS_SPECIES_TYPE_INTERFACIALVOLTAGE + * The chemical potentials for these "species" refer to electrons in + * metal electrodes. They have the following formula + * + * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT + * + * F is Faraday's constant. + * R = gas constant + * T = temperature + * V = potential of the interface = phi_electrode - phi_solution + * + * For these species, the solution vector unknown, z[I], is V, the phase voltage, in volts. + * + * Input + * -------- + * @param ll Determine which group of species gets updated + * ll = 0: Calculate for all species + * < 0: calculate for components and for major non-components + * 1: calculate for components and for minor non-components + * + * @param lbot Restricts the calculation of the chemical potential + * to the species between LBOT <= i < LTOP. Usually + * LBOT and LTOP will be equal to 0 and MR, respectively. + * @param ltop Top value of the loops + * + * + * @param stateCalc Determines whether z is old or new or tentative: + * 1: Use the tentative values for the total number of + * moles in the phases, i.e., use TG1 instead of TG etc. + * 0: Use the base values of the total number of + * moles in each system. + * + * Also needed: + * ff : standard state chemical potentials. These are the + * chemical potentials of the standard states at + * the same T and P as the solution. + * tg : Total Number of moles in the phase. + */ + void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop); + + //! Print out a table of chemical potentials + /*! + * @param vcsState Determines where to get the mole numbers from. + * - VCS_STATECALC_OLD -> from m_molNumSpecies_old + * - VCS_STATECALC_NEW -> from m_molNumSpecies_new + */ + void vcs_printSpeciesChemPot(const int stateCalc) const; + + //! This routine uploads the state of the system into all of the + //! vcs_VolumePhase objects in the current problem. + /*! + * @param vcsState Determines where to get the mole numbers from. + * - VCS_STATECALC_OLD -> from m_molNumSpecies_old + * - VCS_STATECALC_NEW -> from m_molNumSpecies_new + */ + void vcs_updateVP(const int stateCalc); + + //! Utility function that evaluates whether a phase can be popped + //! into existence + /*! + * @param iphasePop id of the phase, which is currently zeroed, + * + * @return Returns true if the phase can come into existence + * and false otherwise. + */ + bool vcs_popPhasePossible(const size_t iphasePop) const; + + + //! Determine the list of problems that need to be checked to see if there are any phases pops + /*! + * This routine evaluates and fills in the following quantities + * phasePopProblemLists_ + * + * @return Returns the number of problems that must be checked. + */ + int vcs_phasePopDeterminePossibleList(); + + + + + + + //! Decision as to whether a phase pops back into existence + /*! + * @param phasePopPhaseIDs Vector containing the phase ids of the phases + * that will be popped this step. + * + * @return returns the phase id of the phase that pops back into + * existence. Returns -1 if there are no phases + */ + size_t vcs_popPhaseID(std::vector &phasePopPhaseIDs); + + //! Calculates the deltas of the reactions due to phases popping + //! into existence + /*! + * @param iphasePop Phase id of the phase that will come into existence + * + * @return Returns an int representing the status of the step + * - 0 : normal return + * - 1 : A single species phase species has been zeroed out + * in this routine. The species is a noncomponent + * - 2 : Same as one but, the zeroed species is a component. + */ + int vcs_popPhaseRxnStepSizes(const size_t iphasePop); + + //! Calculates formation reaction step sizes. + /*! + * This is equation 6.4-16, p. 143 in Smith and Missen. + * + * Output + * ------- + * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers + * to the irxn'th species + * formation reaction. This adjustment is for species + * irxn + M, where M is the number of components. + * + * Special branching occurs sometimes. This causes the component basis + * to be reevaluated + * + * @param forceComponentRecalc integer flagging whether a component recalculation needs + * to be carried out. + * @param kSpecial species number of phase being zeroed. + * + * @return Returns an int representing which phase may need to be zeroed + */ + size_t vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial); + + //! Calculates the total number of moles of species in all phases. + /*! + * Calculates the total number of moles in all phases and updates + * the variable m_totalMolNum. + * Reconciles Phase existence flags with total moles in each phase. + */ + double vcs_tmoles(); +#ifdef DEBUG_MODE + void check_tmoles() const; +#endif + + //! This subroutine calculates reaction free energy changes for + //! all noncomponent formation reactions. + /*! + * Formation reactions are + * reactions which create each noncomponent species from the component + * species. m_stoichCoeffRxnMatrix[irxn][jcomp] are the stoichiometric + * coefficients for these reactions. A stoichiometric coefficient of + * one is assumed for species irxn in this reaction. + * + * INPUT + * @param l + * L < 0 : Calculate reactions corresponding to + * major noncomponent and zeroed species only + * L = 0 : Do all noncomponent reactions, i, between + * 0 <= i < irxnl + * L > 0 : Calculate reactions corresponding to + * minor noncomponent and zeroed species only + * + * @param doDeleted Do deleted species + * @param vcsState Calculate deltaG corresponding to either old or new + * free energies + * @param alterZeroedPhases boolean indicating whether we should + * add in a special section for zeroed phases. + * + * Note we special case one important issue. + * If the component has zero moles, then we do not + * allow deltaG < 0.0 for formation reactions which + * would lead to the loss of more of that same component. + * This dG < 0.0 condition feeds back into the algorithm in several + * places, and leads to a infinite loop in at least one case. + */ + void vcs_deltag(const int l, const bool doDeleted, const int vcsState, + const bool alterZeroedPhases = true); + + void vcs_printDeltaG(const int stateCalc); + + //! Calculate deltag of formation for all species in a single phase. + /*! + * Calculate deltag of formation for all species in a single + * phase. It is assumed that the fe[] is up to date for all species. + * Howevever, if the phase is currently zereoed out, a subproblem + * is calculated to solve for AC[i] and pseudo-X[i] for that + * phase. + * + * @param iphase phase index of the phase to be calculated + * @param doDeleted boolean indicating whether to do deleted + * species or not + * @param stateCalc integer describing which set of free energies + * to use and where to stick the results. + * @param alterZeroedPhases boolean indicating whether we should + * add in a special section for zeroed phases. + * + * NOTE: this is currently not used used anywhere. + * It may be in the future? + */ + void vcs_deltag_Phase(const size_t iphase, const bool doDeleted, + const int stateCalc, const bool alterZeroedPhases = true); + + //! Swaps the indices for all of the global data for two species, k1 + //! and k2. + /*! + * + * @param ifunc: If true, switch the species data and the noncomponent reaction + * data. This must be called for a non-component species only. + * If false, switch the species data only. Typically, we use this + * option when determining the component species and at the + * end of the calculation, when we want to return unscrambled + * results. All rxn data will be out-of-date. + * + * @param k1 First species index + * + * @param k2 Second species index + */ + void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2); + + + //! Birth guess returns the number of moles of a species + //! that is coming back to life. + /*! + * Birth guess returns the number of moles of a species + * that is coming back to life. + * Note, this routine is not applicable if the whole phase is coming + * back to life, not just one species in that phase. + * + * Do a minor alt calculation. But, cap the mole numbers at + * 1.0E-15. + * For SS phases use VCS_DELETE_SPECIES_CUTOFF * 100. + * + * The routine makes sure the guess doesn't reduce the concentration + * of a component by more than 1/3. Note this may mean that + * the vlaue coming back from this routine is zero or a + * very small number. + * + * + * @param kspec Species number that is coming back to life + * + * @return Returns the number of kmol that the species should + * have. + */ + double vcs_birthGuess(const int kspec); + + int vcs_solve_phaseStability(const int iphase, int ifunc, double& funcval, int print_lvl); + + //! Main program to test whether a deleted phase should be brought + //! back into existence + /*! + * + * @param iph Phase id of the deleted phase + */ + double vcs_phaseStabilityTest(const size_t iph); + + //! Solve an equilibrium problem at a particular fixed temperature + //! and pressure + /*! + * The actual problem statement is assumed to be in the structure + * already. This is a wrapper around the solve_TP() function. + * In this wrapper, we nondimensionalize the system + * we calculate the standard state gibbs free energies of the + * species, and we decide whether to we need to use the + * initial guess algorithm. + * + * @param ipr = 1 -> Print results to standard output + * 0 -> don't report on anything + * @param ip1 = 1 -> Print intermediate results. + * 0 -> Dont print any intermediate results + * @param maxit Maximum number of iterations for the algorithm + * @param T Value of the Temperature (Kelvin) + * Param pres Value of the Pressure (units given by m_VCS_UnitsFormat variable + * + * @return Returns an integer representing the success of the algorithm + * 0 = Equilibrium Achieved + * 1 = Range space error encountered. The element abundance criteria are + * only partially satisfied. Specifically, the first NC= (number of + * components) conditions are satisfied. However, the full NE + * (number of elements) conditions are not satisfied. The equilibrirum + * condition is returned. + * -1 = Maximum number of iterations is exceeded. Convergence was not + * found. + */ + int vcs_TP(int ipr, int ip1, int maxit, double T, double pres); + + int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres); + + //! Initialize the chemical potential of single species phases + /*! + * For single species phases, initialize the chemical + * potential with the value of the standard state chemical + * potential. This value doesn't change during the calculation + */ + void vcs_fePrep_TP(); + + //! Calculation of the total volume and the partial molar volumes + /*! + * This function calculates the partial molar volume + * for all species, kspec, in the thermo problem + * at the temperature TKelvin and pressure, Pres, pres is in atm. + * And, it calculates the total volume of the combined system. + * + * Input + * --------------- + * @param tkelvin Temperature in kelvin() + * @param pres Pressure in Pascal + * @param w w[] is thevector containing the current mole numbers + * in units of kmol. + * + * Output + * ---------------- + * @param volPM[] For species in all phase, the entries are the + * partial molar volumes units of M**3 / kmol. + * + * @return The return value is the total volume of + * the entire system in units of m**3. + */ + double vcs_VolTotal(const double tkelvin, const double pres, + const double w[], double volPM[]); + + //! This routine is mostly concerned with changing the private data + //! to be consistent with what's needed for solution. It is called one + //! time for each new problem structure definition. + /*! + * This routine is always followed by vcs_prep(). Therefore, tasks + * that need to be done for every call to vcsc() should be placed in + * vcs_prep() and not in this routine. + * + * The problem structure refers to: + * + * the number and identity of the species. + * the formula matrix and thus the number of components. + * the number and identity of the phases. + * the equation of state + * the method and parameters for determining the standard state + * The method and parameters for determining the activity coefficients. + * + * Tasks: + * 0) Fill in the SSPhase[] array. + * 1) Check to see if any multispecies phases actually have only one + * species in that phase. If true, reassign that phase and species + * to be a single-species phase. + * 2) Determine the number of components in the problem if not already + * done so. During this process the order of the species is changed + * in the private data structure. All references to the species + * properties must employ the ind[] index vector. + * + * @param printLvl Print level of the routine + * + * @return the return code + * VCS_SUCCESS = everything went OK + * + */ + int vcs_prep_oneTime(int printLvl); + + //! Prepare the object for solution + /*! + * This routine is mostly concerned with changing the private data + * to be consistent with that needed for solution. It is called for + * every invocation of the vcs_solve() except for the cleanup invocation. + * + * Tasks: + * 1) Initialization of arrays to zero. + * + * return code + * VCS_SUCCESS = everything went OK + * VCS_PUB_BAD = There is an irreconcilable difference in the + * public data structure from when the problem was + * initially set up. + */ + int vcs_prep(); + + //! In this routine, we check for things that will cause the algorithm + //! to fail. + /*! + * We check to see if the problem is well posed. If it is not, we return + * false and print out error conditions. + * + * Current there is one condition. If all the element abundances are + * zero, the algorithm will fail. + * + * @param vprob VCS_PROB pointer to the definition of the equilibrium + * problem + * + * @return If true, the problem is well-posed. If false, the problem + * is not well posed. + */ + bool vcs_wellPosed(VCS_PROB* vprob); + + //! Rearrange the constraint equations represented by the Formula + //! Matrix so that the operational ones are in the front + /*! + * + * This subroutine handles the rearrangement of the constraint + * equations represented by the Formula Matrix. Rearrangement is only + * necessary when the number of components is less than the number of + * elements. For this case, some constraints can never be satisfied + * exactly, because the range space represented by the Formula + * Matrix of the components can't span the extra space. These + * constraints, which are out of the range space of the component + * Formula matrix entries, are migrated to the back of the Formula + * matrix. + * + * A prototypical example is an extra element column in + * FormulaMatrix[], + * which is identically zero. For example, let's say that argon is + * has an element column in FormulaMatrix[], but no species in the + * mechanism + * actually contains argon. Then, nc < ne. Also, without perturbation + * of FormulaMatrix[] vcs_basopt[] would produce a zero pivot + * because the matrix + * would be singular (unless the argon element column was already the + * last column of FormulaMatrix[]. + * This routine borrows heavily from vcs_basopt's algorithm. It + * finds nc constraints which span the range space of the Component + * Formula matrix, and assigns them as the first nc components in the + * formula matrix. This guarantees that vcs_basopt[] has a + * nonsingular matrix to invert. + * + * Other Variables + * @param aw aw[i[ Mole fraction work space (ne in length) + * @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length) + * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) + * @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length) + * + */ + int vcs_elem_rearrange(double* const aw, double* const sa, + double* const sm, double* const ss); + + //! Swaps the indices for all of the global data for two elements, ipos + //! and jpos. + /*! + * This function knows all of the element information with VCS_SOLVE, and + * can therefore switch element positions + * + * @param ipos first global element index + * @param jpos second global element index + */ + void vcs_switch_elem_pos(size_t ipos, size_t jpos); + + //! Calculates reaction adjustments using a full Hessian approximation + /*! + * Calculates reaction adjustments. This does what equation 6.4-16, p. 143 + * in Smith and Missen is suppose to do. However, a full matrix is + * formed and then solved via a conjugate gradient algorithm. No + * preconditioning is done. + * + * If special branching is warranted, then the program bails out. + * + * Output + * ------- + * DS(I) : reaction adjustment, where I refers to the Ith species + * Special branching occurs sometimes. This causes the component basis + * to be reevaluated + * return = 0 : normal return + * 1 : A single species phase species has been zeroed out + * in this routine. The species is a noncomponent + * 2 : Same as one but, the zeroed species is a component. + * + * Special attention is taken to flag cases where the direction of the + * update is contrary to the steepest descent rule. This is an important + * attribute of the regular vcs algorithm. We don't want to violate this. + * + * NOTE: currently this routine is not used. + */ + int vcs_rxn_adj_cg(void); + + //! Calculates the diagonal contribution to the Hessian due to + //! the dependence of the activity coefficients on the mole numbers. + /*! + * (See framemaker notes, Eqn. 20 - VCS Equations document) + * + * We allow the diagonal to be increased positively to any degree. + * We allow the diagonal to be decreased to 1/3 of the ideal solution + * value, but no more -> it must remain positive. + * + * NOTE: currently this routine is not used + */ + double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal); + + //! Calculates the diagonal contribution to the Hessian due to + //! the dependence of the activity coefficients on the mole numbers. + /*! + * (See framemaker notes, Eqn. 20 - VCS Equations document) + * + * NOTE: currently this routine is not used + */ + double vcs_Hessian_actCoeff_diag(size_t irxn); + + void vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS); + +#ifdef DEBUG_MODE + //! A line search algorithm is carried out on one reaction + /*! + * In this routine we carry out a rough line search algorithm + * to make sure that the m_deltaGRxn_new doesn't switch signs prematurely. + * + * @param irxn Reaction number + * @param dx_orig Original step length + * + * @param ANOTE Output character string stating the conclusions of the + * line search + * + */ + double vcs_line_search(const size_t irxn, const double dx_orig, + char* const ANOTE); +#else + double vcs_line_search(const size_t irxn, const double dx_orig); +#endif + + + //! Print out a report on the state of the equilibrium problem to + //! standard output. + /*! + * @param iconv Indicator of convergence, to be printed out in the report: + * - 0 converged + * - 1 range space error + * - -1 not converged + */ + int vcs_report(int iconv); + + //! Switch all species data back to the original order. + /*! + * This destroys the data based on reaction ordering. + */ + int vcs_rearrange(); + + //! Returns the multiplier for electric charge terms + /* + * This is basically equal to F/RT + * + * @param mu_units integer representing the dimensional units system + * @param TKelvin double Temperature in Kelvin + * + * @return Returns the value of F/RT + */ + double vcs_nondim_Farad(int mu_units, double TKelvin) const; + + //! Returns the multiplier for the nondimensionalization of the equations + /*! + * This is basically equal to RT + * + * @param mu_units integer representing the dimensional units system + * @param TKelvin double Temperature in Kelvin + * + * @return Returns the value of RT + */ + double vcs_nondimMult_TP(int mu_units, double TKelvin) const; + + //! Nondimensionalize the problem data + /*! + * Nondimensionalize the free energies using the divisor, R * T + * + * Essentially the internal data can either be in dimensional form + * or in nondimensional form. This routine switches the data from + * dimensional form into nondimensional form. + * + * What we do is to divide by RT. + * + * @todo Add a scale factor based on the total mole numbers. + * The algorithm contains hard coded numbers based on the + * total mole number. If we ever were faced with a problem + * with significantly different total kmol numbers than one + * the algorithm would have problems. + */ + void vcs_nondim_TP(); + + //! Redimensionalize the problem data + /*! + * Reddimensionalize the free energies using the multiplier R * T + * + * Essentially the internal data can either be in dimensional form + * or in nondimensional form. This routine switches the data from + * nondimensional form into dimensional form. + * + * What we do is to multiply by RT. + */ + void vcs_redim_TP(); + + //! Print the string representing the Chemical potential units + /*! + * This gets printed using plogf() + * + * @param unitsFormat Integer representing the units system + */ + void vcs_printChemPotUnits(int unitsFormat) const; + + //! Computes the current elemental abundances vector + /*! + * Computes the elemental abundances vector, m_elemAbundances[], and stores it + * back into the global structure + */ + void vcs_elab(); + + bool vcs_elabcheck(int ibound); + void vcs_elabPhase(size_t iphase, double* const elemAbundPhase); + int vcs_elcorr(double aa[], double x[]); + + + //! Create an initial estimate of the solution to the thermodynamic + //! equilibrium problem. + /*! + * @return Return value indicates success: + * - 0: successful initial guess + * - -1: Unsuccessful initial guess; the elemental abundances aren't + * satisfied. + */ + int vcs_inest_TP(); + +#ifdef ALTLINPROG + //! Extimate the initial mole numbers by constrained linear programming + /*! + * This is done by running + * each reaction as far forward or backward as possible, subject + * to the constraint that all mole numbers remain + * non-negative. Reactions for which \f$ \Delta \mu^0 \f$ are + * positive are run in reverse, and ones for which it is negative + * are run in the forward direction. The end result is equivalent + * to solving the linear programming problem of minimizing the + * linear Gibbs function subject to the element and + * non-negativity constraints. + */ + int vcs_setMolesLinProg(); +#endif + + double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles); + + //! Calculate the total dimensionless Gibbs free energy of a single phase + /*! + * -> Inert species are handled as if they had a standard free + * energy of zero and if they obeyed ideal solution/gas theory + * + * @param iphase ID of the phase + * @param w Species mole number vector for all species + * @param fe vector of partial molar free energies of all of the + * species + */ + double vcs_GibbsPhase(size_t iphase, const double* const w, + const double* const fe); + + //! Transfer the results of the equilibrium calculation back to VCS_PROB + /*! + * The VCS_PUB structure is returned to the user. + * + * @param pub Pointer to VCS_PROB object that will get the results of the + * equilibrium calculation transfered to it. + */ + int vcs_prob_update(VCS_PROB* pub); + + //! Fully specify the problem to be solved using VCS_PROB + /*! + * Use the contents of the VCS_PROB to specify the contents of the + * private data, VCS_SOLVE. + * + * @param pub Pointer to VCS_PROB that will be used to + * initialize the current equilibrium problem + */ + int vcs_prob_specifyFully(const VCS_PROB* pub); + + //! Specify the problem to be solved using VCS_PROB, incrementally + /*! + * Use the contents of the VCS_PROB to specify the contents of the + * private data, VCS_SOLVE. + * + * It's assumed we are solving the same problem. + * + * @param pub Pointer to VCS_PROB that will be used to + * initialize the current equilibrium problem + */ + int vcs_prob_specify(const VCS_PROB* pub); + +private: + + //! Zero out the concentration of a species. + /*! + * Zero out the concentration of a species. Make sure to conserve + * elements and keep track of the total moles in all phases. + * w[] + * m_tPhaseMoles_old[] + * + * @param kspec Species index + * + * @return: + * 1: succeeded + * 0: failed. + */ + int vcs_zero_species(const size_t kspec); + + //! Change a single species from active to inactive status + /*! + * Rearrange data when species is added or removed. The Lth species is + * moved to the back of the species vector. The back of the species + * vector is indicated by the value of MR, the current number of + * active species in the mechanism. + * + * @param kspec Species Index + * @return + * Returns 0 unless. + * The return is 1 when the current number of + * noncomponent species is equal to zero. A recheck of deleted species + * is carried out in the main code. + */ + int vcs_delete_species(const size_t kspec); + + //! This routine handles the bookkeepking involved with the + //! deletion of multiphase phases from the problem. + /*! + * When they are deleted, all of their species become active + * species, even though their mole numbers are set to zero. + * The routine does not make the decision to eliminate multiphases. + * + * Note, species in phases with zero mole numbers are still + * considered active. Whether the phase pops back into + * existence or not is checked as part of the main iteration + * loop. + * + * @param iph Phase to be deleted + * + * @return Returns whether the operation was successful or not + */ + bool vcs_delete_multiphase(const size_t iph); + + //! Change the concentration of a species by delta moles. + /*! + * Make sure to conserve elements and keep track of the total kmoles in all phases. + * + * @param kspec The species index + * @delta_ptr pointer to the delta for the species. This may change during + * the calculation + * + * @return + * 1: succeeded without change of dx + * 0: Had to adjust dx, perhaps to zero, in order to do the delta. + */ + int delta_species(const size_t kspec, double* const delta_ptr); + + //! Provide an estimate for the deleted species in phases that + //! are not zeroed out + /*! + * Try to add back in all deleted species. An estimate of the kmol numbers + * are obtained and the species is added back into the equation system, + * into the old state vector. + * + * This routine is called at the end of the calculation, just before + * returning to the user. + */ + size_t vcs_add_all_deleted(); + + //! Recheck deleted species in multispecies phases. + /*! + * We are checking the equation: + * + * sum_u = sum_j_comp [ sigma_i_j * u_j ] + * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) + * + * by first evaluating: + * + * DG_i_O = u_i_O - sum_u. + * + * Then, if TL is zero, the phase pops into existence if DG_i_O < 0. + * Also, if the phase exists, then we check to see if the species + * can have a mole number larger than VCS_DELETE_SPECIES_CUTOFF + * (default value = 1.0E-32). + * + */ + int vcs_recheck_deleted(); + + //! Recheck deletion condition for multispecies phases. + /*! + * We assume here that DG_i_0 has been calculated for deleted species correctly + * + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I]) + * - ln(Mnaught * m_units) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * sum_u = sum_j_comp [ sigma_i_j * u_j ] + * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) + * + * DG_i_0 = m_feSpecies(I) - sum_m{ a_i_m DG_m } + * + * + * by first evaluating: + * + * DG_i_O = u_i_O - sum_u. + * + * Then, the phase pops into existence iff + * + * phaseDG = 1.0 - sum_i{exp(-DG_i_O)} < 0.0 + * + * This formula works for both single species phases and for multispecies + * phases. It's an overkill for single species phases. + * + * @param iphase Phase index number + * + * @return Returns true if the phase is currently deleted + * but should be reinstated. Returns false otherwise. + * + * NOTE: this routine is currently not used in the code, and + * contains some basic changes that are incompatible. + * + * assumptions: + * 1) Vphase Existence is up to date + * 2) Vphase->IndSpecies is up to date + * 3) m_deltaGRxn_old[irxn] is up to date + */ + bool recheck_deleted_phase(const int iph); + + //! Minor species alternative calculation + /*! + * This is based upon the following approximation: + * The mole fraction changes due to these reactions don't affect + * the mole numbers of the component species. Therefore the following + * approximation is valid for a small component of an ideal phase: + * + * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I)) + * + * m_deltaGRxn_old contains the contribution from + * + * m_feSpecies_old(I) = + * m_SSfeSpecies(I) + + * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph)) + * Thus, + * + * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I)) + * + * Most of this section is mainly restricting the update to reasonable + * values. + * We restrict the update a factor of 1.0E10 up and 1.0E-10 down + * because we run into trouble with the addition operator due to roundoff + * if we go larger than ~1.0E15. Roundoff will then sometimes produce + * zero mole fractions. + * + * Note: This routine was generalized to incorporate + * nonideal phases and phases on the molality basis + * + * Input: + * ------ + * @param kspec The current species and corresponding formation + * reaction number. + * @param irxn The current species and corresponding formation + * reaction number. + * + * Output: + * --------- + * @param do_delete: BOOLEAN which if true on return, then we branch + * to the section that deletes a species from the + * current set of active species. + * + * @param dx The change in mole number + */ + double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete +#ifdef DEBUG_MODE + , char* ANOTE +#endif + ) const; + + //! This routine optimizes the minimization of the total gibbs free + //! energy by making sure the slope of the following functional stays + //! negative: + /*! + * The slope of the following functional is equivalent to the slope + * of the total Gibbs free energy of the system: + * + * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] ) + * + * along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0 0), + * does this code section kick in. It finds the point on the parabola + * where the slope is equal to zero. + * + */ + bool vcs_globStepDamp(); + + //! Switch rows and columns of a sqare matrix + /*! + * Switches the row and column of a matrix. + * So that after + * + * J[k1][j] = J_old[k2][j] and J[j][k1] = J_old[j][k2] + * J[k2][j] = J_old[k1][j] and J[j][k2] = J_old[j][k1] + * + * @param Jac Double pointer to the jacobiam + * @param k1 first row/column value to be switched + * @param k2 second row/column value to be switched + */ + void vcs_switch2D(double* const* const Jac, + const size_t k1, const size_t k2) const; + + //! Calculate the norm of a deltaGibbs free energy vector + /*! + * Positive DG for species which don't exist are ignored. + * + * @param dgLocal Vector of local delta G's. + */ + double l2normdg(double dg[]) const; + +#ifdef DEBUG_MODE + + //! Print out and check the elemental abundance vector + void prneav() const; + + void checkDelta1(double* const ds, double* const delTPhMoles, int kspec); +#endif + + //! Estimate equilibrium compositions + /*! + * Estimates equilibrium compositions. + * Algorithm covered in a section of Smith and Missen's Book. + * + * Linear programming module is based on using dbolm. + * + * @param aw aw[i[ Mole fraction work space (ne in length) + * @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length) + * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) + * @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length) + * @param test This is a small negative number. + */ + void vcs_inest(double* const aw, double* const sa, double* const sm, + double* const ss, double test); + + + //! Calculate the status of single species phases. + void vcs_SSPhase(void); + + //! This function recalculates the deltaG for reaction, irxn + /*! + * This function recalculates the deltaG for reaction irxn, + * given the mole numbers in molNum. It uses the temporary + * space mu_i, to hold the recalculated chemical potentials. + * It only recalculates the chemical potentials for species in phases + * which participate in the irxn reaction. + * + * Input + * ------------ + * @param irxn Reaction number + * @param molNum Current mole numbers of species to be used as + * input to the calculation (units = kmol) + * (length = totalNuMSpecies) + * + * Output + * ------------ + * @param ac output Activity coefficients (length = totalNumSpecies) + * Note this is only partially formed. Only species in + * phases that participate in the reaction will be updated + * @param mu_i diemsionless chemical potentials (length - totalNumSpecies + * Note this is only partially formed. Only species in + * phases that participate in the reaction will be updated + * + * @return Returns the dimensionless deltaG of the reaction + */ + double deltaG_Recalc_Rxn(const int stateCalc, + const size_t irxn, const double* const molNum, + double* const ac, double* const mu_i); + + //! Delete memory that isn't just resizeable STL containers + /*! + * This gets called by the destructor or by InitSizes(). + */ + void vcs_delete_memory(); + + //! Initialize the internal counters + /*! + * Initialize the internal counters containing the subroutine call + * values and times spent in the subroutines. + * + * ifunc = 0 Initialize only those counters appropriate for the top of + * vcs_solve_TP(). + * = 1 Initialize all counters. + */ + void vcs_counters_init(int ifunc); + + //! Create a report on the plog file containing timing and its information + /*! + * @param timing_print_lvl If 0, just report the iteration count. + * If larger than zero, report the timing information + */ + void vcs_TCounters_report(int timing_print_lvl = 1); + + void vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc); + + void vcs_setFlagsVolPhase(const size_t iph, const bool upToDate, const int stateCalc); + + //! Update all underlying vcs_VolPhase objects + /*! + * Update the mole numbers and the phase voltages of all phases in the + * vcs problem + * + * @param stateCalc Location of the update (either VCS_STATECALC_NEW or + * VCS_STATECALC_OLD). + */ + void vcs_updateMolNumVolPhases(const int stateCalc); + + +public: + //! value of the number of species used to malloc data structures + size_t NSPECIES0; + + //! value of the number of phases used to malloc data structures + size_t NPHASE0; + + //! Total number of species in the problems + size_t m_numSpeciesTot; + + //! Number of element constraints in the problem + /*! + * This is typically equal to the number of elements in the problem + */ + size_t m_numElemConstraints; + + //! Number of components calculated for the problem + size_t m_numComponents; + + //! Total number of non-component species in the problem + size_t m_numRxnTot; + + //! Current number of species in the problems + /*! + * Species can be deleted if they aren't + * stable under the current conditions + */ + size_t m_numSpeciesRdc; + + //! Current number of non-component species in the problem + /*! + * Species can be deleted if they aren't + * stable under the current conditions + */ + size_t m_numRxnRdc; + + //! Number of active species which are currently either treated as + //! minor species + size_t m_numRxnMinorZeroed; + + //! Number of Phases in the problem + size_t m_numPhases; + + //! Formula matrix for the problem + /*! + * FormulaMatrix[j][kspec] = Number of elements, j, in the kspec species + * + * Both element and species indecies are swapped. + */ + DoubleStarStar m_formulaMatrix; + + //! Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form. + /*! + * This is the stoichiometric coefficient matrix for the + * reaction which forms species kspec from the component species. A + * stoichiometric coefficient of one is assumed for the species kspec in this mechanism. + * + * NOTE: kspec = irxn + m_numComponents + * + * m_stoichCoeffRxnMatrix[irxn][j] : + * j refers to the component number, and irxn refers to the irxn_th non-component species. + * The stoichiometric coefficents multilpled by the Formula coefficients of the + * component species add up to the negative value of the number of elements in + * the species kspec. + * + * length = [nspecies0][nelements0] + */ + DoubleStarStar m_stoichCoeffRxnMatrix; + + //! Absolute size of the stoichiometric coefficients + /*! + * scSize[irxn] = abs(Size) of the stoichiometric + * coefficients. These are used to determine + * whether a given species should be + * handled by the alt_min treatment or + * should be handled as a major species. + */ + std::vector m_scSize; + + //! total size of the species + /*! + * This is used as a multiplier to the mole number in figuring out which + * species should be components. + */ + std::vector m_spSize; + + //! Standard state chemical potentials for species K at the current + //! temperature and pressure. + /*! + * The first NC entries are for components. The following NR entries are + * for the current non-component species in the mechanism. + */ + std::vector m_SSfeSpecies; + + //! Free energy vector from the start of the current iteration + /*! + * The free energies are saved at the start of the current iteration. + * Length = number of species + */ + std::vector m_feSpecies_old; + + //! Dimensionless new free energy for all the species in the mechanism + //! at the new tentatite T, P, and mole numbers. + /*! + * The first NC entries are for components. The following + * NR entries are for the current non-component species in the mechanism. + * Length = number of species + */ + std::vector m_feSpecies_new; + + //! Setting for whether to do an initial estimate + /*! + * Initial estimate: 0 Do not estimate the solution at all. Use the + * supplied mole numbers as is. + * 1 Only do an estimate if the element abundances + * aren't satisfied. + * -1 Force an estimate of the soln. Throw out the input + * mole numbers. + */ + int m_doEstimateEquil; + + //! Total moles of the species + /*! + * Total number of moles of the kth species. + * Length = Total number of species = m + */ + std::vector m_molNumSpecies_old; + + //! Specifies the species unknown type + /*! + * There are two types. One is the straightforward + * species, with the mole number w[k], as the + * unknown. The second is the an interfacial + * voltage where w[k] refers to the interfacial + * voltage in volts. + * These species types correspond to metalic + * electrons corresponding to electrodes. + * The voltage and other interfacial conditions + * sets up an interfacial current, which is + * set to zero in this initial treatment. + * Later we may have non-zero interfacial currents. + */ + std::vector m_speciesUnknownType; + + //! Change in the number of moles of phase, iphase, due to the noncomponent formation + //! reaction, irxn, for species, k: + /*! + * m_deltaMolNumPhase[irxn][iphase] = k = nc + irxn + */ + DoubleStarStar m_deltaMolNumPhase; + + //! This is 1 if the phase, iphase, participates in the formation reaction + //! irxn, and zero otherwise. PhaseParticipation[irxn][iphase] + IntStarStar m_phaseParticipation; + + //! electric potential of the iph phase + std::vector m_phasePhi; + + //! Tentative value of the mole number vector. It's also used to store the + //! mole fraction vector. + //std::vector wt; + std::vector m_molNumSpecies_new; + + //! Delta G(irxn) for the noncomponent species in the mechanism. + /*! + * Computed by the subroutine deltaG. m_deltaGRxn is the free + * energy change for the reaction which forms species K from the + * component species. This vector has length equal to the number + * of noncomponent species in the mechanism. It starts with + * the first current noncomponent species in the mechanism. + */ + std::vector m_deltaGRxn_new; + + //! Last deltag[irxn] from the previous step + std::vector m_deltaGRxn_old; + + //! Last deltag[irxn] from the previous step with additions for + //! possible births of zeroed phases. + std::vector m_deltaGRxn_Deficient; + + //! Temporary vector of Rxn DeltaG's + /*! + * This is used from time to time, for printing purposes + */ + std::vector m_deltaGRxn_tmp; + + //! Reaction Adjustments for each species during the current step + /*! + * delta Moles for each species during the current step. + * Length = number of species + */ + std::vector m_deltaMolNumSpecies; + + //! Element abundances vector + /*! + * Vector of moles of each element actually in the solution + * vector. Except for certain parts of the algorithm, + * this is a constant. + * Note other constraint conditions are added to this vector. + * This is input from the input file and + * is considered a constant from thereon. + * units = kmoles + */ + std::vector m_elemAbundances; + + //! Element abundances vector Goals + /*! + * Vector of moles of each element that are the goals of the + * simulation. This is a constant in the problem. + * Note other constraint conditions are added to this vector. + * This is input from the input file and + * is considered a constant from thereon. + * units = kmoles + */ + std::vector m_elemAbundancesGoal; + + //! Total number of kmoles in all phases + /*! + * This number includes the inerts. + * -> Don't use this except for scaling + * purposes + */ + double m_totalMolNum; + + //! Total kmols of species in each phase + /*! + * This contains the total number of moles of species in each phase + * + * Length = number of phases + */ + std::vector m_tPhaseMoles_old; + + //! total kmols of species in each phase in the tentative soln vector + /*! + * This contains the total number of moles of species in each phase + * in the tentative solution vector + * + * Length = number of phases + */ + std::vector m_tPhaseMoles_new; + + //! Temporary vector of length NPhase + mutable std::vector m_TmpPhase; + + //! Temporary vector of length NPhase + mutable std::vector m_TmpPhase2; + + //! Change in the total moles in each phase + /*! + * Length number of phases. + */ + std::vector m_deltaPhaseMoles; + + //! Temperature (Kelvin) + double m_temperature; + + //! Pressure (units are determined by m_VCS_UnitsFormat + /*! + * Values units + * -1: atm + * 0: atm + * 1: atm + * 2: atm + * 3: Pa + * Units being changed to Pa + */ + double m_pressurePA; + + //! Total kmoles of inert to add to each phase + /*! + * TPhInertMoles[iph] = Total kmoles of inert to add to each phase + * length = number of phases + */ + std::vector TPhInertMoles; + + //! Tolerance requirement for major species + double m_tolmaj; + + //! Tolerance requirements for minor species + double m_tolmin; + + //! Below this, major species aren't refined any more + double m_tolmaj2; + + //! Below this, minor species aren't refined any more + double m_tolmin2; + + //! Index vector that keeps track of the species vector rearrangement + /*! + * At the end of each run, the species vector and associated data gets put back + * in the original order. + * + * Example + * + * k = m_speciesMapIndex[kspec] + * + * kspec = current order in the vcs_solve object + * k = original order in the vcs_prob object and in the MultiPhase object + */ + std::vector m_speciesMapIndex; + + //! Index that keeps track of the index of the species within the local + //! phase + /*! + * This returns the local index of the species within the phase. Its argument + * is the global species index within the VCS problem. + * + * k = m_speciesLocalPhaseIndex[kspec] + * + * k varies between 0 and the nSpecies in the phase + * + * Length = number of species + */ + std::vector m_speciesLocalPhaseIndex; + + //! Index vector that keeps track of the rearrangement of the elements + /*! + * At the end of each run, the element vector and associated data gets put back + * in the original order. + * + * Example + * + * e = m_elementMapIndex[eNum] + * + * eNum = current order in the vcs_solve object + * e = original order in the vcs_prob object and in the MultiPhase object + */ + std::vector m_elementMapIndex; + + //! Mapping between the species index for noncomponent species and the + //! full species index. + /*! + * ir[irxn] = Mapping between the reaction index for + * noncomponent formation reaction of a species + * and the full species + * index. + * - Initially set to a value of K = NC + I + * This vector has length equal to number + * of noncomponent species in the mechanism. + * It starts with the first current + * noncomponent species in the mechanism. + * kspec = ir[irxn] + */ + std::vector m_indexRxnToSpecies; + + //! Major -Minor status vector for the species in the problem + /*! + * The index for this is species. The reaction that this is referring + * to is + * kspec = irxn + m_numComponents + * + * kspec : 2 -> Component species VCS_SPECIES_COMPONENT + * -> deprecated, want to assign -2 to some + * component species. We can already determine + * whether the species is a component from + * its position in the species vector. + * 1 -> Major species VCS_SPECIES_MAJOR + * 0 -> Minor species VCS_SPECIES_MINOR + * -1 -> The species lies in a multicomponent phase + * that exists. Its concentration is currently + * very low, necessitating a different method + * of calculation. + * - VCS_SPECIES_ZEROEDPHASE + * -2 -> The species lies in a multicomponent phase + * which currently doesn't exist. + * Its concentration is currently zero. + * - VCS_SPECIES_ZEROEDMS + * -3 -> Species lies in a single-species phase which + * is currently zereod out. + * - VCS_SPECIES_ZEREODSS + * -4 -> Species has such a small mole fraction it is + * deleted even though its phase may possibly exist. + * The species is believed to have such a small + * mole fraction that it best to throw the + * calculation of it out. It will be added back in + * at the end of the calculation. + * - VCS_SPECIES_DELETED + * -5 -> Species refers to an electron in the metal + * The unknown is equal to the interfacial voltage + * drop across the interface on the SHE (standard + * hydroogen electrode) scale (volts). + * - VCS_SPECIES_INTERFACIALVOLTAGE + * -6 -> Species lies in a multicomponent phase that + * is zeroed atm and will stay deleted due to a + * choice from a higher level. + * These species will formally always have zero + * mole numbers in the solution vector. + * - VCS_SPECIES_ZEROEDPHASE + * -7 -> The species lies in a multicomponent phase which + * currently does exist. Its concentration is currently + * identically zero, though the phase exists. Note, this + * is a temporary condition that exists at the start + * of an equilibrium problem. + * The species is soon "birthed" or "deleted". + * - VCS_SPECIES_ACTIVEBUTZERO + * -8 -> The species lies in a multicomponent phase which + * currently does exist. Its concentration is currently + * identically zero, though the phase exists. This is + * a permanent condition due to stoich constraints + * - VCS_SPECIES_STOICHZERO + * + */ + std::vector m_speciesStatus; + + //! Mapping from the species number to the phase number + std::vector m_phaseID; + + //! Boolean indicating whether a species belongs to a single-species phase + // vector can't be used here because it doesn't work with std::swap + std::vector m_SSPhase; + + //! Species string name for the kth species + /*! + * Species string name for the kth species + */ + std::vector m_speciesName; + + //! Vector of strings containing the element names + /*! + * ElName[j] = String containing element names + */ + std::vector m_elementName; + + //! Type of the element constraint + /*! + * m_elType[j] = type of the element + * 0 VCS_ELEM_TYPE_ABSPOS Normal element that is positive + * or zero in all species. + * 1 VCS_ELEM_TPYE_ELECTRONCHARGE element dof that corresponds + * to the electronic charge DOF. + * 2 VCS_ELEM_TYPE_CHARGENEUTRALITY element dof that + * corresponds to a required charge + * neutrality constraint on the phase. + * The element abundance is always exactly zero. + * 3 VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may + * mean that a species has neg 0 or pos value + * of that constraint (other than charge) + */ + std::vector m_elType; + + //! Specifies whether an element constraint is active + /*! + * The default is true + * Length = nelements + */ + std::vector m_elementActive; + + //! Array of Phase Structures + /*! + * Length = number of phases + */ + std::vector m_VolPhaseList; + + //! String containing the title of the run + std::string m_title; + + //! This specifies the current state of units for the Gibbs free energy + //! properties in the program. + /*! + *. The default is to have this unitless + */ + char m_unitsState; + + //! Multiplier for the mole numbers within the nondimensionless formulation + /*! + * All numbers within the main routine are on an absolute basis. This + * presents some problems wrt very large and very small mole numbers. + * We get around this by using a multiplier coming into and coming + * out of the equilibrium routines + */ + double m_totalMoleScale; + + //! specifies the activity convention of the phase containing the species + /*! + * SpecActConvention[kspec] + * 0 = molar based + * 1 = molality based + * length = number of species + */ + std::vector m_actConventionSpecies; + + //! specifies the activity convention of the phase. + /*! + * 0 = molar based + * 1 = molality based + * length = number of phases + */ + std::vector m_phaseActConvention; + + //! specifies the ln(Mnaught) used to calculate the chemical potentials + /*! + * For molar based activity conventions + * this will be equal to 0.0 + * length = number of species + */ + std::vector m_lnMnaughtSpecies; + + //! Molar-based Activity Coefficients for Species + /*! + * + * Length = number of species + */ + std::vector m_actCoeffSpecies_new; + + //! Molar-based Activity Coefficients for Species based on old mole numbers + /*! + * These activity coefficients are based on the m_molNumSpecies_old values + * Molar based activity coeffients. + * Length = number of species + */ + std::vector m_actCoeffSpecies_old; + + //! Change in activity coefficient with mole number + /*! + * length = [nspecies][nspecies] + * + * (This is a temporary array that + * gets regenerated every time it's + * needed. It is not swapped wrt species + * (unused atm) + */ + DoubleStarStar m_dLnActCoeffdMolNum; + + //! Molecular weight of each species + /*! + * units = kg/kmol + * length = number of species + * + * note: this is a candidate for removal. I don't think we use it. + */ + std::vector m_wtSpecies; + + //! Charge of each species + /*! + * Length = number of species + */ + std::vector m_chargeSpecies; + + std::vector > phasePopProblemLists_; + + //! Vector of pointers to thermostructures which identify the model + //! and parameters for evaluating the thermodynamic functions for that + //! particular species. + /*! + * SpeciesThermo[k] pointer to the thermo information for the kth species + */ + std::vector m_speciesThermoList; + + //! Choice of Hessians + /*! + * If this is true, then we will use a better approximation to the + * Hessian based on Jacobian of the ln(ActCoeff) with respect to mole + * numbers + */ + int m_useActCoeffJac; + + //! Total volume of all phases + /*! + * units are m^3 + */ + double m_totalVol; + + //! Partial molar volumes of the species + /*! + * units = mks (m^3/kmol) -determined by m_VCS_UnitsFormat + * Length = number of species + */ + std::vector m_PMVolumeSpecies; + + //! dimensionless value of Faraday's constant + /*! + * F / RT (1/volt) + */ + double m_Faraday_dim; + + //! Timing and iteration counters for the vcs object + VCS_COUNTERS* m_VCount; + + + //! Debug printing lvl + /*! + * Levels correspond to the following guidlines + * - 0 No printing at all + * - 1 Serious warnings or fatal errors get one line + * - 2 one line per eacdh successful vcs package call + * - 3 one line per every successful solve_TP calculation + * - 4 one line for every successful operation -> solve_TP gets a summary report + * - 5 each iteration in solve_TP gets a report with one line per species + * - 6 Each decision in solve_TP gets a line per species in addition to 4 + * - 10 Additionally Hessian matrix is printed out + * + * Levels of printing above 4 are only accessible when DEBUG_MODE is turned on + */ + int m_debug_print_lvl; + + //! printing level of timing information + /*! + * 1 allowing printing of timing + * 0 do not allow printing of timing -> everything is printed + * as a NA. + */ + int m_timing_print_lvl; + + //! Units for the chemical potential data: + /*! + * VCS_UnitsFormat = Units for the chemical potential data: + * -1: kcal/mol + * 0: MU/RT + * 1: kJ/mol + * 2: Kelvin + * 3: J / kmol + * and pressure data: + * -1: Pa + * 0: Pa + * 1: Pa + * 2: pa + * 3: Pa + */ + int m_VCS_UnitsFormat; + + friend class vcs_phaseStabilitySolve; + +}; + +#ifdef ALTLINPROG +#else +int linprogmax(double*, double*, double*, double*, size_t, size_t, size_t); +#endif + +} +#endif + diff --git a/src/numerics/BEulerInt.h b/include/cantera/numerics/BEulerInt.h similarity index 100% rename from src/numerics/BEulerInt.h rename to include/cantera/numerics/BEulerInt.h