Cantera  2.5.1
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 // This file is part of Cantera. See License.txt in the top-level directory or
7 // at https://cantera.org/license.txt for license and copyright information.
8 
11 
12 #include "cantera/base/clockWC.h"
13 
14 namespace Cantera
15 {
16 
17 static char pprefix[20] = " --- vcs_inest: ";
18 
19 void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
20  double* const ss, double test)
21 {
22  size_t nrxn = m_numRxnTot;
23 
24  // CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB AND
25  // molNum(I) .GE. 0.0. Note, both of these programs do this.
27 
28  if (m_debug_print_lvl >= 2) {
29  plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
30  pprefix);
31  plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
32  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
33  plogf("%s ", pprefix);
34  plogf("%-12.12s", m_speciesName[kspec]);
35  plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]);
36  }
37  plogf("%s Element Abundance Agreement returned from linear "
38  "programming (vcs_inest initial guess):\n", pprefix);
39  plogf("%s Element Goal Actual\n", pprefix);
40  for (size_t j = 0; j < m_nelem; j++) {
41  if (m_elementActive[j]) {
42  double tmp = 0.0;
43  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
44  tmp += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec];
45  }
46  plogf("%s ", pprefix);
47  plogf(" %-9.9s", m_elementName[j]);
48  plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp);
49  }
50  }
51  writelogendl();
52  }
53 
54  // Make sure all species have positive definite mole numbers Set voltages to
55  // zero for now, until we figure out what to do
56  m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0);
57  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
59  if (m_molNumSpecies_old[kspec] <= 0.0) {
60  // HKM Should eventually include logic here for non SS phases
61  if (!m_SSPhase[kspec]) {
62  m_molNumSpecies_old[kspec] = 1.0e-30;
63  }
64  }
65  } else {
66  m_molNumSpecies_old[kspec] = 0.0;
67  }
68  }
69 
70  // Now find the optimized basis that spans the stoichiometric coefficient
71  // matrix
72  bool conv;
73  vcs_basopt(false, aw, sa, sm, ss, test, &conv);
74 
75  // CALCULATE TOTAL MOLES, CHEMICAL POTENTIALS OF BASIS
76 
77  // Calculate TMoles and m_tPhaseMoles_old[]
78  vcs_tmoles();
79 
80  // m_tPhaseMoles_new[] will consist of just the component moles
81  for (size_t iph = 0; iph < m_numPhases; iph++) {
82  m_tPhaseMoles_new[iph] = TPhInertMoles[iph] + 1.0E-20;
83  }
84  for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
87  }
88  }
89  double TMolesMultiphase = 0.0;
90  for (size_t iph = 0; iph < m_numPhases; iph++) {
91  if (! m_VolPhaseList[iph]->m_singleSpecies) {
92  TMolesMultiphase += m_tPhaseMoles_new[iph];
93  }
94  }
96  for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
98  m_molNumSpecies_new[kspec] = 0.0;
99  }
100  }
102 
103  for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
105  if (! m_SSPhase[kspec]) {
106  size_t iph = m_phaseID[kspec];
107  m_feSpecies_new[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
108  }
109  } else {
110  m_molNumSpecies_new[kspec] = 0.0;
111  }
112  }
113  vcs_deltag(0, true, VCS_STATECALC_NEW);
114  if (m_debug_print_lvl >= 2) {
115  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
116  plogf("%s", pprefix);
117  plogf("%-12.12s", m_speciesName[kspec]);
118  if (kspec < m_numComponents) {
119  plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_new[kspec],
120  m_SSfeSpecies[kspec]);
121  } else {
122  plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n",
124  }
125  }
126  }
127 
128  // ESTIMATE REACTION ADJUSTMENTS
129  vector_fp& xtphMax = m_TmpPhase;
130  vector_fp& xtphMin = m_TmpPhase2;
131  m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
132  for (size_t iph = 0; iph < m_numPhases; iph++) {
133  xtphMax[iph] = log(m_tPhaseMoles_new[iph] * 1.0E32);
134  xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32);
135  }
136  for (size_t irxn = 0; irxn < nrxn; ++irxn) {
137  size_t kspec = m_indexRxnToSpecies[irxn];
138 
139  // For single species phases, we will not estimate the mole numbers. If
140  // the phase exists, it stays. If it doesn't exist in the estimate, it
141  // doesn't come into existence here.
142  if (! m_SSPhase[kspec]) {
143  size_t iph = m_phaseID[kspec];
144  if (m_deltaGRxn_new[irxn] > xtphMax[iph]) {
145  m_deltaGRxn_new[irxn] = 0.8 * xtphMax[iph];
146  }
147  if (m_deltaGRxn_new[irxn] < xtphMin[iph]) {
148  m_deltaGRxn_new[irxn] = 0.8 * xtphMin[iph];
149  }
150 
151  // HKM -> The TMolesMultiphase is a change of mine. It more evenly
152  // distributes the initial moles amongst multiple multispecies
153  // phases according to the relative values of the standard state
154  // free energies. There is no change for problems with one
155  // multispecies phase. It cut diamond4.vin iterations down from 62
156  // to 14.
157  m_deltaMolNumSpecies[kspec] = 0.5 * (m_tPhaseMoles_new[iph] + TMolesMultiphase)
158  * exp(-m_deltaGRxn_new[irxn]);
159 
160  for (size_t k = 0; k < m_numComponents; ++k) {
162  }
163 
164  for (iph = 0; iph < m_numPhases; iph++) {
165  m_deltaPhaseMoles[iph] += m_deltaMolNumPhase(iph,irxn) * m_deltaMolNumSpecies[kspec];
166  }
167  }
168  }
169  if (m_debug_print_lvl >= 2) {
170  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
172  plogf("%sdirection (", pprefix);
173  plogf("%-12.12s", m_speciesName[kspec]);
174  plogf(") = %g", m_deltaMolNumSpecies[kspec]);
175  if (m_SSPhase[kspec]) {
176  if (m_molNumSpecies_old[kspec] > 0.0) {
177  plogf(" (ssPhase exists at w = %g moles)", m_molNumSpecies_old[kspec]);
178  } else {
179  plogf(" (ssPhase doesn't exist -> stability not checked)");
180  }
181  }
182  writelogendl();
183  }
184  }
185  }
186 
187  // KEEP COMPONENT SPECIES POSITIVE
188  double par = 0.5;
189  for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
191  par < -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]) {
192  par = -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec];
193  }
194  }
195  par = 1. / par;
196  if (par <= 1.0 && par > 0.0) {
197  par *= 0.8;
198  } else {
199  par = 1.0;
200  }
201 
202  // CALCULATE NEW MOLE NUMBERS
203  size_t lt = 0;
204  size_t ikl = 0;
205  double s1 = 0.0;
206  while (true) {
207  for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
209  m_molNumSpecies_old[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec];
210  } else {
211  m_deltaMolNumSpecies[kspec] = 0.0;
212  }
213  }
214  for (size_t kspec = m_numComponents; kspec < m_nsp; ++kspec) {
216  m_deltaMolNumSpecies[kspec] != 0.0) {
217  m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par;
218  }
219  }
220 
221  // We have a new w[] estimate, go get the TMoles and m_tPhaseMoles_old[]
222  // values
223  vcs_tmoles();
224  if (lt > 0) {
225  break;
226  }
227 
228  // CONVERGENCE FORCING SECTION
229  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
231  double s = 0.0;
232  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
233  s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec];
234  }
235  if (s == 0.0) {
236  break;
237  }
238  if (s < 0.0 && ikl == 0) {
239  break;
240  }
241 
242  // TRY HALF STEP SIZE
243  if (ikl == 0) {
244  s1 = s;
245  par *= 0.5;
246  ikl = 1;
247  continue;
248  }
249 
250  // FIT PARABOLA THROUGH HALF AND FULL STEPS
251  double xl = (1.0 - s / (s1 - s)) * 0.5;
252  if (xl < 0.0) {
253  // POOR DIRECTION, REDUCE STEP SIZE TO 0.2
254  par *= 0.2;
255  } else {
256  if (xl > 1.0) {
257  // TOO BIG A STEP, TAKE ORIGINAL FULL STEP
258  par *= 2.0;
259  } else {
260  // ACCEPT RESULTS OF FORCER
261  par = par * 2.0 * xl;
262  }
263  }
264  lt = 1;
265  }
266 
267  if (m_debug_print_lvl >= 2) {
268  plogf("%s Final Mole Numbers produced by inest:\n",
269  pprefix);
270  plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
271  for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
272  plogf("%s %-12.12s %g\n",
273  pprefix, m_speciesName[kspec], m_molNumSpecies_old[kspec]);
274  }
275  }
276 }
277 
279 {
280  int retn = 0;
281  clockWC tickTock;
282  if (m_doEstimateEquil > 0) {
283  // Calculate the elemental abundances
284  vcs_elab();
285  if (vcs_elabcheck(0)) {
286  if (m_debug_print_lvl >= 2) {
287  plogf("%s Initial guess passed element abundances on input\n", pprefix);
288  plogf("%s m_doEstimateEquil = 1 so will use the input mole "
289  "numbers as estimates\n", pprefix);
290  }
291  return retn;
292  } else if (m_debug_print_lvl >= 2) {
293  plogf("%s Initial guess failed element abundances on input\n", pprefix);
294  plogf("%s m_doEstimateEquil = 1 so will discard input "
295  "mole numbers and find our own estimate\n", pprefix);
296  }
297  }
298 
299  // temporary space for usage in this routine and in subroutines
300  vector_fp sm(m_nelem*m_nelem, 0.0);
301  vector_fp ss(m_nelem, 0.0);
302  vector_fp sa(m_nelem, 0.0);
303  vector_fp aw(m_nsp + m_nelem, 0.0);
304 
305  // Go get the estimate of the solution
306  if (m_debug_print_lvl >= 2) {
307  plogf("%sGo find an initial estimate for the equilibrium problem\n",
308  pprefix);
309  }
310  double test = -1.0E20;
311  vcs_inest(&aw[0], &sa[0], &sm[0], &ss[0], test);
312 
313  // Calculate the elemental abundances
314  vcs_elab();
315 
316  // If we still fail to achieve the correct elemental abundances, try to fix
317  // the problem again by calling the main elemental abundances fixer routine,
318  // used in the main program. This attempts to tweak the mole numbers of the
319  // component species to satisfy the element abundance constraints.
320  //
321  // Note: We won't do this unless we have to since it involves inverting a
322  // matrix.
323  bool rangeCheck = vcs_elabcheck(1);
324  if (!vcs_elabcheck(0)) {
325  if (m_debug_print_lvl >= 2) {
326  plogf("%sInitial guess failed element abundances\n", pprefix);
327  plogf("%sCall vcs_elcorr to attempt fix\n", pprefix);
328  }
329  vcs_elcorr(&sm[0], &aw[0]);
330  rangeCheck = vcs_elabcheck(1);
331  if (!vcs_elabcheck(0)) {
332  plogf("%sInitial guess still fails element abundance equations\n",
333  pprefix);
334  plogf("%s - Inability to ever satisfy element abundance "
335  "constraints is probable\n", pprefix);
336  retn = -1;
337  } else {
338  if (m_debug_print_lvl >= 2) {
339  if (rangeCheck) {
340  plogf("%sInitial guess now satisfies element abundances\n", pprefix);
341  } else {
342  plogf("%sElement Abundances RANGE ERROR\n", pprefix);
343  plogf("%s - Initial guess satisfies NC=%d element abundances, "
344  "BUT not NE=%d element abundances\n", pprefix,
346  }
347  }
348  }
349  } else {
350  if (m_debug_print_lvl >= 2) {
351  if (rangeCheck) {
352  plogf("%sInitial guess satisfies element abundances\n", pprefix);
353  } else {
354  plogf("%sElement Abundances RANGE ERROR\n", pprefix);
355  plogf("%s - Initial guess satisfies NC=%d element abundances, "
356  "BUT not NE=%d element abundances\n", pprefix,
358  }
359  }
360  }
361 
362  if (m_debug_print_lvl >= 2) {
363  plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E\n", pprefix,
365  &m_tPhaseMoles_old[0]));
366  }
367 
368  // Record time
369  m_VCount->T_Time_inest += tickTock.secondsWC();
371  return retn;
372 }
373 
374 }
int T_Calls_Inest
Current number of times the initial thermo equilibrium estimator has been called.
Definition: vcs_internal.h:54
double T_Time_inest
Time spent in initial estimator.
Definition: vcs_internal.h:72
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1068
int m_doEstimateEquil
Setting for whether to do an initial estimate.
Definition: vcs_solve.h:1145
int vcs_inest_TP()
Create an initial estimate of the solution to the thermodynamic equilibrium problem.
Definition: vcs_inest.cpp:278
void vcs_deltag(const int L, const bool doDeleted, const int vcsState, const bool alterZeroedPhases=true)
This subroutine calculates reaction free energy changes for all noncomponent formation reactions.
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1364
vector_fp m_TmpPhase2
Temporary vector of length NPhase.
Definition: vcs_solve.h:1261
size_t m_nelem
Number of element constraints in the problem.
Definition: vcs_solve.h:1047
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.
int vcs_elcorr(double aa[], double x[])
Definition: vcs_elem.cpp:97
vector_fp m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1152
vector_fp m_molNumSpecies_new
Tentative value of the mole number vector.
Definition: vcs_solve.h:1183
void vcs_inest(double *const aw, double *const sa, double *const sm, double *const ss, double test)
Estimate equilibrium compositions.
Definition: vcs_inest.cpp:19
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1165
void vcs_elab()
Computes the current elemental abundances vector.
Definition: vcs_elem.cpp:17
int vcs_setMolesLinProg()
Estimate the initial mole numbers by constrained linear programming.
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles)
Calculate the total dimensionless Gibbs free energy.
Definition: vcs_Gibbs.cpp:14
vector_fp m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
Definition: vcs_solve.h:1255
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:1135
vector_fp m_tPhaseMoles_old
Total kmols of species in each phase.
Definition: vcs_solve.h:1246
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1498
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_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1076
vector_fp m_deltaPhaseMoles
Change in the total moles in each phase.
Definition: vcs_solve.h:1267
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
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1050
vector_fp m_feSpecies_old
Free energy vector from the start of the current iteration.
Definition: vcs_solve.h:1126
vector_int m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1391
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1361
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
bool vcs_elabcheck(int ibound)
Definition: vcs_elem.cpp:29
vector_fp m_deltaMolNumSpecies
Reaction Adjustments for each species during the current step.
Definition: vcs_solve.h:1213
vector_fp m_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
Definition: vcs_solve.h:1119
Array2D m_deltaMolNumPhase
Change in the number of moles of phase, iphase, due to the noncomponent formation reaction,...
Definition: vcs_solve.h:1172
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
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,...
size_t m_numRxnTot
Total number of non-component species in the problem.
Definition: vcs_solve.h:1053
vector_fp m_TmpPhase
Temporary vector of length NPhase.
Definition: vcs_solve.h:1258
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1370
The class provides the wall clock timer in seconds.
Definition: clockWC.h:45
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition: clockWC.cpp:32
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
void writelogendl()
Write an end of line character to the screen and flush output.
Definition: global.cpp:38
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
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
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_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
Definition: vcs_defs.h:303
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:281
#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...