From 451520d23829c6aa091751e5c672e9b73247ec2c Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 11 Jun 2008 17:14:33 +0000 Subject: [PATCH] Started to reorganize vcs_VolPhase to be more efficient. --- Cantera/src/equil/vcs_Gibbs.cpp | 6 +- Cantera/src/equil/vcs_MultiPhaseEquil.cpp | 10 +- Cantera/src/equil/vcs_TP.cpp | 2 +- Cantera/src/equil/vcs_VolPhase.cpp | 191 ++++-- Cantera/src/equil/vcs_VolPhase.h | 692 ++++++++++++---------- Cantera/src/equil/vcs_prob.cpp | 9 +- Cantera/src/equil/vcs_rxnadj.cpp | 2 +- Cantera/src/equil/vcs_solve.cpp | 13 +- Cantera/src/equil/vcs_solve_TP.cpp | 42 +- 9 files changed, 561 insertions(+), 406 deletions(-) diff --git a/Cantera/src/equil/vcs_Gibbs.cpp b/Cantera/src/equil/vcs_Gibbs.cpp index 494cb2741..eddd0df55 100644 --- a/Cantera/src/equil/vcs_Gibbs.cpp +++ b/Cantera/src/equil/vcs_Gibbs.cpp @@ -44,7 +44,7 @@ namespace VCSnonideal { if ((TPhInertMoles[iph] > 0.0) && (tPhMoles[iph] > 0.0)) { g += TPhInertMoles[iph] * log(TPhInertMoles[iph] / tPhMoles[iph]); - if (Vphase->GasPhase) { + if (Vphase->m_gasPhase) { g += TPhInertMoles[iph] * log(m_pressurePA/(1.01325E5)); } } @@ -85,7 +85,7 @@ namespace VCSnonideal { phaseMols += TPhInertMoles[iphase]; g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols); vcs_VolPhase *Vphase = m_VolPhaseList[iphase]; - if (Vphase->GasPhase == iphase) { + if (Vphase->m_gasPhase == iphase) { g += TPhInertMoles[iphase] * log(m_pressurePA/1.01325E5); } } @@ -94,5 +94,3 @@ namespace VCSnonideal { } } - - diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index e910afcc0..e816eef35 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -1145,7 +1145,7 @@ namespace VCSnonideal { vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; VolPhase->resize(iphase, nSpPhase, phaseName.c_str(), 0.0); - VolPhase->GasPhase = gasPhase; + VolPhase->m_gasPhase = gasPhase; /* * Tell the vcs_VolPhase pointer about cantera */ @@ -1385,7 +1385,7 @@ namespace VCSnonideal { } } - VolPhase->setMolesFromVCS(VCS_DATA_PTR(vprob->w)); + VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w)); /* * Now, calculate a sample naught gibbs free energy calculation * at the specified temperature. @@ -1443,7 +1443,7 @@ namespace VCSnonideal { std::string sEOS = string16_EOSType(VolPhase->EqnState); plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), VolPhase->VP_ID, VolPhase->SingleSpecies, - VolPhase->GasPhase, sEOS.c_str(), + VolPhase->m_gasPhase, sEOS.c_str(), VolPhase->NVolSpecies, VolPhase->TMolesInert ); plogf("%16e\n", VolPhase->TotalMoles()); } @@ -1511,7 +1511,7 @@ namespace VCSnonideal { int kglob = volPhase->IndSpecies[volPhase->m_phiVarIndex]; vprob->w[kglob] = tPhase->electricPotential(); } - volPhase->setMolesFromVCS(VCS_DATA_PTR(vprob->w)); + volPhase->setMolesFromVCS(VCS_STATECALC_OLD, VCS_DATA_PTR(vprob->w)); if (volPhase->TotalMoles() > 0.0) { volPhase->Existence = 1; } else { @@ -1561,7 +1561,7 @@ namespace VCSnonideal { std::string sEOS = string16_EOSType(VolPhase->EqnState); plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), VolPhase->VP_ID, VolPhase->SingleSpecies, - VolPhase->GasPhase, sEOS.c_str(), + VolPhase->m_gasPhase, sEOS.c_str(), VolPhase->NVolSpecies, VolPhase->TMolesInert ); plogf("%16e\n", VolPhase->TotalMoles() ); } diff --git a/Cantera/src/equil/vcs_TP.cpp b/Cantera/src/equil/vcs_TP.cpp index 3204f7c98..03904cc26 100644 --- a/Cantera/src/equil/vcs_TP.cpp +++ b/Cantera/src/equil/vcs_TP.cpp @@ -155,7 +155,7 @@ namespace VCSnonideal { for (int iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* vph = m_VolPhaseList[iph]; vph->setState_TP(m_temperature, m_pressurePA); - vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies)); + vph->sendToVCS_GStar(VCS_DATA_PTR(m_SSfeSpecies)); } if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 1c7efb587..916f0831b 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -13,6 +13,7 @@ #include "vcs_internal.h" #include "vcs_SpeciesProperties.h" #include "vcs_species_thermo.h" +#include "vcs_solve.h" #include "ThermoPhase.h" #include "mix_defs.h" @@ -22,18 +23,18 @@ namespace VCSnonideal { - /***************************************************************************** + /**************************************************************************** * * vcs_VolPhase(): * * Constructor for the VolPhase object. */ - vcs_VolPhase::vcs_VolPhase() : + vcs_VolPhase::vcs_VolPhase(VCS_SOLVE * owningSolverObject) : + m_owningSolverObject(0), VP_ID(-1), Domain_ID(-1), SingleSpecies(true), - GasPhase(false), - LiqPhase(false), + m_gasPhase(false), EqnState(VCS_EOS_CONSTANT), nElemConstraints(0), ChargeNeutralityElement(-1), @@ -53,8 +54,10 @@ namespace VCSnonideal { m_useCanteraCalls(false), TP_ptr(0), TMoles(0.0), - Vol(0.0), + m_totalVol(0.0), + m_vcsStateStatus(VCS_STATECALC_OLD), m_phi(0.0), + m_UpToDate(false), m_UpToDate_AC(false), m_UpToDate_VolStar(false), m_UpToDate_VolPM(false), @@ -63,6 +66,7 @@ namespace VCSnonideal { Pres(1.01325E5), RefPres(1.01325E5) { + m_owningSolverObject = owningSolverObject; } /* @@ -88,11 +92,11 @@ namespace VCSnonideal { * The assignment operator does most of the work. */ vcs_VolPhase::vcs_VolPhase(const vcs_VolPhase& b) : + m_owningSolverObject(b.m_owningSolverObject), VP_ID(b.VP_ID), Domain_ID(b.Domain_ID), SingleSpecies(b.SingleSpecies), - GasPhase(b.GasPhase), - LiqPhase(b.LiqPhase), + m_gasPhase(b.m_gasPhase), EqnState(b.EqnState), nElemConstraints(b.nElemConstraints), ChargeNeutralityElement(b.ChargeNeutralityElement), @@ -110,8 +114,10 @@ namespace VCSnonideal { TP_ptr(b.TP_ptr), TMoles(b.TMoles), m_phiVarIndex(-1), - Vol(b.Vol), + m_totalVol(b.m_totalVol), + m_vcsStateStatus(VCS_STATECALC_OLD), m_phi(b.m_phi), + m_UpToDate(false), m_UpToDate_AC(false), m_UpToDate_VolStar(false), m_UpToDate_VolPM(false), @@ -137,11 +143,11 @@ namespace VCSnonideal { if (&b != this) { int old_num = NVolSpecies; + m_owningSolverObject = b.m_owningSolverObject; VP_ID = b.VP_ID; Domain_ID = b.Domain_ID; SingleSpecies = b.SingleSpecies; - GasPhase = b.GasPhase; - LiqPhase = b.LiqPhase; + m_gasPhase = b.m_gasPhase; EqnState = b.EqnState; NVolSpecies = b.NVolSpecies; @@ -195,7 +201,7 @@ namespace VCSnonideal { } m_VCS_UnitsFormat = b.m_VCS_UnitsFormat; - m_useCanteraCalls = b.m_useCanteraCalls; + m_useCanteraCalls = b.m_useCanteraCalls; /* * Do a shallow copy of the ThermoPhase object pointer. * We don't duplicate the object. @@ -220,6 +226,8 @@ namespace VCSnonideal { dLnActCoeffdMolNumber = b.dLnActCoeffdMolNumber; + m_UpToDate = false; + m_vcsStateStatus = b.m_vcsStateStatus; m_UpToDate_AC = false; m_UpToDate_VolStar = false; m_UpToDate_VolPM = false; @@ -312,6 +320,8 @@ namespace VCSnonideal { SpeciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM); + m_UpToDate = false; + m_vcsStateStatus = VCS_STATECALC_OLD; m_UpToDate_AC = false; m_UpToDate_VolStar = false; m_UpToDate_VolPM = false; @@ -360,6 +370,7 @@ namespace VCSnonideal { return(ActCoeff[kspec]); } + // Gibbs free energy calculation at a temperature for the reference state // of each species /* @@ -389,6 +400,7 @@ namespace VCSnonideal { } } } + /***********************************************************************/ // Gibbs free energy calculation at a temperature for the reference state // of a species, return a value for one species @@ -402,6 +414,7 @@ namespace VCSnonideal { G0_calc(tkelvin); return SS0ChemicalPotential[kspec]; } + /***********************************************************************/ // Gibbs free energy calculation for standard states /* @@ -411,8 +424,8 @@ namespace VCSnonideal { * @param TKelvin Current temperature * @param pres Current pressure (pascal) */ - void vcs_VolPhase::GStar_calc(double tkelvin, double pres) { - setState_TP(tkelvin, pres); + void vcs_VolPhase::GStar_calc() { + setState_TP(Temp, Pres); if (!m_UpToDate_GStar) { if (m_useCanteraCalls) { TP_ptr->getStandardChemPotentials(VCS_DATA_PTR(StarChemicalPotential)); @@ -423,12 +436,13 @@ namespace VCSnonideal { vcs_SpeciesProperties *sProp = ListSpeciesPtr[k]; VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo; StarChemicalPotential[k] = - R * (sTherm->GStar_R_calc(kglob, tkelvin, pres)); + R * (sTherm->GStar_R_calc(kglob, Temp, Pres)); } } m_UpToDate_GStar = true; } } + /***********************************************************************/ // Gibbs free energy calculation for standard state of one species /* @@ -446,9 +460,11 @@ namespace VCSnonideal { */ double vcs_VolPhase::GStar_calc_one(int kspec, double tkelvin, double pres) { - GStar_calc(tkelvin, pres); + setState_TP(tkelvin, pres); + GStar_calc(); return StarChemicalPotential[kspec]; } + /***********************************************************************/ // Set the moles within the phase /* @@ -457,15 +473,55 @@ namespace VCSnonideal { * a gather routine. * * - * @param molesSpeciesVCS array of mole numbers. Note, the indecises for species in + * @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 vcs_VolPhase::setMolesFromVCS(const double * const molesSpeciesVCS) { + void vcs_VolPhase::setMolesFromVCS(const int stateCalc, + const double * molesSpeciesVCS) { int kglob; double tmp; TMoles = TMolesInert; + + if (molesSpeciesVCS == 0) { +#ifdef DEBUG_MODE + if (m_owningSolverObject == 0) { + printf("shouldn't be here\n"); + std::exit(-1); + } +#endif + if (stateCalc == VCS_STATECALC_OLD) { + molesSpeciesVCS = VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old); + } else if (stateCalc == VCS_STATECALC_NEW) { + molesSpeciesVCS = VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new); + } +#ifdef DEBUG_MODE + else { + printf("shouldn't be here\n"); + std::exit(-1); + } +#endif + } +#ifdef DEBUG_MODE + else { + if (m_owningSolverObject) { + if (stateCalc == VCS_STATECALC_OLD) { + if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_old)) { + printf("shouldn't be here\n"); + std::exit(-1); + } + } else if (stateCalc == VCS_STATECALC_NEW) { + if (molesSpeciesVCS != VCS_DATA_PTR(m_owningSolverObject->m_molNumSpecies_new)) { + printf("shouldn't be here\n"); + std::exit(-1); + } + } + } + } +#endif + for (int k = 0; k < NVolSpecies; k++) { if (SpeciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { kglob = IndSpecies[k]; @@ -509,7 +565,14 @@ namespace VCSnonideal { if (TMolesInert > 0.0) { Existence = 2; } + /* + * Set flags indicating we are up to date with the VCS state vector. + */ + m_UpToDate = true; + m_vcsStateStatus = stateCalc; + } + /***********************************************************************/ // Set the mole fractions from a conventional mole fraction vector /* @@ -529,7 +592,10 @@ namespace VCSnonideal { } } _updateMoleFractionDependencies(); + m_UpToDate = false; + m_vcsStateStatus = VCS_STATECALC_TMP; } + /***********************************************************************/ // Updates the mole fractions in subobjects /* @@ -542,18 +608,22 @@ namespace VCSnonideal { TP_ptr->setState_PX(Pres, VCS_DATA_PTR(Xmol)); } } - m_UpToDate_AC = false; - m_UpToDate_VolPM = false; + if (!m_isIdealSoln) { + m_UpToDate_AC = false; + m_UpToDate_VolStar = false; + m_UpToDate_VolPM = false; + m_UpToDate_GStar = false; + } } // Return a const reference to the mole fraction vector in the phase const std::vector & vcs_VolPhase::moleFractions() const { return Xmol; } + /***********************************************************************/ - - //! Set the moles within the phase - /*! + // 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. @@ -564,10 +634,11 @@ namespace VCSnonideal { * to gather the species into the local contiguous vector * format. */ - void vcs_VolPhase::setMolesFromVCSCheck(const double * const molesSpeciesVCS, + void vcs_VolPhase::setMolesFromVCSCheck(const int stateCalc, + const double * molesSpeciesVCS, const double * const TPhMoles, int iphase) { - setMolesFromVCS(molesSpeciesVCS); + setMolesFromVCS(stateCalc, molesSpeciesVCS); /* * Check for consistency with TPhMoles[] */ @@ -576,14 +647,16 @@ namespace VCSnonideal { if (vcs_doubleEqual(Tcheck, TMoles)) { Tcheck = TMoles; } else { - plogf("We have a consistency problem: %21.16g %21.16g\n", + plogf("vcs_VolPhase::setMolesFromVCSCheck: " + "We have a consistency problem: %21.16g %21.16g\n", Tcheck, TMoles); std::exit(-1); } } } + /***********************************************************************/ - // Fill in an activity coefficients vector for VCS + // 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 @@ -593,7 +666,7 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - void vcs_VolPhase::sendToVCSActCoeff(double * const AC) const { + void vcs_VolPhase::sendToVCS_ActCoeff(double * const AC) const { if (!m_UpToDate_AC) { evaluateActCoeff(); } @@ -603,6 +676,7 @@ namespace VCSnonideal { AC[kglob] = ActCoeff[k]; } } + /***********************************************************************/ // Fill in the partial molar volume vector for VCS /* @@ -614,7 +688,7 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - double vcs_VolPhase::sendToVCSVolPM(double * const VolPM) const { + double vcs_VolPhase::sendToVCS_VolPM(double * const VolPM) const { if (!m_UpToDate_VolPM) { (void) VolPM_calc(); } @@ -623,8 +697,9 @@ namespace VCSnonideal { kglob = IndSpecies[k]; VolPM[kglob] = PartialMolarVol[k]; } - return Vol; + return m_totalVol; } + /***********************************************************************/ // Fill in the partial molar volume vector for VCS /* @@ -636,9 +711,10 @@ namespace VCSnonideal { * in all of the phases in a VCS problem. Only the * entries for the current phase are filled in. */ - void vcs_VolPhase::sendToVCSGStar(double * const gstar){ + void vcs_VolPhase::sendToVCS_GStar(double * const gstar){ if (!m_UpToDate_GStar) { - GStar_calc(Temp, Pres); + setState_TP(Temp, Pres); + GStar_calc(); } int kglob; for (int k = 0; k < NVolSpecies; k++) { @@ -646,7 +722,7 @@ namespace VCSnonideal { gstar[kglob] = StarChemicalPotential[k]; } } - + /***********************************************************************/ void vcs_VolPhase::setElectricPotential(double phi) { @@ -660,10 +736,12 @@ namespace VCSnonideal { m_UpToDate_VolPM = false; m_UpToDate_GStar = false; } + /***********************************************************************/ double vcs_VolPhase::electricPotential() const { return m_phi; } + /***********************************************************************/ // Sets the temperature and pressure in this object and // underlying objects @@ -693,7 +771,7 @@ namespace VCSnonideal { m_UpToDate_VolPM = false; m_UpToDate_GStar = false; } - + /***********************************************************************/ // Molar volume calculation for standard states /* @@ -722,6 +800,30 @@ namespace VCSnonideal { m_UpToDate_VolStar = true; } } + /***********************************************************************/ + + // 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 vcs_VolPhase::updateFromVCS_MoleNumbers(const int stateCalc) { + if (!m_UpToDate) { + if (stateCalc == VCS_STATECALC_OLD || stateCalc == VCS_STATECALC_NEW) { + if (m_owningSolverObject) { + setMolesFromVCS(stateCalc); + } + } + } + } + /***********************************************************************/ // Molar volume calculation for standard state of one species /* @@ -729,20 +831,21 @@ namespace VCSnonideal { * The results are held internally within the object. * Return the molar volume for one species * - * @param kspec Species number (within the phase) + * @param kspec Species number (within the phase) * @param TKelvin Current temperature * @param pres Current pressure (pascal) * * @return molar volume of the kspec species's standard * state */ - double vcs_VolPhase::VolStar_calc_one(int kspec, double tkelvin, double pres) + double vcs_VolPhase::VolStar_calc_one(int kspec, double tkelvin, + double pres) { VolStar_calc(tkelvin, pres); return StarMolarVol[kspec]; } - /****************************************************************************** + /**************************************************************************** * * VolPM_calc */ @@ -763,16 +866,16 @@ namespace VCSnonideal { } } - Vol = 0.0; + m_totalVol = 0.0; for (k = 0; k < NVolSpecies; k++) { - Vol += PartialMolarVol[k] * Xmol[k]; + m_totalVol += PartialMolarVol[k] * Xmol[k]; } - Vol *= TMoles; + m_totalVol *= TMoles; if (TMolesInert > 0.0) { - if (GasPhase) { + if (m_gasPhase) { double volI = TMolesInert * 8314.47215 * Temp / Pres; - Vol += volI; + m_totalVol += volI; } else { printf("unknown situation\n"); std::exit(-1); @@ -780,7 +883,7 @@ namespace VCSnonideal { } } m_UpToDate_VolPM = true; - return Vol; + return m_totalVol; } /* @@ -795,7 +898,7 @@ namespace VCSnonideal { * with the current values of the mole numbers. * -> This sets TMoles and Xmol[] */ - setMolesFromVCS(moleNumbersVCS); + setMolesFromVCS(VCS_STATECALC_OLD, moleNumbersVCS); /* * Evaluate the current base activity coefficients. @@ -871,7 +974,7 @@ namespace VCSnonideal { * j = id of the species mole number * k = id of the species activity coefficient */ - void vcs_VolPhase::sendToVCSLnActCoeffJac(double * const * const LnACJac_VCS) const { + void vcs_VolPhase::sendToVCS_LnActCoeffJac(double * const * const LnACJac_VCS) const { int j, k, jglob, kglob; for (j = 0; j < NVolSpecies; j++) { jglob = IndSpecies[j]; diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 6b49099ba..b221d29d6 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -62,8 +62,9 @@ namespace VCSnonideal { #define VCS_EOS_UNK_CANTERA -1 -struct VCS_SPECIES; -class vcs_SpeciesProperties; + struct VCS_SPECIES; + class vcs_SpeciesProperties; + class VCS_SOLVE; //! Phase information and Phase calculations for vcs. @@ -111,6 +112,343 @@ class vcs_SpeciesProperties; 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. + */ + 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; + + + //! 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(); + + //! 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, 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] + * + * 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 @@ -133,12 +471,10 @@ public: //! If true, this phase is a gas-phase like phase /*! - * A RTlog(p/1atm) term is added onto the chemical potential + * A RTlog(p/1atm) term is added onto the chemical potential for inert + * species if this is true. */ - int GasPhase; - - //! If true, this phase is a liquid-phase like phase*/ - int LiqPhase; + bool m_gasPhase; //! Type of the equation of state /*! @@ -153,7 +489,8 @@ public: */ int nElemConstraints; - //! This is the element number for the charge neutrality condition of the phase + //! 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 @@ -314,7 +651,7 @@ private: * If we are using Cantera, this is the * pointer to the ThermoPhase object. If not, this is null. */ - Cantera::ThermoPhase * TP_ptr; + Cantera::ThermoPhase *TP_ptr; /** * Variables Having to do with Calculated States @@ -340,7 +677,7 @@ public: /*! * units are m**3 */ - mutable double Vol; + mutable double m_totalVol; //! Vector of calculated SS0 chemical potentials for the //! current Temperature. @@ -377,12 +714,15 @@ public: */ mutable std::vector PartialMolarVol; - /** - * Vector of calculated activity coefficients for the current - * state. + + //! 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 /*! @@ -392,11 +732,25 @@ public: */ 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 @@ -434,322 +788,12 @@ private: //! Current value of the pressure for this object, and underlying objects double Pres; + public: //! Reference pressure for the phase double RefPres; - /************************************************************************* - * FUNCTIONS * - ************************************************************************/ - - //! Base constructor for the class - vcs_VolPhase(); - - //! 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. - */ - void resize(int phaseNum, int numSpecies, const char *phaseName, - double molesInert = 0.0); - - - //! 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; - - //! 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 double * const molesSpeciesVCS); - - //! 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 double * const molesSpeciesVCS, - const double * const TPhMoles, - int iphase = -1); - - //! Fill in an activity coefficients vector for VCS - /*! - * 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 sendToVCSActCoeff(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. - * - * @param TKelvin Current temperature - * @param pres Current pressure - */ - void GStar_calc(double TKelvin, double pres); - - //! 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, 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); - - - //! 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); - - //! Calculate the partial molar volumes of all species and return the - //! total volume - /*! - * Calculates these quantitites internally - * - * @return total volume - */ - double VolPM_calc() 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 sendToVCSVolPM(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 sendToVCSGStar(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] - * - * j = id of the species mole number - * k = id of the species activity coefficient - */ - void sendToVCSLnActCoeffJac(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(); }; //! Return a string representing the equation of state diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index b1892580a..8a81cd7f4 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -283,7 +283,7 @@ void VCS_PROB::prob_report(int print_lvl) { 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->GasPhase); + 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()); @@ -321,7 +321,8 @@ void VCS_PROB::prob_report(int print_lvl) { for (iphase = 0; iphase < NPhase; iphase++) { Vphase = VPhaseList[iphase]; Vphase->G0_calc(T); - Vphase->GStar_calc(T, PresPA); + 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()); @@ -491,7 +492,7 @@ void VCS_PROB::reportCSV(const std::string &reportFile) { //const Cantera::ThermoPhase *tptr = volP->ptrThermoPhase(); int nSpeciesPhase = volP->NVolSpecies; volPM.resize(nSpeciesPhase, 0.0); - volP->sendToVCSVolPM(VCS_DATA_PTR(volPM)); + volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); double TMolesPhase = volP->TotalMoles(); double VolPhaseVolumes = 0.0; @@ -519,7 +520,7 @@ void VCS_PROB::reportCSV(const std::string &reportFile) { const Cantera::ThermoPhase *tp = volP->ptrThermoPhase(); string phaseName = volP->PhaseName; int nSpeciesPhase = volP->NVolSpecies; - volP->sendToVCSVolPM(VCS_DATA_PTR(volPM)); + volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); double TMolesPhase = volP->TotalMoles(); //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase)); activity.resize(nSpeciesPhase, 0.0); diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index f3bf12dd3..d7c697676 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -334,7 +334,7 @@ namespace VCSnonideal { * -> This scatter calculation is carried out in the * vcs_VolPhase object. */ - Vphase->sendToVCSLnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr()); + Vphase->sendToVCS_LnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr()); } } } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index bec5e3a89..0e42ad972 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -209,7 +209,7 @@ namespace VCSnonideal { */ m_VolPhaseList.resize(nphase0, 0); for (iph = 0; iph < nphase0; iph++) { - m_VolPhaseList[iph] = new vcs_VolPhase(); + m_VolPhaseList[iph] = new vcs_VolPhase(this); } /* @@ -860,10 +860,10 @@ namespace VCSnonideal { retn = VCS_PUB_BAD; } - if (vPhase->GasPhase != pub_phase_ptr->GasPhase) { + if (vPhase->m_gasPhase != pub_phase_ptr->m_gasPhase) { plogf("%sGasPhase value have changed:%d %d\n", yo.c_str(), - vPhase->GasPhase, - pub_phase_ptr->GasPhase); + vPhase->m_gasPhase, + pub_phase_ptr->m_gasPhase); retn = VCS_PUB_BAD; } @@ -1068,9 +1068,8 @@ namespace VCSnonideal { for (int iphase = 0; iphase < m_numPhases; iphase++) { vcs_VolPhase *Vphase = m_VolPhaseList[iphase]; Vphase->setState_TP(tkelvin, pres); - Vphase->setMolesFromVCS(w); - double Volp = Vphase->VolPM_calc(); - (void) Vphase->sendToVCSVolPM(volPM); + Vphase->setMolesFromVCS(VCS_STATECALC_OLD, w); + double Volp = Vphase->sendToVCS_VolPM(volPM); VolTot += Volp; } return VolTot; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index fad31e504..849fee22c 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -442,8 +442,10 @@ namespace VCSnonideal { /* * Copy the old solution into the new solution as an initial guess */ - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc); - vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), + VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), + VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); /* Go find a new reaction adjustment -> @@ -1016,8 +1018,10 @@ namespace VCSnonideal { /*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/ /*************************************************************************/ /* - * We have a tentative M_DELTAMOLNUMSPECIES(L=1,MR). Now apply other criteria - * to limit it's magnitude. + * We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria + * to limit it's magnitude. + * + * */ par = 0.5; for (k = 0; k < m_numComponents; ++k) { @@ -1057,7 +1061,7 @@ namespace VCSnonideal { m_deltaMolNumSpecies[i] *= par; } for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] *= par; + m_deltaPhaseMoles[iph] *= par; } } else { par = 1.0; @@ -2111,8 +2115,8 @@ namespace VCSnonideal { } /*****************************************************************************/ - //! Change the concentration of a species by delta moles. - /*! + // Change the concentration of a species by delta moles. + /* * Make sure to conserve elements and keep track of the total kmoles in all phases. * * @@ -2267,7 +2271,8 @@ namespace VCSnonideal { /* * Adjust the total moles in a phase downwards. */ - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, + VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old)); /* @@ -2347,7 +2352,8 @@ namespace VCSnonideal { } vcs_VolPhase *Vphase = m_VolPhaseList[iph]; - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, + VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old)); /* * We may have popped a multispecies phase back @@ -2496,7 +2502,8 @@ namespace VCSnonideal { /* * Upload the state to the VP object */ - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, + VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old), iph); } /**********************************************************************************/ @@ -4371,6 +4378,7 @@ namespace VCSnonideal { vcs_VolPhase *Vphase = m_VolPhaseList[iph]; int nkk = Vphase->NVolSpecies; int k, kspec; + int stateCalc = VCS_STATECALC_OLD; #ifdef DEBUG_MODE //if (m_debug_print_lvl >= 2) { @@ -4388,8 +4396,8 @@ namespace VCSnonideal { tlogMoles = log(tMoles); } - Vphase->setMolesFromVCS(molNum); - Vphase->sendToVCSActCoeff(ac); + Vphase->setMolesFromVCS(stateCalc, molNum); + Vphase->sendToVCS_ActCoeff(ac); double phi = Vphase->electricPotential(); double Faraday_phi = m_Faraday_dim * phi; @@ -4669,8 +4677,8 @@ namespace VCSnonideal { if (!m_phaseACAreCurrent[iphase]) { Vphase = m_VolPhaseList[iphase]; if (!Vphase->SingleSpecies) { - Vphase->setMolesFromVCS(molNum); - Vphase->sendToVCSActCoeff(VCS_DATA_PTR(actCoeff_ptr)); + Vphase->setMolesFromVCS(stateCalc, molNum); + Vphase->sendToVCS_ActCoeff(VCS_DATA_PTR(actCoeff_ptr)); } m_phasePhi[iphase] = Vphase->electricPotential(); m_phaseACAreCurrent[iphase] = 1; @@ -4927,10 +4935,12 @@ namespace VCSnonideal { for (int i = 0; i < m_numPhases; i++) { Vphase = m_VolPhaseList[i]; if (vcsState == VCS_STATECALC_OLD) { - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, + VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old), i); } else if (vcsState == VCS_STATECALC_NEW) { - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_new), + Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW, + VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_tPhaseMoles_new), i); } #ifdef DEBUG_MODE