research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Incremental commit of new ability to determine phase stability

Browse source 
without actually solving the equilibrium equations.
This commit is contained in:
Harry Moffat 2010-12-15 20:20:13 +00:00
parent e10fa97af0
commit 32816b1fa7
7 changed files with 519 additions and 25 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/equil/Makefile.in
View file

@ -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 \

202
Cantera/src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -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;
}
//====================================================================================================================
}

47
Cantera/src/equil/vcs_MultiPhaseEquil.h
View file

@ -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;

72
Cantera/src/equil/vcs_equilibrate.cpp
View file

@ -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;
}
//====================================================================================================================
}

2
Cantera/src/equil/vcs_solve.cpp
View file

@ -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;

7
Cantera/src/equil/vcs_solve.h
View file

@ -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

212
Cantera/src/equil/vcs_solve_phaseStability.cpp Normal file
View file

@ -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 <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>
#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<int> phasePopPhaseIDs(0);
int iphasePop;
int iStab = 0;
std::vector<double> sm(m_numElemConstraints*m_numElemConstraints, 0.0);
std::vector<double> ss(m_numElemConstraints, 0.0);
std::vector<double> sa(m_numElemConstraints, 0.0);
std::vector<double> aw(m_numSpeciesTot, 0.0);
std::vector<double> 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;
}
}