diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 916f0831b..9e7008a08 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -23,7 +23,7 @@ namespace VCSnonideal { - /**************************************************************************** + /* * * vcs_VolPhase(): * @@ -68,6 +68,7 @@ namespace VCSnonideal { { m_owningSolverObject = owningSolverObject; } + /************************************************************************************/ /* * @@ -82,6 +83,7 @@ namespace VCSnonideal { sp = 0; } } + /************************************************************************************/ /* * @@ -131,8 +133,9 @@ namespace VCSnonideal { */ *this = b; } - - /***************************************************************************** + /*****************************************************************************/ + + /* * Assignment operator() * * (note, this is used, so keep it current!) @@ -239,7 +242,7 @@ namespace VCSnonideal { } return *this; } - + /************************************************************************************/ void vcs_VolPhase::resize(int phaseNum, int nspecies, const char *phaseName, double molesInert) { @@ -327,7 +330,7 @@ namespace VCSnonideal { m_UpToDate_VolPM = false; m_UpToDate_GStar = false; } - + /*******************************************************************************/ //! Evaluate activity coefficients /*! @@ -357,8 +360,9 @@ namespace VCSnonideal { } m_UpToDate_AC = true; } + /********************************************************************************/ - /****************************************************************************** + /* * * Evaluate one activity coefficients. * @@ -369,7 +373,7 @@ namespace VCSnonideal { evaluateActCoeff(); return(ActCoeff[kspec]); } - + /************************************************************************************/ // Gibbs free energy calculation at a temperature for the reference state // of each species @@ -424,8 +428,7 @@ namespace VCSnonideal { * @param TKelvin Current temperature * @param pres Current pressure (pascal) */ - void vcs_VolPhase::GStar_calc() { - setState_TP(Temp, Pres); + void vcs_VolPhase::GStar_calc() const { if (!m_UpToDate_GStar) { if (m_useCanteraCalls) { TP_ptr->getStandardChemPotentials(VCS_DATA_PTR(StarChemicalPotential)); @@ -452,16 +455,14 @@ namespace VCSnonideal { * The kth species standard state G is returned * * @param kspec Species number (within the phase) - * @param TKelvin Current temperature (kelvin) - * @param pres Current pressure (pascal) * * @return Gstar[kspec] returns the gibbs free energy for the - * standard state of the kth species. + * standard state of the kspec species. */ - double vcs_VolPhase::GStar_calc_one(int kspec, double tkelvin, - double pres) { - setState_TP(tkelvin, pres); - GStar_calc(); + double vcs_VolPhase::GStar_calc_one(int kspec) { + if (!m_UpToDate_GStar) { + GStar_calc(); + } return StarChemicalPotential[kspec]; } /***********************************************************************/ @@ -610,9 +611,7 @@ namespace VCSnonideal { } if (!m_isIdealSoln) { m_UpToDate_AC = false; - m_UpToDate_VolStar = false; m_UpToDate_VolPM = false; - m_UpToDate_GStar = false; } } @@ -713,7 +712,6 @@ namespace VCSnonideal { */ void vcs_VolPhase::sendToVCS_GStar(double * const gstar){ if (!m_UpToDate_GStar) { - setState_TP(Temp, Pres); GStar_calc(); } int kglob; @@ -771,7 +769,7 @@ namespace VCSnonideal { m_UpToDate_VolPM = false; m_UpToDate_GStar = false; } - /***********************************************************************/ + /***********************************************************************/ // Molar volume calculation for standard states /* @@ -783,8 +781,7 @@ namespace VCSnonideal { * * Calculations are in m**3/kmol */ - void vcs_VolPhase::VolStar_calc(double tkelvin, double pres) { - setState_TP(tkelvin, pres); + void vcs_VolPhase::VolStar_calc() const { if (!m_UpToDate_VolStar) { if (m_useCanteraCalls) { TP_ptr->getStandardVolumes(VCS_DATA_PTR(StarMolarVol)); @@ -794,13 +791,13 @@ namespace VCSnonideal { vcs_SpeciesProperties *sProp = ListSpeciesPtr[k]; VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo; StarMolarVol[k] = - (sTherm->VolStar_calc(kglob, tkelvin, pres)); + (sTherm->VolStar_calc(kglob, Temp, Pres)); } } m_UpToDate_VolStar = true; } } - /***********************************************************************/ + /***********************************************************************/ // Update the moles within the phase, if necessary /* @@ -815,7 +812,7 @@ namespace VCSnonideal { * */ void vcs_VolPhase::updateFromVCS_MoleNumbers(const int stateCalc) { - if (!m_UpToDate) { + if (!m_UpToDate || (stateCalc != m_vcsStateStatus)) { if (stateCalc == VCS_STATECALC_OLD || stateCalc == VCS_STATECALC_NEW) { if (m_owningSolverObject) { setMolesFromVCS(stateCalc); @@ -839,13 +836,16 @@ namespace VCSnonideal { * state */ double vcs_VolPhase::VolStar_calc_one(int kspec, double tkelvin, - double pres) - { - VolStar_calc(tkelvin, pres); + double pres) { + setState_TP(tkelvin, pres); + if (!m_UpToDate_VolStar) { + VolStar_calc(); + } return StarMolarVol[kspec]; } + /****************************************************************************/ - /**************************************************************************** + /* * * VolPM_calc */ @@ -885,20 +885,16 @@ namespace VCSnonideal { m_UpToDate_VolPM = true; return m_totalVol; } + /************************************************************************************/ /* * updateLnActCoeffJac(): * */ - void vcs_VolPhase::updateLnActCoeffJac(const double * const moleNumbersVCS) { + void vcs_VolPhase::updateLnActCoeffJac() { int k, j; double deltaMoles_j = 0.0; - /* - * Make sure the base state of this object is fully up to date. - * with the current values of the mole numbers. - * -> This sets TMoles and Xmol[] - */ - setMolesFromVCS(VCS_STATECALC_OLD, moleNumbersVCS); + /* * Evaluate the current base activity coefficients. @@ -961,6 +957,7 @@ namespace VCSnonideal { _updateMoleFractionDependencies(); evaluateActCoeff(); } + /************************************************************************************/ // Downloads the ln ActCoeff jacobian into the VCS version of the // ln ActCoeff jacobian. @@ -986,7 +983,7 @@ namespace VCSnonideal { } } } - /**********************************************************************/ + /************************************************************************************/ // Set the pointer for Cantera's ThermoPhase parameter /* @@ -1043,6 +1040,7 @@ namespace VCSnonideal { m_useCanteraCalls = false; } } + /************************************************************************************/ // Return a const ThermoPhase pointer corresponding to this phase /* @@ -1051,18 +1049,22 @@ namespace VCSnonideal { const Cantera::ThermoPhase *vcs_VolPhase::ptrThermoPhase() const { return TP_ptr; } + /************************************************************************************/ double vcs_VolPhase::TotalMoles() const { return TMoles; } + /************************************************************************************/ double vcs_VolPhase::molefraction(int k) const { return Xmol[k]; } + /************************************************************************************/ void vcs_VolPhase::setTotalMoles(double tmols) { TMoles = tmols; } + /************************************************************************************/ // Return a string representing the equation of state /* diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index b221d29d6..f36970b45 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -28,30 +28,30 @@ namespace Cantera { } namespace VCSnonideal { -/* - * Models for the species activity coefficients - * - */ + /* + * 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_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 - */ + /* + * + * 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 - */ + /* + * 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 @@ -67,742 +67,756 @@ namespace VCSnonideal { 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 + //! Phase information and Phase calculations for vcs. /*! - * @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. - */ - void resize(int phaseNum, int numSpecies, const char *phaseName, - double molesInert = 0.0); - -private: - //! Evaluate activity coefficients - /*! - * 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 evaluateActCoeff() const; -public: - //! 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. + * 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. * - * @param kspec species number - */ - double AC_calc_one(int kspec) const; - - - //! 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. + * 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. * - * @param molesSpeciesVCS array of mole numbers. Note, the indecises 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 * const 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. + * This object is considered not to own the underlying Cantera ThermoPhase + * object for the phase. * - * - * @param molesSpeciesVCS array of mole numbers. Note, the indecises - * 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. + * 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. * - */ - void setMolesFromVCSCheck(const int stateCalc, - const double * molesSpeciesVCS, - const double * const TPhMoles, - int iphase = -1); - - //! 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(double * const AC) const; - - //! set the electric potential of the phase - /*! - * @param phi electric potential (volts) - */ - void setElectricPotential(double phi); - - //! Returns the electric field of the phase - /*! - * Units are potential - */ - double electricPotential() 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. + * 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. * - * @param TKelvin Current temperature - * @param pres Current pressure - */ - void GStar_calc(); - - //! 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 + * 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. * - * @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(int kspec, double TKelvin, double pres); - - //! Gibbs free energy calculation at a temperature for the reference state - //! of each species - /*! - * @param TKelvin temperature - */ - void G0_calc(double TKelvin); - - //! 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(int kspec, double TKelvin); - - //! 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); - - //! 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 VolStar_calc(double TKelvin, double pres); - - //! 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(int kglob, double TKelvin, double pres); - -private: - //! Calculate the partial molar volumes of all species and return the - //! total volume - /*! - * Calculates these quantitites internally - * - * @return total volume - */ - double VolPM_calc() const; -public: - //! 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); - - //! 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(double temperature_Kelvin, double pressure_PA); - - //! 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 - * - * @param moleNumbers Mole numbers are input. - */ - void updateLnActCoeffJac(const double * const moleNumbers); - - // 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] + * 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. * - * j = id of the species mole number - * k = id of the species activity coefficient - */ - void sendToVCS_LnActCoeffJac(double * const * const LnACJac_VCS) const; - - //! 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 - /*! - * Returns the mole fraction of the kspec species - * - */ - double molefraction(int kspec) const; - - //! Sets the total moles in the phase - /*! - * - */ - void setTotalMoles(double tmols); - - //! 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); - - //! Return a const reference to the mole fractions - const std::vector & moleFractions() const; - - //! Returns whether the phase is an ideal solution phase - bool isIdealSoln() const; - - //! Returns whether the object is using cantera calls. - bool usingCanteraCalls() const; - -private: - - //! Updates the mole fractions in subobjects - /*! - * Whenever the mole fractions change, this routine - * should be called. - */ - void _updateMoleFractionDependencies(); - - /******************************************************************************/ -public: - //! 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; - - //! 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. - */ - int 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. - */ - int Domain_ID; - - //! If true, this phase consists of a single species - int 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 EqnState; - - //! Number of element constraints within the problem - /*! - * This is usually equal to the number of elements. - * - */ - int nElemConstraints; - - //! 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 - */ - int ChargeNeutralityElement; - - //! vector of strings containing the element names - /*! - * Length = nElemConstraints - */ - std::vector ElName; - - //! boolean indicating whether element constraint is active - //! for the current problem - std::vector ElActive; - - //! Type of the element - /*! - * 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_elType; - - //! Formula Matrix for the phase - /*! - * FormulaMatrix[j][kspec] - * = Formula Matrix for the species - * Number of elements, j, - * in the kspec species - */ - DoubleStarStar 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 SpeciesUnknownType; - - //! Index of the element number in the global list of elements - //! storred in VCS_PROB or VCS_SOLVE - std::vector ElGlobalIndex; - - //! Number of species in the phase - int NVolSpecies; - - //! String name for the phase - std::string PhaseName; - - //! Total moles of inert in the phase - double TMolesInert; - - //! molar volume of the inert species if present - /*! - * units m**3 / kg - */ - double m_molarVolInert; - - //! 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 ActivityConvention; - - //! Boolean indicating whether the phase is an ideal solution - //! and therefore it's molar-based activity coefficients are - //! uniformly equal to one. - bool m_isIdealSoln; - - //! Current state of existence: - /*! - * 0 : Doesn't exist currently - * 1 : Does exist currently - * 2 : Always exists because it contains - * inerts which can't exist in any other - * phase - */ - int Existence; - - //! Index of the species which is special in - //! with respect to the thermo treatment. - /*! - * For water models this index will point to - * the index for water. - * defaults to 0 - */ - int IndexSpecialSpecies; - - //! Integer representing the activity coefficient model - /*! - * The known models are listed at the top of this page - */ - int Activity_Coeff_Model; - - //! General pointer for hanging stuff off of - /*! - * Currently, not implemented very well - */ - void *Activity_Coeff_Params; - - //! 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; - - //! Boolean indicating whether IndSpecies is contiguous - bool IndSpeciesContig; - - //! 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; - - //! Units for the chemical potential data, pressure data, volume, - //! and species amounts - /*! - * All internally storred 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. - */ - int m_VCS_UnitsFormat; - -private: - //! 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; - - /** - * Variables Having to do with Calculated States - */ - - //! Total mols in the phase - /*! - * units vary - */ - double TMoles; - - //! Vector of the current mole fractions for species - //! in the phase - std::vector Xmol; - -public: - - //! 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 - int 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 - * - */ - int m_vcsStateStatus; - -private: - - //! Value of the potential for the phase (Volts) - double m_phi; - - //! Boolean indicating whether the object has an uptodate mole number vector - //! and potential with respect to the current vcs state calc status - bool m_UpToDate; - - //! Boolean indicating whether activity coefficients are uptodate. - /*! * Activity coefficients and volume calculations are lagged. They are only * called when they are needed (and when the state has changed so that they * need to be recalculated). */ - mutable bool m_UpToDate_AC; + class vcs_VolPhase { + public: - //! Boolean indicating whether Star volumes are uptodate. + /************************************************************************* + * 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(int phaseNum, int numSpecies, const char *phaseName, + double molesInert = 0.0); + + private: + //! Evaluate activity coefficients + /*! + * 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 evaluateActCoeff() const; + + public: + //! 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(int kspec) const; + + //! 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 indecises 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 * const 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 indecises + * 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, + int iphase = -1); + + //! 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(double * const AC) const; + + //! set the electric potential of the phase + /*! + * @param phi electric potential (volts) + */ + void setElectricPotential(double phi); + + //! Returns the electric field of the phase + /*! + * Units are potential + */ + double electricPotential() const; + + private: + //! 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 GStar_calc() const; + + public: + //! 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(int kspec); + + //! Gibbs free energy calculation at a temperature for the reference state + //! of each species + /*! + * @param TKelvin temperature + */ + void G0_calc(double TKelvin); + + //! 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(int kspec, double TKelvin); + + //! 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); + + private: + //! 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 VolStar_calc() const; + + public: + //! 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(int kglob, double TKelvin, double pres); + + private: + //! Calculate the partial molar volumes of all species and return the + //! total volume + /*! + * Calculates these quantitites internally + * + * @return total volume + */ + double VolPM_calc() const; + + public: + //! 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); + + //! 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(double temperature_Kelvin, double pressure_PA); + + //! 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(); + + // 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) const; + + //! 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(int kspec) const; + + //! Sets the total moles in the phase + /*! + * @param totalMols Total moles in the phase (kmol) + */ + void setTotalMoles(double totalMols); + + //! 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); + + //! Return a const reference to the mole fractions + const std::vector & moleFractions() const; + + //! Returns whether the phase is an ideal solution phase + bool isIdealSoln() const; + + //! Returns whether the object is using cantera calls. + bool usingCanteraCalls() const; + + private: + + //! Updates the mole fraction depenpencies + /*! + * Whenever the mole fractions change, this routine + * should be called. + */ + void _updateMoleFractionDependencies(); + + /************************************************************************* + * MEMBER DATA * + ************************************************************************/ + + public: + //! 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; + + //! 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. + */ + int 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. + */ + int Domain_ID; + + //! If true, this phase consists of a single species + int 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 EqnState; + + //! Number of element constraints within the problem + /*! + * This is usually equal to the number of elements. + */ + int nElemConstraints; + + //! 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 + */ + int ChargeNeutralityElement; + + //! vector of strings containing the element names + /*! + * Length = nElemConstraints + */ + std::vector ElName; + + //! boolean indicating whether element constraint is active + //! for the current problem + std::vector ElActive; + + //! Type of the element + /*! + * 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_elType; + + //! Formula Matrix for the phase + /*! + * FormulaMatrix[j][kspec] + * = Formula Matrix for the species + * Number of elements, j, + * in the kspec species + */ + DoubleStarStar 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 SpeciesUnknownType; + + //! Index of the element number in the global list of elements + //! storred in VCS_PROB or VCS_SOLVE + std::vector ElGlobalIndex; + + //! Number of species in the phase + int NVolSpecies; + + //! String name for the phase + std::string PhaseName; + + //! Total moles of inert in the phase + double TMolesInert; + + //! molar volume of the inert species if present + /*! + * units m**3 / kg + */ + double m_molarVolInert; + + //! 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 ActivityConvention; + + //! Boolean indicating whether the phase is an ideal solution + //! and therefore it's molar-based activity coefficients are + //! uniformly equal to one. + bool m_isIdealSoln; + + //! Current state of existence: + /*! + * 0 : Doesn't exist currently + * 1 : Does exist currently + * 2 : Always exists because it contains + * inerts which can't exist in any other + * phase + */ + int Existence; + + //! Index of the species which is special in + //! with respect to the thermo treatment. + /*! + * For water models this index will point to + * the index for water. + * defaults to 0 + */ + int IndexSpecialSpecies; + + //! Integer representing the activity coefficient model + /*! + * The known models are listed at the top of this page + */ + int Activity_Coeff_Model; + + //! General pointer for hanging stuff off of + /*! + * Currently, not implemented very well + */ + void *Activity_Coeff_Params; + + //! 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; + + //! Boolean indicating whether IndSpecies is contiguous + bool IndSpeciesContig; + + //! 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; + + //! Units for the chemical potential data, pressure data, volume, + //! and species amounts + /*! + * All internally storred 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. + */ + int m_VCS_UnitsFormat; + + private: + //! 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 TMoles; + + //! Vector of the current mole fractions for species + //! in the phase + std::vector Xmol; + + public: + + //! 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 + int 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; + + private: + //! 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; + + public: + //! 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 + */ + int m_vcsStateStatus; + + private: + + //! Value of the potential for the phase (Volts) + double m_phi; + + public: + //! Boolean indicating whether the object has an uptodate mole number vector + //! and potential with respect to the current vcs state calc status + bool m_UpToDate; + + private: + //! Boolean indicating whether activity coefficients are uptodate. + /*! + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + */ + mutable bool m_UpToDate_AC; + + //! Boolean indicating whether Star volumes are uptodate. + /*! + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + * Star volumes are sensitive to temperature and pressure + */ + mutable bool m_UpToDate_VolStar; + + //! Boolean indicating whether partial molar volumes are uptodate. + /*! + * Activity coefficients and volume calculations are lagged. They are only + * called when they are needed (and when the state has changed so that they + * need to be recalculated). + * partial molar volumes are sensitive to everything + */ + mutable bool m_UpToDate_VolPM; + + //! Boolean indicating whether GStar is uptodate. + /*! + * GStar is sensitive to the temperature and the pressure, only + */ + mutable bool m_UpToDate_GStar; + + //! 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; + + public: + //! Reference pressure for the phase + double RefPres; + + }; + + //! Return a string representing the equation of state /*! - * 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 + * @param EOSType : integer value of the equation of state + * + * @return returns a string representing the EOS */ - mutable bool m_UpToDate_VolStar; - - //! Boolean indicating whether partial molar volumes are uptodate. - /*! - * Activity coefficients and volume calculations are lagged. They are only - * called when they are needed (and when the state has changed so that they - * need to be recalculated). - * partial molar volumes are sensitive to everything - */ - mutable bool m_UpToDate_VolPM; - - //! Boolean indicating whether GStar is uptodate. - /*! - * GStar is sensitive to the temperature and the pressure, only - */ - mutable bool m_UpToDate_GStar; - - //! 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; - -public: - - //! Reference pressure for the phase - double RefPres; - -}; - -//! 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); + std::string string16_EOSType(int EOSType); } diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index 8a81cd7f4..e0a69458d 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -27,13 +27,13 @@ using namespace std; namespace VCSnonideal { -/* - * VCS_PROB: constructor - * - * We initialize the arrays in the structure to the appropriate sizes. - * And, we initialize all of the elements of the arrays to defaults. - */ -VCS_PROB::VCS_PROB(int nsp, int nel, int nph) : + /* + * VCS_PROB: constructor + * + * We initialize the arrays in the structure to the appropriate sizes. + * And, we initialize all of the elements of the arrays to defaults. + */ + VCS_PROB::VCS_PROB(int nsp, int nel, int nph) : prob_type(VCS_PROBTYPE_TP), nspecies(nsp), NSPECIES0(0), @@ -55,361 +55,360 @@ VCS_PROB::VCS_PROB(int nsp, int nel, int nph) : m_NumBasisOptimizations(0), m_printLvl(0), vcs_debug_print_lvl(0) -{ - NSPECIES0 = nspecies; - if (nspecies <= 0) { - plogf("number of species is zero or neg\n"); - exit(-1); - } - NE0 = ne; - if (ne <= 0) { - plogf("number of elements is zero or neg\n"); - exit(-1); - } - NPHASE0 = NPhase; - if (NPhase <= 0) { - plogf("number of phases is zero or neg\n"); - exit(-1); - } - if (nspecies < NPhase) { - plogf("number of species is less than number of phases\n"); - exit(-1); - } - - m_gibbsSpecies.resize(nspecies, 0.0); - w.resize(nspecies, 0.0); - mf.resize(nspecies, 0.0); - gai.resize(ne, 0.0); - FormulaMatrix.resize(ne, nspecies, 0.0); - SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); - VolPM.resize(nspecies, 0.0); - PhaseID.resize(nspecies, -1); - SpName.resize(nspecies, ""); - ElName.resize(ne, ""); - m_elType.resize(ne, VCS_ELEM_TYPE_ABSPOS); - ElActive.resize(ne, 1); - WtSpecies.resize(nspecies, 0.0); - Charge.resize(nspecies, 0.0); - SpeciesThermo.resize(nspecies,0); - for (int kspec = 0; kspec < nspecies; kspec++) { - VCS_SPECIES_THERMO *ts_tmp = new VCS_SPECIES_THERMO(0, 0); - if (ts_tmp == 0) { - plogf("Failed to init a ts struct\n"); + { + NSPECIES0 = nspecies; + if (nspecies <= 0) { + plogf("number of species is zero or neg\n"); exit(-1); } - SpeciesThermo[kspec] = ts_tmp; - } - VPhaseList.resize(nph, 0); - for (int iphase = 0; iphase < NPhase; iphase++) { - VPhaseList[iphase] = new vcs_VolPhase(); - } -} -/**************************************************************************/ -/**************************************************************************/ -/**************************************************************************/ -/* - * VCS_PROB_INPUT:destructor - * - * We need to manually free all of the arrays. - */ -VCS_PROB::~VCS_PROB() { - for (int i = 0; i < nspecies; i++) { - delete SpeciesThermo[i]; - SpeciesThermo[i] = 0; - } - for (int iph = 0; iph < NPhase; iph++) { - delete VPhaseList[iph]; - VPhaseList[iph] = 0; - } -} - -// 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 VCS_PROB::resizePhase(int nPhase, int force) { - if (force || nPhase > NPHASE0) { - NPHASE0 = nPhase; - } -} - -// 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 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 VCS_PROB::resizeSpecies(int nsp, int force) { - if (force || nsp > NSPECIES0) { - m_gibbsSpecies.resize(nsp, 0.0); - w.resize(nsp, 0.0); - mf.resize(nsp, 0.0); - FormulaMatrix.resize(NE0, nsp, 0.0); - SpeciesUnknownType.resize(nsp, VCS_SPECIES_TYPE_MOLNUM); - VolPM.resize(nsp, 0.0); - PhaseID.resize(nsp, 0); - SpName.resize(nsp, ""); - WtSpecies.resize(nsp, 0.0); - Charge.resize(nsp, 0.0); - NSPECIES0 = nsp; - if (nspecies > NSPECIES0) { - nspecies = NSPECIES0; - plogf("shouldn't be here\n"); + NE0 = ne; + if (ne <= 0) { + plogf("number of elements is zero or neg\n"); + exit(-1); + } + NPHASE0 = NPhase; + if (NPhase <= 0) { + plogf("number of phases is zero or neg\n"); + exit(-1); + } + if (nspecies < NPhase) { + plogf("number of species is less than number of phases\n"); exit(-1); } - } -} -// 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 VCS_PROB::resizeElements(int nel, int force) { - if (force || nel > NE0) { - gai.resize(nel, 0.0); - FormulaMatrix.resize(nel, NSPECIES0, 0.0); - ElName.resize(nel, ""); - m_elType.resize(nel, VCS_ELEM_TYPE_ABSPOS); - ElActive.resize(nel, 1); - NE0 = nel; - if (ne > NE0) ne = NE0; - } -} - -// Calculate the element abundance vector -/* - * Calculates the element abundance vectors from the mole - * numbers - */ -void VCS_PROB::set_gai () -{ - int kspec, j; - double *ElemAbund = VCS_DATA_PTR(gai); - double *const *const fm = FormulaMatrix.baseDataAddr(); - vcs_dzero(ElemAbund, ne); - - for (j = 0; j < ne; j++) { - for (kspec = 0; kspec < nspecies; kspec++) { - ElemAbund[j] += fm[j][kspec] * w[kspec]; + m_gibbsSpecies.resize(nspecies, 0.0); + w.resize(nspecies, 0.0); + mf.resize(nspecies, 0.0); + gai.resize(ne, 0.0); + FormulaMatrix.resize(ne, nspecies, 0.0); + SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); + VolPM.resize(nspecies, 0.0); + PhaseID.resize(nspecies, -1); + SpName.resize(nspecies, ""); + ElName.resize(ne, ""); + m_elType.resize(ne, VCS_ELEM_TYPE_ABSPOS); + ElActive.resize(ne, 1); + WtSpecies.resize(nspecies, 0.0); + Charge.resize(nspecies, 0.0); + SpeciesThermo.resize(nspecies,0); + for (int kspec = 0; kspec < nspecies; kspec++) { + VCS_SPECIES_THERMO *ts_tmp = new VCS_SPECIES_THERMO(0, 0); + if (ts_tmp == 0) { + plogf("Failed to init a ts struct\n"); + exit(-1); + } + SpeciesThermo[kspec] = ts_tmp; + } + VPhaseList.resize(nph, 0); + for (int iphase = 0; iphase < NPhase; iphase++) { + VPhaseList[iphase] = new vcs_VolPhase(); } } -} - -/*****************************************************************************/ -static void print_space(int num) { - for (int j = 0; j < num; j++) (void) plogf(" "); -} -/*****************************************************************************/ -static void print_char(const char letter, const int num) { - for (int i = 0; i < num; i++) plogf("%c", letter); -} - -/***************************************************************************** - * prob_report(): - * - * Print out the problem specification in all generality - * as it currently exists in the VCS_PROB object - * - */ -void VCS_PROB::prob_report(int print_lvl) { - m_printLvl = print_lvl; - int i, iphase; - vcs_VolPhase *Vphase = 0; + /**************************************************************************/ + /**************************************************************************/ + /**************************************************************************/ /* - * Printout the species information: PhaseID's and mole nums + * VCS_PROB_INPUT:destructor + * + * We need to manually free all of the arrays. */ - if (m_printLvl > 0) { - plogf("\n"); print_char('=', 80); plogf("\n"); - print_char('=', 20); - plogf(" VCS_PROB: PROBLEM STATEMENT "); - print_char('=', 31); plogf("\n"); - print_char('=', 80); plogf("\n"); - - plogf("\n"); - if (prob_type == 0) { - plogf("\tSolve a constant T, P problem:\n"); - plogf("\t\tT = %g K\n", T); - double pres_atm = PresPA / 1.01325E5; - - plogf("\t\tPres = %g atm\n", pres_atm); - } else { - plogf("\tUnknown problem type\n"); - exit(-1); + VCS_PROB::~VCS_PROB() { + for (int i = 0; i < nspecies; i++) { + delete SpeciesThermo[i]; + SpeciesThermo[i] = 0; } - plogf("\n"); - plogf(" Phase IDs of species\n"); - plogf(" species phaseID phaseName "); - plogf(" Initial_Estimated_Moles Species_Type\n"); - for (i = 0; i < nspecies; i++) { - Vphase = VPhaseList[PhaseID[i]]; - plogf("%16s %5d %16s", SpName[i].c_str(), PhaseID[i], - Vphase->PhaseName.c_str()); - if (iest >= 0) plogf(" %-10.5g", w[i]); - else plogf(" N/A"); - if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { - plogf(" Mol_Num"); - } else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Voltage"); + for (int iph = 0; iph < NPhase; iph++) { + delete VPhaseList[iph]; + VPhaseList[iph] = 0; + } + } + + // 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 VCS_PROB::resizePhase(int nPhase, int force) { + if (force || nPhase > NPHASE0) { + NPHASE0 = nPhase; + } + } + + // 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 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 VCS_PROB::resizeSpecies(int nsp, int force) { + if (force || nsp > NSPECIES0) { + m_gibbsSpecies.resize(nsp, 0.0); + w.resize(nsp, 0.0); + mf.resize(nsp, 0.0); + FormulaMatrix.resize(NE0, nsp, 0.0); + SpeciesUnknownType.resize(nsp, VCS_SPECIES_TYPE_MOLNUM); + VolPM.resize(nsp, 0.0); + PhaseID.resize(nsp, 0); + SpName.resize(nsp, ""); + WtSpecies.resize(nsp, 0.0); + Charge.resize(nsp, 0.0); + NSPECIES0 = nsp; + if (nspecies > NSPECIES0) { + nspecies = NSPECIES0; + plogf("shouldn't be here\n"); + exit(-1); + } + } + } + + // 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 VCS_PROB::resizeElements(int nel, int force) { + if (force || nel > NE0) { + gai.resize(nel, 0.0); + FormulaMatrix.resize(nel, NSPECIES0, 0.0); + ElName.resize(nel, ""); + m_elType.resize(nel, VCS_ELEM_TYPE_ABSPOS); + ElActive.resize(nel, 1); + NE0 = nel; + if (ne > NE0) ne = NE0; + } + } + + // Calculate the element abundance vector + /* + * Calculates the element abundance vectors from the mole + * numbers + */ + void VCS_PROB::set_gai () + { + int kspec, j; + double *ElemAbund = VCS_DATA_PTR(gai); + double *const *const fm = FormulaMatrix.baseDataAddr(); + vcs_dzero(ElemAbund, ne); + + for (j = 0; j < ne; j++) { + for (kspec = 0; kspec < nspecies; kspec++) { + ElemAbund[j] += fm[j][kspec] * w[kspec]; + } + } + } + + /*****************************************************************************/ + static void print_space(int num) { + for (int j = 0; j < num; j++) (void) plogf(" "); + } + /*****************************************************************************/ + static void print_char(const char letter, const int num) { + for (int i = 0; i < num; i++) plogf("%c", letter); + } + + /***************************************************************************** + * prob_report(): + * + * Print out the problem specification in all generality + * as it currently exists in the VCS_PROB object + * + */ + void VCS_PROB::prob_report(int print_lvl) { + m_printLvl = print_lvl; + int i, iphase; + vcs_VolPhase *Vphase = 0; + /* + * Printout the species information: PhaseID's and mole nums + */ + if (m_printLvl > 0) { + plogf("\n"); print_char('=', 80); plogf("\n"); + print_char('=', 20); + plogf(" VCS_PROB: PROBLEM STATEMENT "); + print_char('=', 31); plogf("\n"); + print_char('=', 80); plogf("\n"); + + plogf("\n"); + if (prob_type == 0) { + plogf("\tSolve a constant T, P problem:\n"); + plogf("\t\tT = %g K\n", T); + double pres_atm = PresPA / 1.01325E5; + + plogf("\t\tPres = %g atm\n", pres_atm); } else { - plogf(" "); + plogf("\tUnknown problem type\n"); + exit(-1); } plogf("\n"); - } - - /* - * Printout of the Phase structure information - */ - plogf("\n"); print_char('-', 80); plogf("\n"); - plogf(" Information about phases\n"); - plogf(" PhaseName PhaseNum SingSpec GasPhase " - " EqnState NumSpec"); - plogf(" TMolesInert TKmoles\n"); - - for (iphase = 0; iphase < NPhase; iphase++) { - Vphase = VPhaseList[iphase]; - std::string EOS_cstr = string16_EOSType(Vphase->EqnState); - plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(), - Vphase->VP_ID, Vphase->SingleSpecies, Vphase->m_gasPhase); - plogf("%16s %8d %16e ", EOS_cstr.c_str(), - Vphase->NVolSpecies, Vphase->TMolesInert); - if (iest >= 0) plogf("%16e\n", Vphase->TotalMoles()); - else plogf(" N/A\n"); - } - - plogf("\nElemental Abundances: "); - plogf(" Target_kmol ElemType ElActive\n"); - double fac = 1.0; - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - //fac = 1.0E3; - fac = 1.0; - } - for (i = 0; i < ne; ++i) { - print_space(26); plogf("%-2.2s", ElName[i].c_str()); - plogf("%20.12E ", fac * gai[i]); - plogf("%3d %3d\n", m_elType[i], ElActive[i]); - } - - plogf("\nChemical Potentials: "); - if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - plogf("(unitless)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) { - plogf("(kcal/gmol)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) { - plogf("(kJ/gmol)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) { - plogf("(Kelvin)"); - } else if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - plogf("(J/kmol)"); - } - plogf("\n"); - plogf(" Species (phase) " - " SS0ChemPot StarChemPot\n"); - for (iphase = 0; iphase < NPhase; iphase++) { - Vphase = VPhaseList[iphase]; - Vphase->G0_calc(T); - Vphase->setState_TP(T, PresPA); - Vphase->GStar_calc(); - for (int kindex = 0; kindex < Vphase->NVolSpecies; kindex++) { - int kglob = Vphase->IndSpecies[kindex]; - plogf("%16s ", SpName[kglob].c_str()); - if (kindex == 0) { - plogf("%16s", Vphase->PhaseName.c_str()); + plogf(" Phase IDs of species\n"); + plogf(" species phaseID phaseName "); + plogf(" Initial_Estimated_Moles Species_Type\n"); + for (i = 0; i < nspecies; i++) { + Vphase = VPhaseList[PhaseID[i]]; + plogf("%16s %5d %16s", SpName[i].c_str(), PhaseID[i], + Vphase->PhaseName.c_str()); + if (iest >= 0) plogf(" %-10.5g", w[i]); + else plogf(" N/A"); + if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { + plogf(" Mol_Num"); + } else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" Voltage"); } else { - plogf(" "); + plogf(" "); } - - plogf("%16g %16g\n", - Vphase->SS0ChemicalPotential[kindex], - Vphase->StarChemicalPotential[kindex] ); + plogf("\n"); } + + /* + * Printout of the Phase structure information + */ + plogf("\n"); print_char('-', 80); plogf("\n"); + plogf(" Information about phases\n"); + plogf(" PhaseName PhaseNum SingSpec GasPhase " + " EqnState NumSpec"); + plogf(" TMolesInert TKmoles\n"); + + for (iphase = 0; iphase < NPhase; iphase++) { + Vphase = VPhaseList[iphase]; + std::string EOS_cstr = string16_EOSType(Vphase->EqnState); + plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(), + Vphase->VP_ID, Vphase->SingleSpecies, Vphase->m_gasPhase); + plogf("%16s %8d %16e ", EOS_cstr.c_str(), + Vphase->NVolSpecies, Vphase->TMolesInert); + if (iest >= 0) plogf("%16e\n", Vphase->TotalMoles()); + else plogf(" N/A\n"); + } + + plogf("\nElemental Abundances: "); + plogf(" Target_kmol ElemType ElActive\n"); + double fac = 1.0; + if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { + //fac = 1.0E3; + fac = 1.0; + } + for (i = 0; i < ne; ++i) { + print_space(26); plogf("%-2.2s", ElName[i].c_str()); + plogf("%20.12E ", fac * gai[i]); + plogf("%3d %3d\n", m_elType[i], ElActive[i]); + } + + plogf("\nChemical Potentials: "); + if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { + plogf("(unitless)"); + } else if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) { + plogf("(kcal/gmol)"); + } else if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) { + plogf("(kJ/gmol)"); + } else if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) { + plogf("(Kelvin)"); + } else if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { + plogf("(J/kmol)"); + } + plogf("\n"); + plogf(" Species (phase) " + " SS0ChemPot StarChemPot\n"); + for (iphase = 0; iphase < NPhase; iphase++) { + Vphase = VPhaseList[iphase]; + Vphase->G0_calc(T); + Vphase->setState_TP(T, PresPA); + for (int kindex = 0; kindex < Vphase->NVolSpecies; kindex++) { + int kglob = Vphase->IndSpecies[kindex]; + plogf("%16s ", SpName[kglob].c_str()); + if (kindex == 0) { + plogf("%16s", Vphase->PhaseName.c_str()); + } else { + plogf(" "); + } + + plogf("%16g %16g\n", + Vphase->SS0ChemicalPotential[kindex], + Vphase->GStar_calc_one(kindex)); + } + } + plogf("\n"); print_char('=', 80); plogf("\n"); + print_char('=', 20); plogf(" VCS_PROB: END OF PROBLEM STATEMENT "); + print_char('=', 24); plogf("\n"); + print_char('=', 80); plogf("\n\n"); } - plogf("\n"); print_char('=', 80); plogf("\n"); - print_char('=', 20); plogf(" VCS_PROB: END OF PROBLEM STATEMENT "); - print_char('=', 24); plogf("\n"); - print_char('=', 80); plogf("\n\n"); } -} -// 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 - * storred 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 VCS_PROB::addPhaseElements(vcs_VolPhase *volPhase) { - int e, eVP; - int foundPos = -1; - int neVP = volPhase->nElemConstraints; - std::string en; - std::string enVP; + // Add elements to the local element list /* - * Loop through the elements in the vol phase object + * 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 + * storred 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 */ - for (eVP = 0; eVP < neVP; eVP++) { - foundPos = -1; - enVP = volPhase->ElName[eVP]; + void VCS_PROB::addPhaseElements(vcs_VolPhase *volPhase) { + int e, eVP; + int foundPos = -1; + int neVP = volPhase->nElemConstraints; + std::string en; + std::string enVP; /* - * Search for matches with the existing elements. - * If found, then fill in the entry in the global - * mapping array. + * Loop through the elements in the vol phase object */ - for (e = 0; e < ne; e++) { - en = ElName[e]; - if (!strcmp(enVP.c_str(), en.c_str())) { + for (eVP = 0; eVP < neVP; eVP++) { + foundPos = -1; + enVP = volPhase->ElName[eVP]; + /* + * Search for matches with the existing elements. + * If found, then fill in the entry in the global + * mapping array. + */ + for (e = 0; e < ne; e++) { + en = ElName[e]; + if (!strcmp(enVP.c_str(), en.c_str())) { + volPhase->ElGlobalIndex[eVP] = e; + foundPos = e; + } + } + if (foundPos == -1) { + int elType = volPhase->m_elType[eVP]; + int elactive = volPhase->ElActive[eVP]; + e = addElement(enVP.c_str(), elType, elactive); volPhase->ElGlobalIndex[eVP] = e; - foundPos = e; } } - if (foundPos == -1) { - int elType = volPhase->m_elType[eVP]; - int elactive = volPhase->ElActive[eVP]; - e = addElement(enVP.c_str(), elType, elactive); - volPhase->ElGlobalIndex[eVP] = e; - } } -} -// 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 - */ -int VCS_PROB::addElement(const char *elNameNew, int elType, int elactive) { + // 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 + */ + int VCS_PROB::addElement(const char *elNameNew, int elType, int elactive) { if (!elNameNew) { plogf("error: element must have a name\n"); exit(-1); @@ -421,207 +420,207 @@ int VCS_PROB::addElement(const char *elNameNew, int elType, int elactive) { m_elType[ne-1] = elType; ElActive[ne-1] = elactive; return (ne - 1); -} - -// 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 - * - */ -int VCS_PROB::addOnePhaseSpecies(vcs_VolPhase *volPhase, int k, int kT) { - int e, eVP; - if (kT > nspecies) { - /* - * Need to expand the number of species here - */ - plogf("Shouldn't be here\n"); - exit(-1); } - double *const *const fm = volPhase->FormulaMatrix.baseDataAddr(); - for (eVP = 0; eVP < volPhase->nElemConstraints; eVP++) { - e = volPhase->ElGlobalIndex[eVP]; - if (e < 0) { + + // 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 + * + */ + int VCS_PROB::addOnePhaseSpecies(vcs_VolPhase *volPhase, int k, int kT) { + int e, eVP; + if (kT > nspecies) { + /* + * Need to expand the number of species here + */ + plogf("Shouldn't be here\n"); exit(-1); } - FormulaMatrix[e][kT] = fm[eVP][k]; - } - /* - * Tell the phase object about the current position of the - * species within the global species vector - */ - volPhase->IndSpecies[k] = kT; - return kT; -} - -void VCS_PROB::reportCSV(const std::string &reportFile) { - int k; - int istart; - - double vol = 0.0; - string sName; - int nphase = NPhase; - - FILE * FP = fopen(reportFile.c_str(), "w"); - if (!FP) { - plogf("Failure to open file\n"); - exit(-1); - } - double Temp = T; - - std::vector volPM(nspecies, 0.0); - std::vector activity(nspecies, 0.0);; - std::vector ac(nspecies, 0.0);; - std::vector mu(nspecies, 0.0);; - std::vector mu0(nspecies, 0.0);; - std::vector molalities(nspecies, 0.0);; - - - vol = 0.0; - int iK = 0; - for (int iphase = 0; iphase < nphase; iphase++) { - istart = iK; - vcs_VolPhase *volP = VPhaseList[iphase]; - //const Cantera::ThermoPhase *tptr = volP->ptrThermoPhase(); - int nSpeciesPhase = volP->NVolSpecies; - volPM.resize(nSpeciesPhase, 0.0); - volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); - - double TMolesPhase = volP->TotalMoles(); - double VolPhaseVolumes = 0.0; - for (k = 0; k < nSpeciesPhase; k++) { - iK++; - VolPhaseVolumes += volPM[istart + k] * mf[istart + k]; - } - VolPhaseVolumes *= TMolesPhase; - vol += VolPhaseVolumes; - } - - fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT" - " -----------------------------\n"); - fprintf(FP,"Temperature = %11.5g kelvin\n", Temp); - fprintf(FP,"Pressure = %11.5g Pascal\n", PresPA); - fprintf(FP,"Total Volume = %11.5g m**3\n", vol); - fprintf(FP,"Number Basis optimizations = %d\n", m_NumBasisOptimizations); - fprintf(FP,"Number VCS iterations = %d\n", m_Iterations); - - iK = 0; - for (int iphase = 0; iphase < nphase; iphase++) { - istart = iK; - - vcs_VolPhase *volP = VPhaseList[iphase]; - const Cantera::ThermoPhase *tp = volP->ptrThermoPhase(); - string phaseName = volP->PhaseName; - int nSpeciesPhase = volP->NVolSpecies; - volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); - double TMolesPhase = volP->TotalMoles(); - //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase)); - activity.resize(nSpeciesPhase, 0.0); - ac.resize(nSpeciesPhase, 0.0); - - mu0.resize(nSpeciesPhase, 0.0); - mu.resize(nSpeciesPhase, 0.0); - volPM.resize(nSpeciesPhase, 0.0); - molalities.resize(nSpeciesPhase, 0.0); - - int actConvention = tp->activityConvention(); - tp->getActivities(VCS_DATA_PTR(activity)); - tp->getActivityCoefficients(VCS_DATA_PTR(ac)); - tp->getStandardChemPotentials(VCS_DATA_PTR(mu0)); - - tp->getPartialMolarVolumes(VCS_DATA_PTR(volPM)); - tp->getChemPotentials(VCS_DATA_PTR(mu)); - double VolPhaseVolumes = 0.0; - for (k = 0; k < nSpeciesPhase; k++) { - VolPhaseVolumes += volPM[k] * mf[istart + k]; - } - VolPhaseVolumes *= TMolesPhase; - vol += VolPhaseVolumes; - - - if (actConvention == 1) { - const Cantera::MolalityVPSSTP *mTP = static_cast(tp); - tp->getChemPotentials(VCS_DATA_PTR(mu)); - mTP->getMolalities(VCS_DATA_PTR(molalities)); - tp->getChemPotentials(VCS_DATA_PTR(mu)); - - if (iphase == 0) { - fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, " - "Molalities, ActCoeff, Activity," - "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n"); - - fprintf(FP," , , (kmol), , " - " , , ," - " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n"); - } - for (k = 0; k < nSpeciesPhase; k++) { - sName = tp->speciesName(k); - fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e," - "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n", - sName.c_str(), - phaseName.c_str(), TMolesPhase, - mf[istart + k], molalities[k], ac[k], activity[k], - mu0[k]*1.0E-6, mu[k]*1.0E-6, - mf[istart + k] * TMolesPhase, - volPM[k], VolPhaseVolumes ); - } - - } else { - if (iphase == 0) { - fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, " - "Molalities, ActCoeff, Activity," - " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n"); - - fprintf(FP," , , (kmol), , " - " , , ," - " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n"); - } - for (k = 0; k < nSpeciesPhase; k++) { - molalities[k] = 0.0; - } - for (k = 0; k < nSpeciesPhase; k++) { - sName = tp->speciesName(k); - fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, " - "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n", - sName.c_str(), - phaseName.c_str(), TMolesPhase, - mf[istart + k], molalities[k], ac[k], - activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6, - mf[istart + k] * TMolesPhase, - volPM[k], VolPhaseVolumes ); - } - } - -#ifdef DEBUG_MODE - /* - * Check consistency: These should be equal - */ - tp->getChemPotentials(VCS_DATA_PTR(m_gibbsSpecies)+istart); - for (k = 0; k < nSpeciesPhase; k++) { - if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) { - fprintf(FP,"ERROR: incompatibility!\n"); - fclose(FP); - plogf("ERROR: incompatibility!\n"); + double *const *const fm = volPhase->FormulaMatrix.baseDataAddr(); + for (eVP = 0; eVP < volPhase->nElemConstraints; eVP++) { + e = volPhase->ElGlobalIndex[eVP]; + if (e < 0) { exit(-1); } + FormulaMatrix[e][kT] = fm[eVP][k]; } -#endif - iK += nSpeciesPhase; + /* + * Tell the phase object about the current position of the + * species within the global species vector + */ + volPhase->IndSpecies[k] = kT; + return kT; } - fclose(FP); -} + + void VCS_PROB::reportCSV(const std::string &reportFile) { + int k; + int istart; + + double vol = 0.0; + string sName; + int nphase = NPhase; + + FILE * FP = fopen(reportFile.c_str(), "w"); + if (!FP) { + plogf("Failure to open file\n"); + exit(-1); + } + double Temp = T; + + std::vector volPM(nspecies, 0.0); + std::vector activity(nspecies, 0.0);; + std::vector ac(nspecies, 0.0);; + std::vector mu(nspecies, 0.0);; + std::vector mu0(nspecies, 0.0);; + std::vector molalities(nspecies, 0.0);; -void VCS_PROB::setDebugPrintLvl(int lvl) { - vcs_debug_print_lvl = lvl; -} + vol = 0.0; + int iK = 0; + for (int iphase = 0; iphase < nphase; iphase++) { + istart = iK; + vcs_VolPhase *volP = VPhaseList[iphase]; + //const Cantera::ThermoPhase *tptr = volP->ptrThermoPhase(); + int nSpeciesPhase = volP->NVolSpecies; + volPM.resize(nSpeciesPhase, 0.0); + volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); + + double TMolesPhase = volP->TotalMoles(); + double VolPhaseVolumes = 0.0; + for (k = 0; k < nSpeciesPhase; k++) { + iK++; + VolPhaseVolumes += volPM[istart + k] * mf[istart + k]; + } + VolPhaseVolumes *= TMolesPhase; + vol += VolPhaseVolumes; + } + + fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT" + " -----------------------------\n"); + fprintf(FP,"Temperature = %11.5g kelvin\n", Temp); + fprintf(FP,"Pressure = %11.5g Pascal\n", PresPA); + fprintf(FP,"Total Volume = %11.5g m**3\n", vol); + fprintf(FP,"Number Basis optimizations = %d\n", m_NumBasisOptimizations); + fprintf(FP,"Number VCS iterations = %d\n", m_Iterations); + + iK = 0; + for (int iphase = 0; iphase < nphase; iphase++) { + istart = iK; + + vcs_VolPhase *volP = VPhaseList[iphase]; + const Cantera::ThermoPhase *tp = volP->ptrThermoPhase(); + string phaseName = volP->PhaseName; + int nSpeciesPhase = volP->NVolSpecies; + volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); + double TMolesPhase = volP->TotalMoles(); + //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase)); + activity.resize(nSpeciesPhase, 0.0); + ac.resize(nSpeciesPhase, 0.0); + + mu0.resize(nSpeciesPhase, 0.0); + mu.resize(nSpeciesPhase, 0.0); + volPM.resize(nSpeciesPhase, 0.0); + molalities.resize(nSpeciesPhase, 0.0); + + int actConvention = tp->activityConvention(); + tp->getActivities(VCS_DATA_PTR(activity)); + tp->getActivityCoefficients(VCS_DATA_PTR(ac)); + tp->getStandardChemPotentials(VCS_DATA_PTR(mu0)); + + tp->getPartialMolarVolumes(VCS_DATA_PTR(volPM)); + tp->getChemPotentials(VCS_DATA_PTR(mu)); + double VolPhaseVolumes = 0.0; + for (k = 0; k < nSpeciesPhase; k++) { + VolPhaseVolumes += volPM[k] * mf[istart + k]; + } + VolPhaseVolumes *= TMolesPhase; + vol += VolPhaseVolumes; + + + if (actConvention == 1) { + const Cantera::MolalityVPSSTP *mTP = static_cast(tp); + tp->getChemPotentials(VCS_DATA_PTR(mu)); + mTP->getMolalities(VCS_DATA_PTR(molalities)); + tp->getChemPotentials(VCS_DATA_PTR(mu)); + + if (iphase == 0) { + fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, " + "Molalities, ActCoeff, Activity," + "ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume\n"); + + fprintf(FP," , , (kmol), , " + " , , ," + " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n"); + } + for (k = 0; k < nSpeciesPhase; k++) { + sName = tp->speciesName(k); + fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e," + "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n", + sName.c_str(), + phaseName.c_str(), TMolesPhase, + mf[istart + k], molalities[k], ac[k], activity[k], + mu0[k]*1.0E-6, mu[k]*1.0E-6, + mf[istart + k] * TMolesPhase, + volPM[k], VolPhaseVolumes ); + } + + } else { + if (iphase == 0) { + fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, " + "Molalities, ActCoeff, Activity," + " ChemPotSS0, ChemPot, mole_num, PMVol, Phase_Volume\n"); + + fprintf(FP," , , (kmol), , " + " , , ," + " (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n"); + } + for (k = 0; k < nSpeciesPhase; k++) { + molalities[k] = 0.0; + } + for (k = 0; k < nSpeciesPhase; k++) { + sName = tp->speciesName(k); + fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, " + "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n", + sName.c_str(), + phaseName.c_str(), TMolesPhase, + mf[istart + k], molalities[k], ac[k], + activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6, + mf[istart + k] * TMolesPhase, + volPM[k], VolPhaseVolumes ); + } + } + +#ifdef DEBUG_MODE + /* + * Check consistency: These should be equal + */ + tp->getChemPotentials(VCS_DATA_PTR(m_gibbsSpecies)+istart); + for (k = 0; k < nSpeciesPhase; k++) { + if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) { + fprintf(FP,"ERROR: incompatibility!\n"); + fclose(FP); + plogf("ERROR: incompatibility!\n"); + exit(-1); + } + } +#endif + iK += nSpeciesPhase; + } + fclose(FP); + } + + + void VCS_PROB::setDebugPrintLvl(int lvl) { + vcs_debug_print_lvl = lvl; + } } diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index d7c697676..4593feec5 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -324,11 +324,16 @@ namespace VCSnonideal { * We don't need to call single species phases; */ if (!Vphase->SingleSpecies && !Vphase->isIdealSoln()) { + + /* + * update the mole numbers + */ + Vphase->setMolesFromVCS(VCS_STATECALC_OLD, moleSpeciesVCS); /* * update the Ln Act Coeff jacobian entries with respect to the * mole number of species in the phase */ - Vphase->updateLnActCoeffJac(moleSpeciesVCS); + Vphase->updateLnActCoeffJac(); /* * Download the resulting calculation into the full vector * -> This scatter calculation is carried out in the diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 0e42ad972..bee5b21b4 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -971,9 +971,7 @@ namespace VCSnonideal { pubPhase->setMoleFractions(VCS_DATA_PTR(vPhase->moleFractions())); for (int k = 0; k < pubPhase->NVolSpecies; k++) { kT = pubPhase->IndSpecies[k]; - //pubPhase->Xmol[k] = vPhase->Xmol[k]; pubPhase->SS0ChemicalPotential[k] = vPhase->SS0ChemicalPotential[k]; - pubPhase->StarChemicalPotential[k] = vPhase->StarChemicalPotential[k]; pubPhase->StarMolarVol[k] = vPhase->StarMolarVol[k]; pubPhase->PartialMolarVol[k] = vPhase->PartialMolarVol[k]; pubPhase->ActCoeff[k] = vPhase->ActCoeff[k]; diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 58760f599..5d65034e1 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -1323,6 +1323,8 @@ private: void vcs_TCounters_report(int timing_print_lvl = 1); + void vcs_setMoleNumVolPhases(bool upToDate, int stateCalc); + public: //! value of the number of species used to malloc data structures int NSPECIES0; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 849fee22c..4931c59f8 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -1099,6 +1099,7 @@ namespace VCSnonideal { * solution values. We only calculate a subset of these, because * we have only updated a subset of the W(). */ + vcs_setMoleNumVolPhases(false, VCS_STATECALC_NEW); vcs_updateVP(VCS_STATECALC_NEW); vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); @@ -3273,7 +3274,6 @@ namespace VCSnonideal { } /*****************************************************************************/ - // This subroutine calculates reaction free energy changes for // all noncomponent formation reactions. /* @@ -5361,4 +5361,24 @@ namespace VCSnonideal { return dx; } /*******************************************************************************/ + + + void VCS_SOLVE::vcs_setMoleNumVolPhases(bool upToDate, int stateCalc) { + int iph; + vcs_VolPhase *Vphase; + if (!upToDate) { + for (iph = 0; iph < m_numPhases; iph++) { + Vphase = m_VolPhaseList[iph]; + Vphase->m_UpToDate = false; + } + } else { + for (iph = 0; iph < m_numPhases; iph++) { + Vphase = m_VolPhaseList[iph]; + Vphase->m_UpToDate = true; + Vphase->m_vcsStateStatus = stateCalc; + } + } + } + + } diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index 0258dd270..1ecb5244d 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -187,7 +187,8 @@ double VCS_SPECIES_THERMO::GStar_R_calc(int kglob, double TKelvin, if (UseCanteraCalls) { AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency"); int kspec = IndexSpeciesPhase; - fe = OwningPhase->GStar_calc_one(kspec, TKelvin, pres); + OwningPhase->setState_TP(TKelvin, pres); + fe = OwningPhase->GStar_calc_one(kspec); double R = vcsUtil_gasConstant(m_VCS_UnitsFormat); fe /= R; } else {