From 458fbf579da0e67ef4df84036a833634a817c40f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 5 Apr 2012 19:17:01 +0000 Subject: [PATCH] removed include files that were suppose to have been moved. --- src/equil/vcs_DoubleStarStar.h | 127 - src/equil/vcs_IntStarStar.cpp | 4 +- src/equil/vcs_IntStarStar.h | 101 - src/equil/vcs_VolPhase.h | 1038 -------- src/equil/vcs_internal.h | 528 ----- src/equil/vcs_prob.h | 372 --- src/equil/vcs_solve.h | 2083 ----------------- .../VCSnonideal/NaCl_equil/nacl_equil.cpp | 2 +- 8 files changed, 3 insertions(+), 4252 deletions(-) delete mode 100644 src/equil/vcs_DoubleStarStar.h delete mode 100644 src/equil/vcs_IntStarStar.h delete mode 100644 src/equil/vcs_VolPhase.h delete mode 100644 src/equil/vcs_internal.h delete mode 100644 src/equil/vcs_prob.h delete mode 100644 src/equil/vcs_solve.h diff --git a/src/equil/vcs_DoubleStarStar.h b/src/equil/vcs_DoubleStarStar.h deleted file mode 100644 index b15630681..000000000 --- a/src/equil/vcs_DoubleStarStar.h +++ /dev/null @@ -1,127 +0,0 @@ -/** - * @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/src/equil/vcs_IntStarStar.cpp b/src/equil/vcs_IntStarStar.cpp index 626561dd0..76cfef9e7 100644 --- a/src/equil/vcs_IntStarStar.cpp +++ b/src/equil/vcs_IntStarStar.cpp @@ -3,12 +3,12 @@ * * Header file for class IntStarStar */ -#include "vcs_IntStarStar.h" +#include "cantera/equil/vcs_IntStarStar.h" namespace VCSnonideal { -//!Default constructor. Create an empty array. +//Default constructor. Create an empty array. IntStarStar::IntStarStar() : m_nrows(0), m_ncols(0) diff --git a/src/equil/vcs_IntStarStar.h b/src/equil/vcs_IntStarStar.h deleted file mode 100644 index 0f328b4b7..000000000 --- a/src/equil/vcs_IntStarStar.h +++ /dev/null @@ -1,101 +0,0 @@ -/** - * @file vcs_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 - IntStarStar(const IntStarStar& y); - - /// assignment operator - 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 - std::vector m_data; - - std::vector m_colAddr; - - //! number of rows - size_t m_nrows; - - //! number of columns - size_t m_ncols; -}; - -} - -#endif - - diff --git a/src/equil/vcs_VolPhase.h b/src/equil/vcs_VolPhase.h deleted file mode 100644 index 3013b3b1b..000000000 --- a/src/equil/vcs_VolPhase.h +++ /dev/null @@ -1,1038 +0,0 @@ -/** - * @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 "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 maximizes 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 up-to-date mole number vector - //! and potential with respect to the current vcs state calc status - bool m_UpToDate; - - //! Boolean indicating whether activity coefficients are up to date. - /*! - * 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 up to date. - /*! - * 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 up to date. - /*! - * 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 up to date. - /*! - * GStar is sensitive to the temperature and the pressure, only - */ - mutable bool m_UpToDate_GStar; - - - //! Boolean indicating whether G0 is up to date. - /*! - * 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/src/equil/vcs_internal.h b/src/equil/vcs_internal.h deleted file mode 100644 index 7d9da7f04..000000000 --- a/src/equil/vcs_internal.h +++ /dev/null @@ -1,528 +0,0 @@ -/** - * @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 "vcs_DoubleStarStar.h" -#include "vcs_Exception.h" - -#include "cantera/base/global.h" - -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 singular - * - 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 singular - * - 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 parameter - * @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 multiple occurrences 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/src/equil/vcs_prob.h b/src/equil/vcs_prob.h deleted file mode 100644 index b2fe90fe7..000000000 --- a/src/equil/vcs_prob.h +++ /dev/null @@ -1,372 +0,0 @@ -/** - * @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 - * quantities 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/src/equil/vcs_solve.h b/src/equil/vcs_solve.h deleted file mode 100644 index d5242fba1..000000000 --- a/src/equil/vcs_solve.h +++ /dev/null @@ -1,2083 +0,0 @@ -/** - * @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 "vcs_DoubleStarStar.h" -#include "vcs_IntStarStar.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 coefficients multiplied 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 tentative 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 metallic - * 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/test_problems/VCSnonideal/NaCl_equil/nacl_equil.cpp b/test_problems/VCSnonideal/NaCl_equil/nacl_equil.cpp index a0a84e0a4..204d1eadf 100644 --- a/test_problems/VCSnonideal/NaCl_equil/nacl_equil.cpp +++ b/test_problems/VCSnonideal/NaCl_equil/nacl_equil.cpp @@ -4,7 +4,7 @@ #include "cantera/equilibrium.h" #include "cantera/equil/vcs_MultiPhaseEquil.h" -#include "equil/vcs_internal.h" +#include "cantera/equil/vcs_internal.h" #include "cantera/thermo/ThermoFactory.h" #include "cantera/thermo/IdealGasPhase.h"