diff --git a/Cantera/src/equil/MultiPhase.cpp b/Cantera/src/equil/MultiPhase.cpp index f746a035d..31f024995 100644 --- a/Cantera/src/equil/MultiPhase.cpp +++ b/Cantera/src/equil/MultiPhase.cpp @@ -433,8 +433,12 @@ namespace Cantera { } m_moles[ip] = phasemoles; if (nsp > 1) { + if (phasemoles > 0.0) { p->setState_TPX(m_temp, m_press, n + loc); p->getMoleFractions(DATA_PTR(m_moleFractions) + loc); + } else { + p->getMoleFractions(DATA_PTR(m_moleFractions) + loc); + } } else { m_moleFractions[loc] = 1.0; diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index 6359debf0..f578b74de 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -887,177 +887,7 @@ namespace VCSnonideal { fclose(FP); } - /*! - * This function decides whether a phase has charged species - * or not. - */ - static bool hasChargedSpecies(Cantera::ThermoPhase *tPhase) { - int nSpPhase = tPhase->nSpecies(); - for (int k = 0; k < nSpPhase; k++) { - if (tPhase->charge(k) != 0.0) { - return true; - } - } - return false; - } - - /********************************************************************** - * - * chargeNeutralityElement(): - * - * This utility routine decides whether a Cantera ThermoPhase needs - * a constraint equation representing the charge neutrality of the - * phase. It does this by searching for charged species. If it - * finds one, and if the phase needs one, then it returns true. - */ - static bool chargeNeutralityElement(Cantera::ThermoPhase *tPhase) { - int hasCharge = hasChargedSpecies(tPhase); - if (tPhase->chargeNeutralityNecessary()) { - if (hasCharge) { - return true; - } - } - return false; - } - - /************************************************************************* - * - * - * - */ - static int setUpElements(vcs_VolPhase *VolPhase, Cantera::ThermoPhase *tPhase) { - int e, k, eT; - string ename; - int eFound = -2; - /* - * - */ - int nebase = tPhase->nElements(); - int ne = nebase; - int ns = tPhase->nSpecies(); - - /* - * Decide whether we need an extra element constraint for charge - * neutrality of the phase - */ - bool cne = chargeNeutralityElement(tPhase); - if (cne) { - VolPhase->ChargeNeutralityElement = ne; - ne++; - } - - /* - * Assign and malloc structures - */ - VolPhase->elemResize(ne); - - - if (VolPhase->ChargeNeutralityElement >= 0) { - VolPhase->m_elType[VolPhase->ChargeNeutralityElement] = - VCS_ELEM_TYPE_CHARGENEUTRALITY; - } - - if (hasChargedSpecies(tPhase)) { - if (cne) { - /* - * We need a charge neutrality constraint. - * We also have an Electron Element. These are - * duplicates of each other. To avoid trouble with - * possible range error conflicts, sometimes we eliminate - * the Electron condition. Flag that condition for elimination - * by toggling the ElActive variable. If we find we need it - * later, we will retoggle ElActive to true. - */ - for (eT = 0; eT < nebase; eT++) { - ename = tPhase->elementName(eT); - if (ename == "E") { - eFound = eT; - VolPhase->ElActive[eT] = 0; - VolPhase->m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE; - } - } - } else { - for (eT = 0; eT < nebase; eT++) { - ename = tPhase->elementName(eT); - if (ename == "E") { - eFound = eT; - VolPhase->m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE; - } - } - } - if (eFound == -2) { - eFound = ne; - VolPhase->m_elType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE; - VolPhase->ElActive[ne] = 0; - string ename = "E"; - VolPhase->ElName[ne] = ename; - ne++; - VolPhase->elemResize(ne); - } - - } - - VolPhase->FormulaMatrix.resize(ne, ns, 0.0); - - VolPhase->m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM); - - VolPhase->elemResize(ne); - //ElGlobalIndex.resize(ne, -1); - - - e = 0; - for (eT = 0; eT < nebase; eT++) { - ename = tPhase->elementName(eT); - VolPhase->ElName[e] = ename; - e++; - } - - if (cne) { - string pname = tPhase->id(); - if (pname == "") { - char sss[50]; - sprintf(sss, "phase%d", VolPhase->VP_ID); - pname = sss; - } - ename = "cn_" + pname; - e = VolPhase->ChargeNeutralityElement; - VolPhase->ElName[e] = ename; - } - double * const * const fm = VolPhase->FormulaMatrix.baseDataAddr(); - for (k = 0; k < ns; k++) { - e = 0; - for (eT = 0; eT < nebase; eT++) { - fm[e][k] = tPhase->nAtoms(k, eT); - e++; - } - if (eFound >= 0) { - fm[eFound][k] = - tPhase->charge(k); - } - } - - if (cne) { - for (k = 0; k < ns; k++) { - fm[VolPhase->ChargeNeutralityElement][k] = tPhase->charge(k); - } - } - - - /* - * Here, we figure out what is the species types are - * The logic isn't set in stone, and is just for a particular type - * of problem that I'm solving first. - */ - if (ns == 1) { - if (tPhase->charge(0) != 0.0) { - VolPhase->m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE; - VolPhase->setPhiVarIndex(0); - } - } - - return ne; - } - static void print_char(const char letter, const int num) { for (int i = 0; i < num; i++) plogf("%c", letter); } @@ -1207,7 +1037,7 @@ namespace VCSnonideal { * element in the phase to enforce a charge neutrality * constraint. */ - setUpElements(VolPhase, tPhase); + VolPhase->transferElementsFM(tPhase); /* * Combine the element information in the vcs_VolPhase diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 206fe7b57..c9ca477c5 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -18,6 +18,7 @@ #include "ThermoPhase.h" #include "mix_defs.h" +#include #include #include @@ -1322,7 +1323,8 @@ namespace VCSnonideal { //! Returns the global index of the local element index for the phase void vcs_VolPhase::setElemGlobalIndex(const int eLocal, const int eGlobal) { DebugAssertThrowVCS(eLocal >= 0, "vcs_VolPhase::setElemGlobalIndex"); - DebugAssertThrowVCS(eLocal < m_numElemConstraints, "vcs_VolPhase::setElemGlobalIndex"); + DebugAssertThrowVCS(eLocal < m_numElemConstraints, + "vcs_VolPhase::setElemGlobalIndex"); m_elemGlobalIndex[eLocal] = eGlobal; } @@ -1330,5 +1332,173 @@ namespace VCSnonideal { return m_numElemConstraints; } + std::string vcs_VolPhase::elementName(const int e) const { + return ElName[e]; + } + + /*! + * This function decides whether a phase has charged species + * or not. + */ + static bool hasChargedSpecies(const Cantera::ThermoPhase * const tPhase) { + int nSpPhase = tPhase->nSpecies(); + for (int k = 0; k < nSpPhase; k++) { + if (tPhase->charge(k) != 0.0) { + return true; + } + } + return false; + } + /********************************************************************** + * + * chargeNeutralityElement(): + * + * This utility routine decides whether a Cantera ThermoPhase needs + * a constraint equation representing the charge neutrality of the + * phase. It does this by searching for charged species. If it + * finds one, and if the phase needs one, then it returns true. + */ + static bool chargeNeutralityElement(const Cantera::ThermoPhase * const tPhase) { + int hasCharge = hasChargedSpecies(tPhase); + if (tPhase->chargeNeutralityNecessary()) { + if (hasCharge) { + return true; + } + } + return false; + } + + int vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase * const tPhase) { + int e, k, eT; + std::string ename; + int eFound = -2; + /* + * + */ + int nebase = tPhase->nElements(); + int ne = nebase; + int ns = tPhase->nSpecies(); + + /* + * Decide whether we need an extra element constraint for charge + * neutrality of the phase + */ + bool cne = chargeNeutralityElement(tPhase); + if (cne) { + ChargeNeutralityElement = ne; + ne++; + } + + /* + * Assign and malloc structures + */ + elemResize(ne); + + + if (ChargeNeutralityElement >= 0) { + m_elType[ChargeNeutralityElement] = + VCS_ELEM_TYPE_CHARGENEUTRALITY; + } + + if (hasChargedSpecies(tPhase)) { + if (cne) { + /* + * We need a charge neutrality constraint. + * We also have an Electron Element. These are + * duplicates of each other. To avoid trouble with + * possible range error conflicts, sometimes we eliminate + * the Electron condition. Flag that condition for elimination + * by toggling the ElActive variable. If we find we need it + * later, we will retoggle ElActive to true. + */ + for (eT = 0; eT < nebase; eT++) { + ename = tPhase->elementName(eT); + if (ename == "E") { + eFound = eT; + ElActive[eT] = 0; + m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE; + } + } + } else { + for (eT = 0; eT < nebase; eT++) { + ename = tPhase->elementName(eT); + if (ename == "E") { + eFound = eT; + m_elType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE; + } + } + } + if (eFound == -2) { + eFound = ne; + m_elType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE; + ElActive[ne] = 0; + std::string ename = "E"; + ElName[ne] = ename; + ne++; + elemResize(ne); + } + + } + + FormulaMatrix.resize(ne, ns, 0.0); + + m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM); + + elemResize(ne); + //ElGlobalIndex.resize(ne, -1); + + + e = 0; + for (eT = 0; eT < nebase; eT++) { + ename = tPhase->elementName(eT); + ElName[e] = ename; + e++; + } + + if (cne) { + std::string pname = tPhase->id(); + if (pname == "") { + char sss[50]; + sprintf(sss, "phase%d", VP_ID); + pname = sss; + } + ename = "cn_" + pname; + e = ChargeNeutralityElement; + ElName[e] = ename; + } + + double * const * const fm = FormulaMatrix.baseDataAddr(); + for (k = 0; k < ns; k++) { + e = 0; + for (eT = 0; eT < nebase; eT++) { + fm[e][k] = tPhase->nAtoms(k, eT); + e++; + } + if (eFound >= 0) { + fm[eFound][k] = - tPhase->charge(k); + } + } + + if (cne) { + for (k = 0; k < ns; k++) { + fm[ChargeNeutralityElement][k] = tPhase->charge(k); + } + } + + + /* + * Here, we figure out what is the species types are + * The logic isn't set in stone, and is just for a particular type + * of problem that I'm solving first. + */ + if (ns == 1) { + if (tPhase->charge(0) != 0.0) { + m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE; + setPhiVarIndex(0); + } + } + + return ne; + } } diff --git a/Cantera/src/equil/vcs_VolPhase.h b/Cantera/src/equil/vcs_VolPhase.h index 6878a6be5..66032676b 100644 --- a/Cantera/src/equil/vcs_VolPhase.h +++ b/Cantera/src/equil/vcs_VolPhase.h @@ -480,6 +480,22 @@ namespace VCSnonideal { int nElemConstraints() const; + std::string elementName(const int e) const; + + + + //! Transfer all of the element information from the + //! ThermoPhase object to the vcs_VolPhase object. + /*! + * Also decide whether we need a new charge neutrality + * element in the phase to enforce a charge neutrality + * constraint. + * + * @param tPhase Pointer to the thermophase object + */ + int transferElementsFM(const Cantera::ThermoPhase * const tPhase); + + private: //! Evaluate the activity coefficients at the current conditions @@ -549,11 +565,12 @@ namespace VCSnonideal { */ void _updateMoleFractionDependencies(); + /************************************************************************* * MEMBER DATA * ************************************************************************/ - public: + private: //! Backtrack value of VCS_SOLVE * /*! * Note the default for this is 0. That's a valid value too, since @@ -561,6 +578,7 @@ namespace VCSnonideal { */ VCS_SOLVE *m_owningSolverObject; + public: //! Original ID of the phase in the problem. /*! * If a non-ideal phase splits into two due to a @@ -610,12 +628,13 @@ namespace VCSnonideal { */ int ChargeNeutralityElement; + private: //! vector of strings containing the element names /*! * Length = nElemConstraints */ std::vector ElName; - + public: //! boolean indicating whether an element constraint is active //! for the current problem std::vector ElActive; diff --git a/Cantera/src/equil/vcs_prob.cpp b/Cantera/src/equil/vcs_prob.cpp index 3222dc1a5..8408c3807 100644 --- a/Cantera/src/equil/vcs_prob.cpp +++ b/Cantera/src/equil/vcs_prob.cpp @@ -369,7 +369,7 @@ namespace VCSnonideal { */ for (eVP = 0; eVP < neVP; eVP++) { foundPos = -1; - enVP = volPhase->ElName[eVP]; + enVP = volPhase->elementName(eVP); /* * Search for matches with the existing elements. * If found, then fill in the entry in the global diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 97ae896de..58e9d897d 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -1059,8 +1059,10 @@ private: * loop. * * @param iph Phase to be deleted + * + * @return Returns whether the operation was successful or not */ - void vcs_delete_multiphase(const int iph); + bool vcs_delete_multiphase(const int iph); //! Change the concentration of a species by delta moles. /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 04d38a174..e518846c4 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -2116,8 +2116,8 @@ namespace VCSnonideal { * * * @param kspec The species index - * @delta_ptr pointer to the delta for the species. This may change during - * the calculation + * @param delta_ptr pointer to the delta for the species. This may change during + * the calculation * * @return * 1: succeeded without change of dx @@ -2129,6 +2129,13 @@ namespace VCSnonideal { int j; double tmp; double delta = *delta_ptr; +#ifdef DEBUG_MODE + if (irxn < 0) { + plogf(" --- delete_species() ERROR: called for a component %d", kspec); + plogendl(); + std::exit(-1); + } +#endif if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { /* * Attempt the given dx. If it doesn't work, try to see if a smaller @@ -2406,10 +2413,11 @@ namespace VCSnonideal { * * @param iph Phase to be deleted */ - void VCS_SOLVE::vcs_delete_multiphase(const int iph) { - int kspec, j, irxn; + bool VCS_SOLVE::vcs_delete_multiphase(const int iph) { + int kspec, irxn; double dx; vcs_VolPhase *Vphase = m_VolPhaseList[iph]; + bool successful = true; /* * set the phase existence flag to dead */ @@ -2426,6 +2434,8 @@ namespace VCSnonideal { m_tPhaseMoles_new[iph] = 0.0; m_deltaPhaseMoles[iph] = 0.0; + + /* * Loop over all of the active species in the phase. */ @@ -2433,33 +2443,71 @@ namespace VCSnonideal { if (m_phaseID[kspec] == iph) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { irxn = kspec - m_numComponents; - /* - * calculate an extent of rxn, dx, that zeroes out the species. - */ - dx = - (m_molNumSpecies_old[kspec]); - /* - * Set the mole numbers of that species to zero. - */ - m_molNumSpecies_old[kspec] = 0.0; - m_molNumSpecies_new[kspec] = 0.0; - m_deltaMolNumSpecies[kspec] = 0.0; - /* - * Change the status flag of the species to that of an - * zeroed phase - */ - m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDPHASE; - /* - * Changed the component mole numbers to account for the - * final extent of reaction. Make sure to keep component - * mole numbers above zero - * - */ - for (j = 0; j < m_numComponents; ++j) { - m_molNumSpecies_old[j] += m_stoichCoeffRxnMatrix[irxn][j] * dx; - if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) { - if (m_molNumSpecies_old[j] < 0.0) { - m_molNumSpecies_old[j] = 0.0; + if (irxn >= 0) { + /* + * calculate an extent of rxn, dx, that zeroes out the species. + */ + dx = - (m_molNumSpecies_old[kspec]); + double dxTent = dx; + + int retn = delta_species(kspec, &dxTent); + if (retn != 1) { + successful = false; +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- delete_multiphase %d, %s ERROR problems deleting species %s\n", + iph, Vphase->PhaseName.c_str(), m_speciesName[kspec].c_str() ); + plogf(" --- delta attempted: %g achieved: %g " + " Zeroing it manually\n", dx, dxTent); } +#endif + m_molNumSpecies_old[kspec] = 0.0; + m_molNumSpecies_new[kspec] = 0.0; + m_deltaMolNumSpecies[kspec] = 0.0; + // recover the total phase moles. + vcs_tmoles(); + } else { + /* + * Set the mole number of that species to zero. + */ + m_molNumSpecies_old[kspec] = 0.0; + m_molNumSpecies_new[kspec] = 0.0; + m_deltaMolNumSpecies[kspec] = 0.0; + } + /* + * Change the status flag of the species to that of an + * zeroed phase + */ + m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDPHASE; + /* + * Changed the component mole numbers to account for the + * final extent of reaction. Make sure to keep component + * mole numbers above zero + * + */ + // for (j = 0; j < m_numComponents; ++j) { + // m_molNumSpecies_old[j] += m_stoichCoeffRxnMatrix[irxn][j] * dx; + //if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) { + //if (m_molNumSpecies_old[j] < 0.0) { + // m_molNumSpecies_old[j] = 0.0; + //} + //} + //} + } + else { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n", + kspec, m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec]); + } +#endif + if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF * VCS_DELETE_PHASE_CUTOFF ) { + plogf(" --- delete_multiphase unknown situation error exit"); + plogendl(); + std::exit(-1); + } else { + m_molNumSpecies_old[kspec] = 0.0; + m_molNumSpecies_new[kspec] = 0.0; } } } @@ -2505,6 +2553,7 @@ namespace VCSnonideal { Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old)); + return successful; } /**********************************************************************************/ @@ -4647,7 +4696,7 @@ namespace VCSnonideal { if (! vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph])) { plogf("phase Moles may be off, iph = %d, %20.14g %20.14g \n", iph, tlogMoles[iph], tPhMoles_ptr[iph]); - exit(0); + std::exit(0); } } #endif