diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index 9d8863917..b142a3078 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -1280,7 +1280,7 @@ namespace VCSnonideal { vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; std::string sEOS = string16_EOSType(VolPhase->m_eqnState); plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), - VolPhase->VP_ID, VolPhase->m_singleSpecies, + VolPhase->VP_ID_, VolPhase->m_singleSpecies, VolPhase->m_gasPhase, sEOS.c_str(), VolPhase->nSpecies(), VolPhase->totalMolesInert() ); plogf("%16e\n", VolPhase->totalMoles()); @@ -1402,7 +1402,7 @@ namespace VCSnonideal { vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase]; std::string sEOS = string16_EOSType(VolPhase->m_eqnState); plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(), - VolPhase->VP_ID, VolPhase->m_singleSpecies, + VolPhase->VP_ID_, VolPhase->m_singleSpecies, VolPhase->m_gasPhase, sEOS.c_str(), VolPhase->nSpecies(), VolPhase->totalMolesInert() ); plogf("%16e\n", VolPhase->totalMoles() ); diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 760f919ff..47a3f1d7c 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -32,7 +32,7 @@ namespace VCSnonideal { */ vcs_VolPhase::vcs_VolPhase(VCS_SOLVE * owningSolverObject) : m_owningSolverObject(0), - VP_ID(-1), + VP_ID_(-1), Domain_ID(-1), m_singleSpecies(true), m_gasPhase(false), @@ -94,7 +94,7 @@ namespace VCSnonideal { */ vcs_VolPhase::vcs_VolPhase(const vcs_VolPhase& b) : m_owningSolverObject(b.m_owningSolverObject), - VP_ID(b.VP_ID), + VP_ID_(b.VP_ID_), Domain_ID(b.Domain_ID), m_singleSpecies(b.m_singleSpecies), m_gasPhase(b.m_gasPhase), @@ -147,7 +147,7 @@ namespace VCSnonideal { // operator but is true for a copy constructor // m_owningSolverObject = b.m_owningSolverObject; - VP_ID = b.VP_ID; + VP_ID_ = b.VP_ID_; Domain_ID = b.Domain_ID; m_singleSpecies = b.m_singleSpecies; m_gasPhase = b.m_gasPhase; @@ -261,17 +261,17 @@ namespace VCSnonideal { m_phi = 0.0; m_phiVarIndex = -1; - if (phaseNum == VP_ID) { + if (phaseNum == VP_ID_) { if (strcmp(PhaseName.c_str(), phaseName)) { plogf("Strings are different: %s %s :unknown situation\n", PhaseName.c_str(), phaseName); exit(EXIT_FAILURE); } } else { - VP_ID = phaseNum; + VP_ID_ = phaseNum; if (!phaseName) { char itmp[40]; - sprintf(itmp, "Phase_%d", VP_ID); + sprintf(itmp, "Phase_%d", VP_ID_); PhaseName = itmp; } else { PhaseName = phaseName; @@ -711,7 +711,7 @@ namespace VCSnonideal { /* * Check for consistency with TPhMoles[] */ - double Tcheck = TPhMoles[VP_ID]; + double Tcheck = TPhMoles[VP_ID_]; if (Tcheck != v_totalMoles) { if (vcs_doubleEqual(Tcheck, v_totalMoles)) { Tcheck = v_totalMoles; @@ -1104,7 +1104,7 @@ namespace VCSnonideal { if (m_numSpecies != 0) { plogf("Warning Nsp != NVolSpeces: %d %d \n", nsp, m_numSpecies); } - resize(VP_ID, nsp, nelem, PhaseName.c_str()); + resize(VP_ID_, nsp, nelem, PhaseName.c_str()); } TP_ptr->getMoleFractions(VCS_DATA_PTR(Xmol)); fractionCreationDelta_ = Xmol; @@ -1567,7 +1567,7 @@ namespace VCSnonideal { std::string pname = tPhase->id(); if (pname == "") { char sss[50]; - sprintf(sss, "phase%d", VP_ID); + sprintf(sss, "phase%d", VP_ID_); pname = sss; } ename = "cn_" + pname; diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 808558a16..2ea6349c2 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -673,7 +673,7 @@ namespace VCSnonideal { * miscibility gap, these numbers will stay the * same after the split. */ - int VP_ID; + int VP_ID_; //! ID of the surface or volume domain in which the //! this phase exists diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp index e92dd7f54..3a370809b 100644 --- a/Cantera/src/equil/vcs_phaseStability.cpp +++ b/Cantera/src/equil/vcs_phaseStability.cpp @@ -1,8 +1,17 @@ -/* ======================================================================= */ -/* $RCSfile: vcs_phaseStability.cpp,v $ */ -/* $Date$ */ -/* $Revision$ */ -/* ======================================================================= */ +/** + * @file vcs_phaseStability.cpp + * Implementation class for functions associated with determining the stability of a phase + * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). + */ + +/* + * $Id$ + */ +/* + * Copywrite (2005) Sandia Corporation. Under the terms of + * Contract DE-AC04-94AL85000 with Sandia Corporation, the + * U.S. Government retains certain rights in this software. + */ #include "vcs_solve.h" @@ -138,7 +147,7 @@ namespace VCSnonideal { * @return returns the phase id of the phase that pops back into * existence. Returns -1 if there are no phases */ - int VCS_SOLVE::vcs_popPhaseID() { + int VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) { int iphasePop = -1; int iph; int irxn, kspec; @@ -151,8 +160,8 @@ namespace VCSnonideal { char anote[128]; if (m_debug_print_lvl >= 2) { plogf(" --- vcs_popPhaseID() called\n"); - plogf(" --- Phase Status F_e MoleNum\n"); - plogf(" --------------------------------------------------------------\n"); + plogf(" --- Phase Status F_e MoleNum\n"); + plogf(" --------------------------------------------------------------------------\n"); } #endif for (iph = 0; iph < m_numPhases; iph++) { @@ -165,7 +174,7 @@ namespace VCSnonideal { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - plogf(" --- %18s %5d NA %11.3e\n", + plogf(" --- %18s %5d NA %11.3e\n", Vphase->PhaseName.c_str(), existence, m_tPhaseMoles_old[iph]); @@ -194,10 +203,19 @@ namespace VCSnonideal { #endif } } +#ifdef DEBUG_MODE + if (Fephase < 0.0) { + strcpy(anote," (not stable)"); + if (m_tPhaseMoles_old[iph] > 0.0) { + printf("shouldn't be here\n"); + exit(-1); + } + } +#endif #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - plogf(" --- %18s %5d NA %11.3g %s\n", + plogf(" --- %18s %5d %10.3g %10.3g %s\n", Vphase->PhaseName.c_str(), existence, Fephase, m_tPhaseMoles_old[iph], anote); @@ -242,11 +260,19 @@ namespace VCSnonideal { } } } + phasePopPhaseIDs.resize(0); + if (iphasePop >= 0) { + phasePopPhaseIDs.push_back(iphasePop); + } + /* + * Insert logic here to figure out if phase pops are linked together. Only do one linked + * pop at a time. + */ #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - plogf(" --------------------------------------------------------------\n"); + plogf(" ---------------------------------------------------------------------\n"); } #endif return iphasePop; @@ -280,6 +306,7 @@ namespace VCSnonideal { int kspec = Vphase->spGlobalIndexVCS(0); // Identify the formation reaction for that species int irxn = kspec - m_numComponents; + std::vector creationGlobalRxnNumbers; doublereal s; int j, k; @@ -378,8 +405,10 @@ namespace VCSnonideal { } else { vector fracDelta(Vphase->nSpecies()); vector X_est(Vphase->nSpecies()); + + // fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); fracDelta = Vphase->fractionCreationDeltas(); - + double sumFrac = 0.0; for (k = 0; k < Vphase->nSpecies(); k++) { sumFrac += fracDelta[k]; @@ -490,6 +519,8 @@ namespace VCSnonideal { vector fracDelta_new(Vphase->nSpecies(), 0.0); vector fracDelta_old(Vphase->nSpecies(), 0.0); vector fracDelta_raw(Vphase->nSpecies(), 0.0); + vector creationGlobalRxnNumbers(Vphase->nSpecies(), -1); + vector m_feSpecies_Deficient(m_numComponents, 0.0); doublereal damp = 1.0; @@ -505,8 +536,8 @@ namespace VCSnonideal { // Get the storred estimate for the composition of the phase if // it gets created + // fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); fracDelta_new = Vphase->fractionCreationDeltas(); - bool oneIsComponent = false; std::vector componentList; @@ -539,6 +570,12 @@ namespace VCSnonideal { plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph); } #endif + + for (k = 0; k < Vphase->nSpecies(); k++) { + if (fracDelta_new[k] < 1.0E-13) { + fracDelta_new[k] = 1.0E-13; + } + } bool converged = false; for (int its = 0; its < 200 && (!converged); its++) { @@ -547,6 +584,8 @@ namespace VCSnonideal { fracDelta_old = fracDelta_new; dirProdOld = dirProd; + + // Given a set of fracDelta's, we calculate the fracDelta's // for the component species, if any for (i = 0; i < (int) componentList.size(); i++) { @@ -753,7 +792,9 @@ namespace VCSnonideal { if (converged) { Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), VCS_STATECALC_PHASESTABILITY); - Vphase->setFractionCreationDeltas( VCS_DATA_PTR(fracDelta_new)); + // Vphase->setCreationMoleNumbers(VCS_DATA_PTR(fracDelta_new), creationGlobalRxnNumbers); + Vphase->setFractionCreationDeltas(VCS_DATA_PTR(fracDelta_new)); + } diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index 6d1891b4c..7fe4d6542 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -284,7 +284,7 @@ namespace VCSnonideal { Vphase = VPhaseList[iphase]; std::string EOS_cstr = string16_EOSType(Vphase->m_eqnState); plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(), - Vphase->VP_ID, Vphase->m_singleSpecies, Vphase->m_gasPhase); + Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase); plogf("%16s %8d %16e ", EOS_cstr.c_str(), Vphase->nSpecies(), Vphase->totalMolesInert()); if (iest >= 0) plogf("%16e\n", Vphase->totalMoles()); diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 108d1dcbc..6331ef942 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -869,9 +869,9 @@ namespace VCSnonideal { vcs_VolPhase *vPhase = m_VolPhaseList[iph]; vcs_VolPhase *pub_phase_ptr = pub->VPhaseList[iph]; - if (vPhase->VP_ID != pub_phase_ptr->VP_ID) { + if (vPhase->VP_ID_ != pub_phase_ptr->VP_ID_) { plogf("%sPhase numbers have changed:%d %d\n", yo.c_str(), - vPhase->VP_ID, pub_phase_ptr->VP_ID); + vPhase->VP_ID_, pub_phase_ptr->VP_ID_); retn = VCS_PUB_BAD; } diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 053e681c8..6100de9e5 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -1,6 +1,7 @@ /** * @file vcs_solve.h - * Header file for the internal object that holds the problem + * Header file for the internal object that holds the vcs equilibrium problem + * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). */ /* * $Id$ @@ -501,7 +502,6 @@ public: */ void vcs_updateVP(const int stateCalc); - //! Utility function that evaluates whether a phase can be popped //! into existence /*! @@ -514,10 +514,13 @@ public: //! Decision as to whether a phase pops back into existence /*! + * @param phasePopPhaseIDs Vector containing the phase ids of the phases + * that will be popped this step. + * * @return returns the phase id of the phase that pops back into * existence. Returns -1 if there are no phases */ - int vcs_popPhaseID(); + int vcs_popPhaseID(std::vector &phasePopPhaseIDs); //! Calculates the deltas of the reactions due to phases popping //! into existence diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 5d8bae98b..d6de12ad7 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -122,6 +122,7 @@ namespace VCSnonideal { double *dnPhase_irxn; double atomComp; int iphasePop; + std::vector phasePopPhaseIDs(0); #ifdef DEBUG_MODE char ANOTE[128]; /* @@ -399,7 +400,8 @@ namespace VCSnonideal { * First step is a major branch in the algorithm. * We first determine if a phase pops into existence. */ - iphasePop = vcs_popPhaseID(); + phasePopPhaseIDs.clear(); + iphasePop = vcs_popPhaseID(phasePopPhaseIDs); /* * */