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