Cantera  2.0
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  */
8 #include "vcs_species_thermo.h"
10 #include "vcs_Exception.h"
11 
12 #include <cstdio>
13 #include <cstdlib>
14 #include <cmath>
15 #include <vector>
16 #include <cstring>
17 #include <algorithm>
18 
19 using namespace std;
20 
21 namespace VCSnonideal
22 {
23 
24 //====================================================================================================================
25 // Utility function that evaluates whether a phase can be popped into existence
26 /*
27  * A phase can be popped iff the stoichiometric coefficients for the
28  * component species, whose concentrations will be lowered during the
29  * process, are positive by at least a small degree.
30  *
31  * If one of the phase species is a zeroed component, then the phase can
32  * be popped if the component increases in mole number as the phase moles
33  * are increased.
34  *
35  * @param iphasePop id of the phase, which is currently zeroed,
36  *
37  * @return Returns true if the phase can come into existence
38  * and false otherwise.
39  */
40 bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
41 {
42 
43  vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
44 
45 #ifdef DEBUG_MODE
46  int existence = Vphase->exists();
47  if (existence > 0) {
48  printf("ERROR vcs_popPhasePossible called for a phase that exists!");
49  std::exit(-1);
50  }
51 #endif
52 
53  /*
54  * Loop through all of the species in the phase. We say the phase
55  * can be popped, if there is one species in the phase that can be
56  * popped.
57  */
58  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
59  size_t kspec = Vphase->spGlobalIndexVCS(k);
60 #ifdef DEBUG_MODE
61  if (m_molNumSpecies_old[kspec] > 0.0) {
62  printf("ERROR vcs_popPhasePossible we shouldn't be here %lu %g > 0.0",
63  kspec, m_molNumSpecies_old[kspec]);
64  exit(-1);
65  }
66 #endif
67  size_t irxn = kspec - m_numComponents;
68  if (kspec >= m_numComponents) {
69  bool iPopPossible = true;
70  for (size_t j = 0; j < m_numComponents; ++j) {
71  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
72  double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
73  if (stoicC != 0.0) {
74  double negChangeComp = - stoicC * 1.0;
75  if (negChangeComp > 0.0) {
76  // TODO: We may have to come up with a tolerance here
77  if (m_molNumSpecies_old[j] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
78  iPopPossible = false;
79  }
80  }
81  }
82  }
83  }
84  if (iPopPossible) {
85  return true;
86  }
87  } else {
88  /*
89  * We are here when the species in the phase is a component. Its mole number is zero.
90  * We loop through the regular reaction looking for a reaction that can pop the
91  * component.
92  */
93  //printf("WE are here at new logic - CHECK\n");
94  for (size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) {
95  bool foundJrxn = false;
96  // First, if the component is a product of the reaction
97  if (m_stoichCoeffRxnMatrix[jrxn][kspec] > 0.0) {
98  foundJrxn = true;
99  for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
100  if (m_stoichCoeffRxnMatrix[jrxn][kcomp] < 0.0) {
101  if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
102  foundJrxn = false;
103  }
104  }
105  }
106  if (foundJrxn) {
107  //printf("We have found a component phase pop! CHECK1 \n");
108  return true;
109  }
110  }
111  // Second we are here if the component is a reactant in the reaction, and the reaction goes backwards.
112  else if (m_stoichCoeffRxnMatrix[jrxn][kspec] < 0.0) {
113  foundJrxn = true;
114  size_t jspec = jrxn + m_numComponents;
115  if (m_molNumSpecies_old[jspec] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
116  foundJrxn = false;
117  continue;
118  }
119  for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
120  if (m_stoichCoeffRxnMatrix[jrxn][kcomp] > 0.0) {
121  if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
122  foundJrxn = false;
123  }
124  }
125  }
126  if (foundJrxn) {
127  //printf("We have found a component phase pop! CHECK2 \n");
128  return true;
129  }
130  }
131  }
132  }
133  }
134  return false;
135 }
136 
137 //====================================================================================================================
138 // Determine the list of problems that need to be checked to see if there are any phases pops
139 /*
140  * This routine evaluates and fills in the following quantities
141  * phasePopProblemLists_
142  *
143  * Need to work in species that are zeroed by element constraints
144  *
145  * @return Returns the number of problems that must be checked.
146  */
147 int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
148 {
149  int nfound = 0;
150  vcs_VolPhase* Vphase = 0;
151  double stoicC;
152  double molComp;
153  std::vector<int> linkedPhases;
154  phasePopProblemLists_.clear();
155 
156  /*
157  * This is a vector over each component.
158  * For zeroed components it lists the phases, which are currently zeroed,
159  * which have a species with a positive stoichiometric value wrt the component.
160  * Therefore, we could pop the component species and pop that phase at the same time
161  * if we considered no other factors than keeping the component mole number positve.
162  *
163  * It does not count species with positive stoichiometric values if that species
164  * already has a positive mole number. The phase is already popped.
165  */
166  std::vector< std::vector<size_t> > zeroedComponentLinkedPhasePops(m_numComponents);
167  /*
168  * The logic below calculates zeroedComponentLinkedPhasePops
169  */
170  for (size_t j = 0; j < m_numComponents; j++) {
171  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
172  molComp = m_molNumSpecies_old[j];
173  if (molComp <= 0.0) {
174  std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
175  size_t iph = m_phaseID[j];
176  jList.push_back(iph);
177  for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
178  size_t kspec = irxn + m_numComponents;
179  iph = m_phaseID[kspec];
180  Vphase = m_VolPhaseList[iph];
181  int existence = Vphase->exists();
182  if (existence < 0) {
183  stoicC = m_stoichCoeffRxnMatrix[irxn][j];
184  if (stoicC > 0.0) {
185  if (std::find(jList.begin(), jList.end(), iph) != jList.end()) {
186  jList.push_back(iph);
187  }
188  }
189  }
190  }
191  }
192  }
193  }
194  /*
195  * This is a vector over each zeroed phase
196  * For zeroed phases, it lists the components, which are currently zereoed,
197  * which have a species with a negative stoichiometric value wrt one or more species in the phase.
198  * Cut out components which have a pos stoichiometric value with another species in the phase.
199  */
200  std::vector< std::vector<size_t> > zeroedPhaseLinkedZeroComponents(m_numPhases);
201  /*
202  * The logic below calculates zeroedPhaseLinkedZeroComponents
203  */
204  for (size_t iph = 0; iph < m_numPhases; iph++) {
205  std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
206  iphList.clear();
207  Vphase = m_VolPhaseList[iph];
208  int existence = Vphase->exists();
209  if (existence < 0) {
210 
211  linkedPhases.clear();
212  size_t nsp = Vphase->nSpecies();
213  for (size_t k = 0; k < nsp; k++) {
214  size_t kspec = Vphase->spGlobalIndexVCS(k);
215  size_t irxn = kspec - m_numComponents;
216 
217  for (size_t j = 0; j < m_numComponents; j++) {
218  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
219  molComp = m_molNumSpecies_old[j];
220  if (molComp <= 0.0) {
221  stoicC = m_stoichCoeffRxnMatrix[irxn][j];
222  if (stoicC < 0.0) {
223  bool foundPos = false;
224  for (size_t kk = 0; kk < nsp; kk++) {
225  size_t kkspec = Vphase->spGlobalIndexVCS(kk);
226  if (kkspec >= m_numComponents) {
227  size_t iirxn = kkspec - m_numComponents;
228  if (m_stoichCoeffRxnMatrix[iirxn][j] > 0.0) {
229  foundPos = true;
230  }
231  }
232  }
233  if (!foundPos) {
234  if (std::find(iphList.begin(), iphList.end(), j) != iphList.end()) {
235  iphList.push_back(j);
236  }
237  }
238  }
239  }
240  }
241  }
242  }
243  }
244  }
245 
246  /*
247  * Now fill in the phasePopProblemLists_ list.
248  *
249  */
250  for (size_t iph = 0; iph < m_numPhases; iph++) {
251  Vphase = m_VolPhaseList[iph];
252  int existence = Vphase->exists();
253  if (existence < 0) {
254  std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
255  std::vector<size_t> popProblem(0);
256  popProblem.push_back(iph);
257  for (size_t i = 0; i < iphList.size(); i++) {
258  size_t j = iphList[i];
259  std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
260  for (size_t jjl = 0; jjl < jList.size(); jjl++) {
261  size_t jph = jList[jjl];
262  if (std::find(popProblem.begin(), popProblem.end(), jph) != popProblem.end()) {
263  popProblem.push_back(jph);
264  }
265  }
266  }
267  phasePopProblemLists_.push_back(popProblem);
268  }
269  }
270 
271  return nfound;
272 }
273 
274 
275 //====================================================================================================================
276 // Decision as to whether a phase pops back into existence
277 /*
278  * @return returns the phase id of the phases that pops back into
279  * existence. Returns npos if there are no phases
280  */
281 size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
282 {
283  size_t iphasePop = npos;
284  size_t irxn, kspec;
285  doublereal FephaseMax = -1.0E30;
286  doublereal Fephase = -1.0E30;
287  vcs_VolPhase* Vphase = 0;
288 
289 
290 #ifdef DEBUG_MODE
291  char anote[128];
292  if (m_debug_print_lvl >= 2) {
293  plogf(" --- vcs_popPhaseID() called\n");
294  plogf(" --- Phase Status F_e MoleNum\n");
295  plogf(" --------------------------------------------------------------------------\n");
296  }
297 #endif
298  for (size_t iph = 0; iph < m_numPhases; iph++) {
299  Vphase = m_VolPhaseList[iph];
300  int existence = Vphase->exists();
301 #ifdef DEBUG_MODE
302  strcpy(anote, "");
303 #endif
304  if (existence > 0) {
305 
306 #ifdef DEBUG_MODE
307  if (m_debug_print_lvl >= 2) {
308  plogf(" --- %18s %5d NA %11.3e\n",
309  Vphase->PhaseName.c_str(),
310  existence,
311  m_tPhaseMoles_old[iph]);
312  }
313 #endif
314  } else {
315  if (Vphase->m_singleSpecies) {
316  /***********************************************************************
317  *
318  * Single Phase Stability Resolution
319  *
320  ***********************************************************************/
321  kspec = Vphase->spGlobalIndexVCS(0);
322  irxn = kspec - m_numComponents;
323  doublereal deltaGRxn = m_deltaGRxn_old[irxn];
324  Fephase = exp(-deltaGRxn) - 1.0;
325  if (Fephase > 0.0) {
326 #ifdef DEBUG_MODE
327  strcpy(anote," (ready to be birthed)");
328 #endif
329  if (Fephase > FephaseMax) {
330  iphasePop = iph;
331  FephaseMax = Fephase;
332 #ifdef DEBUG_MODE
333  strcpy(anote," (chosen to be birthed)");
334 #endif
335  }
336  }
337 #ifdef DEBUG_MODE
338  if (Fephase < 0.0) {
339  strcpy(anote," (not stable)");
340  if (m_tPhaseMoles_old[iph] > 0.0) {
341  printf("shouldn't be here\n");
342  exit(-1);
343  }
344  }
345 #endif
346 
347 #ifdef DEBUG_MODE
348  if (m_debug_print_lvl >= 2) {
349  plogf(" --- %18s %5d %10.3g %10.3g %s\n",
350  Vphase->PhaseName.c_str(),
351  existence, Fephase,
352  m_tPhaseMoles_old[iph], anote);
353  }
354 #endif
355 
356  } else {
357  /***********************************************************************
358  *
359  * MultiSpecies Phase Stability Resolution
360  *
361  ***********************************************************************/
362  if (vcs_popPhasePossible(iph)) {
363  Fephase = vcs_phaseStabilityTest(iph);
364  if (Fephase > 0.0) {
365  if (Fephase > FephaseMax) {
366  iphasePop = iph;
367  FephaseMax = Fephase;
368  }
369  } else {
370  if (Fephase > FephaseMax) {
371  FephaseMax = Fephase;
372  }
373  }
374 #ifdef DEBUG_MODE
375  if (m_debug_print_lvl >= 2) {
376  plogf(" --- %18s %5d %11.3g %11.3g\n",
377  Vphase->PhaseName.c_str(),
378  existence, Fephase,
379  m_tPhaseMoles_old[iph]);
380  }
381 #endif
382  } else {
383 #ifdef DEBUG_MODE
384  if (m_debug_print_lvl >= 2) {
385  plogf(" --- %18s %5d blocked %11.3g\n",
386  Vphase->PhaseName.c_str(),
387  existence, m_tPhaseMoles_old[iph]);
388  }
389 #endif
390  }
391  }
392  }
393  }
394  phasePopPhaseIDs.resize(0);
395  if (iphasePop != npos) {
396  phasePopPhaseIDs.push_back(iphasePop);
397  }
398 
399  /*
400  * Insert logic here to figure out if phase pops are linked together. Only do one linked
401  * pop at a time.
402  */
403 
404 #ifdef DEBUG_MODE
405  if (m_debug_print_lvl >= 2) {
406  plogf(" ---------------------------------------------------------------------\n");
407  }
408 #endif
409  return iphasePop;
410 }
411 //====================================================================================================================
412 // Calculates the deltas of the reactions due to phases popping
413 // into existence
414 /*
415  * @param iphasePop Phase id of the phase that will come into existence
416  *
417  * Output
418  * -------
419  * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers
420  * to the irxn'th species
421  * formation reaction. This adjustment
422  * is for species
423  * irxn + M, where M is the number
424  * of components.
425  *
426  * @return Returns an int representing the status of the step
427  * - 0 : normal return
428  * - 1 : A single species phase species has been zeroed out
429  * in this routine. The species is a noncomponent
430  * - 2 : Same as one but, the zeroed species is a component.
431  * - 3 : Nothing was done because the phase couldn't be birthed
432  * because a needed component is zero.
433  */
434 int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
435 {
436  vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
437  // Identify the first species in the phase
438  size_t kspec = Vphase->spGlobalIndexVCS(0);
439  // Identify the formation reaction for that species
440  size_t irxn = kspec - m_numComponents;
441  std::vector<size_t> creationGlobalRxnNumbers;
442 
443  doublereal s;
444  // Calculate the initial moles of the phase being born.
445  // Here we set it to 10x of the value which would cause the phase to be
446  // zeroed out within the algorithm. We may later adjust the value.
447  doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
448 
449 
450 #ifdef DEBUG_MODE
451  int existence = Vphase->exists();
452  if (existence > 0) {
453  printf("ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!");
454  exit(-1);
455  }
456  char anote[256];
457  if (m_debug_print_lvl >= 2) {
458  plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
459  Vphase->PhaseName.c_str(), iphasePop);
460  }
461 #endif
462  // Section for a single-species phase
463  if (Vphase->m_singleSpecies) {
464  s = 0.0;
465  double* dnPhase_irxn = m_deltaMolNumPhase[irxn];
466  for (size_t j = 0; j < m_numComponents; ++j) {
467  if (!m_SSPhase[j]) {
468  if (m_molNumSpecies_old[j] > 0.0) {
469  s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
470  }
471  }
472  }
473  for (size_t j = 0; j < m_numPhases; j++) {
474  Vphase = m_VolPhaseList[j];
475  if (! Vphase->m_singleSpecies) {
476  if (m_tPhaseMoles_old[j] > 0.0) {
477  s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
478  }
479  }
480  }
481  if (s != 0.0) {
482  double s_old = s;
483  s = vcs_Hessian_diag_adj(irxn, s_old);
484 #ifdef DEBUG_MODE
485  if (s_old != s) {
486  sprintf(anote, "Normal calc: diag adjusted from %g "
487  "to %g due to act coeff", s_old, s);
488  }
489 #endif
490  m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
491  } else {
492  // Ok, s is equal to zero. We can not apply a sophisticated theory
493  // to birth the phase. Just pick a small delta and go with it.
494  m_deltaMolNumSpecies[kspec] = tPhaseMoles;
495  }
496 
497  /*
498  * section to do damping of the m_deltaMolNumSpecies[]
499  */
500  for (size_t j = 0; j < m_numComponents; ++j) {
501  double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
502  if (stoicC != 0.0) {
503  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
504  double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
505  if (negChangeComp > m_molNumSpecies_old[j]) {
506  if (m_molNumSpecies_old[j] > 0.0) {
507 #ifdef DEBUG_MODE
508  sprintf(anote, "Delta damped from %g "
509  "to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec],
510  -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
511 #endif
512  m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
513  } else {
514 #ifdef DEBUG_MODE
515  sprintf(anote, "Delta damped from %g "
516  "to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec],
517  -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
518 #endif
519  m_deltaMolNumSpecies[kspec] = 0.0;
520  }
521  }
522  }
523  }
524  }
525  // Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number
526  if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
527 #ifdef DEBUG_MODE
528  sprintf(anote, "Delta damped from %g "
529  "to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
530  -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
531 #endif
532  m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
533  }
534 
535 
536  } else {
537  vector<doublereal> fracDelta(Vphase->nSpecies());
538  vector<doublereal> X_est(Vphase->nSpecies());
539  fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
540 
541  double sumFrac = 0.0;
542  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
543  sumFrac += fracDelta[k];
544  }
545  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
546  X_est[k] = fracDelta[k] / sumFrac;
547  }
548 
549  doublereal deltaMolNumPhase = tPhaseMoles;
550  doublereal damp = 1.0;
551  m_deltaGRxn_tmp = m_molNumSpecies_old;
552  double* molNumSpecies_tmp = DATA_PTR(m_deltaGRxn_tmp);
553 
554 
555  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
556  kspec = Vphase->spGlobalIndexVCS(k);
557  double delmol = deltaMolNumPhase * X_est[k];
558  if (kspec >= m_numComponents) {
559  irxn = kspec - m_numComponents;
560  for (size_t j = 0; j < m_numComponents; ++j) {
561  double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
562  if (stoicC != 0.0) {
563  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
564  molNumSpecies_tmp[j] += stoicC * delmol;
565  }
566  }
567  }
568  }
569  }
570 
571  doublereal ratioComp = 0.0;
572  for (size_t j = 0; j < m_numComponents; ++j) {
573  double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j];
574  if (molNumSpecies_tmp[j] < 0.0) {
575  ratioComp = 1.0;
576  if (deltaJ > 0.0) {
577  double delta0 = m_molNumSpecies_old[j];
578  double dampj = delta0 / deltaJ * 0.9;
579  if (dampj < damp) {
580  damp = dampj;
581  }
582  }
583  } else {
584  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
585  size_t jph = m_phaseID[j];
586  if ((jph != iphasePop) && (!m_SSPhase[j])) {
587  double fdeltaJ = fabs(deltaJ);
588  if (m_molNumSpecies_old[j] > 0.0) {
589  ratioComp = std::max(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]);
590  }
591  }
592  }
593  }
594  }
595 
596  // We may have greatly underestimated the deltaMoles for the phase pop
597  // Here we create a damp > 1 to account for this possibility.
598  // We adjust upwards to make sure that a component in an existing multispecies
599  // phase is modified by a factor of 1/1000.
600  if (ratioComp > 1.0E-30) {
601  if (ratioComp < 0.001) {
602  damp = 0.001 / ratioComp;
603  }
604  }
605 
606 
607  if (damp <= 1.0E-6) {
608  return 3;
609  }
610 
611  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
612  kspec = Vphase->spGlobalIndexVCS(k);
613  if (kspec < m_numComponents) {
614  m_speciesStatus[kspec] = VCS_SPECIES_COMPONENT;
615  } else {
616  m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp;
617  if (X_est[k] > 1.0E-3) {
618  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
619  } else {
620  m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
621  }
622  }
623  }
624 
625  }
626 
627  return 0;
628 }
629 
630 
631 double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
632 {
633 
634  /*
635  * We will use the _new state calc here
636  */
637  size_t kspec, irxn, k, i, kc, kc_spec;
638  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
639  doublereal deltaGRxn;
640 
641  // We will do a full newton calculation later, but for now, ...
642  bool doSuccessiveSubstitution = true;
643  double funcPhaseStability;
644  vector<doublereal> X_est(Vphase->nSpecies(), 0.0);
645  vector<doublereal> delFrac(Vphase->nSpecies(), 0.0);
646  vector<doublereal> E_phi(Vphase->nSpecies(), 0.0);
647  vector<doublereal> fracDelta_new(Vphase->nSpecies(), 0.0);
648  vector<doublereal> fracDelta_old(Vphase->nSpecies(), 0.0);
649  vector<doublereal> fracDelta_raw(Vphase->nSpecies(), 0.0);
650  vector<size_t> creationGlobalRxnNumbers(Vphase->nSpecies(), npos);
651  vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
652 
653  vector<doublereal> m_feSpecies_Deficient(m_numComponents, 0.0);
654  doublereal damp = 1.0;
655  doublereal dampOld = 1.0;
656  doublereal normUpdate = 1.0;
657  doublereal normUpdateOld = 1.0;
658  doublereal sum = 0.0;
659  doublereal dirProd = 0.0;
660  doublereal dirProdOld = 0.0;
661 
662  // get the activity coefficients
663  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
664 
665  // Get the stored estimate for the composition of the phase if
666  // it gets created
667  fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
668 
669 
670  std::vector<size_t> componentList;
671 
672  for (k = 0; k < Vphase->nSpecies(); k++) {
673  kspec = Vphase->spGlobalIndexVCS(k);
674  if (kspec < m_numComponents) {
675  componentList.push_back(k);
676  }
677  }
678 
679  for (k = 0; k < m_numComponents; k++) {
680  m_feSpecies_Deficient[k] = m_feSpecies_old[k];
681  }
682  normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
683  damp = 1.0E-2;
684 
685  if (doSuccessiveSubstitution) {
686 
687 #ifdef DEBUG_MODE
688  int KP = 0;
689  if (m_debug_print_lvl >= 2) {
690  plogf(" --- vcs_phaseStabilityTest() called\n");
691  plogf(" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
692  " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
693  plogf(" --------------------------------------------------------------"
694  "--------------------------------------------------------\n");
695  } else if (m_debug_print_lvl == 1) {
696  plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
697  }
698 #endif
699 
700  for (k = 0; k < Vphase->nSpecies(); k++) {
701  if (fracDelta_new[k] < 1.0E-13) {
702  fracDelta_new[k] = 1.0E-13;
703  }
704  }
705  bool converged = false;
706  for (int its = 0; its < 200 && (!converged); its++) {
707 
708  dampOld = damp;
709  normUpdateOld = normUpdate;
710  fracDelta_old = fracDelta_new;
711  dirProdOld = dirProd;
712 
713  // Given a set of fracDelta's, we calculate the fracDelta's
714  // for the component species, if any
715  for (i = 0; i < componentList.size(); i++) {
716  kc = componentList[i];
717  kc_spec = Vphase->spGlobalIndexVCS(kc);
718  fracDelta_old[kc] = 0.0;
719  for (k = 0; k < Vphase->nSpecies(); k++) {
720  kspec = Vphase->spGlobalIndexVCS(k);
721  if (kspec >= m_numComponents) {
722  irxn = kspec - m_numComponents;
723  fracDelta_old[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_old[k];
724  }
725  }
726  }
727 
728  // Now, calculate the predicted mole fractions, X_est[k]
729  double sumFrac = 0.0;
730  for (k = 0; k < Vphase->nSpecies(); k++) {
731  sumFrac += fracDelta_old[k];
732  }
733  // Necessary because this can be identically zero. -> we need to fix this algorithm!
734  if (sumFrac <= 0.0) {
735  sumFrac = 1.0;
736  }
737  double sum_Xcomp = 0.0;
738  for (k = 0; k < Vphase->nSpecies(); k++) {
739  X_est[k] = fracDelta_old[k] / sumFrac;
740  kc_spec = Vphase->spGlobalIndexVCS(k);
741  if (kc_spec < m_numComponents) {
742  sum_Xcomp += X_est[k];
743  }
744  }
745 
746 
747  /*
748  * Feed the newly formed estimate of the mole fractions back into the
749  * ThermoPhase object
750  */
752 
753  /*
754  * get the activity coefficients
755  */
756  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
757 
758  /*
759  * First calculate altered chemical potentials for component species
760  * belonging to this phase.
761  */
762  for (i = 0; i < componentList.size(); i++) {
763  kc = componentList[i];
764  kc_spec = Vphase->spGlobalIndexVCS(kc);
765  if (X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
766  m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
767  + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
768  } else {
769  m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
770  + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
771  }
772  }
773 
774  for (i = 0; i < componentList.size(); i++) {
775  kc = componentList[i];
776  kc_spec = Vphase->spGlobalIndexVCS(kc);
777 
778  for (k = 0; k < Vphase->nSpecies(); k++) {
779  kspec = Vphase->spGlobalIndexVCS(k);
780  if (kspec >= m_numComponents) {
781  irxn = kspec - m_numComponents;
782  if (i == 0) {
783  m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn];
784  }
785  double* dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
786  if (dtmp_ptr[kc_spec] != 0.0) {
787  m_deltaGRxn_Deficient[irxn] +=
788  dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
789  }
790  }
791 
792  }
793  }
794 
795  /*
796  * Calculate the E_phi's
797  */
798  sum = 0.0;
799  funcPhaseStability = sum_Xcomp - 1.0;
800  for (k = 0; k < Vphase->nSpecies(); k++) {
801  kspec = Vphase->spGlobalIndexVCS(k);
802  if (kspec >= m_numComponents) {
803  irxn = kspec - m_numComponents;
804  deltaGRxn = m_deltaGRxn_Deficient[irxn];
805  if (deltaGRxn > 50.0) {
806  deltaGRxn = 50.0;
807  }
808  if (deltaGRxn < -50.0) {
809  deltaGRxn = -50.0;
810  }
811  E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec];
812  sum += E_phi[k];
813  funcPhaseStability += E_phi[k];
814  } else {
815  E_phi[k] = 0.0;
816  }
817  }
818 
819  /*
820  * Calculate the raw estimate of the new fracs
821  */
822  for (k = 0; k < Vphase->nSpecies(); k++) {
823  kspec = Vphase->spGlobalIndexVCS(k);
824  double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
825  if (kspec >= m_numComponents) {
826  irxn = kspec - m_numComponents;
827  fracDelta_raw[k] = b;
828  }
829  }
830 
831 
832  // Given a set of fracDelta's, we calculate the fracDelta's
833  // for the component species, if any
834  for (i = 0; i < componentList.size(); i++) {
835  kc = componentList[i];
836  kc_spec = Vphase->spGlobalIndexVCS(kc);
837  fracDelta_raw[kc] = 0.0;
838  for (k = 0; k < Vphase->nSpecies(); k++) {
839  kspec = Vphase->spGlobalIndexVCS(k);
840  if (kspec >= m_numComponents) {
841  irxn = kspec - m_numComponents;
842  fracDelta_raw[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_raw[k];
843  }
844  }
845  }
846 
847 
848 
849  /*
850  * Now possibly dampen the estimate.
851  */
852  doublereal sumADel = 0.0;
853  for (k = 0; k < Vphase->nSpecies(); k++) {
854  delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
855  sumADel += fabs(delFrac[k]);
856  }
857  normUpdate = vcs_l2norm(delFrac);
858 
859  dirProd = 0.0;
860  for (k = 0; k < Vphase->nSpecies(); k++) {
861  dirProd += fracDelta_old[k] * delFrac[k];
862  }
863  bool crossedSign = false;
864  if (dirProd * dirProdOld < 0.0) {
865  crossedSign = true;
866  }
867 
868 
869  damp = 0.5;
870  if (dampOld < 0.25) {
871  damp = 2.0 * dampOld;
872  }
873  if (crossedSign) {
874  if (normUpdate *1.5 > normUpdateOld) {
875  damp = 0.5 * dampOld;
876  } else if (normUpdate *2.0 > normUpdateOld) {
877  damp = 0.8 * dampOld;
878  }
879  } else {
880  if (normUpdate > normUpdateOld * 2.0) {
881  damp = 0.6 * dampOld;
882  } else if (normUpdate > normUpdateOld * 1.2) {
883  damp = 0.9 * dampOld;
884  }
885  }
886 
887  for (k = 0; k < Vphase->nSpecies(); k++) {
888  if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
889  damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
890  1.0E-8/fabs(delFrac[k]));
891  }
892  if (delFrac[k] < 0.0) {
893  if (2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
894  damp = fracDelta_old[k] / (2.0 * (-delFrac[k]));
895  }
896  }
897  if (delFrac[k] > 0.0) {
898  if (2.0 * damp * delFrac[k] > fracDelta_old[k]) {
899  damp = fracDelta_old[k] / (2.0 * delFrac[k]);
900  }
901  }
902  }
903  if (damp < 0.000001) {
904  damp = 0.000001;
905  }
906 
907  for (k = 0; k < Vphase->nSpecies(); k++) {
908  fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]);
909  }
910 
911 #ifdef DEBUG_MODE
912  if (m_debug_print_lvl >= 2) {
913  plogf(" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
914  delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
915  }
916 #endif
917 
918  if (normUpdate < 1.0E-5) {
919  converged = true;
920  }
921 
922  }
923 
924  if (converged) {
925  Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est),
927  Vphase->setCreationMoleNumbers(VCS_DATA_PTR(fracDelta_new), creationGlobalRxnNumbers);
928  }
929 
930 
931  } else {
932  printf("not done yet\n");
933  exit(-1);
934  }
935 #ifdef DEBUG_MODE
936  if (m_debug_print_lvl >= 2) {
937  plogf(" ------------------------------------------------------------"
938  "-------------------------------------------------------------\n");
939  } else if (m_debug_print_lvl == 1) {
940  if (funcPhaseStability > 0.0) {
941  plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
942  } else {
943  plogf(" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability);
944  }
945  }
946 #endif
947  return funcPhaseStability;
948 }
949 
950 }
951