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