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

[Equil] Simplify initialization of VCS_SOLVE

Browse source
This commit is contained in:
Ray Speth 2017-08-17 17:39:00 -04:00
parent 8522095dea
commit bba0d8edf0
16 changed files with 865 additions and 1330 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

24
include/cantera/equil/vcs_MultiPhaseEquil.h
View file

@ -15,26 +15,6 @@
namespace Cantera namespace Cantera
{ {
//! Translate a MultiPhase object into a VCS_SOLVE problem definition object
/*!
* @param mphase MultiPhase object that is the source for all of the information
* @param vsolve VCS_SOLVE problem definition that gets all of the information
*
* Note, both objects share the underlying ThermoPhase objects. So, neither can
* be const objects.
*/
int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve);
//! Translate a MultiPhase information into a VCS_SOLVE problem definition object
/*!
* This version updates the problem statement information only. All species and
* phase definitions remain the same.
*
* @param mphase MultiPhase object that is the source for all of the information
* @param vsolve VCS_SOLVE problem definition that gets all of the information
*/
int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve);
//! %Cantera's Interface to the Multiphase chemical equilibrium solver. //! %Cantera's Interface to the Multiphase chemical equilibrium solver.
/*! /*!
* Class vcs_MultiPhaseEquil is designed to be used to set a mixture containing * Class vcs_MultiPhaseEquil is designed to be used to set a mixture containing
@ -243,10 +223,6 @@ public:
*/ */
void reportCSV(const std::string& reportFile); void reportCSV(const std::string& reportFile);
// Friend functions
friend int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve);
friend int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve);
protected: protected:
//! Vector that takes into account of the current sorting of the species //! Vector that takes into account of the current sorting of the species
/*! /*!

9
include/cantera/equil/vcs_defs.h
View file

@ -25,15 +25,6 @@ namespace Cantera
#define VCS_MP_FAIL -6 #define VCS_MP_FAIL -6
//@} //@}
/*!
* @name Type of the underlying equilibrium solve
* @{
*/
//! Current, it is always done holding T and P constant.
#define VCS_PROBTYPE_TP 0
//@}
/*! /*!
* @name Sizes of Phases and Cutoff Mole Numbers * @name Sizes of Phases and Cutoff Mole Numbers
* *

144
include/cantera/equil/vcs_solve.h
View file

@ -46,20 +46,9 @@ class VCS_SOLVE
public: public:
//! Initialize the sizes within the VCS_SOLVE object //! Initialize the sizes within the VCS_SOLVE object
/*! /*!
* This resizes all of the internal arrays within the object. This routine * This resizes all of the internal arrays within the object.
* operates in two modes. If all of the parameters are the same as it
* currently exists in the object, nothing is done by this routine; a quick
* exit is carried out and all of the data in the object persists.
*
* If any of the parameters are different than currently exists in the
* object, then all of the data in the object must be redone. It may not be
* zeroed, but it must be redone.
*
* @param nspecies0 Number of species within the object
* @param nelements Number of element constraints within the problem
* @param nphase0 Number of phases defined within the problem.
*/ */
VCS_SOLVE(size_t nspecies0, size_t nelements, size_t nphase0); VCS_SOLVE(MultiPhase* mphase, int printLvl=0);
~VCS_SOLVE(); ~VCS_SOLVE();
@ -67,16 +56,6 @@ public:
/*! /*!
* This is the main interface routine to the equilibrium solver * This is the main interface routine to the equilibrium solver
* *
* @param ifunc Determines the operation to be done: Valid values:
* 0 -> Solve a new problem by initializing structures first. An
* initial estimate may or may not have been already
* determined. This is indicated in the VCS_SOLVE structure.
* 1 -> The problem has already been initialized and set up. We
* call this routine to resolve it using the problem
* statement and solution estimate contained in the VCS_SOLVE
* structure.
* 2 -> Don't solve a problem. Destroy all the private
* structures.
* @param ipr Printing of results * @param ipr Printing of results
* ipr = 1 -> Print problem statement and final results to * ipr = 1 -> Print problem statement and final results to
* standard output * standard output
@ -87,7 +66,7 @@ public:
* @return nonzero value: failure to solve the problem at hand. zero : * @return nonzero value: failure to solve the problem at hand. zero :
* success * success
*/ */
int vcs(int ifunc, int ipr, int ip1, int maxit); int vcs(int ipr, int ip1, int maxit);
//! Main routine that solves for equilibrium at constant T and P using a //! Main routine that solves for equilibrium at constant T and P using a
//! variant of the VCS method //! variant of the VCS method
@ -530,10 +509,6 @@ public:
//! consistent with what's needed for solution. It is called one time for //! consistent with what's needed for solution. It is called one time for
//! each new problem structure definition. //! 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 problem structure refers to:
* *
* - the number and identity of the species. * - the number and identity of the species.
@ -552,42 +527,18 @@ public:
* done so. During this process the order of the species is changed * done so. During this process the order of the species is changed
* in the private data structure. All references to the species * in the private data structure. All references to the species
* properties must employ the ind[] index vector. * properties must employ the ind[] index vector.
* 4. Initialization of arrays to zero.
* 5. Check to see if the problem is well posed (If all the element
* abundances are zero, the algorithm will fail)
* *
* @param printLvl Print level of the routine * @param printLvl Print level of the routine
* @return VCS_SUCCESS = everything went OK
*/
int vcs_prep_oneTime(int printLvl);
//! 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 every
* invocation of the vcs_solve() except for the cleanup invocation.
*
* Tasks:
* 1. Initialization of arrays to zero.
*
* @return * @return
* VCS_SUCCESS = everything went OK; * VCS_SUCCESS = everything went OK;
* VCS_PUB_BAD = There is an irreconcilable difference in the * VCS_PUB_BAD = There is an irreconcilable difference in the
* public data structure from when the problem was * public data structure from when the problem was
* initially set up. * initially set up.
*/ */
int vcs_prep(); int vcs_prep(int printLvl);
//! 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.
*
* @return If true, the problem is well-posed. If false, the problem
* is not well posed.
*/
bool vcs_wellPosed();
//! Rearrange the constraint equations represented by the Formula //! Rearrange the constraint equations represented by the Formula
//! Matrix so that the operational ones are in the front //! Matrix so that the operational ones are in the front
@ -817,12 +768,6 @@ public:
//! Fully specify the problem to be solved //! Fully specify the problem to be solved
int vcs_prob_specifyFully(); int vcs_prob_specifyFully();
//! Specify the problem to be solved, incrementally
/*!
* It's assumed we are solving the same problem.
*/
int vcs_prob_specify();
private: private:
//! Zero out the concentration of a species. //! Zero out the concentration of a species.
/*! /*!
@ -1063,20 +1008,6 @@ private:
public: public:
//! @{ Variables moved from VCS_PROB //! @{ Variables moved from VCS_PROB
//! Problem type. I.e., the identity of what is held constant. Currently, T
//! and P are held constant, and this input is ignored
int prob_type;
//! Total number of species in the problems
size_t nspecies;
//! Number of element constraints in the equilibrium problem
size_t ne;
//! Number of element constraints used to size data structures
//! involving elements
size_t NE0;
//! Number of phases in the problem
size_t NPhase;
//! Vector of chemical potentials of the species. This is a calculated //! Vector of chemical potentials of the species. This is a calculated
//! output quantity. length = number of species. //! output quantity. length = number of species.
vector_fp m_gibbsSpecies; vector_fp m_gibbsSpecies;
@ -1119,42 +1050,6 @@ public:
*/ */
vector_int SpeciesUnknownType; vector_int SpeciesUnknownType;
//! Temperature (Kelvin)
/*!
* Specification of the temperature for the equilibrium problem
*/
double T;
//! Pressure
double PresPA;
//! Volume of the entire system
/*!
* Note, this is an output variable atm
*/
double Vol;
//! Partial Molar Volumes of species
/*!
* This is a calculated vector, calculated from w[].
* length number of species.
*/
vector_fp VolPM;
//! Specification of the initial estimate method
/*!
* * 0: user estimate
* * 1: user estimate if satisifies elements
* * -1: machine estimate
*/
int iest;
//! Tolerance requirement for major species
double tolmaj;
//! Tolerance requirement for minor species
double tolmin;
//! Mapping between the species and the phases //! Mapping between the species and the phases
std::vector<size_t> PhaseID; std::vector<size_t> PhaseID;
@ -1177,20 +1072,11 @@ public:
//! Array of phase structures //! Array of phase structures
std::vector<vcs_VolPhase*> VPhaseList; std::vector<vcs_VolPhase*> VPhaseList;
// String containing the title of the run
std::string Title;
//! Vector of pointers to thermo structures which identify the model and //! Vector of pointers to thermo structures which identify the model and
//! parameters for evaluating the thermodynamic functions for that //! parameters for evaluating the thermodynamic functions for that
//! particular species //! particular species
std::vector<VCS_SPECIES_THERMO*> SpeciesThermo; std::vector<VCS_SPECIES_THERMO*> SpeciesThermo;
//! Number of iterations. This is an output variable
int m_Iterations;
//! Number of basis optimizations used. This is an output variable.
int m_NumBasisOptimizations;
//! Print level for print routines //! Print level for print routines
int m_printLvl; int m_printLvl;
@ -1258,20 +1144,11 @@ public:
void reportCSV(const std::string& reportFile); void reportCSV(const std::string& reportFile);
//! value of the number of species used to size data structures
size_t NSPECIES0;
//! value of the number of phases used to size data structures
size_t NPHASE0;
//! Total number of species in the problems //! Total number of species in the problems
size_t m_numSpeciesTot; size_t m_nsp;
//! Number of element constraints in the problem //! Number of element constraints in the problem
/*! size_t m_nelem;
* This is typically equal to the number of elements in the problem
*/
size_t m_numElemConstraints;
//! Number of components calculated for the problem //! Number of components calculated for the problem
size_t m_numComponents; size_t m_numComponents;
@ -1620,9 +1497,6 @@ public:
//! Array of Phase Structures. Length = number of phases. //! Array of Phase Structures. Length = number of phases.
std::vector<vcs_VolPhase*> m_VolPhaseList; std::vector<vcs_VolPhase*> m_VolPhaseList;
//! String containing the title of the run
std::string m_title;
//! This specifies the current state of units for the Gibbs free energy //! This specifies the current state of units for the Gibbs free energy
//! properties in the program. //! properties in the program.
/*! /*!

444
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -22,20 +22,10 @@ using namespace std;
namespace Cantera namespace Cantera
{ {
vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(MultiPhase* mix, int printLvl) : vcs_MultiPhaseEquil::vcs_MultiPhaseEquil(MultiPhase* mix, int printLvl) :
m_mix(0), m_mix(mix),
m_printLvl(printLvl), m_printLvl(printLvl),
m_vsolve(mix->nSpecies(), mix->nElements(), mix->nPhases()) m_vsolve(mix, printLvl)
{ {
m_mix = mix;
m_vsolve.m_printLvl = m_printLvl;
m_vsolve.m_mix = m_mix;
// Work out the details of the VCS_SOLVE construction and Transfer the
// current problem to the VCS_SOLVE object
int res = vcs_Cantera_to_vprob(mix, &m_vsolve);
if (res != 0) {
plogf("problems\n");
}
} }
int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
@ -430,20 +420,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
clockWC tickTock; clockWC tickTock;
m_printLvl = printLvl; m_printLvl = printLvl;
m_vsolve.m_printLvl = printLvl; m_vsolve.m_printLvl = printLvl;
m_vsolve.m_doEstimateEquil = estimateEquil;
// 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_vsolve);
if (res != 0) {
plogf("problems\n");
}
// Set the estimation technique
if (estimateEquil) {
m_vsolve.iest = estimateEquil;
} else {
m_vsolve.iest = 0;
}
// Check obvious bounds on the temperature and pressure NOTE, we may want to // 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. // do more here with the real bounds given by the ThermoPhase objects.
@ -456,10 +433,6 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
"Pressure less than zero on input"); "Pressure less than zero on input");
} }
// Print out the problem specification from the point of
// view of the vprob object.
m_vsolve.prob_report(m_printLvl);
//! Call the thermo Program //! Call the thermo Program
int ip1 = m_printLvl; int ip1 = m_printLvl;
int ipr = std::max(0, m_printLvl-1); int ipr = std::max(0, m_printLvl-1);
@ -468,7 +441,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
} else { } else {
ip1 = 0; ip1 = 0;
} }
int iSuccess = m_vsolve.vcs(0, ipr, ip1, maxit); int iSuccess = m_vsolve.vcs(ipr, ip1, maxit);
// Transfer the information back to the MultiPhase object. Note we don't // Transfer the information back to the MultiPhase object. Note we don't
// just call setMoles, because some multispecies solution phases may be // just call setMoles, because some multispecies solution phases may be
@ -477,7 +450,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
// about likely reemergent states. // about likely reemergent states.
m_mix->uploadMoleFractionsFromPhases(); m_mix->uploadMoleFractionsFromPhases();
size_t kGlob = 0; size_t kGlob = 0;
for (size_t ip = 0; ip < m_vsolve.NPhase; ip++) { for (size_t ip = 0; ip < m_vsolve.m_numPhases; ip++) {
double phaseMole = 0.0; double phaseMole = 0.0;
ThermoPhase& tref = m_mix->phase(ip); ThermoPhase& tref = m_mix->phase(ip);
for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) { for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) {
@ -493,8 +466,8 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
plogf("\nVCS FAILED TO CONVERGE!\n"); plogf("\nVCS FAILED TO CONVERGE!\n");
} }
plogf("\n"); plogf("\n");
plogf("Temperature = %g Kelvin\n", m_vsolve.T); plogf("Temperature = %g Kelvin\n", m_vsolve.m_temperature);
plogf("Pressure = %g Pa\n", m_vsolve.PresPA); plogf("Pressure = %g Pa\n", m_vsolve.m_pressurePA);
plogf("\n"); plogf("\n");
plogf("----------------------------------------" plogf("----------------------------------------"
"---------------------\n"); "---------------------\n");
@ -502,7 +475,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
plogf(" Mole_Fraction Chem_Potential (J/kmol)\n"); plogf(" Mole_Fraction Chem_Potential (J/kmol)\n");
plogf("--------------------------------------------------" plogf("--------------------------------------------------"
"-----------\n"); "-----------\n");
for (size_t i = 0; i < m_vsolve.nspecies; i++) { for (size_t i = 0; i < m_mix->nSpecies(); i++) {
plogf("%-12s", m_mix->speciesName(i)); plogf("%-12s", m_mix->speciesName(i));
if (m_vsolve.SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_vsolve.SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" %15.3e %15.3e ", 0.0, m_vsolve.mf[i]); plogf(" %15.3e %15.3e ", 0.0, m_vsolve.mf[i]);
@ -533,7 +506,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
{ {
size_t nphase = m_vsolve.NPhase; size_t nphase = m_vsolve.m_numPhases;
FILE* FP = fopen(reportFile.c_str(), "w"); FILE* FP = fopen(reportFile.c_str(), "w");
if (!FP) { if (!FP) {
@ -572,8 +545,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature()); fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature());
fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure()); fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure());
fprintf(FP,"Total Volume = %11.5g m**3\n", vol); fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
fprintf(FP,"Number Basis optimizations = %d\n", m_vsolve.m_NumBasisOptimizations); fprintf(FP,"Number Basis optimizations = %d\n", m_vsolve.m_VCount->Basis_Opts);
fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_Iterations); fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_VCount->Its);
for (size_t iphase = 0; iphase < nphase; iphase++) { for (size_t iphase = 0; iphase < nphase; iphase++) {
size_t istart = m_mix->speciesIndex(0, iphase); size_t istart = m_mix->speciesIndex(0, iphase);
@ -666,399 +639,4 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
fclose(FP); fclose(FP);
} }
// HKM -> Work on transferring the current value of the voltages into the
// equilibrium problem.
int vcs_Cantera_to_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve)
{
VCS_SPECIES_THERMO* ts_ptr = 0;
// Calculate the total number of species and phases in the problem
size_t totNumPhases = mphase->nPhases();
size_t totNumSpecies = mphase->nSpecies();
// Problem type has yet to be worked out.
vsolve->prob_type = 0;
vsolve->nspecies = totNumSpecies;
vsolve->ne = 0;
vsolve->NPhase = totNumPhases;
// Set the initial estimate to a machine generated estimate for now
// We will work out the details later.
vsolve->iest = -1;
vsolve->T = mphase->temperature();
vsolve->PresPA = mphase->pressure();
vsolve->Vol = mphase->volume();
vsolve->Title = "MultiPhase Object";
int printLvl = vsolve->m_printLvl;
// Loop over the phases, transferring pertinent information
int kT = 0;
for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
// Get the ThermoPhase object - assume volume phase
ThermoPhase* tPhase = &mphase->phase(iphase);
size_t nelem = tPhase->nElements();
// Query Cantera for the equation of state type of the current phase.
std::string eos = tPhase->type();
bool gasPhase = (eos == "IdealGas");
// Find out the number of species in the phase
size_t nSpPhase = tPhase->nSpecies();
// Find out the name of the phase
string phaseName = tPhase->name();
// Call the basic vcs_VolPhase creation routine.
// Properties set here:
// ->PhaseNum = phase number in the thermo problem
// ->GasPhase = Boolean indicating whether it is a gas phase
// ->NumSpecies = number of species in the phase
// ->TMolesInert = Inerts in the phase = 0.0 for cantera
// ->PhaseName = Name of the phase
vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase];
VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
VolPhase->m_gasPhase = gasPhase;
// Tell the vcs_VolPhase pointer about cantera
VolPhase->setPtrThermoPhase(tPhase);
VolPhase->setTotalMoles(0.0);
// Set the electric potential of the volume phase from the
// ThermoPhase object's value.
VolPhase->setElectricPotential(tPhase->electricPotential());
// Query the ThermoPhase object to find out what convention
// it uses for the specification of activity and Standard State.
VolPhase->p_activityConvention = tPhase->activityConvention();
// Assign the value of eqn of state. Handle conflicts here.
if (eos == "IdealGas") {
VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS;
} else if (eos == "ConstDensity") {
VolPhase->m_eqnState = VCS_EOS_CONSTANT;
} else if (eos == "StoichSubstance") {
VolPhase->m_eqnState = VCS_EOS_STOICH_SUB;
} else if (eos == "IdealSolidSoln") {
VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN;
} else if (eos == "Surf" || eos == "Edge") {
throw CanteraError("VCSnonideal",
"Surface/edge phase not handled yet.");
} else {
if (printLvl > 1) {
writelog("Unknown Cantera EOS to VCSnonideal: '{}'\n", eos);
}
VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
}
// Transfer all of the element information from the ThermoPhase object
// to the vcs_VolPhase object. Also decide whether we need a new charge
// neutrality element in the phase to enforce a charge neutrality
// constraint. We also decide whether this is a single species phase
// with the voltage being the independent variable setting the chemical
// potential of the electrons.
VolPhase->transferElementsFM(tPhase);
// Combine the element information in the vcs_VolPhase
// object into the vprob object.
vsolve->addPhaseElements(VolPhase);
VolPhase->setState_TP(vsolve->T, vsolve->PresPA);
vector_fp muPhase(tPhase->nSpecies(),0.0);
tPhase->getChemPotentials(&muPhase[0]);
double tMoles = 0.0;
// Loop through each species in the current phase
for (size_t k = 0; k < nSpPhase; k++) {
// Obtain the molecular weight of the species from the
// ThermoPhase object
vsolve->WtSpecies[kT] = tPhase->molecularWeight(k);
// Obtain the charges of the species from the ThermoPhase object
vsolve->Charge[kT] = tPhase->charge(k);
// Set the phaseid of the species
vsolve->PhaseID[kT] = iphase;
// Transfer the type of unknown
vsolve->SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k);
if (vsolve->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) {
// Set the initial number of kmoles of the species
// and the mole fraction vector
vsolve->w[kT] = mphase->speciesMoles(kT);
tMoles += vsolve->w[kT];
vsolve->mf[kT] = mphase->moleFraction(kT);
} else if (vsolve->SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
vsolve->w[kT] = tPhase->electricPotential();
vsolve->mf[kT] = mphase->moleFraction(kT);
} else {
throw CanteraError(" vcs_Cantera_to_vsolve() ERROR",
"Unknown species type: {}", vsolve->SpeciesUnknownType[kT]);
}
// transfer chemical potential vector
vsolve->m_gibbsSpecies[kT] = muPhase[k];
// Transfer the species information from the
// volPhase structure to the VPROB structure
// This includes:
// FormulaMatrix[][]
// VolPhase->IndSpecies[]
vsolve->addOnePhaseSpecies(VolPhase, k, kT);
// Get a pointer to the thermo object
ts_ptr = vsolve->SpeciesThermo[kT];
// Fill in the vcs_SpeciesProperty structure
vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
sProp->NumElements = vsolve->ne;
sProp->SpName = mphase->speciesName(kT);
sProp->SpeciesThermo = ts_ptr;
sProp->WtSpecies = tPhase->molecularWeight(k);
sProp->FormulaMatrixCol.resize(vsolve->ne, 0.0);
for (size_t e = 0; e < vsolve->ne; e++) {
sProp->FormulaMatrixCol[e] = vsolve->FormulaMatrix(kT,e);
}
sProp->Charge = tPhase->charge(k);
sProp->SurfaceSpecies = false;
sProp->VolPM = 0.0;
// Transfer the thermo specification of the species
// vsolve->SpeciesThermo[]
// Add lookback connectivity into the thermo object first
ts_ptr->IndexPhase = iphase;
ts_ptr->IndexSpeciesPhase = k;
ts_ptr->OwningPhase = VolPhase;
// get a reference to the Cantera species thermo.
MultiSpeciesThermo& sp = tPhase->speciesThermo();
int spType = sp.reportType(k);
if (spType == SIMPLE) {
double c[4];
double minTemp, maxTemp, refPressure;
sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure);
ts_ptr->SS0_Model = VCS_SS0_CONSTANT;
ts_ptr->SS0_T0 = c[0];
ts_ptr->SS0_H0 = c[1];
ts_ptr->SS0_S0 = c[2];
ts_ptr->SS0_Cp0 = c[3];
if (gasPhase) {
ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS;
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS;
} else {
ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT;
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
}
} else {
if (vsolve->m_printLvl > 2) {
plogf("vcs_Cantera_convert: Species Type %d not known \n",
spType);
}
ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED;
ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED;
}
// Transfer the Volume Information -> NEEDS WORK
if (gasPhase) {
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_IDEALGAS;
ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0;
} else {
vector_fp phaseTermCoeff(nSpPhase, 0.0);
int nCoeff;
tPhase->getParameters(nCoeff, &phaseTermCoeff[0]);
ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
ts_ptr->SSStar_Vol0 = phaseTermCoeff[k];
}
kT++;
}
// Now go back through the species in the phase and assign a valid mole
// fraction to all phases, even if the initial estimate of the total
// number of moles is zero.
if (tMoles > 0.0) {
for (size_t k = 0; k < nSpPhase; k++) {
size_t kTa = VolPhase->spGlobalIndexVCS(k);
vsolve->mf[kTa] = vsolve->w[kTa] / tMoles;
}
} else {
// Perhaps, we could do a more sophisticated treatment below.
// But, will start with this.
for (size_t k = 0; k < nSpPhase; k++) {
size_t kTa = VolPhase->spGlobalIndexVCS(k);
vsolve->mf[kTa]= 1.0 / (double) nSpPhase;
}
}
VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vsolve->w[0]);
// Now, calculate a sample naught Gibbs free energy calculation
// at the specified temperature.
for (size_t k = 0; k < nSpPhase; k++) {
vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
ts_ptr = sProp->SpeciesThermo;
ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant;
ts_ptr->SS0_TSave = vsolve->T;
}
}
// Transfer initial element abundances to the vprob object.
// We have to find the mapping index from one to the other
vsolve->gai.resize(vsolve->ne, 0.0);
vsolve->set_gai();
// Printout the species information: PhaseID's and mole nums
if (vsolve->m_printLvl > 1) {
writeline('=', 80, true, true);
writeline('=', 16, false);
plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
writeline('=', 20);
writeline('=', 80);
plogf(" Phase IDs of species\n");
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_kMols\n");
for (size_t i = 0; i < vsolve->nspecies; i++) {
size_t iphase = vsolve->PhaseID[i];
vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase];
plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase,
VolPhase->PhaseName.c_str());
if (vsolve->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Volts = %-10.5g\n", vsolve->w[i]);
} else {
plogf(" %-10.5g\n", vsolve->w[i]);
}
}
// Printout of the Phase structure information
writeline('-', 80, true, true);
plogf(" Information about phases\n");
plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
plogf(" TMolesInert Tmoles(kmol)\n");
for (size_t iphase = 0; iphase < vsolve->NPhase; iphase++) {
vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase];
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
VolPhase->m_gasPhase, VolPhase->eos_name(),
VolPhase->nSpecies(), VolPhase->totalMolesInert());
plogf("%16e\n", VolPhase->totalMoles());
}
writeline('=', 80, true, true);
writeline('=', 16, false);
plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
writeline('=', 20);
writeline('=', 80);
plogf("\n");
}
return VCS_SUCCESS;
}
int vcs_Cantera_update_vprob(MultiPhase* mphase, VCS_SOLVE* vsolve)
{
size_t totNumPhases = mphase->nPhases();
size_t kT = 0;
vector_fp tmpMoles;
// Problem type has yet to be worked out.
vsolve->prob_type = 0;
// Whether we have an estimate or not gets overwritten on
// the call to the equilibrium solver.
vsolve->iest = -1;
vsolve->T = mphase->temperature();
vsolve->PresPA = mphase->pressure();
vsolve->Vol = mphase->volume();
for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
ThermoPhase* tPhase = &mphase->phase(iphase);
vcs_VolPhase* volPhase = vsolve->VPhaseList[iphase];
// Set the electric potential of the volume phase from the
// ThermoPhase object's value.
volPhase->setElectricPotential(tPhase->electricPotential());
volPhase->setState_TP(vsolve->T, vsolve->PresPA);
vector_fp muPhase(tPhase->nSpecies(),0.0);
tPhase->getChemPotentials(&muPhase[0]);
// Loop through each species in the current phase
size_t nSpPhase = tPhase->nSpecies();
tmpMoles.resize(nSpPhase);
for (size_t k = 0; k < nSpPhase; k++) {
tmpMoles[k] = mphase->speciesMoles(kT);
vsolve->w[kT] = mphase->speciesMoles(kT);
vsolve->mf[kT] = mphase->moleFraction(kT);
// transfer chemical potential vector
vsolve->m_gibbsSpecies[kT] = muPhase[k];
kT++;
}
if (volPhase->phiVarIndex() != npos) {
size_t kphi = volPhase->phiVarIndex();
size_t kglob = volPhase->spGlobalIndexVCS(kphi);
vsolve->w[kglob] = tPhase->electricPotential();
}
volPhase->setMolesFromVCS(VCS_STATECALC_OLD, &vsolve->w[0]);
if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) {
volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS);
} else if (volPhase->totalMoles() > 0.0) {
volPhase->setExistence(VCS_PHASE_EXIST_YES);
} else {
volPhase->setExistence(VCS_PHASE_EXIST_NO);
}
}
// Transfer initial element abundances to the vprob object. Put them in the
// front of the object. There may be more constraints than there are
// elements. But, we know the element abundances are in the front of the
// vector.
vsolve->set_gai();
// Printout the species information: PhaseID's and mole nums
if (vsolve->m_printLvl > 1) {
writeline('=', 80, true, true);
writeline('=', 20, false);
plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
writeline('=', 20);
writeline('=', 80);
plogf("\n");
plogf(" Phase IDs of species\n");
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_kMols\n");
for (size_t i = 0; i < vsolve->nspecies; i++) {
size_t iphase = vsolve->PhaseID[i];
vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase];
plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase,
VolPhase->PhaseName.c_str());
if (vsolve->SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Volts = %-10.5g\n", vsolve->w[i]);
} else {
plogf(" %-10.5g\n", vsolve->w[i]);
}
}
// Printout of the Phase structure information
writeline('-', 80, true, true);
plogf(" Information about phases\n");
plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
plogf(" TMolesInert Tmoles(kmol)\n");
for (size_t iphase = 0; iphase < vsolve->NPhase; iphase++) {
vcs_VolPhase* VolPhase = vsolve->VPhaseList[iphase];
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
VolPhase->m_gasPhase, VolPhase->eos_name(),
VolPhase->nSpecies(), VolPhase->totalMolesInert());
plogf("%16e\n", VolPhase->totalMoles());
}
writeline('=', 80, true, true);
writeline('=', 20, false);
plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
writeline('=', 20);
writeline('=', 80);
plogf("\n");
}
return VCS_SUCCESS;
}
} }

2
src/equil/vcs_TP.cpp
View file

@ -60,7 +60,7 @@ int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
void VCS_SOLVE::vcs_fePrep_TP() void VCS_SOLVE::vcs_fePrep_TP()
{ {
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
// For single species phases, initialize the chemical potential with the // For single species phases, initialize the chemical potential with the
// value of the standard state chemical potential. This value doesn't // value of the standard state chemical potential. This value doesn't
// change during the calculation // change during the calculation

48
src/equil/vcs_elem.cpp
View file

@ -16,9 +16,9 @@ namespace Cantera
{ {
void VCS_SOLVE::vcs_elab() void VCS_SOLVE::vcs_elab()
{ {
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
m_elemAbundances[j] = 0.0; m_elemAbundances[j] = 0.0;
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_elemAbundances[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; m_elemAbundances[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i];
} }
@ -30,7 +30,7 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound)
{ {
size_t top = m_numComponents; size_t top = m_numComponents;
if (ibound) { if (ibound) {
top = m_numElemConstraints; top = m_nelem;
} }
for (size_t i = 0; i < top; ++i) { for (size_t i = 0; i < top; ++i) {
@ -50,7 +50,7 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound)
// addition of terms. We are limited to 13 digits of finite // addition of terms. We are limited to 13 digits of finite
// arithmetic accuracy. // arithmetic accuracy.
bool multisign = false; bool multisign = false;
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
double eval = m_formulaMatrix(kspec,i); double eval = m_formulaMatrix(kspec,i);
if (eval < 0.0) { if (eval < 0.0) {
multisign = true; multisign = true;
@ -84,9 +84,9 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound)
void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase) void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase)
{ {
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
elemAbundPhase[j] = 0.0; elemAbundPhase[j] = 0.0;
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && m_phaseID[i] == iphase) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && m_phaseID[i] == iphase) {
elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i]; elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i];
} }
@ -101,29 +101,29 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
vector_fp ga_save(m_elemAbundances); vector_fp ga_save(m_elemAbundances);
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
plogf(" --- vcsc_elcorr: Element abundances correction routine"); plogf(" --- vcsc_elcorr: Element abundances correction routine");
if (m_numElemConstraints != m_numComponents) { if (m_nelem != m_numComponents) {
plogf(" (m_numComponents != m_numElemConstraints)"); plogf(" (m_numComponents != m_nelem)");
} }
plogf("\n"); plogf("\n");
} }
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
} }
double l2before = 0.0; double l2before = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
l2before += x[i] * x[i]; l2before += x[i] * x[i];
} }
l2before = sqrt(l2before/m_numElemConstraints); l2before = sqrt(l2before/m_nelem);
// Special section to take out single species, single component, // Special section to take out single species, single component,
// moles. These are species which have non-zero entries in the // moles. These are species which have non-zero entries in the
// formula matrix, and no other species have zero values either. // formula matrix, and no other species have zero values either.
bool changed = false; bool changed = false;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
int numNonZero = 0; int numNonZero = 0;
bool multisign = false; bool multisign = false;
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double eval = m_formulaMatrix(kspec,i); double eval = m_formulaMatrix(kspec,i);
if (eval < 0.0) { if (eval < 0.0) {
@ -136,7 +136,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
} }
if (!multisign) { if (!multisign) {
if (numNonZero < 2) { if (numNonZero < 2) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double eval = m_formulaMatrix(kspec,i); double eval = m_formulaMatrix(kspec,i);
if (eval > 0.0) { if (eval > 0.0) {
@ -159,7 +159,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
} }
if (numCompNonZero == 1) { if (numCompNonZero == 1) {
double diff = m_elemAbundancesGoal[i]; double diff = m_elemAbundancesGoal[i];
for (size_t kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = m_numComponents; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double eval = m_formulaMatrix(kspec,i); double eval = m_formulaMatrix(kspec,i);
diff -= eval * m_molNumSpecies_old[kspec]; diff -= eval * m_molNumSpecies_old[kspec];
@ -183,10 +183,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
// Note, also we can do this over ne, the number of elements, not just the // Note, also we can do this over ne, the number of elements, not just the
// number of components. // number of components.
changed = false; changed = false;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
int elType = m_elType[i]; int elType = m_elType[i];
if (elType == VCS_ELEM_TYPE_ABSPOS) { if (elType == VCS_ELEM_TYPE_ABSPOS) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double atomComp = m_formulaMatrix(kspec,i); double atomComp = m_formulaMatrix(kspec,i);
if (atomComp > 0.0) { if (atomComp > 0.0) {
@ -238,7 +238,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
} }
} }
solve(A, x, 1, m_numElemConstraints); solve(A, x, 1, m_nelem);
// Now apply the new direction without creating negative species. // Now apply the new direction without creating negative species.
double par = 0.5; double par = 0.5;
@ -288,7 +288,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
// Try some ad hoc procedures for fixing the problem // Try some ad hoc procedures for fixing the problem
if (retn >= 2) { if (retn >= 2) {
// First find a species whose adjustment is a win-win situation. // First find a species whose adjustment is a win-win situation.
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
continue; continue;
} }
@ -348,7 +348,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
goto L_CLEANUP; goto L_CLEANUP;
} }
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY || if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY ||
(m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) {
for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) { for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
@ -374,7 +374,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
// For electron charges element types, we try positive deltas in the species // For electron charges element types, we try positive deltas in the species
// concentrations to match the desired electron charge exactly. // concentrations to match the desired electron charge exactly.
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i]; double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i];
if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) { if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) {
bool useZeroed = true; bool useZeroed = true;
@ -418,14 +418,14 @@ L_CLEANUP:
; ;
vcs_tmoles(); vcs_tmoles();
double l2after = 0.0; double l2after = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2); l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2);
} }
l2after = sqrt(l2after/m_numElemConstraints); l2after = sqrt(l2after/m_nelem);
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
plogf(" --- Elem_Abund: Correct Initial " plogf(" --- Elem_Abund: Correct Initial "
" Final\n"); " Final\n");
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
plogf(" --- "); plogf(" --- ");
plogf("%-2.2s", m_elementName[i]); plogf("%-2.2s", m_elementName[i]);
plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]); plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]);

12
src/equil/vcs_elem_rearrange.cpp
View file

@ -33,7 +33,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
double test = -1.0E10; double test = -1.0E10;
while (lindep) { while (lindep) {
lindep = false; lindep = false;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
test -= 1.0; test -= 1.0;
aw[i] = m_elemAbundancesGoal[i]; aw[i] = m_elemAbundancesGoal[i];
if (test == aw[i]) { if (test == aw[i]) {
@ -53,14 +53,14 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
while (true) { while (true) {
// Search the remaining part of the mole fraction vector, AW, for // Search the remaining part of the mole fraction vector, AW, for
// the largest remaining species. Return its identity in K. // the largest remaining species. Return its identity in K.
k = m_numElemConstraints; k = m_nelem;
for (size_t ielem = jr; ielem < m_numElemConstraints; ielem++) { for (size_t ielem = jr; ielem < m_nelem; ielem++) {
if (m_elementActive[ielem] && aw[ielem] != test) { if (m_elementActive[ielem] && aw[ielem] != test) {
k = ielem; k = ielem;
break; break;
} }
} }
if (k == m_numElemConstraints) { if (k == m_nelem) {
throw CanteraError("vcs_elem_rearrange", throw CanteraError("vcs_elem_rearrange",
"Shouldn't be here. Algorithm misfired."); "Shouldn't be here. Algorithm misfired.");
} }
@ -136,7 +136,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
if (ipos == jpos) { if (ipos == jpos) {
return; return;
} }
AssertThrowMsg(ipos < m_numElemConstraints && jpos < m_numElemConstraints, AssertThrowMsg(ipos < m_nelem && jpos < m_nelem,
"vcs_switch_elem_pos", "vcs_switch_elem_pos",
"inappropriate args: {} {}", ipos, jpos); "inappropriate args: {} {}", ipos, jpos);
@ -158,7 +158,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]); std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]);
std::swap(m_elType[ipos], m_elType[jpos]); std::swap(m_elType[ipos], m_elType[jpos]);
std::swap(m_elementActive[ipos], m_elementActive[jpos]); std::swap(m_elementActive[ipos], m_elementActive[jpos]);
for (size_t j = 0; j < m_numSpeciesTot; ++j) { for (size_t j = 0; j < m_nsp; ++j) {
std::swap(m_formulaMatrix(j,ipos), m_formulaMatrix(j,jpos)); std::swap(m_formulaMatrix(j,ipos), m_formulaMatrix(j,jpos));
} }
std::swap(m_elementName[ipos], m_elementName[jpos]); std::swap(m_elementName[ipos], m_elementName[jpos]);

33
src/equil/vcs_inest.cpp
View file

@ -19,7 +19,6 @@ static char pprefix[20] = " --- vcs_inest: ";
void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
double* const ss, double test) double* const ss, double test)
{ {
size_t nspecies = m_numSpeciesTot;
size_t nrxn = m_numRxnTot; size_t nrxn = m_numRxnTot;
// CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB AND // CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB AND
@ -30,7 +29,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n", plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
pprefix); pprefix);
plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix); plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
plogf("%s ", pprefix); plogf("%s ", pprefix);
plogf("%-12.12s", m_speciesName[kspec]); plogf("%-12.12s", m_speciesName[kspec]);
plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]); plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]);
@ -38,10 +37,10 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
plogf("%s Element Abundance Agreement returned from linear " plogf("%s Element Abundance Agreement returned from linear "
"programming (vcs_inest initial guess):\n", pprefix); "programming (vcs_inest initial guess):\n", pprefix);
plogf("%s Element Goal Actual\n", pprefix); plogf("%s Element Goal Actual\n", pprefix);
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
if (m_elementActive[j]) { if (m_elementActive[j]) {
double tmp = 0.0; double tmp = 0.0;
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
tmp += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; tmp += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec];
} }
plogf("%s ", pprefix); plogf("%s ", pprefix);
@ -55,7 +54,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
// Make sure all species have positive definite mole numbers Set voltages to // Make sure all species have positive definite mole numbers Set voltages to
// zero for now, until we figure out what to do // zero for now, until we figure out what to do
m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0); m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0);
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_molNumSpecies_old[kspec] <= 0.0) { if (m_molNumSpecies_old[kspec] <= 0.0) {
// HKM Should eventually include logic here for non SS phases // HKM Should eventually include logic here for non SS phases
@ -113,7 +112,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
} }
vcs_deltag(0, true, VCS_STATECALC_NEW); vcs_deltag(0, true, VCS_STATECALC_NEW);
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
plogf("%s", pprefix); plogf("%s", pprefix);
plogf("%-12.12s", m_speciesName[kspec]); plogf("%-12.12s", m_speciesName[kspec]);
if (kspec < m_numComponents) { if (kspec < m_numComponents) {
@ -168,7 +167,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
} }
} }
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf("%sdirection (", pprefix); plogf("%sdirection (", pprefix);
plogf("%-12.12s", m_speciesName[kspec]); plogf("%-12.12s", m_speciesName[kspec]);
@ -212,7 +211,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
m_deltaMolNumSpecies[kspec] = 0.0; m_deltaMolNumSpecies[kspec] = 0.0;
} }
} }
for (size_t kspec = m_numComponents; kspec < nspecies; ++kspec) { for (size_t kspec = m_numComponents; kspec < m_nsp; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE && if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE &&
m_deltaMolNumSpecies[kspec] != 0.0) { m_deltaMolNumSpecies[kspec] != 0.0) {
m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par; m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par;
@ -228,9 +227,9 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
// CONVERGENCE FORCING SECTION // CONVERGENCE FORCING SECTION
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_dfe(VCS_STATECALC_OLD, 0, 0, nspecies); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp);
double s = 0.0; double s = 0.0;
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec]; s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec];
} }
if (s == 0.0) { if (s == 0.0) {
@ -269,7 +268,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
plogf("%s Final Mole Numbers produced by inest:\n", plogf("%s Final Mole Numbers produced by inest:\n",
pprefix); pprefix);
plogf("%s SPECIES MOLE_NUMBER\n", pprefix); plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
for (size_t kspec = 0; kspec < nspecies; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
plogf("%s %-12.12s %g\n", plogf("%s %-12.12s %g\n",
pprefix, m_speciesName[kspec], m_molNumSpecies_old[kspec]); pprefix, m_speciesName[kspec], m_molNumSpecies_old[kspec]);
} }
@ -298,10 +297,10 @@ int VCS_SOLVE::vcs_inest_TP()
} }
// temporary space for usage in this routine and in subroutines // temporary space for usage in this routine and in subroutines
vector_fp sm(m_numElemConstraints*m_numElemConstraints, 0.0); vector_fp sm(m_nelem*m_nelem, 0.0);
vector_fp ss(m_numElemConstraints, 0.0); vector_fp ss(m_nelem, 0.0);
vector_fp sa(m_numElemConstraints, 0.0); vector_fp sa(m_nelem, 0.0);
vector_fp aw(m_numSpeciesTot+ m_numElemConstraints, 0.0); vector_fp aw(m_nsp + m_nelem, 0.0);
// Go get the estimate of the solution // Go get the estimate of the solution
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
@ -343,7 +342,7 @@ int VCS_SOLVE::vcs_inest_TP()
plogf("%sElement Abundances RANGE ERROR\n", pprefix); plogf("%sElement Abundances RANGE ERROR\n", pprefix);
plogf("%s - Initial guess satisfies NC=%d element abundances, " plogf("%s - Initial guess satisfies NC=%d element abundances, "
"BUT not NE=%d element abundances\n", pprefix, "BUT not NE=%d element abundances\n", pprefix,
m_numComponents, m_numElemConstraints); m_numComponents, m_nelem);
} }
} }
} }
@ -355,7 +354,7 @@ int VCS_SOLVE::vcs_inest_TP()
plogf("%sElement Abundances RANGE ERROR\n", pprefix); plogf("%sElement Abundances RANGE ERROR\n", pprefix);
plogf("%s - Initial guess satisfies NC=%d element abundances, " plogf("%s - Initial guess satisfies NC=%d element abundances, "
"BUT not NE=%d element abundances\n", pprefix, "BUT not NE=%d element abundances\n", pprefix,
m_numComponents, m_numElemConstraints); m_numComponents, m_nelem);
} }
} }
} }

14
src/equil/vcs_nondim.cpp
View file

@ -19,7 +19,7 @@ void VCS_SOLVE::vcs_nondim_TP()
if (m_unitsState == VCS_DIMENSIONAL_G) { if (m_unitsState == VCS_DIMENSIONAL_G) {
m_unitsState = VCS_NONDIMENSIONAL_G; m_unitsState = VCS_NONDIMENSIONAL_G;
double tf = 1.0 / (GasConstant * m_temperature); double tf = 1.0 / (GasConstant * m_temperature);
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
// Modify the standard state and total chemical potential data, // Modify the standard state and total chemical potential data,
// FF(I), to make it dimensionless, i.e., mu / RT. Thus, we may // FF(I), to make it dimensionless, i.e., mu / RT. Thus, we may
// divide it by the temperature. // divide it by the temperature.
@ -37,7 +37,7 @@ void VCS_SOLVE::vcs_nondim_TP()
// Then add in the total moles of elements that are goals. Either one or // Then add in the total moles of elements that are goals. Either one or
// the other is specified here. // the other is specified here.
double esum = 0.0; double esum = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
esum += fabs(m_elemAbundancesGoal[i]); esum += fabs(m_elemAbundancesGoal[i]);
} }
@ -66,12 +66,12 @@ void VCS_SOLVE::vcs_nondim_TP()
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice\n", m_totalMoleScale); plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice\n", m_totalMoleScale);
} }
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale); m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale);
} }
} }
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale); m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale);
} }
@ -92,7 +92,7 @@ void VCS_SOLVE::vcs_redim_TP()
if (m_unitsState != VCS_DIMENSIONAL_G) { if (m_unitsState != VCS_DIMENSIONAL_G) {
m_unitsState = VCS_DIMENSIONAL_G; m_unitsState = VCS_DIMENSIONAL_G;
double tf = m_temperature * GasConstant; double tf = m_temperature * GasConstant;
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
// Modify the standard state and total chemical potential data, // Modify the standard state and total chemical potential data,
// FF(I), to make it have units, i.e. mu = RT * mu_star // FF(I), to make it have units, i.e. mu = RT * mu_star
@ -107,12 +107,12 @@ void VCS_SOLVE::vcs_redim_TP()
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g\n", m_totalMoleScale); plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g\n", m_totalMoleScale);
} }
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= m_totalMoleScale; m_molNumSpecies_old[i] *= m_totalMoleScale;
} }
} }
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
m_elemAbundancesGoal[i] *= m_totalMoleScale; m_elemAbundancesGoal[i] *= m_totalMoleScale;
} }

60
src/equil/vcs_prep.cpp
View file

@ -15,7 +15,7 @@ namespace Cantera
void VCS_SOLVE::vcs_SSPhase() void VCS_SOLVE::vcs_SSPhase()
{ {
vector_int numPhSpecies(m_numPhases, 0); vector_int numPhSpecies(m_numPhases, 0);
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
numPhSpecies[m_phaseID[kspec]]++; numPhSpecies[m_phaseID[kspec]]++;
} }
@ -36,7 +36,7 @@ void VCS_SOLVE::vcs_SSPhase()
// Fill in some useful arrays here that have to do with the static // Fill in some useful arrays here that have to do with the static
// information concerning the phase ID of species. SSPhase = Boolean // information concerning the phase ID of species. SSPhase = Boolean
// indicating whether a species is in a single species phase or not. // indicating whether a species is in a single species phase or not.
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
size_t iph = m_phaseID[kspec]; size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph]; vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (Vphase->m_singleSpecies) { if (Vphase->m_singleSpecies) {
@ -47,7 +47,7 @@ void VCS_SOLVE::vcs_SSPhase()
} }
} }
int VCS_SOLVE::vcs_prep_oneTime(int printLvl) int VCS_SOLVE::vcs_prep(int printLvl)
{ {
int retn = VCS_SUCCESS; int retn = VCS_SUCCESS;
m_debug_print_lvl = printLvl; m_debug_print_lvl = printLvl;
@ -58,18 +58,18 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
// Set an initial estimate for the number of noncomponent species equal to // Set an initial estimate for the number of noncomponent species equal to
// nspecies - nelements. This may be changed below // nspecies - nelements. This may be changed below
if (m_numElemConstraints > m_numSpeciesTot) { if (m_nelem > m_nsp) {
m_numRxnTot = 0; m_numRxnTot = 0;
} else { } else {
m_numRxnTot = m_numSpeciesTot - m_numElemConstraints; m_numRxnTot = m_nsp - m_nelem;
} }
m_numRxnRdc = m_numRxnTot; m_numRxnRdc = m_numRxnTot;
m_numSpeciesRdc = m_numSpeciesTot; m_numSpeciesRdc = m_nsp;
for (size_t i = 0; i < m_numRxnRdc; ++i) { for (size_t i = 0; i < m_numRxnRdc; ++i) {
m_indexRxnToSpecies[i] = m_numElemConstraints + i; m_indexRxnToSpecies[i] = m_nelem + i;
} }
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
size_t pID = m_phaseID[kspec]; size_t pID = m_phaseID[kspec];
size_t spPhIndex = m_speciesLocalPhaseIndex[kspec]; size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
vcs_VolPhase* vPhase = m_VolPhaseList[pID]; vcs_VolPhase* vPhase = m_VolPhaseList[pID];
@ -104,7 +104,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
bool modifiedSoln = false; bool modifiedSoln = false;
if (m_doEstimateEquil < 0) { if (m_doEstimateEquil < 0) {
double sum = 0.0; double sum = 0.0;
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
sum += fabs(m_molNumSpecies_old[kspec]); sum += fabs(m_molNumSpecies_old[kspec]);
} }
@ -113,7 +113,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
modifiedSoln = true; modifiedSoln = true;
double pres = (m_pressurePA <= 0.0) ? 1.01325E5 : m_pressurePA; double pres = (m_pressurePA <= 0.0) ? 1.01325E5 : m_pressurePA;
retn = vcs_evalSS_TP(0, 0, m_temperature, pres); retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec]; m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
} else { } else {
@ -126,15 +126,15 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
// NC = number of components is in the vcs.h common block. This call to // NC = number of components is in the vcs.h common block. This call to
// BASOPT doesn't calculate the stoichiometric reaction matrix. // BASOPT doesn't calculate the stoichiometric reaction matrix.
vector_fp awSpace(m_numSpeciesTot + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0); vector_fp awSpace(m_nsp + (m_nelem + 2)*(m_nelem), 0.0);
double* aw = &awSpace[0]; double* aw = &awSpace[0];
if (aw == NULL) { if (aw == NULL) {
plogf("vcs_prep_oneTime: failed to get memory: global bailout\n"); plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
return VCS_NOMEMORY; return VCS_NOMEMORY;
} }
double* sa = aw + m_numSpeciesTot; double* sa = aw + m_nsp;
double* sm = sa + m_numElemConstraints; double* sm = sa + m_nelem;
double* ss = sm + (m_numElemConstraints)*(m_numElemConstraints); double* ss = sm + m_nelem * m_nelem;
bool conv; bool conv;
retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv); retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv);
if (retn != VCS_SUCCESS) { if (retn != VCS_SUCCESS) {
@ -145,8 +145,8 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
return retn; return retn;
} }
if (m_numSpeciesTot >= m_numComponents) { if (m_nsp >= m_numComponents) {
m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents; m_numRxnTot = m_numRxnRdc = m_nsp - m_numComponents;
for (size_t i = 0; i < m_numRxnRdc; ++i) { for (size_t i = 0; i < m_numRxnRdc; ++i) {
m_indexRxnToSpecies[i] = m_numComponents + i; m_indexRxnToSpecies[i] = m_numComponents + i;
} }
@ -155,11 +155,11 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
} }
// The elements might need to be rearranged. // The elements might need to be rearranged.
awSpace.resize(m_numElemConstraints + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0); awSpace.resize(m_nelem + (m_nelem + 2)*m_nelem, 0.0);
aw = &awSpace[0]; aw = &awSpace[0];
sa = aw + m_numElemConstraints; sa = aw + m_nelem;
sm = sa + m_numElemConstraints; sm = sa + m_nelem;
ss = sm + (m_numElemConstraints)*(m_numElemConstraints); ss = sm + m_nelem * m_nelem;
retn = vcs_elem_rearrange(aw, sa, sm, ss); retn = vcs_elem_rearrange(aw, sa, sm, ss);
if (retn != VCS_SUCCESS) { if (retn != VCS_SUCCESS) {
plogf("vcs_prep_oneTime:"); plogf("vcs_prep_oneTime:");
@ -172,15 +172,11 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
// If we mucked up the solution unknowns because they were all // If we mucked up the solution unknowns because they were all
// zero to start with, set them back to zero here // zero to start with, set them back to zero here
if (modifiedSoln) { if (modifiedSoln) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
m_molNumSpecies_old[kspec] = 0.0; m_molNumSpecies_old[kspec] = 0.0;
} }
} }
return VCS_SUCCESS;
}
int VCS_SOLVE::vcs_prep()
{
// Initialize various arrays in the data to zero // Initialize various arrays in the data to zero
m_feSpecies_old.assign(m_feSpecies_old.size(), 0.0); m_feSpecies_old.assign(m_feSpecies_old.size(), 0.0);
m_feSpecies_new.assign(m_feSpecies_new.size(), 0.0); m_feSpecies_new.assign(m_feSpecies_new.size(), 0.0);
@ -192,20 +188,18 @@ int VCS_SOLVE::vcs_prep()
// Calculate the total number of moles in all phases. // Calculate the total number of moles in all phases.
vcs_tmoles(); vcs_tmoles();
return VCS_SUCCESS;
}
bool VCS_SOLVE::vcs_wellPosed() // Check to see if the current problem is well posed.
{
double sum = 0.0; double sum = 0.0;
for (size_t e = 0; e < ne; e++) { for (size_t e = 0; e < m_nelem; e++) {
sum += m_mix->elementMoles(e); sum += m_mix->elementMoles(e);
} }
if (sum < 1.0E-20) { if (sum < 1.0E-20) {
plogf("vcs_wellPosed: Element abundance is close to zero\n"); // Check to see if the current problem is well posed.
return false; plogf("vcs has determined the problem is not well posed: Bailing\n");
return VCS_PUB_BAD;
} }
return true; return VCS_SUCCESS;
} }
} }

46
src/equil/vcs_prob.cpp
View file

@ -25,8 +25,8 @@ namespace Cantera
void VCS_SOLVE::set_gai() void VCS_SOLVE::set_gai()
{ {
gai.assign(gai.size(), 0.0); gai.assign(gai.size(), 0.0);
for (size_t j = 0; j < ne; j++) { for (size_t j = 0; j < m_nelem; j++) {
for (size_t kspec = 0; kspec < nspecies; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
gai[j] += FormulaMatrix(kspec,j) * w[kspec]; gai[j] += FormulaMatrix(kspec,j) * w[kspec];
} }
@ -46,24 +46,20 @@ void VCS_SOLVE::prob_report(int print_lvl)
writeline('=', 31); writeline('=', 31);
writeline('=', 80); writeline('=', 80);
plogf("\n"); plogf("\n");
if (prob_type == 0) { plogf("\tSolve a constant T, P problem:\n");
plogf("\tSolve a constant T, P problem:\n"); plogf("\t\tT = %g K\n", m_temperature);
plogf("\t\tT = %g K\n", T); double pres_atm = m_pressurePA / 1.01325E5;
double pres_atm = PresPA / 1.01325E5;
plogf("\t\tPres = %g atm\n", pres_atm); plogf("\t\tPres = %g atm\n", pres_atm);
} else {
throw CanteraError("VCS_PROB::prob_report", "Unknown problem type");
}
plogf("\n"); plogf("\n");
plogf(" Phase IDs of species\n"); plogf(" Phase IDs of species\n");
plogf(" species phaseID phaseName "); plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_Moles Species_Type\n"); plogf(" Initial_Estimated_Moles Species_Type\n");
for (size_t i = 0; i < nspecies; i++) { for (size_t i = 0; i < m_nsp; i++) {
vcs_VolPhase* Vphase = VPhaseList[PhaseID[i]]; vcs_VolPhase* Vphase = VPhaseList[PhaseID[i]];
plogf("%16s %5d %16s", m_mix->speciesName(i), PhaseID[i], plogf("%16s %5d %16s", m_mix->speciesName(i), PhaseID[i],
Vphase->PhaseName); Vphase->PhaseName);
if (iest >= 0) { if (m_doEstimateEquil >= 0) {
plogf(" %-10.5g", w[i]); plogf(" %-10.5g", w[i]);
} else { } else {
plogf(" N/A"); plogf(" N/A");
@ -85,13 +81,13 @@ void VCS_SOLVE::prob_report(int print_lvl)
" EqnState NumSpec"); " EqnState NumSpec");
plogf(" TMolesInert TKmoles\n"); plogf(" TMolesInert TKmoles\n");
for (size_t iphase = 0; iphase < NPhase; iphase++) { for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = VPhaseList[iphase]; vcs_VolPhase* Vphase = VPhaseList[iphase];
plogf("%16s %5d %5d %8d ", Vphase->PhaseName, plogf("%16s %5d %5d %8d ", Vphase->PhaseName,
Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase); Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase);
plogf("%16s %8d %16e ", Vphase->eos_name(), plogf("%16s %8d %16e ", Vphase->eos_name(),
Vphase->nSpecies(), Vphase->totalMolesInert()); Vphase->nSpecies(), Vphase->totalMolesInert());
if (iest >= 0) { if (m_doEstimateEquil >= 0) {
plogf("%16e\n", Vphase->totalMoles()); plogf("%16e\n", Vphase->totalMoles());
} else { } else {
plogf(" N/A\n"); plogf(" N/A\n");
@ -100,7 +96,7 @@ void VCS_SOLVE::prob_report(int print_lvl)
plogf("\nElemental Abundances: "); plogf("\nElemental Abundances: ");
plogf(" Target_kmol ElemType ElActive\n"); plogf(" Target_kmol ElemType ElActive\n");
for (size_t i = 0; i < ne; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
writeline(' ', 26, false); writeline(' ', 26, false);
plogf("%-2.2s", m_elementName[i]); plogf("%-2.2s", m_elementName[i]);
plogf("%20.12E ", gai[i]); plogf("%20.12E ", gai[i]);
@ -110,9 +106,9 @@ void VCS_SOLVE::prob_report(int print_lvl)
plogf("\nChemical Potentials: (J/kmol)\n"); plogf("\nChemical Potentials: (J/kmol)\n");
plogf(" Species (phase) " plogf(" Species (phase) "
" SS0ChemPot StarChemPot\n"); " SS0ChemPot StarChemPot\n");
for (size_t iphase = 0; iphase < NPhase; iphase++) { for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = VPhaseList[iphase]; vcs_VolPhase* Vphase = VPhaseList[iphase];
Vphase->setState_TP(T, PresPA); Vphase->setState_TP(m_temperature, m_pressurePA);
for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) { for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) {
size_t kglob = Vphase->spGlobalIndexVCS(kindex); size_t kglob = Vphase->spGlobalIndexVCS(kindex);
plogf("%16s ", m_mix->speciesName(kglob)); plogf("%16s ", m_mix->speciesName(kglob));
@ -145,7 +141,7 @@ void VCS_SOLVE::addPhaseElements(vcs_VolPhase* volPhase)
// Search for matches with the existing elements. If found, then fill in // Search for matches with the existing elements. If found, then fill in
// the entry in the global mapping array. // the entry in the global mapping array.
for (size_t e = 0; e < ne; e++) { for (size_t e = 0; e < m_nelem; e++) {
std::string en = m_elementName[e]; std::string en = m_elementName[e];
if (!strcmp(enVP.c_str(), en.c_str())) { if (!strcmp(enVP.c_str(), en.c_str())) {
volPhase->setElemGlobalIndex(eVP, e); volPhase->setElemGlobalIndex(eVP, e);
@ -167,14 +163,13 @@ size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive)
throw CanteraError("VCS_SOLVE::addElement", throw CanteraError("VCS_SOLVE::addElement",
"error: element must have a name"); "error: element must have a name");
} }
ne++; m_nelem++;
m_numElemConstraints++;
m_numComponents++; m_numComponents++;
gai.push_back(0.0); gai.push_back(0.0);
FormulaMatrix.resize(NSPECIES0, ne, 0.0); FormulaMatrix.resize(m_nsp, m_nelem, 0.0);
m_formulaMatrix.resize(NSPECIES0, ne); m_formulaMatrix.resize(m_nsp, m_nelem);
m_stoichCoeffRxnMatrix.resize(ne, NSPECIES0, 0.0); m_stoichCoeffRxnMatrix.resize(m_nelem, m_nsp, 0.0);
m_elType.push_back(elType); m_elType.push_back(elType);
ElActive.push_back(elactive); ElActive.push_back(elactive);
m_elementActive.push_back(elactive); m_elementActive.push_back(elactive);
@ -182,13 +177,12 @@ size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive)
m_elemAbundancesGoal.push_back(0.0); m_elemAbundancesGoal.push_back(0.0);
m_elementMapIndex.push_back(0); m_elementMapIndex.push_back(0);
m_elementName.push_back(elNameNew); m_elementName.push_back(elNameNew);
NE0 = ne; return m_nelem - 1;
return ne - 1;
} }
size_t VCS_SOLVE::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT) size_t VCS_SOLVE::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
{ {
if (kT > nspecies) { if (kT > m_nsp) {
// Need to expand the number of species here // Need to expand the number of species here
throw CanteraError("VCS_PROB::addOnePhaseSpecies", "Shouldn't be here"); throw CanteraError("VCS_PROB::addOnePhaseSpecies", "Shouldn't be here");
} }

39
src/equil/vcs_report.cpp
View file

@ -12,13 +12,12 @@ namespace Cantera
int VCS_SOLVE::vcs_report(int iconv) int VCS_SOLVE::vcs_report(int iconv)
{ {
bool printActualMoles = true, inertYes = false; bool printActualMoles = true, inertYes = false;
size_t nspecies = m_numSpeciesTot;
char originalUnitsState = m_unitsState; char originalUnitsState = m_unitsState;
std::vector<size_t> sortindex(nspecies,0); std::vector<size_t> sortindex(m_nsp, 0);
vector_fp xy(nspecies,0.0); vector_fp xy(m_nsp, 0.0);
// SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES // SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES
for (size_t i = 0; i < nspecies; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
sortindex[i] = i; sortindex[i] = i;
xy[i] = m_molNumSpecies_old[i]; xy[i] = m_molNumSpecies_old[i];
} }
@ -45,7 +44,7 @@ int VCS_SOLVE::vcs_report(int iconv)
} }
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp);
// PRINT OUT RESULTS // PRINT OUT RESULTS
plogf("\n\n\n\n"); plogf("\n\n\n\n");
@ -117,9 +116,9 @@ int VCS_SOLVE::vcs_report(int iconv)
TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0); TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
} }
} }
if (m_numSpeciesRdc != nspecies) { if (m_numSpeciesRdc != m_nsp) {
plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n"); plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
for (size_t kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) { for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
plogf(" %-12.12s", m_speciesName[kspec]); plogf(" %-12.12s", m_speciesName[kspec]);
// Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
plogf(" %14.7E %14.7E %12.4E", plogf(" %14.7E %14.7E %12.4E",
@ -171,30 +170,30 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\n"); plogf("\n");
// TABLE OF PHASE INFORMATION // TABLE OF PHASE INFORMATION
vector_fp gaPhase(m_numElemConstraints, 0.0); vector_fp gaPhase(m_nelem, 0.0);
vector_fp gaTPhase(m_numElemConstraints, 0.0); vector_fp gaTPhase(m_nelem, 0.0);
double totalMoles = 0.0; double totalMoles = 0.0;
double gibbsPhase = 0.0; double gibbsPhase = 0.0;
double gibbsTotal = 0.0; double gibbsTotal = 0.0;
plogf("\n\n"); plogf("\n\n");
plogf("\n"); plogf("\n");
writeline('-', m_numElemConstraints*10 + 58); writeline('-', m_nelem*10 + 58);
plogf(" | ElementID |"); plogf(" | ElementID |");
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %3d", j); plogf(" %3d", j);
} }
plogf(" | |\n"); plogf(" | |\n");
plogf(" | Element |"); plogf(" | Element |");
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %10.10s", m_elementName[j]); plogf(" %10.10s", m_elementName[j]);
} }
plogf(" | |\n"); plogf(" | |\n");
plogf(" PhaseName |KMolTarget |"); plogf(" PhaseName |KMolTarget |");
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %10.3g", m_elemAbundancesGoal[j]); plogf(" %10.3g", m_elemAbundancesGoal[j]);
} }
plogf(" | Gibbs Total |\n"); plogf(" | Gibbs Total |\n");
writeline('-', m_numElemConstraints*10 + 58); writeline('-', m_nelem*10 + 58);
for (size_t iphase = 0; iphase < m_numPhases; iphase++) { for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
plogf(" %3d ", iphase); plogf(" %3d ", iphase);
vcs_VolPhase* VPhase = m_VolPhaseList[iphase]; vcs_VolPhase* VPhase = m_VolPhaseList[iphase];
@ -206,7 +205,7 @@ int VCS_SOLVE::vcs_report(int iconv)
throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem"); throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
} }
vcs_elabPhase(iphase, &gaPhase[0]); vcs_elabPhase(iphase, &gaPhase[0]);
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %10.3g", gaPhase[j]); plogf(" %10.3g", gaPhase[j]);
gaTPhase[j] += gaPhase[j]; gaTPhase[j] += gaPhase[j];
} }
@ -215,14 +214,14 @@ int VCS_SOLVE::vcs_report(int iconv)
gibbsTotal += gibbsPhase; gibbsTotal += gibbsPhase;
plogf(" | %18.11E |\n", gibbsPhase); plogf(" | %18.11E |\n", gibbsPhase);
} }
writeline('-', m_numElemConstraints*10 + 58); writeline('-', m_nelem*10 + 58);
plogf(" TOTAL |%10.3e |", totalMoles); plogf(" TOTAL |%10.3e |", totalMoles);
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %10.3g", gaTPhase[j]); plogf(" %10.3g", gaTPhase[j]);
} }
plogf(" | %18.11E |\n", gibbsTotal); plogf(" | %18.11E |\n", gibbsTotal);
writeline('-', m_numElemConstraints*10 + 58); writeline('-', m_nelem*10 + 58);
plogf("\n"); plogf("\n");
// GLOBAL SATISFACTION INFORMATION // GLOBAL SATISFACTION INFORMATION
@ -238,7 +237,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\nElemental Abundances (kmol): "); plogf("\nElemental Abundances (kmol): ");
plogf(" Actual Target Type ElActive\n"); plogf(" Actual Target Type ElActive\n");
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
writeline(' ', 26, false); writeline(' ', 26, false);
plogf("%-2.2s", m_elementName[i]); plogf("%-2.2s", m_elementName[i]);
plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale); plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale);
@ -255,7 +254,7 @@ int VCS_SOLVE::vcs_report(int iconv)
" ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)"); " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
plogf("| (MolNum ChemPot)|"); plogf("| (MolNum ChemPot)|");
writeline('-', 147, true, true); writeline('-', 147, true, true);
for (size_t i = 0; i < nspecies; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
size_t j = sortindex[i]; size_t j = sortindex[i];
size_t pid = m_phaseID[j]; size_t pid = m_phaseID[j];
plogf(" %-12.12s", m_speciesName[j]); plogf(" %-12.12s", m_speciesName[j]);

2
src/equil/vcs_rxnadj.cpp
View file

@ -242,7 +242,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
} }
// Delete the single species phase // Delete the single species phase
for (size_t j = 0; j < m_numSpeciesTot; j++) { for (size_t j = 0; j < m_nsp; j++) {
m_deltaMolNumSpecies[j] = 0.0; m_deltaMolNumSpecies[j] = 0.0;
} }
m_deltaMolNumSpecies[kspec] = dss; m_deltaMolNumSpecies[kspec] = dss;

14
src/equil/vcs_setMolesLinProg.cpp
View file

@ -40,13 +40,13 @@ int VCS_SOLVE::vcs_setMolesLinProg()
int iter = 0; int iter = 0;
bool abundancesOK = true; bool abundancesOK = true;
bool usedZeroedSpecies; bool usedZeroedSpecies;
vector_fp sm(m_numElemConstraints*m_numElemConstraints, 0.0); vector_fp sm(m_nelem * m_nelem, 0.0);
vector_fp ss(m_numElemConstraints, 0.0); vector_fp ss(m_nelem, 0.0);
vector_fp sa(m_numElemConstraints, 0.0); vector_fp sa(m_nelem, 0.0);
vector_fp wx(m_numElemConstraints, 0.0); vector_fp wx(m_nelem, 0.0);
vector_fp aw(m_numSpeciesTot, 0.0); vector_fp aw(m_nsp, 0.0);
for (size_t ik = 0; ik < m_numSpeciesTot; ik++) { for (size_t ik = 0; ik < m_nsp; ik++) {
if (m_speciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]); m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]);
} }
@ -98,7 +98,7 @@ int VCS_SOLVE::vcs_setMolesLinProg()
double dg_rt = m_SSfeSpecies[ik]; double dg_rt = m_SSfeSpecies[ik];
dxi_min = 1.0e10; dxi_min = 1.0e10;
const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
for (size_t jcomp = 0; jcomp < m_numElemConstraints; jcomp++) { for (size_t jcomp = 0; jcomp < m_nelem; jcomp++) {
dg_rt += m_SSfeSpecies[jcomp] * sc_irxn[jcomp]; dg_rt += m_SSfeSpecies[jcomp] * sc_irxn[jcomp];
} }
// fwd or rev direction. // fwd or rev direction.

1162
src/equil/vcs_solve.cpp
View file

File diff suppressed because it is too large Load diff

142
src/equil/vcs_solve_TP.cpp
View file

@ -66,11 +66,11 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
clockWC ticktock; clockWC ticktock;
// temporary space for usage in this routine and in subroutines // temporary space for usage in this routine and in subroutines
m_sm.assign(m_numElemConstraints*m_numElemConstraints, 0.0); m_sm.assign(m_nelem * m_nelem, 0.0);
m_ss.assign(m_numElemConstraints, 0.0); m_ss.assign(m_nelem, 0.0);
m_sa.assign(m_numElemConstraints, 0.0); m_sa.assign(m_nelem, 0.0);
m_aw.assign(m_numSpeciesTot, 0.0); m_aw.assign(m_nsp, 0.0);
m_wx.assign(m_numElemConstraints, 0.0); m_wx.assign(m_nelem, 0.0);
int solveFail = false; int solveFail = false;
@ -80,8 +80,8 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
// Printout the initial conditions for problem // Printout the initial conditions for problem
if (print_lvl != 0) { if (print_lvl != 0) {
plogf("VCS CALCULATION METHOD\n\n "); plogf("VCS CALCULATION METHOD\n\n ");
plogf("%s\n", m_title); plogf("MultiPhase Object\n");
plogf("\n\n%5d SPECIES\n%5d ELEMENTS\n", m_numSpeciesTot, m_numElemConstraints); plogf("\n\n%5d SPECIES\n%5d ELEMENTS\n", m_nsp, m_nelem);
plogf("%5d COMPONENTS\n", m_numComponents); plogf("%5d COMPONENTS\n", m_numComponents);
plogf("%5d PHASES\n", m_numPhases); plogf("%5d PHASES\n", m_numPhases);
plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa "); plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa ");
@ -95,7 +95,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
} }
plogf("\n ELEMENTAL ABUNDANCES CORRECT"); plogf("\n ELEMENTAL ABUNDANCES CORRECT");
plogf(" FROM ESTIMATE Type\n\n"); plogf(" FROM ESTIMATE Type\n\n");
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
writeline(' ', 26, false); writeline(' ', 26, false);
plogf("%-2.2s", m_elementName[i]); plogf("%-2.2s", m_elementName[i]);
plogf("%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i], plogf("%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i],
@ -114,18 +114,18 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
writeline(' ', 41, false); writeline(' ', 41, false);
plogf(" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n"); plogf(" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n");
writeline(' ', 20, false); writeline(' ', 20, false);
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
plogf("%-4.4s ", m_elementName[i]); plogf("%-4.4s ", m_elementName[i]);
} }
plogf(" PhaseID\n"); plogf(" PhaseID\n");
double RT = GasConstant * m_temperature; double RT = GasConstant * m_temperature;
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
plogf(" %-18.18s", m_speciesName[i]); plogf(" %-18.18s", m_speciesName[i]);
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
plogf("% -7.3g ", m_formulaMatrix(i,j)); plogf("% -7.3g ", m_formulaMatrix(i,j));
} }
plogf(" %3d ", m_phaseID[i]); plogf(" %3d ", m_phaseID[i]);
writeline(' ', std::max(55-int(m_numElemConstraints)*8, 0), false); writeline(' ', std::max(55-int(m_nelem)*8, 0), false);
plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]); plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]);
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" Mol_Num\n"); plogf(" Mol_Num\n");
@ -137,7 +137,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
} }
} }
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
if (m_molNumSpecies_old[i] < 0.0) { if (m_molNumSpecies_old[i] < 0.0) {
plogf("On Input species %-12s has a negative MF, setting it small\n", plogf("On Input species %-12s has a negative MF, setting it small\n",
m_speciesName[i]); m_speciesName[i]);
@ -250,7 +250,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
// Evaluate the final mole fractions storing them in wt[] // Evaluate the final mole fractions storing them in wt[]
m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0); m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0);
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
if (m_SSPhase[kspec]) { if (m_SSPhase[kspec]) {
m_molNumSpecies_new[kspec] = 1.0; m_molNumSpecies_new[kspec] = 1.0;
} else { } else {
@ -424,13 +424,13 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers // m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers
if (iphaseDelete != npos) { if (iphaseDelete != npos) {
debuglog(" --- Main Loop Treatment -> Circumvented due to Phase Deletion\n", m_debug_print_lvl >= 2); debuglog(" --- Main Loop Treatment -> Circumvented due to Phase Deletion\n", m_debug_print_lvl >= 2);
for (size_t k = 0; k < m_numSpeciesTot; k++) { for (size_t k = 0; k < m_nsp; k++) {
m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k]; m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k];
size_t iph = m_phaseID[k]; size_t iph = m_phaseID[k];
m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k]; m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k];
} }
if (kspec >= m_numComponents) { if (kspec >= m_numComponents) {
if (m_molNumSpecies_new[m_numSpeciesTot] != 0.0) { if (m_molNumSpecies_new[m_nsp] != 0.0) {
throw CanteraError("VCS_SOLVE::solve_tp_inner", throw CanteraError("VCS_SOLVE::solve_tp_inner",
"we shouldn't be here!"); "we shouldn't be here!");
} }
@ -451,7 +451,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// Calculate the new chemical potentials using the tentative solution // Calculate the new chemical potentials using the tentative solution
// values. We only calculate a subset of these, because we have only // values. We only calculate a subset of these, because we have only
// updated a subset of the W(). // updated a subset of the W().
vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_nsp);
// Evaluate DeltaG for all components if ITI=0, and for major components // Evaluate DeltaG for all components if ITI=0, and for major components
// only if ITI NE 0 // only if ITI NE 0
@ -515,7 +515,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
} }
} }
} else { } else {
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
int elType = m_elType[j]; int elType = m_elType[j];
if (elType == VCS_ELEM_TYPE_ABSPOS) { if (elType == VCS_ELEM_TYPE_ABSPOS) {
double atomComp = m_formulaMatrix(kspec,j); double atomComp = m_formulaMatrix(kspec,j);
@ -734,7 +734,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
++m_numRxnMinorZeroed; ++m_numRxnMinorZeroed;
allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
for (size_t kk = 0; kk < m_numSpeciesTot; kk++) { for (size_t kk = 0; kk < m_nsp; kk++) {
m_deltaMolNumSpecies[kk] = 0.0; m_deltaMolNumSpecies[kk] = 0.0;
m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk]; m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk];
} }
@ -837,7 +837,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
plogf(" --- Reduction in step size due to component %s going negative = %11.3E\n", plogf(" --- Reduction in step size due to component %s going negative = %11.3E\n",
m_speciesName[ll], par); m_speciesName[ll], par);
} }
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
m_deltaMolNumSpecies[i] *= par; m_deltaMolNumSpecies[i] *= par;
} }
for (size_t iph = 0; iph < m_numPhases; iph++) { for (size_t iph = 0; iph < m_numPhases; iph++) {
@ -847,14 +847,14 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
par = 1.0; par = 1.0;
} }
checkDelta1(&m_deltaMolNumSpecies[0], checkDelta1(&m_deltaMolNumSpecies[0],
&m_deltaPhaseMoles[0], m_numSpeciesTot); &m_deltaPhaseMoles[0], m_nsp);
// Now adjust the wt[kspec]'s so that the reflect the decrease in the // Now adjust the wt[kspec]'s so that the reflect the decrease in the
// overall length of m_deltaMolNumSpecies[kspec] just calculated. At the // overall length of m_deltaMolNumSpecies[kspec] just calculated. At the
// end of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and // end of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and
// tPhMoles1 should all be consistent with a new estimate of the state // tPhMoles1 should all be consistent with a new estimate of the state
// of the system. // of the system.
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec] if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec]
!= VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
@ -876,7 +876,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// Calculate the new chemical potentials using the tentative solution // Calculate the new chemical potentials using the tentative solution
// values. We only calculate a subset of these, because we have only // values. We only calculate a subset of these, because we have only
// updated a subset of the W(). // updated a subset of the W().
vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_nsp);
// Evaluate DeltaG for all components if ITI=0, and for major components // Evaluate DeltaG for all components if ITI=0, and for major components
// only if ITI NE 0 // only if ITI NE 0
@ -962,7 +962,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i], m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i],
m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]); m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
} }
for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
size_t l1 = kspec - m_numComponents; size_t l1 = kspec - m_numComponents;
plogf(" --- %-12.12s", m_speciesName[kspec]); plogf(" --- %-12.12s", m_speciesName[kspec]);
plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
@ -1280,7 +1280,7 @@ void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
// If we have deleted a species then we need to recheck the the deleted // If we have deleted a species then we need to recheck the the deleted
// species, before exiting // species, before exiting
if (m_numSpeciesRdc != m_numSpeciesTot) { if (m_numSpeciesRdc != m_nsp) {
stage = RECHECK_DELETED; stage = RECHECK_DELETED;
return; return;
} }
@ -1600,7 +1600,7 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec)
if (! m_SSPhase[kspec]) { if (! m_SSPhase[kspec]) {
if (Vphase->exists() == VCS_PHASE_EXIST_NO) { if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
Vphase->setExistence(VCS_PHASE_EXIST_YES); Vphase->setExistence(VCS_PHASE_EXIST_YES);
for (size_t k = 0; k < m_numSpeciesTot; k++) { for (size_t k = 0; k < m_nsp; k++) {
if (m_phaseID[k] == iph && m_speciesStatus[k] != VCS_SPECIES_DELETED) { if (m_phaseID[k] == iph && m_speciesStatus[k] != VCS_SPECIES_DELETED) {
m_speciesStatus[k] = VCS_SPECIES_MINOR; m_speciesStatus[k] = VCS_SPECIES_MINOR;
} }
@ -1723,7 +1723,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
// reinstate the "major ones" in the future. Note, species in phases with // reinstate the "major ones" in the future. Note, species in phases with
// zero mole numbers are still considered active. Whether the phase pops // zero mole numbers are still considered active. Whether the phase pops
// back into existence or not is checked as part of the main iteration loop. // back into existence or not is checked as part of the main iteration loop.
for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
if (m_phaseID[kspec] == iph) { if (m_phaseID[kspec] == iph) {
m_molNumSpecies_old[kspec] = 0.0; m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0; m_molNumSpecies_new[kspec] = 0.0;
@ -1760,14 +1760,14 @@ int VCS_SOLVE::vcs_recheck_deleted()
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
plogf(" --- Start rechecking deleted species in multispec phases\n"); plogf(" --- Start rechecking deleted species in multispec phases\n");
} }
if (m_numSpeciesRdc == m_numSpeciesTot) { if (m_numSpeciesRdc == m_nsp) {
return 0; return 0;
} }
// Use the standard chemical potentials for the chemical potentials of // Use the standard chemical potentials for the chemical potentials of
// deleted species. Then, calculate Delta G for for formation reactions. // deleted species. Then, calculate Delta G for for formation reactions.
// Note: fe[] here includes everything except for the ln(x[i]) term // Note: fe[] here includes everything except for the ln(x[i]) term
for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
size_t iph = m_phaseID[kspec]; size_t iph = m_phaseID[kspec];
m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
- m_lnMnaughtSpecies[kspec] - m_lnMnaughtSpecies[kspec]
@ -1830,7 +1830,7 @@ int VCS_SOLVE::vcs_recheck_deleted()
size_t VCS_SOLVE::vcs_add_all_deleted() size_t VCS_SOLVE::vcs_add_all_deleted()
{ {
if (m_numSpeciesRdc == m_numSpeciesTot) { if (m_numSpeciesRdc == m_nsp) {
return 0; return 0;
} }
@ -1841,7 +1841,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
// of the routine. // of the routine.
m_molNumSpecies_new = m_molNumSpecies_old; m_molNumSpecies_new = m_molNumSpecies_old;
for (int cits = 0; cits < 3; cits++) { for (int cits = 0; cits < 3; cits++) {
for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
size_t iph = m_phaseID[kspec]; size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph]; vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (m_molNumSpecies_new[kspec] == 0.0) { if (m_molNumSpecies_new[kspec] == 0.0) {
@ -1898,7 +1898,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
} }
} }
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp);
vcs_deltag(0, true, VCS_STATECALC_OLD); vcs_deltag(0, true, VCS_STATECALC_OLD);
size_t retn = 0; size_t retn = 0;
@ -2050,18 +2050,18 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
plogf("\n"); plogf("\n");
plogf(" --- Formula Matrix used in BASOPT calculation\n"); plogf(" --- Formula Matrix used in BASOPT calculation\n");
plogf(" --- Active | "); plogf(" --- Active | ");
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %1d ", m_elementActive[j]); plogf(" %1d ", m_elementActive[j]);
} }
plogf("\n"); plogf("\n");
plogf(" --- Species | "); plogf(" --- Species | ");
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
writelog(" {:>8.8s}", m_elementName[j]); writelog(" {:>8.8s}", m_elementName[j]);
} }
plogf("\n"); plogf("\n");
for (k = 0; k < m_numSpeciesTot; k++) { for (k = 0; k < m_nsp; k++) {
writelog(" --- {:>11.11s} | ", m_speciesName[k]); writelog(" --- {:>11.11s} | ", m_speciesName[k]);
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
plogf(" %8.2g", m_formulaMatrix(k,j)); plogf(" %8.2g", m_formulaMatrix(k,j));
} }
plogf("\n"); plogf("\n");
@ -2073,17 +2073,17 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
// Calculate the maximum value of the number of components possible. It's // Calculate the maximum value of the number of components possible. It's
// equal to the minimum of the number of elements and the number of total // equal to the minimum of the number of elements and the number of total
// species. // species.
size_t ncTrial = std::min(m_numElemConstraints, m_numSpeciesTot); size_t ncTrial = std::min(m_nelem, m_nsp);
m_numComponents = ncTrial; m_numComponents = ncTrial;
*usedZeroedSpecies = false; *usedZeroedSpecies = false;
vector_int ipiv(ncTrial); vector_int ipiv(ncTrial);
// Use a temporary work array for the mole numbers, aw[] // Use a temporary work array for the mole numbers, aw[]
std::copy(m_molNumSpecies_old.begin(), std::copy(m_molNumSpecies_old.begin(),
m_molNumSpecies_old.begin() + m_numSpeciesTot, aw); m_molNumSpecies_old.begin() + m_nsp, aw);
// Take out the Voltage unknowns from consideration // Take out the Voltage unknowns from consideration
for (k = 0; k < m_numSpeciesTot; k++) { for (k = 0; k < m_nsp; k++) {
if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
aw[k] = test; aw[k] = test;
} }
@ -2101,7 +2101,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
// the largest remaining species. Return its identity in K. The // the largest remaining species. Return its identity in K. The
// first search criteria is always the largest positive magnitude of // first search criteria is always the largest positive magnitude of
// the mole number. // the mole number.
k = vcs_basisOptMax(aw, jr, m_numSpeciesTot); k = vcs_basisOptMax(aw, jr, m_nsp);
// The fun really starts when you have run out of species that have // The fun really starts when you have run out of species that have
// a significant concentration. It becomes extremely important to // a significant concentration. It becomes extremely important to
@ -2147,11 +2147,11 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
size_t kfound = npos; size_t kfound = npos;
int minNonZeroes = 100000; int minNonZeroes = 100000;
int nonZeroesKspec = 0; int nonZeroesKspec = 0;
for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = ncTrial; kspec < m_nsp; kspec++) {
if (aw[kspec] >= 0.0 && m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (aw[kspec] >= 0.0 && m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
maxConcPossKspec = 1.0E10; maxConcPossKspec = 1.0E10;
nonZeroesKspec = 0; nonZeroesKspec = 0;
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
if (m_elementActive[j] && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { if (m_elementActive[j] && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
double nu = m_formulaMatrix(kspec,j); double nu = m_formulaMatrix(kspec,j);
if (nu != 0.0) { if (nu != 0.0) {
@ -2180,7 +2180,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
if (kfound == npos) { if (kfound == npos) {
double gmin = 0.0; double gmin = 0.0;
kfound = k; kfound = k;
for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = ncTrial; kspec < m_nsp; kspec++) {
if (aw[kspec] >= 0.0) { if (aw[kspec] >= 0.0) {
size_t irxn = kspec - ncTrial; size_t irxn = kspec - ncTrial;
if (m_deltaGRxn_new[irxn] < gmin) { if (m_deltaGRxn_new[irxn] < gmin) {
@ -2197,18 +2197,18 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
m_numComponents = jr; m_numComponents = jr;
ncTrial = m_numComponents; ncTrial = m_numComponents;
size_t numPreDeleted = m_numRxnTot - m_numRxnRdc; size_t numPreDeleted = m_numRxnTot - m_numRxnRdc;
if (numPreDeleted != (m_numSpeciesTot - m_numSpeciesRdc)) { if (numPreDeleted != (m_nsp - m_numSpeciesRdc)) {
throw CanteraError("VCS_SOLVE::vcs_basopt", "we shouldn't be here"); throw CanteraError("VCS_SOLVE::vcs_basopt", "we shouldn't be here");
} }
m_numRxnTot = m_numSpeciesTot - ncTrial; m_numRxnTot = m_nsp - ncTrial;
m_numRxnRdc = m_numRxnTot - numPreDeleted; m_numRxnRdc = m_numRxnTot - numPreDeleted;
m_numSpeciesRdc = m_numSpeciesTot - numPreDeleted; m_numSpeciesRdc = m_nsp - numPreDeleted;
for (size_t i = 0; i < m_numSpeciesTot; ++i) { for (size_t i = 0; i < m_nsp; ++i) {
m_indexRxnToSpecies[i] = ncTrial + i; m_indexRxnToSpecies[i] = ncTrial + i;
} }
if (m_debug_print_lvl >= 2) { if (m_debug_print_lvl >= 2) {
plogf(" --- Total number of components found = %3d (ne = %d)\n ", plogf(" --- Total number of components found = %3d (ne = %d)\n ",
ncTrial, m_numElemConstraints); ncTrial, m_nelem);
} }
goto L_END_LOOP; goto L_END_LOOP;
} }
@ -2222,8 +2222,8 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
// Modified Gram-Schmidt Method, p. 202 Dalquist // Modified Gram-Schmidt Method, p. 202 Dalquist
// QR factorization of a matrix without row pivoting. // QR factorization of a matrix without row pivoting.
size_t jl = jr; size_t jl = jr;
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
sm[j + jr*m_numElemConstraints] = m_formulaMatrix(k,j); sm[j + jr*m_nelem] = m_formulaMatrix(k,j);
} }
if (jl > 0) { if (jl > 0) {
// Compute the coefficients of JA column of the the upper // Compute the coefficients of JA column of the the upper
@ -2231,16 +2231,16 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
// different than Dalquist) R_JA_JA = 1 // different than Dalquist) R_JA_JA = 1
for (size_t j = 0; j < jl; ++j) { for (size_t j = 0; j < jl; ++j) {
ss[j] = 0.0; ss[j] = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
ss[j] += sm[i + jr*m_numElemConstraints] * sm[i + j*m_numElemConstraints]; ss[j] += sm[i + jr*m_nelem] * sm[i + j*m_nelem];
} }
ss[j] /= sa[j]; ss[j] /= sa[j];
} }
// Now make the new column, (*,JR), orthogonal to the previous // Now make the new column, (*,JR), orthogonal to the previous
// columns // columns
for (size_t j = 0; j < jl; ++j) { for (size_t j = 0; j < jl; ++j) {
for (size_t i = 0; i < m_numElemConstraints; ++i) { for (size_t i = 0; i < m_nelem; ++i) {
sm[i + jr*m_numElemConstraints] -= ss[j] * sm[i + j*m_numElemConstraints]; sm[i + jr*m_nelem] -= ss[j] * sm[i + j*m_nelem];
} }
} }
} }
@ -2248,8 +2248,8 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
// Find the new length of the new column in Q. It will be used in // Find the new length of the new column in Q. It will be used in
// the denominator in future row calcs. // the denominator in future row calcs.
sa[jr] = 0.0; sa[jr] = 0.0;
for (size_t ml = 0; ml < m_numElemConstraints; ++ml) { for (size_t ml = 0; ml < m_nelem; ++ml) {
sa[jr] += pow(sm[ml + jr*m_numElemConstraints], 2); sa[jr] += pow(sm[ml + jr*m_nelem], 2);
} }
// IF NORM OF NEW ROW .LT. 1E-3 REJECT // IF NORM OF NEW ROW .LT. 1E-3 REJECT
@ -2337,23 +2337,23 @@ L_END_LOOP:
} }
// Solve the linear system to calculate the reaction matrix, // Solve the linear system to calculate the reaction matrix,
// m_stoichCoeffRxnMatrix. // m_stoichCoeffRxnMatrix.
solve(C, m_stoichCoeffRxnMatrix.ptrColumn(0), m_numRxnTot, m_numElemConstraints); solve(C, m_stoichCoeffRxnMatrix.ptrColumn(0), m_numRxnTot, m_nelem);
// NOW, if we have interfacial voltage unknowns, what we did was just wrong // NOW, if we have interfacial voltage unknowns, what we did was just wrong
// -> hopefully it didn't blow up. Redo the problem. Search for inactive E // -> hopefully it didn't blow up. Redo the problem. Search for inactive E
juse = npos; juse = npos;
jlose = npos; jlose = npos;
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
if (!m_elementActive[j] && !strcmp(m_elementName[j].c_str(), "E")) { if (!m_elementActive[j] && !strcmp(m_elementName[j].c_str(), "E")) {
juse = j; juse = j;
} }
} }
for (size_t j = 0; j < m_numElemConstraints; j++) { for (size_t j = 0; j < m_nelem; j++) {
if (m_elementActive[j] && !strncmp((m_elementName[j]).c_str(), "cn_", 3)) { if (m_elementActive[j] && !strncmp((m_elementName[j]).c_str(), "cn_", 3)) {
jlose = j; jlose = j;
} }
} }
for (k = 0; k < m_numSpeciesTot; k++) { for (k = 0; k < m_nsp; k++) {
if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
for (size_t j = 0; j < ncTrial; ++j) { for (size_t j = 0; j < ncTrial; ++j) {
for (size_t i = 0; i < ncTrial; ++i) { for (size_t i = 0; i < ncTrial; ++i) {
@ -2375,7 +2375,7 @@ L_END_LOOP:
} }
} }
solve(C, aw, 1, m_numElemConstraints); solve(C, aw, 1, m_nelem);
size_t i = k - ncTrial; size_t i = k - ncTrial;
for (size_t j = 0; j < ncTrial; j++) { for (size_t j = 0; j < ncTrial; j++) {
m_stoichCoeffRxnMatrix(j,i) = aw[j]; m_stoichCoeffRxnMatrix(j,i) = aw[j];
@ -2580,7 +2580,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
// see if the species has an element which is so low that species will // see if the species has an element which is so low that species will
// always be zero // always be zero
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
double atomComp = m_formulaMatrix(kspec,j); double atomComp = m_formulaMatrix(kspec,j);
if (atomComp > 0.0) { if (atomComp > 0.0) {
@ -2767,7 +2767,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
tlogMoles[iph] = tPhInertMoles[iph]; tlogMoles[iph] = tPhInertMoles[iph];
} }
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
size_t iph = m_phaseID[kspec]; size_t iph = m_phaseID[kspec];
tlogMoles[iph] += molNum[kspec]; tlogMoles[iph] += molNum[kspec];
@ -2964,7 +2964,7 @@ double VCS_SOLVE::vcs_tmoles()
for (size_t i = 0; i < m_numPhases; i++) { for (size_t i = 0; i < m_numPhases; i++) {
m_tPhaseMoles_old[i] = TPhInertMoles[i]; m_tPhaseMoles_old[i] = TPhInertMoles[i];
} }
for (size_t i = 0; i < m_numSpeciesTot; i++) { for (size_t i = 0; i < m_nsp; i++) {
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
m_tPhaseMoles_old[m_phaseID[i]] += m_molNumSpecies_old[i]; m_tPhaseMoles_old[m_phaseID[i]] += m_molNumSpecies_old[i];
} }
@ -2989,7 +2989,7 @@ void VCS_SOLVE::check_tmoles() const
for (size_t i = 0; i < m_numPhases; i++) { for (size_t i = 0; i < m_numPhases; i++) {
double m_tPhaseMoles_old_a = TPhInertMoles[i]; double m_tPhaseMoles_old_a = TPhInertMoles[i];
for (size_t k = 0; k < m_numSpeciesTot; k++) { for (size_t k = 0; k < m_nsp; k++) {
if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM && m_phaseID[k] == i) { if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM && m_phaseID[k] == i) {
m_tPhaseMoles_old_a += m_molNumSpecies_old[k]; m_tPhaseMoles_old_a += m_molNumSpecies_old[k];
} }
@ -3031,7 +3031,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType()
} else if (m_debug_print_lvl >= 5) { } else if (m_debug_print_lvl >= 5) {
plogf(" --- Species Status decision is reevaluated\n"); plogf(" --- Species Status decision is reevaluated\n");
} }
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) { for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
m_speciesStatus[kspec] = vcs_species_type(kspec); m_speciesStatus[kspec] = vcs_species_type(kspec);
if (m_debug_print_lvl >= 5) { if (m_debug_print_lvl >= 5) {
plogf(" --- %-16s: ", m_speciesName[kspec]); plogf(" --- %-16s: ", m_speciesName[kspec]);
@ -3326,7 +3326,7 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc)
writelog(" "); writelog(" ");
writeline('-', 132); writeline('-', 132);
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) {
size_t irxn = npos; size_t irxn = npos;
if (kspec >= m_numComponents) { if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents; irxn = kspec - m_numComponents;
@ -3402,7 +3402,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
if (k1 == k2) { if (k1 == k2) {
return; return;
} }
if (k1 >= m_numSpeciesTot || k2 >= m_numSpeciesTot) { if (k1 >= m_nsp || k2 >= m_nsp) {
plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n", plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n",
k1, k2); k1, k2);
} }
@ -3440,14 +3440,14 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
std::swap(m_speciesThermoList[k1], m_speciesThermoList[k2]); std::swap(m_speciesThermoList[k1], m_speciesThermoList[k2]);
std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]); std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]);
for (size_t j = 0; j < m_numElemConstraints; ++j) { for (size_t j = 0; j < m_nelem; ++j) {
std::swap(m_formulaMatrix(k1,j), m_formulaMatrix(k2,j)); std::swap(m_formulaMatrix(k1,j), m_formulaMatrix(k2,j));
} }
if (m_useActCoeffJac && k1 != k2) { if (m_useActCoeffJac && k1 != k2) {
for (size_t i = 0; i < m_numSpeciesTot; i++) { for (size_t i = 0; i < m_nsp; i++) {
std::swap(m_np_dLnActCoeffdMolNum(k1,i), m_np_dLnActCoeffdMolNum(k2,i)); std::swap(m_np_dLnActCoeffdMolNum(k1,i), m_np_dLnActCoeffdMolNum(k2,i));
} }
for (size_t i = 0; i < m_numSpeciesTot; i++) { for (size_t i = 0; i < m_nsp; i++) {
std::swap(m_np_dLnActCoeffdMolNum(i,k1), m_np_dLnActCoeffdMolNum(i,k2)); std::swap(m_np_dLnActCoeffdMolNum(i,k1), m_np_dLnActCoeffdMolNum(i,k2));
} }
} }