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

Changed the interface to allow values of 0, 1, and -1 for

Browse source 
estimateEquil.
  0 don't estimate
  1 estimate only if element abundances aren't satisfied
 -1 always estimate.
This commit is contained in:
Harry Moffat 2008-05-06 16:53:37 +00:00
parent 7b4e4ec1c4
commit 9da8402543
11 changed files with 504 additions and 352 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -127,7 +127,7 @@ namespace Cantera {
if (solver >= 2) {
#ifdef WITH_VCSNONIDEAL
int printLvlSub = 0;
bool estimateEquil = false;
int estimateEquil = 0;
m = new MultiPhase;
try {
m->addPhase(&s, 1.0);

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

@ -91,7 +91,7 @@ namespace VCSnonideal {
}
int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
bool estimateEquil,
int estimateEquil,
int printLvl, doublereal err,
int maxsteps, int loglevel) {
@ -111,7 +111,7 @@ namespace VCSnonideal {
m_mix->setTemperature(xtarget);
}
double Pnew;
bool strt = estimateEquil;
int strt = estimateEquil;
double P1 = 0.0;
double V1 = 0.0;
double V2 = 0.0;
@ -212,7 +212,7 @@ namespace VCSnonideal {
int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
int XY, double Tlow, double Thigh,
bool estimateEquil,
int estimateEquil,
int printLvl, doublereal err,
int maxsteps, int loglevel) {
int maxiter = 100;
@ -221,7 +221,7 @@ namespace VCSnonideal {
throw CanteraError("vcs_MultiPhaseEquil::equilibrate_HP",
"Wrong XP" + XY);
}
bool strt = estimateEquil;
int strt = estimateEquil;
// Lower bound on T. This will change as we progress in the calculation
if (Tlow <= 0.0) {
@ -252,7 +252,7 @@ namespace VCSnonideal {
try {
Tnow = m_mix->temperature();
iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
strt = false;
strt = 0;
if (XY == UP) {
Hnow = m_mix->IntEnergy();
} else {
@ -339,7 +339,7 @@ namespace VCSnonideal {
if (!estimateEquil) {
addLogEntry("no convergence",
"try estimating composition at the start");
strt = true;
strt = -1;
}
else {
Tnew = 0.5*(Tnow + Thigh);
@ -362,12 +362,12 @@ namespace VCSnonideal {
int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
double Tlow, double Thigh,
bool estimateEquil,
int estimateEquil,
int printLvl, doublereal err,
int maxsteps, int loglevel) {
int maxiter = 100;
int iSuccess;
bool strt = estimateEquil;
int strt = estimateEquil;
// Lower bound on T. This will change as we progress in the calculation
if (Tlow <= 0.0) {
@ -403,7 +403,7 @@ namespace VCSnonideal {
try {
Tnow = m_mix->temperature();
iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
strt = false;
strt = 0;
Snow = m_mix->entropy();
double pmoles[10];
pmoles[0] = m_mix->phaseMoles(0);
@ -495,7 +495,7 @@ namespace VCSnonideal {
if (!estimateEquil) {
addLogEntry("no convergence",
"try estimating composition at the start");
strt = true;
strt = -1;
}
else {
Tnew = 0.5*(Tnow + Thigh);
@ -520,7 +520,7 @@ namespace VCSnonideal {
/*
* Equilibrate the solution using the current element abundances
*/
int vcs_MultiPhaseEquil::equilibrate(int XY, bool estimateEquil,
int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil,
int printLvl, doublereal err,
int maxsteps, int loglevel) {
int iSuccess;
@ -570,7 +570,7 @@ namespace VCSnonideal {
/*
* Equilibrate the solution using the current element abundances
*/
int vcs_MultiPhaseEquil::equilibrate_TP(bool estimateEquil,
int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
int printLvl, doublereal err,
int maxsteps, int loglevel) {
// Debugging level
@ -598,7 +598,7 @@ namespace VCSnonideal {
// Set the estimation technique
if (estimateEquil) {
m_vprob->iest = -1;
m_vprob->iest = estimateEquil;
} else {
m_vprob->iest = 0;
}

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

@ -34,11 +34,17 @@ namespace Cantera {
* @param XY An integer specifying the two properties to be held
* constant.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -73,7 +79,7 @@ namespace Cantera {
* @ingroup equilfunctions
*/
int vcs_equilibrate(thermo_t& s, const char* XY,
bool estimateEquil = false, int printLvl = 0,
int estimateEquil = 0, int printLvl = 0,
int solver = -1, doublereal rtol = 1.0e-9,
int maxsteps = 5000,
int maxiter = 100, int loglevel = -99);
@ -93,11 +99,17 @@ namespace Cantera {
* @param XY A character string representing the unknowns
* to be held constant
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -128,7 +140,7 @@ namespace Cantera {
* @ingroup equilfunctions
*/
int vcs_equilibrate(MultiPhase& s, const char* XY,
bool estimateEquil = false, int printLvl = 0,
int estimateEquil = 0, int printLvl = 0,
int solver = 2,
doublereal rtol = 1.0e-9, int maxsteps = 5000,
int maxiter = 100, int loglevel = -99);
@ -147,11 +159,17 @@ namespace Cantera {
* @param ixy An integer specifying the two properties to be held
* constant.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -182,7 +200,7 @@ namespace Cantera {
* @ingroup equilfunctions
*/
int vcs_equilibrate_1(MultiPhase& s, int ixy,
bool estimateEquil = false, int printLvl = 0,
int estimateEquil = 0, int printLvl = 0,
int solver = 2,
doublereal rtol = 1.0e-9, int maxsteps = 5000,
int maxiter = 100, int loglevel = -99);
@ -315,11 +333,17 @@ namespace VCSnonideal {
* @param XY Integer representing what two thermo quantities
* are held constant during the equilibration
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -329,7 +353,7 @@ namespace VCSnonideal {
* @param maxsteps max steps allowed.
* @param loglevel for
*/
int equilibrate(int XY, bool estimateEquil = false,
int equilibrate(int XY, int estimateEquil = 0,
int printLvl= 0, doublereal err = 1.0e-6,
int maxsteps = 5000, int loglevel=-99);
@ -339,11 +363,17 @@ namespace VCSnonideal {
* Use the vcs algorithm to equilibrate the current multiphase
* mixture.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -353,7 +383,7 @@ namespace VCSnonideal {
* @param maxsteps max steps allowed.
* @param loglevel for
*/
int equilibrate_TP(bool estimateEquil = false,
int equilibrate_TP(int estimateEquil = 0,
int printLvl= 0, doublereal err = 1.0e-6,
int maxsteps = 5000, int loglevel=-99);
@ -381,11 +411,17 @@ namespace VCSnonideal {
* error condition if the temperature goes
* higher than Thigh.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -401,7 +437,7 @@ namespace VCSnonideal {
* output file.
*/
int equilibrate_HP(doublereal Htarget, int XY, double Tlow, double Thigh,
bool estimateEquil = false,
int estimateEquil = 0,
int printLvl = 0, doublereal err = 1.0E-6,
int maxsteps = 5000, int loglevel=-99);
@ -426,11 +462,17 @@ namespace VCSnonideal {
* error condition if the temperature goes
* higher than Thigh.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -446,7 +488,7 @@ namespace VCSnonideal {
* output file.
*/
int equilibrate_SP(doublereal Starget, double Tlow, double Thigh,
bool estimateEquil = false,
int estimateEquil = 0,
int printLvl = 0, doublereal err = 1.0E-6,
int maxsteps = 5000, int loglevel=-99);
@ -468,12 +510,17 @@ namespace VCSnonideal {
* Note, except for T, this must be an extensive
* quantity. units = Joules/K or Joules
*
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -489,7 +536,7 @@ namespace VCSnonideal {
* output file.
*/
int equilibrate_TV(int XY, doublereal xtarget,
bool estimateEquil = false,
int estimateEquil = 0,
int printLvl = 0, doublereal err = 1.0E-6,
int maxsteps = 5000, int loglevel = -99);

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

@ -76,7 +76,7 @@ namespace VCSnonideal {
* Decide whether we need an initial estimate of the solution
* If so, go get one. If not, then
*/
if (iest == -1) {
if (m_doEstimateEquil) {
retn = vcs_inest_TP();
if (retn != VCS_SUCCESS) {
plogf("vcs_inest_TP returned a failure flag\n");

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

@ -50,11 +50,17 @@ namespace Cantera {
* @param XY An integer specifying the two properties to be held
* constant.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -84,7 +90,7 @@ namespace Cantera {
* in a web browser. @see HTML_logs
*/
int vcs_equilibrate(thermo_t& s, const char* XY,
bool estimateEquil, int printLvl,
int estimateEquil, int printLvl,
int solver,
doublereal rtol, int maxsteps, int maxiter,
int loglevel) {
@ -203,11 +209,17 @@ namespace Cantera {
* @param XY A character string specifying the two properties to
* be held constant
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -231,7 +243,7 @@ namespace Cantera {
* @ingroup equilfunctions
*/
int vcs_equilibrate(MultiPhase& s, const char* XY,
bool estimateEquil, int printLvl, int solver,
int estimateEquil, int printLvl, int solver,
doublereal tol, int maxsteps, int maxiter,
int loglevel) {
int ixy = _equilflag(XY);
@ -255,11 +267,17 @@ namespace Cantera {
* @param XY An integer specifying the two properties to be held
* constant.
*
* @param estimateEquil Boolean indicating whether the solver
* @param estimateEquil integer indicating whether the solver
* should estimate its own initial condition.
* If false, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* If 0, the initial mole fraction vector
* in the %ThermoPhase object is used as the
* initial condition.
* If 1, the initial mole fraction vector
* is used if the element abundances are
* satisfied.
* if -1, the initial mole fraction vector
* is thrown out, and an estimate is
* formulated.
*
* @param printLvl Determines the amount of printing that
* gets sent to stdout from the vcs package
@ -283,7 +301,7 @@ namespace Cantera {
* @ingroup equilfunctions
*/
int vcs_equilibrate_1(MultiPhase& s, int ixy,
bool estimateEquil, int printLvl, int solver,
int estimateEquil, int printLvl, int solver,
doublereal tol, int maxsteps, int maxiter, int loglevel) {
static int counter = 0;
int retn = 1;

88
Cantera/src/equil/vcs_inest.cpp
View file

@ -382,32 +382,48 @@ namespace VCSnonideal {
}
}
#endif
} /* inest() *****************************************************************/
}
/***************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
int VCS_SOLVE::vcs_inest_TP(void)
/**************************************************************************
*
* vcs_inest_TP:
*
* Create an initial estimate of the solution to the thermodynamic
* equilibrium problem.
*
* Return value:
*
* 0: successful initial guess
* -1: Unsuccessful initial guess, the elemental abundances aren't
* satisfied.
***************************************************************************/
{
// Create an initial estimate of the solution to the thermodynamic
// equilibrium problem.
/*
* @return Return value indicates success:
* - 0: successful initial guess
* - -1: Unsuccessful initial guess; the elemental abundances aren't
* satisfied.
*/
int VCS_SOLVE::vcs_inest_TP() {
int retn = 0;
double test;
Cantera::clockWC tickTock;
test = -1.0E20;
if (m_doEstimateEquil > 0) {
/*
* Calculate the elemental abundances
*/
vcs_elab();
if (vcs_elabcheck(0)) {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%s Initial guess passed element abundances on input\n", pprefix);
plogf("%s m_doEstimateEquil = 1 so will use the input mole numbers as estimates", pprefix);
plogendl();
}
#endif
return retn;
#ifdef DEBUG_MODE
} else {
if (vcs_debug_print_lvl >= 2) {
plogf("%s Initial guess failed element abundances on input\n", pprefix);
plogf("%s m_doEstimateEquil = 1 so will discard input mole numbers and find our own estimate", pprefix);
plogendl();
}
#endif
}
}
/*
* Malloc temporary space for usage in this routine and in
* subroutines
@ -416,8 +432,6 @@ namespace VCSnonideal {
* sa[ne]
* aw[m]
*/
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);
@ -428,22 +442,23 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%sGo find an initial estimate for the equilibrium problem",
pprefix);
pprefix);
plogendl();
}
#endif
inest(VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm),
VCS_DATA_PTR(ss), test);
/*
* Calculate the elemental abundances
* Calculate the elemental abundances
*/
vcs_elab();
/*
* If we still fail to achieve the correct elemental abundances,
* try to fix the problem again by calling the main elemental abundances
* fixer routine, used in the main program. What this does, is that it
* fixer routine, used in the main program. This
* attempts to tweak the mole numbers of the component species to
* satisfy the element abundance constraints.
* satisfy the element abundance constraints.
*
* Note: We won't do this unless we have to since it involves inverting
* a matrix.
@ -461,9 +476,9 @@ namespace VCSnonideal {
rangeCheck = vcs_elabcheck(1);
if (!vcs_elabcheck(0)) {
plogf("%sInitial guess still fails element abundance equations\n",
pprefix);
pprefix);
plogf("%s - Inability to ever satisfy element abundance "
"constraints is probable", pprefix);
"constraints is probable", pprefix);
plogendl();
retn = -1;
} else {
@ -475,8 +490,8 @@ namespace VCSnonideal {
} else {
plogf("%sElement Abundances RANGE ERROR\n", pprefix);
plogf("%s - Initial guess satisfies NC=%d element abundances, "
"BUT not NE=%d element abundances", pprefix,
m_numComponents, m_numElemConstraints);
"BUT not NE=%d element abundances", pprefix,
m_numComponents, m_numElemConstraints);
plogendl();
}
}
@ -492,8 +507,8 @@ namespace VCSnonideal {
} else {
plogf("%sElement Abundances RANGE ERROR\n", pprefix);
plogf("%s - Initial guess satisfies NC=%d element abundances, "
"BUT not NE=%d element abundances", pprefix,
m_numComponents, m_numElemConstraints);
"BUT not NE=%d element abundances", pprefix,
m_numComponents, m_numElemConstraints);
plogendl();
}
}
@ -503,14 +518,14 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(m_tPhaseMoles_old)));
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(m_tPhaseMoles_old)));
plogendl();
}
#endif
/*
* Free malloced memory
* Record time
*/
double tsecond = tickTock.secondsWC();
m_VCount->T_Time_inest += tsecond;
@ -519,4 +534,3 @@ namespace VCSnonideal {
}
}

494
Cantera/src/equil/vcs_prep.cpp
View file

@ -24,73 +24,66 @@
namespace VCSnonideal {
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::vcs_SSPhase(void)
/**************************************************************************
*
* vcs_SSPhase:
*
* Calculate the status of single species phases.
*
*************************************************************************/
{
int kspec, iph;
vcs_VolPhase *Vphase;
void VCS_SOLVE::vcs_SSPhase(void)
/**************************************************************************
*
* vcs_SSPhase:
*
* Calculate the status of single species phases.
*
*************************************************************************/
{
int kspec, iph;
vcs_VolPhase *Vphase;
std::vector<int> numPhSpecies(NPhase, 0);
std::vector<int> numPhSpecies(NPhase, 0);
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
numPhSpecies[PhaseID[kspec]]++;
}
/*
* Handle the special case of a single species in a phase that
* has been earmarked as a multispecies phase.
* Treat that species as a single-species phase
*/
for (iph = 0; iph < NPhase; iph++) {
Vphase = VPhaseList[iph];
Vphase->SingleSpecies = false;
if (TPhInertMoles[iph] > 0.0) {
Vphase->Existence = 2;
}
if (numPhSpecies[iph] <= 1) {
if (TPhInertMoles[iph] == 0.0) {
Vphase->SingleSpecies = true;
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
numPhSpecies[PhaseID[kspec]]++;
}
/*
* Handle the special case of a single species in a phase that
* has been earmarked as a multispecies phase.
* Treat that species as a single-species phase
*/
for (iph = 0; iph < NPhase; iph++) {
Vphase = VPhaseList[iph];
Vphase->SingleSpecies = false;
if (TPhInertMoles[iph] > 0.0) {
Vphase->Existence = 2;
}
}
Vphase->NVolSpecies = numPhSpecies[iph];
}
if (numPhSpecies[iph] <= 1) {
if (TPhInertMoles[iph] == 0.0) {
Vphase->SingleSpecies = true;
}
}
Vphase->NVolSpecies = numPhSpecies[iph];
}
/*
* Fill in some useful arrays here that have to do with the
* static information concerning the phase ID of species.
* SSPhase = Boolean indicating whether a species is in a
* single species phase or not.
*/
for (kspec = 0; kspec < m_numSpeciesTot; kspec++) {
iph = PhaseID[kspec];
Vphase = VPhaseList[iph];
if (Vphase->SingleSpecies) SSPhase[kspec] = TRUE;
else SSPhase[kspec] = FALSE;
/*
* Fill in some useful arrays here that have to do with the
* static information concerning the phase ID of species.
* SSPhase = Boolean indicating whether a species is in a
* single species phase or not.
*/
for (kspec = 0; kspec < m_numSpeciesTot; kspec++) {
iph = PhaseID[kspec];
Vphase = VPhaseList[iph];
if (Vphase->SingleSpecies) SSPhase[kspec] = TRUE;
else SSPhase[kspec] = FALSE;
}
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
/**************************************************************************
*
* vcs_prep_oneTime:
*
* This routine is mostly concerned with changing the private data
* to be consistent with what's needed for solution. It is called one
* time for each new problem structure definition.
*
/*****************************************************************************/
// This routine is mostly concerned with changing the private data
// to be consistent with what's needed for solution. It is called one
// time for each new problem structure definition.
/*
* This routine is always followed by vcs_prep(). Therefore, tasks
* that need to be done for every call to vcsc() should be placed in
* vcs_prep() and not in this routine.
@ -114,213 +107,220 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
* in the private data structure. All references to the species
* properties must employ the ind[] index vector.
*
* return code
* VCS_SUCCESS = everything went OK
* @param printLvl Print level of the routine
*
**************************************************************************/
{
int kspec, i, conv, retn = VCS_SUCCESS;
double pres, test;
double *aw, *sa, *sm, *ss;
bool modifiedSoln = false;
* @return the return code
* VCS_SUCCESS = everything went OK
*
*/
int VCS_SOLVE::vcs_prep_oneTime(int printLvl) {
int kspec, i, conv, retn = VCS_SUCCESS;
double pres, test;
double *aw, *sa, *sm, *ss;
bool modifiedSoln = false;
#ifdef DEBUG_MODE
vcs_debug_print_lvl = printLvl;
vcs_debug_print_lvl = printLvl;
#endif
/*
* Calculate the Single Species status of phases
* Also calculate the number of species per phase
*/
vcs_SSPhase();
/*
* Calculate the Single Species status of phases
* Also calculate the number of species per phase
*/
vcs_SSPhase();
/*
* Set an initial estimate for the number of noncomponent species
* equal to nspecies - nelements. This may be changed below
*/
m_numRxnTot = m_numSpeciesTot - m_numElemConstraints;
m_numRxnRdc = m_numRxnTot;
m_numSpeciesRdc = m_numSpeciesTot;
for (i = 0; i < m_numRxnRdc; ++i) {
ir[i] = m_numElemConstraints + i;
}
/*
* Set an initial estimate for the number of noncomponent species
* equal to nspecies - nelements. This may be changed below
*/
m_numRxnTot = m_numSpeciesTot - m_numElemConstraints;
m_numRxnRdc = m_numRxnTot;
m_numSpeciesRdc = m_numSpeciesTot;
for (i = 0; i < m_numRxnRdc; ++i) {
ir[i] = m_numElemConstraints + i;
}
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
int pID = PhaseID[kspec];
int spPhIndex = indPhSp[kspec];
vcs_VolPhase *vPhase = VPhaseList[pID];
vcs_SpeciesProperties *spProp = vPhase->ListSpeciesPtr[spPhIndex];
double sz = 0.0;
int eSize = spProp->FormulaMatrixCol.size();
for (int e = 0; e < eSize; e++) {
sz += fabs(spProp->FormulaMatrixCol[e]);
}
if (sz > 0.0) {
m_spSize[kspec] = sz;
} else {
m_spSize[kspec] = 1.0;
}
}
/* ***************************************************** */
/* **** DETERMINE THE NUMBER OF COMPONENTS ************* */
/* ***************************************************** */
/*
* Obtain a valid estimate of the mole fraction. This will
* be used as an initial ordering vector for prioritizing
* which species are defined as components.
*
* If a mole number estimate was supplied from the
* input file, use that mole number estimate.
*
* If a solution estimate wasn't supplied from the input file,
* supply an initial estimate for the mole fractions
* based on the relative reverse ordering of the
* chemical potentials.
*
* For voltage unknowns, set these to zero for the moment.
*/
test = -1.0e-10;
if (iest < 0) {
double sum = 0.0;
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
sum += fabs(m_molNumSpecies_old[kspec]);
int pID = PhaseID[kspec];
int spPhIndex = indPhSp[kspec];
vcs_VolPhase *vPhase = VPhaseList[pID];
vcs_SpeciesProperties *spProp = vPhase->ListSpeciesPtr[spPhIndex];
double sz = 0.0;
int eSize = spProp->FormulaMatrixCol.size();
for (int e = 0; e < eSize; e++) {
sz += fabs(spProp->FormulaMatrixCol[e]);
}
if (sz > 0.0) {
m_spSize[kspec] = sz;
} else {
m_spSize[kspec] = 1.0;
}
}
if (fabs(sum) < 1.0E-6) {
modifiedSoln = true;
if (m_pressurePA <= 0.0) pres = 1.01325E5;
else pres = m_pressurePA;
retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
/* ***************************************************** */
/* **** DETERMINE THE NUMBER OF COMPONENTS ************* */
/* ***************************************************** */
/*
* Obtain a valid estimate of the mole fraction. This will
* be used as an initial ordering vector for prioritizing
* which species are defined as components.
*
* If a mole number estimate was supplied from the
* input file, use that mole number estimate.
*
* If a solution estimate wasn't supplied from the input file,
* supply an initial estimate for the mole fractions
* based on the relative reverse ordering of the
* chemical potentials.
*
* For voltage unknowns, set these to zero for the moment.
*/
test = -1.0e-10;
if (m_doEstimateEquil < 0) {
double sum = 0.0;
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
} else {
m_molNumSpecies_old[kspec] = 0.0;
sum += fabs(m_molNumSpecies_old[kspec]);
}
}
if (fabs(sum) < 1.0E-6) {
modifiedSoln = true;
if (m_pressurePA <= 0.0) pres = 1.01325E5;
else pres = m_pressurePA;
retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
} else {
m_molNumSpecies_old[kspec] = 0.0;
}
}
}
test = -1.0e20;
}
test = -1.0e20;
}
/*
* NC = number of components is in the vcs.h common block
* This call to BASOPT doesn't calculate the stoichiometric
* reaction matrix.
*/
std::vector<double> awSpace( m_numSpeciesTot + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
aw = VCS_DATA_PTR(awSpace);
if (aw == NULL) {
plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
return VCS_NOMEMORY;
}
sa = aw + m_numSpeciesTot;
sm = sa + m_numElemConstraints;
ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
retn = vcs_basopt(TRUE, aw, sa, sm, ss, test, &conv);
if (retn != VCS_SUCCESS) {
plogf("vcs_prep_oneTime:");
plogf(" Determination of number of components failed: %d\n",
retn);
plogf(" Global Bailout!\n");
return retn;
}
if (m_numElemConstraints != m_numComponents) {
m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents;
for (i = 0; i < m_numRxnRdc; ++i) {
ir[i] = m_numComponents + i;
/*
* NC = number of components is in the vcs.h common block
* This call to BASOPT doesn't calculate the stoichiometric
* reaction matrix.
*/
std::vector<double> awSpace( m_numSpeciesTot + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
aw = VCS_DATA_PTR(awSpace);
if (aw == NULL) {
plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
return VCS_NOMEMORY;
}
sa = aw + m_numSpeciesTot;
sm = sa + m_numElemConstraints;
ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
retn = vcs_basopt(TRUE, aw, sa, sm, ss, test, &conv);
if (retn != VCS_SUCCESS) {
plogf("vcs_prep_oneTime:");
plogf(" Determination of number of components failed: %d\n",
retn);
plogf(" Global Bailout!\n");
return retn;
}
if (m_numElemConstraints != m_numComponents) {
m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents;
for (i = 0; i < m_numRxnRdc; ++i) {
ir[i] = m_numComponents + i;
}
}
}
/*
* The elements might need to be rearranged.
*/
awSpace.resize(m_numElemConstraints + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
aw = VCS_DATA_PTR(awSpace);
sa = aw + m_numElemConstraints;
sm = sa + m_numElemConstraints;
ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
retn = vcs_elem_rearrange(aw, sa, sm, ss);
if (retn != VCS_SUCCESS) {
plogf("vcs_prep_oneTime:");
plogf(" Determination of element reordering failed: %d\n",
retn);
plogf(" Global Bailout!\n");
return retn;
}
// If we mucked up the solution unknowns because they were all
// zero to start with, set them back to zero here
if (modifiedSoln) {
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_molNumSpecies_old[kspec] = 0.0;
/*
* The elements might need to be rearranged.
*/
awSpace.resize(m_numElemConstraints + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
aw = VCS_DATA_PTR(awSpace);
sa = aw + m_numElemConstraints;
sm = sa + m_numElemConstraints;
ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
retn = vcs_elem_rearrange(aw, sa, sm, ss);
if (retn != VCS_SUCCESS) {
plogf("vcs_prep_oneTime:");
plogf(" Determination of element reordering failed: %d\n",
retn);
plogf(" Global Bailout!\n");
return retn;
}
}
return VCS_SUCCESS;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
// Prepare the object for re-solution
/*
* This routine is mostly concerned with changing the private data
* to be consistent with that needed for solution. It is called for
* every invocation of the vcs_solve() except for the cleanup invocation.
*
* Tasks:
* 1) Initialization of arrays to zero.
* 2) Calculate total number of moles in all phases
*
* return code
* VCS_SUCCESS = everything went OK
* VCS_PUB_BAD = There is an irreconcilable difference in the
* public data structure from when the problem was
* initially set up.
*/
int VCS_SOLVE::vcs_prep(void) {
/*
* Initialize various arrays in the data to zero
*/
vcs_dzero(VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesTot);
vcs_vdzero(m_feSpecies_old, m_numSpeciesTot);
vcs_vdzero(m_molNumSpecies_new, m_numSpeciesTot);
vcs_dzero(&(DnPhase[0][0]), m_numSpeciesTot*NPhase);
vcs_izero(&(PhaseParticipation[0][0]), m_numSpeciesTot*NPhase);
vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), NPhase);
vcs_dzero(VCS_DATA_PTR(m_tPhaseMoles_new), NPhase);
/*
* Calculate the total number of moles in all phases.
*/
vcs_tmoles();
return VCS_SUCCESS;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
bool VCS_SOLVE::vcs_wellPosed(VCS_PROB *vprob)
/**************************************************************************
*
* vcs_wellPosed:
*
* In this routine, we check for things that will cause the algorithm
* to fail.
*
**************************************************************************/
{
double sum = 0.0;
for (int e = 0; e < vprob->ne; e++) {
sum = sum + vprob->gai[e];
// If we mucked up the solution unknowns because they were all
// zero to start with, set them back to zero here
if (modifiedSoln) {
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_molNumSpecies_old[kspec] = 0.0;
}
}
return VCS_SUCCESS;
}
if (sum < 1.0E-20) {
plogf("vcs_wellPosed: Element abundance is close to zero\n");
return false;
/*****************************************************************************/
// Prepare the object for re-solution
/*
* This routine is mostly concerned with changing the private data
* to be consistent with that needed for solution. It is called for
* every invocation of the vcs_solve() except for the cleanup invocation.
*
* Tasks:
* 1) Initialization of arrays to zero.
* 2) Calculate total number of moles in all phases
*
* return code
* VCS_SUCCESS = everything went OK
* VCS_PUB_BAD = There is an irreconcilable difference in the
* public data structure from when the problem was
* initially set up.
*/
int VCS_SOLVE::vcs_prep() {
/*
* Initialize various arrays in the data to zero
*/
vcs_dzero(VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesTot);
vcs_vdzero(m_feSpecies_old, m_numSpeciesTot);
vcs_vdzero(m_molNumSpecies_new, m_numSpeciesTot);
vcs_dzero(&(DnPhase[0][0]), m_numSpeciesTot*NPhase);
vcs_izero(&(PhaseParticipation[0][0]), m_numSpeciesTot*NPhase);
vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), NPhase);
vcs_dzero(VCS_DATA_PTR(m_tPhaseMoles_new), NPhase);
/*
* Calculate the total number of moles in all phases.
*/
vcs_tmoles();
return VCS_SUCCESS;
}
return true;
}
/*****************************************************************************/
/*****************************************************************************/
// In this routine, we check for things that will cause the algorithm
// to fail.
/*
* We check to see if the problem is well posed. If it is not, we return
* false and print out error conditions.
*
* Current there is one condition. If all the element abundances are
* zero, the algorithm will fail
*
* @param vprob VCS_PROB pointer to the definition of the equilibrium
* problem
*
* @return If true, the problem is well-posed. If false, the problem
* is not well posed.
*/
bool VCS_SOLVE::vcs_wellPosed(VCS_PROB *vprob) {
double sum = 0.0;
for (int e = 0; e < vprob->ne; e++) {
sum = sum + vprob->gai[e];
}
if (sum < 1.0E-20) {
plogf("vcs_wellPosed: Element abundance is close to zero\n");
return false;
}
return true;
}
/*****************************************************************************/
}

1
Cantera/src/equil/vcs_prob.h
View file

@ -177,6 +177,7 @@ namespace VCSnonideal {
//! Specification of the initial estimate method
/*!
* iest = Initial estimate: 0 user estimate
* 1 user estimate if satisifies elements
* -1 machine estimate
*/
int iest;

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

@ -41,7 +41,7 @@ namespace VCSnonideal {
m_numSpeciesRdc(0),
m_numRxnMinorZeroed(0),
NPhase(0),
iest(0),
m_doEstimateEquil(0),
TMoles(0.0),
m_temperature(0.0),
m_pressurePA(0.0),
@ -536,7 +536,7 @@ namespace VCSnonideal {
/*
* iest => Do we have an initial estimate of the species mole numbers ?
*/
iest = pub->iest;
m_doEstimateEquil = pub->iest;
/*
* w[] -> Copy the equilibrium mole number estimate if it exists.
@ -544,7 +544,7 @@ namespace VCSnonideal {
if (pub->w.size() != 0) {
vcs_vdcopy(m_molNumSpecies_old, pub->w, nspecies);
} else {
iest = -1;
m_doEstimateEquil = -1;
vcs_dzero(VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
}
@ -554,7 +554,7 @@ namespace VCSnonideal {
if (pub->gai.size() != 0) {
for (i = 0; i < nelements; i++) m_elemAbundancesGoal[i] = pub->gai[i];
} else {
if (iest == 0) {
if (m_doEstimateEquil == 0) {
for (j = 0; j < nelements; j++) {
m_elemAbundancesGoal[j] = 0.0;
for (kspec = 0; kspec < nspecies; kspec++) {
@ -775,7 +775,7 @@ namespace VCSnonideal {
m_temperature = pub->T;
m_pressurePA = pub->PresPA;
m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
iest = pub->iest;
m_doEstimateEquil = pub->iest;
Vol = pub->Vol;

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

@ -176,9 +176,43 @@ public:
void vcs_fePrep_TP(void);
double vcs_VolTotal(double, double, double [], double []);
//! This routine is mostly concerned with changing the private data
//! to be consistent with what's needed for solution. It is called one
//! time for each new problem structure definition.
/*!
* This routine is always followed by vcs_prep(). Therefore, tasks
* that need to be done for every call to vcsc() should be placed in
* vcs_prep() and not in this routine.
*
* The problem structure refers to:
*
* the number and identity of the species.
* the formula matrix and thus the number of components.
* the number and identity of the phases.
* the equation of state
* the method and parameters for determining the standard state
* The method and parameters for determining the activity coefficients.
*
* Tasks:
* 0) Fill in the SSPhase[] array.
* 1) Check to see if any multispecies phases actually have only one
* species in that phase. If true, reassign that phase and species
* to be a single-species phase.
* 2) Determine the number of components in the problem if not already
* done so. During this process the order of the species is changed
* in the private data structure. All references to the species
* properties must employ the ind[] index vector.
*
* @param printLvl Print level of the routine
*
* @return the return code
* VCS_SUCCESS = everything went OK
*
*/
int vcs_prep_oneTime(int printLvl);
//! Prepare the object for resolution
//! Prepare the object for solution
/*!
* This routine is mostly concerned with changing the private data
* to be consistent with that needed for solution. It is called for
@ -193,8 +227,23 @@ public:
* public data structure from when the problem was
* initially set up.
*/
int vcs_prep(void);
int vcs_prep();
//! In this routine, we check for things that will cause the algorithm
//! to fail.
/*!
* We check to see if the problem is well posed. If it is not, we return
* false and print out error conditions.
*
* Current there is one condition. If all the element abundances are
* zero, the algorithm will fail.
*
* @param vprob VCS_PROB pointer to the definition of the equilibrium
* problem
*
* @return If true, the problem is well-posed. If false, the problem
* is not well posed.
*/
bool vcs_wellPosed(VCS_PROB *vprob);
//! Rearrange the constraint equations represented by the Formula
@ -258,7 +307,16 @@ public:
double vcs_line_search(int irxn, double dx_orig);
#endif
int vcs_report(int);
//! Print out a report on the state of the equilibrium problem to
//! standard output.
/*!
* @param iconv Indicator of convergence, to be printed out in the report:
* - 0 converged
* - 1 range space error
* - -1 not converged
*/
int vcs_report(int iconv);
int vcs_rearrange(void);
@ -281,7 +339,16 @@ public:
void vcs_elabPhase(int iphase, double * const elemAbundPhase);
int vcs_elcorr(double aa[], double x[]);
int vcs_inest_TP(void);
//! Create an initial estimate of the solution to the thermodynamic
//! equilibrium problem.
/*!
* @return Return value indicates success:
* - 0: successful initial guess
* - -1: Unsuccessful initial guess; the elemental abundances aren't
* satisfied.
*/
int vcs_inest_TP();
#ifdef ALTLINPROG
//! Extimate the initial mole numbers by constrained linear programming
@ -535,12 +602,15 @@ public:
*/
std::vector<double> m_feSpecies_new;
//! Setting for the initial estimate
//! Setting for whether to do an initial estimate
/*!
* Initial estimate: 0 user estimate
* -1 machine estimate
* Initial estimate: 0 Do not estimate the solution at all. Use the supplied
* mole numbers as is.
* 1 Only do an estimate if the element abundances aren't satisfied.
* -1 Force an estimate of the soln. Throw out the input
* mole numbers.
*/
int iest;
int m_doEstimateEquil;
//! Total moles of the species
/*!
@ -573,7 +643,7 @@ public:
DoubleStarStar DnPhase;
//! This is 1 if the phase, iphase, participates in the formation reaction
//! irxn, and zero otherwise. PhaseParticipation[irxn][iphase]
//! irxn, and zero otherwise. PhaseParticipation[irxn][iphase]
IntStarStar PhaseParticipation;
//! electric potential of the iph phase

8
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -257,10 +257,12 @@ namespace VCSnonideal {
plogf("%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i],
m_elType[i]);
}
if (iest < 0) {
plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM\n");
if (m_doEstimateEquil < 0) {
plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - forced\n");
} else if (m_doEstimateEquil > 0) {
plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - where necessary\n");
}
if (iest >= 0) {
if (m_doEstimateEquil == 0) {
plogf("\n USER ESTIMATE OF EQUILIBRIUM\n");
}
if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {