diff --git a/Cantera/src/equil/Makefile.in b/Cantera/src/equil/Makefile.in index 363c44348..0cdec0ecd 100644 --- a/Cantera/src/equil/Makefile.in +++ b/Cantera/src/equil/Makefile.in @@ -74,7 +74,7 @@ VCSNONIDEAL_OBJ = vcs_solve_TP.o vcs_VolPhase.o vcs_solve.o vcs_prob.o \ vcs_root1d.o vcs_rxnadj.o \ vcs_SpeciesProperties.o vcs_equilibrate.o \ vcs_prep.o vcs_species_thermo.o vcs_Gibbs.o vcs_phaseStability.o \ - $(DALT_OBJ) + vcs_solve_phaseStability.o $(DALT_OBJ) VCSNONIDEAL_H = vcs_internal.h vcs_VolPhase.h vcs_solve.h vcs_prob.h \ vcs_IntStarStar.h vcs_DoubleStarStar.h vcs_defs.h \ diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index 387a83f1e..f4b0a3b7c 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -42,7 +42,7 @@ using namespace std; //using namespace VCSnonideal; namespace VCSnonideal { - + //==================================================================================================================== vcs_MultiPhaseEquil::vcs_MultiPhaseEquil() : m_vprob(0), @@ -51,7 +51,7 @@ namespace VCSnonideal { m_vsolvePtr(0) { } - + //==================================================================================================================== vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(mix_t* mix, int printLvl) : m_vprob(0), m_mix(0), @@ -89,7 +89,7 @@ namespace VCSnonideal { m_vsolvePtr = 0; } } - + //==================================================================================================================== int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, int estimateEquil, int printLvl, doublereal err, @@ -209,7 +209,7 @@ namespace VCSnonideal { return iSuccess; } - + //==================================================================================================================== int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget, int XY, double Tlow, double Thigh, int estimateEquil, @@ -359,7 +359,7 @@ namespace VCSnonideal { done:; return iSuccess; } - + //==================================================================================================================== int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget, double Tlow, double Thigh, int estimateEquil, @@ -515,7 +515,7 @@ namespace VCSnonideal { done:; return iSuccess; } - + //==================================================================================================================== /* * Equilibrate the solution using the current element abundances @@ -566,7 +566,7 @@ namespace VCSnonideal { } return iSuccess; } - + //==================================================================================================================== /* * Equilibrate the solution using the current element abundances */ @@ -733,7 +733,7 @@ namespace VCSnonideal { } - + //==================================================================================================================== /************************************************************************** * * @@ -905,7 +905,7 @@ namespace VCSnonideal { static void print_char(const char letter, const int num) { for (int i = 0; i < num; i++) plogf("%c", letter); } - + //==================================================================================================================== /* * * @@ -1299,7 +1299,7 @@ namespace VCSnonideal { return VCS_SUCCESS; } - + //==================================================================================================================== // Transfer the current state of mphase into the VCS_PROB object /* * The basic problem has already been set up. @@ -1421,7 +1421,7 @@ namespace VCSnonideal { return VCS_SUCCESS; } - + //==================================================================================================================== // This routine hasn't been checked yet void vcs_MultiPhaseEquil::getStoichVector(index_t rxn, Cantera::vector_fp& nu) { int nsp = m_vsolvePtr->m_numSpeciesTot; @@ -1442,7 +1442,7 @@ namespace VCSnonideal { } } - + //==================================================================================================================== int vcs_MultiPhaseEquil::numComponents() const { int nc = -1; if (m_vsolvePtr) { @@ -1450,7 +1450,7 @@ namespace VCSnonideal { } return nc; } - + //==================================================================================================================== int vcs_MultiPhaseEquil::numElemConstraints() const { int nec = -1; if (m_vsolvePtr) { @@ -1459,11 +1459,185 @@ namespace VCSnonideal { return nec; } - + //==================================================================================================================== int vcs_MultiPhaseEquil::component(int m) const { int nc = numComponents(); if (m < nc) return m_vsolvePtr->m_speciesMapIndex[m]; else return -1; } + //==================================================================================================================== + // Determine the phase stability of a phase at the current conditions + /* + * Equilibration of the solution is not done before the determination is made. + * + * @param iph Phase number to determine the equilibrium. If the phase + * has a non-zero mole number.... + * + * @param funcStab Value of the phase pop function + * + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * + * @param loglevel Determines the amount of printing to the HTML + * output file. + */ + int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double &funcStab, int printLvl, int loglevel) { + + + clockWC tickTock; + int nsp = m_mix->nSpecies(); + int nel = m_mix->nElements(); + int nph = m_mix->nPhases(); + if (m_vprob == 0) { + m_vprob = new VCS_PROB(nsp, nel, nph); + } + m_printLvl = printLvl; + m_vprob->m_printLvl = printLvl; + + /* + * Extract the current state information + * from the MultiPhase object and + * Transfer it to VCS_PROB object. + */ + int res = vcs_Cantera_update_vprob(m_mix, m_vprob); + if (res != 0) { + plogf("problems\n"); + } + + + + // Check obvious bounds on the temperature and pressure + // NOTE, we may want to do more here with the real bounds + // given by the ThermoPhase objects. + double T = m_mix->temperature(); + if (T <= 0.0) { + throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability", + "Temperature less than zero on input"); + } + double pres = m_mix->pressure(); + if (pres <= 0.0) { + throw CanteraError("vcs_MultiPhaseEquil::determine_PhaseStability", + "Pressure less than zero on input"); + } + + beginLogGroup("vcs_MultiPhaseEquil::determine_PhaseStability", loglevel); + addLogEntry("problem type", "fixed T,P"); + addLogEntry("Temperature", T); + addLogEntry("Pressure", pres); + + + /* + * Print out the problem specification from the point of + * view of the vprob object. + */ + m_vprob->prob_report(m_printLvl); + + /* + * Call the thermo Program + */ + int ip1 = m_printLvl; + if (m_printLvl >= 3) { + ip1 = m_printLvl - 2; + } else { + ip1 = 0; + } + if (!m_vsolvePtr) { + m_vsolvePtr = new VCS_SOLVE(); + } + double feStable; + int iStable = m_vsolvePtr->vcs_PS(m_vprob, iph, printLvl, feStable); + + /* + * Transfer the information back to the MultiPhase object. + * Note we don't just call setMoles, because some multispecies + * solution phases may be zeroed out, and that would cause a problem + * for that routine. Also, the mole fractions of such zereod out + * phases actually contain information about likely reemergent + * states. + */ + m_mix->uploadMoleFractionsFromPhases(); + int kGlob = 0; + for (int ip = 0; ip < m_vprob->NPhase; ip++) { + double phaseMole = 0.0; + Cantera::ThermoPhase &tref = m_mix->phase(ip); + int nspPhase = tref.nSpecies(); + for (int k = 0; k < nspPhase; k++, kGlob++) { + phaseMole += m_vprob->w[kGlob]; + } + //phaseMole *= 1.0E-3; + m_mix->setPhaseMoles(ip, phaseMole); + } + + double te = tickTock.secondsWC(); + if (printLvl > 0) { + plogf("\n Results from vcs_PS:\n"); + + plogf("\n"); + plogf("Temperature = %g Kelvin\n", m_vprob->T); + plogf("Pressure = %g Pa\n", m_vprob->PresPA); + std::string sss = m_mix->phaseName(iph); + if (iStable) { + plogf("Phase %d named %s is stable, function value = %g > 0\n", iph, sss.c_str(), feStable); + } else { + plogf("Phase %d named %s is not stable + function value = %g < 0\n", iph, sss.c_str(), feStable); + } + plogf("\n"); + plogf("----------------------------------------" + "---------------------\n"); + plogf(" Name Mole_Number"); + if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) { + plogf("(kmol)"); + } else { + plogf("(gmol)"); + } + plogf(" Mole_Fraction Chem_Potential"); + if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) + plogf(" (kcal/mol)\n"); + else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) + plogf(" (Dimensionless)\n"); + else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KJMOL) + plogf(" (kJ/mol)\n"); + else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_KELVIN) + plogf(" (Kelvin)\n"); + else if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) + plogf(" (J/kmol)\n"); + plogf("-------------------------------------------------------------\n"); + for (int i = 0; i < m_vprob->nspecies; i++) { + plogf("%-12s", m_vprob->SpName[i].c_str()); + if (m_vprob->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + plogf(" %15.3e %15.3e ", 0.0, m_vprob->mf[i]); + plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]); + } else { + plogf(" %15.3e %15.3e ", m_vprob->w[i], m_vprob->mf[i]); + if (m_vprob->w[i] <= 0.0) { + int iph = m_vprob->PhaseID[i]; + vcs_VolPhase *VPhase = m_vprob->VPhaseList[iph]; + if (VPhase->nSpecies() > 1) { + plogf(" -1.000e+300\n"); + } else { + plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]); + } + } else { + plogf("%15.3e\n", m_vprob->m_gibbsSpecies[i]); + } + } + } + plogf("------------------------------------------" + "-------------------\n"); + if (printLvl > 2) { + if (m_vsolvePtr->m_timing_print_lvl > 0) { + plogf("Total time = %12.6e seconds\n", te); + } + } + } + if (loglevel > 0) { + endLogGroup(); + } + return iStable; + } +//==================================================================================================================== + } diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.h b/Cantera/src/equil/vcs_MultiPhaseEquil.h index c090c5a7c..699b6e46a 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.h +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.h @@ -154,7 +154,7 @@ namespace Cantera { * themselves. Two other thermodynamic quantities, determined by the * XY string, are held constant during the equilibration. * - * @param s The object to set to an equilibrium state + * @param s The MultiPhase object to be set to an equilibrium state * * @param ixy An integer specifying the two properties to be held * constant. @@ -205,6 +205,30 @@ namespace Cantera { doublereal rtol = 1.0e-9, int maxsteps = VCS_MAXSTEPS, int maxiter = 100, int loglevel = -99); + //! Determine the phase stability of a single phase given the current conditions + //! in a MultiPhase object + /*! + * + * @param s The MultiPhase object to be set to an equilibrium state + * @param iphase Phase index within the multiphase object to be + * tested for stability. + * @param funcStab Function value that tests equilibrium. > 0 indicates stable + * < 0 indicates unstable + * + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * + * @param loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose + * messages are written as loglevel increases. The + * messages are written to a file in HTML format for viewing + * in a web browser. @see HTML_logs + */ + int vcs_determine_PhaseStability(MultiPhase& s, int iphase, + double &funcStab, int printLvl, int loglevel); + } namespace VCSnonideal { @@ -540,6 +564,22 @@ namespace VCSnonideal { int printLvl = 0, doublereal err = 1.0E-6, int maxsteps = VCS_MAXSTEPS, int loglevel = -99); + //! Determine the phase stability of a phase at the current conditions + /*! + * Equilibration of the solution is not done before the determination is made. + * + * @param iph Phase number to determine the equilibrium. If the phase + * has a non-zero mole number.... + * @param funcStab Value of the phase pop function + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * @param loglevel Determines the amount of printing to the HTML + * output file. + */ + int determine_PhaseStability(int iph, double &funcStab, int printLvl= 0, int logLevel = -99); + //! Report the equilibrium answer in a comma separated table format /*! * This routine is used for in the test suite. @@ -563,6 +603,8 @@ namespace VCSnonideal { */ int numElemConstraints() const; + + // Friend functions friend int vcs_Cantera_to_vprob(Cantera::MultiPhase *mphase, @@ -638,8 +680,7 @@ namespace VCSnonideal { //! Pointer to the object that does all of the equilibration work. /*! * VCS_SOLVE will have different ordering for species and element constraints - * than this object or the VCS_PROB object. - * This object owns the pointer. + * than this object or the VCS_PROB object. This object owns the pointer. */ VCSnonideal::VCS_SOLVE *m_vsolvePtr; diff --git a/Cantera/src/equil/vcs_equilibrate.cpp b/Cantera/src/equil/vcs_equilibrate.cpp index ed46e6c14..03692c6b5 100644 --- a/Cantera/src/equil/vcs_equilibrate.cpp +++ b/Cantera/src/equil/vcs_equilibrate.cpp @@ -327,10 +327,8 @@ namespace Cantera { if (solver == 2) { try { - VCSnonideal::vcs_MultiPhaseEquil *eqsolve = - new VCSnonideal::vcs_MultiPhaseEquil(&s, printLvlSub); - int err = eqsolve->equilibrate(ixy, estimateEquil, printLvlSub, - tol, maxsteps, loglevel); + VCSnonideal::vcs_MultiPhaseEquil *eqsolve = new VCSnonideal::vcs_MultiPhaseEquil(&s, printLvlSub); + int err = eqsolve->equilibrate(ixy, estimateEquil, printLvlSub, tol, maxsteps, loglevel); if (err != 0) { retn = -1; addLogEntry("vcs_equilibrate Error - ", err); @@ -386,4 +384,70 @@ namespace Cantera { } return retn; } + + //==================================================================================================================== + // Determine the phase stability of a single phase given the current conditions + // in a MultiPhase object + /* + * + * @param s The MultiPhase object to be set to an equilibrium state + * @param iphase Phase index within the multiphase object to be + * tested for stability. + * @param funcStab Function value that tests equilibrium. > 0 indicates stable + * < 0 indicates unstable + * + * @param printLvl Determines the amount of printing that + * gets sent to stdout from the vcs package + * (Note, you may have to compile with debug + * flags to get some printing). + * + * @param loglevel Controls amount of diagnostic output. loglevel + * = 0 suppresses diagnostics, and increasingly-verbose + * messages are written as loglevel increases. The + * messages are written to a file in HTML format for viewing + * in a web browser. @see HTML_logs + */ + int vcs_determine_PhaseStability(MultiPhase& s, int iphase, + double &funcStab, int printLvl, int loglevel) + { + int iStab = 0; + static int counter = 0; + beginLogGroup("PhaseStability",loglevel); + addLogEntry("multiphase phase stability function"); + beginLogGroup("arguments"); + addLogEntry("iphase",iphase); + addLogEntry("loglevel",loglevel); + endLogGroup("arguments"); + + int printLvlSub = MAX(0, printLvl-1); + + s.init(); + try { + VCSnonideal::vcs_MultiPhaseEquil *eqsolve = new VCSnonideal::vcs_MultiPhaseEquil(&s, printLvlSub); + iStab = eqsolve->determine_PhaseStability(iphase, funcStab, printLvlSub, loglevel); + if (iStab != 0) { + addLogEntry("Phase is stable - ", iphase); + } else { + addLogEntry("Phase is not stable - ", iphase); + } + endLogGroup("PhaseStability"); + // hard code a csv output file. + if (printLvl > 0) { + string reportFile = "vcs_phaseStability.csv"; + if (counter > 0) { + reportFile = "vcs_phaseStability_" + int2str(counter) + ".csv"; + } + eqsolve->reportCSV(reportFile); + counter++; + } + delete eqsolve; + } + catch (CanteraError &e) { + addLogEntry("Failure.", lastErrorMessage()); + endLogGroup("equilibrate"); + throw e; + } + return iStab; + } + //==================================================================================================================== } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 6331ef942..dc790bb28 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -411,9 +411,7 @@ namespace VCSnonideal { * a 2x2 Newton's method, using loops over vcs_TP() to * calculate the residual and Jacobian) */ - iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->PresPA); - /* * If requested to print anything out, go ahead and do so; diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 6da5b2358..cc47b2061 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -147,6 +147,9 @@ public: */ int vcs_solve_TP(int print_lvl, int printDetails, int maxit); + + int vcs_PS(VCS_PROB *vprob, int iph, int printLvl, double &feStable); + void vcs_reinsert_deleted(int kspec); //! Choose the optimum species basis for the calculations @@ -550,6 +553,7 @@ public: */ int vcs_popPhaseRxnStepSizes(const int iphasePop); + //! Calculates formation reaction step sizes. /*! * This is equation 6.4-16, p. 143 in Smith and Missen. @@ -681,7 +685,8 @@ public: * have. */ double vcs_birthGuess(const int kspec); - + + int vcs_solve_phaseStability(const int iphase, int ifunc, double &funcval, int print_lvl); //! Main program to test whether a deleted phase should be brought //! back into existence diff --git a/Cantera/src/equil/vcs_solve_phaseStability.cpp b/Cantera/src/equil/vcs_solve_phaseStability.cpp new file mode 100644 index 000000000..5244414c8 --- /dev/null +++ b/Cantera/src/equil/vcs_solve_phaseStability.cpp @@ -0,0 +1,212 @@ +/** + * @file vcs_solve_TP.cpp Implementation file that contains the + * main algorithm for finding an equilibrium + */ +/* + * $Id: vcs_solve_TP.cpp 626 2010-10-28 01:33:54Z hkmoffa $ + */ +/* + * 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 +#include +#include +#include + +#include "vcs_solve.h" +#include "vcs_internal.h" +#include "vcs_VolPhase.h" +#include "vcs_species_thermo.h" +#include "vcs_prob.h" + +#include "clockWC.h" + +#ifdef WIN32 +#pragma warning(disable:4996) +#endif + +using namespace std; + +namespace VCSnonideal { + + + int VCS_SOLVE::vcs_PS(VCS_PROB *vprob, int iphase, int printLvl, double &feStable) { + + /* + * ifunc determines the problem type + */ + int ifunc = 0; + int iStab = 0; + + + + if (ifunc < 0 || ifunc > 2) { + plogf("vcs: Unrecognized value of ifunc, %d: bailing!\n", + ifunc); + return VCS_PUB_BAD; + } + + /* + * This function is called to copy the public data + * and the current problem specification + * into the current object's data structure. + */ + int retn = vcs_prob_specifyFully(vprob); + if (retn != 0) { + plogf("vcs_pub_to_priv returned a bad status, %d: bailing!\n", + retn); + return retn; + } + /* + * Prep the problem data + * - adjust the identity of any phases + * - determine the number of components in the problem + */ + retn = vcs_prep_oneTime(printLvl); + if (retn != 0) { + plogf("vcs_prep_oneTime returned a bad status, %d: bailing!\n", + retn); + return retn; + } + + + /* + * This function is called to copy the current problem + * into the current object's data structure. + */ + retn = vcs_prob_specify(vprob); + if (retn != 0) { + plogf("vcs_prob_specify returned a bad status, %d: bailing!\n", + retn); + return retn; + } + + + /* + * Prep the problem data for this particular instantiation of + * the problem + */ + retn = vcs_prep(); + if (retn != VCS_SUCCESS) { + plogf("vcs_prep returned a bad status, %d: bailing!\n", retn); + return retn; + } + /* + * Check to see if the current problem is well posed. + */ + if (!vcs_wellPosed(vprob)) { + plogf("vcs has determined the problem is not well posed: Bailing\n"); + return VCS_PUB_BAD; + } + + + int iconv; + /* + * Store the temperature and pressure in the private global variables + */ + m_temperature = vprob->T; + m_pressurePA = vprob->PresPA; + /* + * Evaluate the standard state free energies + * at the current temperatures and pressures. + */ + iconv = vcs_evalSS_TP(printLvl, printLvl, m_temperature, m_pressurePA); + + /* + * Prepare the problem data: + * ->nondimensionalize the free energies using + * the divisor, R * T + */ + vcs_nondim_TP(); + /* + * Prep the fe field + */ + vcs_fePrep_TP(); + + /* + * Solve the problem at a fixed Temperature and Pressure + * (all information concerning Temperature and Pressure has already + * been derived. The free energies are now in dimensionless form.) + */ + double funcVal; + iStab = vcs_solve_phaseStability(iphase, ifunc, funcVal, printLvl); + + + /* + * Redimensionalize the free energies using + * the reverse of vcs_nondim to add back units. + */ + vcs_redim_TP(); + /* + * Return the convergence success flag. + */ + return iStab; + + + } + //==================================================================================================================== + // Routine that independently determines whether a phase should be popped + // under the current conditions. + /* + * This is the main routine that solves for equilibrium at constant T and P + * using a variant of the VCS method. Nonideal phases can be accommodated + * as well. + * + * Any number of single-species phases and multi-species phases + * can be handled by the present version. + * + * Input + * ------------ + * @param print_lvl 1 -> Print results to standard output + * 0 -> don't report on anything + * + * @param printDetails 1 -> Print intermediate results. + * + * @param maxit Maximum number of iterations for the algorithm + * + * @return 0 = Equilibrium Achieved + * 1 = Range space error encountered. The element abundance criteria are + * only partially satisfied. Specifically, the first NC= (number of + * components) conditions are satisfied. However, the full NE + * (number of elements) conditions are not satisfied. The equilibrirum + * condition is returned. + * -1 = Maximum number of iterations is exceeded. Convergence was not + * found. + */ + int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, + double &funcVal, + int printLv) { + int retn = 0; + double test = -1.0E-10; + int usedZeroedSpecies; + std::vector phasePopPhaseIDs(0); + int iphasePop; + int iStab = 0; + + std::vector sm(m_numElemConstraints*m_numElemConstraints, 0.0); + std::vector ss(m_numElemConstraints, 0.0); + std::vector sa(m_numElemConstraints, 0.0); + + std::vector aw(m_numSpeciesTot, 0.0); + std::vector wx(m_numElemConstraints, 0.0); + + + retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), + VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), + test, &usedZeroedSpecies); + phasePopPhaseIDs.clear(); + iphasePop = vcs_popPhaseID(phasePopPhaseIDs); + funcVal = vcs_phaseStabilityTest(iph); + if (funcVal > 0.0) { + iStab = 1; + } else { + iStab = 0; + } + + return iStab; + } + +}