Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
vcs_solve_TP.cpp
Go to the documentation of this file.
1 /**
2  * @file vcs_solve_TP.cpp Implementation file that contains the
3  * main algorithm for finding an equilibrium
4  */
5 /*
6  * Copyright (2005) Sandia Corporation. Under the terms of
7  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
8  * U.S. Government retains certain rights in this software.
9  */
10 
11 #include "cantera/equil/vcs_solve.h"
12 #include "cantera/equil/vcs_VolPhase.h"
13 
14 #include "cantera/base/ctexceptions.h"
15 #include "cantera/base/clockWC.h"
16 #include "cantera/base/stringUtils.h"
17 #include "cantera/numerics/ctlapack.h"
18 
19 #include <cstdio>
20 
21 using namespace std;
22 namespace {
23 enum stages {MAIN, EQUILIB_CHECK, ELEM_ABUND_CHECK,
24  RECHECK_DELETED, RETURN_A, RETURN_B};
25 }
26 
27 namespace Cantera
28 {
29 
30 void VCS_SOLVE::checkDelta1(double* const dsLocal,
31  double* const delTPhMoles, size_t kspec)
32 {
33  std::vector<double> dchange(m_numPhases, 0.0);
34  for (size_t k = 0; k < kspec; k++) {
35  if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
36  size_t iph = m_phaseID[k];
37  dchange[iph] += dsLocal[k];
38  }
39  }
40  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
41  double denom = max(m_totalMolNum, 1.0E-4);
42  if (!vcs_doubleEqual(dchange[iphase]/denom, delTPhMoles[iphase]/denom)) {
43  throw CanteraError("VCS_SOLVE::checkDelta1",
44  "we have found a problem");
45  }
46  }
47 }
48 
49 int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
50 {
51  int stage = MAIN;
52  int solveFail;
53  bool allMinorZeroedSpecies = false;
54  size_t it1 = 0;
55  size_t npb, iti;
56  int rangeErrorFound = 0;
57  bool giveUpOnElemAbund = false;
58  int finalElemAbundAttempts = 0;
59  bool uptodate_minors = true;
60  int forceComponentCalc = 1;
61 
62 #ifdef DEBUG_MODE
63  char ANOTE[128];
64  /*
65  * Set the debug print lvl to the same as the print lvl.
66  */
67  m_debug_print_lvl = printDetails;
68 #else
69  char* ANOTE = 0;
70 #endif
71  if (printDetails > 0 && print_lvl == 0) {
72  print_lvl = 1;
73  }
74  /*
75  * Initialize and set up all counters
76  */
77  vcs_counters_init(0);
78  clockWC ticktock;
79 
80  /*
81  * Malloc temporary space for usage in this routine and in
82  * subroutines
83  * sm[ne*ne]
84  * ss[ne]
85  * sa[ne]
86  * aw[m]
87  * wx[ne]
88  * xy[m]
89  */
90 
91  m_sm.assign(m_numElemConstraints*m_numElemConstraints, 0.0);
92  m_ss.assign(m_numElemConstraints, 0.0);
93  m_sa.assign(m_numElemConstraints, 0.0);
94  m_aw.assign(m_numSpeciesTot, 0.0);
95  m_wx.assign(m_numElemConstraints, 0.0);
96 
97  solveFail = false;
98 
99  /* ****************************************************** */
100  /* **** Evaluate the elemental composition ****** */
101  /* ****************************************************** */
102  vcs_elab();
103 
104  /* ******************************************************* */
105  /* **** Printout the initial conditions for problem ****** */
106  /* ******************************************************* */
107 
108  if (print_lvl != 0) {
109  plogf("VCS CALCULATION METHOD\n\n ");
110  plogf("%s\n", m_title.c_str());
111  plogf("\n\n%5d SPECIES\n%5d ELEMENTS\n", m_numSpeciesTot, m_numElemConstraints);
112  plogf("%5d COMPONENTS\n", m_numComponents);
113  plogf("%5d PHASES\n", m_numPhases);
114 
115  plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa ");
116  plogf(" TEMPERATURE%19.3f K\n", m_temperature);
117  vcs_VolPhase* Vphase = m_VolPhaseList[0];
118  if (Vphase->nSpecies() > 0) {
119  plogf(" PHASE1 INERTS%17.3f\n", TPhInertMoles[0]);
120  }
121  if (m_numPhases > 1) {
122  plogf(" PHASE2 INERTS%17.3f\n", TPhInertMoles[1]);
123  }
124  plogf("\n ELEMENTAL ABUNDANCES CORRECT");
125  plogf(" FROM ESTIMATE Type\n\n");
126  for (size_t i = 0; i < m_numElemConstraints; ++i) {
127  writeline(' ', 26, false);
128  plogf("%-2.2s", (m_elementName[i]).c_str());
129  plogf("%20.12E%20.12E %3d\n", m_elemAbundancesGoal[i], m_elemAbundances[i],
130  m_elType[i]);
131  }
132  if (m_doEstimateEquil < 0) {
133  plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - forced\n");
134  } else if (m_doEstimateEquil > 0) {
135  plogf("\n MODIFIED LINEAR PROGRAMMING ESTIMATE OF EQUILIBRIUM - where necessary\n");
136  }
137  if (m_doEstimateEquil == 0) {
138  plogf("\n USER ESTIMATE OF EQUILIBRIUM\n");
139  }
140  if (m_VCS_UnitsFormat == VCS_UNITS_KCALMOL) {
141  plogf(" Stan. Chem. Pot. in kcal/mole\n");
142  }
143  if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
144  plogf(" Stan. Chem. Pot. is MU/RT\n");
145  }
146  if (m_VCS_UnitsFormat == VCS_UNITS_KJMOL) {
147  plogf(" Stan. Chem. Pot. in KJ/mole\n");
148  }
149  if (m_VCS_UnitsFormat == VCS_UNITS_KELVIN) {
150  plogf(" Stan. Chem. Pot. in Kelvin\n");
151  }
152  if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
153  plogf(" Stan. Chem. Pot. in J/kmol\n");
154  }
155  plogf("\n SPECIES FORMULA VECTOR ");
156  writeline(' ', 41, false);
157  plogf(" STAN_CHEM_POT EQUILIBRIUM_EST. Species_Type\n\n");
158  writeline(' ', 20, false);
159  for (size_t i = 0; i < m_numElemConstraints; ++i) {
160  plogf("%-4.4s ", m_elementName[i].c_str());
161  }
162  plogf(" PhaseID\n");
163  double RT = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
164  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
165  plogf(" %-18.18s", m_speciesName[i].c_str());
166  for (size_t j = 0; j < m_numElemConstraints; ++j) {
167  plogf("% -7.3g ", m_formulaMatrix(i,j));
168  }
169  plogf(" %3d ", m_phaseID[i]);
170  writeline(' ', std::max(55-int(m_numElemConstraints)*8, 0), false);
171  plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]);
172  if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
173  plogf(" Mol_Num");
174  } else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
175  plogf(" Voltage");
176  } else {
177  plogf(" Unknown");
178  }
179  plogendl();
180  }
181  }
182 
183  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
184  if (m_molNumSpecies_old[i] < 0.0) {
185  plogf("On Input species %-12s has a "
186  "negative MF, setting it small",
187  m_speciesName[i].c_str());
188  plogendl();
189  size_t iph = m_phaseID[i];
190  double tmp = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10;
191  tmp = std::max(tmp, VCS_DELETE_MINORSPECIES_CUTOFF*10.);
192  m_molNumSpecies_old[i] = tmp;
193  }
194  }
195 
196  /*
197  * Evaluate the total moles of species in the problem
198  */
199  vcs_tmoles();
200 
201  /*
202  * Evaluate all chemical potentials at the old mole numbers at the
203  * outset of the calculation.
204  */
205  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
206  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
207 
208  bool lec;
209  while (true) {
210  if (stage == MAIN) {
211  /* *********************************************************** */
212  /* **** DETERMINE BASIS SPECIES, EVALUATE STOICHIOMETRY ****** */
213  /* *********************************************************** */
214  if (forceComponentCalc) {
215  int retn = solve_tp_component_calc(allMinorZeroedSpecies);
216  if (retn != VCS_SUCCESS) {
217  return retn;
218  }
219 
220  it1 = 1;
221  forceComponentCalc = 0;
222  iti = 0;
223  }
224  /*
225  * Check on too many iterations.
226  * If we have too many iterations, Clean up and exit code even though we haven't
227  * converged. -> we have run out of iterations!
228  */
229  if (m_VCount->Its > maxit) {
230  return -1;
231  }
232  solve_tp_inner(iti, it1, uptodate_minors, allMinorZeroedSpecies,
233  forceComponentCalc, stage, printDetails, ANOTE);
234  lec = false;
235 
236  } else if (stage == EQUILIB_CHECK) {
237  /*************************************************************************/
238  /***************** EQUILIBRIUM CHECK FOR MAJOR SPECIES *******************/
239  /*************************************************************************/
240  solve_tp_equilib_check(allMinorZeroedSpecies, uptodate_minors,
241  giveUpOnElemAbund, solveFail, iti, it1,
242  maxit, stage, lec);
243 
244  } else if (stage == ELEM_ABUND_CHECK) {
245  /* *************************************************** */
246  /* **** CORRECT ELEMENTAL ABUNDANCES ***************** */
247  /* *************************************************** */
248  solve_tp_elem_abund_check(iti, stage, lec, giveUpOnElemAbund,
249  finalElemAbundAttempts, rangeErrorFound);
250  } else if (stage == RECHECK_DELETED) {
251  /* *************************************************** */
252  /* **** RECHECK DELETED SPECIES ********************** */
253  /* *************************************************** */
254  /*
255  * We are here for two reasons. One is if we have
256  * achieved convergence, but some species have been eliminated
257  * from the problem because they were in multispecies phases
258  * and their mole fractions drifted less than
259  * VCS_DELETE_SPECIES_CUTOFF .
260  * The other reason why we are here is because all of the
261  * non-component species in the problem have been eliminated
262  * for one reason or another.
263  */
264  npb = vcs_recheck_deleted();
265  /*
266  * If we haven't found any species that needed adding we are done.
267  */
268  if (npb <= 0) {
269  stage = RETURN_B;
270  } else {
271  /*
272  * If we have found something to add, recalculate everything
273  * for minor species and go back to do a full iteration
274  */
275  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
276  vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
277  vcs_deltag(0, false, VCS_STATECALC_OLD);
278  iti = 0;
279  stage = MAIN;
280  }
281  } else if (stage == RETURN_A) {
282  /*************************************************************************/
283  /******************** CLEANUP AND RETURN BLOCK ***************************/
284  /*************************************************************************/
285  npb = vcs_recheck_deleted();
286  /*
287  * If we haven't found any species that needed adding we are done.
288  */
289  if (npb > 0) {
290  /*
291  * If we have found something to add, recalculate everything
292  * for minor species and go back to do a full iteration
293  */
294  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
295  vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
296  vcs_deltag(0, false, VCS_STATECALC_OLD);
297  iti = 0;
298  stage = MAIN;
299  } else {
300  stage = RETURN_B;
301  }
302  } else if (stage == RETURN_B) {
303  /*
304  * Add back deleted species in non-zeroed phases. Estimate their
305  * mole numbers.
306  */
307  npb = vcs_add_all_deleted();
308  if (npb > 0) {
309  iti = 0;
310  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) {
311  plogf(" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!");
312  plogendl();
313  }
314  stage = MAIN;
315  } else {
316  break;
317  }
318  }
319  }
320  /*
321  * Make sure the volume phase objects hold the same state and
322  * information as the vcs object. This also update the Cantera objects
323  * with this information.
324  */
325  vcs_updateVP(VCS_STATECALC_OLD);
326  /*
327  * Evaluate the final mole fractions
328  * storing them in wt[]
329  */
330  m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0);
331  for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
332  if (m_SSPhase[kspec]) {
333  m_molNumSpecies_new[kspec] = 1.0;
334  } else {
335  size_t iph = m_phaseID[kspec];
336  if (m_tPhaseMoles_old[iph] != 0.0) {
337  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / m_tPhaseMoles_old[iph];
338  } else {
339  /*
340  * For MultiSpecies phases that are zeroed out,
341  * return the mole fraction vector from the VolPhase object.
342  * This contains the mole fraction that would be true if
343  * the phase just pops into existence.
344  */
345  size_t i = m_speciesLocalPhaseIndex[kspec];
346  m_molNumSpecies_new[kspec] = m_VolPhaseList[iph]->molefraction(i);
347  }
348  }
349  }
350  // Return an error code if a Range Space Error is thought to have occurred.
351  if (rangeErrorFound) {
352  solveFail = 1;
353  }
354  /*
355  * Calculate counters
356  */
357  double tsecond = ticktock.secondsWC();
358  m_VCount->Time_vcs_TP = tsecond;
359  m_VCount->T_Time_vcs_TP += m_VCount->Time_vcs_TP;
360  (m_VCount->T_Calls_vcs_TP)++;
361  m_VCount->T_Its += m_VCount->Its;
362  m_VCount->T_Basis_Opts += m_VCount->Basis_Opts;
363  m_VCount->T_Time_basopt += m_VCount->Time_basopt;
364  /*
365  * Return a Flag indicating whether convergence occurred
366  */
367  return solveFail;
368 }
369 
370 int VCS_SOLVE::solve_tp_component_calc(bool& allMinorZeroedSpecies)
371 {
372  double test = -1.0e-10;
373  bool usedZeroedSpecies;
374  int retn = vcs_basopt(false, VCS_DATA_PTR(m_aw), VCS_DATA_PTR(m_sa),
375  VCS_DATA_PTR(m_sm), VCS_DATA_PTR(m_ss),
376  test, &usedZeroedSpecies);
377  if (retn != VCS_SUCCESS) {
378  return retn;
379  }
380 
381  // Update the phase objects with the contents of the soln vector
382  vcs_updateVP(VCS_STATECALC_OLD);
383  vcs_deltag(0, false, VCS_STATECALC_OLD);
384 
385  /*************************************************************************/
386  /************** EVALUATE INITIAL SPECIES STATUS VECTOR *******************/
387  /*************************************************************************/
388  allMinorZeroedSpecies = vcs_evaluate_speciesType();
389 
390  /*************************************************************************/
391  /************** EVALUATE THE ELELEMT ABUNDANCE CHECK ******************/
392  /*************************************************************************/
393  if (! vcs_elabcheck(0)) {
394  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
395  plogf(" --- Element Abundance check failed");
396  plogendl();
397  }
398  vcs_elcorr(VCS_DATA_PTR(m_sm), VCS_DATA_PTR(m_wx));
399  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
400  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
401  // Update the phase objects with the contents of the soln vector
402  vcs_updateVP(VCS_STATECALC_OLD);
403  vcs_deltag(0, false, VCS_STATECALC_OLD);
404  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
405  plogf(" --- Element Abundance check passed");
406  plogendl();
407  }
408  return retn;
409 }
410 
411 void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
412  bool& uptodate_minors,
413  bool& allMinorZeroedSpecies,
414  int& forceComponentCalc,
415  int& stage, bool printDetails, char* ANOTE)
416 {
417  if (iti == 0) {
418  /* ********************************************************* */
419  /* **** SET INITIAL VALUES FOR ITERATION ******************* */
420  /* **** EVALUATE REACTION ADJUSTMENTS ******************* */
421  /* ********************************************************* */
422  /*
423  * Evaluate the minor non-component species chemical
424  * potentials and delta G for their formation reactions
425  * We have already evaluated the major non-components
426  */
427  if (!uptodate_minors) {
428  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
429  vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
430  vcs_deltag(1, false, VCS_STATECALC_OLD);
431  }
432  uptodate_minors = true;
433  } else {
434  uptodate_minors = false;
435  }
436 
437  if (printDetails) {
438  plogf("\n");
439  writeline('=', 110);
440  plogf(" Iteration = %3d, Iterations since last evaluation of "
441  "optimal basis = %3d",
442  m_VCount->Its, it1 - 1);
443  if (iti == 0) {
444  plogf(" (all species)\n");
445  } else {
446  plogf(" (only major species)\n");
447  }
448  }
449  /*
450  * Calculate the total moles in each phase -> old solution
451  * -> Needed for numerical stability when phases disappear.
452  * -> the phase moles tend to drift off without this step.
453  */
454 #ifdef DEBUG_MODE
455  check_tmoles();
456 #endif
457  vcs_tmoles();
458  /*************************************************************************/
459  /************** COPY OLD into NEW and ZERO VECTORS ***********************/
460  /*************************************************************************/
461  /*
462  * Copy the old solution into the new solution as an initial guess
463  */
464  m_feSpecies_new = m_feSpecies_old;
465  m_actCoeffSpecies_new = m_actCoeffSpecies_old;
466  m_deltaGRxn_new = m_deltaGRxn_old;
467  m_deltaGRxn_Deficient = m_deltaGRxn_old;
468  m_tPhaseMoles_new = m_tPhaseMoles_old;
469 
470  /*
471  * Zero out the entire vector of updates. We sometimes would
472  * query these values below, and we want to be sure that no
473  * information is left from previous iterations.
474  */
475  m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0);
476 
477  /*************************************************************************/
478  /************** DETERMINE IF DEAD PHASES POP INTO EXISTENCE **************/
479  /*************************************************************************/
480  /*
481  * First step is a major branch in the algorithm.
482  * We first determine if a phase pops into existence.
483  */
484  std::vector<size_t> phasePopPhaseIDs(0);
485  size_t iphasePop = vcs_popPhaseID(phasePopPhaseIDs);
486  if (iphasePop != npos) {
487  int soldel = vcs_popPhaseRxnStepSizes(iphasePop);
488  if (soldel == 3) {
489  iphasePop = npos;
490  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
491  plogf(" --- vcs_popPhaseRxnStepSizes() was called but stoich "
492  "prevented phase %d popping\n");
493  }
494  }
495  }
496 
497  /*************************************************************************/
498  /* DETERMINE THE REACTION STEP SIZES FOR MAIN STEP AND IF PHASES DIE *****/
499  /*************************************************************************/
500  /*
501  * Don't do this step if there is a phase pop
502  */
503  size_t iphaseDelete = npos;
504  size_t kspec;
505  if (iphasePop == npos) {
506  /*
507  * Figure out the new reaction step sizes
508  * for the major species (do minor species in the future too)
509  */
510  kspec = npos;
511  iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec);
512  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
513  plogf(" --- vcs_RxnStepSizes not called because alternative"
514  "phase creation delta was used instead\n");
515  }
516  size_t doPhaseDeleteKspec = npos;
517  size_t doPhaseDeleteIph = npos;
518  /*
519  * Zero out the net change in moles of multispecies phases
520  */
521  m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
522 
523  /* ********************************************************************** */
524  /* ***************** MAIN LOOP IN CALCULATION *************************** */
525  /* ***************** LOOP OVER IRXN TO DETERMINE STEP SIZE ************** */
526  /* ********************************************************************** */
527  /*
528  * Loop through all of the reactions, irxn, pertaining to the
529  * formation reaction for species kspec in canonical form.
530  *
531  * At the end of this loop, we will have a new estimate for the
532  * mole numbers for all species consistent with an extent
533  * of reaction for all noncomponent species formation
534  * reactions. We will have also ensured that all predicted
535  * non-component mole numbers are greater than zero.
536  *
537  * Old_Solution New_Solution Description
538  * -----------------------------------------------------------------------------
539  * m_molNumSpecies_old[kspec] m_molNumSpecies_new[kspec] Species Mole Numbers
540  * m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers
541  */
542  if (iphaseDelete != npos) {
543  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
544  plogf(" --- Main Loop Treatment -> Circumvented due to Phase Deletion ");
545  plogendl();
546  }
547  for (size_t k = 0; k < m_numSpeciesTot; k++) {
548  m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k];
549  size_t iph = m_phaseID[k];
550  m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k];
551  }
552  if (kspec >= m_numComponents) {
553  if (m_molNumSpecies_new[m_numSpeciesTot] != 0.0) {
554  throw CanteraError("VCS_SOLVE::solve_tp_inner",
555  "we shouldn't be here!");
556  }
557  if (m_SSPhase[kspec] == 1) {
558  m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
559  } else {
560  throw CanteraError("VCS_SOLVE::solve_tp_inner",
561  "we shouldn't be here!");
562  }
563  ++m_numRxnMinorZeroed;
564  allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
565  }
566  /*
567  * Set the flags indicating the mole numbers in the vcs_VolPhase
568  * objects are out of date.
569  */
570  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
571 
572  /*
573  * Calculate the new chemical potentials using the tentative
574  * solution values. We only calculate a subset of these, because
575  * we have only updated a subset of the W().
576  */
577  vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot);
578 
579  /*
580  * Evaluate DeltaG for all components if ITI=0, and for
581  * major components only if ITI NE 0
582  */
583  vcs_deltag(0, false, VCS_STATECALC_NEW);
584  } else {
585  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
586  plogf(" --- Main Loop Treatment of each non-component species ");
587  if (iti == 0) {
588  plogf("- Full Calculation:\n");
589  } else {
590  plogf("- Major Components Calculation:\n");
591  }
592  plogf(" --- Species IC ");
593  plogf(" KMoles Tent_KMoles Rxn_Adj | Comment \n");
594  }
595  for (size_t irxn = 0; irxn < m_numRxnRdc; irxn++) {
596  size_t kspec = m_indexRxnToSpecies[irxn];
597  double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
598  size_t iph = m_phaseID[kspec];
599  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
600  if (DEBUG_MODE_ENABLED) {
601  ANOTE[0] = '\0';
602  }
603  double dx;
604 
605  if (iphasePop != npos) {
606  if (iph == iphasePop) {
607  dx = m_deltaMolNumSpecies[kspec];
608  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
609  if (DEBUG_MODE_ENABLED) {
610  sprintf(ANOTE, "Phase pop");
611  }
612  } else {
613  dx = 0.0;
614  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
615  }
616  } else if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) {
617  /********************************************************************/
618  /************************ VOLTAGE SPECIES ***************************/
619  /********************************************************************/
620  bool soldel_ret;
621  dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
622  m_deltaMolNumSpecies[kspec] = dx;
623  } else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) {
624  /********************************************************************/
625  /********************** ZEROED OUT SPECIES **************************/
626  /********************************************************************/
627  bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0);
628  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 3) {
629  plogf(" --- %s currently zeroed (SpStatus=%-2d):",
630  m_speciesName[kspec].c_str(), m_speciesStatus[kspec]);
631  plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
632  irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec],
633  m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]);
634  }
635  if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) {
636  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
637  m_deltaMolNumSpecies[kspec] = 0.0;
638  resurrect = false;
639  if (DEBUG_MODE_ENABLED) {
640  sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]);
641  if (m_deltaGRxn_new[irxn] < 0.0) {
642  if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
643  sprintf(ANOTE, "Species stays zeroed even though dg neg due to "
644  "STOICH/PHASEPOP constraint: DG = %11.4E",
645  m_deltaGRxn_new[irxn]);
646  } else {
647  sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed",
648  m_deltaGRxn_new[irxn]);
649  }
650  }
651  }
652  } else {
653  for (size_t j = 0; j < m_numElemConstraints; ++j) {
654  int elType = m_elType[j];
655  if (elType == VCS_ELEM_TYPE_ABSPOS) {
656  double atomComp = m_formulaMatrix(kspec,j);
657  if (atomComp > 0.0) {
658  double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
659  if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
660  if (DEBUG_MODE_ENABLED) {
661  sprintf(ANOTE, "Species stays zeroed even though dG "
662  "neg, because of %s elemAbund",
663  m_elementName[j].c_str());
664  }
665  resurrect = false;
666  break;
667  }
668  }
669  }
670  }
671  }
672  /*
673  * Resurrect the species
674  */
675  if (resurrect) {
676  if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
677  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
678  plogf(" --- Zeroed species changed to major: ");
679  plogf("%-12s\n", m_speciesName[kspec].c_str());
680  }
681  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
682  allMinorZeroedSpecies = false;
683  } else {
684  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
685  plogf(" --- Zeroed species changed to minor: ");
686  plogf("%-12s\n", m_speciesName[kspec].c_str());
687  }
688  m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
689  }
690  if (m_deltaMolNumSpecies[kspec] > 0.0) {
691  dx = m_deltaMolNumSpecies[kspec] * 0.01;
692  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
693  } else {
694  m_molNumSpecies_new[kspec] = m_totalMolNum * VCS_DELETE_PHASE_CUTOFF * 10.;
695  dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
696  }
697  m_deltaMolNumSpecies[kspec] = dx;
698  if (DEBUG_MODE_ENABLED) {
699  sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]);
700  }
701  } else {
702  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
703  m_deltaMolNumSpecies[kspec] = 0.0;
704  dx = 0.0;
705  }
706  } else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
707  /********************************************************************/
708  /***************************** MINOR SPECIES ************************/
709  /********************************************************************/
710  /*
711  * Unless ITI isn't equal to zero we zero out changes
712  * to minor species.
713  */
714  if (iti != 0) {
715  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
716  m_deltaMolNumSpecies[kspec] = 0.0;
717  dx = 0.0;
718  if (DEBUG_MODE_ENABLED) {
719  sprintf(ANOTE,"minor species not considered");
720  if (m_debug_print_lvl >= 2) {
721  plogf(" --- ");
722  plogf("%-12s", m_speciesName[kspec].c_str());
723  plogf("%3d%11.4E%11.4E%11.4E | %s",
724  m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
725  m_deltaMolNumSpecies[kspec], ANOTE);
726  plogendl();
727  }
728  }
729  continue;
730  }
731  /*
732  * Minor species alternative calculation
733  * ---------------------------------------
734  * This is based upon the following approximation:
735  * The mole fraction changes due to these reactions don't affect
736  * the mole numbers of the component species. Therefore the
737  * following approximation is valid for an ideal solution
738  * 0 = DG(I) + log(WT(I)/W(I))
739  * (DG contains the contribution from FF(I) + log(W(I)/TL) )
740  * Thus,
741  * WT(I) = W(I) EXP(-DG(I))
742  * If soldel is true on return, then we branch to the section
743  * that deletes a species from the current set of active species.
744  */
745  bool soldel_ret;
746  dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE);
747  m_deltaMolNumSpecies[kspec] = dx;
748  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
749  if (soldel_ret) {
750  /*******************************************************************/
751  /***** DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF */
752  /***** MOLE NUMBER */
753  /*******************************************************************/
754  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
755  plogf(" --- Delete minor species in multispec phase: %-12s",
756  m_speciesName[kspec].c_str());
757  plogendl();
758  }
759  m_deltaMolNumSpecies[kspec] = 0.0;
760  /*
761  * Delete species, kspec. The alternate return is for the case
762  * where all species become deleted. Then, we need to
763  * branch to the code where we reevaluate the deletion
764  * of all species.
765  */
766  size_t lnospec = vcs_delete_species(kspec);
767 
768  if (lnospec) {
769  stage = RECHECK_DELETED;
770  break;
771  }
772  /*
773  * Go back to consider the next species in the list.
774  * Note, however, that the next species in the list is now
775  * in slot l. In deleting the previous species L, We have
776  * exchanged slot MR with slot l, and then have
777  * decremented MR.
778  * Therefore, we will decrement the species counter, here.
779  */
780  --irxn;
781  continue;
782  }
783  } else {
784  /********************************************************************/
785  /*********************** MAJOR SPECIES ******************************/
786  /********************************************************************/
787  if (DEBUG_MODE_ENABLED) {
788  sprintf(ANOTE, "Normal Major Calc");
789  }
790  /*
791  * Check for superconvergence of the formation reaction. Do
792  * nothing if it is superconverged. Skip to the end of the
793  * irxn loop if it is superconverged.
794  */
795  if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
796  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
797  m_deltaMolNumSpecies[kspec] = 0.0;
798  dx = 0.0;
799  if (DEBUG_MODE_ENABLED) {
800  sprintf(ANOTE, "major species is converged");
801  if (m_debug_print_lvl >= 2) {
802  plogf(" --- ");
803  plogf("%-12s", m_speciesName[kspec].c_str());
804  plogf("%3d%11.4E%11.4E%11.4E | %s",
805  m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
806  m_deltaMolNumSpecies[kspec], ANOTE);
807  plogendl();
808  }
809  }
810  continue;
811  }
812  /*
813  * Set the initial step size, dx, equal to the value produced
814  * by the routine, vcs_RxnStepSize().
815  *
816  * Note the multiplication logic is to make sure that
817  * dg[] didn't change sign due to w[] changing in the
818  * middle of the iteration. (it can if a single species
819  * phase goes out of existence).
820  */
821  if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) {
822  dx = m_deltaMolNumSpecies[kspec];
823  } else {
824  dx = 0.0;
825  m_deltaMolNumSpecies[kspec] = 0.0;
826  if (DEBUG_MODE_ENABLED) {
827  sprintf(ANOTE, "dx set to 0, DG flipped sign due to "
828  "changed initial point");
829  }
830  }
831  /*
832  * Form a tentative value of the new species moles
833  */
834  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
835 
836  /*
837  * Check for non-positive mole fraction of major species.
838  * If we find one, we branch to a section below. Then,
839  * depending upon the outcome, we branch to sections below,
840  * or we restart the entire iteration.
841  */
842  if (m_molNumSpecies_new[kspec] <= 0.0) {
843  if (DEBUG_MODE_ENABLED) {
844  sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
845  m_molNumSpecies_new[kspec]);
846  }
847  /* ************************************************* */
848  /* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */
849  /* ************************************************* */
850  /*
851  * We are here when a tentative value of a mole fraction
852  * created by a tentative value of M_DELTAMOLNUMSPECIES(*) is negative.
853  * We branch from here depending upon whether this
854  * species is in a single species phase or in
855  * a multispecies phase.
856  */
857  if (!(m_SSPhase[kspec])) {
858  /*
859  * Section for multispecies phases:
860  * - Cut reaction adjustment for positive kmoles of
861  * major species in multispecies phases.
862  * Decrease its concentration by a factor of 10.
863  */
864  dx = -0.9 * m_molNumSpecies_old[kspec];
865  m_deltaMolNumSpecies[kspec] = dx;
866  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
867  } else {
868  /*
869  * Section for single species phases:
870  * Calculate a dx that will wipe out the
871  * moles in the phase.
872  */
873  dx = -m_molNumSpecies_old[kspec];
874  /*
875  * Calculate an update that doesn't create a negative mole
876  * number for a component species. Actually, restrict this
877  * a little more so that the component values can only be
878  * reduced by two 99%,
879  */
880  for (size_t j = 0; j < m_numComponents; ++j) {
881  if (sc_irxn[j] != 0.0) {
882  m_wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx;
883  if (m_wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) {
884  dx = std::max(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]);
885  }
886  } else {
887  m_wx[j] = m_molNumSpecies_old[j];
888  }
889  }
890  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
891  if (m_molNumSpecies_new[kspec] > 0.0) {
892  m_deltaMolNumSpecies[kspec] = dx;
893  if (DEBUG_MODE_ENABLED) {
894  sprintf(ANOTE,
895  "zeroing SS phase created a neg component species "
896  "-> reducing step size instead");
897  }
898  } else {
899  /*
900  * We are going to zero the single species phase.
901  * Set the existence flag
902  */
903  iph = m_phaseID[kspec];
904  Vphase = m_VolPhaseList[iph];
905  if (DEBUG_MODE_ENABLED) {
906  sprintf(ANOTE, "zeroing out SS phase: ");
907  }
908  /*
909  * Change the base mole numbers for the iteration.
910  * We need to do this here, because we have decided
911  * to eliminate the phase in this special section
912  * outside the main loop.
913  */
914  m_molNumSpecies_new[kspec] = 0.0;
915  doPhaseDeleteIph = iph;
916 
917  if (DEBUG_MODE_ENABLED) {
918  doPhaseDeleteKspec = kspec;
919  if (m_debug_print_lvl >= 2) {
920  if (m_speciesStatus[kspec] >= 0) {
921  plogf(" --- SS species changed to zeroedss: ");
922  plogf("%-12s", m_speciesName[kspec].c_str());
923  plogendl();
924  }
925  }
926  }
927  m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS;
928  ++m_numRxnMinorZeroed;
929  allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
930 
931  for (size_t kk = 0; kk < m_numSpeciesTot; kk++) {
932  m_deltaMolNumSpecies[kk] = 0.0;
933  m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk];
934  }
935  m_deltaMolNumSpecies[kspec] = dx;
936  m_molNumSpecies_new[kspec] = 0.0;
937 
938  for (size_t k = 0; k < m_numComponents; ++k) {
939  m_deltaMolNumSpecies[k] = 0.0;
940  }
941  for (iph = 0; iph < m_numPhases; iph++) {
942  m_deltaPhaseMoles[iph] = 0.0;
943  }
944  }
945  }
946  }
947 
948 #ifdef VCS_LINE_SEARCH
949  /*********************************************************************/
950  /*** LINE SEARCH ALGORITHM FOR MAJOR SPECIES IN NON-IDEAL PHASES *****/
951  /*********************************************************************/
952  /*
953  * Skip the line search if we are birthing a species
954  */
955  if ((dx != 0.0) &&
956  (m_molNumSpecies_old[kspec] > 0.0) &&
957  (doPhaseDeleteIph == -1) &&
958  (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
959  double dx_old = dx;
960 
961  dx = vcs_line_search(irxn, dx_old, ANOTE);
962  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
963  }
964  m_deltaMolNumSpecies[kspec] = dx;
965 #endif
966  } /* End of Loop on ic[irxn] -> the type of species */
967  /***********************************************************************/
968  /****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/
969  /***********************************************************************/
970  if (dx != 0.0 && (m_speciesUnknownType[kspec] !=
971  VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
972  /*
973  * Change the amount of the component compounds according
974  * to the reaction delta that we just computed.
975  * This should keep the amount of material constant.
976  */
977  AssertThrowMsg(fabs(m_deltaMolNumSpecies[kspec] -dx) <
978  1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32),
979  "VCS_SOLVE::solve_tp_inner",
980  "ds[kspec] = " + fp2str(m_deltaMolNumSpecies[kspec]) +
981  " dx = " + fp2str(dx) + " , kspec = " + int2str(kspec) +
982  "\nwe have a problem!");
983  for (size_t k = 0; k < m_numComponents; ++k) {
984  m_deltaMolNumSpecies[k] += sc_irxn[k] * dx;
985  }
986  /*
987  * Calculate the tentative change in the total number of
988  * moles in all of the phases
989  */
990  for (iph = 0; iph < m_numPhases; iph++) {
991  m_deltaPhaseMoles[iph] += dx * m_deltaMolNumPhase(iph,irxn);
992  }
993  }
994 
995  if (DEBUG_MODE_ENABLED) {
996  checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies),
997  VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1);
998  }
999  /*
1000  * Branch point for returning -
1001  */
1002 
1003  if (m_debug_print_lvl >= 2) {
1004  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
1005  plogf(" --- ");
1006  plogf("%-12.12s", m_speciesName[kspec].c_str());
1007  plogf("%3d%11.4E%11.4E%11.4E | %s",
1008  m_speciesStatus[kspec], m_molNumSpecies_old[kspec],
1009  m_molNumSpecies_new[kspec],
1010  m_deltaMolNumSpecies[kspec], ANOTE);
1011  plogendl();
1012  }
1013 
1014  if (doPhaseDeleteIph != npos) {
1015  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1016  plogf(" --- ");
1017  plogf("%-12.12s Main Loop Special Case deleting phase with species: ",
1018  m_speciesName[doPhaseDeleteKspec].c_str());
1019  plogendl();
1020  }
1021  break;
1022  }
1023  } /**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/
1024  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1025  for (size_t k = 0; k < m_numComponents; k++) {
1026  plogf(" --- ");
1027  plogf("%-12.12s", m_speciesName[k].c_str());
1028  plogf(" c%11.4E%11.4E%11.4E |\n",
1029  m_molNumSpecies_old[k],
1030  m_molNumSpecies_old[k]+m_deltaMolNumSpecies[k], m_deltaMolNumSpecies[k]);
1031  }
1032  plogf(" ");
1033  writeline('-', 80);
1034  plogf(" --- Finished Main Loop");
1035  plogendl();
1036  }
1037 
1038  /*************************************************************************/
1039  /*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/
1040  /*************************************************************************/
1041  /*
1042  * We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria
1043  * to limit its magnitude.
1044  */
1045  double par = 0.5;
1046  size_t ll; // only used in DEBUG_MODE
1047  for (size_t k = 0; k < m_numComponents; ++k) {
1048  if (m_molNumSpecies_old[k] > 0.0) {
1049  double xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k];
1050  if (par < xx) {
1051  par = xx;
1052  ll = k;
1053  }
1054  } else if (m_deltaMolNumSpecies[k] < 0.0) {
1055  /*
1056  * If we are here, we then do a step which violates element
1057  * conservation.
1058  */
1059  size_t iph = m_phaseID[k];
1060  m_deltaPhaseMoles[iph] -= m_deltaMolNumSpecies[k];
1061  m_deltaMolNumSpecies[k] = 0.0;
1062  }
1063  }
1064  par = 1.0 / par;
1065  if (par <= 1.01 && par > 0.0) {
1066  /* Reduce the size of the step by the multiplicative factor, par */
1067  par *= 0.99;
1068  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1069  plogf(" --- Reduction in step size due to component ");
1070  plogf("%s", m_speciesName[ll].c_str());
1071  plogf(" going negative = %11.3E", par);
1072  plogendl();
1073  }
1074  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
1075  m_deltaMolNumSpecies[i] *= par;
1076  }
1077  for (size_t iph = 0; iph < m_numPhases; iph++) {
1078  m_deltaPhaseMoles[iph] *= par;
1079  }
1080  } else {
1081  par = 1.0;
1082  }
1083  if (DEBUG_MODE_ENABLED) {
1084  checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies),
1085  VCS_DATA_PTR(m_deltaPhaseMoles), m_numSpeciesTot);
1086  }
1087 
1088  /*
1089  * Now adjust the wt[kspec]'s so that the reflect the decrease in
1090  * the overall length of m_deltaMolNumSpecies[kspec] just calculated. At the end
1091  * of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and tPhMoles1 should all be
1092  * consistent with a new estimate of the state of the system.
1093  */
1094  for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
1095  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec];
1096  if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec]
1097  != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
1098  throw CanteraError("VCS_SOLVE::solve_tp_inner",
1099  "vcs_solve_TP: ERROR on step change wt[" + int2str(kspec) + ":" +
1100  m_speciesName[kspec] + "]: " +
1101  fp2str(m_molNumSpecies_new[kspec]) + " < 0.0");
1102  }
1103  }
1104 
1105  /*
1106  * Calculate the tentative total mole numbers for each phase
1107  */
1108  for (size_t iph = 0; iph < m_numPhases; iph++) {
1109  m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + m_deltaPhaseMoles[iph];
1110  }
1111 
1112  /*
1113  * Set the flags indicating the mole numbers in the vcs_VolPhase
1114  * objects are out of date.
1115  */
1116  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
1117 
1118  /*
1119  * Calculate the new chemical potentials using the tentative
1120  * solution values. We only calculate a subset of these, because
1121  * we have only updated a subset of the W().
1122  */
1123  vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot);
1124 
1125  /*
1126  * Evaluate DeltaG for all components if ITI=0, and for
1127  * major components only if ITI NE 0
1128  */
1129  vcs_deltag(0, false, VCS_STATECALC_NEW);
1130 
1131  /* *************************************************************** */
1132  /* **** CONVERGENCE FORCER SECTION ******************************* */
1133  /* *************************************************************** */
1134  if (printDetails) {
1135  plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
1136  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
1137  VCS_DATA_PTR(m_tPhaseMoles_old)));
1138  plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E",
1139  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
1140  VCS_DATA_PTR(m_tPhaseMoles_new)));
1141  plogendl();
1142  }
1143 
1144  bool forced = vcs_globStepDamp();
1145 
1146  /*
1147  * Print out the changes to the solution that FORCER produced
1148  */
1149  if (printDetails && forced) {
1150 
1151  plogf(" -----------------------------------------------------\n");
1152  plogf(" --- FORCER SUBROUTINE changed the solution:\n");
1153  plogf(" --- SPECIES Status INIT MOLES TENT_MOLES");
1154  plogf(" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n");
1155  for (size_t i = 0; i < m_numComponents; ++i) {
1156  plogf(" --- %-12.12s", m_speciesName[i].c_str());
1157  plogf(" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i],
1158  m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]);
1159  }
1160  for (size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
1161  size_t irxn = kspec - m_numComponents;
1162  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
1163  plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_speciesStatus[kspec],
1164  m_molNumSpecies_old[kspec],
1165  m_molNumSpecies_old[kspec]+m_deltaMolNumSpecies[kspec],
1166  m_molNumSpecies_new[kspec], m_deltaGRxn_old[irxn],
1167  m_deltaGRxn_tmp[irxn], m_deltaGRxn_new[irxn]);
1168  }
1169  writeline(' ', 26, false);
1170  plogf("Norms of Delta G():%14.6E%14.6E\n",
1171  l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)),
1172  l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
1173  plogf(" Total kmoles of gas = %15.7E\n", m_tPhaseMoles_old[0]);
1174  if ((m_numPhases > 1) && (!(m_VolPhaseList[1])->m_singleSpecies)) {
1175  plogf(" Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]);
1176  } else {
1177  plogf(" Total kmoles of liquid = %15.7E\n", 0.0);
1178  }
1179  plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n",
1180  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
1181  VCS_DATA_PTR(m_tPhaseMoles_new)));
1182  plogf(" -----------------------------------------------------");
1183  plogendl();
1184  }
1185 
1186  }
1187  /* *************************************************************** */
1188  /* **** ITERATION SUMMARY PRINTOUT SECTION *********************** */
1189  /* *************************************************************** */
1190 
1191  if (printDetails) {
1192  plogf(" ");
1193  writeline('-', 103);
1194  plogf(" --- Summary of the Update ");
1195  if (iti == 0) {
1196  plogf(" (all species):");
1197  } else {
1198  plogf(" (only major species):");
1199  }
1200  if (m_totalMoleScale != 1.0) {
1201  plogf(" (Total Mole Scale = %g)", m_totalMoleScale);
1202  }
1203  plogf("\n");
1204  plogf(" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT");
1205  plogf(" Mu/RT Init_Del_G/RT Delta_G/RT\n");
1206  for (size_t i = 0; i < m_numComponents; ++i) {
1207  plogf(" --- %-12.12s", m_speciesName[i].c_str());
1208  plogf(" ");
1209  plogf("%14.6E%14.6E%14.6E%14.6E\n", m_molNumSpecies_old[i],
1210  m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i]);
1211  }
1212  for (size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
1213  size_t l1 = i - m_numComponents;
1214  plogf(" --- %-12.12s", m_speciesName[i].c_str());
1215  plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
1216  m_speciesStatus[i], m_molNumSpecies_old[i],
1217  m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_new[i],
1218  m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
1219  }
1220  for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
1221  size_t l1 = kspec - m_numComponents;
1222  plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
1223  plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
1224  m_speciesStatus[kspec], m_molNumSpecies_old[kspec],
1225  m_molNumSpecies_new[kspec], m_feSpecies_old[kspec], m_feSpecies_new[kspec],
1226  m_deltaGRxn_old[l1], m_deltaGRxn_new[l1]);
1227  }
1228  plogf(" ---");
1229  writeline(' ', 56, false);
1230  plogf("Norms of Delta G():%14.6E%14.6E",
1231  l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)),
1232  l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
1233  plogendl();
1234 
1235  plogf(" --- Phase_Name KMoles(after update)\n");
1236  plogf(" --- ");
1237  writeline('-', 50);
1238  for (size_t iph = 0; iph < m_numPhases; iph++) {
1239  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
1240  plogf(" --- %18s = %15.7E\n", Vphase->PhaseName.c_str(), m_tPhaseMoles_new[iph]);
1241  }
1242  plogf(" ");
1243  writeline('-', 103);
1244  plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
1245  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
1246  VCS_DATA_PTR(m_tPhaseMoles_old)));
1247  plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E",
1248  vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
1249  VCS_DATA_PTR(m_tPhaseMoles_new)));
1250  plogendl();
1251  if (DEBUG_MODE_ENABLED && m_VCount->Its > 550) {
1252  plogf(" --- Troublesome solve");
1253  plogendl();
1254  }
1255  }
1256  /*************************************************************************/
1257  /******************* RESET VALUES AT END OF ITERATION ********************/
1258  /******************* UPDATE MOLE NUMBERS *********************************/
1259  /*************************************************************************/
1260  /*
1261  * If the solution wasn't changed in the forcer routine,
1262  * then copy the tentative mole numbers and Phase moles
1263  * into the actual mole numbers and phase moles.
1264  * We will consider this current step to be completed.
1265  *
1266  * Accept the step. -> the tentative solution now becomes
1267  * the real solution. If FORCED is true, then
1268  * we have already done this inside the FORCED
1269  * loop.
1270  */
1271  vcs_updateMolNumVolPhases(VCS_STATECALC_NEW);
1272  m_tPhaseMoles_old = m_tPhaseMoles_new;
1273  m_molNumSpecies_old = m_molNumSpecies_new;
1274  m_actCoeffSpecies_old = m_actCoeffSpecies_new;
1275  m_deltaGRxn_old = m_deltaGRxn_new;
1276  m_feSpecies_old = m_feSpecies_new;
1277 
1278  vcs_setFlagsVolPhases(true, VCS_STATECALC_OLD);
1279  /*
1280  * Increment the iteration counters
1281  */
1282  ++(m_VCount->Its);
1283  ++it1;
1284  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1285  plogf(" --- Increment counter increased, step is accepted: %4d",
1286  m_VCount->Its);
1287  plogendl();
1288  }
1289  /*************************************************************************/
1290  /******************* HANDLE DELETION OF MULTISPECIES PHASES **************/
1291  /*************************************************************************/
1292  /*
1293  * We delete multiphases, when the total moles in the multiphase
1294  * is reduced below a relative threshold.
1295  * Set microscopic multispecies phases with total relative
1296  * number of moles less than VCS_DELETE_PHASE_CUTOFF to
1297  * absolute zero.
1298  */
1299  bool justDeletedMultiPhase = false;
1300  for (size_t iph = 0; iph < m_numPhases; iph++) {
1301  if (!m_VolPhaseList[iph]->m_singleSpecies && m_tPhaseMoles_old[iph] != 0.0 &&
1302  m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) {
1303  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 1) {
1304  plogf(" --- Setting microscopic phase %d to zero", iph);
1305  plogendl();
1306  }
1307  justDeletedMultiPhase = true;
1308  vcs_delete_multiphase(iph);
1309  }
1310  }
1311  /*
1312  * If we have deleted a multispecies phase because the
1313  * equilibrium moles decreased, then we will update all
1314  * the component basis calculation, and therefore all
1315  * of the thermo functions just to be safe.
1316  */
1317  if (justDeletedMultiPhase) {
1318  bool usedZeroedSpecies;
1319  double test = -1.0e-10;
1320  int retn = vcs_basopt(false, VCS_DATA_PTR(m_aw), VCS_DATA_PTR(m_sa),
1321  VCS_DATA_PTR(m_sm), VCS_DATA_PTR(m_ss), test,
1322  &usedZeroedSpecies);
1323  if (retn != VCS_SUCCESS) {
1324  throw CanteraError("VCS_SOLVE::solve_tp_inner",
1325  "BASOPT returned with an error condition");
1326  }
1327  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
1328  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
1329 
1330  vcs_deltag(0, true, VCS_STATECALC_OLD);
1331  iti = 0;
1332  return;
1333  }
1334  /*************************************************************************/
1335  /***************** CHECK FOR ELEMENT ABUNDANCE****************************/
1336  /*************************************************************************/
1337  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1338  plogf(" --- Normal element abundance check");
1339  }
1340  vcs_elab();
1341  if (! vcs_elabcheck(0)) {
1342  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1343  plogf(" - failed -> redoing element abundances.");
1344  plogendl();
1345  }
1346  vcs_elcorr(VCS_DATA_PTR(m_sm), VCS_DATA_PTR(m_wx));
1347  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
1348  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
1349  vcs_deltag(0, true, VCS_STATECALC_OLD);
1350  uptodate_minors = true;
1351  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1352  plogf(" - passed");
1353  plogendl();
1354  }
1355  /*************************************************************************/
1356  /***************** CHECK FOR OPTIMUM BASIS *******************************/
1357  /*************************************************************************/
1358 
1359  for (size_t i = 0; i < m_numRxnRdc; ++i) {
1360  size_t l = m_indexRxnToSpecies[i];
1361  if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
1362  continue;
1363  }
1364  for (size_t j = 0; j < m_numComponents; ++j) {
1365  bool doSwap = false;
1366  if (m_SSPhase[j]) {
1367  doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
1368  (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
1369  if (!m_SSPhase[l] && doSwap) {
1370  doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
1371  }
1372  } else {
1373  if (m_SSPhase[l]) {
1374  doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
1375  (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
1376  if (!doSwap) {
1377  doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
1378  }
1379  } else {
1380  doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) >
1381  (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
1382  }
1383  }
1384  if (doSwap && m_stoichCoeffRxnMatrix(j,i) != 0.0) {
1385  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1386  plogf(" --- Get a new basis because ");
1387  plogf("%s", m_speciesName[l].c_str());
1388  plogf(" is better than comp ");
1389  plogf("%s", m_speciesName[j].c_str());
1390  plogf(" and share nonzero stoic: %-9.1f",
1391  m_stoichCoeffRxnMatrix(j,i));
1392  plogendl();
1393  }
1394  forceComponentCalc = 1;
1395  return;
1396  }
1397 #ifdef DEBUG_NOT
1398  if (m_speciesStatus[l] == VCS_SPECIES_ZEROEDMS && m_molNumSpecies_old[j] == 0.0 && m_stoichCoeffRxnMatrix(j,i) != 0.0 && dg[i] < 0.0) {
1399  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1400  plogf(" --- Get a new basis because %s", m_speciesName[l].c_str());
1401  plogf(" has dg < 0.0 and comp %s has zero mole num",
1402  m_speciesName[j].c_str());
1403  plogf(" and share nonzero stoic: %-9.1f",
1404  m_stoichCoeffRxnMatrix(j,i));
1405  plogendl();
1406  }
1407  return;
1408  }
1409 #endif
1410  }
1411  }
1412  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1413  plogf(" --- Check for an optimum basis passed");
1414  plogendl();
1415  }
1416  stage = EQUILIB_CHECK;
1417  /*************************************************************************/
1418  /********************** RE-EVALUATE MAJOR-MINOR VECTOR IF NECESSARY ******/
1419  /*************************************************************************/
1420  /*
1421  * Skip this section if we haven't done a full calculation.
1422  * Go right to the check equilibrium section
1423  */
1424  if (iti != 0) {
1425  return;
1426  }
1427  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1428  plogf(" --- Reevaluate major-minor status of noncomponents:\n");
1429  }
1430  m_numRxnMinorZeroed = 0;
1431  for (size_t irxn = 0; irxn < m_numRxnRdc; irxn++) {
1432  size_t kspec = m_indexRxnToSpecies[irxn];
1433  int speciesType = vcs_species_type(kspec);
1434  if (speciesType < VCS_SPECIES_MINOR) {
1435  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1436  if (m_speciesStatus[kspec] >= VCS_SPECIES_MINOR) {
1437  plogf(" --- major/minor species is now zeroed out: %s\n",
1438  m_speciesName[kspec].c_str());
1439  }
1440  }
1441  ++m_numRxnMinorZeroed;
1442  } else if (speciesType == VCS_SPECIES_MINOR) {
1443  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1444  if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
1445  if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) {
1446  plogf(" --- Noncomponent turned from major to minor: ");
1447  } else if (kspec < m_numComponents) {
1448  plogf(" --- Component turned into a minor species: ");
1449  } else {
1450  plogf(" --- Zeroed Species turned into a "
1451  "minor species: ");
1452  }
1453  plogf("%s\n", m_speciesName[kspec].c_str());
1454  }
1455  }
1456  ++m_numRxnMinorZeroed;
1457  } else if (speciesType == VCS_SPECIES_MAJOR) {
1458  if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
1459  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1460  if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
1461  plogf(" --- Noncomponent turned from minor to major: ");
1462  } else if (kspec < m_numComponents) {
1463  plogf(" --- Component turned into a major: ");
1464  } else {
1465  plogf(" --- Noncomponent turned from zeroed to major: ");
1466  }
1467  plogf("%s\n", m_speciesName[kspec].c_str());
1468  }
1469  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
1470  }
1471  }
1472  m_speciesStatus[kspec] = speciesType;
1473  }
1474  /*
1475  * This logical variable indicates whether all current
1476  * non-component species are minor or nonexistent
1477  */
1478  allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc);
1479 }
1480 
1481 void VCS_SOLVE::solve_tp_equilib_check(bool& allMinorZeroedSpecies,
1482  bool& uptodate_minors,
1483  bool& giveUpOnElemAbund, int& solveFail,
1484  size_t& iti, size_t& it1, int maxit,
1485  int& stage, bool& lec)
1486 {
1487  if (! allMinorZeroedSpecies) {
1488  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1489  plogf(" --- Equilibrium check for major species: ");
1490  }
1491  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
1492  size_t kspec = irxn + m_numComponents;
1493  if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmaj)) {
1494  if (m_VCount->Its >= maxit) {
1495  solveFail = -1;
1496  /*
1497  * Clean up and exit code even though we haven't
1498  * converged. -> we have run out of iterations!
1499  */
1500  stage = RETURN_A;
1501  } else {
1502  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1503  plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str());
1504  }
1505  // Convergence amongst major species has not been achieved
1506  /*
1507  * Go back and do another iteration with variable ITI
1508  */
1509  iti = ((it1/4) *4) - it1;
1510  stage = MAIN;
1511  }
1512  return;
1513  }
1514  }
1515  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1516  plogf(" MAJOR SPECIES CONVERGENCE achieved");
1517  plogendl();
1518  }
1519  }
1520  else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1521  plogf(" MAJOR SPECIES CONVERGENCE achieved "
1522  "(because there are no major species)");
1523  plogendl();
1524  }
1525  // Convergence amongst major species has been achieved
1526 
1527  /*************************************************************************/
1528  /*************** EQUILIBRIUM CHECK FOR MINOR SPECIES *********************/
1529  /*************************************************************************/
1530 
1531  if (m_numRxnMinorZeroed != 0) {
1532  /*
1533  * Calculate the chemical potential and reaction DeltaG
1534  * for minor species, if needed.
1535  */
1536  if (iti != 0) {
1537  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
1538  vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
1539  vcs_deltag(1, false, VCS_STATECALC_OLD);
1540  uptodate_minors = true;
1541  }
1542  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1543  plogf(" --- Equilibrium check for minor species: ");
1544  }
1545  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
1546  size_t kspec = irxn + m_numComponents;
1547  if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR && (fabs(m_deltaGRxn_new[irxn]) > m_tolmin)) {
1548  if (m_VCount->Its >= maxit) {
1549  solveFail = -1;
1550  /*
1551  * Clean up and exit code. -> Even though we have not
1552  * converged, we have run out of iterations !
1553  */
1554  stage = RETURN_A;
1555  return;
1556  }
1557  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1558  plogf("%s failed\n", m_speciesName[m_indexRxnToSpecies[irxn]].c_str());
1559  }
1560  /*
1561  * Set iti to zero to force a full calculation, and go back
1562  * to the main loop to do another iteration.
1563  */
1564  iti = 0;
1565  stage = MAIN;
1566  return;
1567  }
1568  }
1569  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1570  plogf(" CONVERGENCE achieved\n");
1571  }
1572  }
1573  /*************************************************************************/
1574  /*********************** FINAL ELEMENTAL ABUNDANCE CHECK *****************/
1575  /*************************************************************************/
1576  /*
1577  * Recalculate the element abundance vector again
1578  */
1579  vcs_updateVP(VCS_STATECALC_OLD);
1580  vcs_elab();
1581 
1582  /* LEC is only true when we are near the end game */
1583  if (lec) {
1584  if (!giveUpOnElemAbund) {
1585  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1586  plogf(" --- Check the Full Element Abundances: ");
1587  }
1588  /*
1589  * Final element abundance check:
1590  * If we fail then we need to go back and correct
1591  * the element abundances, and then go do a major step
1592  */
1593  if (! vcs_elabcheck(1)) {
1594  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1595  if (! vcs_elabcheck(0)) {
1596  plogf(" failed\n");
1597  } else {
1598  plogf(" passed for NC but failed for NE: RANGE ERROR\n");
1599  }
1600  }
1601  // delete?
1602  stage = ELEM_ABUND_CHECK;
1603  return;
1604  }
1605  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1606  plogf(" passed\n");
1607  }
1608  }
1609  /*
1610  * If we have deleted a species then we need to recheck the
1611  * the deleted species, before exiting
1612  */
1613  if (m_numSpeciesRdc != m_numSpeciesTot) {
1614  stage = RECHECK_DELETED;
1615  return;
1616  }
1617  /* - Final checks are passed -> go check out */
1618  stage = RETURN_A;
1619  }
1620  lec = true;
1621 }
1622 
1623 void VCS_SOLVE::solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec,
1624  bool& giveUpOnElemAbund,
1625  int& finalElemAbundAttempts,
1626  int& rangeErrorFound)
1627 {
1628  /*
1629  * HKM - Put in an element abundance check. The element abundances
1630  * were being corrected even if they were perfectly OK to
1631  * start with. This is actually an expensive operation, so
1632  * I took it out. Also vcs_dfe() doesn't need to be called if
1633  * no changes were made.
1634  */
1635  rangeErrorFound = 0;
1636  if (! vcs_elabcheck(1)) {
1637  bool ncBefore = vcs_elabcheck(0);
1638  vcs_elcorr(VCS_DATA_PTR(m_sm), VCS_DATA_PTR(m_wx));
1639  bool ncAfter = vcs_elabcheck(0);
1640  bool neAfter = vcs_elabcheck(1);
1641  /*
1642  * Go back to evaluate the total moles of gas and liquid.
1643  */
1644  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
1645  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
1646  vcs_deltag(0, false, VCS_STATECALC_OLD);
1647  if (!ncBefore) {
1648  if (ncAfter) {
1649  /*
1650  * We have breathed new life into the old problem. Now the
1651  * element abundances up to NC agree. Go back and
1652  * restart the main loop calculation, resetting the
1653  * end conditions.
1654  */
1655  lec = false;
1656  iti = 0;
1657  stage = MAIN;
1658  } else {
1659  /*
1660  * We are still hosed
1661  */
1662  if (finalElemAbundAttempts >= 3) {
1663  giveUpOnElemAbund = true;
1664  stage = EQUILIB_CHECK;
1665  } else {
1666  finalElemAbundAttempts++;
1667  lec = false;
1668  iti = 0;
1669  stage = MAIN;
1670  }
1671  }
1672  return;
1673  } else if (ncAfter) {
1674  if (!neAfter) {
1675  /*
1676  * Probably an unrecoverable range error
1677  */
1678  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1679  plogf(" --- vcs_solve_tp: RANGE SPACE ERROR ENCOUNTERED\n");
1680  plogf(" --- vcs_solve_tp: - Giving up on NE Element Abundance satisfaction \n");
1681  plogf(" --- vcs_solve_tp: - However, NC Element Abundance criteria is satisfied \n");
1682  plogf(" --- vcs_solve_tp: - Returning the calculated equilibrium condition ");
1683  plogendl();
1684  }
1685  rangeErrorFound = 1;
1686  giveUpOnElemAbund = true;
1687  }
1688  /*
1689  * Recovery of end element abundances
1690  * -> go do equilibrium check again and then
1691  * check out.
1692  */
1693  stage = EQUILIB_CHECK;
1694  return;
1695  }
1696  }
1697  // Calculate delta g's
1698  vcs_deltag(0, false, VCS_STATECALC_OLD);
1699  // Go back to equilibrium check as a prep to eventually checking out
1700  stage = EQUILIB_CHECK;
1701 }
1702 
1703 double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
1704  char* ANOTE) const
1705 {
1706  double dx = 0.0, a;
1707  double w_kspec = m_molNumSpecies_old[kspec];
1708  double molNum_kspec_new;
1709  double wTrial, tmp;
1710  double dg_irxn = m_deltaGRxn_old[irxn];
1711  doublereal s;
1712  size_t iph = m_phaseID[kspec];
1713 
1714  *do_delete = false;
1715  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
1716  if (w_kspec <= 0.0) {
1717  w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
1718  }
1719  dg_irxn = std::max(dg_irxn, -200.0);
1720  if (DEBUG_MODE_ENABLED && ANOTE) {
1721  sprintf(ANOTE,"minor species alternative calc");
1722  }
1723  if (dg_irxn >= 23.0) {
1724  molNum_kspec_new = w_kspec * 1.0e-10;
1725  if (w_kspec < VCS_DELETE_MINORSPECIES_CUTOFF) {
1726  goto L_ZERO_SPECIES;
1727  }
1728  return molNum_kspec_new - w_kspec;
1729  } else {
1730  if (fabs(dg_irxn) <= m_tolmin2) {
1731  molNum_kspec_new = w_kspec;
1732  return 0.0;
1733  }
1734  }
1735 
1736  /*
1737  * get the diagonal of the activity coefficient Jacobian
1738  */
1739  s = m_np_dLnActCoeffdMolNum(kspec,kspec) / (m_tPhaseMoles_old[iph]);
1740  /*
1741  * We fit it to a power law approximation of the activity coefficient
1742  *
1743  * gamma = gamma_0 * ( x / x0)**a
1744  *
1745  * where a is forced to be a little bit greater than -1.
1746  * We do this so that the resulting expression is always nonnegative
1747  *
1748  * We then solve the resulting calculation:
1749  *
1750  * gamma * x = gamma_0 * x0 exp (-deltaG/RT);
1751  *
1752  *
1753  */
1754  a = clip(w_kspec * s, -1.0+1e-8, 100.0);
1755  tmp = clip(-dg_irxn / (1.0 + a), -200.0, 200.0);
1756  wTrial = w_kspec * exp(tmp);
1757 
1758  molNum_kspec_new = wTrial;
1759 
1760  if (wTrial > 100. * w_kspec) {
1761  double molNumMax = 0.0001 * m_tPhaseMoles_old[iph];
1762  if (molNumMax < 100. * w_kspec) {
1763  molNumMax = 100. * w_kspec;
1764  }
1765  if (wTrial > molNumMax) {
1766  molNum_kspec_new = molNumMax;
1767  } else {
1768  molNum_kspec_new = wTrial;
1769  }
1770 
1771  } else if (1.0E10 * wTrial < w_kspec) {
1772  molNum_kspec_new= 1.0E-10 * w_kspec;
1773  } else {
1774  molNum_kspec_new = wTrial;
1775  }
1776 
1777 
1778  if ((molNum_kspec_new) < VCS_DELETE_MINORSPECIES_CUTOFF) {
1779  goto L_ZERO_SPECIES;
1780  }
1781  return molNum_kspec_new - w_kspec;
1782  /*
1783  *
1784  * Alternate return based for cases where we need to delete the species
1785  * from the current list of active species, because its concentration
1786  * has gotten too small.
1787  */
1788 L_ZERO_SPECIES:
1789  ;
1790  *do_delete = true;
1791  return - w_kspec;
1792  } else {
1793  /*
1794  * Voltage calculation
1795  * Need to check the sign -> This is good for electrons
1796  */
1797  dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim;
1798  if (DEBUG_MODE_ENABLED && ANOTE) {
1799  sprintf(ANOTE,"voltage species alternative calc");
1800  }
1801  }
1802  return dx;
1803 }
1804 
1805 int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr)
1806 {
1807  size_t irxn = kspec - m_numComponents;
1808  int retn = 1;
1809  double delta = *delta_ptr;
1810  AssertThrowMsg(kspec >= m_numComponents, "VCS_SOLVE::delta_species",
1811  "delete_species() ERROR: called for a component " + int2str(kspec));
1812  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
1813  /*
1814  * Attempt the given dx. If it doesn't work, try to see if a smaller
1815  * one would work,
1816  */
1817  double dx = delta;
1818  double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
1819  for (size_t j = 0; j < m_numComponents; ++j) {
1820  if (m_molNumSpecies_old[j] > 0.0) {
1821  double tmp = sc_irxn[j] * dx;
1822  if (-tmp > m_molNumSpecies_old[j]) {
1823  retn = 0;
1824  dx = std::min(dx, - m_molNumSpecies_old[j] / sc_irxn[j]);
1825  }
1826  }
1827  /*
1828  * If the component has a zero concentration and is a reactant
1829  * in the formation reaction, then dx == 0.0, and we just return.
1830  */
1831  if (m_molNumSpecies_old[j] <= 0.0) {
1832  if (sc_irxn[j] < 0.0) {
1833  *delta_ptr = 0.0;
1834  return 0;
1835  }
1836  }
1837  }
1838  /*
1839  * ok, we found a positive dx. implement it.
1840  */
1841  *delta_ptr = dx;
1842  m_molNumSpecies_old[kspec] += dx;
1843  size_t iph = m_phaseID[kspec];
1844  m_tPhaseMoles_old[iph] += dx;
1845  vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD);
1846 
1847  for (size_t j = 0; j < m_numComponents; ++j) {
1848  double tmp = sc_irxn[j] * dx;
1849  if (tmp != 0.0) {
1850  iph = m_phaseID[j];
1851  m_molNumSpecies_old[j] += tmp;
1852  m_tPhaseMoles_old[iph] += tmp;
1853  vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD);
1854  m_molNumSpecies_old[j] = std::max(m_molNumSpecies_old[j], 0.0);
1855  }
1856  }
1857  }
1858  return retn;
1859 }
1860 
1861 int VCS_SOLVE::vcs_zero_species(const size_t kspec)
1862 {
1863  int retn = 1;
1864  /*
1865  * Calculate a delta that will eliminate the species.
1866  */
1867  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
1868  double dx = -(m_molNumSpecies_old[kspec]);
1869  if (dx != 0.0) {
1870  retn = delta_species(kspec, &dx);
1871  if (DEBUG_MODE_ENABLED && !retn && m_debug_print_lvl >= 1) {
1872  plogf("vcs_zero_species: Couldn't zero the species %d, "
1873  "did delta of %g. orig conc of %g",
1874  kspec, dx, m_molNumSpecies_old[kspec] + dx);
1875  plogendl();
1876  }
1877  }
1878  }
1879  return retn;
1880 }
1881 
1882 int VCS_SOLVE::vcs_delete_species(const size_t kspec)
1883 {
1884  const size_t klast = m_numSpeciesRdc - 1;
1885  const size_t iph = m_phaseID[kspec];
1886  vcs_VolPhase* const Vphase = m_VolPhaseList[iph];
1887  const size_t irxn = kspec - m_numComponents;
1888  /*
1889  * Zero the concentration of the species.
1890  * -> This zeroes w[kspec] and modifies m_tPhaseMoles_old[]
1891  */
1892  const int retn = vcs_zero_species(kspec);
1893  if (!retn) {
1894  throw CanteraError("VCS_SOLVE::vcs_delete_species",
1895  "Failed to delete a species!");
1896  }
1897  /*
1898  * Decrement the minor species counter if the current species is
1899  * a minor species
1900  */
1901  if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
1902  --(m_numRxnMinorZeroed);
1903  }
1904  m_speciesStatus[kspec] = VCS_SPECIES_DELETED;
1905  m_deltaGRxn_new[irxn] = 0.0;
1906  m_deltaGRxn_old[irxn] = 0.0;
1907  m_feSpecies_new[kspec] = 0.0;
1908  m_feSpecies_old[kspec] = 0.0;
1909  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
1910  /*
1911  * Rearrange the data if the current species isn't the last active
1912  * species.
1913  */
1914  if (kspec != klast) {
1915  vcs_switch_pos(true, klast, kspec);
1916  }
1917  /*
1918  * Adjust the total moles in a phase downwards.
1919  */
1920  Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
1921  VCS_DATA_PTR(m_molNumSpecies_old),
1922  VCS_DATA_PTR(m_tPhaseMoles_old));
1923 
1924  /*
1925  * Adjust the current number of active species and reactions counters
1926  */
1927  --(m_numRxnRdc);
1928  --(m_numSpeciesRdc);
1929 
1930  /*
1931  * Check to see whether we have just annihilated a multispecies phase.
1932  * If it is extinct, call the delete_multiphase() function.
1933  */
1934  if (! m_SSPhase[klast]) {
1935  if (Vphase->exists() != VCS_PHASE_EXIST_ALWAYS) {
1936  bool stillExists = false;
1937  for (size_t k = 0; k < m_numSpeciesRdc; k++) {
1938  if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
1939  if (m_phaseID[k] == iph) {
1940  if (m_molNumSpecies_old[k] > 0.0) {
1941  stillExists = true;
1942  break;
1943  }
1944  }
1945  }
1946  }
1947  if (!stillExists) {
1948  vcs_delete_multiphase(iph);
1949  }
1950  }
1951  }
1952  /*
1953  * When the total number of noncomponent species is zero, we
1954  * have to signal the calling code
1955  */
1956  return (m_numRxnRdc == 0);
1957 }
1958 
1959 void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec)
1960 {
1961  size_t iph = m_phaseID[kspec];
1962  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
1963  plogf(" --- Add back a deleted species: %-12s\n", m_speciesName[kspec].c_str());
1964  }
1965  /*
1966  * Set the species back to minor species status
1967  * this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[]
1968  * HKM -> make this a relative mole number!
1969  */
1970  double dx = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10.;
1971  delta_species(kspec, &dx);
1972  m_speciesStatus[kspec] = VCS_SPECIES_MINOR;
1973 
1974  if (m_SSPhase[kspec]) {
1975  m_speciesStatus[kspec] = VCS_SPECIES_MAJOR;
1976  --(m_numRxnMinorZeroed);
1977  }
1978 
1979  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
1980  Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
1981  VCS_DATA_PTR(m_molNumSpecies_old),
1982  VCS_DATA_PTR(m_tPhaseMoles_old));
1983  /*
1984  * We may have popped a multispecies phase back
1985  * into existence. If we did, we have to check
1986  * the other species in that phase.
1987  * Take care of the m_speciesStatus[] flag.
1988  * The value of m_speciesStatus[] must change from
1989  * VCS_SPECIES_ZEROEDPHASE to VCS_SPECIES_ZEROEDMS
1990  * for those other species.
1991  */
1992  if (! m_SSPhase[kspec]) {
1993  if (Vphase->exists() == VCS_PHASE_EXIST_NO) {
1994  Vphase->setExistence(VCS_PHASE_EXIST_YES);
1995  for (size_t k = 0; k < m_numSpeciesTot; k++) {
1996  if (m_phaseID[k] == iph) {
1997  if (m_speciesStatus[k] != VCS_SPECIES_DELETED) {
1998  m_speciesStatus[k] = VCS_SPECIES_MINOR;
1999  }
2000  }
2001  }
2002  }
2003  } else {
2004  Vphase->setExistence(VCS_PHASE_EXIST_YES);
2005  }
2006 
2007  ++(m_numRxnRdc);
2008  ++(m_numSpeciesRdc);
2009  ++(m_numRxnMinorZeroed);
2010 
2011  if (kspec != (m_numSpeciesRdc - 1)) {
2012  /*
2013  * Rearrange both the species and the non-component global data
2014  */
2015  vcs_switch_pos(true, (m_numSpeciesRdc - 1), kspec);
2016  }
2017 }
2018 
2019 bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
2020 {
2021  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
2022  bool successful = true;
2023  /*
2024  * set the phase existence flag to dead
2025  */
2026  Vphase->setTotalMoles(0.0);
2027  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2028  plogf(" --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName.c_str());
2029  }
2030 
2031  /*
2032  * Loop over all of the species in the phase.
2033  */
2034  for (size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
2035  if (m_phaseID[kspec] == iph) {
2036  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2037  /*
2038  * calculate an extent of rxn, dx, that zeroes out the species.
2039  */
2040  double dx = - (m_molNumSpecies_old[kspec]);
2041  double dxTent = dx;
2042 
2043  int retn = delta_species(kspec, &dxTent);
2044  if (retn != 1) {
2045  successful = false;
2046  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2047  plogf(" --- delete_multiphase %d, %s ERROR problems deleting species %s\n",
2048  iph, Vphase->PhaseName.c_str(), m_speciesName[kspec].c_str());
2049  plogf(" --- delta attempted: %g achieved: %g "
2050  " Zeroing it manually\n", dx, dxTent);
2051  }
2052  m_molNumSpecies_old[kspec] = 0.0;
2053  m_molNumSpecies_new[kspec] = 0.0;
2054  m_deltaMolNumSpecies[kspec] = 0.0;
2055  // recover the total phase moles.
2056  vcs_tmoles();
2057  } else {
2058  /*
2059  * Set the mole number of that species to zero.
2060  */
2061  m_molNumSpecies_old[kspec] = 0.0;
2062  m_molNumSpecies_new[kspec] = 0.0;
2063  m_deltaMolNumSpecies[kspec] = 0.0;
2064  }
2065  /*
2066  * Change the status flag of the species to that of an
2067  * zeroed phase
2068  */
2069  m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS;
2070  }
2071  }
2072  }
2073 
2074  double dj, dxWant, dxPerm = 0.0, dxPerm2 = 0.0;
2075  for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
2076  if (m_phaseID[kcomp] == iph) {
2077  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2078  plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n",
2079  kcomp, m_speciesName[kcomp].c_str(), m_molNumSpecies_old[kcomp]);
2080  }
2081  if (m_molNumSpecies_old[kcomp] != 0.0) {
2082  for (size_t kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
2083  size_t irxn = kspec - m_numComponents;
2084  if (m_phaseID[kspec] != iph) {
2085  if (m_stoichCoeffRxnMatrix(kcomp,irxn) != 0.0) {
2086  dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn);
2087  if (dxWant + m_molNumSpecies_old[kspec] < 0.0) {
2088  dxPerm = -m_molNumSpecies_old[kspec];
2089  }
2090  for (size_t jcomp = 0; kcomp < m_numComponents; ++kcomp) {
2091  if (jcomp != kcomp) {
2092  if (m_phaseID[jcomp] == iph) {
2093  dxPerm = 0.0;
2094  } else {
2095  dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn);
2096  if (dj + m_molNumSpecies_old[kcomp] < 0.0) {
2097  dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(jcomp,irxn);
2098  }
2099  if (fabs(dxPerm2) < fabs(dxPerm)) {
2100  dxPerm = dxPerm2;
2101  }
2102  }
2103  }
2104  }
2105  }
2106  if (dxPerm != 0.0) {
2107  delta_species(kspec, &dxPerm);
2108  }
2109  }
2110  }
2111 
2112  }
2113  if (m_molNumSpecies_old[kcomp] != 0.0) {
2114  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2115  plogf(" --- delete_multiphase One of the species is a component %d - %s still with mole number %g\n",
2116  kcomp, m_speciesName[kcomp].c_str(), m_molNumSpecies_old[kcomp]);
2117  plogf(" --- zeroing it \n");
2118  }
2119  m_molNumSpecies_old[kcomp] = 0.0;
2120  }
2121  m_speciesStatus[kcomp] = VCS_SPECIES_ZEROEDMS;
2122  }
2123  }
2124 
2125 
2126  /*
2127  * Loop over all of the inactive species in the phase:
2128  * Right now we reinstate all species in a deleted multiphase.
2129  * We may only want to reinstate the "major ones" in the future.
2130  * Note, species in phases with zero mole numbers are still
2131  * considered active. Whether the phase pops back into
2132  * existence or not is checked as part of the main iteration
2133  * loop.
2134  */
2135  for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
2136  if (m_phaseID[kspec] == iph) {
2137  m_molNumSpecies_old[kspec] = 0.0;
2138  m_molNumSpecies_new[kspec] = 0.0;
2139  m_deltaMolNumSpecies[kspec] = 0.0;
2140  m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDMS;
2141 
2142  ++(m_numRxnRdc);
2143  ++(m_numSpeciesRdc);
2144  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2145  plogf(" --- Make %s", m_speciesName[kspec].c_str());
2146  plogf(" an active but zeroed species because its phase "
2147  "was zeroed\n");
2148  }
2149  if (kspec != (m_numSpeciesRdc - 1)) {
2150  /*
2151  * Rearrange both the species and the non-component global data
2152  */
2153  vcs_switch_pos(true, (m_numSpeciesRdc - 1), kspec);
2154  }
2155  }
2156  }
2157 
2158  /*
2159  * Zero out the total moles counters for the phase
2160  */
2161  m_tPhaseMoles_old[iph] = 0.0;
2162  m_tPhaseMoles_new[iph] = 0.0;
2163  m_deltaPhaseMoles[iph] = 0.0;
2164  /*
2165  * Upload the state to the VP object
2166  */
2167  Vphase->setTotalMoles(0.0);
2168 
2169  return successful;
2170 }
2171 
2172 int VCS_SOLVE::vcs_recheck_deleted()
2173 {
2174  double* xtcutoff = VCS_DATA_PTR(m_TmpPhase);
2175  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2176  plogf(" --- Start rechecking deleted species in multispec phases\n");
2177  }
2178  if (m_numSpeciesRdc == m_numSpeciesTot) {
2179  return 0;
2180  }
2181  /*
2182  * Use the standard chemical potentials for the chemical potentials
2183  * of deleted species. Then, calculate Delta G for
2184  * for formation reactions.
2185  * Note: fe[] here includes everything except for the ln(x[i]) term
2186  */
2187  for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
2188  size_t iph = m_phaseID[kspec];
2189  m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
2190  - m_lnMnaughtSpecies[kspec]
2191  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
2192  }
2193 
2194  /*
2195  * Recalculate the DeltaG's of the formation reactions for the
2196  * deleted species in the mechanism
2197  */
2198  vcs_deltag(0, true, VCS_STATECALC_NEW);
2199 
2200  for (size_t iph = 0; iph < m_numPhases; iph++) {
2201  if (m_tPhaseMoles_old[iph] > 0.0) {
2202  xtcutoff[iph] = log(1.0 / VCS_RELDELETE_SPECIES_CUTOFF);
2203  } else {
2204  xtcutoff[iph] = 0.0;
2205  }
2206  }
2207 
2208  /*
2209  *
2210  * We are checking the equation:
2211  *
2212  * sum_u = sum_j_comp [ sigma_i_j * u_j ]
2213  * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
2214  *
2215  * by first evaluating:
2216  *
2217  * DG_i_O = u_i_O - sum_u.
2218  *
2219  * Then, if TL is zero, the phase pops into existence if DG_i_O < 0.
2220  * Also, if the phase exists, then we check to see if the species
2221  * can have a mole number larger than VCS_DELETE_SPECIES_CUTOFF
2222  * (default value = 1.0E-32).
2223  *
2224  * HKM:
2225  * This seems to be an inconsistency in the algorithm here that needs
2226  * correcting. The requirement above may bypass some multiphases which
2227  * should exist. The real requirement for the phase to exist is:
2228  *
2229  * sum_i_in_phase [ exp(-DG_i_O) ] >= 1.0
2230  *
2231  * Thus, we need to amend the code. Also nonideal solutions will tend to
2232  * complicate matters severely also.
2233  */
2234  int npb = 0;
2235  for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
2236  size_t kspec = m_indexRxnToSpecies[irxn];
2237  size_t iph = m_phaseID[kspec];
2238  if (m_tPhaseMoles_old[iph] == 0.0) {
2239  if (m_deltaGRxn_new[irxn] < 0.0) {
2240  vcs_reinsert_deleted(kspec);
2241  npb++;
2242  } else {
2243  m_molNumSpecies_old[kspec] = 0.0;
2244  }
2245  } else if (m_tPhaseMoles_old[iph] > 0.0) {
2246  if (m_deltaGRxn_new[irxn] < xtcutoff[iph]) {
2247  vcs_reinsert_deleted(kspec);
2248  npb++;
2249  }
2250  }
2251  }
2252  return npb;
2253 }
2254 
2255 bool VCS_SOLVE::recheck_deleted_phase(const int iphase)
2256 {
2257 
2258  // Check first to see if the phase is in fact deleted
2259  const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
2260  if (Vphase->exists() != VCS_PHASE_EXIST_NO) {
2261  return false;
2262  }
2263  if (Vphase->exists() == VCS_PHASE_EXIST_ZEROEDPHASE) {
2264  return false;
2265  }
2266  if (Vphase->m_singleSpecies) {
2267  size_t kspec = Vphase->spGlobalIndexVCS(0);
2268  size_t irxn = kspec + m_numComponents;
2269  if (m_deltaGRxn_old[irxn] < 0.0) {
2270  return true;
2271  }
2272  return false;
2273  }
2274 
2275  double phaseDG = 1.0;
2276  for (size_t kk = 0; kk < Vphase->nSpecies(); kk++) {
2277  size_t kspec = Vphase->spGlobalIndexVCS(kk);
2278  size_t irxn = kspec + m_numComponents;
2279  m_deltaGRxn_old[irxn] = clip(m_deltaGRxn_old[irxn], -50.0, 50.0);
2280  phaseDG -= exp(-m_deltaGRxn_old[irxn]);
2281  }
2282 
2283  if (phaseDG < 0.0) {
2284  return true;
2285  }
2286  return false;
2287 }
2288 
2289 size_t VCS_SOLVE::vcs_add_all_deleted()
2290 {
2291  size_t retn;
2292  if (m_numSpeciesRdc == m_numSpeciesTot) {
2293  return 0;
2294  }
2295  /*
2296  * Use the standard chemical potentials for the chemical potentials of deleted species. Then, calculate Delta G for
2297  * for formation reactions.
2298  * We are relying here on a old saved value of m_actCoeffSpecies_old[kspec]
2299  * being sufficiently good. Note, we will recalculate everything at the end of the routine.
2300  */
2301  m_molNumSpecies_new = m_molNumSpecies_old;
2302 
2303  for (int cits = 0; cits < 3; cits++) {
2304  for (size_t kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
2305  size_t iph = m_phaseID[kspec];
2306  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
2307  if (m_molNumSpecies_new[kspec] == 0.0) {
2308  m_molNumSpecies_new[kspec] = VCS_DELETE_MINORSPECIES_CUTOFF * 1.0E-10;
2309  }
2310  if (!Vphase->m_singleSpecies) {
2311  Vphase->sendToVCS_ActCoeff(VCS_STATECALC_NEW, VCS_DATA_PTR(m_actCoeffSpecies_new));
2312  }
2313  m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_new[kspec]) - m_lnMnaughtSpecies[kspec]
2314  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
2315  }
2316  /*
2317  * Recalculate the DeltaG's of the formation reactions for the deleted species in the mechanism
2318  */
2319  vcs_deltag(0, true, VCS_STATECALC_NEW);
2320  for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
2321  size_t kspec = m_indexRxnToSpecies[irxn];
2322  size_t iph = m_phaseID[kspec];
2323  if (m_tPhaseMoles_old[iph] > 0.0) {
2324  double maxDG = std::min(m_deltaGRxn_new[irxn], 690.0);
2325  double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
2326  m_molNumSpecies_new[kspec] = dx;
2327  }
2328  }
2329  }
2330  for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
2331  size_t kspec = m_indexRxnToSpecies[irxn];
2332  size_t iph = m_phaseID[kspec];
2333  if (m_tPhaseMoles_old[iph] > 0.0) {
2334  double dx = m_molNumSpecies_new[kspec];
2335  retn = delta_species(kspec, &dx);
2336  if (retn == 0) {
2337  if (DEBUG_MODE_ENABLED && m_debug_print_lvl) {
2338  plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
2339  m_speciesName[kspec].c_str(), kspec, dx);
2340  }
2341  if (dx > 1.0E-50) {
2342  dx = 1.0E-50;
2343  retn = delta_species(kspec, &dx);
2344  if (DEBUG_MODE_ENABLED && retn == 0 && m_debug_print_lvl) {
2345  plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
2346  m_speciesName[kspec].c_str(), kspec, dx);
2347  }
2348  }
2349  }
2350  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2351  if (retn != 0) {
2352  plogf(" --- add_deleted(): species %s added back in with mol number %g",
2353  m_speciesName[kspec].c_str(), dx);
2354  plogendl();
2355  } else {
2356  plogf(" --- add_deleted(): species %s failed to be added back in");
2357  plogendl();
2358  }
2359  }
2360  }
2361  }
2362  vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
2363  vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
2364  vcs_deltag(0, true, VCS_STATECALC_OLD);
2365 
2366  retn = 0;
2367  for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
2368  size_t kspec = m_indexRxnToSpecies[irxn];
2369  size_t iph = m_phaseID[kspec];
2370  if (m_tPhaseMoles_old[iph] > 0.0) {
2371  if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) {
2372  if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) >
2373  VCS_DELETE_MINORSPECIES_CUTOFF) ||
2374  (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) {
2375  retn++;
2376  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2377  plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g",
2378  m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
2379  m_deltaGRxn_old[irxn]);
2380  plogendl();
2381  }
2382  }
2383  }
2384  }
2385  }
2386  return retn;
2387 }
2388 
2389 bool VCS_SOLVE::vcs_globStepDamp()
2390 {
2391  double* dptr = VCS_DATA_PTR(m_deltaGRxn_new);
2392 
2393  /* *************************************************** */
2394  /* **** CALCULATE SLOPE AT END OF THE STEP ********** */
2395  /* *************************************************** */
2396  double s2 = 0.0;
2397  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2398  size_t kspec = irxn + m_numComponents;
2399  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2400  s2 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
2401  }
2402  }
2403 
2404 
2405  /* *************************************************** */
2406  /* **** CALCULATE ORIGINAL SLOPE ********************* */
2407  /* ************************************************** */
2408  double s1 = 0.0;
2409  dptr = VCS_DATA_PTR(m_deltaGRxn_old);
2410  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2411  size_t kspec = irxn + m_numComponents;
2412  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2413  s1 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
2414  }
2415  }
2416 
2417  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2418  plogf(" --- subroutine FORCE: Beginning Slope = %g\n", s1);
2419  plogf(" --- subroutine FORCE: End Slope = %g\n", s2);
2420  }
2421 
2422  if (s1 > 0.0) {
2423  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2424  plogf(" --- subroutine FORCE produced no adjustments,");
2425  if (s1 < 1.0E-40) {
2426  plogf(" s1 positive but really small");
2427  } else {
2428  plogf(" failed s1 test");
2429  }
2430  plogendl();
2431  }
2432  return false;
2433  }
2434 
2435  if (s2 <= 0.0) {
2436  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2437  plogf(" --- subroutine FORCE produced no adjustments, s2 < 0");
2438  plogendl();
2439  }
2440  return false;
2441  }
2442 
2443  /* *************************************************** */
2444  /* **** FIT PARABOLA ********************************* */
2445  /* *************************************************** */
2446  double al = 1.0;
2447  if (fabs(s1 -s2) > 1.0E-200) {
2448  al = s1 / (s1 - s2);
2449  }
2450  if (al >= 0.95 || al < 0.0) {
2451  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2452  plogf(" --- subroutine FORCE produced no adjustments (al = %g)\n", al);
2453  }
2454  return false;
2455  }
2456  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2457  plogf(" --- subroutine FORCE produced a damping factor = %g\n", al);
2458  }
2459 
2460  /* *************************************************** */
2461  /* **** ADJUST MOLE NUMBERS, CHEM. POT *************** */
2462  /* *************************************************** */
2463  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2464  m_deltaGRxn_tmp = m_deltaGRxn_new;
2465  }
2466 
2467  dptr = VCS_DATA_PTR(m_molNumSpecies_new);
2468  for (size_t kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
2469  m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] +
2470  al * m_deltaMolNumSpecies[kspec];
2471  }
2472  for (size_t iph = 0; iph < m_numPhases; iph++) {
2473  m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + al * m_deltaPhaseMoles[iph];
2474  }
2475  vcs_updateVP(VCS_STATECALC_NEW);
2476 
2477  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2478  plogf(" --- subroutine FORCE adjusted the mole "
2479  "numbers, AL = %10.3f\n", al);
2480  }
2481  /*
2482  * Because we changed the mole numbers, we need to
2483  * calculate the chemical potentials again. If a major-
2484  * only step is being carried out, then we don't need to
2485  * update the minor noncomponents.
2486  */
2487  vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
2488  vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesRdc);
2489 
2490  /*
2491  * Evaluate DeltaG for all components if ITI=0, and for
2492  * major components only if ITI NE 0
2493  */
2494  vcs_deltag(0, false, VCS_STATECALC_NEW);
2495 
2496  dptr = VCS_DATA_PTR(m_deltaGRxn_new);
2497  s2 = 0.0;
2498  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
2499  size_t kspec = irxn + m_numComponents;
2500  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2501  s2 += dptr[irxn] * m_deltaMolNumSpecies[kspec];
2502  }
2503  }
2504 
2505  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2506  plogf(" --- subroutine FORCE: Adj End Slope = %g", s2);
2507  plogendl();
2508  }
2509  return true;
2510 }
2511 
2512 int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[],
2513  double ss[], double test, bool* const usedZeroedSpecies)
2514 {
2515  size_t k;
2516  bool lindep;
2517  size_t juse = npos;
2518  size_t jlose = npos;
2519  double* scrxn_ptr;
2520  clockWC tickTock;
2521  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2522  plogf(" ");
2523  for (size_t i=0; i<77; i++) {
2524  plogf("-");
2525  }
2526  plogf("\n");
2527  plogf(" --- Subroutine BASOPT called to ");
2528  if (doJustComponents) {
2529  plogf("calculate the number of components\n");
2530  } else {
2531  plogf("reevaluate the components\n");
2532  }
2533  if (m_debug_print_lvl >= 2) {
2534  plogf("\n");
2535  plogf(" --- Formula Matrix used in BASOPT calculation\n");
2536  plogf(" --- Active | ");
2537  for (size_t j = 0; j < m_numElemConstraints; j++) {
2538  plogf(" %1d ", m_elementActive[j]);
2539  }
2540  plogf("\n");
2541  plogf(" --- Species | ");
2542  for (size_t j = 0; j < m_numElemConstraints; j++) {
2543  plogf(" ");
2544  vcs_print_stringTrunc(m_elementName[j].c_str(), 8, 1);
2545  }
2546  plogf("\n");
2547  for (k = 0; k < m_numSpeciesTot; k++) {
2548  plogf(" --- ");
2549  vcs_print_stringTrunc(m_speciesName[k].c_str(), 11, 1);
2550  plogf(" | ");
2551  for (size_t j = 0; j < m_numElemConstraints; j++) {
2552  plogf(" %8.2g", m_formulaMatrix(k,j));
2553  }
2554  plogf("\n");
2555  }
2556  plogendl();
2557  }
2558  }
2559 
2560  /*
2561  * Calculate the maximum value of the number of components possible
2562  * It's equal to the minimum of the number of elements and the
2563  * number of total species.
2564  */
2565  size_t ncTrial = std::min(m_numElemConstraints, m_numSpeciesTot);
2566  m_numComponents = ncTrial;
2567  *usedZeroedSpecies = false;
2568  vector_int ipiv(ncTrial);
2569  int info;
2570 
2571  /*
2572  * Use a temporary work array for the mole numbers, aw[]
2573  */
2574  std::copy(m_molNumSpecies_old.begin(),
2575  m_molNumSpecies_old.begin() + m_numSpeciesTot, aw);
2576  /*
2577  * Take out the Voltage unknowns from consideration
2578  */
2579  for (k = 0; k < m_numSpeciesTot; k++) {
2580  if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2581  aw[k] = test;
2582  }
2583  }
2584 
2585  size_t jr = npos;
2586  /*
2587  * Top of a loop of some sort based on the index JR. JR is the
2588  * current number of component species found.
2589  */
2590  do {
2591  ++jr;
2592  /* - Top of another loop point based on finding a linearly */
2593  /* - independent species */
2594  do {
2595  /*
2596  * Search the remaining part of the mole fraction vector, AW,
2597  * for the largest remaining species. Return its identity in K.
2598  * The first search criteria is always the largest positive
2599  * magnitude of the mole number.
2600  */
2601  k = vcs_basisOptMax(aw, jr, m_numSpeciesTot);
2602 
2603  /*
2604  * The fun really starts when you have run out of species that have a significant
2605  * concentration. It becomes extremely important to make a good choice of which
2606  * species you want to pick to fill out the basis. Basically, you don't want to
2607  * use species with elements abundances which aren't pegged to zero. This means
2608  * that those modes will never be allowed to grow. You want to have the
2609  * best chance that the component will grow positively.
2610  *
2611  * Suppose you start with CH4, N2, as the only species with nonzero compositions.
2612  * You have the following abundances:
2613  *
2614  * Abundances:
2615  * ----------------
2616  * C 2.0
2617  * N 2.0
2618  * H 4.0
2619  * O 0.0
2620  *
2621  * For example, Make the following choice:
2622  *
2623  * CH4 N2 O choose -> OH
2624  * or
2625  * CH4 N2 O choose -> H2
2626  *
2627  * OH and H2 both fill out the basis. They will pass the algorithm. However,
2628  * choosing OH as the next species will create a situation where H2 can not
2629  * grow in concentration. This happened in practice, btw. The reason is that
2630  * the formation reaction for H2 will cause one of the component species
2631  * to go negative.
2632  *
2633  * The basic idea here is to pick a simple species whose mole number
2634  * can grow according to the element compositions. Candidates are still
2635  * filtered according to their linear independence.
2636  *
2637  * Note, if there is electronic charge and the electron species,
2638  * you should probably pick the electron as a component, if it
2639  * linearly independent. The algorithm below will do this automagically.
2640  *
2641  */
2642  if ((aw[k] != test) && aw[k] < VCS_DELETE_MINORSPECIES_CUTOFF) {
2643  *usedZeroedSpecies = true;
2644 
2645  double maxConcPossKspec = 0.0;
2646  double maxConcPoss = 0.0;
2647  size_t kfound = npos;
2648  int minNonZeroes = 100000;
2649  int nonZeroesKspec = 0;
2650  for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) {
2651  if (aw[kspec] >= 0.0) {
2652  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2653  maxConcPossKspec = 1.0E10;
2654  nonZeroesKspec = 0;
2655  for (size_t j = 0; j < m_numElemConstraints; ++j) {
2656  if (m_elementActive[j]) {
2657  if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
2658  double nu = m_formulaMatrix(kspec,j);
2659  if (nu != 0.0) {
2660  nonZeroesKspec++;
2661  maxConcPossKspec = std::min(m_elemAbundancesGoal[j] / nu, maxConcPossKspec);
2662  }
2663  }
2664  }
2665  }
2666  if ((maxConcPossKspec >= maxConcPoss) || (maxConcPossKspec > 1.0E-5)) {
2667  if (nonZeroesKspec <= minNonZeroes) {
2668  if (kfound == npos || nonZeroesKspec < minNonZeroes) {
2669  kfound = kspec;
2670  } else {
2671  // ok we are sitting pretty equal here decide on the raw ss Gibbs energy
2672  if (m_SSfeSpecies[kspec] <= m_SSfeSpecies[kfound]) {
2673  kfound = kspec;
2674  }
2675  }
2676  }
2677  minNonZeroes = std::min(minNonZeroes, nonZeroesKspec);
2678  maxConcPoss = std::max(maxConcPoss, maxConcPossKspec);
2679  }
2680  }
2681  }
2682  }
2683  if (kfound == npos) {
2684  double gmin = 0.0;
2685  kfound = k;
2686  for (size_t kspec = ncTrial; kspec < m_numSpeciesTot; kspec++) {
2687  if (aw[kspec] >= 0.0) {
2688  size_t irxn = kspec - ncTrial;
2689  if (m_deltaGRxn_new[irxn] < gmin) {
2690  gmin = m_deltaGRxn_new[irxn];
2691  kfound = kspec;
2692  }
2693  }
2694  }
2695  }
2696  k = kfound;
2697  }
2698 
2699 
2700  if (aw[k] == test) {
2701  m_numComponents = jr;
2702  ncTrial = m_numComponents;
2703  size_t numPreDeleted = m_numRxnTot - m_numRxnRdc;
2704  if (numPreDeleted != (m_numSpeciesTot - m_numSpeciesRdc)) {
2705  throw CanteraError("VCS_SOLVE::vcs_basopt", "we shouldn't be here");
2706  }
2707  m_numRxnTot = m_numSpeciesTot - ncTrial;
2708  m_numRxnRdc = m_numRxnTot - numPreDeleted;
2709  m_numSpeciesRdc = m_numSpeciesTot - numPreDeleted;
2710  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
2711  m_indexRxnToSpecies[i] = ncTrial + i;
2712  }
2713  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2714  plogf(" --- Total number of components found = %3d (ne = %d)\n ",
2715  ncTrial, m_numElemConstraints);
2716  }
2717  goto L_END_LOOP;
2718  }
2719  /*
2720  * Assign a small negative number to the component that we have
2721  * just found, in order to take it out of further consideration.
2722  */
2723  aw[k] = test;
2724  /* *********************************************************** */
2725  /* **** CHECK LINEAR INDEPENDENCE WITH PREVIOUS SPECIES ****** */
2726  /* *********************************************************** */
2727  /*
2728  * Modified Gram-Schmidt Method, p. 202 Dalquist
2729  * QR factorization of a matrix without row pivoting.
2730  */
2731  size_t jl = jr;
2732  for (size_t j = 0; j < m_numElemConstraints; ++j) {
2733  sm[j + jr*m_numElemConstraints] = m_formulaMatrix(k,j);
2734  }
2735  if (jl > 0) {
2736  /*
2737  * Compute the coefficients of JA column of the
2738  * the upper triangular R matrix, SS(J) = R_J_JR
2739  * (this is slightly different than Dalquist)
2740  * R_JA_JA = 1
2741  */
2742  for (size_t j = 0; j < jl; ++j) {
2743  ss[j] = 0.0;
2744  for (size_t i = 0; i < m_numElemConstraints; ++i) {
2745  ss[j] += sm[i + jr*m_numElemConstraints] * sm[i + j*m_numElemConstraints];
2746  }
2747  ss[j] /= sa[j];
2748  }
2749  /*
2750  * Now make the new column, (*,JR), orthogonal to the
2751  * previous columns
2752  */
2753  for (size_t j = 0; j < jl; ++j) {
2754  for (size_t l = 0; l < m_numElemConstraints; ++l) {
2755  sm[l + jr*m_numElemConstraints] -= ss[j] * sm[l + j*m_numElemConstraints];
2756  }
2757  }
2758  }
2759  /*
2760  * Find the new length of the new column in Q.
2761  * It will be used in the denominator in future row calcs.
2762  */
2763  sa[jr] = 0.0;
2764  for (size_t ml = 0; ml < m_numElemConstraints; ++ml) {
2765  sa[jr] += pow(sm[ml + jr*m_numElemConstraints], 2);
2766  }
2767  /* **************************************************** */
2768  /* **** IF NORM OF NEW ROW .LT. 1E-3 REJECT ********** */
2769  /* **************************************************** */
2770  lindep = (sa[jr] < 1.0e-6);
2771  } while (lindep);
2772  /* ****************************************** */
2773  /* **** REARRANGE THE DATA ****************** */
2774  /* ****************************************** */
2775  if (jr != k) {
2776  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2777  plogf(" --- %-12.12s", (m_speciesName[k]).c_str());
2778  if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2779  plogf("(Volts = %9.2g)", m_molNumSpecies_old[k]);
2780  } else {
2781  plogf("(%9.2g)", m_molNumSpecies_old[k]);
2782  }
2783  plogf(" replaces %-12.12s", m_speciesName[jr].c_str());
2784  if (m_speciesUnknownType[jr] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2785  plogf("(Volts = %9.2g)", m_molNumSpecies_old[jr]);
2786  } else {
2787  plogf("(%9.2g)", m_molNumSpecies_old[jr]);
2788  }
2789  plogf(" as component %3d\n", jr);
2790  }
2791  vcs_switch_pos(false, jr, k);
2792  std::swap(aw[jr], aw[k]);
2793  }
2794  else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2795  plogf(" --- %-12.12s", m_speciesName[k].c_str());
2796  if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2797  plogf("(Volts = %9.2g) remains ", m_molNumSpecies_old[k]);
2798  } else {
2799  plogf("(%9.2g) remains ", m_molNumSpecies_old[k]);
2800  }
2801  plogf(" as component %3d\n", jr);
2802  }
2803  /* - entry point from up above */
2804 L_END_LOOP:
2805  ;
2806  /*
2807  * If we haven't found enough components, go back
2808  * and find some more. (nc -1 is used below, because
2809  * jr is counted from 0, via the C convention.
2810  */
2811  } while (jr < (ncTrial-1));
2812 
2813  if (doJustComponents) {
2814  goto L_CLEANUP;
2815  }
2816  /* ****************************************************** */
2817  /* **** EVALUATE THE STOICHIOMETRY ********************** */
2818  /* ****************************************************** */
2819  /*
2820  * Formulate the matrix problem for the stoichiometric
2821  * coefficients. CX + B = 0
2822  * C will be an nc x nc matrix made up of the formula
2823  * vectors for the components.
2824  * n RHS's will be solved for. Thus, B is an nc x n
2825  * matrix.
2826  *
2827  * BIG PROBLEM 1/21/99:
2828  *
2829  * This algorithm makes the assumption that the
2830  * first nc rows of the formula matrix aren't rank deficient.
2831  * However, this might not be the case. For example, assume
2832  * that the first element in m_formulaMatrix[] is argon. Assume that
2833  * no species in the matrix problem actually includes argon.
2834  * Then, the first row in sm[], below will be identically
2835  * zero. bleh.
2836  * What needs to be done is to perform a rearrangement
2837  * of the ELEMENTS -> i.e. rearrange, m_formulaMatrix, sp,
2838  * and m_elemAbundancesGoal, such
2839  * that the first nc elements form in combination with the
2840  * nc components create an invertible sm[]. not a small
2841  * project, but very doable.
2842  * An alternative would be to turn the matrix problem
2843  * below into an ne x nc problem, and do QR elimination instead
2844  * of Gauss-Jordan elimination.
2845  * Note the rearrangement of elements need only be done once
2846  * in the problem. It's actually very similar to the top of
2847  * this program with ne being the species and nc being the
2848  * elements!!
2849  */
2850  for (size_t j = 0; j < ncTrial; ++j) {
2851  for (size_t i = 0; i < ncTrial; ++i) {
2852  sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,i);
2853  }
2854  }
2855  for (size_t i = 0; i < m_numRxnTot; ++i) {
2856  k = m_indexRxnToSpecies[i];
2857  for (size_t j = 0; j < ncTrial; ++j) {
2858  m_stoichCoeffRxnMatrix(j,i) = - m_formulaMatrix(k,j);
2859  }
2860  }
2861  // Solve the linear system to calculate the reaction matrix,
2862  // m_stoichCoeffRxnMatrix.
2863  ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info);
2864  if (info) {
2865  plogf("vcs_solve_TP ERROR: Error factorizing stoichiometric coefficient matrix\n");
2866  return VCS_FAILED_CONVERGENCE;
2867  }
2868  ct_dgetrs(ctlapack::NoTranspose, ncTrial, m_numRxnTot, sm, m_numElemConstraints,
2869  &ipiv[0], m_stoichCoeffRxnMatrix.ptrColumn(0), m_numElemConstraints, info);
2870 
2871  /*
2872  * NOW, if we have interfacial voltage unknowns, what we did
2873  * was just wrong -> hopefully it didn't blow up. Redo the problem.
2874  * Search for inactive E
2875  */
2876  juse = npos;
2877  jlose = npos;
2878  for (size_t j = 0; j < m_numElemConstraints; j++) {
2879  if (!(m_elementActive[j])) {
2880  if (!strcmp((m_elementName[j]).c_str(), "E")) {
2881  juse = j;
2882  }
2883  }
2884  }
2885  for (size_t j = 0; j < m_numElemConstraints; j++) {
2886  if (m_elementActive[j]) {
2887  if (!strncmp((m_elementName[j]).c_str(), "cn_", 3)) {
2888  jlose = j;
2889  }
2890  }
2891  }
2892  for (k = 0; k < m_numSpeciesTot; k++) {
2893  if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2894 
2895  for (size_t j = 0; j < ncTrial; ++j) {
2896  for (size_t i = 0; i < ncTrial; ++i) {
2897  if (i == jlose) {
2898  sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,juse);
2899  } else {
2900  sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,i);
2901  }
2902  }
2903  }
2904  for (size_t i = 0; i < m_numRxnTot; ++i) {
2905  k = m_indexRxnToSpecies[i];
2906  for (size_t j = 0; j < ncTrial; ++j) {
2907  if (j == jlose) {
2908  aw[j] = - m_formulaMatrix(k,juse);
2909  } else {
2910  aw[j] = - m_formulaMatrix(k,j);
2911  }
2912  }
2913  }
2914 
2915  ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info);
2916  if (info) {
2917  plogf("vcs_solve_TP ERROR: Error factorizing matrix\n");
2918  return VCS_FAILED_CONVERGENCE;
2919  }
2920  ct_dgetrs(ctlapack::NoTranspose, ncTrial, 1, sm, m_numElemConstraints,
2921  &ipiv[0], aw, m_numElemConstraints, info);
2922  size_t i = k - ncTrial;
2923  for (size_t j = 0; j < ncTrial; j++) {
2924  m_stoichCoeffRxnMatrix(j,i) = aw[j];
2925  }
2926  }
2927  }
2928 
2929 
2930  /*
2931  * Calculate the szTmp array for each formation reaction
2932  */
2933  for (size_t i = 0; i < m_numRxnTot; i++) {
2934  double szTmp = 0.0;
2935  for (size_t j = 0; j < ncTrial; j++) {
2936  szTmp += fabs(m_stoichCoeffRxnMatrix(j,i));
2937  }
2938  m_scSize[i] = szTmp;
2939  }
2940 
2941  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
2942  plogf(" --- Components:");
2943  for (size_t j = 0; j < ncTrial; j++) {
2944  plogf(" %3d", j);
2945  }
2946  plogf("\n --- Components Moles:");
2947  for (size_t j = 0; j < ncTrial; j++) {
2948  if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2949  plogf(" % -10.3E", 0.0);
2950  } else {
2951  plogf(" % -10.3E", m_molNumSpecies_old[j]);
2952  }
2953  }
2954  plogf("\n --- NonComponent| Moles |");
2955  for (size_t j = 0; j < ncTrial; j++) {
2956  plogf(" %10.10s", m_speciesName[j].c_str());
2957  }
2958  plogf("\n");
2959  for (size_t i = 0; i < m_numRxnTot; i++) {
2960  plogf(" --- %3d ", m_indexRxnToSpecies[i]);
2961  plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str());
2962  if (m_speciesUnknownType[m_indexRxnToSpecies[i]] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2963  plogf("|% -10.3E|", 0.0);
2964  } else {
2965  plogf("|% -10.3E|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
2966  }
2967  for (size_t j = 0; j < ncTrial; j++) {
2968  plogf(" %+7.3f", m_stoichCoeffRxnMatrix(j,i));
2969  }
2970  plogf("\n");
2971  }
2972 
2973 
2974  /*
2975  * Manual check on the satisfaction of the reaction matrix's ability
2976  * to conserve elements
2977  */
2978  double sumMax = -1.0;
2979  size_t iMax = npos;
2980  size_t jMax = npos;
2981  for (size_t i = 0; i < m_numRxnTot; ++i) {
2982  k = m_indexRxnToSpecies[i];
2983  double sum;
2984  for (size_t j = 0; j < ncTrial; ++j) {
2985  if (j == jlose) {
2986  sum = m_formulaMatrix(k,juse);
2987  for (size_t n = 0; n < ncTrial; n++) {
2988  double numElements = m_formulaMatrix(n,juse);
2989  double coeff = m_stoichCoeffRxnMatrix(n,i);
2990  sum += coeff * numElements;
2991  }
2992  } else {
2993  sum = m_formulaMatrix(k,j);
2994  for (size_t n = 0; n < ncTrial; n++) {
2995  double numElements = m_formulaMatrix(n,j);
2996  double coeff = m_stoichCoeffRxnMatrix(n,i);
2997  sum += coeff * numElements;
2998  }
2999  }
3000  if (fabs(sum) > sumMax) {
3001  sumMax = fabs(sum);
3002  iMax = i;
3003  jMax = j;
3004  if (j == jlose) {
3005  jMax = juse;
3006  }
3007  }
3008  if (fabs(sum) > 1.0E-6) {
3009  throw CanteraError("VCS_SOLVE::vcs_basopt", "we have a prob");
3010  }
3011  }
3012  }
3013  plogf(" --- largest error in Stoich coeff = %g at rxn = %d ", sumMax, iMax);
3014  plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[iMax]].c_str());
3015  plogf(" element = %d ", jMax);
3016  plogf("%-5.5s", m_elementName[jMax].c_str());
3017  plogf("\n");
3018  plogf(" ");
3019  for (size_t i=0; i<77; i++) {
3020  plogf("-");
3021  }
3022  plogf("\n");
3023  }
3024  /* **************************************************** */
3025  /* **** EVALUATE DELTA N VALUES *********************** */
3026  /* **************************************************** */
3027  /*
3028  * Evaluate the change in gas and liquid total moles
3029  * due to reaction vectors, DNG and DNL.
3030  */
3031 
3032  /*
3033  * Zero out the change of Phase Moles array
3034  */
3035  m_deltaMolNumPhase.zero();
3036  m_phaseParticipation.zero();
3037  /*
3038  * Loop over each reaction, creating the change in Phase Moles
3039  * array, m_deltaMolNumPhase(iphase,irxn),
3040  * and the phase participation array, PhaseParticipation[irxn][iphase]
3041  */
3042  for (size_t irxn = 0; irxn < m_numRxnTot; ++irxn) {
3043  scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3044  size_t kspec = m_indexRxnToSpecies[irxn];
3045  size_t iph = m_phaseID[kspec];
3046  m_deltaMolNumPhase(iph,irxn) = 1.0;
3047  m_phaseParticipation(iph,irxn)++;
3048  for (size_t j = 0; j < ncTrial; ++j) {
3049  iph = m_phaseID[j];
3050  if (fabs(scrxn_ptr[j]) <= 1.0e-6) {
3051  scrxn_ptr[j] = 0.0;
3052  } else {
3053  m_deltaMolNumPhase(iph,irxn) += scrxn_ptr[j];
3054  m_phaseParticipation(iph,irxn)++;
3055  }
3056  }
3057  }
3058 
3059 L_CLEANUP:
3060  ;
3061  double tsecond = tickTock.secondsWC();
3062  m_VCount->Time_basopt += tsecond;
3063  (m_VCount->Basis_Opts)++;
3064  return VCS_SUCCESS;
3065 }
3066 
3067 size_t VCS_SOLVE::vcs_basisOptMax(const double* const molNum, const size_t j,
3068  const size_t n)
3069 {
3070  /*
3071  * The factors of 1.01 and 1.001 are placed in this routine for a purpose.
3072  * The purpose is to ensure that roundoff errors don't influence major
3073  * decisions. This means that the optimized and non-optimized versions of
3074  * the code remain close to each other.
3075  *
3076  * (we try to avoid the logic: a = b
3077  * if (a > b) { do this }
3078  * else { do something else }
3079  * because roundoff error makes a difference in the inequality evaluation)
3080  *
3081  * Mole numbers are frequently equal to each other in equilibrium problems
3082  * due to constraints. Swaps are only done if there are a 1% difference in
3083  * the mole numbers. Of course this logic isn't foolproof.
3084  */
3085  size_t largest = j;
3086  double big = molNum[j] * m_spSize[j] * 1.01;
3087  if (m_spSize[j] <= 0.0) {
3088  throw CanteraError("VCS_SOLVE::vcs_basisOptMax",
3089  "spSize is nonpositive");
3090  }
3091  for (size_t i = j + 1; i < n; ++i) {
3092  if (m_spSize[i] <= 0.0) {
3093  throw CanteraError("VCS_SOLVE::vcs_basisOptMax",
3094  "spSize is nonpositive");
3095  }
3096  bool doSwap = false;
3097  if (m_SSPhase[j]) {
3098  doSwap = (molNum[i] * m_spSize[i]) > (big);
3099  if (!m_SSPhase[i]) {
3100  if (doSwap) {
3101  doSwap = (molNum[i]) > (molNum[largest] * 1.001);
3102  }
3103  }
3104  } else {
3105  if (m_SSPhase[i]) {
3106  doSwap = (molNum[i] * m_spSize[i]) > (big);
3107  if (!doSwap) {
3108  doSwap = (molNum[i]) > (molNum[largest] * 1.001);
3109  }
3110  } else {
3111  doSwap = (molNum[i] * m_spSize[i]) > (big);
3112  }
3113  }
3114  if (doSwap) {
3115  largest = i;
3116  big = molNum[i] * m_spSize[i] * 1.01;
3117  }
3118  }
3119  return largest;
3120 }
3121 
3122 int VCS_SOLVE::vcs_species_type(const size_t kspec) const
3123 {
3124 
3125  // ---------- Treat special cases first ---------------------
3126 
3127 
3128  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3129  return VCS_SPECIES_INTERFACIALVOLTAGE;
3130  }
3131 
3132  size_t iph = m_phaseID[kspec];
3133  int irxn = int(kspec) - int(m_numComponents);
3134  vcs_VolPhase* VPhase = m_VolPhaseList[iph];
3135  int phaseExist = VPhase->exists();
3136 
3137  // ---------- Treat zeroed out species first ----------------
3138 
3139  if (m_molNumSpecies_old[kspec] <= 0.0) {
3140 
3141  if (m_tPhaseMoles_old[iph] <= 0.0) {
3142  if (!m_SSPhase[kspec]) {
3143  return VCS_SPECIES_ZEROEDMS;
3144  }
3145  }
3146 
3147  /*
3148  * see if the species has an element
3149  * which is so low that species will always be zero
3150  *
3151  */
3152  for (size_t j = 0; j < m_numElemConstraints; ++j) {
3153  int elType = m_elType[j];
3154  if (elType == VCS_ELEM_TYPE_ABSPOS) {
3155  double atomComp = m_formulaMatrix(kspec,j);
3156  if (atomComp > 0.0) {
3157  double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
3158  if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
3159  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3160  plogf(" --- %s can not be nonzero because"
3161  " needed element %s is zero\n",
3162  m_speciesName[kspec].c_str(), (m_elementName[j]).c_str());
3163  }
3164  if (m_SSPhase[kspec]) {
3165  return VCS_SPECIES_ZEROEDSS;
3166  } else {
3167  return VCS_SPECIES_STOICHZERO;
3168  }
3169  }
3170  }
3171  }
3172  }
3173 
3174  /*
3175  * The Gibbs free energy for this species is such that
3176  * it will pop back into existence.
3177  */
3178  /*
3179  * -> An exception to this is if a needed regular element
3180  * is also zeroed out. Then, don't pop the phase or the species back into
3181  * existence.
3182  */
3183  if (irxn >= 0) {
3184  for (size_t j = 0; j < m_numComponents; ++j) {
3185  double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
3186  if (stoicC != 0.0) {
3187  double negChangeComp = - stoicC;
3188  if (negChangeComp > 0.0) {
3189  if (m_molNumSpecies_old[j] < 1.0E-60) {
3190  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3191  plogf(" --- %s is prevented from popping into existence because"
3192  " a needed component to be consumed, %s, has a zero mole number\n",
3193  m_speciesName[kspec].c_str(), m_speciesName[j].c_str());
3194  }
3195  if (m_SSPhase[kspec]) {
3196  return VCS_SPECIES_ZEROEDSS;
3197  } else {
3198  return VCS_SPECIES_STOICHZERO;
3199  }
3200  }
3201  } else if (negChangeComp < 0.0) {
3202  size_t jph = m_phaseID[j];
3203  vcs_VolPhase* jVPhase = m_VolPhaseList[jph];
3204  if (jVPhase->exists() <= 0) {
3205  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3206  plogf(" --- %s is prevented from popping into existence because"
3207  " a needed component %s is in a zeroed-phase that would be "
3208  "popped into existence at the same time\n",
3209  m_speciesName[kspec].c_str(), m_speciesName[j].c_str());
3210  }
3211  if (m_SSPhase[kspec]) {
3212  return VCS_SPECIES_ZEROEDSS;
3213  } else {
3214  return VCS_SPECIES_STOICHZERO;
3215  }
3216  }
3217  }
3218  }
3219  }
3220  }
3221 
3222  if (irxn >= 0) {
3223  if (m_deltaGRxn_old[irxn] >= 0.0) {
3224  /*
3225  * We are here when the species is or should remain zeroed out
3226  */
3227  if (m_SSPhase[kspec]) {
3228  return VCS_SPECIES_ZEROEDSS;
3229  } else {
3230  if (phaseExist >= VCS_PHASE_EXIST_YES) {
3231  return VCS_SPECIES_ACTIVEBUTZERO;
3232  } else if (phaseExist == VCS_PHASE_EXIST_ZEROEDPHASE) {
3233  return VCS_SPECIES_ZEROEDPHASE;
3234  } else {
3235  return VCS_SPECIES_ZEROEDMS;
3236  }
3237  }
3238  }
3239  }
3240  /*
3241  * If the current phase already exists,
3242  */
3243  if (m_tPhaseMoles_old[iph] > 0.0) {
3244  if (m_SSPhase[kspec]) {
3245  return VCS_SPECIES_MAJOR;
3246  } else {
3247  return VCS_SPECIES_ACTIVEBUTZERO;
3248  }
3249  }
3250 
3251  /*
3252  * The Gibbs free energy for this species is such that
3253  * it will pop back into existence.
3254  *
3255  * -> Set it to a major species in anticipation.
3256  * -> note, if we had an estimate for the emerging mole
3257  * fraction of the species in the phase, we could do
3258  * better here.
3259  */
3260  if (m_tPhaseMoles_old[iph] <= 0.0) {
3261  if (m_SSPhase[kspec]) {
3262  return VCS_SPECIES_MAJOR;
3263  } else {
3264  return VCS_SPECIES_ZEROEDMS;
3265  }
3266  }
3267  }
3268 
3269  // ---------- Treat species with non-zero mole numbers next ------------
3270 
3271  /*
3272  * Always treat species in single species phases as majors if the
3273  * phase exists.
3274  */
3275  if (m_SSPhase[kspec]) {
3276  return VCS_SPECIES_MAJOR;
3277  }
3278 
3279  /*
3280  * Check to see whether the current species is a major component
3281  * of its phase. If it is, it is a major component. This is consistent
3282  * with the above rule about single species phases. A major component
3283  * (i.e., a species with a high mole fraction)
3284  * in any phase is always treated as a major species
3285  */
3286  if (m_molNumSpecies_old[kspec] > (m_tPhaseMoles_old[iph] * 0.001)) {
3287  return VCS_SPECIES_MAJOR;
3288  }
3289 
3290  /*
3291  * Main check in the loop:
3292  * Check to see if there is a component with a mole number that is
3293  * within a factor of 100 of the current species.
3294  * If there is and that component is not part of a single species
3295  * phase and shares a non-zero stoichiometric coefficient, then
3296  * the current species is a major species.
3297  */
3298  if (irxn < 0) {
3299  return VCS_SPECIES_MAJOR;
3300  } else {
3301  double szAdj = m_scSize[irxn] * std::sqrt((double)m_numRxnTot);
3302  for (size_t k = 0; k < m_numComponents; ++k) {
3303  if (!(m_SSPhase[k])) {
3304  if (m_stoichCoeffRxnMatrix(k,irxn) != 0.0) {
3305  if (m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) {
3306  return VCS_SPECIES_MAJOR;
3307  }
3308  }
3309  }
3310  }
3311  }
3312  return VCS_SPECIES_MINOR;
3313 }
3314 
3315 void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc,
3316  const size_t iph, const double* const molNum,
3317  double* const ac, double* const mu_i,
3318  const bool do_deleted)
3319 {
3320  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
3321  size_t nkk = Vphase->nSpecies();
3322  double tMoles = TPhInertMoles[iph];
3323  for (size_t k = 0; k < nkk; k++) {
3324  size_t kspec = Vphase->spGlobalIndexVCS(k);
3325  tMoles += molNum[kspec];
3326  }
3327  double tlogMoles = 0.0;
3328  if (tMoles > 0.0) {
3329  tlogMoles = log(tMoles);
3330  }
3331 
3332  Vphase->setMolesFromVCS(stateCalc, molNum);
3333  Vphase->sendToVCS_ActCoeff(stateCalc, ac);
3334 
3335  double phi = Vphase->electricPotential();
3336  double Faraday_phi = m_Faraday_dim * phi;
3337 
3338  for (size_t k = 0; k < nkk; k++) {
3339  size_t kspec = Vphase->spGlobalIndexVCS(k);
3340  if (kspec >= m_numComponents) {
3341  if (!do_deleted &&
3342  (m_speciesStatus[kspec] == VCS_SPECIES_DELETED)) {
3343  continue;
3344  }
3345  }
3346  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3347  AssertThrowMsg(molNum[kspec] == phi, "VCS_SOLVE::vcs_chemPotPhase",
3348  "We have an inconsistency!");
3349  AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_chemPotPhase",
3350  "We have an unexpected situation!");
3351  mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
3352  } else {
3353  if (m_SSPhase[kspec]) {
3354  mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
3355  } else if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
3356  mu_i[kspec] = m_SSfeSpecies[kspec] + log(ac[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
3357  - tlogMoles - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
3358  } else {
3359  mu_i[kspec] = m_SSfeSpecies[kspec] + log(ac[kspec] * molNum[kspec])
3360  - tlogMoles - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
3361  }
3362  }
3363  }
3364 }
3365 
3366 void VCS_SOLVE::vcs_dfe(const int stateCalc,
3367  const int ll, const size_t lbot, const size_t ltop)
3368 {
3369  double* tPhMoles_ptr=0;
3370  double* actCoeff_ptr=0;
3371  double* feSpecies=0;
3372  double* molNum=0;
3373  if (stateCalc == VCS_STATECALC_OLD) {
3374  feSpecies = VCS_DATA_PTR(m_feSpecies_old);
3375  tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
3376  actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
3377  molNum = VCS_DATA_PTR(m_molNumSpecies_old);
3378  } else if (stateCalc == VCS_STATECALC_NEW) {
3379  feSpecies = VCS_DATA_PTR(m_feSpecies_new);
3380  tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
3381  actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
3382  molNum = VCS_DATA_PTR(m_molNumSpecies_new);
3383  } else if (DEBUG_MODE_ENABLED) {
3384  throw CanteraError("VCS_SOLVE::vcs_dfe",
3385  "Subroutine vcs_dfe called with bad stateCalc value: "+
3386  int2str(stateCalc));
3387  }
3388 
3389  AssertThrowMsg(m_unitsState != VCS_DIMENSIONAL_G, "VCS_SOLVE::vcs_dfe",
3390  "called with wrong units state");
3391 
3392  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3393  if (ll == 0) {
3394  if (lbot != 0) {
3395  plogf(" --- Subroutine vcs_dfe called for one species: ");
3396  plogf("%-12.12s", m_speciesName[lbot].c_str());
3397  } else {
3398  plogf(" --- Subroutine vcs_dfe called for all species");
3399  }
3400  } else if (ll > 0) {
3401  plogf(" --- Subroutine vcs_dfe called for components and minors");
3402  } else {
3403  plogf(" --- Subroutine vcs_dfe called for components and majors");
3404  }
3405  if (stateCalc == VCS_STATECALC_NEW) {
3406  plogf(" using tentative solution");
3407  }
3408  plogendl();
3409  }
3410 
3411  double* tlogMoles = VCS_DATA_PTR(m_TmpPhase);
3412  /*
3413  * Might as well recalculate the phase mole vector
3414  * and compare to the stored one. They should be correct.
3415  */
3416  double* tPhInertMoles = VCS_DATA_PTR(TPhInertMoles);
3417  for (size_t iph = 0; iph < m_numPhases; iph++) {
3418  tlogMoles[iph] = tPhInertMoles[iph];
3419 
3420  }
3421  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
3422  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3423  size_t iph = m_phaseID[kspec];
3424  tlogMoles[iph] += molNum[kspec];
3425  }
3426  }
3427  for (size_t iph = 0; iph < m_numPhases; iph++) {
3428  AssertThrowMsg(vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph]),
3429  "VCS_SOLVE::vcs_dfe",
3430  "phase Moles may be off, iph = " + int2str(iph) + ", " +
3431  fp2str(tlogMoles[iph]) + " " + fp2str(tPhMoles_ptr[iph]));
3432  }
3433  m_TmpPhase.assign(m_TmpPhase.size(), 0.0);
3434  for (size_t iph = 0; iph < m_numPhases; iph++) {
3435  if (tPhMoles_ptr[iph] > 0.0) {
3436  tlogMoles[iph] = log(tPhMoles_ptr[iph]);
3437  }
3438  }
3439 
3440  size_t l1, l2;
3441  if (ll != 0) {
3442  l1 = lbot;
3443  l2 = m_numComponents;
3444  } else {
3445  l1 = lbot;
3446  l2 = ltop;
3447  }
3448 
3449  /*
3450  * Calculate activity coefficients for all phases that are
3451  * not current. Here we also trigger an update check for each
3452  * VolPhase to see if its mole numbers are current with vcs
3453  */
3454  for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
3455  vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
3456  Vphase->updateFromVCS_MoleNumbers(stateCalc);
3457  if (!Vphase->m_singleSpecies) {
3458  Vphase->sendToVCS_ActCoeff(stateCalc, VCS_DATA_PTR(actCoeff_ptr));
3459  }
3460  m_phasePhi[iphase] = Vphase->electricPotential();
3461  }
3462  /* ************************************************************** */
3463  /* **** ALL SPECIES, OR COMPONENTS ****************************** */
3464  /* ************************************************************** */
3465  /*
3466  * Do all of the species when LL = 0. Then we are done for the routine
3467  * When LL ne 0., just do the initial components. We will then
3468  * finish up below with loops over either the major noncomponent
3469  * species or the minor noncomponent species.
3470  */
3471  for (size_t kspec = l1; kspec < l2; ++kspec) {
3472  size_t iphase = m_phaseID[kspec];
3473  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3474  AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase], "VCS_SOLVE::vcs_dfe",
3475  "We have an inconsistency!");
3476  AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_dfe",
3477  "We have an unexpected situation!");
3478  feSpecies[kspec] = m_SSfeSpecies[kspec]
3479  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3480  } else {
3481  if (m_SSPhase[kspec]) {
3482  feSpecies[kspec] = m_SSfeSpecies[kspec]
3483  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3484  } else if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
3485  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
3486  feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
3487  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3488  } else {
3489  if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
3490  size_t iph = m_phaseID[kspec];
3491  if (tPhMoles_ptr[iph] > 0.0) {
3492  feSpecies[kspec] = m_SSfeSpecies[kspec]
3493  + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
3494  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
3495  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3496  } else {
3497  feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
3498  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3499  }
3500  } else {
3501  feSpecies[kspec] = m_SSfeSpecies[kspec]
3502  + log(actCoeff_ptr[kspec] * molNum[kspec])
3503  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
3504  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3505  }
3506  }
3507  }
3508  }
3509  /* ************************************************ */
3510  /* **** MAJORS ONLY ******************************* */
3511  /* ************************************************ */
3512  if (ll < 0) {
3513  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3514  size_t kspec = m_indexRxnToSpecies[irxn];
3515  if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
3516  size_t iphase = m_phaseID[kspec];
3517  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3518  AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase], "VCS_SOLVE::vcs_dfe",
3519  "We have an inconsistency!");
3520  AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_dfe",
3521  "We have an unexpected situation!");
3522  feSpecies[kspec] = m_SSfeSpecies[kspec]
3523  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3524  } else {
3525  if (m_SSPhase[kspec]) {
3526  feSpecies[kspec] = m_SSfeSpecies[kspec]
3527  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3528  } else if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
3529  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
3530  feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
3531  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3532  } else {
3533  if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
3534  size_t iph = m_phaseID[kspec];
3535  if (tPhMoles_ptr[iph] > 0.0) {
3536  feSpecies[kspec] = m_SSfeSpecies[kspec]
3537  + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
3538  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
3539  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3540  } else {
3541  feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
3542  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3543  }
3544  } else {
3545  feSpecies[kspec] = m_SSfeSpecies[kspec]
3546  + log(actCoeff_ptr[kspec] * molNum[kspec])
3547  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
3548  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3549  }
3550  }
3551  }
3552  }
3553  }
3554  /* ************************************************ */
3555  /* **** MINORS ONLY ******************************* */
3556  /* ************************************************ */
3557  } else if (ll > 0) {
3558  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3559  size_t kspec = m_indexRxnToSpecies[irxn];
3560  if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
3561  size_t iphase = m_phaseID[kspec];
3562  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3563  AssertThrowMsg(molNum[kspec] == m_phasePhi[iphase], "VCS_SOLVE::vcs_dfe",
3564  "We have an inconsistency!");
3565  AssertThrowMsg(m_chargeSpecies[kspec] == -1.0, "VCS_SOLVE::vcs_dfe",
3566  "We have an unexpected situation!");
3567  feSpecies[kspec] = m_SSfeSpecies[kspec]
3568  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ;
3569  } else {
3570  if (m_SSPhase[kspec]) {
3571  feSpecies[kspec] = m_SSfeSpecies[kspec]
3572  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3573  } else if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
3574  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
3575  feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
3576  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3577  } else {
3578  if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
3579  size_t iph = m_phaseID[kspec];
3580  if (tPhMoles_ptr[iph] > 0.0) {
3581  feSpecies[kspec] = m_SSfeSpecies[kspec]
3582  + log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
3583  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
3584  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3585  } else {
3586  feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
3587  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3588  }
3589  } else {
3590  feSpecies[kspec] = m_SSfeSpecies[kspec]
3591  + log(actCoeff_ptr[kspec] * molNum[kspec])
3592  - tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
3593  + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
3594  }
3595  }
3596  }
3597  }
3598  }
3599  }
3600 
3601 
3602 }
3603 
3604 void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const
3605 {
3606  double mfValue = 1.0;
3607  bool zeroedPhase = false;
3608 
3609  const double* molNum = VCS_DATA_PTR(m_molNumSpecies_old);
3610  const double* actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
3611  if (stateCalc == VCS_STATECALC_NEW) {
3612  actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
3613  molNum = VCS_DATA_PTR(m_molNumSpecies_new);
3614  }
3615 
3616  double* tMoles = VCS_DATA_PTR(m_TmpPhase);
3617  const double* tPhInertMoles = VCS_DATA_PTR(TPhInertMoles);
3618  for (size_t iph = 0; iph < m_numPhases; iph++) {
3619  tMoles[iph] = tPhInertMoles[iph];
3620  }
3621  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
3622  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3623  size_t iph = m_phaseID[kspec];
3624  tMoles[iph] += molNum[kspec];
3625  }
3626  }
3627 
3628  double RT = m_temperature * GasConstant;
3629  printf(" --- CHEMICAL POT TABLE (J/kmol) Name PhID MolFR ChemoSS "
3630  " logMF Gamma Elect extra ElectrChem\n");
3631  printf(" ");
3632  writeline('-', 132);
3633 
3634  for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
3635  mfValue = 1.0;
3636  size_t iphase = m_phaseID[kspec];
3637  const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
3638  if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
3639  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
3640  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS)) {
3641  zeroedPhase = true;
3642  } else {
3643  zeroedPhase = false;
3644  }
3645  if (tMoles[iphase] > 0.0) {
3646  if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
3647  mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tMoles[iphase];
3648  } else {
3649  mfValue = molNum[kspec]/tMoles[iphase];
3650  }
3651  } else {
3652  size_t klocal = m_speciesLocalPhaseIndex[kspec];
3653  mfValue = Vphase->moleFraction(klocal);
3654  }
3655  double volts = Vphase->electricPotential();
3656  double elect = m_chargeSpecies[kspec] * m_Faraday_dim * volts;
3657  double comb = - m_lnMnaughtSpecies[kspec];
3658  double total = (m_SSfeSpecies[kspec] + log(mfValue) + elect + log(actCoeff_ptr[kspec]) + comb);
3659 
3660  if (zeroedPhase) {
3661  printf(" --- ** zp *** ");
3662  } else {
3663  printf(" --- ");
3664  }
3665  printf("%-24.24s", m_speciesName[kspec].c_str());
3666  printf(" %-3s", int2str(iphase).c_str());
3667  printf(" % -12.4e", mfValue);
3668  printf(" % -12.4e", m_SSfeSpecies[kspec] * RT);
3669  printf(" % -12.4e", log(mfValue) * RT);
3670  printf(" % -12.4e", log(actCoeff_ptr[kspec]) * RT);
3671  printf(" % -12.4e", elect * RT);
3672  printf(" % -12.4e", comb * RT);
3673  printf(" % -12.4e\n", total *RT);
3674  }
3675  printf(" ");
3676  writeline('-', 132);
3677 }
3678 
3679 #ifdef DEBUG_MODE
3680 void VCS_SOLVE::prneav() const
3681 {
3682  std::vector<double> eav(m_numElemConstraints, 0.0);
3683 
3684  for (size_t j = 0; j < m_numElemConstraints; ++j) {
3685  for (size_t i = 0; i < m_numSpeciesTot; ++i) {
3686  if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
3687  eav[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i];
3688  }
3689  }
3690  }
3691  bool kerr = false;
3692  plogf("--------------------------------------------------");
3693  plogf("ELEMENT ABUNDANCE VECTOR:\n");
3694  plogf(" Element Now Orignal Deviation Type\n");
3695  for (size_t j = 0; j < m_numElemConstraints; ++j) {
3696  plogf(" ");
3697  plogf("%-2.2s", (m_elementName[j]).c_str());
3698  plogf(" = %15.6E %15.6E %15.6E %3d\n",
3699  eav[j], m_elemAbundancesGoal[j], eav[j] - m_elemAbundancesGoal[j], m_elType[j]);
3700  if (m_elemAbundancesGoal[j] != 0.) {
3701  if (fabs(eav[j] - m_elemAbundancesGoal[j]) > m_elemAbundancesGoal[j] * 5.0e-9) {
3702  kerr = true;
3703  }
3704  } else {
3705  if (fabs(eav[j]) > 1.0e-10) {
3706  kerr = true;
3707  }
3708  }
3709  }
3710  if (kerr) {
3711  plogf("Element abundance check failure\n");
3712  }
3713  plogf("--------------------------------------------------");
3714  plogendl();
3715 }
3716 #endif
3717 
3718 double VCS_SOLVE::l2normdg(double dgLocal[]) const
3719 {
3720  if (m_numRxnRdc <= 0) {
3721  return 0.0;
3722  }
3723  double tmp = 0;
3724  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3725  size_t kspec = irxn + m_numComponents;
3726  if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR || m_speciesStatus[kspec] == VCS_SPECIES_MINOR ||
3727  dgLocal[irxn] < 0.0) {
3728  if (m_speciesStatus[kspec] != VCS_SPECIES_ZEROEDMS) {
3729  tmp += dgLocal[irxn] * dgLocal[irxn];
3730  }
3731  }
3732  }
3733  return std::sqrt(tmp / m_numRxnRdc);
3734 }
3735 
3736 double VCS_SOLVE::vcs_tmoles()
3737 {
3738  for (size_t i = 0; i < m_numPhases; i++) {
3739  m_tPhaseMoles_old[i] = TPhInertMoles[i];
3740  }
3741  for (size_t i = 0; i < m_numSpeciesTot; i++) {
3742  if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
3743  m_tPhaseMoles_old[m_phaseID[i]] += m_molNumSpecies_old[i];
3744  }
3745  }
3746  double sum = 0.0;
3747  for (size_t i = 0; i < m_numPhases; i++) {
3748  sum += m_tPhaseMoles_old[i];
3749  vcs_VolPhase* Vphase = m_VolPhaseList[i];
3750  if (m_tPhaseMoles_old[i] == 0.0) {
3751  Vphase->setTotalMoles(0.0);
3752  } else {
3753  Vphase->setTotalMoles(m_tPhaseMoles_old[i]);
3754  }
3755  }
3756  m_totalMolNum = sum;
3757  return m_totalMolNum;
3758 }
3759 
3760 #ifdef DEBUG_MODE
3761 void VCS_SOLVE::check_tmoles() const
3762 {
3763  double sum = 0.0;
3764  for (size_t i = 0; i < m_numPhases; i++) {
3765  double m_tPhaseMoles_old_a = TPhInertMoles[i];
3766 
3767  for (size_t k = 0; k < m_numSpeciesTot; k++) {
3768  if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM) {
3769  if (m_phaseID[k] == i) {
3770  m_tPhaseMoles_old_a += m_molNumSpecies_old[k];
3771  }
3772  }
3773  }
3774  sum += m_tPhaseMoles_old_a;
3775 
3776  double denom = m_tPhaseMoles_old[i]+ m_tPhaseMoles_old_a + 1.0E-19;
3777  if (!vcs_doubleEqual(m_tPhaseMoles_old[i]/denom, m_tPhaseMoles_old_a/denom)) {
3778  plogf("check_tmoles: we have found a problem with phase %d: %20.15g, %20.15g\n",
3779  i, m_tPhaseMoles_old[i], m_tPhaseMoles_old_a);
3780  }
3781  }
3782 }
3783 #endif
3784 
3785 void VCS_SOLVE::vcs_updateVP(const int vcsState)
3786 {
3787  for (size_t i = 0; i < m_numPhases; i++) {
3788  vcs_VolPhase* Vphase = m_VolPhaseList[i];
3789  if (vcsState == VCS_STATECALC_OLD) {
3790  Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
3791  VCS_DATA_PTR(m_molNumSpecies_old),
3792  VCS_DATA_PTR(m_tPhaseMoles_old));
3793  } else if (vcsState == VCS_STATECALC_NEW) {
3794  Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW,
3795  VCS_DATA_PTR(m_molNumSpecies_new),
3796  VCS_DATA_PTR(m_tPhaseMoles_new));
3797  } else if (DEBUG_MODE_ENABLED) {
3798  throw CanteraError("VCS_SOLVE::vcs_updateVP",
3799  "wrong stateCalc value: " + int2str(vcsState));
3800  }
3801  }
3802 }
3803 
3804 bool VCS_SOLVE::vcs_evaluate_speciesType()
3805 {
3806  m_numRxnMinorZeroed = 0;
3807  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3808  plogf(" --- Species Status decision is reevaluated: All species are minor except for:\n");
3809  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) {
3810  plogf(" --- Species Status decision is reevaluated");
3811  plogendl();
3812  }
3813  for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
3814  m_speciesStatus[kspec] = vcs_species_type(kspec);
3815  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 5) {
3816  plogf(" --- %-16s: ", m_speciesName[kspec].c_str());
3817  if (kspec < m_numComponents) {
3818  plogf("(COMP) ");
3819  } else {
3820  plogf(" ");
3821  }
3822  plogf(" %10.3g ", m_molNumSpecies_old[kspec]);
3823  const char* sString = vcs_speciesType_string(m_speciesStatus[kspec], 100);
3824  plogf("%s\n", sString);
3825 
3826  } else if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3827  if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
3828  switch (m_speciesStatus[kspec]) {
3829  case VCS_SPECIES_COMPONENT:
3830  break;
3831  case VCS_SPECIES_MAJOR:
3832  plogf(" --- Major Species : %-s\n", m_speciesName[kspec].c_str());
3833  break;
3834  case VCS_SPECIES_ZEROEDPHASE:
3835  plogf(" --- Purposely Zeroed-Phase Species (not in problem): %-s\n",
3836  m_speciesName[kspec].c_str());
3837  break;
3838  case VCS_SPECIES_ZEROEDMS:
3839  plogf(" --- Zeroed-MS Phase Species: %-s\n", m_speciesName[kspec].c_str());
3840  break;
3841  case VCS_SPECIES_ZEROEDSS:
3842  plogf(" --- Zeroed-SS Phase Species: %-s\n", m_speciesName[kspec].c_str());
3843  break;
3844  case VCS_SPECIES_DELETED:
3845  plogf(" --- Deleted-Small Species : %-s\n", m_speciesName[kspec].c_str());
3846  break;
3847  case VCS_SPECIES_ACTIVEBUTZERO:
3848  plogf(" --- Zeroed Species in an active MS phase (tmp): %-s\n",
3849  m_speciesName[kspec].c_str());
3850  break;
3851  case VCS_SPECIES_STOICHZERO:
3852  plogf(" --- Zeroed Species in an active MS phase (Stoich Constraint): %-s\n",
3853  m_speciesName[kspec].c_str());
3854  break;
3855  case VCS_SPECIES_INTERFACIALVOLTAGE:
3856  plogf(" --- InterfaceVoltage Species: %-s\n", m_speciesName[kspec].c_str());
3857  break;
3858  default:
3859  throw CanteraError("VCS_SOLVE::vcs_evaluate_speciesType",
3860  "Unknown type: " + int2str(m_speciesStatus[kspec]));
3861  }
3862  }
3863  }
3864  if (kspec >= m_numComponents) {
3865  if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
3866  ++m_numRxnMinorZeroed;
3867  }
3868  }
3869  }
3870  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3871  plogf(" ---");
3872  plogendl();
3873  }
3874  return (m_numRxnMinorZeroed >= m_numRxnRdc);
3875 }
3876 
3877 void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted,
3878  const int vcsState, const bool alterZeroedPhases)
3879 {
3880  int icase = 0;
3881  size_t irxnl = m_numRxnRdc;
3882  if (doDeleted) {
3883  irxnl = m_numRxnTot;
3884  }
3885 
3886  double* deltaGRxn;
3887  double* feSpecies;
3888  double* molNumSpecies;
3889  double* actCoeffSpecies;
3890  if (vcsState == VCS_STATECALC_NEW) {
3891  deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new);
3892  feSpecies = VCS_DATA_PTR(m_feSpecies_new);
3893  molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new);
3894  actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new);
3895  } else if (vcsState == VCS_STATECALC_OLD) {
3896  deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
3897  feSpecies = VCS_DATA_PTR(m_feSpecies_old);
3898  molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
3899  actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old);
3900  } else {
3901  throw CanteraError("VCS_SOLVE::vcs_deltag", "bad vcsState");
3902  }
3903 
3904  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
3905  plogf(" --- Subroutine vcs_deltag called for ");
3906  if (l < 0) {
3907  plogf("major noncomponents\n");
3908  } else if (l == 0) {
3909  plogf("all noncomponents\n");
3910  } else {
3911  plogf("minor noncomponents\n");
3912  }
3913  }
3914  /* ************************************************* */
3915  /* **** MAJORS and ZEROED SPECIES ONLY ************* */
3916  /* ************************************************* */
3917  if (l < 0) {
3918  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3919  size_t kspec = irxn + m_numComponents;
3920  if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
3921  icase = 0;
3922  deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
3923  double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3924  for (kspec = 0; kspec < m_numComponents; ++kspec) {
3925  deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3926  if (molNumSpecies[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
3927  icase = 1;
3928  }
3929  }
3930  if (icase) {
3931  deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3932  }
3933  }
3934  }
3935  } else if (l == 0) {
3936  /* ************************************************* */
3937  /* **** ALL REACTIONS ****************************** */
3938  /* ************************************************* */
3939  for (size_t irxn = 0; irxn < irxnl; ++irxn) {
3940  icase = 0;
3941  deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
3942  double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3943  for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3944  deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3945  if (molNumSpecies[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF &&
3946  dtmp_ptr[kspec] < 0.0) {
3947  icase = 1;
3948  }
3949  }
3950  if (icase) {
3951  deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3952  }
3953  }
3954  } else {
3955  /* ************************************************* */
3956  /* **** MINORS AND ZEROED SPECIES ****************** */
3957  /* ************************************************* */
3958  for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
3959  size_t kspec = irxn + m_numComponents;
3960  if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR) {
3961  icase = 0;
3962  deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
3963  double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
3964  for (kspec = 0; kspec < m_numComponents; ++kspec) {
3965  deltaGRxn[irxn] += dtmp_ptr[kspec] * feSpecies[kspec];
3966  if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF &&
3967  dtmp_ptr[kspec] < 0.0) {
3968  icase = 1;
3969  }
3970  }
3971  if (icase) {
3972  deltaGRxn[irxn] = std::max(0.0, deltaGRxn[irxn]);
3973  }
3974  }
3975  }
3976  }
3977  /* ************************************************* */
3978  /* **** MULTISPECIES PHASES WITH ZERO MOLES ******** */
3979  /* ************************************************* */
3980  /*
3981  * Massage the free energies for species with zero mole fractions
3982  * in multispecies phases. This section implements the
3983  * Equation 3.8-5 in Smith and Missen, p.59.
3984  * A multispecies phase will exist iff
3985  * 1 < sum_i(exp(-dg_i)/AC_i)
3986  * If DG is negative then that species wants to be reintroduced into
3987  * the calculation.
3988  * For small dg_i, the expression below becomes:
3989  * 1 - sum_i(exp(-dg_i)/AC_i) ~ sum_i((dg_i-1)/AC_i) + 1
3990  *
3991  * So, what we are doing here is equalizing all DG's in a multispecies
3992  * phase whose total mole number has already been zeroed out.
3993  * It must have to do with the case where a complete multispecies
3994  * phase is currently zeroed out. In that case, when one species
3995  * in that phase has a negative DG, then the phase should kick in.
3996  * This code section will cause that to happen, because a negative
3997  * DG will dominate the calculation of SDEL. Then, DG(I) for all
3998  * species in that phase will be forced to be equal and negative.
3999  * Thus, all species in that phase will come into being at the
4000  * same time.
4001  *
4002  * HKM -> The ratio of mole fractions at the reinstatement
4003  * time should be equal to the normalized weighting
4004  * of exp(-dg_i) / AC_i. This should be implemented.
4005  *
4006  * HKM -> There is circular logic here. ActCoeff depends on the
4007  * mole fractions of a phase that does not exist. In actuality
4008  * the proto-mole fractions should be selected from the
4009  * solution of a nonlinear problem with NsPhase unknowns
4010  *
4011  * X_i = exp(-dg[irxn]) / ActCoeff_i / denom
4012  *
4013  * where
4014  * denom = sum_i[ exp(-dg[irxn]) / ActCoeff_i ]
4015  *
4016  * This can probably be solved by successive iteration.
4017  * This should be implemented.
4018  */
4019  if (alterZeroedPhases && false) {
4020  for (size_t iph = 0; iph < m_numPhases; iph++) {
4021  bool lneed = false;
4022  vcs_VolPhase* Vphase = m_VolPhaseList[iph];
4023  if (! Vphase->m_singleSpecies) {
4024  double sum = 0.0;
4025  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
4026  size_t kspec = Vphase->spGlobalIndexVCS(k);
4027  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
4028  sum += molNumSpecies[kspec];
4029  }
4030  if (sum > 0.0) {
4031  break;
4032  }
4033  }
4034  if (sum == 0.0) {
4035  lneed = true;
4036  }
4037  }
4038 
4039  if (lneed) {
4040  double poly = 0.0;
4041  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
4042  size_t kspec = Vphase->spGlobalIndexVCS(k);
4043  // We may need to look at deltaGRxn for components!
4044  if (kspec >= m_numComponents) {
4045  size_t irxn = kspec - m_numComponents;
4046  deltaGRxn[irxn] = clip(deltaGRxn[irxn], -50.0, 50.0);
4047  poly += exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec];
4048  }
4049  }
4050  /*
4051  * Calculate deltaGRxn[] for each species in a zeroed multispecies phase.
4052  * All of the m_deltaGRxn_new[]'s will be equal. If deltaGRxn[] is
4053  * negative, then the phase will come back into existence.
4054  */
4055  for (size_t k = 0; k < Vphase->nSpecies(); k++) {
4056  size_t kspec = Vphase->spGlobalIndexVCS(k);
4057  if (kspec >= m_numComponents) {
4058  size_t irxn = kspec - m_numComponents;
4059  deltaGRxn[irxn] = 1.0 - poly;
4060  }
4061  }
4062 
4063  }
4064  }
4065  }
4066 
4067 
4068 #ifdef DEBUG_NOT
4069  for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
4070  checkFinite(deltaGRxn[irxn]);
4071  }
4072 #endif
4073 }
4074 
4075 void VCS_SOLVE::vcs_printDeltaG(const int stateCalc)
4076 {
4077  double* deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
4078  double* feSpecies = VCS_DATA_PTR(m_feSpecies_old);
4079  double* molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
4080  const double* tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
4081  const double* actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
4082  if (stateCalc == VCS_STATECALC_NEW) {
4083  deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new);
4084  feSpecies = VCS_DATA_PTR(m_feSpecies_new);
4085  molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new);
4086  actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
4087  tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
4088  }
4089  double RT = m_temperature * GasConstant;
4090  bool zeroedPhase = false;
4091  if (m_debug_print_lvl >= 2) {
4092  plogf(" --- DELTA_G TABLE Components:");
4093  for (size_t j = 0; j < m_numComponents; j++) {
4094  plogf(" %3d ", j);
4095  }
4096  plogf("\n --- Components Moles:");
4097  for (size_t j = 0; j < m_numComponents; j++) {
4098  plogf("%10.3g", m_molNumSpecies_old[j]);
4099  }
4100  plogf("\n --- NonComponent| Moles | ");
4101  for (size_t j = 0; j < m_numComponents; j++) {
4102  plogf("%-10.10s", m_speciesName[j].c_str());
4103  }
4104  plogf("\n");
4105  for (size_t i = 0; i < m_numRxnTot; i++) {
4106  plogf(" --- %3d ", m_indexRxnToSpecies[i]);
4107  plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str());
4108  if (m_speciesUnknownType[m_indexRxnToSpecies[i]] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
4109  plogf("| NA |");
4110  } else {
4111  plogf("|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
4112  }
4113  for (size_t j = 0; j < m_numComponents; j++) {
4114  plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,i));
4115  }
4116  plogf("\n");
4117  }
4118  plogf(" ");
4119  for (int i=0; i<77; i++) {
4120  plogf("-");
4121  }
4122  plogf("\n");
4123  }
4124 
4125  printf(" --- DeltaG Table (J/kmol) Name PhID MoleNum MolFR "
4126  " ElectrChemStar ElectrChem DeltaGStar DeltaG(Pred) Stability\n");
4127  printf(" ");
4128  writeline('-', 132);
4129 
4130  for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
4131 
4132  size_t irxn = npos;
4133  if (kspec >= m_numComponents) {
4134  irxn = kspec - m_numComponents;
4135  }
4136 
4137  double mfValue = 1.0;
4138  size_t iphase = m_phaseID[kspec];
4139  const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
4140  if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
4141  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
4142  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS)) {
4143  zeroedPhase = true;
4144  } else {
4145  zeroedPhase = false;
4146  }
4147  if (tPhMoles_ptr[iphase] > 0.0) {
4148  if (molNumSpecies[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
4149  mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tPhMoles_ptr[iphase];
4150  } else {
4151  mfValue = molNumSpecies[kspec] / tPhMoles_ptr[iphase];
4152  }
4153  } else {
4154  size_t klocal = m_speciesLocalPhaseIndex[kspec];
4155  mfValue = Vphase->moleFraction(klocal);
4156  }
4157  if (zeroedPhase) {
4158  printf(" --- ** zp *** ");
4159  } else {
4160  printf(" --- ");
4161  }
4162  double feFull = feSpecies[kspec];
4163  if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
4164  (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE)) {
4165  feFull += log(actCoeff_ptr[kspec]) + log(mfValue);
4166  }
4167  printf("%-24.24s", m_speciesName[kspec].c_str());
4168  printf(" %-3s", int2str(iphase).c_str());
4169  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
4170  printf(" NA ");
4171  } else {
4172  printf(" % -12.4e", molNumSpecies[kspec]);
4173  }
4174  printf(" % -12.4e", mfValue);
4175  printf(" % -12.4e", feSpecies[kspec] * RT);
4176  printf(" % -12.4e", feFull * RT);
4177  if (irxn != npos) {
4178  printf(" % -12.4e", deltaGRxn[irxn] * RT);
4179  printf(" % -12.4e", (deltaGRxn[irxn] + feFull - feSpecies[kspec]) * RT);
4180 
4181  if (deltaGRxn[irxn] < 0.0) {
4182  if (molNumSpecies[kspec] > 0.0) {
4183  printf(" growing");
4184  } else {
4185  printf(" stable");
4186  }
4187  } else if (deltaGRxn[irxn] > 0.0) {
4188  if (molNumSpecies[kspec] > 0.0) {
4189  printf(" shrinking");
4190  } else {
4191  printf(" unstable");
4192  }
4193  } else {
4194  printf(" balanced");
4195  }
4196  }
4197 
4198  printf(" \n");
4199  }
4200 
4201  printf(" ");
4202  writeline('-', 132);
4203 
4204 }
4205 
4206 void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
4207  const int stateCalc, const bool alterZeroedPhases)
4208 {
4209  double* feSpecies=0;
4210  double* deltaGRxn=0;
4211  double* actCoeffSpecies=0;
4212  if (stateCalc == VCS_STATECALC_NEW) {
4213  feSpecies = VCS_DATA_PTR(m_feSpecies_new);
4214  deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new);
4215  actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new);
4216  } else if (stateCalc == VCS_STATECALC_OLD) {
4217  feSpecies = VCS_DATA_PTR(m_feSpecies_old);
4218  deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
4219  actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_old);
4220  } else if (DEBUG_MODE_ENABLED) {
4221  throw CanteraError("VCS_SOLVE::vcs_deltag_Phase", "bad stateCalc");
4222  }
4223 
4224  size_t irxnl = m_numRxnRdc;
4225  if (doDeleted) {
4226  irxnl = m_numRxnTot;
4227  }
4228  vcs_VolPhase* vPhase = m_VolPhaseList[iphase];
4229 
4230  if (DEBUG_MODE_ENABLED && m_debug_print_lvl >= 2) {
4231  plogf(" --- Subroutine vcs_deltag_Phase called for phase %d\n",
4232  iphase);
4233  }
4234 
4235  /*
4236  * Single species Phase
4237  */
4238  if (vPhase->m_singleSpecies) {
4239  size_t kspec = vPhase->spGlobalIndexVCS(0);
4240  AssertThrowMsg(iphase == m_phaseID[kspec], "VCS_SOLVE::vcs_deltag_Phase",
4241  "index error");
4242  if (kspec >= m_numComponents) {
4243  size_t irxn = kspec - m_numComponents;
4244  deltaGRxn[irxn] = feSpecies[kspec];
4245  for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
4246  deltaGRxn[irxn] += m_stoichCoeffRxnMatrix(kcomp,irxn) * feSpecies[kcomp];
4247  }
4248  }
4249  }
4250  /*
4251  * Multispecies Phase
4252  */
4253  else {
4254  bool zeroedPhase = true;
4255 
4256  for (size_t irxn = 0; irxn < irxnl; ++irxn) {
4257  size_t kspec = m_indexRxnToSpecies[irxn];
4258  if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
4259  if (m_phaseID[kspec] == iphase) {
4260  if (m_molNumSpecies_old[kspec] > 0.0) {
4261  zeroedPhase = false;
4262  }
4263  deltaGRxn[irxn] = feSpecies[kspec];
4264  for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
4265  deltaGRxn[irxn] += m_stoichCoeffRxnMatrix(kcomp,irxn) * feSpecies[kcomp];
4266  }
4267  }
4268  }
4269  }
4270 
4271  /*
4272  * special section for zeroed phases
4273  */
4274  /* ************************************************* */
4275  /* **** MULTISPECIES PHASES WITH ZERO MOLES************ */
4276  /* ************************************************* */
4277  /*
4278  * Massage the free energies for species with zero mole fractions
4279  * in multispecies phases. This section implements the
4280  * Equation 3.8-5 in Smith and Missen, p.59.
4281  * A multispecies phase will exist iff
4282  * 1 < sum_i(exp(-dg_i)/AC_i)
4283  * If DG is negative then that species wants to be reintroduced into
4284  * the calculation.
4285  * For small dg_i, the expression below becomes:
4286  * 1 - sum_i(exp(-dg_i)/AC_i) ~ sum_i((dg_i-1)/AC_i) + 1
4287  *
4288  *
4289  * HKM -> The ratio of mole fractions at the reinstatement
4290  * time should be equal to the normalized weighting
4291  * of exp(-dg_i) / AC_i. This should be implemented.
4292  *
4293  * HKM -> There is circular logic here. ActCoeff depends on the
4294  * mole fractions of a phase that does not exist. In actuality
4295  * the proto-mole fractions should be selected from the
4296  * solution of a nonlinear problem with NsPhase unknowns
4297  *
4298  * X_i = exp(-dg[irxn]) / ActCoeff_i / denom
4299  *
4300  * where
4301  * denom = sum_i[ exp(-dg[irxn]) / ActCoeff_i ]
4302  *
4303  * This can probably be solved by successive iteration.
4304  * This should be implemented.
4305  */
4306  /*
4307  * Calculate dg[] for each species in a zeroed multispecies phase.
4308  * All of the dg[]'s will be equal. If dg[] is negative, then
4309  * the phase will come back into existence.
4310  */
4311  if (alterZeroedPhases) {
4312  if (zeroedPhase) {
4313  double phaseDG = 1.0;
4314  for (size_t irxn = 0; irxn < irxnl; ++irxn) {
4315  size_t kspec = m_indexRxnToSpecies[irxn];
4316  if (m_phaseID[kspec] == iphase) {
4317  deltaGRxn[irxn] = clip(deltaGRxn[irxn], -50.0, 50.0);
4318  phaseDG -= exp(-deltaGRxn[irxn])/actCoeffSpecies[kspec];
4319  }
4320  }
4321  /*
4322  * Overwrite the individual dg's with the phase DG.
4323  */
4324  for (size_t irxn = 0; irxn < irxnl; ++irxn) {
4325  size_t kspec = m_indexRxnToSpecies[irxn];
4326  if (m_phaseID[kspec] == iphase) {
4327  deltaGRxn[irxn] = 1.0 - phaseDG;
4328  }
4329  }
4330  }
4331  }
4332  }
4333 }
4334 
4335 void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2)
4336 {
4337  if (k1 == k2) {
4338  return;
4339  }
4340  if (DEBUG_MODE_ENABLED && (k1 >= m_numSpeciesTot ||
4341  k2 >= m_numSpeciesTot)) {
4342  plogf("vcs_switch_pos: ifunc = 0: inappropriate args: %d %d\n",
4343  k1, k2);
4344  }
4345  /*
4346  * Handle the index pointer in the phase structures first
4347  */
4348  vcs_VolPhase* pv1 = m_VolPhaseList[m_phaseID[k1]];
4349  vcs_VolPhase* pv2 = m_VolPhaseList[m_phaseID[k2]];
4350 
4351  size_t kp1 = m_speciesLocalPhaseIndex[k1];
4352  size_t kp2 = m_speciesLocalPhaseIndex[k2];
4353  AssertThrowMsg(pv1->spGlobalIndexVCS(kp1) == k1, "VCS_SOLVE::vcs_switch_pos",
4354  "Indexing error");
4355  AssertThrowMsg(pv2->spGlobalIndexVCS(kp2) == k2, "VCS_SOLVE::vcs_switch_pos",
4356  "Indexing error");
4357  pv1->setSpGlobalIndexVCS(kp1, k2);
4358  pv2->setSpGlobalIndexVCS(kp2, k1);
4359  std::swap(m_speciesName[k1], m_speciesName[k2]);
4360  std::swap(m_molNumSpecies_old[k1], m_molNumSpecies_old[k2]);
4361  std::swap(m_speciesUnknownType[k1], m_speciesUnknownType[k2]);
4362  std::swap(m_molNumSpecies_new[k1], m_molNumSpecies_new[k2]);
4363  std::swap(m_SSfeSpecies[k1], m_SSfeSpecies[k2]);
4364  std::swap(m_spSize[k1], m_spSize[k2]);
4365  std::swap(m_deltaMolNumSpecies[k1], m_deltaMolNumSpecies[k2]);
4366  std::swap(m_feSpecies_old[k1], m_feSpecies_old[k2]);
4367  std::swap(m_feSpecies_new[k1], m_feSpecies_new[k2]);
4368  std::swap(m_SSPhase[k1], m_SSPhase[k2]);
4369  std::swap(m_phaseID[k1], m_phaseID[k2]);
4370  std::swap(m_speciesMapIndex[k1], m_speciesMapIndex[k2]);
4371  std::swap(m_speciesLocalPhaseIndex[k1], m_speciesLocalPhaseIndex[k2]);
4372  std::swap(m_actConventionSpecies[k1], m_actConventionSpecies[k2]);
4373  std::swap(m_lnMnaughtSpecies[k1], m_lnMnaughtSpecies[k2]);
4374  std::swap(m_actCoeffSpecies_new[k1], m_actCoeffSpecies_new[k2]);
4375  std::swap(m_actCoeffSpecies_old[k1], m_actCoeffSpecies_old[k2]);
4376  std::swap(m_wtSpecies[k1], m_wtSpecies[k2]);
4377  std::swap(m_chargeSpecies[k1], m_chargeSpecies[k2]);
4378  std::swap(m_speciesThermoList[k1], m_speciesThermoList[k2]);
4379  std::swap(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2]);
4380 
4381  for (size_t j = 0; j < m_numElemConstraints; ++j) {
4382  std::swap(m_formulaMatrix(k1,j), m_formulaMatrix(k2,j));
4383  }
4384  if (m_useActCoeffJac && k1 != k2) {
4385  for (size_t i = 0; i < m_numSpeciesTot; i++) {
4386  std::swap(m_np_dLnActCoeffdMolNum(k1,i), m_np_dLnActCoeffdMolNum(k2,i));
4387  }
4388  for (size_t i = 0; i < m_numSpeciesTot; i++) {
4389  std::swap(m_np_dLnActCoeffdMolNum(i,k1), m_np_dLnActCoeffdMolNum(i,k2));
4390  }
4391  }
4392  std::swap(m_speciesStatus[k1], m_speciesStatus[k2]);
4393  /*
4394  * Handle the index pointer in the phase structures
4395  */
4396 
4397 
4398  if (ifunc) {
4399  /*
4400  * Find the Rxn indices corresponding to the two species
4401  */
4402  size_t i1 = k1 - m_numComponents;
4403  size_t i2 = k2 - m_numComponents;
4404  if (DEBUG_MODE_ENABLED && (i1 > m_numRxnTot || i2 >= m_numRxnTot)) {
4405  plogf("switch_pos: ifunc = 1: inappropriate noncomp values: %d %d\n",
4406  i1 , i2);
4407  }
4408  for (size_t j = 0; j < m_numComponents; ++j) {
4409  std::swap(m_stoichCoeffRxnMatrix(j,i1), m_stoichCoeffRxnMatrix(j,i2));
4410  }
4411  std::swap(m_scSize[i1], m_scSize[i2]);
4412  for (size_t iph = 0; iph < m_numPhases; iph++) {
4413  std::swap(m_deltaMolNumPhase(iph,i1), m_deltaMolNumPhase(iph,i2));
4414  std::swap(m_phaseParticipation(iph,i1),
4415  m_phaseParticipation(iph,i2));
4416  }
4417  std::swap(m_deltaGRxn_new[i1], m_deltaGRxn_new[i2]);
4418  std::swap(m_deltaGRxn_old[i1], m_deltaGRxn_old[i2]);
4419  std::swap(m_deltaGRxn_tmp[i1], m_deltaGRxn_tmp[i2]);
4420  }
4421 }
4422 
4423 double VCS_SOLVE::vcs_birthGuess(const int kspec)
4424 {
4425  size_t irxn = kspec - m_numComponents;
4426  double dx = 0.0;
4427  if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
4428  return dx;
4429  }
4430  double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
4431 
4432  // Check to make sure that species is zero in the solution vector
4433  // If it isn't, we don't know what's happening
4434  AssertThrowMsg(m_molNumSpecies_old[kspec] == 0.0,
4435  "VCS_SOLVE::vcs_birthGuess", "we shouldn't be here");
4436  int ss = m_SSPhase[kspec];
4437  if (!ss) {
4438  /*
4439  * Logic to handle species in multiple species phases
4440  * we cap the moles here at 1.0E-15 kmol.
4441  */
4442  bool soldel_ret;
4443  double dxm = vcs_minor_alt_calc(kspec, irxn, &soldel_ret);
4444  dx = std::min(w_kspec + dxm, 1e-15);
4445  } else {
4446  /*
4447  * Logic to handle single species phases
4448  * There is no real way to estimate the moles. So
4449  * we set it to a small number.
4450  */
4451  dx = 1.0E-30;
4452  }
4453 
4454  /*
4455  * Check to see if the current value of the components
4456  * allow the dx just estimated.
4457  * If we are in danger of zeroing a component,
4458  * only go 1/3 the way to zeroing the component with
4459  * this dx. Note, this may mean that dx= 0 coming
4460  * back from this routine. This evaluation should
4461  * be respected.
4462  */
4463  double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
4464  for (size_t j = 0; j < m_numComponents; ++j) {
4465  // Only loop over element constraints that involve positive def. constraints
4466  if (m_speciesUnknownType[j] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
4467  if (m_molNumSpecies_old[j] > 0.0) {
4468  double tmp = sc_irxn[j] * dx;
4469  if (3.0*(-tmp) > m_molNumSpecies_old[j]) {
4470  dx = std::min(dx, - 0.3333* m_molNumSpecies_old[j] / sc_irxn[j]);
4471  }
4472  }
4473  if (m_molNumSpecies_old[j] <= 0.0) {
4474  if (sc_irxn[j] < 0.0) {
4475  dx = 0.0;
4476  }
4477  }
4478  }
4479  }
4480  return dx;
4481 }
4482 
4483 void VCS_SOLVE::vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc)
4484 {
4485  if (!upToDate) {
4486  for (size_t iph = 0; iph < m_numPhases; iph++) {
4487  m_VolPhaseList[iph]->setMolesOutOfDate(stateCalc);
4488  }
4489  } else {
4490  for (size_t iph = 0; iph < m_numPhases; iph++) {
4491  m_VolPhaseList[iph]->setMolesCurrent(stateCalc);
4492  }
4493  }
4494 }
4495 
4496 void VCS_SOLVE::vcs_setFlagsVolPhase(const size_t iph, const bool upToDate,
4497  const int stateCalc)
4498 {
4499  if (!upToDate) {
4500  m_VolPhaseList[iph]->setMolesOutOfDate(stateCalc);
4501  } else {
4502  m_VolPhaseList[iph]->setMolesCurrent(stateCalc);
4503  }
4504 }
4505 
4506 void VCS_SOLVE::vcs_updateMolNumVolPhases(const int stateCalc)
4507 {
4508  for (size_t iph = 0; iph < m_numPhases; iph++) {
4509  m_VolPhaseList[iph]->updateFromVCS_MoleNumbers(stateCalc);
4510  }
4511 }
4512 
4513 }
VCS_PHASE_EXIST_NO
#define VCS_PHASE_EXIST_NO
Phase doesn't currently exist in the mixture.
Definition: vcs_defs.h:225
Cantera::vcs_doubleEqual
bool vcs_doubleEqual(double d1, double d2)
Simple routine to check whether two doubles are equal up to roundoff error.
Definition: vcs_util.cpp:185
Cantera::int2str
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
Cantera::vcs_VolPhase::m_singleSpecies
bool m_singleSpecies
If true, this phase consists of a single species.
Definition: vcs_VolPhase.h:584
Cantera::checkFinite
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
Definition: checkFinite.cpp:43
Cantera::vcs_VolPhase::updateFromVCS_MoleNumbers
void updateFromVCS_MoleNumbers(const int stateCalc)
Update the moles within the phase, if necessary.
Definition: vcs_VolPhase.cpp:578
VCS_SPECIES_ZEROEDMS
#define VCS_SPECIES_ZEROEDMS
Species lies in a multicomponent phase with concentration zero.
Definition: vcs_defs.h:150
VCS_SPECIES_MINOR
#define VCS_SPECIES_MINOR
Species is a major species.
Definition: vcs_defs.h:135
clockWC.h
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
VCS_PHASE_EXIST_ALWAYS
#define VCS_PHASE_EXIST_ALWAYS
Always exists because it contains inerts which can't exist in any other phase.
Definition: vcs_defs.h:212
Cantera::vcs_VolPhase::setMolesFromVCSCheck
void setMolesFromVCSCheck(const int vcsStateStatus, const double *molesSpeciesVCS, const double *const TPhMoles)
Set the moles within the phase.
Definition: vcs_VolPhase.cpp:558
VCS_PHASE_EXIST_ZEROEDPHASE
#define VCS_PHASE_EXIST_ZEROEDPHASE
Phase currently is zeroed due to a programmatic issue.
Definition: vcs_defs.h:231
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
VCS_DATA_PTR
#define VCS_DATA_PTR(vvv)
Points to the data in a std::vector<> object.
Definition: vcs_internal.h:18
Cantera::vcs_VolPhase::nSpecies
size_t nSpecies() const
Return the number of species in the phase.
Definition: vcs_VolPhase.cpp:1296
Cantera::clip
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:270
Cantera::clockWC
The class provides the wall clock timer in seconds.
Definition: clockWC.h:51
VCS_SPECIES_INTERFACIALVOLTAGE
#define VCS_SPECIES_INTERFACIALVOLTAGE
Species refers to an electron in the metal.
Definition: vcs_defs.h:173
Cantera::vcs_VolPhase::PhaseName
std::string PhaseName
String name for the phase.
Definition: vcs_VolPhase.h:696
VCS_PHASE_EXIST_YES
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
Definition: vcs_defs.h:215
Cantera::vcs_VolPhase::setExistence
void setExistence(const int existence)
Set the existence flag in the object.
Definition: vcs_VolPhase.cpp:1024
Cantera::vcs_VolPhase::electricPotential
double electricPotential() const
Returns the electric field of the phase.
Definition: vcs_VolPhase.cpp:641
AssertThrowMsg
#define AssertThrowMsg(expr, procedure, message)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:301
Cantera::vcs_VolPhase::setMolesFromVCS
void setMolesFromVCS(const int stateCalc, const double *molesSpeciesVCS=0)
Set the moles within the phase.
Definition: vcs_VolPhase.cpp:471
Cantera::vector_int
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:159
vcs_solve.h
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
vcs_VolPhase.h
Header for the object representing each phase within vcs.
VCS_SPECIES_ZEROEDSS
#define VCS_SPECIES_ZEROEDSS
Species is a SS phase, that is currently zeroed out.
Definition: vcs_defs.h:157
VCS_SPECIES_TYPE_MOLNUM
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:343
VCS_ELEM_TYPE_ABSPOS
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:290
VCS_SPECIES_ZEROEDPHASE
#define VCS_SPECIES_ZEROEDPHASE
Species lies in a multicomponent phase that is zeroed atm.
Definition: vcs_defs.h:182
VCS_DELETE_MINORSPECIES_CUTOFF
#define VCS_DELETE_MINORSPECIES_CUTOFF
Cutoff relative mole number value, below which species are deleted from the equilibrium problem...
Definition: vcs_defs.h:57
VCS_SUCCESS
#define VCS_SUCCESS
Definition: vcs_defs.h:21
Cantera::clockWC::secondsWC
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition: clockWC.cpp:35
Cantera::fp2str
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:28
VCS_DELETE_PHASE_CUTOFF
#define VCS_DELETE_PHASE_CUTOFF
Cutoff relative moles below which a phase is deleted from the equilibrium problem.
Definition: vcs_defs.h:69
VCS_STATECALC_NEW
#define VCS_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
Definition: vcs_defs.h:368
Cantera::vcs_print_stringTrunc
void vcs_print_stringTrunc(const char *str, size_t space, int alignment)
Print a string within a given space limit.
Definition: vcs_util.cpp:154
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
VCS_SPECIES_ACTIVEBUTZERO
#define VCS_SPECIES_ACTIVEBUTZERO
Species lies in a multicomponent phase that is active, but species concentration is zero...
Definition: vcs_defs.h:191
Cantera::vcs_VolPhase::sendToVCS_ActCoeff
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
Definition: vcs_VolPhase.cpp:589
Cantera::vcs_VolPhase::exists
int exists() const
int indicating whether the phase exists or not
Definition: vcs_VolPhase.cpp:1019
VCS_SPECIES_MAJOR
#define VCS_SPECIES_MAJOR
Species is a major species.
Definition: vcs_defs.h:128
VCS_RELDELETE_SPECIES_CUTOFF
#define VCS_RELDELETE_SPECIES_CUTOFF
Cutoff relative mole fraction value, below which species are deleted from the equilibrium problem...
Definition: vcs_defs.h:51
VCS_SPECIES_STOICHZERO
#define VCS_SPECIES_STOICHZERO
Species lies in a multicomponent phase that is active, but species concentration is zero due to stoic...
Definition: vcs_defs.h:202
Cantera::vcs_VolPhase
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:86
VCS_DIMENSIONAL_G
#define VCS_DIMENSIONAL_G
dimensioned
Definition: vcs_defs.h:102
VCS_STATECALC_OLD
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:365
plogendl
#define plogendl()
define this Cantera function to replace cout << endl;
Definition: vcs_internal.h:31
Cantera::vcs_VolPhase::setSpGlobalIndexVCS
void setSpGlobalIndexVCS(const size_t spIndex, const size_t spGlobalIndex)
set the Global VCS index of the kth species in the phase
Definition: vcs_VolPhase.cpp:1059
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
Cantera::vcs_speciesType_string
const char * vcs_speciesType_string(int speciesStatus, int length)
Returns a const char string representing the type of the species given by the first argument...
Definition: vcs_util.cpp:96
plogf
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:24
VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:352
VCS_SPECIES_COMPONENT
#define VCS_SPECIES_COMPONENT
Species is a component which can be nonzero.
Definition: vcs_defs.h:121
Cantera::vcs_VolPhase::spGlobalIndexVCS
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
Definition: vcs_VolPhase.cpp:1054
VCS_SPECIES_DELETED
#define VCS_SPECIES_DELETED
Species has such a small mole fraction it is deleted even though its phase may possibly exist...
Definition: vcs_defs.h:166
ctexceptions.h
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
Cantera::vcs_VolPhase::setTotalMoles
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.
Definition: vcs_VolPhase.cpp:918
Generated by   doxygen 1.8.6