Cantera  2.5.1
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  */
6 
7 // This file is part of Cantera. See License.txt in the top-level directory or
8 // at https://cantera.org/license.txt for license and copyright information.
9 
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 
20 bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
21 {
22  vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop].get();
23  AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhasePossible",
24  "called for a phase that exists!");
25 
26  // Loop through all of the species in the phase. We say the phase can be
27  // popped, if there is one species in the phase that can be popped. This
28  // does not mean that the phase will be popped or that it leads to a lower
29  // Gibbs free energy.
30  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
31  size_t kspec = Vphase->spGlobalIndexVCS(k);
32  AssertThrowMsg(m_molNumSpecies_old[kspec] <= 0.0,
33  "VCS_SOLVE::vcs_popPhasePossible",
34  "we shouldn't be here {}: {} > 0.0", kspec,
35  m_molNumSpecies_old[kspec]);
36  size_t irxn = kspec - m_numComponents;
37  if (kspec >= m_numComponents) {
38  bool iPopPossible = true;
39 
40  // Note one case is if the component is a member of the popping
41  // phase. This component will be zeroed and the logic here will
42  // negate the current species from causing a positive if this
43  // component is consumed.
44  for (size_t j = 0; j < m_numComponents; ++j) {
45  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
46  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
47  if (stoicC != 0.0) {
48  double negChangeComp = - stoicC;
49  if (negChangeComp > 0.0) {
50  // If there is no component to give, then the
51  // species can't be created
52  if (m_molNumSpecies_old[j] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
53  iPopPossible = false;
54  }
55  }
56  }
57  }
58  }
59  // We are here when the species can be popped because all its needed
60  // components have positive mole numbers
61  if (iPopPossible) {
62  return true;
63  }
64  } else {
65  // We are here when the species, k, in the phase is a component. Its
66  // mole number is zero. We loop through the regular reaction looking
67  // for a reaction that can pop the component.
68  for (size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) {
69  bool foundJrxn = false;
70  // First, if the component is a product of the reaction
71  if (m_stoichCoeffRxnMatrix(kspec,jrxn) > 0.0) {
72  foundJrxn = true;
73  // We can do the reaction if all other reactant components
74  // have positive mole fractions
75  for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
76  if (m_stoichCoeffRxnMatrix(kcomp,jrxn) < 0.0 && m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
77  foundJrxn = false;
78  }
79  }
80  if (foundJrxn) {
81  return true;
82  }
83  } else if (m_stoichCoeffRxnMatrix(kspec,jrxn) < 0.0) {
84  // Second we are here if the component is a reactant in the
85  // reaction, and the reaction goes backwards.
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
93  // components species are positive
94  for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
95  if (m_stoichCoeffRxnMatrix(kcomp,jrxn) > 0.0 && m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
96  foundJrxn = false;
97  }
98  }
99  if (foundJrxn) {
100  return true;
101  }
102  }
103  }
104  }
105  }
106  return false;
107 }
108 
109 size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
110 {
111  size_t iphasePop = npos;
112  doublereal FephaseMax = -1.0E30;
113  doublereal Fephase = -1.0E30;
114 
115  char anote[128];
116  if (m_debug_print_lvl >= 2) {
117  plogf(" --- vcs_popPhaseID() called\n");
118  plogf(" --- Phase Status F_e MoleNum\n");
119  plogf(" --------------------------------------------------------------------------\n");
120  }
121  for (size_t iph = 0; iph < m_numPhases; iph++) {
122  vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
123  int existence = Vphase->exists();
124  strcpy(anote, "");
125  if (existence > 0) {
126  if (m_debug_print_lvl >= 2) {
127  plogf(" --- %18s %5d NA %11.3e\n",
128  Vphase->PhaseName, existence, m_tPhaseMoles_old[iph]);
129  }
130  } else {
131  if (Vphase->m_singleSpecies) {
132  // Single Phase Stability Resolution
133  size_t kspec = Vphase->spGlobalIndexVCS(0);
134  size_t irxn = kspec - m_numComponents;
135  if (irxn > m_deltaGRxn_old.size()) {
136  throw CanteraError("VCS_SOLVE::vcs_popPhaseID",
137  "Index out of bounds due to logic error.");
138  }
139  doublereal deltaGRxn = m_deltaGRxn_old[irxn];
140  Fephase = exp(-deltaGRxn) - 1.0;
141  if (Fephase > 0.0) {
142  strcpy(anote," (ready to be birthed)");
143  if (Fephase > FephaseMax) {
144  iphasePop = iph;
145  FephaseMax = Fephase;
146  strcpy(anote," (chosen to be birthed)");
147  }
148  }
149  if (Fephase < 0.0) {
150  strcpy(anote," (not stable)");
151  AssertThrowMsg(m_tPhaseMoles_old[iph] <= 0.0,
152  "VCS_SOLVE::vcs_popPhaseID", "shouldn't be here");
153  }
154 
155  if (m_debug_print_lvl >= 2) {
156  plogf(" --- %18s %5d %10.3g %10.3g %s\n",
157  Vphase->PhaseName, existence, Fephase,
158  m_tPhaseMoles_old[iph], anote);
159  }
160  } else {
161  // MultiSpecies Phase Stability Resolution
162  if (vcs_popPhasePossible(iph)) {
163  Fephase = vcs_phaseStabilityTest(iph);
164  if (Fephase > 0.0) {
165  if (Fephase > FephaseMax) {
166  iphasePop = iph;
167  FephaseMax = Fephase;
168  }
169  } else {
170  FephaseMax = std::max(FephaseMax, Fephase);
171  }
172  if (m_debug_print_lvl >= 2) {
173  plogf(" --- %18s %5d %11.3g %11.3g\n",
174  Vphase->PhaseName, existence, Fephase,
175  m_tPhaseMoles_old[iph]);
176  }
177  } else {
178  if (m_debug_print_lvl >= 2) {
179  plogf(" --- %18s %5d blocked %11.3g\n",
180  Vphase->PhaseName,
181  existence, m_tPhaseMoles_old[iph]);
182  }
183  }
184  }
185  }
186  }
187  phasePopPhaseIDs.resize(0);
188  if (iphasePop != npos) {
189  phasePopPhaseIDs.push_back(iphasePop);
190  }
191 
192  // Insert logic here to figure out if phase pops are linked together. Only
193  // do one linked pop at a time.
194  if (m_debug_print_lvl >= 2) {
195  plogf(" ---------------------------------------------------------------------\n");
196  }
197  return iphasePop;
198 }
199 
200 int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
201 {
202  vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop].get();
203  // Identify the first species in the phase
204  size_t kspec = Vphase->spGlobalIndexVCS(0);
205  // Identify the formation reaction for that species
206  size_t irxn = kspec - m_numComponents;
207  std::vector<size_t> creationGlobalRxnNumbers;
208 
209  // Calculate the initial moles of the phase being born.
210  // Here we set it to 10x of the value which would cause the phase to be
211  // zeroed out within the algorithm. We may later adjust the value.
212  doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
213 
214  AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhaseRxnStepSizes",
215  "called for a phase that exists!");
216  if (m_debug_print_lvl >= 2) {
217  plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
218  Vphase->PhaseName, iphasePop);
219  }
220  // Section for a single-species phase
221  if (Vphase->m_singleSpecies) {
222  double s = 0.0;
223  for (size_t j = 0; j < m_numComponents; ++j) {
224  if (!m_SSPhase[j] && m_molNumSpecies_old[j] > 0.0) {
225  s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j];
226  }
227  }
228  for (size_t j = 0; j < m_numPhases; j++) {
229  Vphase = m_VolPhaseList[j].get();
230  if (! Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) {
231  s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j];
232  }
233  }
234  if (s != 0.0) {
235  double s_old = s;
236  s = vcs_Hessian_diag_adj(irxn, s_old);
237  m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
238  } else {
239  // Ok, s is equal to zero. We can not apply a sophisticated theory
240  // to birth the phase. Just pick a small delta and go with it.
241  m_deltaMolNumSpecies[kspec] = tPhaseMoles;
242  }
243 
244  // section to do damping of the m_deltaMolNumSpecies[]
245  for (size_t j = 0; j < m_numComponents; ++j) {
246  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
247  if (stoicC != 0.0 && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
248  double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
249  if (negChangeComp > m_molNumSpecies_old[j]) {
250  if (m_molNumSpecies_old[j] > 0.0) {
251  m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
252  } else {
253  m_deltaMolNumSpecies[kspec] = 0.0;
254  }
255  }
256  }
257  }
258  // Implement a damping term that limits m_deltaMolNumSpecies to the size
259  // of the mole number
260  if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
261  m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
262  }
263  } else {
264  vector_fp fracDelta(Vphase->nSpecies());
265  vector_fp X_est(Vphase->nSpecies());
266  fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
267 
268  double sumFrac = 0.0;
269  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
270  sumFrac += fracDelta[k];
271  }
272  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
273  X_est[k] = fracDelta[k] / sumFrac;
274  }
275 
276  doublereal deltaMolNumPhase = tPhaseMoles;
277  doublereal damp = 1.0;
278  m_deltaGRxn_tmp = m_molNumSpecies_old;
279  double* molNumSpecies_tmp = m_deltaGRxn_tmp.data();
280 
281  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
282  kspec = Vphase->spGlobalIndexVCS(k);
283  double delmol = deltaMolNumPhase * X_est[k];
284  if (kspec >= m_numComponents) {
285  irxn = kspec - m_numComponents;
286  if (irxn > m_stoichCoeffRxnMatrix.nColumns()) {
287  throw CanteraError("VCS_SOLVE::vcs_popPhaseRxnStepSizes",
288  "Index out of bounds due to logic error.");
289  }
290  for (size_t j = 0; j < m_numComponents; ++j) {
291  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
292  if (stoicC != 0.0 && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
293  molNumSpecies_tmp[j] += stoicC * delmol;
294  }
295  }
296  }
297  }
298 
299  doublereal ratioComp = 0.0;
300  for (size_t j = 0; j < m_numComponents; ++j) {
301  double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j];
302  if (molNumSpecies_tmp[j] < 0.0) {
303  ratioComp = 1.0;
304  if (deltaJ > 0.0) {
305  double delta0 = m_molNumSpecies_old[j];
306  damp = std::min(damp, delta0 / deltaJ * 0.9);
307  }
308  } else {
309  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
310  size_t jph = m_phaseID[j];
311  if ((jph != iphasePop) && (!m_SSPhase[j])) {
312  double fdeltaJ = fabs(deltaJ);
313  if (m_molNumSpecies_old[j] > 0.0) {
314  ratioComp = std::max(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]);
315  }
316  }
317  }
318  }
319  }
320 
321  // We may have greatly underestimated the deltaMoles for the phase pop
322  // Here we create a damp > 1 to account for this possibility. We adjust
323  // upwards to make sure that a component in an existing multispecies
324  // phase is modified by a factor of 1/1000.
325  if (ratioComp > 1.0E-30 && ratioComp < 0.001) {
326  damp = 0.001 / ratioComp;
327  }
328  if (damp <= 1.0E-6) {
329  return 3;
330  }
331 
332  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
333  kspec = Vphase->spGlobalIndexVCS(k);
334  if (kspec < m_numComponents) {
335  m_speciesStatus[kspec] = VCS_SPECIES_COMPONENT;
336  } else {
337  m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp;
338  if (X_est[k] > 1.0E-3) {
339  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
340  } else {
341  m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
342  }
343  }
344  }
345  }
346  return 0;
347 }
348 
349 double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
350 {
351  // We will use the _new state calc here
352  vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
353  const size_t nsp = Vphase->nSpecies();
354  int minNumberIterations = 3;
355  if (nsp <= 1) {
356  minNumberIterations = 1;
357  }
358 
359  // We will do a full Newton calculation later, but for now, ...
360  bool doSuccessiveSubstitution = true;
361  double funcPhaseStability;
362  vector_fp X_est(nsp, 0.0);
363  vector_fp delFrac(nsp, 0.0);
364  vector_fp E_phi(nsp, 0.0);
365  vector_fp fracDelta_old(nsp, 0.0);
366  vector_fp fracDelta_raw(nsp, 0.0);
367  vector<size_t> creationGlobalRxnNumbers(nsp, npos);
368  m_deltaGRxn_Deficient = m_deltaGRxn_old;
369  vector_fp feSpecies_Deficient = m_feSpecies_old;
370 
371  // get the activity coefficients
372  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, &m_actCoeffSpecies_new[0]);
373 
374  // Get the stored estimate for the composition of the phase if
375  // it gets created
376  vector_fp fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
377 
378  std::vector<size_t> componentList;
379  for (size_t k = 0; k < nsp; k++) {
380  size_t kspec = Vphase->spGlobalIndexVCS(k);
381  if (kspec < m_numComponents) {
382  componentList.push_back(k);
383  }
384  }
385 
386  double normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
387  double damp = 1.0E-2;
388 
389  if (doSuccessiveSubstitution) {
390  int KP = 0;
391  if (m_debug_print_lvl >= 2) {
392  plogf(" --- vcs_phaseStabilityTest() called\n");
393  plogf(" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
394  " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
395  plogf(" --------------------------------------------------------------"
396  "--------------------------------------------------------\n");
397  } else if (m_debug_print_lvl == 1) {
398  plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
399  }
400 
401  for (size_t k = 0; k < nsp; k++) {
402  if (fracDelta_new[k] < 1.0E-13) {
403  fracDelta_new[k] =1.0E-13;
404  }
405  }
406  bool converged = false;
407  double dirProd = 0.0;
408  for (int its = 0; its < 200 && (!converged); its++) {
409  double dampOld = damp;
410  double normUpdateOld = normUpdate;
411  fracDelta_old = fracDelta_new;
412  double dirProdOld = dirProd;
413 
414  // Given a set of fracDelta's, we calculate the fracDelta's
415  // for the component species, if any
416  for (size_t i = 0; i < componentList.size(); i++) {
417  size_t kc = componentList[i];
418  size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
419  fracDelta_old[kc] = 0.0;
420  for (size_t k = 0; k < nsp; k++) {
421  size_t kspec = Vphase->spGlobalIndexVCS(k);
422  if (kspec >= m_numComponents) {
423  size_t irxn = kspec - m_numComponents;
424  fracDelta_old[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_old[k];
425  }
426  }
427  }
428 
429  // Now, calculate the predicted mole fractions, X_est[k]
430  double sumFrac = 0.0;
431  for (size_t k = 0; k < nsp; k++) {
432  sumFrac += fracDelta_old[k];
433  }
434  // Necessary because this can be identically zero. -> we need to fix
435  // this algorithm!
436  if (sumFrac <= 0.0) {
437  sumFrac = 1.0;
438  }
439  double sum_Xcomp = 0.0;
440  for (size_t k = 0; k < nsp; k++) {
441  X_est[k] = fracDelta_old[k] / sumFrac;
442  if (Vphase->spGlobalIndexVCS(k) < m_numComponents) {
443  sum_Xcomp += X_est[k];
444  }
445  }
446 
447  // Feed the newly formed estimate of the mole fractions back into the
448  // ThermoPhase object
449  Vphase->setMoleFractionsState(0.0, &X_est[0], VCS_STATECALC_PHASESTABILITY);
450 
451  // get the activity coefficients
452  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, &m_actCoeffSpecies_new[0]);
453 
454  // First calculate altered chemical potentials for component species
455  // belonging to this phase.
456  for (size_t i = 0; i < componentList.size(); i++) {
457  size_t kc = componentList[i];
458  size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
459  if (X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
460  feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
461  + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
462  } else {
463  feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
464  + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
465  }
466  }
467 
468  for (size_t i = 0; i < componentList.size(); i++) {
469  size_t kc_spec = Vphase->spGlobalIndexVCS(componentList[i]);
470  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
471  size_t kspec = Vphase->spGlobalIndexVCS(k);
472  if (kspec >= m_numComponents) {
473  size_t irxn = kspec - m_numComponents;
474  if (i == 0) {
475  m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn];
476  }
477  if (m_stoichCoeffRxnMatrix(kc_spec,irxn) != 0.0) {
478  m_deltaGRxn_Deficient[irxn] +=
479  m_stoichCoeffRxnMatrix(kc_spec,irxn) * (feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
480  }
481  }
482  }
483  }
484 
485  // Calculate the E_phi's
486  double sum = 0.0;
487  funcPhaseStability = sum_Xcomp - 1.0;
488  for (size_t k = 0; k < nsp; k++) {
489  size_t kspec = Vphase->spGlobalIndexVCS(k);
490  if (kspec >= m_numComponents) {
491  size_t irxn = kspec - m_numComponents;
492  double deltaGRxn = clip(m_deltaGRxn_Deficient[irxn], -50.0, 50.0);
493  E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec];
494  sum += E_phi[k];
495  funcPhaseStability += E_phi[k];
496  } else {
497  E_phi[k] = 0.0;
498  }
499  }
500 
501  // Calculate the raw estimate of the new fracs
502  for (size_t k = 0; k < nsp; k++) {
503  size_t kspec = Vphase->spGlobalIndexVCS(k);
504  double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
505  if (kspec >= m_numComponents) {
506  fracDelta_raw[k] = b;
507  }
508  }
509 
510  // Given a set of fracDelta's, we calculate the fracDelta's
511  // for the component species, if any
512  for (size_t i = 0; i < componentList.size(); i++) {
513  size_t kc = componentList[i];
514  size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
515  fracDelta_raw[kc] = 0.0;
516  for (size_t k = 0; k < nsp; k++) {
517  size_t kspec = Vphase->spGlobalIndexVCS(k);
518  if (kspec >= m_numComponents) {
519  size_t irxn = kspec - m_numComponents;
520  fracDelta_raw[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_raw[k];
521  }
522  }
523  }
524 
525  // Now possibly dampen the estimate.
526  doublereal sumADel = 0.0;
527  for (size_t k = 0; k < nsp; k++) {
528  delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
529  sumADel += fabs(delFrac[k]);
530  }
531  normUpdate = vcs_l2norm(delFrac);
532 
533  dirProd = 0.0;
534  for (size_t k = 0; k < nsp; k++) {
535  dirProd += fracDelta_old[k] * delFrac[k];
536  }
537  bool crossedSign = false;
538  if (dirProd * dirProdOld < 0.0) {
539  crossedSign = true;
540  }
541 
542  damp = 0.5;
543  if (dampOld < 0.25) {
544  damp = 2.0 * dampOld;
545  }
546  if (crossedSign) {
547  if (normUpdate *1.5 > normUpdateOld) {
548  damp = 0.5 * dampOld;
549  } else if (normUpdate *2.0 > normUpdateOld) {
550  damp = 0.8 * dampOld;
551  }
552  } else {
553  if (normUpdate > normUpdateOld * 2.0) {
554  damp = 0.6 * dampOld;
555  } else if (normUpdate > normUpdateOld * 1.2) {
556  damp = 0.9 * dampOld;
557  }
558  }
559 
560  for (size_t k = 0; k < nsp; k++) {
561  if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
562  damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
563  1.0E-8/fabs(delFrac[k]));
564  }
565  if (delFrac[k] < 0.0 && 2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
566  damp = fracDelta_old[k] / (2.0 * -delFrac[k]);
567  }
568  if (delFrac[k] > 0.0 && 2.0 * damp * delFrac[k] > fracDelta_old[k]) {
569  damp = fracDelta_old[k] / (2.0 * delFrac[k]);
570  }
571  }
572  damp = std::max(damp, 0.000001);
573  for (size_t k = 0; k < nsp; k++) {
574  fracDelta_new[k] = fracDelta_old[k] + damp * delFrac[k];
575  }
576 
577  if (m_debug_print_lvl >= 2) {
578  plogf(" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
579  delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
580  }
581 
582  if (normUpdate < 1.0E-5 * damp) {
583  converged = true;
584  if (its < minNumberIterations) {
585  converged = false;
586  }
587  }
588  }
589 
590  if (converged) {
591  // Save the final optimized stated back into the VolPhase object for later use
592  Vphase->setMoleFractionsState(0.0, &X_est[0], VCS_STATECALC_PHASESTABILITY);
593 
594  // Save fracDelta for later use to initialize the problem better
595  // @TODO creationGlobalRxnNumbers needs to be calculated here and stored.
596  Vphase->setCreationMoleNumbers(&fracDelta_new[0], creationGlobalRxnNumbers);
597  }
598  } else {
599  throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet");
600  }
601  if (m_debug_print_lvl >= 2) {
602  plogf(" ------------------------------------------------------------"
603  "-------------------------------------------------------------\n");
604  } else if (m_debug_print_lvl == 1) {
605  if (funcPhaseStability > 0.0) {
606  plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
607  } else {
608  plogf(" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability);
609  }
610  }
611  return funcPhaseStability;
612 }
613 
614 }
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:82
size_t nSpecies() const
Return the number of species in the phase.
void setMoleFractionsState(const double molNum, const double *const moleFracVec, const int vcsStateStatus)
Set the moles and/or mole fractions within the phase.
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:541
int exists() const
int indicating whether the phase exists or not
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:627
void setCreationMoleNumbers(const double *const n_k, const std::vector< size_t > &creationGlobalRxnNumbers)
Sets the creationMoleNum's within the phase object.
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
const vector_fp & creationMoleNumbers(std::vector< size_t > &creationGlobalRxnNumbers) const
Return a const reference to the creationMoleNumbers stored in the object.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:265
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:300
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:188
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
double vcs_l2norm(const vector_fp &vec)
determine the l2 norm of a vector of doubles
Definition: vcs_util.cpp:21
Contains declarations for string manipulation functions within Cantera.
Header for the object representing each phase within vcs.
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:300
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem.
Definition: vcs_defs.h:44
#define VCS_SPECIES_COMPONENT
Species is a component which can be nonzero.
Definition: vcs_defs.h:95
#define VCS_SPECIES_MAJOR
Species is a major species.
Definition: vcs_defs.h:102
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:228
#define VCS_DELETE_PHASE_CUTOFF
Cutoff relative moles below which a phase is deleted from the equilibrium problem.
Definition: vcs_defs.h:56
#define VCS_SPECIES_MINOR
Species is a major species.
Definition: vcs_defs.h:109
#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:70
#define VCS_STATECALC_PHASESTABILITY
State Calculation based on tentative mole numbers for a phase which is currently zeroed,...
Definition: vcs_defs.h:308
#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...