Cantera  2.0
vcs_inest.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_inest.cpp
3  * Implementation methods for obtaining a good initial guess
4  */
5 /*
6  * Copyright (2005) Sandia Corporation. Under the terms of
7  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
8  * U.S. Government retains certain rights in this software.
9  */
10 
11 #include "cantera/equil/vcs_solve.h"
12 #include "cantera/equil/vcs_internal.h"
13 #include "cantera/equil/vcs_VolPhase.h"
14 
15 #include "cantera/base/clockWC.h"
16 
17 #include <cstdio>
18 #include <cstdlib>
19 #include <cmath>
20 
21 namespace VCSnonideal
22 {
23 
24 static char pprefix[20] = " --- vcs_inest: ";
25 
26 // Estimate equilibrium compositions
27 /*
28  * Estimates equilibrium compositions.
29  * Algorithm covered in a section of Smith and Missen's Book.
30  *
31  * Linear programming module is based on using dbolm.
32  *
33  * @param aw aw[i[ Mole fraction work space (ne in length)
34  * @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length)
35  * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length)
36  * @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length)
37  * @param test This is a small negative number.
38  */
39 void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
40  double* const ss, double test)
41 {
42  size_t lt, ikl, kspec, iph, irxn;
43  double s;
44  double s1 = 0.0;
45  double xl, par;
46  bool finished, conv;
47  size_t nspecies = m_numSpeciesTot;
48  size_t nrxn = m_numRxnTot;
49 
50  // double *molNum = VCS_DATA_PTR(m_molNumSpecies_old);
51  double TMolesMultiphase;
52  double* xtphMax = VCS_DATA_PTR(m_TmpPhase);
53  double* xtphMin = VCS_DATA_PTR(m_TmpPhase2);
54 
55  ikl = 0;
56  lt = 0;
57 
58 
59  /*
60  * CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB
61  * AND molNum(I) .GE. 0.0
62  *
63  * Note, both of these programs do this.
64  */
65 #ifdef ALTLINPROG
66  vcs_setMolesLinProg();
67 #else
68  std::vector<double> ax(m_numElemConstraints*nspecies, 0.0);
69  std::vector<double> bb(m_numElemConstraints, 0.0);
70  std::vector<double> cc(nspecies, 0.0);
71 
72  int neActive = 0;
73  size_t jj = 0;
74  for (size_t j = 0; j < m_numElemConstraints; j++) {
75  if (m_elementActive[j]) {
76  neActive++;
77  bb[jj] = m_elemAbundancesGoal[j];
78  jj++;
79  }
80  }
81  for (kspec = 0; kspec < nspecies; ++kspec) {
82  cc[kspec] = -m_SSfeSpecies[kspec];
83  jj = 0;
84  for (size_t j = 0; j < m_numElemConstraints; ++j) {
85  if (m_elementActive[j]) {
86  ax[jj + kspec * neActive] = m_formulaMatrix[j][kspec];
87  jj++;
88  }
89  }
90  }
91  linprogmax(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(cc), VCS_DATA_PTR(ax),
92  VCS_DATA_PTR(bb), neActive, nspecies, neActive);
93 #endif
94 
95 #ifdef DEBUG_MODE
96  if (m_debug_print_lvl >= 2) {
97  plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
98  pprefix);
99  plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
100  for (kspec = 0; kspec < nspecies; ++kspec) {
101  plogf("%s ", pprefix);
102  plogf("%-12.12s", m_speciesName[kspec].c_str());
103  plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]);
104  }
105  plogf("%s Element Abundance Agreement returned from linear "
106  "programming (vcs_inest initial guess):",
107  pprefix);
108  plogendl();
109  plogf("%s Element Goal Actual\n", pprefix);
110  int jj = 0;
111  for (size_t j = 0; j < m_numElemConstraints; j++) {
112  if (m_elementActive[j]) {
113  double tmp = 0.0;
114  for (kspec = 0; kspec < nspecies; ++kspec) {
115  tmp += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec];
116  }
117  plogf("%s ", pprefix);
118  plogf(" %-9.9s", (m_elementName[j]).c_str());
119  plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp);
120  jj++;
121  }
122  }
123  plogendl();
124  }
125 #endif
126 
127 
128  /*
129  * Make sure all species have positive definite mole numbers
130  * Set voltages to zero for now, until we figure out what to do
131  */
132  vcs_dzero(VCS_DATA_PTR(m_deltaMolNumSpecies), nspecies);
133  for (kspec = 0; kspec < nspecies; ++kspec) {
134  iph = m_phaseID[kspec];
135  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
136  if (m_molNumSpecies_old[kspec] <= 0.0) {
137  /*
138  * HKM Should eventually include logic here for non SS phases
139  */
140  if (!m_SSPhase[kspec]) {
141  m_molNumSpecies_old[kspec] = 1.0e-30;
142  }
143  }
144  } else {
145  m_molNumSpecies_old[kspec] = 0.0;
146  }
147  }
148 
149  /*
150  * Now find the optimized basis that spans the stoichiometric
151  * coefficient matrix
152  */
153  (void) vcs_basopt(false, aw, sa, sm, ss, test, &conv);
154 
155  /* ***************************************************************** */
156  /* **** CALCULATE TOTAL MOLES, ****************** */
157  /* **** CHEMICAL POTENTIALS OF BASIS ****************** */
158  /* ***************************************************************** */
159  /*
160  * Calculate TMoles and m_tPhaseMoles_old[]
161  */
162  vcs_tmoles();
163  /*
164  * m_tPhaseMoles_new[] will consist of just the component moles
165  */
166  for (iph = 0; iph < m_numPhases; iph++) {
167  m_tPhaseMoles_new[iph] = TPhInertMoles[iph] + 1.0E-20;
168  }
169  for (kspec = 0; kspec < m_numComponents; ++kspec) {
170  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
171  m_tPhaseMoles_new[m_phaseID[kspec]] += m_molNumSpecies_old[kspec];
172  }
173  }
174  TMolesMultiphase = 0.0;
175  for (iph = 0; iph < m_numPhases; iph++) {
176  if (! m_VolPhaseList[iph]->m_singleSpecies) {
177  TMolesMultiphase += m_tPhaseMoles_new[iph];
178  }
179  }
180  vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
181  for (kspec = 0; kspec < m_numComponents; ++kspec) {
182  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_MOLNUM) {
183  m_molNumSpecies_new[kspec] = 0.0;
184  }
185  }
186  vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_SSfeSpecies),
187  nspecies);
188 
189  for (kspec = 0; kspec < m_numComponents; ++kspec) {
190  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
191  if (! m_SSPhase[kspec]) {
192  iph = m_phaseID[kspec];
193  m_feSpecies_new[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
194  }
195  } else {
196  m_molNumSpecies_new[kspec] = 0.0;
197  }
198  }
199  vcs_deltag(0, true, VCS_STATECALC_NEW);
200 #ifdef DEBUG_MODE
201  if (m_debug_print_lvl >= 2) {
202  for (kspec = 0; kspec < nspecies; ++kspec) {
203  plogf("%s", pprefix);
204  plogf("%-12.12s", m_speciesName[kspec].c_str());
205  if (kspec < m_numComponents)
206  plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_new[kspec],
207  m_SSfeSpecies[kspec]);
208  else
209  plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n",
210  m_feSpecies_new[kspec], m_SSfeSpecies[kspec], m_deltaGRxn_new[kspec-m_numComponents]);
211  }
212  }
213 #endif
214  /* ********************************************************** */
215  /* **** ESTIMATE REACTION ADJUSTMENTS *********************** */
216  /* ********************************************************** */
217  vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), m_numPhases);
218  for (iph = 0; iph < m_numPhases; iph++) {
219  xtphMax[iph] = log(m_tPhaseMoles_new[iph] * 1.0E32);
220  xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32);
221  }
222  for (irxn = 0; irxn < nrxn; ++irxn) {
223  kspec = m_indexRxnToSpecies[irxn];
224  /*
225  * For single species phases, we will not estimate the
226  * mole numbers. If the phase exists, it stays. If it
227  * doesn't exist in the estimate, it doesn't come into
228  * existence here.
229  */
230  if (! m_SSPhase[kspec]) {
231  iph = m_phaseID[kspec];
232  if (m_deltaGRxn_new[irxn] > xtphMax[iph]) {
233  m_deltaGRxn_new[irxn] = 0.8 * xtphMax[iph];
234  }
235  if (m_deltaGRxn_new[irxn] < xtphMin[iph]) {
236  m_deltaGRxn_new[irxn] = 0.8 * xtphMin[iph];
237  }
238  /*
239  * HKM -> The TMolesMultiphase is a change of mine.
240  * It more evenly distributes the initial moles amongst
241  * multiple multispecies phases according to the
242  * relative values of the standard state free energies.
243  * There is no change for problems with one multispecies
244  * phase.
245  * It cut diamond4.vin iterations down from 62 to 14.
246  */
247  m_deltaMolNumSpecies[kspec] = 0.5 * (m_tPhaseMoles_new[iph] + TMolesMultiphase)
248  * exp(-m_deltaGRxn_new[irxn]);
249 
250  for (size_t k = 0; k < m_numComponents; ++k) {
251  m_deltaMolNumSpecies[k] += m_stoichCoeffRxnMatrix[irxn][k] * m_deltaMolNumSpecies[kspec];
252  }
253 
254  for (iph = 0; iph < m_numPhases; iph++) {
255  m_deltaPhaseMoles[iph] += m_deltaMolNumPhase[irxn][iph] * m_deltaMolNumSpecies[kspec];
256  }
257  }
258  }
259 #ifdef DEBUG_MODE
260  if (m_debug_print_lvl >= 2) {
261  for (kspec = 0; kspec < nspecies; ++kspec) {
262  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
263  plogf("%sdirection (", pprefix);
264  plogf("%-12.12s", m_speciesName[kspec].c_str());
265  plogf(") = %g", m_deltaMolNumSpecies[kspec]);
266  if (m_SSPhase[kspec]) {
267  if (m_molNumSpecies_old[kspec] > 0.0) {
268  plogf(" (ssPhase exists at w = %g moles)", m_molNumSpecies_old[kspec]);
269  } else {
270  plogf(" (ssPhase doesn't exist -> stability not checked)");
271  }
272  }
273  plogendl();
274  }
275  }
276  }
277 #endif
278  /* *********************************************************** */
279  /* **** KEEP COMPONENT SPECIES POSITIVE ********************** */
280  /* *********************************************************** */
281  par = 0.5;
282  for (kspec = 0; kspec < m_numComponents; ++kspec) {
283  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
284  if (par < -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]) {
285  par = -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec];
286  }
287  }
288  }
289  par = 1. / par;
290  if (par <= 1.0 && par > 0.0) {
291  par *= 0.8;
292  } else {
293  par = 1.0;
294  }
295  /* ******************************************** */
296  /* **** CALCULATE NEW MOLE NUMBERS ************ */
297  /* ******************************************** */
298  finished = false;
299  do {
300  for (kspec = 0; kspec < m_numComponents; ++kspec) {
301  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
302  m_molNumSpecies_old[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec];
303  } else {
304  m_deltaMolNumSpecies[kspec] = 0.0;
305  }
306  }
307  for (kspec = m_numComponents; kspec < nspecies; ++kspec) {
308  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
309  if (m_deltaMolNumSpecies[kspec] != 0.0) {
310  m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par;
311  }
312  }
313  }
314  /*
315  * We have a new w[] estimate, go get the
316  * TMoles and m_tPhaseMoles_old[] values
317  */
318  vcs_tmoles();
319  if (lt > 0) {
320  goto finished;
321  }
322  /* ******************************************* */
323  /* **** CONVERGENCE FORCING SECTION ********** */
324  /* ******************************************* */
325  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
326  vcs_dfe(VCS_STATECALC_OLD, 0, 0, nspecies);
327  for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) {
328  s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec];
329  }
330  if (s == 0.0) {
331  finished = true;
332  continue;
333  }
334  if (s < 0.0) {
335  if (ikl == 0) {
336  finished = true;
337  continue;
338  }
339  }
340  /* ***************************************** */
341  /* *** TRY HALF STEP SIZE ****************** */
342  /* ***************************************** */
343  if (ikl == 0) {
344  s1 = s;
345  par *= 0.5;
346  ikl = 1;
347  continue;
348  }
349  /* **************************************************** */
350  /* **** FIT PARABOLA THROUGH HALF AND FULL STEPS ****** */
351  /* **************************************************** */
352  xl = (1.0 - s / (s1 - s)) * 0.5;
353  if (xl < 0.0) {
354  /* *************************************************** */
355  /* *** POOR DIRECTION, REDUCE STEP SIZE TO 0.2 ******* */
356  /* *************************************************** */
357  par *= 0.2;
358  } else {
359  if (xl > 1.0) {
360  /* *************************************************** */
361  /* **** TOO BIG A STEP, TAKE ORIGINAL FULL STEP ****** */
362  /* *************************************************** */
363  par *= 2.0;
364  } else {
365  /* *************************************************** */
366  /* **** ACCEPT RESULTS OF FORCER ********************* */
367  /* *************************************************** */
368  par = par * 2.0 * xl;
369  }
370  }
371  lt = 1;
372  } while (!finished);
373 finished:
374  ;
375 #ifdef DEBUG_MODE
376  if (m_debug_print_lvl >= 2) {
377  plogf("%s Final Mole Numbers produced by inest:\n",
378  pprefix);
379  plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
380  for (kspec = 0; kspec < nspecies; ++kspec) {
381  plogf("%s ", pprefix);
382  plogf("%-12.12s", m_speciesName[kspec].c_str());
383  plogf(" %g", m_molNumSpecies_old[kspec]);
384  plogendl();
385  }
386  }
387 #endif
388 }
389 /***************************************************************************/
390 
391 // Create an initial estimate of the solution to the thermodynamic
392 // equilibrium problem.
393 /*
394  * @return Return value indicates success:
395  * - 0: successful initial guess
396  * - -1: Unsuccessful initial guess; the elemental abundances aren't
397  * satisfied.
398  */
399 int VCS_SOLVE::vcs_inest_TP()
400 {
401  int retn = 0;
402  double test;
403  Cantera::clockWC tickTock;
404  test = -1.0E20;
405 
406  if (m_doEstimateEquil > 0) {
407  /*
408  * Calculate the elemental abundances
409  */
410  vcs_elab();
411  if (vcs_elabcheck(0)) {
412 #ifdef DEBUG_MODE
413  if (m_debug_print_lvl >= 2) {
414  plogf("%s Initial guess passed element abundances on input\n", pprefix);
415  plogf("%s m_doEstimateEquil = 1 so will use the input mole "
416  "numbers as estimates", pprefix);
417  plogendl();
418  }
419 #endif
420  return retn;
421 #ifdef DEBUG_MODE
422  } else {
423  if (m_debug_print_lvl >= 2) {
424  plogf("%s Initial guess failed element abundances on input\n", pprefix);
425  plogf("%s m_doEstimateEquil = 1 so will discard input "
426  "mole numbers and find our own estimate", pprefix);
427  plogendl();
428  }
429 #endif
430  }
431  }
432 
433  /*
434  * Malloc temporary space for usage in this routine and in
435  * subroutines
436  * sm[ne*ne]
437  * ss[ne]
438  * sa[ne]
439  * aw[m]
440  */
441  std::vector<double> sm(m_numElemConstraints*m_numElemConstraints, 0.0);
442  std::vector<double> ss(m_numElemConstraints, 0.0);
443  std::vector<double> sa(m_numElemConstraints, 0.0);
444  std::vector<double> aw(m_numSpeciesTot+ m_numElemConstraints, 0.0);
445  /*
446  * Go get the estimate of the solution
447  */
448 #ifdef DEBUG_MODE
449  if (m_debug_print_lvl >= 2) {
450  plogf("%sGo find an initial estimate for the equilibrium problem",
451  pprefix);
452  plogendl();
453  }
454 #endif
455  vcs_inest(VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm),
456  VCS_DATA_PTR(ss), test);
457  /*
458  * Calculate the elemental abundances
459  */
460  vcs_elab();
461 
462  /*
463  * If we still fail to achieve the correct elemental abundances,
464  * try to fix the problem again by calling the main elemental abundances
465  * fixer routine, used in the main program. This
466  * attempts to tweak the mole numbers of the component species to
467  * satisfy the element abundance constraints.
468  *
469  * Note: We won't do this unless we have to since it involves inverting
470  * a matrix.
471  */
472  bool rangeCheck = vcs_elabcheck(1);
473  if (!vcs_elabcheck(0)) {
474  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
475  plogf("%sInitial guess failed element abundances\n", pprefix);
476  plogf("%sCall vcs_elcorr to attempt fix", pprefix);
477  plogendl();
478  }
479  vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(aw));
480  rangeCheck = vcs_elabcheck(1);
481  if (!vcs_elabcheck(0)) {
482  plogf("%sInitial guess still fails element abundance equations\n",
483  pprefix);
484  plogf("%s - Inability to ever satisfy element abundance "
485  "constraints is probable", pprefix);
486  plogendl();
487  retn = -1;
488  } else {
489  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
490  if (rangeCheck) {
491  plogf("%sInitial guess now satisfies element abundances", pprefix);
492  plogendl();
493  } else {
494  plogf("%sElement Abundances RANGE ERROR\n", pprefix);
495  plogf("%s - Initial guess satisfies NC=%d element abundances, "
496  "BUT not NE=%d element abundances", pprefix,
497  m_numComponents, m_numElemConstraints);
498  plogendl();
499  }
500  }
501  }
502  } else {
503  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
504  if (rangeCheck) {
505  plogf("%sInitial guess satisfies element abundances", pprefix);
506  plogendl();
507  } else {
508  plogf("%sElement Abundances RANGE ERROR\n", pprefix);
509  plogf("%s - Initial guess satisfies NC=%d element abundances, "
510  "BUT not NE=%d element abundances", pprefix,
511  m_numComponents, m_numElemConstraints);
512  plogendl();
513  }
514  }
515  }
516 
517 #ifdef DEBUG_MODE
518  if (m_debug_print_lvl >= 2) {
519  plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
520  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_new),
521  VCS_DATA_PTR(m_tPhaseMoles_old)));
522  plogendl();
523  }
524 #endif
525 
526  /*
527  * Record time
528  */
529  double tsecond = tickTock.secondsWC();
530  m_VCount->T_Time_inest += tsecond;
531  (m_VCount->T_Calls_Inest)++;
532  return retn;
533 }
534 
535 }
Generated by   doxygen 1.8.2