Cantera  2.5.1
vcs_prep.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_prep.cpp
3  * This file contains some prepatory functions.
4  */
5 
6 // This file is part of Cantera. See License.txt in the top-level directory or
7 // at https://cantera.org/license.txt for license and copyright information.
8 
12 
13 namespace Cantera
14 {
16 {
17  vector_int numPhSpecies(m_numPhases, 0);
18  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
19  numPhSpecies[m_phaseID[kspec]]++;
20  }
21 
22  // Handle the special case of a single species in a phase that has been
23  // earmarked as a multispecies phase. Treat that species as a single-species
24  // phase
25  for (size_t iph = 0; iph < m_numPhases; iph++) {
26  vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
27  Vphase->m_singleSpecies = false;
28  if (TPhInertMoles[iph] > 0.0) {
29  Vphase->setExistence(2);
30  }
31  if (numPhSpecies[iph] <= 1 && TPhInertMoles[iph] == 0.0) {
32  Vphase->m_singleSpecies = true;
33  }
34  }
35 
36  // Fill in some useful arrays here that have to do with the static
37  // information concerning the phase ID of species. SSPhase = Boolean
38  // indicating whether a species is in a single species phase or not.
39  for (size_t kspec = 0; kspec < m_nsp; kspec++) {
40  size_t iph = m_phaseID[kspec];
41  vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
42  if (Vphase->m_singleSpecies) {
43  m_SSPhase[kspec] = true;
44  } else {
45  m_SSPhase[kspec] = false;
46  }
47  }
48 }
49 
50 int VCS_SOLVE::vcs_prep(int printLvl)
51 {
52  int retn = VCS_SUCCESS;
53  m_debug_print_lvl = printLvl;
54 
55  // Calculate the Single Species status of phases
56  // Also calculate the number of species per phase
57  vcs_SSPhase();
58 
59  // Set an initial estimate for the number of noncomponent species equal to
60  // nspecies - nelements. This may be changed below
61  if (m_nelem > m_nsp) {
62  m_numRxnTot = 0;
63  } else {
65  }
68  for (size_t i = 0; i < m_numRxnRdc; ++i) {
70  }
71 
72  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
73  size_t pID = m_phaseID[kspec];
74  size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
75  vcs_VolPhase* vPhase = m_VolPhaseList[pID].get();
76  vcs_SpeciesProperties* spProp = vPhase->speciesProperty(spPhIndex);
77  double sz = 0.0;
78  size_t eSize = spProp->FormulaMatrixCol.size();
79  for (size_t e = 0; e < eSize; e++) {
80  sz += fabs(spProp->FormulaMatrixCol[e]);
81  }
82  if (sz > 0.0) {
83  m_spSize[kspec] = sz;
84  } else {
85  m_spSize[kspec] = 1.0;
86  }
87  }
88 
89  // DETERMINE THE NUMBER OF COMPONENTS
90  //
91  // Obtain a valid estimate of the mole fraction. This will be used as an
92  // initial ordering vector for prioritizing which species are defined as
93  // components.
94  //
95  // If a mole number estimate was supplied from the input file, use that mole
96  // number estimate.
97  //
98  // If a solution estimate wasn't supplied from the input file, supply an
99  // initial estimate for the mole fractions based on the relative reverse
100  // ordering of the chemical potentials.
101  //
102  // For voltage unknowns, set these to zero for the moment.
103  double test = -1.0e-10;
104  bool modifiedSoln = false;
105  if (m_doEstimateEquil < 0) {
106  double sum = 0.0;
107  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
109  sum += fabs(m_molNumSpecies_old[kspec]);
110  }
111  }
112  if (fabs(sum) < 1.0E-6) {
113  modifiedSoln = true;
114  double pres = (m_pressurePA <= 0.0) ? 1.01325E5 : m_pressurePA;
115  retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
116  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
118  m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
119  } else {
120  m_molNumSpecies_old[kspec] = 0.0;
121  }
122  }
123  }
124  test = -1.0e20;
125  }
126 
127  // NC = number of components is in the vcs.h common block. This call to
128  // BASOPT doesn't calculate the stoichiometric reaction matrix.
129  vector_fp awSpace(m_nsp + (m_nelem + 2)*(m_nelem), 0.0);
130  double* aw = &awSpace[0];
131  if (aw == NULL) {
132  plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
133  return VCS_NOMEMORY;
134  }
135  double* sa = aw + m_nsp;
136  double* sm = sa + m_nelem;
137  double* ss = sm + m_nelem * m_nelem;
138  bool conv;
139  retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv);
140  if (retn != VCS_SUCCESS) {
141  plogf("vcs_prep_oneTime:");
142  plogf(" Determination of number of components failed: %d\n",
143  retn);
144  plogf(" Global Bailout!\n");
145  return retn;
146  }
147 
148  if (m_nsp >= m_numComponents) {
150  for (size_t i = 0; i < m_numRxnRdc; ++i) {
152  }
153  } else {
154  m_numRxnTot = m_numRxnRdc = 0;
155  }
156 
157  // The elements might need to be rearranged.
158  awSpace.resize(m_nelem + (m_nelem + 2)*m_nelem, 0.0);
159  aw = &awSpace[0];
160  sa = aw + m_nelem;
161  sm = sa + m_nelem;
162  ss = sm + m_nelem * m_nelem;
163  retn = vcs_elem_rearrange(aw, sa, sm, ss);
164  if (retn != VCS_SUCCESS) {
165  plogf("vcs_prep_oneTime:");
166  plogf(" Determination of element reordering failed: %d\n",
167  retn);
168  plogf(" Global Bailout!\n");
169  return retn;
170  }
171 
172  // If we mucked up the solution unknowns because they were all
173  // zero to start with, set them back to zero here
174  if (modifiedSoln) {
175  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
176  m_molNumSpecies_old[kspec] = 0.0;
177  }
178  }
179 
180  // Initialize various arrays in the data to zero
181  m_feSpecies_old.assign(m_feSpecies_old.size(), 0.0);
182  m_feSpecies_new.assign(m_feSpecies_new.size(), 0.0);
183  m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0);
186  m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
187  m_tPhaseMoles_new.assign(m_tPhaseMoles_new.size(), 0.0);
188 
189  // Calculate the total number of moles in all phases.
190  vcs_tmoles();
191 
192  // Check to see if the current problem is well posed.
193  double sum = 0.0;
194  for (size_t e = 0; e < m_mix->nElements(); e++) {
195  sum += m_mix->elementMoles(e);
196  }
197  if (sum < 1.0E-20) {
198  // Check to see if the current problem is well posed.
199  plogf("vcs has determined the problem is not well posed: Bailing\n");
200  return VCS_PUB_BAD;
201  }
202  return VCS_SUCCESS;
203 }
204 
205 }
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Equilfu...
void zero()
Set all of the entries to zero.
Definition: Array.h:198
doublereal elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
Definition: MultiPhase.cpp:190
size_t nElements() const
Number of elements.
Definition: MultiPhase.h:98
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1068
int m_doEstimateEquil
Setting for whether to do an initial estimate.
Definition: vcs_solve.h:1145
int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
Definition: vcs_TP.cpp:43
int vcs_elem_rearrange(double *const aw, double *const sa, double *const sm, double *const ss)
Rearrange the constraint equations represented by the Formula Matrix so that the operational ones are...
size_t m_nelem
Number of element constraints in the problem.
Definition: vcs_solve.h:1047
int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], double ss[], double test, bool *const usedZeroedSpecies)
Choose the optimum species basis for the calculations.
int vcs_prep(int printLvl)
This routine is mostly concerned with changing the private data to be consistent with what's needed f...
Definition: vcs_prep.cpp:50
vector_fp m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1152
Array2D m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise.
Definition: vcs_solve.h:1176
vector_fp m_molNumSpecies_new
Tentative value of the mole number vector.
Definition: vcs_solve.h:1183
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1165
vector_fp m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
Definition: vcs_solve.h:1255
vector_fp m_feSpecies_new
Dimensionless new free energy for all the species in the mechanism at the new tentative T,...
Definition: vcs_solve.h:1135
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1057
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1498
vector_fp TPhInertMoles
Total kmoles of inert to add to each phase.
Definition: vcs_solve.h:1280
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1357
std::vector< size_t > m_indexRxnToSpecies
Mapping between the species index for noncomponent species and the full species index.
Definition: vcs_solve.h:1346
vector_fp m_deltaPhaseMoles
Change in the total moles in each phase.
Definition: vcs_solve.h:1267
std::vector< size_t > m_speciesLocalPhaseIndex
Index that keeps track of the index of the species within the local phase.
Definition: vcs_solve.h:1320
size_t m_numRxnRdc
Current number of non-component species in the problem.
Definition: vcs_solve.h:1061
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1050
double m_pressurePA
Pressure.
Definition: vcs_solve.h:1273
void vcs_SSPhase()
Calculate the status of single species phases.
Definition: vcs_prep.cpp:15
vector_fp m_feSpecies_old
Free energy vector from the start of the current iteration.
Definition: vcs_solve.h:1126
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1361
size_t m_nsp
Total number of species in the problems.
Definition: vcs_solve.h:1044
std::vector< std::unique_ptr< vcs_VolPhase > > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1394
vector_fp m_spSize
total size of the species
Definition: vcs_solve.h:1111
vector_fp m_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
Definition: vcs_solve.h:1119
Array2D m_deltaMolNumPhase
Change in the number of moles of phase, iphase, due to the noncomponent formation reaction,...
Definition: vcs_solve.h:1172
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
size_t m_numRxnTot
Total number of non-component species in the problem.
Definition: vcs_solve.h:1053
double m_temperature
Temperature (Kelvin)
Definition: vcs_solve.h:1270
Properties of a single species.
vector_fp FormulaMatrixCol
Column of the formula matrix, comprising the element composition of the species.
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:82
vcs_SpeciesProperties * speciesProperty(const size_t kindex)
Retrieve the kth Species structure for the species belonging to this phase.
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:541
void setExistence(const int existence)
Set the existence flag in the object.
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:182
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:180
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:281
#define VCS_SUCCESS
Definition: vcs_defs.h:18
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...