Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcs_phaseStability.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_phaseStability.cpp
3  * Implementation class for functions associated with determining the stability of a phase
4  * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ).
5  */
10 
11 using namespace std;
12 
13 namespace Cantera
14 {
15 
16 bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
17 {
18  vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
19 
20  AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhasePossible",
21  "called for a phase that exists!");
22 
23  /*
24  * Loop through all of the species in the phase. We say the phase
25  * can be popped, if there is one species in the phase that can be
26  * popped. This does not mean that the phase will be popped or that it
27  * leads to a lower Gibbs free energy.
28  */
29  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
30  size_t kspec = Vphase->spGlobalIndexVCS(k);
31  AssertThrowMsg(m_molNumSpecies_old[kspec] <= 0.0,
32  "VCS_SOLVE::vcs_popPhasePossible",
33  "we shouldn't be here " + int2str(kspec) + " "+
34  fp2str(m_molNumSpecies_old[kspec]) + " > 0.0");
35  size_t irxn = kspec - m_numComponents;
36  if (kspec >= m_numComponents) {
37  bool iPopPossible = true;
38  /*
39  * Note one case is if the component is a member of the popping phase.
40  * This component will be zeroed and the logic here will negate the current
41  * species from causing a positive if this component is consumed.
42  */
43  for (size_t j = 0; j < m_numComponents; ++j) {
44  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
45  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
46  if (stoicC != 0.0) {
47  double negChangeComp = - stoicC;
48  if (negChangeComp > 0.0) {
49  // If there is no component to give, then the species can't be created
50  if (m_molNumSpecies_old[j] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
51  iPopPossible = false;
52  }
53  }
54  }
55  }
56  }
57  // We are here when the species can be popped because all its needed components have positive mole numbers
58  if (iPopPossible) {
59  return true;
60  }
61  } else {
62  /*
63  * We are here when the species, k, in the phase is a component. Its mole number is zero.
64  * We loop through the regular reaction looking for a reaction that can pop the
65  * component.
66  */
67  for (size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) {
68  bool foundJrxn = false;
69  // First, if the component is a product of the reaction
70  if (m_stoichCoeffRxnMatrix(kspec,jrxn) > 0.0) {
71  foundJrxn = true;
72  // We can do the reaction if all other reactant components have positive mole fractions
73  for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
74  if (m_stoichCoeffRxnMatrix(kcomp,jrxn) < 0.0) {
75  if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
76  foundJrxn = false;
77  }
78  }
79  }
80  if (foundJrxn) {
81  return true;
82  }
83  }
84  // Second we are here if the component is a reactant in the reaction, and the reaction goes backwards.
85  else if (m_stoichCoeffRxnMatrix(kspec,jrxn) < 0.0) {
86  foundJrxn = true;
87  size_t jspec = jrxn + m_numComponents;
88  if (m_molNumSpecies_old[jspec] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
89  foundJrxn = false;
90  continue;
91  }
92  // We can do the backwards reaction if all of the product components species are positive
93  for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
94  if (m_stoichCoeffRxnMatrix(kcomp,jrxn) > 0.0) {
95  if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
96  foundJrxn = false;
97  }
98  }
99  }
100  if (foundJrxn) {
101  return true;
102  }
103  }
104  }
105  }
106  }
107  return false;
108 }
109 
110 int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
111 {
112  int nfound = 0;
113  phasePopProblemLists_.clear();
114 
115  /*
116  * This is a vector over each component.
117  * For zeroed components it lists the phases, which are currently zeroed,
118  * which have a species with a positive stoichiometric value wrt the component.
119  * Therefore, we could pop the component species and pop that phase at the same time
120  * if we considered no other factors than keeping the component mole number positive.
121  *
122  * It does not count species with positive stoichiometric values if that species
123  * already has a positive mole number. The phase is already popped.
124  */
125  std::vector< std::vector<size_t> > zeroedComponentLinkedPhasePops(m_numComponents);
126  /*
127  * The logic below calculates zeroedComponentLinkedPhasePops
128  */
129  for (size_t j = 0; j < m_numComponents; j++) {
130  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
131  if (m_molNumSpecies_old[j] <= 0.0) {
132  std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
133  size_t iph = m_phaseID[j];
134  jList.push_back(iph);
135  for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
136  size_t kspec = irxn + m_numComponents;
137  iph = m_phaseID[kspec];
138  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
139  int existence = Vphase->exists();
140  if (existence < 0) {
141  if (m_stoichCoeffRxnMatrix(j,irxn) > 0.0) {
142  if (std::find(jList.begin(), jList.end(), iph) != jList.end()) {
143  jList.push_back(iph);
144  }
145  }
146  }
147  }
148  }
149  }
150  }
151  /*
152  * This is a vector over each zeroed phase
153  * For zeroed phases, it lists the components, which are currently zeroed,
154  * which have a species with a negative stoichiometric value wrt one or more species in the phase.
155  * Cut out components which have a pos stoichiometric value with another species in the phase.
156  */
157  std::vector< std::vector<size_t> > zeroedPhaseLinkedZeroComponents(m_numPhases);
158  std::vector<int> linkedPhases;
159  /*
160  * The logic below calculates zeroedPhaseLinkedZeroComponents
161  */
162  for (size_t iph = 0; iph < m_numPhases; iph++) {
163  std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
164  iphList.clear();
165  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
166  if (Vphase->exists() < 0) {
167 
168  linkedPhases.clear();
169  size_t nsp = Vphase->nSpecies();
170  for (size_t k = 0; k < nsp; k++) {
171  size_t kspec = Vphase->spGlobalIndexVCS(k);
172  size_t irxn = kspec - m_numComponents;
173 
174  for (size_t j = 0; j < m_numComponents; j++) {
175  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
176  if (m_molNumSpecies_old[j] <= 0.0) {
177  if (m_stoichCoeffRxnMatrix(j,irxn) < 0.0) {
178  bool foundPos = false;
179  for (size_t kk = 0; kk < nsp; kk++) {
180  size_t kkspec = Vphase->spGlobalIndexVCS(kk);
181  if (kkspec >= m_numComponents) {
182  size_t iirxn = kkspec - m_numComponents;
183  if (m_stoichCoeffRxnMatrix(j,iirxn) > 0.0) {
184  foundPos = true;
185  }
186  }
187  }
188  if (!foundPos) {
189  if (std::find(iphList.begin(), iphList.end(), j) != iphList.end()) {
190  iphList.push_back(j);
191  }
192  }
193  }
194  }
195  }
196  }
197  }
198  }
199  }
200 
201  /*
202  * Now fill in the phasePopProblemLists_ list.
203  *
204  */
205  for (size_t iph = 0; iph < m_numPhases; iph++) {
206  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
207  if (Vphase->exists() < 0) {
208  std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
209  std::vector<size_t> popProblem(0);
210  popProblem.push_back(iph);
211  for (size_t i = 0; i < iphList.size(); i++) {
212  size_t j = iphList[i];
213  std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
214  for (size_t jjl = 0; jjl < jList.size(); jjl++) {
215  size_t jph = jList[jjl];
216  if (std::find(popProblem.begin(), popProblem.end(), jph) != popProblem.end()) {
217  popProblem.push_back(jph);
218  }
219  }
220  }
221  phasePopProblemLists_.push_back(popProblem);
222  }
223  }
224 
225  return nfound;
226 }
227 
228 size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
229 {
230  size_t iphasePop = npos;
231  doublereal FephaseMax = -1.0E30;
232  doublereal Fephase = -1.0E30;
233 
234 #ifdef DEBUG_MODE
235  char anote[128];
236  if (m_debug_print_lvl >= 2) {
237  plogf(" --- vcs_popPhaseID() called\n");
238  plogf(" --- Phase Status F_e MoleNum\n");
239  plogf(" --------------------------------------------------------------------------\n");
240  }
241 #else
242  char* anote;
243 #endif
244  for (size_t iph = 0; iph < m_numPhases; iph++) {
245  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
246  int existence = Vphase->exists();
247  if (DEBUG_MODE_ENABLED) {
248  strcpy(anote, "");
249  }
250  if (existence > 0) {
251  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
252  plogf(" --- %18s %5d NA %11.3e\n",
253  Vphase->PhaseName.c_str(),
254  existence,
255  m_tPhaseMoles_old[iph]);
256  }
257  } else {
258  if (Vphase->m_singleSpecies) {
259  /***********************************************************************
260  *
261  * Single Phase Stability Resolution
262  *
263  ***********************************************************************/
264  size_t kspec = Vphase->spGlobalIndexVCS(0);
265  size_t irxn = kspec - m_numComponents;
266  doublereal deltaGRxn = m_deltaGRxn_old[irxn];
267  Fephase = exp(-deltaGRxn) - 1.0;
268  if (Fephase > 0.0) {
269  if (DEBUG_MODE_ENABLED) {
270  strcpy(anote," (ready to be birthed)");
271  }
272  if (Fephase > FephaseMax) {
273  iphasePop = iph;
274  FephaseMax = Fephase;
275  if (DEBUG_MODE_ENABLED) {
276  strcpy(anote," (chosen to be birthed)");
277  }
278  }
279  }
280  if (DEBUG_MODE_ENABLED && Fephase < 0.0) {
281  strcpy(anote," (not stable)");
282  AssertThrowMsg(m_tPhaseMoles_old[iph] <= 0.0,
283  "VCS_SOLVE::vcs_popPhaseID", "shouldn't be here");
284  }
285 
286  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
287  plogf(" --- %18s %5d %10.3g %10.3g %s\n",
288  Vphase->PhaseName.c_str(),
289  existence, Fephase,
290  m_tPhaseMoles_old[iph], anote);
291  }
292 
293  } else {
294  /***********************************************************************
295  *
296  * MultiSpecies Phase Stability Resolution
297  *
298  ***********************************************************************/
299  if (vcs_popPhasePossible(iph)) {
300  Fephase = vcs_phaseStabilityTest(iph);
301  if (Fephase > 0.0) {
302  if (Fephase > FephaseMax) {
303  iphasePop = iph;
304  FephaseMax = Fephase;
305  }
306  } else {
307  FephaseMax = std::max(FephaseMax, Fephase);
308  }
309  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
310  plogf(" --- %18s %5d %11.3g %11.3g\n",
311  Vphase->PhaseName.c_str(),
312  existence, Fephase,
313  m_tPhaseMoles_old[iph]);
314  }
315  } else {
316  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
317  plogf(" --- %18s %5d blocked %11.3g\n",
318  Vphase->PhaseName.c_str(),
319  existence, m_tPhaseMoles_old[iph]);
320  }
321  }
322  }
323  }
324  }
325  phasePopPhaseIDs.resize(0);
326  if (iphasePop != npos) {
327  phasePopPhaseIDs.push_back(iphasePop);
328  }
329 
330  /*
331  * Insert logic here to figure out if phase pops are linked together. Only do one linked
332  * pop at a time.
333  */
334 
335  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
336  plogf(" ---------------------------------------------------------------------\n");
337  }
338  return iphasePop;
339 }
340 
341 int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
342 {
343  vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
344  // Identify the first species in the phase
345  size_t kspec = Vphase->spGlobalIndexVCS(0);
346  // Identify the formation reaction for that species
347  size_t irxn = kspec - m_numComponents;
348  std::vector<size_t> creationGlobalRxnNumbers;
349 
350  // Calculate the initial moles of the phase being born.
351  // Here we set it to 10x of the value which would cause the phase to be
352  // zeroed out within the algorithm. We may later adjust the value.
353  doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
354 
355  AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhaseRxnStepSizes",
356  "called for a phase that exists!");
357  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
358  plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
359  Vphase->PhaseName.c_str(), iphasePop);
360  }
361  // Section for a single-species phase
362  if (Vphase->m_singleSpecies) {
363  double s = 0.0;
364  for (size_t j = 0; j < m_numComponents; ++j) {
365  if (!m_SSPhase[j]) {
366  if (m_molNumSpecies_old[j] > 0.0) {
367  s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j];
368  }
369  }
370  }
371  for (size_t j = 0; j < m_numPhases; j++) {
372  Vphase = m_VolPhaseList[j];
373  if (! Vphase->m_singleSpecies) {
374  if (m_tPhaseMoles_old[j] > 0.0) {
375  s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j];
376  }
377  }
378  }
379  if (s != 0.0) {
380  double s_old = s;
381  s = vcs_Hessian_diag_adj(irxn, s_old);
382  m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
383  } else {
384  // Ok, s is equal to zero. We can not apply a sophisticated theory
385  // to birth the phase. Just pick a small delta and go with it.
386  m_deltaMolNumSpecies[kspec] = tPhaseMoles;
387  }
388 
389  /*
390  * section to do damping of the m_deltaMolNumSpecies[]
391  */
392  for (size_t j = 0; j < m_numComponents; ++j) {
393  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
394  if (stoicC != 0.0) {
395  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
396  double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
397  if (negChangeComp > m_molNumSpecies_old[j]) {
398  if (m_molNumSpecies_old[j] > 0.0) {
399  m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
400  } else {
401  m_deltaMolNumSpecies[kspec] = 0.0;
402  }
403  }
404  }
405  }
406  }
407  // Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number
408  if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
409  m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
410  }
411 
412 
413  } else {
414  vector<doublereal> fracDelta(Vphase->nSpecies());
415  vector<doublereal> X_est(Vphase->nSpecies());
416  fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
417 
418  double sumFrac = 0.0;
419  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
420  sumFrac += fracDelta[k];
421  }
422  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
423  X_est[k] = fracDelta[k] / sumFrac;
424  }
425 
426  doublereal deltaMolNumPhase = tPhaseMoles;
427  doublereal damp = 1.0;
428  m_deltaGRxn_tmp = m_molNumSpecies_old;
429  double* molNumSpecies_tmp = DATA_PTR(m_deltaGRxn_tmp);
430 
431 
432  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
433  kspec = Vphase->spGlobalIndexVCS(k);
434  double delmol = deltaMolNumPhase * X_est[k];
435  if (kspec >= m_numComponents) {
436  irxn = kspec - m_numComponents;
437  for (size_t j = 0; j < m_numComponents; ++j) {
438  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
439  if (stoicC != 0.0) {
440  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
441  molNumSpecies_tmp[j] += stoicC * delmol;
442  }
443  }
444  }
445  }
446  }
447 
448  doublereal ratioComp = 0.0;
449  for (size_t j = 0; j < m_numComponents; ++j) {
450  double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j];
451  if (molNumSpecies_tmp[j] < 0.0) {
452  ratioComp = 1.0;
453  if (deltaJ > 0.0) {
454  double delta0 = m_molNumSpecies_old[j];
455  damp = std::min(damp, delta0 / deltaJ * 0.9);
456  }
457  } else {
458  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
459  size_t jph = m_phaseID[j];
460  if ((jph != iphasePop) && (!m_SSPhase[j])) {
461  double fdeltaJ = fabs(deltaJ);
462  if (m_molNumSpecies_old[j] > 0.0) {
463  ratioComp = std::max(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]);
464  }
465  }
466  }
467  }
468  }
469 
470  // We may have greatly underestimated the deltaMoles for the phase pop
471  // Here we create a damp > 1 to account for this possibility.
472  // We adjust upwards to make sure that a component in an existing multispecies
473  // phase is modified by a factor of 1/1000.
474  if (ratioComp > 1.0E-30) {
475  if (ratioComp < 0.001) {
476  damp = 0.001 / ratioComp;
477  }
478  }
479 
480 
481  if (damp <= 1.0E-6) {
482  return 3;
483  }
484 
485  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
486  kspec = Vphase->spGlobalIndexVCS(k);
487  if (kspec < m_numComponents) {
488  m_speciesStatus[kspec] = VCS_SPECIES_COMPONENT;
489  } else {
490  m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp;
491  if (X_est[k] > 1.0E-3) {
492  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
493  } else {
494  m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
495  }
496  }
497  }
498 
499  }
500 
501  return 0;
502 }
503 
504 double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
505 {
506  /*
507  * We will use the _new state calc here
508  */
509  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
510  const size_t nsp = Vphase->nSpecies();
511  int minNumberIterations = 3;
512  if (nsp <= 1) {
513  minNumberIterations = 1;
514  }
515 
516  // We will do a full Newton calculation later, but for now, ...
517  bool doSuccessiveSubstitution = true;
518  double funcPhaseStability;
519  vector<doublereal> X_est(nsp, 0.0);
520  vector<doublereal> delFrac(nsp, 0.0);
521  vector<doublereal> E_phi(nsp, 0.0);
522  vector<doublereal> fracDelta_new(nsp, 0.0);
523  vector<doublereal> fracDelta_old(nsp, 0.0);
524  vector<doublereal> fracDelta_raw(nsp, 0.0);
525  vector<size_t> creationGlobalRxnNumbers(nsp, npos);
526  m_deltaGRxn_Deficient = m_deltaGRxn_old;
527 
528  vector<doublereal> m_feSpecies_Deficient(m_numComponents, 0.0);
529  doublereal damp = 1.0;
530  doublereal dampOld = 1.0;
531  doublereal normUpdate = 1.0;
532  doublereal normUpdateOld = 1.0;
533  doublereal sum = 0.0;
534  doublereal dirProd = 0.0;
535  doublereal dirProdOld = 0.0;
536 
537  // get the activity coefficients
538  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
539 
540  // Get the stored estimate for the composition of the phase if
541  // it gets created
542  fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
543 
544 
545  std::vector<size_t> componentList;
546 
547  for (size_t k = 0; k < nsp; k++) {
548  size_t kspec = Vphase->spGlobalIndexVCS(k);
549  if (kspec < m_numComponents) {
550  componentList.push_back(k);
551  }
552  }
553 
554  for (size_t k = 0; k < m_numComponents; k++) {
555  m_feSpecies_Deficient[k] = m_feSpecies_old[k];
556  }
557  normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
558  damp = 1.0E-2;
559 
560  if (doSuccessiveSubstitution) {
561  int KP = 0;
562  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
563  plogf(" --- vcs_phaseStabilityTest() called\n");
564  plogf(" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
565  " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
566  plogf(" --------------------------------------------------------------"
567  "--------------------------------------------------------\n");
568  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) {
569  plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
570  }
571 
572  for (size_t k = 0; k < nsp; k++) {
573  if (fracDelta_new[k] < 1.0E-13) {
574  fracDelta_new[k] =1.0E-13;
575  }
576  }
577  bool converged = false;
578  for (int its = 0; its < 200 && (!converged); its++) {
579 
580  dampOld = damp;
581  normUpdateOld = normUpdate;
582  fracDelta_old = fracDelta_new;
583  dirProdOld = dirProd;
584 
585  // Given a set of fracDelta's, we calculate the fracDelta's
586  // for the component species, if any
587  for (size_t i = 0; i < componentList.size(); i++) {
588  size_t kc = componentList[i];
589  size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
590  fracDelta_old[kc] = 0.0;
591  for (size_t k = 0; k < nsp; k++) {
592  size_t kspec = Vphase->spGlobalIndexVCS(k);
593  if (kspec >= m_numComponents) {
594  size_t irxn = kspec - m_numComponents;
595  fracDelta_old[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_old[k];
596  }
597  }
598  }
599 
600  // Now, calculate the predicted mole fractions, X_est[k]
601  double sumFrac = 0.0;
602  for (size_t k = 0; k < nsp; k++) {
603  sumFrac += fracDelta_old[k];
604  }
605  // Necessary because this can be identically zero. -> we need to fix this algorithm!
606  if (sumFrac <= 0.0) {
607  sumFrac = 1.0;
608  }
609  double sum_Xcomp = 0.0;
610  for (size_t k = 0; k < nsp; k++) {
611  X_est[k] = fracDelta_old[k] / sumFrac;
612  if (Vphase->spGlobalIndexVCS(k) < m_numComponents) {
613  sum_Xcomp += X_est[k];
614  }
615  }
616 
617 
618  /*
619  * Feed the newly formed estimate of the mole fractions back into the
620  * ThermoPhase object
621  */
623 
624  /*
625  * get the activity coefficients
626  */
627  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
628 
629  /*
630  * First calculate altered chemical potentials for component species
631  * belonging to this phase.
632  */
633  for (size_t i = 0; i < componentList.size(); i++) {
634  size_t kc = componentList[i];
635  size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
636  if (X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
637  m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
638  + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
639  } else {
640  m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
641  + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
642  }
643  }
644 
645  for (size_t i = 0; i < componentList.size(); i++) {
646  size_t kc_spec = Vphase->spGlobalIndexVCS(componentList[i]);
647  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
648  size_t kspec = Vphase->spGlobalIndexVCS(k);
649  if (kspec >= m_numComponents) {
650  size_t irxn = kspec - m_numComponents;
651  if (i == 0) {
652  m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn];
653  }
654  if (m_stoichCoeffRxnMatrix(kc_spec,irxn) != 0.0) {
655  m_deltaGRxn_Deficient[irxn] +=
656  m_stoichCoeffRxnMatrix(kc_spec,irxn) * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
657  }
658  }
659 
660  }
661  }
662 
663  /*
664  * Calculate the E_phi's
665  */
666  sum = 0.0;
667  funcPhaseStability = sum_Xcomp - 1.0;
668  for (size_t k = 0; k < nsp; k++) {
669  size_t kspec = Vphase->spGlobalIndexVCS(k);
670  if (kspec >= m_numComponents) {
671  size_t irxn = kspec - m_numComponents;
672  double deltaGRxn = clip(m_deltaGRxn_Deficient[irxn], -50.0, 50.0);
673  E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec];
674  sum += E_phi[k];
675  funcPhaseStability += E_phi[k];
676  } else {
677  E_phi[k] = 0.0;
678  }
679  }
680 
681  /*
682  * Calculate the raw estimate of the new fracs
683  */
684  for (size_t k = 0; k < nsp; k++) {
685  size_t kspec = Vphase->spGlobalIndexVCS(k);
686  double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
687  if (kspec >= m_numComponents) {
688  fracDelta_raw[k] = b;
689  }
690  }
691 
692 
693  // Given a set of fracDelta's, we calculate the fracDelta's
694  // for the component species, if any
695  for (size_t i = 0; i < componentList.size(); i++) {
696  size_t kc = componentList[i];
697  size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
698  fracDelta_raw[kc] = 0.0;
699  for (size_t k = 0; k < nsp; k++) {
700  size_t kspec = Vphase->spGlobalIndexVCS(k);
701  if (kspec >= m_numComponents) {
702  size_t irxn = kspec - m_numComponents;
703  fracDelta_raw[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_raw[k];
704  }
705  }
706  }
707 
708 
709 
710  /*
711  * Now possibly dampen the estimate.
712  */
713  doublereal sumADel = 0.0;
714  for (size_t k = 0; k < nsp; k++) {
715  delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
716  sumADel += fabs(delFrac[k]);
717  }
718  normUpdate = vcs_l2norm(delFrac);
719 
720  dirProd = 0.0;
721  for (size_t k = 0; k < nsp; k++) {
722  dirProd += fracDelta_old[k] * delFrac[k];
723  }
724  bool crossedSign = false;
725  if (dirProd * dirProdOld < 0.0) {
726  crossedSign = true;
727  }
728 
729 
730  damp = 0.5;
731  if (dampOld < 0.25) {
732  damp = 2.0 * dampOld;
733  }
734  if (crossedSign) {
735  if (normUpdate *1.5 > normUpdateOld) {
736  damp = 0.5 * dampOld;
737  } else if (normUpdate *2.0 > normUpdateOld) {
738  damp = 0.8 * dampOld;
739  }
740  } else {
741  if (normUpdate > normUpdateOld * 2.0) {
742  damp = 0.6 * dampOld;
743  } else if (normUpdate > normUpdateOld * 1.2) {
744  damp = 0.9 * dampOld;
745  }
746  }
747 
748  for (size_t k = 0; k < nsp; k++) {
749  if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
750  damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
751  1.0E-8/fabs(delFrac[k]));
752  }
753  if (delFrac[k] < 0.0) {
754  if (2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
755  damp = fracDelta_old[k] / (2.0 * (-delFrac[k]));
756  }
757  }
758  if (delFrac[k] > 0.0) {
759  if (2.0 * damp * delFrac[k] > fracDelta_old[k]) {
760  damp = fracDelta_old[k] / (2.0 * delFrac[k]);
761  }
762  }
763  }
764  damp = std::max(damp, 0.000001);
765  for (size_t k = 0; k < nsp; k++) {
766  fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]);
767  }
768 
769  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
770  plogf(" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
771  delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
772  }
773 
774  if (normUpdate < 1.0E-5 * damp) {
775  converged = true;
776  if (its < minNumberIterations) {
777  converged = false;
778  }
779  }
780 
781  }
782 
783  if (converged) {
784  /*
785  * Save the final optimized stated back into the VolPhase object for later use
786  */
788  /*
789  * Save fracDelta for later use to initialize the problem better
790  * @TODO creationGlobalRxnNumbers needs to be calculated here and stored.
791  */
792  Vphase->setCreationMoleNumbers(VCS_DATA_PTR(fracDelta_new), creationGlobalRxnNumbers);
793  }
794 
795 
796  } else {
797  throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet");
798  }
799  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
800  plogf(" ------------------------------------------------------------"
801  "-------------------------------------------------------------\n");
802  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl == 1) {
803  if (funcPhaseStability > 0.0) {
804  plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
805  } else {
806  plogf(" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability);
807  }
808  }
809  return funcPhaseStability;
810 }
811 
812 }
void setCreationMoleNumbers(const double *const n_k, const std::vector< size_t > &creationGlobalRxnNumbers)
Sets the creationMoleNum's within the phase object.
const std::vector< double > & creationMoleNumbers(std::vector< size_t > &creationGlobalRxnNumbers) const
Return a const reference to the creationMoleNumbers stored in the object.
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:584
void setMoleFractionsState(const double molNum, const double *const moleFracVec, const int vcsStateStatus)
Set the moles and/or mole fractions within the phase.
#define VCS_SPECIES_MINOR
Species is a major species.
Definition: vcs_defs.h:135
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
#define VCS_DELETE_ELEMENTABS_CUTOFF
Cutoff moles below which a phase or species which comprises the bulk of an element's total concentrat...
Definition: vcs_defs.h:84
#define VCS_DATA_PTR(vvv)
Points to the data in a std::vector<> object.
Definition: vcs_internal.h:18
size_t nSpecies() const
Return the number of species in the phase.
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:270
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:696
#define VCS_STATECALC_PHASESTABILITY
State Calculation based on tentative mole numbers for a phase which is currently zeroed, but is being evaluated for whether it should pop back into existence.
Definition: vcs_defs.h:373
#define AssertThrowMsg(expr, procedure, message)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:301
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_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:290
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem...
Definition: vcs_defs.h:57
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:28
double vcs_l2norm(const std::vector< double > vec)
determine the l2 norm of a vector of doubles
Definition: vcs_util.cpp:23
#define VCS_DELETE_PHASE_CUTOFF
Cutoff relative moles below which a phase is deleted from the equilibrium problem.
Definition: vcs_defs.h:69
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
int exists() const
int indicating whether the phase exists or not
#define VCS_SPECIES_MAJOR
Species is a major species.
Definition: vcs_defs.h:128
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:86
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:365
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:24
#define VCS_SPECIES_COMPONENT
Species is a component which can be nonzero.
Definition: vcs_defs.h:121
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...