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