Cantera  2.5.1
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
4 // at https://cantera.org/license.txt for license and copyright information.
5 
9 
10 namespace Cantera
11 {
12 int VCS_SOLVE::vcs_report(int iconv)
13 {
14  bool inertYes = false;
15 
16  // SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES
17  std::vector<std::pair<double, size_t>> x_order;
18  for (size_t i = 0; i < m_nsp; i++) {
19  x_order.push_back({-m_molNumSpecies_old[i], i});
20  }
21  std::sort(x_order.begin() + m_numComponents,
22  x_order.begin() + m_numSpeciesRdc);
23 
24  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
26 
27  // PRINT OUT RESULTS
28  plogf("\n\n\n\n");
29  writeline('-', 80);
30  writeline('-', 80);
31  plogf("\t\t VCS_TP REPORT\n");
32  writeline('-', 80);
33  writeline('-', 80);
34  if (iconv < 0) {
35  plogf(" ERROR: CONVERGENCE CRITERION NOT SATISFIED.\n");
36  } else if (iconv == 1) {
37  plogf(" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
38  }
39 
40  // Calculate some quantities that may need updating
41  vcs_tmoles();
44 
45  plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
46  plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
47  plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol);
48 
49  // TABLE OF SPECIES IN DECREASING MOLE NUMBERS
50  plogf("\n\n");
51  writeline('-', 80);
52  plogf(" Species Equilibrium kmoles ");
53  plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
54  writeline('-', 80);
55  for (size_t i = 0; i < m_numComponents; ++i) {
56  plogf(" %-12.12s", m_speciesName[i]);
57  writeline(' ', 13, false);
58  plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i],
60  plogf(" %3d", m_speciesUnknownType[i]);
61  plogf("\n");
62  }
63  for (size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
64  size_t j = x_order[i].second;
65  plogf(" %-12.12s", m_speciesName[j]);
66  writeline(' ', 13, false);
67 
69  plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j],
71  plogf(" KMolNum ");
73  plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[j]);
74  plogf(" Voltage = %14.7E", m_molNumSpecies_old[j]);
75  } else {
76  throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
77  }
78  plogf("\n");
79  }
80  for (size_t i = 0; i < m_numPhases; i++) {
81  if (TPhInertMoles[i] > 0.0) {
82  inertYes = true;
83  if (i == 0) {
84  plogf(" Inert Gas Species ");
85  } else {
86  plogf(" Inert Species in phase %16s ",
87  m_VolPhaseList[i]->PhaseName);
88  }
89  plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i],
90  TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
91  }
92  }
93  if (m_numSpeciesRdc != m_nsp) {
94  plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
95  for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
96  plogf(" %-12.12s", m_speciesName[kspec]);
97  // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
98  plogf(" %14.7E %14.7E %12.4E",
99  m_molNumSpecies_old[kspec],
100  m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]);
102  plogf(" KMol_Num");
104  plogf(" Voltage");
105  } else {
106  plogf(" Unknown");
107  }
108  plogf("\n");
109  }
110  }
111  writeline('-', 80);
112  plogf("\n");
113 
114  // TABLE OF SPECIES FORMATION REACTIONS
115  writeline('-', m_numComponents*10 + 45, true, true);
116  plogf(" |ComponentID|");
117  for (size_t j = 0; j < m_numComponents; j++) {
118  plogf(" %3d", j);
119  }
120  plogf(" | |\n");
121  plogf(" | Components|");
122  for (size_t j = 0; j < m_numComponents; j++) {
123  plogf(" %10.10s", m_speciesName[j]);
124  }
125  plogf(" | |\n");
126  plogf(" NonComponent | Moles |");
127  for (size_t j = 0; j < m_numComponents; j++) {
128  plogf(" %10.3g", m_molNumSpecies_old[j]);
129  }
130  plogf(" | DG/RT Rxn |\n");
131  writeline('-', m_numComponents*10 + 45);
132  for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
133  size_t kspec = m_indexRxnToSpecies[irxn];
134  plogf(" %3d ", kspec);
135  plogf("%-10.10s", m_speciesName[kspec]);
136  plogf("|%10.3g |", m_molNumSpecies_old[kspec]);
137  for (size_t j = 0; j < m_numComponents; j++) {
138  plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,irxn));
139  }
140  plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
141  plogf("\n");
142  }
143  writeline('-', m_numComponents*10 + 45);
144  plogf("\n");
145 
146  // TABLE OF PHASE INFORMATION
147  vector_fp gaPhase(m_nelem, 0.0);
148  vector_fp gaTPhase(m_nelem, 0.0);
149  double totalMoles = 0.0;
150  double gibbsPhase = 0.0;
151  double gibbsTotal = 0.0;
152  plogf("\n\n");
153  plogf("\n");
154  writeline('-', m_nelem*10 + 58);
155  plogf(" | ElementID |");
156  for (size_t j = 0; j < m_nelem; j++) {
157  plogf(" %3d", j);
158  }
159  plogf(" | |\n");
160  plogf(" | Element |");
161  for (size_t j = 0; j < m_nelem; j++) {
162  plogf(" %10.10s", m_elementName[j]);
163  }
164  plogf(" | |\n");
165  plogf(" PhaseName |KMolTarget |");
166  for (size_t j = 0; j < m_nelem; j++) {
167  plogf(" %10.3g", m_elemAbundancesGoal[j]);
168  }
169  plogf(" | Gibbs Total |\n");
170  writeline('-', m_nelem*10 + 58);
171  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
172  plogf(" %3d ", iphase);
173  vcs_VolPhase* VPhase = m_VolPhaseList[iphase].get();
174  plogf("%-12.12s |",VPhase->PhaseName);
175  plogf("%10.3e |", m_tPhaseMoles_old[iphase]);
176  totalMoles += m_tPhaseMoles_old[iphase];
177  if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles() &&
178  !vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) {
179  throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
180  }
181  vcs_elabPhase(iphase, &gaPhase[0]);
182  for (size_t j = 0; j < m_nelem; j++) {
183  plogf(" %10.3g", gaPhase[j]);
184  gaTPhase[j] += gaPhase[j];
185  }
186  gibbsPhase = vcs_GibbsPhase(iphase, &m_molNumSpecies_old[0],
187  &m_feSpecies_old[0]);
188  gibbsTotal += gibbsPhase;
189  plogf(" | %18.11E |\n", gibbsPhase);
190  }
191  writeline('-', m_nelem*10 + 58);
192  plogf(" TOTAL |%10.3e |", totalMoles);
193  for (size_t j = 0; j < m_nelem; j++) {
194  plogf(" %10.3g", gaTPhase[j]);
195  }
196  plogf(" | %18.11E |\n", gibbsTotal);
197 
198  writeline('-', m_nelem*10 + 58);
199  plogf("\n");
200 
201  // GLOBAL SATISFACTION INFORMATION
202 
203  // Calculate the total dimensionless Gibbs Free Energy. Inert species are
204  // handled as if they had a standard free energy of zero
206  &m_tPhaseMoles_old[0]);
207  plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
208  if (inertYes) {
209  plogf("\t\t(Inert species have standard free energy of zero)\n");
210  }
211 
212  plogf("\nElemental Abundances (kmol): ");
213  plogf(" Actual Target Type ElActive\n");
214  for (size_t i = 0; i < m_nelem; ++i) {
215  writeline(' ', 26, false);
216  plogf("%-2.2s", m_elementName[i]);
217  plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]);
218  plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]);
219  }
220  plogf("\n");
221 
222  // TABLE OF SPECIES CHEM POTS
223  writeline('-', 93, true, true);
224  plogf("Chemical Potentials of the Species: (dimensionless)\n");
225 
226  plogf("\t\t(RT = %g J/kmol)\n", GasConstant * m_temperature);
227  plogf(" Name TKMoles StandStateChemPot "
228  " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
229  plogf("| (MolNum ChemPot)|");
230  writeline('-', 147, true, true);
231  for (size_t i = 0; i < m_nsp; ++i) {
232  size_t j = x_order[i].second;
233  size_t pid = m_phaseID[j];
234  plogf(" %-12.12s", m_speciesName[j]);
235  plogf(" %14.7E ", m_molNumSpecies_old[j]);
236  plogf("%14.7E ", m_SSfeSpecies[j]);
237  plogf("%14.7E ", log(m_actCoeffSpecies_old[j]));
238  double tpmoles = m_tPhaseMoles_old[pid];
239  double phi = m_phasePhi[pid];
240  double eContrib = phi * m_chargeSpecies[j] * m_Faraday_dim;
241  double lx = 0.0;
243  lx = 0.0;
244  } else {
245  if (tpmoles > 0.0 && m_molNumSpecies_old[j] > 0.0) {
246  double tmp = std::max(VCS_DELETE_MINORSPECIES_CUTOFF, m_molNumSpecies_old[j]);
247  lx = log(tmp) - log(tpmoles);
248  } else {
249  lx = m_feSpecies_old[j] - m_SSfeSpecies[j]
251  }
252  }
253  plogf("%14.7E |", lx);
254  plogf("%14.7E | ", eContrib);
255  double tmp = m_SSfeSpecies[j] + log(m_actCoeffSpecies_old[j])
256  + lx - m_lnMnaughtSpecies[j] + eContrib;
257  if (fabs(m_feSpecies_old[j] - tmp) > 1.0E-7) {
258  throw CanteraError("VCS_SOLVE::vcs_report",
259  "we have a problem - doesn't add up");
260  }
261  plogf(" %12.4E |", m_feSpecies_old[j]);
262  if (m_lnMnaughtSpecies[j] != 0.0) {
263  plogf("(%11.5E)", - m_lnMnaughtSpecies[j]);
264  } else {
265  plogf(" ");
266  }
267 
268  plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j]);
269  plogf("\n");
270  }
271  for (size_t i = 0; i < 125; i++) {
272  plogf(" ");
273  }
274  plogf(" %20.9E\n", g);
275  writeline('-', 147);
276 
277  // TABLE OF SOLUTION COUNTERS
278  plogf("\n");
279  plogf("\nCounters: Iterations Time (seconds)\n");
280  if (m_timing_print_lvl > 0) {
281  plogf(" vcs_basopt: %5d %11.5E\n",
283  plogf(" vcs_TP: %5d %11.5E\n",
285  } else {
286  plogf(" vcs_basopt: %5d %11s\n",
287  m_VCount->Basis_Opts," NA ");
288  plogf(" vcs_TP: %5d %11s\n",
289  m_VCount->Its," NA ");
290  }
291  writeline('-', 80);
292  writeline('-', 80);
293 
294  // Return a successful completion flag
295  return VCS_SUCCESS;
296 }
297 
298 void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl)
299 {
300  plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
301  if (timing_print_lvl > 0) {
302  plogf(" vcs_basopt: %5d %5d %11.5E\n",
305  plogf(" vcs_TP: %5d %5d %11.5E\n",
308  plogf(" vcs_inest: %5d %11.5E\n",
310  plogf(" vcs_TotalTime: %11.5E\n",
312  } else {
313  plogf(" vcs_basopt: %5d %5d %11s\n",
315  plogf(" vcs_TP: %5d %5d %11s\n",
317  plogf(" vcs_inest: %5d %11s\n",
318  m_VCount->T_Calls_Inest, " NA ");
319  plogf(" vcs_TotalTime: %11s\n",
320  " NA ");
321  }
322 }
323 
324 }
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
int Its
Current number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
Definition: vcs_internal.h:44
int T_Calls_Inest
Current number of times the initial thermo equilibrium estimator has been called.
Definition: vcs_internal.h:54
double T_Time_basopt
Total Time spent in basopt.
Definition: vcs_internal.h:66
double T_Time_inest
Time spent in initial estimator.
Definition: vcs_internal.h:72
int T_Basis_Opts
Total number of optimizations of the components basis set done.
Definition: vcs_internal.h:47
double T_Time_vcs_TP
Current time spent in vcs_TP.
Definition: vcs_internal.h:60
int Basis_Opts
number of optimizations of the components basis set done
Definition: vcs_internal.h:50
double T_Time_vcs
Time spent in the vcs suite of programs.
Definition: vcs_internal.h:75
int T_Its
Total number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
Definition: vcs_internal.h:40
int T_Calls_vcs_TP
Current number of calls to vcs_TP.
Definition: vcs_internal.h:57
double Time_basopt
Current Time spent in basopt.
Definition: vcs_internal.h:69
double Time_vcs_TP
Current time spent in vcs_TP.
Definition: vcs_internal.h:63
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1068
vector_fp m_phasePhi
electric potential of the iph phase
Definition: vcs_solve.h:1179
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:659
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1364
vector_int m_elType
Type of the element constraint.
Definition: vcs_solve.h:1385
size_t m_nelem
Number of element constraints in the problem.
Definition: vcs_solve.h:1047
double m_totalVol
Total volume of all phases. Units are m^3.
Definition: vcs_solve.h:1471
vector_fp m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1152
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:298
vector_fp m_molNumSpecies_new
Tentative value of the mole number vector.
Definition: vcs_solve.h:1183
vector_fp m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1222
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1165
int m_timing_print_lvl
printing level of timing information
Definition: vcs_solve.h:1505
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles)
Calculate the total dimensionless Gibbs free energy.
Definition: vcs_Gibbs.cpp:14
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Definition: vcs_elem.cpp:85
vector_fp m_tPhaseMoles_old
Total kmols of species in each phase.
Definition: vcs_solve.h:1246
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1057
vector_fp m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1231
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
Array2D m_stoichCoeffRxnMatrix
Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form.
Definition: vcs_solve.h:1095
VCS_COUNTERS * m_VCount
Timing and iteration counters for the vcs object.
Definition: vcs_solve.h:1484
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
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
vector_fp m_chargeSpecies
Charge of each species. Length = number of species.
Definition: vcs_solve.h:1451
vector_fp m_feSpecies_old
Free energy vector from the start of the current iteration.
Definition: vcs_solve.h:1126
int vcs_report(int iconv)
Print out a report on the state of the equilibrium problem to standard output.
Definition: vcs_report.cpp:12
vector_int m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1391
double m_Faraday_dim
dimensionless value of Faraday's constant, F / RT (1/volt)
Definition: vcs_solve.h:1481
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_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
Definition: vcs_solve.h:1119
vector_fp m_actCoeffSpecies_old
Molar-based Activity Coefficients for Species based on old mole numbers.
Definition: vcs_solve.h:1430
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
vector_fp m_PMVolumeSpecies
Partial molar volumes of the species.
Definition: vcs_solve.h:1478
vector_fp m_deltaGRxn_new
Delta G(irxn) for the noncomponent species in the mechanism.
Definition: vcs_solve.h:1193
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,...
vector_fp m_lnMnaughtSpecies
specifies the ln(Mnaught) used to calculate the chemical potentials
Definition: vcs_solve.h:1419
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
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1370
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:82
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:627
double totalMoles() const
Return the total moles in the phase.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
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
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:109
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
Definition: vcs_util.cpp:91
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:289
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:300
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem.
Definition: vcs_defs.h:44
#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...