diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 9e7008a08..30eee5b42 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -467,6 +467,53 @@ namespace VCSnonideal { } /***********************************************************************/ + // 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 vcs_VolPhase::setMoleFractions(const double * const xmol) { + double sum = -1.0; + for (int k = 0; k < NVolSpecies; k++) { + Xmol[k] = xmol[k]; + sum+= xmol[k]; + } + if (std::fabs(sum) > 1.0E-13) { + for (int k = 0; k < NVolSpecies; k++) { + Xmol[k] /= sum; + } + } + _updateMoleFractionDependencies(); + m_UpToDate = false; + m_vcsStateStatus = VCS_STATECALC_TMP; + } + /***********************************************************************/ + + // Updates the mole fractions in subobjects + /* + * Whenever the mole fractions change, this routine + * should be called. + */ + void vcs_VolPhase::_updateMoleFractionDependencies() { + if (m_useCanteraCalls) { + if (TP_ptr) { + TP_ptr->setState_PX(Pres, VCS_DATA_PTR(Xmol)); + } + } + if (!m_isIdealSoln) { + m_UpToDate_AC = false; + m_UpToDate_VolPM = 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 /* * This function takes as input the mole numbers in vcs format, and @@ -575,52 +622,6 @@ namespace VCSnonideal { } /***********************************************************************/ - // 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 vcs_VolPhase::setMoleFractions(const double * const xmol) { - double sum = -1.0; - for (int k = 0; k < NVolSpecies; k++) { - Xmol[k] = xmol[k]; - sum+= xmol[k]; - } - if (std::fabs(sum) > 1.0E-13) { - for (int k = 0; k < NVolSpecies; k++) { - Xmol[k] /= sum; - } - } - _updateMoleFractionDependencies(); - m_UpToDate = false; - m_vcsStateStatus = VCS_STATECALC_TMP; - } - /***********************************************************************/ - - // Updates the mole fractions in subobjects - /* - * Whenever the mole fractions change, this routine - * should be called. - */ - void vcs_VolPhase::_updateMoleFractionDependencies() { - if (m_useCanteraCalls) { - if (TP_ptr) { - TP_ptr->setState_PX(Pres, VCS_DATA_PTR(Xmol)); - } - } - if (!m_isIdealSoln) { - m_UpToDate_AC = false; - m_UpToDate_VolPM = 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 /* * This function takes as input the mole numbers in vcs format, and @@ -635,8 +636,7 @@ namespace VCSnonideal { */ void vcs_VolPhase::setMolesFromVCSCheck(const int stateCalc, const double * molesSpeciesVCS, - const double * const TPhMoles, - int iphase) { + const double * const TPhMoles) { setMolesFromVCS(stateCalc, molesSpeciesVCS); /* * Check for consistency with TPhMoles[] @@ -655,6 +655,29 @@ namespace VCSnonideal { } /***********************************************************************/ + // 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 || (stateCalc != m_vcsStateStatus)) { + if (stateCalc == VCS_STATECALC_OLD || stateCalc == VCS_STATECALC_NEW) { + if (m_owningSolverObject) { + setMolesFromVCS(stateCalc); + } + } + } + } + /***********************************************************************/ + // Fill in an activity coefficients vector within a VCS_SOLVE object /* * This routine will calculate the activity coefficients for the @@ -799,28 +822,6 @@ namespace VCSnonideal { } /***********************************************************************/ - // 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 || (stateCalc != m_vcsStateStatus)) { - if (stateCalc == VCS_STATECALC_OLD || stateCalc == VCS_STATECALC_NEW) { - if (m_owningSolverObject) { - setMolesFromVCS(stateCalc); - } - } - } - } - /***********************************************************************/ // Molar volume calculation for standard state of one species /* @@ -971,7 +972,16 @@ namespace VCSnonideal { * j = id of the species mole number * k = id of the species activity coefficient */ - void vcs_VolPhase::sendToVCS_LnActCoeffJac(double * const * const LnACJac_VCS) const { + void vcs_VolPhase::sendToVCS_LnActCoeffJac(double * const * const LnACJac_VCS) { + /* + * update the Ln Act Coeff jacobian entries with respect to the + * mole number of species in the phase -> we always assume that + * they are out of date. + */ + updateLnActCoeffJac(); + /* + * Now copy over the values + */ 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 f36970b45..29d6a8fb0 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -175,7 +175,8 @@ namespace VCSnonideal { * to gather the species into the local contiguous vector * format. */ - void setMolesFromVCS(const int stateCalc, const double * const molesSpeciesVCS = 0); + void setMolesFromVCS(const int stateCalc, + const double * const molesSpeciesVCS = 0); //! Set the moles within the phase /*! @@ -199,8 +200,21 @@ namespace VCSnonideal { */ void setMolesFromVCSCheck(const int stateCalc, const double * molesSpeciesVCS, - const double * const TPhMoles, - int iphase = -1); + const double * const TPhMoles); + + //! Update the moles within the phase, if necessary + /*! + * This function takes as input the stateCalc value, which + * determines where within VCS_SOLVE to fetch the mole numbers. + * It then updates this object with their values. This is essentially + * a gather routine. + * + * @param stateCalc State calc value either VCS_STATECALC_OLD + * or VCS_STATECALC_NEW. With any other value + * nothing is done. + * + */ + void updateFromVCS_MoleNumbers(const int stateCalc); //! Fill in an activity coefficients vector within a VCS_SOLVE object /*! @@ -271,20 +285,7 @@ namespace VCSnonideal { */ 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 /*! @@ -361,6 +362,7 @@ namespace VCSnonideal { */ void setState_TP(double temperature_Kelvin, double pressure_PA); + private: //! Evaluation of Activity Coefficient Jacobians /*! * This is the derivative of the ln of the activity coefficient @@ -376,6 +378,7 @@ namespace VCSnonideal { */ void updateLnActCoeffJac(); + public: // Downloads the ln ActCoeff jacobian into the VCS version of the // ln ActCoeff jacobian. /* @@ -389,7 +392,7 @@ namespace VCSnonideal { * j = id of the species mole number * k = id of the species activity coefficient */ - void sendToVCS_LnActCoeffJac(double * const * const LnACJac_VCS) const; + void sendToVCS_LnActCoeffJac(double * const * const LnACJac_VCS); //! Set the pointer for Cantera's ThermoPhase parameter /*! @@ -736,7 +739,6 @@ namespace VCSnonideal { */ mutable std::vector ActCoeff; - //! Vector of the derivatives of the ln activity coefficient wrt to the //! current mole number /*! diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index feeea0d35..74cfae72b 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -325,6 +325,7 @@ namespace VCSnonideal { /* ******************************************* */ /* **** CONVERGENCE FORCING SECTION ********** */ /* ******************************************* */ + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, nspecies); for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) { s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec]; diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 08fa14528..96edd0725 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -86,6 +86,7 @@ namespace VCSnonideal { if (m_unitsState == VCS_DIMENSIONAL_G) { vcs_nondim_TP(); } + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); /* ******************************************************** */ /* *** PRINT OUT RESULTS ********************************** */ diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 4593feec5..8e33bc259 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -324,16 +324,10 @@ 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(); /* * Download the resulting calculation into the full vector * -> This scatter calculation is carried out in the diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 5d65034e1..0ddf84cb4 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -1323,7 +1323,12 @@ private: void vcs_TCounters_report(int timing_print_lvl = 1); - void vcs_setMoleNumVolPhases(bool upToDate, int stateCalc); + void vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc); + + void vcs_setFlagsVolPhase(const int iph, const bool upToDate, const int stateCalc); + + void vcs_forceMolUpdateVolPhase(const int stateCalc); + public: //! value of the number of species used to malloc data structures diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 4931c59f8..a3ae00397 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -284,6 +284,7 @@ namespace VCSnonideal { /* ***************************************************************************** */ /* **** EVALUATE ALL CHEMICAL POTENTIALS AT THE OLD (CURRENT) MOLE NUMBERS ***** */ /* ***************************************************************************** */ + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); /* @@ -382,6 +383,7 @@ namespace VCSnonideal { } #endif vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx)); + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); } #ifdef DEBUG_MODE @@ -419,6 +421,7 @@ namespace VCSnonideal { * We have already evaluated the major non-components */ if (uptodate_minors == FALSE) { + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); vcs_deltag(1, false, VCS_STATECALC_NEW); } @@ -608,7 +611,8 @@ namespace VCSnonideal { double maxPermissible = m_elemAbundancesGoal[j] / atomComp; if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { #ifdef DEBUG_MODE - sprintf(ANOTE, "Species stays zeroed even though dG neg, because of %s elemAbund", + sprintf(ANOTE, "Species stays zeroed even though dG " + "neg, because of %s elemAbund", m_elementName[j].c_str()); #endif resurrect = false; @@ -1099,7 +1103,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_setFlagsVolPhases(false, VCS_STATECALC_NEW); vcs_updateVP(VCS_STATECALC_NEW); vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); @@ -1241,6 +1245,7 @@ namespace VCSnonideal { * we have already done this inside the FORCED * loop. */ + vcs_forceMolUpdateVolPhase(VCS_STATECALC_NEW); vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_old), VCS_DATA_PTR(m_tPhaseMoles_new), m_numPhases); vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new), m_numSpeciesRdc); @@ -1249,7 +1254,8 @@ namespace VCSnonideal { vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc); vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_new), m_numSpeciesRdc); - vcs_updateVP(VCS_STATECALC_OLD); + //vcs_updateVP(VCS_STATECALC_OLD); + vcs_setFlagsVolPhases(true, VCS_STATECALC_OLD); /* * Increment the iteration counters */ @@ -1325,6 +1331,7 @@ namespace VCSnonideal { VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, &usedZeroedSpecies); if (retn != VCS_SUCCESS) return retn; + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); vcs_deltag(0, true, VCS_STATECALC_OLD); uptodate_minors = TRUE; @@ -1358,6 +1365,7 @@ namespace VCSnonideal { } #endif vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx)); + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); vcs_deltag(0, true, VCS_STATECALC_OLD); uptodate_minors = TRUE; @@ -1507,7 +1515,8 @@ namespace VCSnonideal { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { plogf(" --- Get a new basis because %s", m_speciesName[l].c_str()); - plogf(" has dg < 0.0 and comp %s has zero mole num", m_speciesName[j].c_str()); + plogf(" has dg < 0.0 and comp %s has zero mole num", + m_speciesName[j].c_str()); plogf(" and share nonzero stoic: %-9.1f", m_stoichCoeffRxnMatrix[i][j]); plogendl(); @@ -1593,6 +1602,7 @@ namespace VCSnonideal { * For this special case, we must reevaluate thermo functions */ if (iti != 0) { + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, kspec, kspec+1); vcs_deltag(0, false, VCS_STATECALC_OLD); } @@ -1674,6 +1684,7 @@ namespace VCSnonideal { * for minor species, if needed. */ if (iti != 0) { + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); vcs_deltag(1, false, VCS_STATECALC_OLD); uptodate_minors = TRUE; @@ -1784,6 +1795,7 @@ namespace VCSnonideal { /* * Go back to evaluate the total moles of gas and liquid. */ + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); vcs_deltag(0, false, VCS_STATECALC_OLD); /* @@ -1874,6 +1886,7 @@ namespace VCSnonideal { * for minor species and go back to do a full iteration */ MajorSpeciesHaveConverged = true; + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); vcs_deltag(0, false, VCS_STATECALC_OLD); iti = 0; @@ -1893,6 +1906,7 @@ namespace VCSnonideal { * for minor species and go back to do a full iteration */ MajorSpeciesHaveConverged = true; + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); vcs_deltag(0, false, VCS_STATECALC_OLD); iti = 0; @@ -2505,7 +2519,7 @@ namespace VCSnonideal { */ Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, VCS_DATA_PTR(m_molNumSpecies_old), - VCS_DATA_PTR(m_tPhaseMoles_old), iph); + VCS_DATA_PTR(m_tPhaseMoles_old)); } /**********************************************************************************/ @@ -2757,7 +2771,7 @@ namespace VCSnonideal { #endif } } - + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); vcs_deltag(0, true, VCS_STATECALC_OLD); @@ -2767,12 +2781,14 @@ namespace VCSnonideal { iph = m_phaseID[kspec]; if (m_tPhaseMoles_old[iph] > 0.0) { if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) { - if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) > VCS_DELETE_MINORSPECIES_CUTOFF) || + if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) > + VCS_DELETE_MINORSPECIES_CUTOFF) || (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) { retn++; #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", + plogf(" --- add_deleted(): species %s " + "with mol number %g not converged: DG = %g", m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], m_deltaGRxn_old[irxn]); plogendl(); @@ -2920,6 +2936,7 @@ namespace VCSnonideal { * only step is being carried out, then we don't need to * update the minor noncomponents. */ + vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesRdc); /* @@ -4937,11 +4954,11 @@ namespace VCSnonideal { if (vcsState == VCS_STATECALC_OLD) { Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, VCS_DATA_PTR(m_molNumSpecies_old), - VCS_DATA_PTR(m_tPhaseMoles_old), i); + VCS_DATA_PTR(m_tPhaseMoles_old)); } else if (vcsState == VCS_STATECALC_NEW) { Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW, VCS_DATA_PTR(m_molNumSpecies_new), - VCS_DATA_PTR(m_tPhaseMoles_new), i); + VCS_DATA_PTR(m_tPhaseMoles_new)); } #ifdef DEBUG_MODE else { @@ -5362,8 +5379,7 @@ namespace VCSnonideal { } /*******************************************************************************/ - - void VCS_SOLVE::vcs_setMoleNumVolPhases(bool upToDate, int stateCalc) { + void VCS_SOLVE::vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc) { int iph; vcs_VolPhase *Vphase; if (!upToDate) { @@ -5379,6 +5395,29 @@ namespace VCSnonideal { } } } + /*******************************************************************************/ + void VCS_SOLVE::vcs_setFlagsVolPhase(const int iph, const bool upToDate, + const int stateCalc) { + vcs_VolPhase *Vphase; + if (!upToDate) { + Vphase = m_VolPhaseList[iph]; + Vphase->m_UpToDate = false; + } else { + Vphase = m_VolPhaseList[iph]; + Vphase->m_UpToDate = true; + Vphase->m_vcsStateStatus = stateCalc; + } + } + /*******************************************************************************/ + + void VCS_SOLVE::vcs_forceMolUpdateVolPhase(const int stateCalc) { + int iph; + vcs_VolPhase *Vphase; + for (iph = 0; iph < m_numPhases; iph++) { + Vphase = m_VolPhaseList[iph]; + Vphase->updateFromVCS_MoleNumbers(stateCalc); + } + } }