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