Cantera  2.3.0
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 http://www.cantera.org/license.txt for license and copyright information.
8 
10 #include "cantera/equil/vcs_prob.h"
12 
13 namespace Cantera
14 {
16 {
17  vector_int numPhSpecies(m_numPhases, 0);
18  for (size_t kspec = 0; kspec < m_numSpeciesTot; ++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];
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_numSpeciesTot; kspec++) {
40  size_t iph = m_phaseID[kspec];
41  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
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_oneTime(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
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_numSpeciesTot; ++kspec) {
73  size_t pID = m_phaseID[kspec];
74  size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
75  vcs_VolPhase* vPhase = m_VolPhaseList[pID];
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_numSpeciesTot; ++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_numSpeciesTot; ++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.
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_numSpeciesTot;
136  double* sm = sa + m_numElemConstraints;
137  double* ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
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 
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_numElemConstraints + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
159  aw = &awSpace[0];
160  sa = aw + m_numElemConstraints;
161  sm = sa + m_numElemConstraints;
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_numSpeciesTot; ++kspec) {
176  m_molNumSpecies_old[kspec] = 0.0;
177  }
178  }
179  return VCS_SUCCESS;
180 }
181 
183 {
184  // Initialize various arrays in the data to zero
185  m_feSpecies_old.assign(m_feSpecies_old.size(), 0.0);
186  m_feSpecies_new.assign(m_feSpecies_new.size(), 0.0);
187  m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0);
190  m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
191  m_tPhaseMoles_new.assign(m_tPhaseMoles_new.size(), 0.0);
192 
193  // Calculate the total number of moles in all phases.
194  vcs_tmoles();
195  return VCS_SUCCESS;
196 }
197 
199 {
200  double sum = 0.0;
201  for (size_t e = 0; e < vprob->ne; e++) {
202  sum += vprob->gai[e];
203  }
204  if (sum < 1.0E-20) {
205  plogf("vcs_wellPosed: Element abundance is close to zero\n");
206  return false;
207  }
208  return true;
209 }
210 
211 }
void vcs_SSPhase()
Calculate the status of single species phases.
Definition: vcs_prep.cpp:15
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...
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1703
vector_fp FormulaMatrixCol
Column of the formula matrix, comprising the element composition of the species.
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1511
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:566
vector_fp m_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
Definition: vcs_solve.h:1465
vector_fp m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1498
std::vector< vcs_VolPhase * > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1740
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
bool vcs_wellPosed(VCS_PROB *vprob)
In this routine, we check for things that will cause the algorithm to fail.
Definition: vcs_prep.cpp:198
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:1481
vector_fp m_spSize
total size of the species
Definition: vcs_solve.h:1457
size_t m_numRxnRdc
Current number of non-component species in the problem.
Definition: vcs_solve.h:1407
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1403
void setExistence(const int existence)
Set the existence flag in the object.
int vcs_prep()
Prepare the object for solution.
Definition: vcs_prep.cpp:182
vector_fp m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
Definition: vcs_solve.h:1601
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
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:301
int m_doEstimateEquil
Setting for whether to do an initial estimate.
Definition: vcs_solve.h:1491
vector_fp gai
Element abundances for jth element.
Definition: vcs_prob.h:72
#define VCS_SUCCESS
Definition: vcs_defs.h:18
Properties of a single species.
Header for the Interface class for the vcs thermo equilibrium solver package,.
size_t m_numElemConstraints
Number of element constraints in the problem.
Definition: vcs_solve.h:1393
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1387
vector_fp m_deltaPhaseMoles
Change in the total moles in each phase.
Definition: vcs_solve.h:1613
vector_fp m_feSpecies_old
Free energy vector from the start of the current iteration.
Definition: vcs_solve.h:1472
vector_fp TPhInertMoles
Total kmoles of inert to add to each phase.
Definition: vcs_solve.h:1626
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1863
int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
Definition: vcs_TP.cpp:50
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1396
int vcs_prep_oneTime(int printLvl)
This routine is mostly concerned with changing the private data to be consistent with what&#39;s needed f...
Definition: vcs_prep.cpp:50
double m_temperature
Temperature (Kelvin)
Definition: vcs_solve.h:1616
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1707
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:81
size_t ne
Number of element constraints in the equilibrium problem.
Definition: vcs_prob.h:36
Array2D m_deltaMolNumPhase
Change in the number of moles of phase, iphase, due to the noncomponent formation reaction...
Definition: vcs_solve.h:1518
double m_pressurePA
Pressure.
Definition: vcs_solve.h:1619
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:157
vcs_SpeciesProperties * speciesProperty(const size_t kindex)
Retrieve the kth Species structure for the species belonging to this phase.
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.
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
void zero()
Set all of the entries to zero.
Definition: Array.h:220
Array2D m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise...
Definition: vcs_solve.h:1522
Namespace for the Cantera kernel.
Definition: application.cpp:29
std::vector< size_t > m_indexRxnToSpecies
Mapping between the species index for noncomponent species and the full species index.
Definition: vcs_solve.h:1692
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1414
Interface class for the vcs thermo equilibrium solver package, which generally describes the problem ...
Definition: vcs_prob.h:22
size_t m_numRxnTot
Total number of non-component species in the problem.
Definition: vcs_solve.h:1399
vector_fp m_molNumSpecies_new
Tentative value of the mole number vector.
Definition: vcs_solve.h:1529
std::vector< size_t > m_speciesLocalPhaseIndex
Index that keeps track of the index of the species within the local phase.
Definition: vcs_solve.h:1666