From 9da8402543227b239da39ddd93d9f05be05605bf Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 6 May 2008 16:53:37 +0000 Subject: [PATCH] Changed the interface to allow values of 0, 1, and -1 for estimateEquil. 0 don't estimate 1 estimate only if element abundances aren't satisfied -1 always estimate. --- Cantera/src/equil/equilibrate.cpp | 2 +- Cantera/src/equil/vcs_MultiPhaseEquil.cpp | 26 +- Cantera/src/equil/vcs_MultiPhaseEquil.h | 97 +++-- Cantera/src/equil/vcs_TP.cpp | 2 +- Cantera/src/equil/vcs_equilibrate.cpp | 40 +- Cantera/src/equil/vcs_inest.cpp | 88 ++-- Cantera/src/equil/vcs_prep.cpp | 494 +++++++++++----------- Cantera/src/equil/vcs_prob.h | 1 + Cantera/src/equil/vcs_solve.cpp | 10 +- Cantera/src/equil/vcs_solve.h | 88 +++- Cantera/src/equil/vcs_solve_TP.cpp | 8 +- 11 files changed, 504 insertions(+), 352 deletions(-) diff --git a/Cantera/src/equil/equilibrate.cpp b/Cantera/src/equil/equilibrate.cpp index 65e2307f8..f91dedf88 100644 --- a/Cantera/src/equil/equilibrate.cpp +++ b/Cantera/src/equil/equilibrate.cpp @@ -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); diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp index 7a4cedad4..d861094da 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp @@ -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; } diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.h b/Cantera/src/equil/vcs_MultiPhaseEquil.h index 2f3dd4fcb..3b35973a4 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.h +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.h @@ -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); diff --git a/Cantera/src/equil/vcs_TP.cpp b/Cantera/src/equil/vcs_TP.cpp index 19c3df783..26ee17359 100644 --- a/Cantera/src/equil/vcs_TP.cpp +++ b/Cantera/src/equil/vcs_TP.cpp @@ -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"); diff --git a/Cantera/src/equil/vcs_equilibrate.cpp b/Cantera/src/equil/vcs_equilibrate.cpp index d33518b78..f0c7072fe 100644 --- a/Cantera/src/equil/vcs_equilibrate.cpp +++ b/Cantera/src/equil/vcs_equilibrate.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index c33e8bd07..f28cb0220 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -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 sm(m_numElemConstraints*m_numElemConstraints, 0.0); std::vector ss(m_numElemConstraints, 0.0); std::vector 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 { } } - diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index f2b4a3fd0..ba5854e7a 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -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 numPhSpecies(NPhase, 0); + std::vector 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 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 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; + } + + /*****************************************************************************/ } diff --git a/Cantera/src/equil/vcs_prob.h b/Cantera/src/equil/vcs_prob.h index 0b6e09312..70905a50b 100644 --- a/Cantera/src/equil/vcs_prob.h +++ b/Cantera/src/equil/vcs_prob.h @@ -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; diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 09eee6c83..c272ba0e1 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 878c977ab..bf0442b67 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -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 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 diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index caa209a76..a3e226252 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -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) {