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