Cantera  2.0
vcs_report.cpp
1 
2 /*
3  * Copyright (2005) Sandia Corporation. Under the terms of
4  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
5  * U.S. Government retains certain rights in this software.
6  */
7 
8 
9 #include "cantera/equil/vcs_solve.h"
10 #include "cantera/equil/vcs_internal.h"
11 #include "cantera/equil/vcs_VolPhase.h"
12 
13 #include <cstdio>
14 #include <cstdlib>
15 #include <cmath>
16 
17 namespace VCSnonideal
18 {
19 
20 /*****************************************************************************/
21 static void print_space(int num)
22 {
23  for (int j = 0; j < num; j++) {
24  plogf(" ");
25  }
26 }
27 
28 static void print_line(std::string schar, size_t num)
29 {
30  for (size_t j = 0; j < num; j++) {
31  plogf("%s", schar.c_str());
32  }
33  plogf("\n");
34 }
35 
36 /*****************************************************************************/
37 
38 /**************************************************************************
39  *
40  * vcs_report:
41  *
42  * Print out a report on the state of the equilibrium problem to
43  * standard output.
44  * This prints out the current contents of the VCS_SOLVE class, V.
45  * The "old" solution vector is printed out.
46  ***************************************************************************/
47 int VCS_SOLVE::vcs_report(int iconv)
48 {
49  bool printActualMoles = true, inertYes = false;
50  size_t i, j, l, k, kspec;
51  size_t nspecies = m_numSpeciesTot;
52  double g;
53 
54  char originalUnitsState = m_unitsState;
55 
56 
57  std::vector<size_t> sortindex(nspecies,0);
58  std::vector<double> xy(nspecies,0.0);
59 
60  /* ************************************************************** */
61  /* **** SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES ***** */
62  /* ************************************************************** */
63 
64  for (i = 0; i < nspecies; ++i) {
65  sortindex[i] = i;
66  xy[i] = m_molNumSpecies_old[i];
67  }
68  /*
69  * Sort the XY vector, the mole fraction vector,
70  * and the sort index vector, sortindex, according to
71  * the magnitude of the mole fraction vector.
72  */
73  for (l = m_numComponents; l < m_numSpeciesRdc; ++l) {
74  k = vcs_optMax(VCS_DATA_PTR(xy), 0, l, m_numSpeciesRdc);
75  if (k != l) {
76  std::swap(xy[k], xy[l]);
77  std::swap(sortindex[k], sortindex[l]);
78  }
79  }
80 
81  /*
82  * Decide whether we have to nondimensionalize the equations.
83  * -> For the printouts from this routine, we will use nondimensional
84  * representations. This may be expanded in the future.
85  */
86  if (m_unitsState == VCS_DIMENSIONAL_G) {
87  vcs_nondim_TP();
88  }
89  double molScale = 1.0;
90  if (printActualMoles) {
91  molScale = m_totalMoleScale;
92  }
93 
94  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
95  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
96  /* ******************************************************** */
97  /* *** PRINT OUT RESULTS ********************************** */
98  /* ******************************************************** */
99 
100  plogf("\n\n\n\n");
101  print_line("-", 80);
102  print_line("-", 80);
103  plogf("\t\t VCS_TP REPORT\n");
104  print_line("-", 80);
105  print_line("-", 80);
106  if (iconv < 0) {
107  plogf(" ERROR: CONVERGENCE CRITERION NOT SATISFIED.\n");
108  } else if (iconv == 1) {
109  plogf(" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
110  }
111  /*
112  * Calculate some quantities that may need updating
113  */
114  vcs_tmoles();
115  m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA,
116  VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies));
117 
118  plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
119  plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
120  plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol * molScale);
121  if (!printActualMoles) {
122  plogf("\t\tMole Scale = %15.5g kmol (all mole numbers and volumes are scaled by this value)\n",
123  molScale);
124  }
125 
126  /*
127  * -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS --------------
128  */
129  plogf("\n\n");
130  print_line("-", 80);
131  plogf(" Species Equilibrium kmoles ");
132  plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
133  print_line("-", 80);
134  for (i = 0; i < m_numComponents; ++i) {
135  plogf(" %-12.12s", m_speciesName[i].c_str());
136  print_space(13);
137  plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i] * molScale,
138  m_molNumSpecies_new[i] * molScale, m_feSpecies_old[i]);
139  plogf(" %3d", m_speciesUnknownType[i]);
140  plogf("\n");
141  }
142  for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
143  l = sortindex[i];
144  plogf(" %-12.12s", m_speciesName[l].c_str());
145  print_space(13);
146 
147  if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) {
148  plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[l] * molScale,
149  m_molNumSpecies_new[l] * molScale, m_feSpecies_old[l]);
150  plogf(" KMolNum ");
151  } else if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
152  plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[l]);
153  plogf(" Voltage = %14.7E", m_molNumSpecies_old[l] * molScale);
154  } else {
155  plogf("we have a problem\n");
156  exit(EXIT_FAILURE);
157  }
158  plogf("\n");
159  }
160  for (i = 0; i < m_numPhases; i++) {
161  if (TPhInertMoles[i] > 0.0) {
162  inertYes = true;
163  if (i == 0) {
164  plogf(" Inert Gas Species ");
165  } else {
166  plogf(" Inert Species in phase %16s ",
167  (m_VolPhaseList[i])->PhaseName.c_str());
168  }
169  plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i] * molScale,
170  TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
171  }
172  }
173  if (m_numSpeciesRdc != nspecies) {
174  plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
175  for (kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
176  plogf(" %-12.12s", m_speciesName[kspec].c_str());
177  // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
178  plogf(" %14.7E %14.7E %12.4E",
179  m_molNumSpecies_old[kspec]*molScale,
180  m_molNumSpecies_new[kspec]*molScale, m_deltaGRxn_new[kspec]);
181  if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
182  plogf(" KMol_Num");
183  } else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
184  plogf(" Voltage");
185  } else {
186  plogf(" Unknown");
187  }
188 
189  plogf("\n");
190  }
191  }
192  print_line("-", 80);
193  plogf("\n");
194 
195  /*
196  * ---------- TABLE OF SPECIES FORMATION REACTIONS ------------------
197  */
198  plogf("\n");
199  print_line("-", m_numComponents*10 + 45);
200  plogf(" |ComponentID|");
201  for (j = 0; j < m_numComponents; j++) {
202  plogf(" %3d", j);
203  }
204  plogf(" | |\n");
205  plogf(" | Components|");
206  for (j = 0; j < m_numComponents; j++) {
207  plogf(" %10.10s", m_speciesName[j].c_str());
208  }
209  plogf(" | |\n");
210  plogf(" NonComponent | Moles |");
211  for (j = 0; j < m_numComponents; j++) {
212  plogf(" %10.3g", m_molNumSpecies_old[j] * molScale);
213  }
214  plogf(" | DG/RT Rxn |\n");
215  print_line("-", m_numComponents*10 + 45);
216  for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
217  size_t kspec = m_indexRxnToSpecies[irxn];
218  plogf(" %3d ", kspec);
219  plogf("%-10.10s", m_speciesName[kspec].c_str());
220  plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale);
221  for (j = 0; j < m_numComponents; j++) {
222  plogf(" %6.2f", m_stoichCoeffRxnMatrix[irxn][j]);
223  }
224  plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
225  plogf("\n");
226  }
227  print_line("-", m_numComponents*10 + 45);
228  plogf("\n");
229 
230  /*
231  * ------------------ TABLE OF PHASE INFORMATION ---------------------
232  */
233  std::vector<double> gaPhase(m_numElemConstraints, 0.0);
234  std::vector<double> gaTPhase(m_numElemConstraints, 0.0);
235  double totalMoles = 0.0;
236  double gibbsPhase = 0.0;
237  double gibbsTotal = 0.0;
238  plogf("\n\n");
239  plogf("\n");
240  print_line("-", m_numElemConstraints*10 + 58);
241  plogf(" | ElementID |");
242  for (j = 0; j < m_numElemConstraints; j++) {
243  plogf(" %3d", j);
244  }
245  plogf(" | |\n");
246  plogf(" | Element |");
247  for (j = 0; j < m_numElemConstraints; j++) {
248  plogf(" %10.10s", (m_elementName[j]).c_str());
249  }
250  plogf(" | |\n");
251  plogf(" PhaseName |KMolTarget |");
252  for (j = 0; j < m_numElemConstraints; j++) {
253  plogf(" %10.3g", m_elemAbundancesGoal[j]);
254  }
255  plogf(" | Gibbs Total |\n");
256  print_line("-", m_numElemConstraints*10 + 58);
257  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
258  plogf(" %3d ", iphase);
259  vcs_VolPhase* VPhase = m_VolPhaseList[iphase];
260  plogf("%-12.12s |",VPhase->PhaseName.c_str());
261  plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale);
262  totalMoles += m_tPhaseMoles_old[iphase];
263  if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles()) {
264  if (! vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) {
265  plogf("We have a problem\n");
266  exit(EXIT_FAILURE);
267  }
268  }
269  vcs_elabPhase(iphase, VCS_DATA_PTR(gaPhase));
270  for (j = 0; j < m_numElemConstraints; j++) {
271  plogf(" %10.3g", gaPhase[j]);
272  gaTPhase[j] += gaPhase[j];
273  }
274  gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old),
275  VCS_DATA_PTR(m_feSpecies_old));
276  gibbsTotal += gibbsPhase;
277  plogf(" | %18.11E |\n", gibbsPhase);
278  }
279  print_line("-", m_numElemConstraints*10 + 58);
280  plogf(" TOTAL |%10.3e |", totalMoles);
281  for (j = 0; j < m_numElemConstraints; j++) {
282  plogf(" %10.3g", gaTPhase[j]);
283  }
284  plogf(" | %18.11E |\n", gibbsTotal);
285 
286  print_line("-", m_numElemConstraints*10 + 58);
287  plogf("\n");
288 
289  /*
290  * ----------- GLOBAL SATISFACTION INFORMATION -----------------------
291  */
292 
293  /*
294  * Calculate the total dimensionless Gibbs Free Energy
295  * -> Inert species are handled as if they had a standard free
296  * energy of zero
297  */
298 
299  g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
300  VCS_DATA_PTR(m_tPhaseMoles_old));
301  plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
302  if (inertYes) {
303  plogf("\t\t(Inert species have standard free energy of zero)\n");
304  }
305 
306  plogf("\nElemental Abundances (kmol): ");
307  plogf(" Actual Target Type ElActive\n");
308  for (i = 0; i < m_numElemConstraints; ++i) {
309  print_space(26);
310  plogf("%-2.2s", (m_elementName[i]).c_str());
311  plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale);
312  plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]);
313  }
314  plogf("\n");
315 
316  /*
317  * ------------------ TABLE OF SPECIES CHEM POTS ---------------------
318  */
319  plogf("\n");
320  print_line("-", 93);
321  plogf("Chemical Potentials of the Species: (dimensionless)\n");
322 
323  double rt = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
324  plogf("\t\t(RT = %g ", rt);
325  vcs_printChemPotUnits(m_VCS_UnitsFormat);
326  plogf(")\n");
327  plogf(" Name TKMoles StandStateChemPot "
328  " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
329  plogf("| (MolNum ChemPot)|");
330  plogf("\n");
331  print_line("-", 147);
332  for (i = 0; i < nspecies; ++i) {
333  l = sortindex[i];
334  size_t pid = m_phaseID[l];
335  plogf(" %-12.12s", m_speciesName[l].c_str());
336  plogf(" %14.7E ", m_molNumSpecies_old[l]*molScale);
337  plogf("%14.7E ", m_SSfeSpecies[l]);
338  plogf("%14.7E ", log(m_actCoeffSpecies_old[l]));
339  double tpmoles = m_tPhaseMoles_old[pid];
340  double phi = m_phasePhi[pid];
341  double eContrib = phi * m_chargeSpecies[l] * m_Faraday_dim;
342  double lx = 0.0;
343  if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
344  lx = 0.0;
345  } else {
346  if (tpmoles > 0.0 && m_molNumSpecies_old[l] > 0.0) {
347  double tmp = std::max(VCS_DELETE_MINORSPECIES_CUTOFF, m_molNumSpecies_old[l]);
348  lx = log(tmp) - log(tpmoles);
349  } else {
350  lx = m_feSpecies_old[l] - m_SSfeSpecies[l]
351  - log(m_actCoeffSpecies_old[l]) + m_lnMnaughtSpecies[l];
352  }
353  }
354  plogf("%14.7E |", lx);
355  plogf("%14.7E | ", eContrib);
356  double tmp = m_SSfeSpecies[l] + log(m_actCoeffSpecies_old[l])
357  + lx - m_lnMnaughtSpecies[l] + eContrib;
358  if (fabs(m_feSpecies_old[l] - tmp) > 1.0E-7) {
359  plogf("\n\t\twe have a problem - doesn't add up\n");
360  exit(EXIT_FAILURE);
361  }
362  plogf(" %12.4E |", m_feSpecies_old[l]);
363  if (m_lnMnaughtSpecies[l] != 0.0) {
364  plogf("(%11.5E)", - m_lnMnaughtSpecies[l]);
365  } else {
366  plogf(" ");
367  }
368 
369  plogf("| %20.9E |", m_feSpecies_old[l] * m_molNumSpecies_old[l] * molScale);
370  plogf("\n");
371  }
372  for (i = 0; i < 125; i++) {
373  plogf(" ");
374  }
375  plogf(" %20.9E\n", g);
376  print_line("-", 147);
377 
378  /*
379  * ------------- TABLE OF SOLUTION COUNTERS --------------------------
380  */
381  plogf("\n");
382  plogf("\nCounters: Iterations Time (seconds)\n");
383  if (m_timing_print_lvl > 0) {
384  plogf(" vcs_basopt: %5d %11.5E\n",
385  m_VCount->Basis_Opts, m_VCount->Time_basopt);
386  plogf(" vcs_TP: %5d %11.5E\n",
387  m_VCount->Its, m_VCount->Time_vcs_TP);
388  } else {
389  plogf(" vcs_basopt: %5d %11s\n",
390  m_VCount->Basis_Opts," NA ");
391  plogf(" vcs_TP: %5d %11s\n",
392  m_VCount->Its," NA ");
393  }
394  print_line("-", 80);
395  print_line("-", 80);
396 
397  /*
398  * Set the Units state of the system back to where it was when we
399  * entered the program.
400  */
401  if (originalUnitsState != m_unitsState) {
402  if (originalUnitsState == VCS_DIMENSIONAL_G) {
403  vcs_redim_TP();
404  } else {
405  vcs_nondim_TP();
406  }
407  }
408  /*
409  * Return a successful completion flag
410  */
411  return VCS_SUCCESS;
412 } /* vcs_report() ************************************************************/
413 /*****************************************************************************/
414 /*****************************************************************************/
415 /*****************************************************************************/
416 
417 void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl)
418 
419 /**************************************************************************
420  *
421  * vcs_TCounters_report:
422  *
423  * Print out the total Its and time counters to standard output
424  ***************************************************************************/
425 {
426  plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
427  if (timing_print_lvl > 0) {
428  plogf(" vcs_basopt: %5d %5d %11.5E\n",
429  m_VCount->T_Basis_Opts, m_VCount->T_Basis_Opts,
430  m_VCount->T_Time_basopt);
431  plogf(" vcs_TP: %5d %5d %11.5E\n",
432  m_VCount->T_Calls_vcs_TP, m_VCount->T_Its,
433  m_VCount->T_Time_vcs_TP);
434  plogf(" vcs_inest: %5d %11.5E\n",
435  m_VCount->T_Calls_Inest, m_VCount->T_Time_inest);
436  plogf(" vcs_TotalTime: %11.5E\n",
437  m_VCount->T_Time_vcs);
438  } else {
439  plogf(" vcs_basopt: %5d %5d %11s\n",
440  m_VCount->T_Basis_Opts, m_VCount->T_Basis_Opts," NA ");
441  plogf(" vcs_TP: %5d %5d %11s\n",
442  m_VCount->T_Calls_vcs_TP, m_VCount->T_Its," NA ");
443  plogf(" vcs_inest: %5d %11s\n",
444  m_VCount->T_Calls_Inest, " NA ");
445  plogf(" vcs_TotalTime: %11s\n",
446  " NA ");
447  }
448 }
449 
450 /*****************************************************************************/
451 /*****************************************************************************/
452 /*****************************************************************************/
453 }
454 
Generated by   doxygen 1.8.2