Cantera  2.3.0
vcs_report.cpp
Go to the documentation of this file.
1 //! @file vcs_report.cpp
2
3 // This file is part of Cantera. See License.txt in the top-level directory or
5
9
10 namespace Cantera
11 {
12 int VCS_SOLVE::vcs_report(int iconv)
13 {
14  bool printActualMoles = true, inertYes = false;
15  size_t nspecies = m_numSpeciesTot;
16  char originalUnitsState = m_unitsState;
17  std::vector<size_t> sortindex(nspecies,0);
18  vector_fp xy(nspecies,0.0);
19
20  // SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES
21  for (size_t i = 0; i < nspecies; ++i) {
22  sortindex[i] = i;
23  xy[i] = m_molNumSpecies_old[i];
24  }
25
26  // Sort the XY vector, the mole fraction vector, and the sort index vector,
27  // sortindex, according to the magnitude of the mole fraction vector.
28  for (size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
29  size_t k = vcs_optMax(&xy, 0, i, m_numSpeciesRdc);
30  if (k != i) {
31  std::swap(xy[k], xy[i]);
32  std::swap(sortindex[k], sortindex[i]);
33  }
34  }
35
36  // Decide whether we have to nondimensionalize the equations. For the
37  // printouts from this routine, we will use nondimensional representations.
38  // This may be expanded in the future.
40  vcs_nondim_TP();
41  }
42  double molScale = 1.0;
43  if (printActualMoles) {
44  molScale = m_totalMoleScale;
45  }
46
47  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
49
50  // PRINT OUT RESULTS
51  plogf("\n\n\n\n");
52  writeline('-', 80);
53  writeline('-', 80);
54  plogf("\t\t VCS_TP REPORT\n");
55  writeline('-', 80);
56  writeline('-', 80);
57  if (iconv < 0) {
58  plogf(" ERROR: CONVERGENCE CRITERION NOT SATISFIED.\n");
59  } else if (iconv == 1) {
60  plogf(" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
61  }
62
63  // Calculate some quantities that may need updating
64  vcs_tmoles();
67
68  plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
69  plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
70  plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol * molScale);
71  if (!printActualMoles) {
72  plogf("\t\tMole Scale = %15.5g kmol (all mole numbers and volumes are scaled by this value)\n",
73  molScale);
74  }
75
76  // TABLE OF SPECIES IN DECREASING MOLE NUMBERS
77  plogf("\n\n");
78  writeline('-', 80);
79  plogf(" Species Equilibrium kmoles ");
80  plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
81  writeline('-', 80);
82  for (size_t i = 0; i < m_numComponents; ++i) {
83  plogf(" %-12.12s", m_speciesName[i]);
84  writeline(' ', 13, false);
85  plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i] * molScale,
86  m_molNumSpecies_new[i] * molScale, m_feSpecies_old[i]);
87  plogf(" %3d", m_speciesUnknownType[i]);
88  plogf("\n");
89  }
90  for (size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
91  size_t j = sortindex[i];
92  plogf(" %-12.12s", m_speciesName[j]);
93  writeline(' ', 13, false);
94
96  plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j] * molScale,
97  m_molNumSpecies_new[j] * molScale, m_feSpecies_old[j]);
98  plogf(" KMolNum ");
100  plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[j]);
101  plogf(" Voltage = %14.7E", m_molNumSpecies_old[j] * molScale);
102  } else {
103  throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
104  }
105  plogf("\n");
106  }
107  for (size_t i = 0; i < m_numPhases; i++) {
108  if (TPhInertMoles[i] > 0.0) {
109  inertYes = true;
110  if (i == 0) {
111  plogf(" Inert Gas Species ");
112  } else {
113  plogf(" Inert Species in phase %16s ",
114  m_VolPhaseList[i]->PhaseName);
115  }
116  plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i] * molScale,
117  TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
118  }
119  }
120  if (m_numSpeciesRdc != nspecies) {
121  plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
122  for (size_t kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
123  plogf(" %-12.12s", m_speciesName[kspec]);
124  // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
125  plogf(" %14.7E %14.7E %12.4E",
126  m_molNumSpecies_old[kspec]*molScale,
127  m_molNumSpecies_new[kspec]*molScale, m_deltaGRxn_new[kspec]);
129  plogf(" KMol_Num");
131  plogf(" Voltage");
132  } else {
133  plogf(" Unknown");
134  }
135  plogf("\n");
136  }
137  }
138  writeline('-', 80);
139  plogf("\n");
140
141  // TABLE OF SPECIES FORMATION REACTIONS
142  writeline('-', m_numComponents*10 + 45, true, true);
143  plogf(" |ComponentID|");
144  for (size_t j = 0; j < m_numComponents; j++) {
145  plogf(" %3d", j);
146  }
147  plogf(" | |\n");
148  plogf(" | Components|");
149  for (size_t j = 0; j < m_numComponents; j++) {
150  plogf(" %10.10s", m_speciesName[j]);
151  }
152  plogf(" | |\n");
153  plogf(" NonComponent | Moles |");
154  for (size_t j = 0; j < m_numComponents; j++) {
155  plogf(" %10.3g", m_molNumSpecies_old[j] * molScale);
156  }
157  plogf(" | DG/RT Rxn |\n");
158  writeline('-', m_numComponents*10 + 45);
159  for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
160  size_t kspec = m_indexRxnToSpecies[irxn];
161  plogf(" %3d ", kspec);
162  plogf("%-10.10s", m_speciesName[kspec]);
163  plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale);
164  for (size_t j = 0; j < m_numComponents; j++) {
165  plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,irxn));
166  }
167  plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
168  plogf("\n");
169  }
170  writeline('-', m_numComponents*10 + 45);
171  plogf("\n");
172
173  // TABLE OF PHASE INFORMATION
174  vector_fp gaPhase(m_numElemConstraints, 0.0);
175  vector_fp gaTPhase(m_numElemConstraints, 0.0);
176  double totalMoles = 0.0;
177  double gibbsPhase = 0.0;
178  double gibbsTotal = 0.0;
179  plogf("\n\n");
180  plogf("\n");
181  writeline('-', m_numElemConstraints*10 + 58);
182  plogf(" | ElementID |");
183  for (size_t j = 0; j < m_numElemConstraints; j++) {
184  plogf(" %3d", j);
185  }
186  plogf(" | |\n");
187  plogf(" | Element |");
188  for (size_t j = 0; j < m_numElemConstraints; j++) {
189  plogf(" %10.10s", m_elementName[j]);
190  }
191  plogf(" | |\n");
192  plogf(" PhaseName |KMolTarget |");
193  for (size_t j = 0; j < m_numElemConstraints; j++) {
194  plogf(" %10.3g", m_elemAbundancesGoal[j]);
195  }
196  plogf(" | Gibbs Total |\n");
197  writeline('-', m_numElemConstraints*10 + 58);
198  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
199  plogf(" %3d ", iphase);
200  vcs_VolPhase* VPhase = m_VolPhaseList[iphase];
201  plogf("%-12.12s |",VPhase->PhaseName);
202  plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale);
203  totalMoles += m_tPhaseMoles_old[iphase];
204  if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles() &&
205  !vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) {
206  throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
207  }
208  vcs_elabPhase(iphase, &gaPhase);
209  for (size_t j = 0; j < m_numElemConstraints; j++) {
210  plogf(" %10.3g", gaPhase[j]);
211  gaTPhase[j] += gaPhase[j];
212  }
213  gibbsPhase = vcs_GibbsPhase(iphase, &m_molNumSpecies_old,
214  &m_feSpecies_old);
215  gibbsTotal += gibbsPhase;
216  plogf(" | %18.11E |\n", gibbsPhase);
217  }
218  writeline('-', m_numElemConstraints*10 + 58);
219  plogf(" TOTAL |%10.3e |", totalMoles);
220  for (size_t j = 0; j < m_numElemConstraints; j++) {
221  plogf(" %10.3g", gaTPhase[j]);
222  }
223  plogf(" | %18.11E |\n", gibbsTotal);
224
225  writeline('-', m_numElemConstraints*10 + 58);
226  plogf("\n");
227
228  // GLOBAL SATISFACTION INFORMATION
229
230  // Calculate the total dimensionless Gibbs Free Energy. Inert species are
231  // handled as if they had a standard free energy of zero
233  &m_tPhaseMoles_old);
234  plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
235  if (inertYes) {
236  plogf("\t\t(Inert species have standard free energy of zero)\n");
237  }
238
239  plogf("\nElemental Abundances (kmol): ");
240  plogf(" Actual Target Type ElActive\n");
241  for (size_t i = 0; i < m_numElemConstraints; ++i) {
242  writeline(' ', 26, false);
243  plogf("%-2.2s", m_elementName[i]);
244  plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale);
245  plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]);
246  }
247  plogf("\n");
248
249  // TABLE OF SPECIES CHEM POTS
250  writeline('-', 93, true, true);
251  plogf("Chemical Potentials of the Species: (dimensionless)\n");
252
253  plogf("\t\t(RT = %g J/kmol)\n", GasConstant * m_temperature);
254  plogf(" Name TKMoles StandStateChemPot "
255  " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
256  plogf("| (MolNum ChemPot)|");
257  writeline('-', 147, true, true);
258  for (size_t i = 0; i < nspecies; ++i) {
259  size_t j = sortindex[i];
260  size_t pid = m_phaseID[j];
261  plogf(" %-12.12s", m_speciesName[j]);
262  plogf(" %14.7E ", m_molNumSpecies_old[j]*molScale);
263  plogf("%14.7E ", m_SSfeSpecies[j]);
264  plogf("%14.7E ", log(m_actCoeffSpecies_old[j]));
265  double tpmoles = m_tPhaseMoles_old[pid];
266  double phi = m_phasePhi[pid];
267  double eContrib = phi * m_chargeSpecies[j] * m_Faraday_dim;
268  double lx = 0.0;
270  lx = 0.0;
271  } else {
272  if (tpmoles > 0.0 && m_molNumSpecies_old[j] > 0.0) {
273  double tmp = std::max(VCS_DELETE_MINORSPECIES_CUTOFF, m_molNumSpecies_old[j]);
274  lx = log(tmp) - log(tpmoles);
275  } else {
276  lx = m_feSpecies_old[j] - m_SSfeSpecies[j]
278  }
279  }
280  plogf("%14.7E |", lx);
281  plogf("%14.7E | ", eContrib);
282  double tmp = m_SSfeSpecies[j] + log(m_actCoeffSpecies_old[j])
283  + lx - m_lnMnaughtSpecies[j] + eContrib;
284  if (fabs(m_feSpecies_old[j] - tmp) > 1.0E-7) {
285  throw CanteraError("VCS_SOLVE::vcs_report",
286  "we have a problem - doesn't add up");
287  }
288  plogf(" %12.4E |", m_feSpecies_old[j]);
289  if (m_lnMnaughtSpecies[j] != 0.0) {
290  plogf("(%11.5E)", - m_lnMnaughtSpecies[j]);
291  } else {
292  plogf(" ");
293  }
294
295  plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j] * molScale);
296  plogf("\n");
297  }
298  for (size_t i = 0; i < 125; i++) {
299  plogf(" ");
300  }
301  plogf(" %20.9E\n", g);
302  writeline('-', 147);
303
304  // TABLE OF SOLUTION COUNTERS
305  plogf("\n");
306  plogf("\nCounters: Iterations Time (seconds)\n");
307  if (m_timing_print_lvl > 0) {
308  plogf(" vcs_basopt: %5d %11.5E\n",
310  plogf(" vcs_TP: %5d %11.5E\n",
312  } else {
313  plogf(" vcs_basopt: %5d %11s\n",
314  m_VCount->Basis_Opts," NA ");
315  plogf(" vcs_TP: %5d %11s\n",
316  m_VCount->Its," NA ");
317  }
318  writeline('-', 80);
319  writeline('-', 80);
320
321  // Set the Units state of the system back to where it was when we
322  // entered the program.
323  if (originalUnitsState != m_unitsState) {
324  if (originalUnitsState == VCS_DIMENSIONAL_G) {
325  vcs_redim_TP();
326  } else {
327  vcs_nondim_TP();
328  }
329  }
330
331  // Return a successful completion flag
332  return VCS_SUCCESS;
333 }
334
335 void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl)
336 {
337  plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
338  if (timing_print_lvl > 0) {
339  plogf(" vcs_basopt: %5d %5d %11.5E\n",
342  plogf(" vcs_TP: %5d %5d %11.5E\n",
345  plogf(" vcs_inest: %5d %11.5E\n",
347  plogf(" vcs_TotalTime: %11.5E\n",
349  } else {
350  plogf(" vcs_basopt: %5d %5d %11s\n",
352  plogf(" vcs_TP: %5d %5d %11s\n",
354  plogf(" vcs_inest: %5d %11s\n",
355  m_VCount->T_Calls_Inest, " NA ");
356  plogf(" vcs_TotalTime: %11s\n",
357  " NA ");
358  }
359 }
360
361 }
vector_fp m_phasePhi
electric potential of the iph phase
Definition: vcs_solve.h:1525
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1703
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
Definition: vcs_util.cpp:116
int Its
Current number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
Definition: vcs_internal.h:52
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles)
Calculate the total dimensionless Gibbs free energy.
Definition: vcs_Gibbs.cpp:14
int T_Calls_Inest
Current number of times the initial thermo equilibrium estimator has been called. ...
Definition: vcs_internal.h:62
double Time_vcs_TP
Current time spent in vcs_TP.
Definition: vcs_internal.h:71
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1511
int Basis_Opts
number of optimizations of the components basis set done
Definition: vcs_internal.h:58
vector_fp m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1568
double T_Time_basopt
Total Time spent in basopt.
Definition: vcs_internal.h:74
dimensionless value of Faraday&#39;s constant, F / RT (1/volt)
Definition: vcs_solve.h:1846
VCS_COUNTERS * m_VCount
Timing and iteration counters for the vcs object.
Definition: vcs_solve.h:1849
vector_fp m_PMVolumeSpecies
Partial molar volumes of the species.
Definition: vcs_solve.h:1843
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
Array2D m_stoichCoeffRxnMatrix
Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form...
Definition: vcs_solve.h:1441
void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop)
Calculate the dimensionless chemical potentials of all species or of certain groups of species...
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Definition: vcs_elem.cpp:85
std::vector< vcs_VolPhase * > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1740
vector_fp m_lnMnaughtSpecies
specifies the ln(Mnaught) used to calculate the chemical potentials
Definition: vcs_solve.h:1784
void vcs_nondim_TP()
Nondimensionalize the problem data.
Definition: vcs_nondim.cpp:17
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:653
void vcs_redim_TP()
Redimensionalize the problem data.
Definition: vcs_nondim.cpp:91
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1403
int T_Basis_Opts
Total number of optimizations of the components basis set done.
Definition: vcs_internal.h:55
vector_fp m_deltaGRxn_new
Delta G(irxn) for the noncomponent species in the mechanism.
Definition: vcs_solve.h:1539
double vcs_GibbsPhase(size_t iphase, const double *const w, const double *const fe)
Calculate the total dimensionless Gibbs free energy of a single phase.
Definition: vcs_Gibbs.cpp:39
vector_int m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1737
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1716
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 vcs_report(int iconv)
Print out a report on the state of the equilibrium problem to standard output.
Definition: vcs_report.cpp:12
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem...
Definition: vcs_defs.h:53
double vcs_VolTotal(const double tkelvin, const double pres, const double w[], double volPM[])
Calculation of the total volume and the partial molar volumes.
Definition: vcs_solve.cpp:812
#define VCS_SUCCESS
Definition: vcs_defs.h:18
double T_Time_vcs_TP
Current time spent in vcs_TP.
Definition: vcs_internal.h:68
size_t m_numElemConstraints
Number of element constraints in the problem.
Definition: vcs_solve.h:1393
void vcs_TCounters_report(int timing_print_lvl=1)
Create a report on the plog file containing timing and its information.
Definition: vcs_report.cpp:335
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1387
int T_Its
Total number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
Definition: vcs_internal.h:48
double T_Time_inest
Time spent in initial estimator.
Definition: vcs_internal.h:80
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
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
double Time_basopt
Current Time spent in basopt.
Definition: vcs_internal.h:77
int m_timing_print_lvl
printing level of timing information
Definition: vcs_solve.h:1870
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1396
vector_fp m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1577
double totalMoles() const
Return the total moles in the phase.
double m_totalVol
Total volume of all phases. Units are m^3.
Definition: vcs_solve.h:1836
double m_temperature
Temperature (Kelvin)
Definition: vcs_solve.h:1616
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:81
#define VCS_DIMENSIONAL_G
dimensioned
Definition: vcs_defs.h:96
vector_fp m_tPhaseMoles_old
Total kmols of species in each phase.
Definition: vcs_solve.h:1592
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1710
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:320
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
size_t vcs_optMax(const double *x, const double *xSize, size_t j, size_t n)
Finds the location of the maximum component in a double vector.
Definition: vcs_util.cpp:33
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:309
double T_Time_vcs
Time spent in the vcs suite of programs.
Definition: vcs_internal.h:83
vector_int m_elType
Type of the element constraint.
Definition: vcs_solve.h:1731
Namespace for the Cantera kernel.
Definition: application.cpp:29
char m_unitsState
This specifies the current state of units for the Gibbs free energy properties in the program...
Definition: vcs_solve.h:1750
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
double m_totalMoleScale
Multiplier for the mole numbers within the nondimensional formulation.
Definition: vcs_solve.h:1759
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
vector_fp m_actCoeffSpecies_old
Molar-based Activity Coefficients for Species based on old mole numbers.
Definition: vcs_solve.h:1795
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
vector_fp m_chargeSpecies
Charge of each species. Length = number of species.
Definition: vcs_solve.h:1816
int T_Calls_vcs_TP
Current number of calls to vcs_TP.
Definition: vcs_internal.h:65