Cantera  3.0.0
Loading...
Searching...
No Matches
vcs_solve.cpp
Go to the documentation of this file.
1/*!
2 * @file vcs_solve.cpp Implementation file for the internal class that holds
3 * the problem.
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
18
19using namespace std;
20
21namespace Cantera
22{
23
24namespace {
25
26void printProgress(const vector<string>& spName, const vector<double>& soln,
27 const vector<double>& ff)
28{
29 double sum = 0.0;
30 plogf(" --- Summary of current progress:\n");
31 plogf(" --- Name Moles - SSGibbs \n");
32 plogf(" -------------------------------------------------------------------------------------\n");
33 for (size_t k = 0; k < soln.size(); k++) {
34 plogf(" --- %20s %12.4g - %12.4g\n", spName[k], soln[k], ff[k]);
35 sum += soln[k] * ff[k];
36 }
37 plogf(" --- Total sum to be minimized = %g\n", sum);
38}
39
40} // anonymous namespace
41
43
44VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
45 m_printLvl(printLvl),
46 m_mix(mphase),
47 m_nsp(mphase->nSpecies()),
48 m_numSpeciesRdc(mphase->nSpecies()),
49 m_numPhases(mphase->nPhases()),
50 m_temperature(mphase->temperature()),
51 m_pressurePA(mphase->pressure()),
52 m_totalVol(mphase->volume()),
53 m_Faraday_dim(Faraday / (m_temperature * GasConstant))
54{
56 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
57 m_speciesThermoList[kspec] = make_unique<VCS_SPECIES_THERMO>();
58 }
59
60 string ser = "VCS_SOLVE: ERROR:\n\t";
61 if (m_nsp <= 0) {
62 plogf("%s Number of species is nonpositive\n", ser);
63 throw CanteraError("VCS_SOLVE::VCS_SOLVE", ser +
64 " Number of species is nonpositive\n");
65 }
66 if (m_numPhases <= 0) {
67 plogf("%s Number of phases is nonpositive\n", ser);
68 throw CanteraError("VCS_SOLVE::VCS_SOLVE", ser +
69 " Number of phases is nonpositive\n");
70 }
71
72 /*
73 * We will initialize sc[] to note the fact that it needs to be
74 * filled with meaningful information.
75 */
76 m_scSize.resize(m_nsp, 0.0);
77 m_spSize.resize(m_nsp, 1.0);
78 m_SSfeSpecies.resize(m_nsp, 0.0);
79 m_feSpecies_new.resize(m_nsp, 0.0);
80 m_molNumSpecies_old.resize(m_nsp, 0.0);
84 m_phasePhi.resize(m_numPhases, 0.0);
85 m_molNumSpecies_new.resize(m_nsp, 0.0);
86 m_deltaGRxn_new.resize(m_nsp, 0.0);
87 m_deltaGRxn_old.resize(m_nsp, 0.0);
88 m_deltaGRxn_Deficient.resize(m_nsp, 0.0);
89 m_deltaGRxn_tmp.resize(m_nsp, 0.0);
90 m_deltaMolNumSpecies.resize(m_nsp, 0.0);
91 m_feSpecies_old.resize(m_nsp, 0.0);
92 m_tPhaseMoles_old.resize(m_numPhases, 0.0);
93 m_tPhaseMoles_new.resize(m_numPhases, 0.0);
94 m_deltaPhaseMoles.resize(m_numPhases, 0.0);
95 m_TmpPhase.resize(m_numPhases, 0.0);
96 m_TmpPhase2.resize(m_numPhases, 0.0);
97 TPhInertMoles.resize(m_numPhases, 0.0);
98
99 // ind[] is an index variable that keep track of solution vector rotations.
100 m_speciesMapIndex.resize(m_nsp, 0);
102
103 // ir[] is an index vector that keeps track of the irxn to species mapping.
104 // We can't fill it in until we know the number of c components in the
105 // problem
106 m_indexRxnToSpecies.resize(m_nsp, 0);
107
108 // Initialize all species to be major species
110
111 m_SSPhase.resize(2*m_nsp, 0);
112 m_phaseID.resize(m_nsp, 0);
113 m_speciesName.resize(m_nsp);
114
115 // space for activity coefficients for all species. Set it equal to one.
116 m_actConventionSpecies.resize(m_nsp, 0);
118 m_lnMnaughtSpecies.resize(m_nsp, 0.0);
119 m_actCoeffSpecies_new.resize(m_nsp, 1.0);
120 m_actCoeffSpecies_old.resize(m_nsp, 1.0);
121 m_wtSpecies.resize(m_nsp, 0.0);
122 m_chargeSpecies.resize(m_nsp, 0.0);
123
124 // Phase Info
126 for (size_t iph = 0; iph < m_numPhases; iph++) {
127 m_VolPhaseList[iph] = make_unique<vcs_VolPhase>(this);
128 }
129
130 // For Future expansion
131 m_useActCoeffJac = true;
132 if (m_useActCoeffJac) {
134 }
135
136 m_PMVolumeSpecies.resize(m_nsp, 0.0);
137
138 // counters kept within vcs
139 m_VCount = new VCS_COUNTERS();
141
142 if (vcs_timing_print_lvl == 0) {
144 }
145
146 VCS_SPECIES_THERMO* ts_ptr = 0;
147
148 // Loop over the phases, transferring pertinent information
149 int kT = 0;
150 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
151 // Get the ThermoPhase object - assume volume phase
152 ThermoPhase* tPhase = &mphase->phase(iphase);
153 size_t nelem = tPhase->nElements();
154
155 // Query Cantera for the equation of state type of the current phase.
156 string eos = tPhase->type();
157 bool gasPhase = (eos == "ideal-gas");
158
159 // Find out the number of species in the phase
160 size_t nSpPhase = tPhase->nSpecies();
161 // Find out the name of the phase
162 string phaseName = tPhase->name();
163
164 // Call the basic vcs_VolPhase creation routine.
165 // Properties set here:
166 // ->PhaseNum = phase number in the thermo problem
167 // ->GasPhase = Boolean indicating whether it is a gas phase
168 // ->NumSpecies = number of species in the phase
169 // ->TMolesInert = Inerts in the phase = 0.0 for cantera
170 // ->PhaseName = Name of the phase
171 vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
172 VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
173 VolPhase->m_gasPhase = gasPhase;
174
175 // Tell the vcs_VolPhase pointer about cantera
176 VolPhase->setPtrThermoPhase(tPhase);
177 VolPhase->setTotalMoles(0.0);
178
179 // Set the electric potential of the volume phase from the
180 // ThermoPhase object's value.
181 VolPhase->setElectricPotential(tPhase->electricPotential());
182
183 // Query the ThermoPhase object to find out what convention
184 // it uses for the specification of activity and Standard State.
185 VolPhase->p_activityConvention = tPhase->activityConvention();
186
187 // Assign the value of eqn of state. Handle conflicts here.
188 if (eos == "ideal-gas") {
189 VolPhase->m_eqnState = VCS_EOS_IDEAL_GAS;
190 } else if (eos == "fixed-stoichiometry") {
191 VolPhase->m_eqnState = VCS_EOS_STOICH_SUB;
192 } else if (eos == "ideal-condensed") {
193 VolPhase->m_eqnState = VCS_EOS_IDEAL_SOLN;
194 } else if (tPhase->nDim() != 3) {
195 throw CanteraError("VCS_SOLVE::VCS_SOLVE",
196 "Surface/edge phase not handled yet.");
197 } else {
198 if (m_printLvl > 1) {
199 writelog("Unknown Cantera EOS to VCSnonideal: '{}'\n", eos);
200 }
201 VolPhase->m_eqnState = VCS_EOS_UNK_CANTERA;
202 }
203
204 // Transfer all of the element information from the ThermoPhase object
205 // to the vcs_VolPhase object. Also decide whether we need a new charge
206 // neutrality element in the phase to enforce a charge neutrality
207 // constraint. We also decide whether this is a single species phase
208 // with the voltage being the independent variable setting the chemical
209 // potential of the electrons.
210 VolPhase->transferElementsFM(tPhase);
211
212 // Combine the element information in the vcs_VolPhase
213 // object into the vprob object.
214 addPhaseElements(VolPhase);
216
217 // Loop through each species in the current phase
218 for (size_t k = 0; k < nSpPhase; k++) {
219 // Obtain the molecular weight of the species from the
220 // ThermoPhase object
221 m_wtSpecies[kT] = tPhase->molecularWeight(k);
222
223 // Obtain the charges of the species from the ThermoPhase object
224 m_chargeSpecies[kT] = tPhase->charge(k);
225
226 // Set the phaseid of the species
227 m_phaseID[kT] = iphase;
228
229 // Transfer the type of unknown
230 m_speciesUnknownType[kT] = VolPhase->speciesUnknownType(k);
232 // Set the initial number of kmoles of the species
233 m_molNumSpecies_old[kT] = mphase->speciesMoles(kT);
236 } else {
237 throw CanteraError("VCS_SOLVE::VCS_SOLVE",
238 "Unknown species type: {}", m_speciesUnknownType[kT]);
239 }
240
241 // Transfer the species information from the
242 // volPhase structure to the VPROB structure
243 // This includes:
244 // FormulaMatrix[][]
245 // VolPhase->IndSpecies[]
246 addOnePhaseSpecies(VolPhase, k, kT);
247
248 // Get a pointer to the thermo object
249 ts_ptr = m_speciesThermoList[kT].get();
250
251 // Fill in the vcs_SpeciesProperty structure
252 vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
253 sProp->NumElements = m_nelem;
254 sProp->SpName = mphase->speciesName(kT);
255 sProp->SpeciesThermo = ts_ptr;
256 sProp->WtSpecies = tPhase->molecularWeight(k);
257 sProp->FormulaMatrixCol.resize(m_nelem, 0.0);
258 for (size_t e = 0; e < m_nelem; e++) {
259 sProp->FormulaMatrixCol[e] = m_formulaMatrix(kT,e);
260 }
261 sProp->Charge = tPhase->charge(k);
262 sProp->SurfaceSpecies = false;
263 sProp->VolPM = 0.0;
264
265 // Transfer the thermo specification of the species
266 // vsolve->SpeciesThermo[]
267
268 // Add lookback connectivity into the thermo object first
269 ts_ptr->IndexPhase = iphase;
270 ts_ptr->IndexSpeciesPhase = k;
271 ts_ptr->OwningPhase = VolPhase;
272
273 // get a reference to the Cantera species thermo.
274 MultiSpeciesThermo& sp = tPhase->speciesThermo();
275
276 int spType = sp.reportType(k);
277 if (spType == SIMPLE) {
278 double c[4];
279 double minTemp, maxTemp, refPressure;
280 sp.reportParams(k, spType, c, minTemp, maxTemp, refPressure);
281 ts_ptr->SS0_Model = VCS_SS0_CONSTANT;
282 ts_ptr->SS0_T0 = c[0];
283 ts_ptr->SS0_H0 = c[1];
284 ts_ptr->SS0_S0 = c[2];
285 ts_ptr->SS0_Cp0 = c[3];
286 if (gasPhase) {
287 ts_ptr->SSStar_Model = VCS_SSSTAR_IDEAL_GAS;
289 } else {
290 ts_ptr->SSStar_Model = VCS_SSSTAR_CONSTANT;
291 ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
292 }
293 } else {
294 if (m_printLvl > 2) {
295 plogf("vcs_Cantera_convert: Species Type %d not known \n",
296 spType);
297 }
298 ts_ptr->SS0_Model = VCS_SS0_NOTHANDLED;
299 ts_ptr->SSStar_Model = VCS_SSSTAR_NOTHANDLED;
300 }
301
302 // Transfer the Volume Information -> NEEDS WORK
303 if (gasPhase) {
305 ts_ptr->SSStar_Vol0 = 82.05 * 273.15 / 1.0;
306 } else {
307 ts_ptr->SSStar_Vol_Model = VCS_SSVOL_CONSTANT;
308 ts_ptr->SSStar_Vol0 = 0.0;
309 }
310 kT++;
311 }
312
314
315 // Now, calculate a sample naught Gibbs free energy calculation
316 // at the specified temperature.
317 for (size_t k = 0; k < nSpPhase; k++) {
318 vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
319 ts_ptr = sProp->SpeciesThermo;
320 ts_ptr->SS0_feSave = VolPhase->G0_calc_one(k)/ GasConstant;
321 ts_ptr->SS0_TSave = m_temperature;
322 }
323 }
324
325 // Transfer initial element abundances based on the species mole numbers
326 for (size_t j = 0; j < m_nelem; j++) {
327 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
330 }
331 }
332 if (m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10) {
333 m_elemAbundancesGoal[j] = 0.0;
334 }
335 }
336
337 // Printout the species information: PhaseID's and mole nums
338 if (m_printLvl > 1) {
339 writeline('=', 80, true, true);
340 writeline('=', 16, false);
341 plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
342 writeline('=', 20);
343 writeline('=', 80);
344 plogf(" Phase IDs of species\n");
345 plogf(" species phaseID phaseName ");
346 plogf(" Initial_Estimated_kMols\n");
347 for (size_t i = 0; i < m_nsp; i++) {
348 size_t iphase = m_phaseID[i];
349
350 vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
351 plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase,
352 VolPhase->PhaseName.c_str());
354 plogf(" Volts = %-10.5g\n", m_molNumSpecies_old[i]);
355 } else {
356 plogf(" %-10.5g\n", m_molNumSpecies_old[i]);
357 }
358 }
359
360 // Printout of the Phase structure information
361 writeline('-', 80, true, true);
362 plogf(" Information about phases\n");
363 plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
364 plogf(" TMolesInert Tmoles(kmol)\n");
365
366 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
367 vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
368 plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
369 VolPhase->VP_ID_, VolPhase->m_singleSpecies,
370 VolPhase->m_gasPhase, VolPhase->eos_name(),
371 VolPhase->nSpecies(), VolPhase->totalMolesInert());
372 plogf("%16e\n", VolPhase->totalMoles());
373 }
374
375 writeline('=', 80, true, true);
376 writeline('=', 16, false);
377 plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
378 writeline('=', 20);
379 writeline('=', 80);
380 plogf("\n");
381 }
382
383 // TPhInertMoles[] -> must be copied over here
384 for (size_t iph = 0; iph < m_numPhases; iph++) {
385 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
386 TPhInertMoles[iph] = Vphase->totalMolesInert();
387 }
388
389 // m_speciesIndexVector[] is an index variable that keep track of solution
390 // vector rotations.
391 for (size_t i = 0; i < m_nsp; i++) {
392 m_speciesMapIndex[i] = i;
393 }
394
395 // IndEl[] is an index variable that keep track of element vector rotations.
396 for (size_t i = 0; i < m_nelem; i++) {
397 m_elementMapIndex[i] = i;
398 }
399
400 // Fill in the species to phase mapping. Check for bad values at the same
401 // time.
402 vector<size_t> numPhSp(m_numPhases, 0);
403 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
404 size_t iph = m_phaseID[kspec];
405 if (iph >= m_numPhases) {
406 throw CanteraError("VCS_SOLVE::VCS_SOLVE",
407 "Species to Phase Mapping, PhaseID, has a bad value\n"
408 "\tm_phaseID[{}] = {}\n"
409 "Allowed values: 0 to {}", kspec, iph, m_numPhases - 1);
410 }
411 m_phaseID[kspec] = m_phaseID[kspec];
412 m_speciesLocalPhaseIndex[kspec] = numPhSp[iph];
413 numPhSp[iph]++;
414 }
415 for (size_t iph = 0; iph < m_numPhases; iph++) {
416 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
417 if (numPhSp[iph] != Vphase->nSpecies()) {
418 throw CanteraError("VCS_SOLVE::VCS_SOLVE",
419 "Number of species in phase {}, {}, doesn't match ({} != {}) [vphase = {}]",
420 ser, iph, Vphase->PhaseName, numPhSp[iph], Vphase->nSpecies(), (size_t) Vphase);
421 }
422 }
423
424 for (size_t i = 0; i < m_nelem; i++) {
426 if (m_elemAbundancesGoal[i] != 0.0) {
427 if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) {
428 throw CanteraError("VCS_SOLVE::VCS_SOLVE",
429 "Charge neutrality condition {} is signicantly "
430 "nonzero, {}. Giving up",
432 } else {
433 if (m_debug_print_lvl >= 2) {
434 plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n",
436 }
437 m_elemAbundancesGoal[i] = 0.0;
438 }
439 }
440 }
441 }
442
443 // Copy over the species names
444 for (size_t i = 0; i < m_nsp; i++) {
445 m_speciesName[i] = m_mix->speciesName(i);
446 }
447
448 // Specify the Activity Convention information
449 for (size_t iph = 0; iph < m_numPhases; iph++) {
450 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
452 if (Vphase->p_activityConvention != 0) {
453 // We assume here that species 0 is the solvent. The solvent isn't
454 // on a unity activity basis The activity for the solvent assumes
455 // that the it goes to one as the species mole fraction goes to one;
456 // that is, it's really on a molarity framework. So
457 // SpecLnMnaught[iSolvent] = 0.0, and the loop below starts at 1,
458 // not 0.
459 size_t iSolvent = Vphase->spGlobalIndexVCS(0);
460 double mnaught = m_wtSpecies[iSolvent] / 1000.;
461 for (size_t k = 1; k < Vphase->nSpecies(); k++) {
462 size_t kspec = Vphase->spGlobalIndexVCS(k);
464 m_lnMnaughtSpecies[kspec] = log(mnaught);
465 }
466 }
467 }
468
469}
470
471VCS_SOLVE::~VCS_SOLVE()
472{
474}
475
476int VCS_SOLVE::vcs(int ipr, int ip1, int maxit)
477{
478 clockWC tickTock;
479
480 // This function is called to copy the public data and the current
481 // problem specification into the current object's data structure.
483
485
486 // Prep the problem data
487 // - adjust the identity of any phases
488 // - determine the number of components in the problem
489 int retn = vcs_prep(ip1);
490 if (retn != 0) {
491 plogf("vcs_prep_oneTime returned a bad status, %d: bailing!\n",
492 retn);
493 return retn;
494 }
495
496 // Once we have defined the global internal data structure defining the
497 // problem, then we go ahead and solve the problem.
498 //
499 // (right now, all we do is solve fixed T, P problems. Methods for other
500 // problem types will go in at this level. For example, solving for
501 // fixed T, V problems will involve a 2x2 Newton's method, using loops
502 // over vcs_TP() to calculate the residual and Jacobian)
503 int iconv = vcs_TP(ipr, ip1, maxit, m_temperature, m_pressurePA);
504
505 // If requested to print anything out, go ahead and do so;
506 if (ipr > 0) {
507 vcs_report(iconv);
508 }
509
511
512 // Report on the time if requested to do so
513 double te = tickTock.secondsWC();
514 m_VCount->T_Time_vcs += te;
515 if (ipr > 0 || ip1 > 0) {
517 }
518
519 // FILL IN
520 if (iconv < 0) {
521 plogf("ERROR: FAILURE its = %d!\n", m_VCount->Its);
522 } else if (iconv == 1) {
523 plogf("WARNING: RANGE SPACE ERROR encountered\n");
524 }
525 return iconv;
526}
527
528bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
529{
530 vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop].get();
531 AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhasePossible",
532 "called for a phase that exists!");
533
534 // Loop through all of the species in the phase. We say the phase can be
535 // popped, if there is one species in the phase that can be popped. This
536 // does not mean that the phase will be popped or that it leads to a lower
537 // Gibbs free energy.
538 for (size_t k = 0; k < Vphase->nSpecies(); k++) {
539 size_t kspec = Vphase->spGlobalIndexVCS(k);
541 "VCS_SOLVE::vcs_popPhasePossible",
542 "we shouldn't be here {}: {} > 0.0", kspec,
543 m_molNumSpecies_old[kspec]);
544 size_t irxn = kspec - m_numComponents;
545 if (kspec >= m_numComponents) {
546 bool iPopPossible = true;
547
548 // Note one case is if the component is a member of the popping
549 // phase. This component will be zeroed and the logic here will
550 // negate the current species from causing a positive if this
551 // component is consumed.
552 for (size_t j = 0; j < m_numComponents; ++j) {
553 if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
554 double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
555 if (stoicC != 0.0) {
556 double negChangeComp = - stoicC;
557 if (negChangeComp > 0.0) {
558 // If there is no component to give, then the
559 // species can't be created
561 iPopPossible = false;
562 }
563 }
564 }
565 }
566 }
567 // We are here when the species can be popped because all its needed
568 // components have positive mole numbers
569 if (iPopPossible) {
570 return true;
571 }
572 } else {
573 // We are here when the species, k, in the phase is a component. Its
574 // mole number is zero. We loop through the regular reaction looking
575 // for a reaction that can pop the component.
576 for (size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) {
577 bool foundJrxn = false;
578 // First, if the component is a product of the reaction
579 if (m_stoichCoeffRxnMatrix(kspec,jrxn) > 0.0) {
580 foundJrxn = true;
581 // We can do the reaction if all other reactant components
582 // have positive mole fractions
583 for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
584 if (m_stoichCoeffRxnMatrix(kcomp,jrxn) < 0.0 && m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
585 foundJrxn = false;
586 }
587 }
588 if (foundJrxn) {
589 return true;
590 }
591 } else if (m_stoichCoeffRxnMatrix(kspec,jrxn) < 0.0) {
592 // Second we are here if the component is a reactant in the
593 // reaction, and the reaction goes backwards.
594 foundJrxn = true;
595 size_t jspec = jrxn + m_numComponents;
597 foundJrxn = false;
598 continue;
599 }
600 // We can do the backwards reaction if all of the product
601 // components species are positive
602 for (size_t kcomp = 0; kcomp < m_numComponents; kcomp++) {
603 if (m_stoichCoeffRxnMatrix(kcomp,jrxn) > 0.0 && m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) {
604 foundJrxn = false;
605 }
606 }
607 if (foundJrxn) {
608 return true;
609 }
610 }
611 }
612 }
613 }
614 return false;
615}
616
617size_t VCS_SOLVE::vcs_popPhaseID(vector<size_t> & phasePopPhaseIDs)
618{
619 size_t iphasePop = npos;
620 double FephaseMax = -1.0E30;
621 double Fephase = -1.0E30;
622
623 char anote[128];
624 if (m_debug_print_lvl >= 2) {
625 plogf(" --- vcs_popPhaseID() called\n");
626 plogf(" --- Phase Status F_e MoleNum\n");
627 plogf(" --------------------------------------------------------------------------\n");
628 }
629 for (size_t iph = 0; iph < m_numPhases; iph++) {
630 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
631 int existence = Vphase->exists();
632 strcpy(anote, "");
633 if (existence > 0) {
634 if (m_debug_print_lvl >= 2) {
635 plogf(" --- %18s %5d NA %11.3e\n",
636 Vphase->PhaseName, existence, m_tPhaseMoles_old[iph]);
637 }
638 } else {
639 if (Vphase->m_singleSpecies) {
640 // Single Phase Stability Resolution
641 size_t kspec = Vphase->spGlobalIndexVCS(0);
642 size_t irxn = kspec - m_numComponents;
643 if (irxn > m_deltaGRxn_old.size()) {
644 throw CanteraError("VCS_SOLVE::vcs_popPhaseID",
645 "Index out of bounds due to logic error.");
646 }
647 double deltaGRxn = m_deltaGRxn_old[irxn];
648 Fephase = exp(-deltaGRxn) - 1.0;
649 if (Fephase > 0.0) {
650 strcpy(anote," (ready to be birthed)");
651 if (Fephase > FephaseMax) {
652 iphasePop = iph;
653 FephaseMax = Fephase;
654 strcpy(anote," (chosen to be birthed)");
655 }
656 }
657 if (Fephase < 0.0) {
658 strcpy(anote," (not stable)");
660 "VCS_SOLVE::vcs_popPhaseID", "shouldn't be here");
661 }
662
663 if (m_debug_print_lvl >= 2) {
664 plogf(" --- %18s %5d %10.3g %10.3g %s\n",
665 Vphase->PhaseName, existence, Fephase,
666 m_tPhaseMoles_old[iph], anote);
667 }
668 } else {
669 // MultiSpecies Phase Stability Resolution
670 if (vcs_popPhasePossible(iph)) {
671 Fephase = vcs_phaseStabilityTest(iph);
672 if (Fephase > 0.0) {
673 if (Fephase > FephaseMax) {
674 iphasePop = iph;
675 FephaseMax = Fephase;
676 }
677 } else {
678 FephaseMax = std::max(FephaseMax, Fephase);
679 }
680 if (m_debug_print_lvl >= 2) {
681 plogf(" --- %18s %5d %11.3g %11.3g\n",
682 Vphase->PhaseName, existence, Fephase,
683 m_tPhaseMoles_old[iph]);
684 }
685 } else {
686 if (m_debug_print_lvl >= 2) {
687 plogf(" --- %18s %5d blocked %11.3g\n",
688 Vphase->PhaseName,
689 existence, m_tPhaseMoles_old[iph]);
690 }
691 }
692 }
693 }
694 }
695 phasePopPhaseIDs.resize(0);
696 if (iphasePop != npos) {
697 phasePopPhaseIDs.push_back(iphasePop);
698 }
699
700 // Insert logic here to figure out if phase pops are linked together. Only
701 // do one linked pop at a time.
702 if (m_debug_print_lvl >= 2) {
703 plogf(" ---------------------------------------------------------------------\n");
704 }
705 return iphasePop;
706}
707
708int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
709{
710 vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop].get();
711 // Identify the first species in the phase
712 size_t kspec = Vphase->spGlobalIndexVCS(0);
713 // Identify the formation reaction for that species
714 size_t irxn = kspec - m_numComponents;
715 vector<size_t> creationGlobalRxnNumbers;
716
717 // Calculate the initial moles of the phase being born.
718 // Here we set it to 10x of the value which would cause the phase to be
719 // zeroed out within the algorithm. We may later adjust the value.
720 double tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
721
722 AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhaseRxnStepSizes",
723 "called for a phase that exists!");
724 if (m_debug_print_lvl >= 2) {
725 plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n",
726 Vphase->PhaseName, iphasePop);
727 }
728 // Section for a single-species phase
729 if (Vphase->m_singleSpecies) {
730 double s = 0.0;
731 for (size_t j = 0; j < m_numComponents; ++j) {
732 if (!m_SSPhase[j] && m_molNumSpecies_old[j] > 0.0) {
733 s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j];
734 }
735 }
736 for (size_t j = 0; j < m_numPhases; j++) {
737 Vphase = m_VolPhaseList[j].get();
738 if (! Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) {
739 s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j];
740 }
741 }
742 if (s != 0.0) {
743 double s_old = s;
744 s = vcs_Hessian_diag_adj(irxn, s_old);
745 m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
746 } else {
747 // Ok, s is equal to zero. We can not apply a sophisticated theory
748 // to birth the phase. Just pick a small delta and go with it.
749 m_deltaMolNumSpecies[kspec] = tPhaseMoles;
750 }
751
752 // section to do damping of the m_deltaMolNumSpecies[]
753 for (size_t j = 0; j < m_numComponents; ++j) {
754 double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
755 if (stoicC != 0.0 && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
756 double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
757 if (negChangeComp > m_molNumSpecies_old[j]) {
758 if (m_molNumSpecies_old[j] > 0.0) {
759 m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC;
760 } else {
761 m_deltaMolNumSpecies[kspec] = 0.0;
762 }
763 }
764 }
765 }
766 // Implement a damping term that limits m_deltaMolNumSpecies to the size
767 // of the mole number
768 if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
770 }
771 } else {
772 vector<double> fracDelta(Vphase->nSpecies());
773 vector<double> X_est(Vphase->nSpecies());
774 fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
775
776 double sumFrac = 0.0;
777 for (size_t k = 0; k < Vphase->nSpecies(); k++) {
778 sumFrac += fracDelta[k];
779 }
780 for (size_t k = 0; k < Vphase->nSpecies(); k++) {
781 X_est[k] = fracDelta[k] / sumFrac;
782 }
783
784 double deltaMolNumPhase = tPhaseMoles;
785 double damp = 1.0;
787 double* molNumSpecies_tmp = m_deltaGRxn_tmp.data();
788
789 for (size_t k = 0; k < Vphase->nSpecies(); k++) {
790 kspec = Vphase->spGlobalIndexVCS(k);
791 double delmol = deltaMolNumPhase * X_est[k];
792 if (kspec >= m_numComponents) {
793 irxn = kspec - m_numComponents;
794 if (irxn > m_stoichCoeffRxnMatrix.nColumns()) {
795 throw CanteraError("VCS_SOLVE::vcs_popPhaseRxnStepSizes",
796 "Index out of bounds due to logic error.");
797 }
798 for (size_t j = 0; j < m_numComponents; ++j) {
799 double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
800 if (stoicC != 0.0 && m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
801 molNumSpecies_tmp[j] += stoicC * delmol;
802 }
803 }
804 }
805 }
806
807 double ratioComp = 0.0;
808 for (size_t j = 0; j < m_numComponents; ++j) {
809 double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j];
810 if (molNumSpecies_tmp[j] < 0.0) {
811 ratioComp = 1.0;
812 if (deltaJ > 0.0) {
813 double delta0 = m_molNumSpecies_old[j];
814 damp = std::min(damp, delta0 / deltaJ * 0.9);
815 }
816 } else {
817 if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
818 size_t jph = m_phaseID[j];
819 if ((jph != iphasePop) && (!m_SSPhase[j])) {
820 double fdeltaJ = fabs(deltaJ);
821 if (m_molNumSpecies_old[j] > 0.0) {
822 ratioComp = std::max(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]);
823 }
824 }
825 }
826 }
827 }
828
829 // We may have greatly underestimated the deltaMoles for the phase pop
830 // Here we create a damp > 1 to account for this possibility. We adjust
831 // upwards to make sure that a component in an existing multispecies
832 // phase is modified by a factor of 1/1000.
833 if (ratioComp > 1.0E-30 && ratioComp < 0.001) {
834 damp = 0.001 / ratioComp;
835 }
836 if (damp <= 1.0E-6) {
837 return 3;
838 }
839
840 for (size_t k = 0; k < Vphase->nSpecies(); k++) {
841 kspec = Vphase->spGlobalIndexVCS(k);
842 if (kspec < m_numComponents) {
844 } else {
845 m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp;
846 if (X_est[k] > 1.0E-3) {
848 } else {
850 }
851 }
852 }
853 }
854 return 0;
855}
856
857size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
858{
859 size_t iphDel = npos;
860 size_t k = 0;
861 string ANOTE;
862 if (m_debug_print_lvl >= 2) {
863 plogf(" ");
864 for (int j = 0; j < 82; j++) {
865 plogf("-");
866 }
867 plogf("\n");
868 plogf(" --- Subroutine vcs_RxnStepSizes called - Details:\n");
869 plogf(" ");
870 for (int j = 0; j < 82; j++) {
871 plogf("-");
872 }
873 plogf("\n");
874 plogf(" --- Species KMoles Rxn_Adjustment DeltaG"
875 " | Comment\n");
876 }
877
878 // We update the matrix dlnActCoeffdmolNumber[][] at the top of the loop,
879 // when necessary
880 if (m_useActCoeffJac) {
882 }
883
884 // LOOP OVER THE FORMATION REACTIONS
885 for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
886 ANOTE = "Normal Calc";
887
888 size_t kspec = m_indexRxnToSpecies[irxn];
890 m_deltaMolNumSpecies[kspec] = 0.0;
891 ANOTE = "ZeroedPhase: Phase is artificially zeroed";
893 if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
894 // MULTISPECIES PHASE WITH total moles equal to zero
895 //
896 // If dg[irxn] is negative, then the multispecies phase should
897 // come alive again. Add a small positive step size to make it
898 // come alive.
899 if (m_deltaGRxn_new[irxn] < -1.0e-4) {
900 // First decide if this species is part of a multiphase that
901 // is nontrivial in size.
902 size_t iph = m_phaseID[kspec];
903 double tphmoles = m_tPhaseMoles_old[iph];
904 double trphmoles = tphmoles / m_totalMolNum;
905 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
906 if (Vphase->exists() && (trphmoles > VCS_DELETE_PHASE_CUTOFF)) {
909 m_deltaMolNumSpecies[kspec] = 0.0;
910 ANOTE = fmt::sprintf("MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
912 } else {
914 ANOTE = fmt::sprintf("MultSpec (%s): small species born again DG = %11.3E",
916 }
917 } else {
918 ANOTE = fmt::sprintf("MultSpec (%s):still dead, no phase pop, even though DG = %11.3E",
920 m_deltaMolNumSpecies[kspec] = 0.0;
921 if (Vphase->exists() > 0 && trphmoles > 0.0) {
923 ANOTE = fmt::sprintf("MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E",
925 }
926 }
927 } else {
928 ANOTE = fmt::sprintf("MultSpec (%s): still dead DG = %11.3E",
930 m_deltaMolNumSpecies[kspec] = 0.0;
931 }
932 } else {
933 // REGULAR PROCESSING
934 //
935 // First take care of cases where we want to bail out. Don't
936 // bother if superconvergence has already been achieved in this
937 // mode.
938 if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
939 ANOTE = fmt::sprintf("Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
940 if (m_debug_print_lvl >= 2) {
941 plogf(" --- %-12.12s", m_speciesName[kspec]);
942 plogf(" %12.4E %12.4E %12.4E | %s\n",
944 m_deltaGRxn_new[irxn], ANOTE);
945 }
946 continue;
947 }
948
949 // Don't calculate for minor or nonexistent species if their
950 // values are to be decreasing anyway.
951 if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) {
952 ANOTE = fmt::sprintf("Skipped: IC = %3d and DG >0: %11.3E",
953 m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
954 if (m_debug_print_lvl >= 2) {
955 plogf(" --- %-12.12s", m_speciesName[kspec]);
956 plogf(" %12.4E %12.4E %12.4E | %s\n",
958 m_deltaGRxn_new[irxn], ANOTE);
959 }
960 continue;
961 }
962
963 // Start of the regular processing
964 double s;
965 if (m_SSPhase[kspec]) {
966 s = 0.0;
967 } else {
968 s = 1.0 / m_molNumSpecies_old[kspec];
969 }
970 for (size_t j = 0; j < m_numComponents; ++j) {
971 if (!m_SSPhase[j] && m_molNumSpecies_old[j] > 0.0) {
972 s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j];
973 }
974 }
975 for (size_t j = 0; j < m_numPhases; j++) {
976 vcs_VolPhase* Vphase = m_VolPhaseList[j].get();
977 if (!Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) {
978 s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j];
979 }
980 }
981 if (s != 0.0) {
982 // Take into account of the derivatives of the activity
983 // coefficients with respect to the mole numbers, even in
984 // our diagonal approximation.
985 if (m_useActCoeffJac) {
986 double s_old = s;
987 s = vcs_Hessian_diag_adj(irxn, s_old);
988 ANOTE = fmt::sprintf("Normal calc: diag adjusted from %g "
989 "to %g due to act coeff", s_old, s);
990 }
991
992 m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
993 // New section to do damping of the m_deltaMolNumSpecies[]
994 for (size_t j = 0; j < m_numComponents; ++j) {
995 double stoicC = m_stoichCoeffRxnMatrix(j,irxn);
996 if (stoicC != 0.0) {
997 double negChangeComp = -stoicC * m_deltaMolNumSpecies[kspec];
998 if (negChangeComp > m_molNumSpecies_old[j]) {
999 if (m_molNumSpecies_old[j] > 0.0) {
1000 ANOTE = fmt::sprintf("Delta damped from %g "
1001 "to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec],
1002 -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j]);
1003 m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[j] / stoicC;
1004 } else {
1005 ANOTE = fmt::sprintf("Delta damped from %g "
1006 "to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec],
1007 -m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j]);
1008 m_deltaMolNumSpecies[kspec] = 0.0;
1009 }
1010 }
1011 }
1012 }
1013 // Implement a damping term that limits m_deltaMolNumSpecies
1014 // to the size of the mole number
1015 if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
1016 ANOTE = fmt::sprintf("Delta damped from %g "
1017 "to %g due to %s going negative", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[kspec],
1018 m_speciesName[kspec]);
1020 }
1021 } else {
1022 // REACTION IS ENTIRELY AMONGST SINGLE SPECIES PHASES.
1023 // DELETE ONE OF THE PHASES AND RECOMPUTE BASIS.
1024 //
1025 // Either the species L will disappear or one of the
1026 // component single species phases will disappear. The sign
1027 // of DG(I) will indicate which way the reaction will go.
1028 // Then, we need to follow the reaction to see which species
1029 // will zero out first. The species to be zeroed out will be
1030 // "k".
1031 double dss;
1032 if (m_deltaGRxn_new[irxn] > 0.0) {
1033 dss = m_molNumSpecies_old[kspec];
1034 k = kspec;
1035 for (size_t j = 0; j < m_numComponents; ++j) {
1036 if (m_stoichCoeffRxnMatrix(j,irxn) > 0.0) {
1037 double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix(j,irxn);
1038 if (xx < dss) {
1039 dss = xx;
1040 k = j;
1041 }
1042 }
1043 }
1044 dss = -dss;
1045 } else {
1046 dss = 1.0e10;
1047 for (size_t j = 0; j < m_numComponents; ++j) {
1048 if (m_stoichCoeffRxnMatrix(j,irxn) < 0.0) {
1049 double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix(j,irxn);
1050 if (xx < dss) {
1051 dss = xx;
1052 k = j;
1053 }
1054 }
1055 }
1056 }
1057
1058 // Here we adjust the mole fractions according to DSS and
1059 // the stoichiometric array to take into account that we are
1060 // eliminating the kth species. DSS contains the amount of
1061 // moles of the kth species that needs to be added back into
1062 // the component species.
1063 if (dss != 0.0) {
1064 if ((k == kspec) && (m_SSPhase[kspec] != 1)) {
1065 // Found out that we can be in this spot, when
1066 // components of multispecies phases are zeroed,
1067 // leaving noncomponent species of the same phase
1068 // having all of the mole numbers of that phases. it
1069 // seems that we can suggest a zero of the species
1070 // and the code will recover.
1071 ANOTE = fmt::sprintf("Delta damped from %g to %g due to delete %s", m_deltaMolNumSpecies[kspec],
1072 -m_molNumSpecies_old[kspec], m_speciesName[kspec]);
1074 if (m_debug_print_lvl >= 2) {
1075 plogf(" --- %-12.12s", m_speciesName[kspec]);
1076 plogf(" %12.4E %12.4E %12.4E | %s\n",
1078 m_deltaGRxn_new[irxn], ANOTE);
1079 }
1080 continue;
1081 }
1082
1083 // Delete the single species phase
1084 for (size_t j = 0; j < m_nsp; j++) {
1085 m_deltaMolNumSpecies[j] = 0.0;
1086 }
1087 m_deltaMolNumSpecies[kspec] = dss;
1088 for (size_t j = 0; j < m_numComponents; ++j) {
1090 }
1091
1092 iphDel = m_phaseID[k];
1093 kSpecial = k;
1094
1095 if (k != kspec) {
1096 ANOTE = fmt::sprintf("Delete component SS phase %d named %s - SS phases only",
1097 iphDel, m_speciesName[k]);
1098 } else {
1099 ANOTE = fmt::sprintf("Delete this SS phase %d - SS components only", iphDel);
1100 }
1101 if (m_debug_print_lvl >= 2) {
1102 plogf(" --- %-12.12s", m_speciesName[kspec]);
1103 plogf(" %12.4E %12.4E %12.4E | %s\n",
1105 m_deltaGRxn_new[irxn], ANOTE);
1106 plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s\n",
1107 m_speciesName[k]);
1108 }
1109 if (k != kspec) {
1110 forceComponentCalc = 1;
1111 debuglog(" --- Force a component recalculation\n\n", m_debug_print_lvl >= 2);
1112 }
1113 if (m_debug_print_lvl >= 2) {
1114 plogf(" ");
1115 writeline('-', 82);
1116 }
1117 return iphDel;
1118 }
1119 }
1120 } // End of regular processing
1121 if (m_debug_print_lvl >= 2) {
1122 plogf(" --- %-12.12s", m_speciesName[kspec]);
1123 plogf(" %12.4E %12.4E %12.4E | %s\n",
1125 m_deltaGRxn_new[irxn], ANOTE);
1126 }
1127 } // End of loop over m_speciesUnknownType
1128 } // End of loop over non-component stoichiometric formation reactions
1129 if (m_debug_print_lvl >= 2) {
1130 plogf(" ");
1131 writeline('-', 82);
1132 }
1133 return iphDel;
1134}
1135
1137{
1138 // We will use the _new state calc here
1139 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
1140 const size_t nsp = Vphase->nSpecies();
1141 int minNumberIterations = 3;
1142 if (nsp <= 1) {
1143 minNumberIterations = 1;
1144 }
1145
1146 // We will do a full Newton calculation later, but for now, ...
1147 bool doSuccessiveSubstitution = true;
1148 double funcPhaseStability;
1149 vector<double> X_est(nsp, 0.0);
1150 vector<double> delFrac(nsp, 0.0);
1151 vector<double> E_phi(nsp, 0.0);
1152 vector<double> fracDelta_old(nsp, 0.0);
1153 vector<double> fracDelta_raw(nsp, 0.0);
1154 vector<size_t> creationGlobalRxnNumbers(nsp, npos);
1156 vector<double> feSpecies_Deficient = m_feSpecies_old;
1157
1158 // get the activity coefficients
1160
1161 // Get the stored estimate for the composition of the phase if
1162 // it gets created
1163 vector<double> fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
1164
1165 vector<size_t> componentList;
1166 for (size_t k = 0; k < nsp; k++) {
1167 size_t kspec = Vphase->spGlobalIndexVCS(k);
1168 if (kspec < m_numComponents) {
1169 componentList.push_back(k);
1170 }
1171 }
1172
1173 double normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
1174 double damp = 1.0E-2;
1175
1176 if (doSuccessiveSubstitution) {
1177 int KP = 0;
1178 if (m_debug_print_lvl >= 2) {
1179 plogf(" --- vcs_phaseStabilityTest() called\n");
1180 plogf(" --- Its X_old[%2d] FracDel_old[%2d] deltaF[%2d] FracDel_new[%2d]"
1181 " normUpdate damp FuncPhaseStability\n", KP, KP, KP, KP);
1182 plogf(" --------------------------------------------------------------"
1183 "--------------------------------------------------------\n");
1184 } else if (m_debug_print_lvl == 1) {
1185 plogf(" --- vcs_phaseStabilityTest() called for phase %d\n", iph);
1186 }
1187
1188 for (size_t k = 0; k < nsp; k++) {
1189 if (fracDelta_new[k] < 1.0E-13) {
1190 fracDelta_new[k] =1.0E-13;
1191 }
1192 }
1193 bool converged = false;
1194 double dirProd = 0.0;
1195 for (int its = 0; its < 200 && (!converged); its++) {
1196 double dampOld = damp;
1197 double normUpdateOld = normUpdate;
1198 fracDelta_old = fracDelta_new;
1199 double dirProdOld = dirProd;
1200
1201 // Given a set of fracDelta's, we calculate the fracDelta's
1202 // for the component species, if any
1203 for (size_t i = 0; i < componentList.size(); i++) {
1204 size_t kc = componentList[i];
1205 size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
1206 fracDelta_old[kc] = 0.0;
1207 for (size_t k = 0; k < nsp; k++) {
1208 size_t kspec = Vphase->spGlobalIndexVCS(k);
1209 if (kspec >= m_numComponents) {
1210 size_t irxn = kspec - m_numComponents;
1211 fracDelta_old[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_old[k];
1212 }
1213 }
1214 }
1215
1216 // Now, calculate the predicted mole fractions, X_est[k]
1217 double sumFrac = 0.0;
1218 for (size_t k = 0; k < nsp; k++) {
1219 sumFrac += fracDelta_old[k];
1220 }
1221 // Necessary because this can be identically zero. -> we need to fix
1222 // this algorithm!
1223 if (sumFrac <= 0.0) {
1224 sumFrac = 1.0;
1225 }
1226 double sum_Xcomp = 0.0;
1227 for (size_t k = 0; k < nsp; k++) {
1228 X_est[k] = fracDelta_old[k] / sumFrac;
1229 if (Vphase->spGlobalIndexVCS(k) < m_numComponents) {
1230 sum_Xcomp += X_est[k];
1231 }
1232 }
1233
1234 // Feed the newly formed estimate of the mole fractions back into the
1235 // ThermoPhase object
1236 Vphase->setMoleFractionsState(0.0, &X_est[0], VCS_STATECALC_PHASESTABILITY);
1237
1238 // get the activity coefficients
1240
1241 // First calculate altered chemical potentials for component species
1242 // belonging to this phase.
1243 for (size_t i = 0; i < componentList.size(); i++) {
1244 size_t kc = componentList[i];
1245 size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
1246 if (X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
1247 feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
1248 + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
1249 } else {
1250 feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
1252 }
1253 }
1254
1255 for (size_t i = 0; i < componentList.size(); i++) {
1256 size_t kc_spec = Vphase->spGlobalIndexVCS(componentList[i]);
1257 for (size_t k = 0; k < Vphase->nSpecies(); k++) {
1258 size_t kspec = Vphase->spGlobalIndexVCS(k);
1259 if (kspec >= m_numComponents) {
1260 size_t irxn = kspec - m_numComponents;
1261 if (i == 0) {
1263 }
1264 if (m_stoichCoeffRxnMatrix(kc_spec,irxn) != 0.0) {
1265 m_deltaGRxn_Deficient[irxn] +=
1266 m_stoichCoeffRxnMatrix(kc_spec,irxn) * (feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
1267 }
1268 }
1269 }
1270 }
1271
1272 // Calculate the E_phi's
1273 double sum = 0.0;
1274 funcPhaseStability = sum_Xcomp - 1.0;
1275 for (size_t k = 0; k < nsp; k++) {
1276 size_t kspec = Vphase->spGlobalIndexVCS(k);
1277 if (kspec >= m_numComponents) {
1278 size_t irxn = kspec - m_numComponents;
1279 double deltaGRxn = clip(m_deltaGRxn_Deficient[irxn], -50.0, 50.0);
1280 E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec];
1281 sum += E_phi[k];
1282 funcPhaseStability += E_phi[k];
1283 } else {
1284 E_phi[k] = 0.0;
1285 }
1286 }
1287
1288 // Calculate the raw estimate of the new fracs
1289 for (size_t k = 0; k < nsp; k++) {
1290 size_t kspec = Vphase->spGlobalIndexVCS(k);
1291 double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
1292 if (kspec >= m_numComponents) {
1293 fracDelta_raw[k] = b;
1294 }
1295 }
1296
1297 // Given a set of fracDelta's, we calculate the fracDelta's
1298 // for the component species, if any
1299 for (size_t i = 0; i < componentList.size(); i++) {
1300 size_t kc = componentList[i];
1301 size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
1302 fracDelta_raw[kc] = 0.0;
1303 for (size_t k = 0; k < nsp; k++) {
1304 size_t kspec = Vphase->spGlobalIndexVCS(k);
1305 if (kspec >= m_numComponents) {
1306 size_t irxn = kspec - m_numComponents;
1307 fracDelta_raw[kc] += m_stoichCoeffRxnMatrix(kc_spec,irxn) * fracDelta_raw[k];
1308 }
1309 }
1310 }
1311
1312 // Now possibly dampen the estimate.
1313 for (size_t k = 0; k < nsp; k++) {
1314 delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
1315 }
1316 normUpdate = vcs_l2norm(delFrac);
1317
1318 dirProd = 0.0;
1319 for (size_t k = 0; k < nsp; k++) {
1320 dirProd += fracDelta_old[k] * delFrac[k];
1321 }
1322 bool crossedSign = false;
1323 if (dirProd * dirProdOld < 0.0) {
1324 crossedSign = true;
1325 }
1326
1327 damp = 0.5;
1328 if (dampOld < 0.25) {
1329 damp = 2.0 * dampOld;
1330 }
1331 if (crossedSign) {
1332 if (normUpdate *1.5 > normUpdateOld) {
1333 damp = 0.5 * dampOld;
1334 } else if (normUpdate *2.0 > normUpdateOld) {
1335 damp = 0.8 * dampOld;
1336 }
1337 } else {
1338 if (normUpdate > normUpdateOld * 2.0) {
1339 damp = 0.6 * dampOld;
1340 } else if (normUpdate > normUpdateOld * 1.2) {
1341 damp = 0.9 * dampOld;
1342 }
1343 }
1344
1345 for (size_t k = 0; k < nsp; k++) {
1346 if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
1347 damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
1348 1.0E-8/fabs(delFrac[k]));
1349 }
1350 if (delFrac[k] < 0.0 && 2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) {
1351 damp = fracDelta_old[k] / (2.0 * -delFrac[k]);
1352 }
1353 if (delFrac[k] > 0.0 && 2.0 * damp * delFrac[k] > fracDelta_old[k]) {
1354 damp = fracDelta_old[k] / (2.0 * delFrac[k]);
1355 }
1356 }
1357 damp = std::max(damp, 0.000001);
1358 for (size_t k = 0; k < nsp; k++) {
1359 fracDelta_new[k] = fracDelta_old[k] + damp * delFrac[k];
1360 }
1361
1362 if (m_debug_print_lvl >= 2) {
1363 plogf(" --- %3d %12g %12g %12g %12g %12g %12g %12g\n", its, X_est[KP], fracDelta_old[KP],
1364 delFrac[KP], fracDelta_new[KP], normUpdate, damp, funcPhaseStability);
1365 }
1366
1367 if (normUpdate < 1.0E-5 * damp) {
1368 converged = true;
1369 if (its < minNumberIterations) {
1370 converged = false;
1371 }
1372 }
1373 }
1374
1375 if (converged) {
1376 // Save the final optimized stated back into the VolPhase object for later use
1377 Vphase->setMoleFractionsState(0.0, &X_est[0], VCS_STATECALC_PHASESTABILITY);
1378
1379 // Save fracDelta for later use to initialize the problem better
1380 // @todo creationGlobalRxnNumbers needs to be calculated here and stored.
1381 Vphase->setCreationMoleNumbers(&fracDelta_new[0], creationGlobalRxnNumbers);
1382 }
1383 } else {
1384 throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet");
1385 }
1386 if (m_debug_print_lvl >= 2) {
1387 plogf(" ------------------------------------------------------------"
1388 "-------------------------------------------------------------\n");
1389 } else if (m_debug_print_lvl == 1) {
1390 if (funcPhaseStability > 0.0) {
1391 plogf(" --- phase %d with func = %g is to be born\n", iph, funcPhaseStability);
1392 } else {
1393 plogf(" --- phase %d with func = %g stays dead\n", iph, funcPhaseStability);
1394 }
1395 }
1396 return funcPhaseStability;
1397}
1398
1399int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg)
1400{
1401 // Store the temperature and pressure in the private global variables
1402 m_temperature = T_arg;
1403 m_pressurePA = pres_arg;
1405
1406 // Evaluate the standard state free energies
1407 // at the current temperatures and pressures.
1408 int iconv = vcs_evalSS_TP(ipr, ip1, m_temperature, pres_arg);
1409
1410 // Prep the fe field
1411 vcs_fePrep_TP();
1412
1413 // Decide whether we need an initial estimate of the solution If so, go get
1414 // one. If not, then
1415 if (m_doEstimateEquil) {
1416 int retn = vcs_inest_TP();
1417 if (retn != VCS_SUCCESS) {
1418 plogf("vcs_inest_TP returned a failure flag\n");
1419 }
1420 }
1421
1422 // Solve the problem at a fixed Temperature and Pressure (all information
1423 // concerning Temperature and Pressure has already been derived. The free
1424 // energies are now in dimensionless form.)
1425 iconv = vcs_solve_TP(ipr, ip1, maxit);
1426
1427 // Return the convergence success flag.
1428 return iconv;
1429}
1430
1431int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
1432{
1433 for (size_t iph = 0; iph < m_numPhases; iph++) {
1434 vcs_VolPhase* vph = m_VolPhaseList[iph].get();
1437 }
1438 for (size_t k = 0; k < m_nsp; k++) {
1440 }
1441
1442 return VCS_SUCCESS;
1443}
1444
1446{
1447 for (size_t i = 0; i < m_nsp; ++i) {
1448 // For single species phases, initialize the chemical potential with the
1449 // value of the standard state chemical potential. This value doesn't
1450 // change during the calculation
1451 if (m_SSPhase[i]) {
1454 }
1455 }
1456}
1457
1458double VCS_SOLVE::vcs_VolTotal(const double tkelvin, const double pres,
1459 const double w[], double volPM[])
1460{
1461 double VolTot = 0.0;
1462 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
1463 vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
1464 Vphase->setState_TP(tkelvin, pres);
1466 double Volp = Vphase->sendToVCS_VolPM(volPM);
1467 VolTot += Volp;
1468 }
1469 return VolTot;
1470}
1471
1472int VCS_SOLVE::vcs_prep(int printLvl)
1473{
1474 int retn = VCS_SUCCESS;
1475 m_debug_print_lvl = printLvl;
1476
1477 // Calculate the Single Species status of phases
1478 // Also calculate the number of species per phase
1479 vcs_SSPhase();
1480
1481 // Set an initial estimate for the number of noncomponent species equal to
1482 // nspecies - nelements. This may be changed below
1483 if (m_nelem > m_nsp) {
1484 m_numRxnTot = 0;
1485 } else {
1487 }
1490 for (size_t i = 0; i < m_numRxnRdc; ++i) {
1492 }
1493
1494 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
1495 size_t pID = m_phaseID[kspec];
1496 size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
1497 vcs_VolPhase* vPhase = m_VolPhaseList[pID].get();
1498 vcs_SpeciesProperties* spProp = vPhase->speciesProperty(spPhIndex);
1499 double sz = 0.0;
1500 size_t eSize = spProp->FormulaMatrixCol.size();
1501 for (size_t e = 0; e < eSize; e++) {
1502 sz += fabs(spProp->FormulaMatrixCol[e]);
1503 }
1504 if (sz > 0.0) {
1505 m_spSize[kspec] = sz;
1506 } else {
1507 m_spSize[kspec] = 1.0;
1508 }
1509 }
1510
1511 // DETERMINE THE NUMBER OF COMPONENTS
1512 //
1513 // Obtain a valid estimate of the mole fraction. This will be used as an
1514 // initial ordering vector for prioritizing which species are defined as
1515 // components.
1516 //
1517 // If a mole number estimate was supplied from the input file, use that mole
1518 // number estimate.
1519 //
1520 // If a solution estimate wasn't supplied from the input file, supply an
1521 // initial estimate for the mole fractions based on the relative reverse
1522 // ordering of the chemical potentials.
1523 //
1524 // For voltage unknowns, set these to zero for the moment.
1525 double test = -1.0e-10;
1526 bool modifiedSoln = false;
1527 if (m_doEstimateEquil < 0) {
1528 double sum = 0.0;
1529 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
1531 sum += fabs(m_molNumSpecies_old[kspec]);
1532 }
1533 }
1534 if (fabs(sum) < 1.0E-6) {
1535 modifiedSoln = true;
1536 double pres = (m_pressurePA <= 0.0) ? 1.01325E5 : m_pressurePA;
1537 retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
1538 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
1540 m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
1541 } else {
1542 m_molNumSpecies_old[kspec] = 0.0;
1543 }
1544 }
1545 }
1546 test = -1.0e20;
1547 }
1548
1549 // NC = number of components is in the vcs.h common block. This call to
1550 // BASOPT doesn't calculate the stoichiometric reaction matrix.
1551 vector<double> awSpace(m_nsp + (m_nelem + 2)*(m_nelem), 0.0);
1552 double* aw = &awSpace[0];
1553 if (aw == NULL) {
1554 plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
1555 return VCS_NOMEMORY;
1556 }
1557 double* sa = aw + m_nsp;
1558 double* sm = sa + m_nelem;
1559 double* ss = sm + m_nelem * m_nelem;
1560 bool conv;
1561 retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv);
1562 if (retn != VCS_SUCCESS) {
1563 plogf("vcs_prep_oneTime:");
1564 plogf(" Determination of number of components failed: %d\n",
1565 retn);
1566 plogf(" Global Bailout!\n");
1567 return retn;
1568 }
1569
1570 if (m_nsp >= m_numComponents) {
1572 for (size_t i = 0; i < m_numRxnRdc; ++i) {
1574 }
1575 } else {
1577 }
1578
1579 // The elements might need to be rearranged.
1580 awSpace.resize(m_nelem + (m_nelem + 2)*m_nelem, 0.0);
1581 aw = &awSpace[0];
1582 sa = aw + m_nelem;
1583 sm = sa + m_nelem;
1584 ss = sm + m_nelem * m_nelem;
1585 retn = vcs_elem_rearrange(aw, sa, sm, ss);
1586 if (retn != VCS_SUCCESS) {
1587 plogf("vcs_prep_oneTime:");
1588 plogf(" Determination of element reordering failed: %d\n",
1589 retn);
1590 plogf(" Global Bailout!\n");
1591 return retn;
1592 }
1593
1594 // If we mucked up the solution unknowns because they were all
1595 // zero to start with, set them back to zero here
1596 if (modifiedSoln) {
1597 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
1598 m_molNumSpecies_old[kspec] = 0.0;
1599 }
1600 }
1601
1602 // Initialize various arrays in the data to zero
1603 m_feSpecies_old.assign(m_feSpecies_old.size(), 0.0);
1604 m_feSpecies_new.assign(m_feSpecies_new.size(), 0.0);
1605 m_molNumSpecies_new.assign(m_molNumSpecies_new.size(), 0.0);
1608 m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
1609 m_tPhaseMoles_new.assign(m_tPhaseMoles_new.size(), 0.0);
1610
1611 // Calculate the total number of moles in all phases.
1612 vcs_tmoles();
1613
1614 // Check to see if the current problem is well posed.
1615 double sum = 0.0;
1616 for (size_t e = 0; e < m_mix->nElements(); e++) {
1617 sum += m_mix->elementMoles(e);
1618 }
1619 if (sum < 1.0E-20) {
1620 // Check to see if the current problem is well posed.
1621 plogf("vcs has determined the problem is not well posed: Bailing\n");
1622 return VCS_PUB_BAD;
1623 }
1624 return VCS_SUCCESS;
1625}
1626
1627int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
1628 double* const sm, double* const ss)
1629{
1630 size_t ncomponents = m_numComponents;
1631 if (m_debug_print_lvl >= 2) {
1632 plogf(" ");
1633 writeline('-', 77);
1634 plogf(" --- Subroutine elem_rearrange() called to ");
1635 plogf("check stoich. coefficient matrix\n");
1636 plogf(" --- and to rearrange the element ordering once\n");
1637 }
1638
1639 // Use a temporary work array for the element numbers
1640 // Also make sure the value of test is unique.
1641 bool lindep = true;
1642 double test = -1.0E10;
1643 while (lindep) {
1644 lindep = false;
1645 for (size_t i = 0; i < m_nelem; ++i) {
1646 test -= 1.0;
1647 aw[i] = m_elemAbundancesGoal[i];
1648 if (test == aw[i]) {
1649 lindep = true;
1650 }
1651 }
1652 }
1653
1654 // Top of a loop of some sort based on the index JR. JR is the current
1655 // number independent elements found.
1656 size_t jr = 0;
1657 while (jr < ncomponents) {
1658 size_t k;
1659
1660 // Top of another loop point based on finding a linearly independent
1661 // species
1662 while (true) {
1663 // Search the remaining part of the mole fraction vector, AW, for
1664 // the largest remaining species. Return its identity in K.
1665 k = m_nelem;
1666 for (size_t ielem = jr; ielem < m_nelem; ielem++) {
1667 if (m_elementActive[ielem] && aw[ielem] != test) {
1668 k = ielem;
1669 break;
1670 }
1671 }
1672 if (k == m_nelem) {
1673 throw CanteraError("VCS_SOLVE::vcs_elem_rearrange",
1674 "Shouldn't be here. Algorithm misfired.");
1675 }
1676
1677 // Assign a large negative number to the element that we have just
1678 // found, in order to take it out of further consideration.
1679 aw[k] = test;
1680
1681 // CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX LINE WITH
1682 // PREVIOUS LINES OF THE FORMULA MATRIX
1683 //
1684 // Modified Gram-Schmidt Method, p. 202 Dalquist QR factorization of
1685 // a matrix without row pivoting.
1686 size_t jl = jr;
1687
1688 // Fill in the row for the current element, k, under consideration
1689 // The row will contain the Formula matrix value for that element
1690 // from the current component.
1691 for (size_t j = 0; j < ncomponents; ++j) {
1692 sm[j + jr*ncomponents] = m_formulaMatrix(j,k);
1693 }
1694 if (jl > 0) {
1695 // Compute the coefficients of JA column of the the upper
1696 // triangular R matrix, SS(J) = R_J_JR (this is slightly
1697 // different than Dalquist) R_JA_JA = 1
1698 for (size_t j = 0; j < jl; ++j) {
1699 ss[j] = 0.0;
1700 for (size_t i = 0; i < ncomponents; ++i) {
1701 ss[j] += sm[i + jr*ncomponents] * sm[i + j*ncomponents];
1702 }
1703 ss[j] /= sa[j];
1704 }
1705
1706 // Now make the new column, (*,JR), orthogonal to the previous
1707 // columns
1708 for (size_t j = 0; j < jl; ++j) {
1709 for (size_t i = 0; i < ncomponents; ++i) {
1710 sm[i + jr*ncomponents] -= ss[j] * sm[i + j*ncomponents];
1711 }
1712 }
1713 }
1714
1715 // Find the new length of the new column in Q. It will be used in
1716 // the denominator in future row calcs.
1717 sa[jr] = 0.0;
1718 for (size_t ml = 0; ml < ncomponents; ++ml) {
1719 sa[jr] += pow(sm[ml + jr*ncomponents], 2);
1720 }
1721 // IF NORM OF NEW ROW .LT. 1E-6 REJECT
1722 if (sa[jr] > 1.0e-6) {
1723 break;
1724 }
1725 }
1726 // REARRANGE THE DATA
1727 if (jr != k) {
1728 if (m_debug_print_lvl >= 2) {
1729 plogf(" --- %-2.2s(%9.2g) replaces %-2.2s(%9.2g) as element %3d\n",
1731 m_elementName[jr], m_elemAbundancesGoal[jr], jr);
1732 }
1733 vcs_switch_elem_pos(jr, k);
1734 std::swap(aw[jr], aw[k]);
1735 }
1736
1737 // If we haven't found enough components, go back and find some more.
1738 jr++;
1739 }
1740 return VCS_SUCCESS;
1741}
1742
1743void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
1744{
1745 if (ipos == jpos) {
1746 return;
1747 }
1748 AssertThrowMsg(ipos < m_nelem && jpos < m_nelem,
1749 "vcs_switch_elem_pos",
1750 "inappropriate args: {} {}", ipos, jpos);
1751
1752 // Change the element Global Index list in each vcs_VolPhase object
1753 // to reflect the switch in the element positions.
1754 for (size_t iph = 0; iph < m_numPhases; iph++) {
1755 vcs_VolPhase* volPhase = m_VolPhaseList[iph].get();
1756 for (size_t e = 0; e < volPhase->nElemConstraints(); e++) {
1757 if (volPhase->elemGlobalIndex(e) == ipos) {
1758 volPhase->setElemGlobalIndex(e, jpos);
1759 }
1760 if (volPhase->elemGlobalIndex(e) == jpos) {
1761 volPhase->setElemGlobalIndex(e, ipos);
1762 }
1763 }
1764 }
1765 std::swap(m_elemAbundancesGoal[ipos], m_elemAbundancesGoal[jpos]);
1766 std::swap(m_elemAbundances[ipos], m_elemAbundances[jpos]);
1767 std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]);
1768 std::swap(m_elType[ipos], m_elType[jpos]);
1769 std::swap(m_elementActive[ipos], m_elementActive[jpos]);
1770 for (size_t j = 0; j < m_nsp; ++j) {
1771 std::swap(m_formulaMatrix(j,ipos), m_formulaMatrix(j,jpos));
1772 }
1773 std::swap(m_elementName[ipos], m_elementName[jpos]);
1774}
1775
1776double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
1777{
1778 double diag = hessianDiag_Ideal;
1779 double hessActCoef = vcs_Hessian_actCoeff_diag(irxn);
1780 if (hessianDiag_Ideal <= 0.0) {
1781 throw CanteraError("VCS_SOLVE::vcs_Hessian_diag_adj",
1782 "We shouldn't be here");
1783 }
1784 if (hessActCoef >= 0.0) {
1785 diag += hessActCoef;
1786 } else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
1787 diag += hessActCoef;
1788 } else {
1789 diag -= 0.6666 * hessianDiag_Ideal;
1790 }
1791 return diag;
1792}
1793
1795{
1796 size_t kspec = m_indexRxnToSpecies[irxn];
1797 size_t kph = m_phaseID[kspec];
1798 double np_kspec = std::max(m_tPhaseMoles_old[kph], 1e-13);
1799 double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
1800
1801 // First the diagonal term of the Jacobian
1802 double s = m_np_dLnActCoeffdMolNum(kspec,kspec) / np_kspec;
1803
1804 // Next, the other terms. Note this only a loop over the components So, it's
1805 // not too expensive to calculate.
1806 for (size_t j = 0; j < m_numComponents; j++) {
1807 if (!m_SSPhase[j]) {
1808 for (size_t k = 0; k < m_numComponents; ++k) {
1809 if (m_phaseID[k] == m_phaseID[j]) {
1810 double np = m_tPhaseMoles_old[m_phaseID[k]];
1811 if (np > 0.0) {
1812 s += sc_irxn[k] * sc_irxn[j] * m_np_dLnActCoeffdMolNum(j,k) / np;
1813 }
1814 }
1815 }
1816 if (kph == m_phaseID[j]) {
1817 s += sc_irxn[j] * (m_np_dLnActCoeffdMolNum(j,kspec) + m_np_dLnActCoeffdMolNum(kspec,j)) / np_kspec;
1818 }
1819 }
1820 }
1821 return s;
1822}
1823
1824void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS)
1825{
1826 // Loop over all of the phases in the problem
1827 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
1828 vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
1829
1830 // We don't need to call single species phases;
1831 if (!Vphase->m_singleSpecies && !Vphase->isIdealSoln()) {
1832 // update the mole numbers
1833 Vphase->setMolesFromVCS(VCS_STATECALC_OLD, moleSpeciesVCS);
1834
1835 // Download the resulting calculation into the full vector. This
1836 // scatter calculation is carried out in the vcs_VolPhase object.
1838 }
1839 }
1840}
1841
1843{
1844 bool inertYes = false;
1845
1846 // SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES
1847 vector<pair<double, size_t>> x_order;
1848 for (size_t i = 0; i < m_nsp; i++) {
1849 x_order.push_back({-m_molNumSpecies_old[i], i});
1850 }
1851 std::sort(x_order.begin() + m_numComponents,
1852 x_order.begin() + m_numSpeciesRdc);
1853
1854 vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
1856
1857 // PRINT OUT RESULTS
1858 plogf("\n\n\n\n");
1859 writeline('-', 80);
1860 writeline('-', 80);
1861 plogf("\t\t VCS_TP REPORT\n");
1862 writeline('-', 80);
1863 writeline('-', 80);
1864 if (iconv < 0) {
1865 plogf(" ERROR: CONVERGENCE CRITERION NOT SATISFIED.\n");
1866 } else if (iconv == 1) {
1867 plogf(" RANGE SPACE ERROR: Equilibrium Found but not all Element Abundances are Satisfied\n");
1868 }
1869
1870 // Calculate some quantities that may need updating
1871 vcs_tmoles();
1874
1875 plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
1876 plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
1877 plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol);
1878
1879 // TABLE OF SPECIES IN DECREASING MOLE NUMBERS
1880 plogf("\n\n");
1881 writeline('-', 80);
1882 plogf(" Species Equilibrium kmoles ");
1883 plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
1884 writeline('-', 80);
1885 for (size_t i = 0; i < m_numComponents; ++i) {
1886 plogf(" %-12.12s", m_speciesName[i]);
1887 writeline(' ', 13, false);
1888 plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i],
1890 plogf(" %3d", m_speciesUnknownType[i]);
1891 plogf("\n");
1892 }
1893 for (size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
1894 size_t j = x_order[i].second;
1895 plogf(" %-12.12s", m_speciesName[j]);
1896 writeline(' ', 13, false);
1897
1899 plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j],
1901 plogf(" KMolNum ");
1903 plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[j]);
1904 plogf(" Voltage = %14.7E", m_molNumSpecies_old[j]);
1905 } else {
1906 throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
1907 }
1908 plogf("\n");
1909 }
1910 for (size_t i = 0; i < m_numPhases; i++) {
1911 if (TPhInertMoles[i] > 0.0) {
1912 inertYes = true;
1913 if (i == 0) {
1914 plogf(" Inert Gas Species ");
1915 } else {
1916 plogf(" Inert Species in phase %16s ",
1917 m_VolPhaseList[i]->PhaseName);
1918 }
1919 plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i],
1920 TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
1921 }
1922 }
1923 if (m_numSpeciesRdc != m_nsp) {
1924 plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
1925 for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
1926 plogf(" %-12.12s", m_speciesName[kspec]);
1927 // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
1928 plogf(" %14.7E %14.7E %12.4E",
1929 m_molNumSpecies_old[kspec],
1930 m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]);
1932 plogf(" KMol_Num");
1934 plogf(" Voltage");
1935 } else {
1936 plogf(" Unknown");
1937 }
1938 plogf("\n");
1939 }
1940 }
1941 writeline('-', 80);
1942 plogf("\n");
1943
1944 // TABLE OF SPECIES FORMATION REACTIONS
1945 writeline('-', m_numComponents*10 + 45, true, true);
1946 plogf(" |ComponentID|");
1947 for (size_t j = 0; j < m_numComponents; j++) {
1948 plogf(" %3d", j);
1949 }
1950 plogf(" | |\n");
1951 plogf(" | Components|");
1952 for (size_t j = 0; j < m_numComponents; j++) {
1953 plogf(" %10.10s", m_speciesName[j]);
1954 }
1955 plogf(" | |\n");
1956 plogf(" NonComponent | Moles |");
1957 for (size_t j = 0; j < m_numComponents; j++) {
1958 plogf(" %10.3g", m_molNumSpecies_old[j]);
1959 }
1960 plogf(" | DG/RT Rxn |\n");
1961 writeline('-', m_numComponents*10 + 45);
1962 for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
1963 size_t kspec = m_indexRxnToSpecies[irxn];
1964 plogf(" %3d ", kspec);
1965 plogf("%-10.10s", m_speciesName[kspec]);
1966 plogf("|%10.3g |", m_molNumSpecies_old[kspec]);
1967 for (size_t j = 0; j < m_numComponents; j++) {
1968 plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,irxn));
1969 }
1970 plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
1971 plogf("\n");
1972 }
1973 writeline('-', m_numComponents*10 + 45);
1974 plogf("\n");
1975
1976 // TABLE OF PHASE INFORMATION
1977 vector<double> gaPhase(m_nelem, 0.0);
1978 vector<double> gaTPhase(m_nelem, 0.0);
1979 double totalMoles = 0.0;
1980 double gibbsPhase = 0.0;
1981 double gibbsTotal = 0.0;
1982 plogf("\n\n");
1983 plogf("\n");
1984 writeline('-', m_nelem*10 + 58);
1985 plogf(" | ElementID |");
1986 for (size_t j = 0; j < m_nelem; j++) {
1987 plogf(" %3d", j);
1988 }
1989 plogf(" | |\n");
1990 plogf(" | Element |");
1991 for (size_t j = 0; j < m_nelem; j++) {
1992 plogf(" %10.10s", m_elementName[j]);
1993 }
1994 plogf(" | |\n");
1995 plogf(" PhaseName |KMolTarget |");
1996 for (size_t j = 0; j < m_nelem; j++) {
1997 plogf(" %10.3g", m_elemAbundancesGoal[j]);
1998 }
1999 plogf(" | Gibbs Total |\n");
2000 writeline('-', m_nelem*10 + 58);
2001 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
2002 plogf(" %3d ", iphase);
2003 vcs_VolPhase* VPhase = m_VolPhaseList[iphase].get();
2004 plogf("%-12.12s |",VPhase->PhaseName);
2005 plogf("%10.3e |", m_tPhaseMoles_old[iphase]);
2006 totalMoles += m_tPhaseMoles_old[iphase];
2007 if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles() &&
2008 !vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) {
2009 throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
2010 }
2011 vcs_elabPhase(iphase, &gaPhase[0]);
2012 for (size_t j = 0; j < m_nelem; j++) {
2013 plogf(" %10.3g", gaPhase[j]);
2014 gaTPhase[j] += gaPhase[j];
2015 }
2016 gibbsPhase = vcs_GibbsPhase(iphase, &m_molNumSpecies_old[0],
2017 &m_feSpecies_old[0]);
2018 gibbsTotal += gibbsPhase;
2019 plogf(" | %18.11E |\n", gibbsPhase);
2020 }
2021 writeline('-', m_nelem*10 + 58);
2022 plogf(" TOTAL |%10.3e |", totalMoles);
2023 for (size_t j = 0; j < m_nelem; j++) {
2024 plogf(" %10.3g", gaTPhase[j]);
2025 }
2026 plogf(" | %18.11E |\n", gibbsTotal);
2027
2028 writeline('-', m_nelem*10 + 58);
2029 plogf("\n");
2030
2031 // GLOBAL SATISFACTION INFORMATION
2032
2033 // Calculate the total dimensionless Gibbs Free Energy. Inert species are
2034 // handled as if they had a standard free energy of zero
2036 &m_tPhaseMoles_old[0]);
2037 plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
2038 if (inertYes) {
2039 plogf("\t\t(Inert species have standard free energy of zero)\n");
2040 }
2041
2042 plogf("\nElemental Abundances (kmol): ");
2043 plogf(" Actual Target Type ElActive\n");
2044 for (size_t i = 0; i < m_nelem; ++i) {
2045 writeline(' ', 26, false);
2046 plogf("%-2.2s", m_elementName[i]);
2047 plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]);
2048 plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]);
2049 }
2050 plogf("\n");
2051
2052 // TABLE OF SPECIES CHEM POTS
2053 writeline('-', 93, true, true);
2054 plogf("Chemical Potentials of the Species: (dimensionless)\n");
2055
2056 plogf("\t\t(RT = %g J/kmol)\n", GasConstant * m_temperature);
2057 plogf(" Name TKMoles StandStateChemPot "
2058 " ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
2059 plogf("| (MolNum ChemPot)|");
2060 writeline('-', 147, true, true);
2061 for (size_t i = 0; i < m_nsp; ++i) {
2062 size_t j = x_order[i].second;
2063 size_t pid = m_phaseID[j];
2064 plogf(" %-12.12s", m_speciesName[j]);
2065 plogf(" %14.7E ", m_molNumSpecies_old[j]);
2066 plogf("%14.7E ", m_SSfeSpecies[j]);
2067 plogf("%14.7E ", log(m_actCoeffSpecies_old[j]));
2068 double tpmoles = m_tPhaseMoles_old[pid];
2069 double phi = m_phasePhi[pid];
2070 double eContrib = phi * m_chargeSpecies[j] * m_Faraday_dim;
2071 double lx = 0.0;
2073 lx = 0.0;
2074 } else {
2075 if (tpmoles > 0.0 && m_molNumSpecies_old[j] > 0.0) {
2076 double tmp = std::max(VCS_DELETE_MINORSPECIES_CUTOFF, m_molNumSpecies_old[j]);
2077 lx = log(tmp) - log(tpmoles);
2078 } else {
2079 lx = m_feSpecies_old[j] - m_SSfeSpecies[j]
2081 }
2082 }
2083 plogf("%14.7E |", lx);
2084 plogf("%14.7E | ", eContrib);
2085 double tmp = m_SSfeSpecies[j] + log(m_actCoeffSpecies_old[j])
2086 + lx - m_lnMnaughtSpecies[j] + eContrib;
2087 if (fabs(m_feSpecies_old[j] - tmp) > 1.0E-7) {
2088 throw CanteraError("VCS_SOLVE::vcs_report",
2089 "we have a problem - doesn't add up");
2090 }
2091 plogf(" %12.4E |", m_feSpecies_old[j]);
2092 if (m_lnMnaughtSpecies[j] != 0.0) {
2093 plogf("(%11.5E)", - m_lnMnaughtSpecies[j]);
2094 } else {
2095 plogf(" ");
2096 }
2097
2098 plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j]);
2099 plogf("\n");
2100 }
2101 for (size_t i = 0; i < 125; i++) {
2102 plogf(" ");
2103 }
2104 plogf(" %20.9E\n", g);
2105 writeline('-', 147);
2106
2107 // TABLE OF SOLUTION COUNTERS
2108 plogf("\n");
2109 plogf("\nCounters: Iterations Time (seconds)\n");
2110 if (m_timing_print_lvl > 0) {
2111 plogf(" vcs_basopt: %5d %11.5E\n",
2113 plogf(" vcs_TP: %5d %11.5E\n",
2115 } else {
2116 plogf(" vcs_basopt: %5d %11s\n",
2117 m_VCount->Basis_Opts," NA ");
2118 plogf(" vcs_TP: %5d %11s\n",
2119 m_VCount->Its," NA ");
2120 }
2121 writeline('-', 80);
2122 writeline('-', 80);
2123
2124 // Return a successful completion flag
2125 return VCS_SUCCESS;
2126}
2127
2129{
2130 for (size_t j = 0; j < m_nelem; ++j) {
2131 m_elemAbundances[j] = 0.0;
2132 for (size_t i = 0; i < m_nsp; ++i) {
2135 }
2136 }
2137 }
2138}
2139
2141{
2142 size_t top = m_numComponents;
2143 if (ibound) {
2144 top = m_nelem;
2145 }
2146
2147 for (size_t i = 0; i < top; ++i) {
2148 // Require 12 digits of accuracy on non-zero constraints.
2149 if (m_elementActive[i] && fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > fabs(m_elemAbundancesGoal[i]) * 1.0e-12) {
2150 // This logic is for charge neutrality condition
2152 m_elemAbundancesGoal[i] != 0.0) {
2153 throw CanteraError("VCS_SOLVE::vcs_elabcheck",
2154 "Problem with charge neutrality condition");
2155 }
2158
2159 // Find out if the constraint is a multisign constraint. If it
2160 // is, then we have to worry about roundoff error in the
2161 // addition of terms. We are limited to 13 digits of finite
2162 // arithmetic accuracy.
2163 bool multisign = false;
2164 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
2165 double eval = m_formulaMatrix(kspec,i);
2166 if (eval < 0.0) {
2167 multisign = true;
2168 }
2169 if (eval != 0.0) {
2170 scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec]));
2171 }
2172 }
2173 if (multisign) {
2174 if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) {
2175 return false;
2176 }
2177 } else {
2179 return false;
2180 }
2181 }
2182 } else {
2183 // For normal element balances, we require absolute compliance
2184 // even for ridiculously small numbers.
2185 if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
2186 return false;
2187 } else {
2188 return false;
2189 }
2190 }
2191 }
2192 }
2193 return true;
2194}
2195
2196void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase)
2197{
2198 for (size_t j = 0; j < m_nelem; ++j) {
2199 elemAbundPhase[j] = 0.0;
2200 for (size_t i = 0; i < m_nsp; ++i) {
2202 elemAbundPhase[j] += m_formulaMatrix(i,j) * m_molNumSpecies_old[i];
2203 }
2204 }
2205 }
2206}
2207
2208int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
2209{
2210 int retn = 0;
2211
2212 vector<double> ga_save(m_elemAbundances);
2213 if (m_debug_print_lvl >= 2) {
2214 plogf(" --- vcsc_elcorr: Element abundances correction routine");
2215 if (m_nelem != m_numComponents) {
2216 plogf(" (m_numComponents != m_nelem)");
2217 }
2218 plogf("\n");
2219 }
2220
2221 for (size_t i = 0; i < m_nelem; ++i) {
2223 }
2224 double l2before = 0.0;
2225 for (size_t i = 0; i < m_nelem; ++i) {
2226 l2before += x[i] * x[i];
2227 }
2228 l2before = sqrt(l2before/m_nelem);
2229
2230 // Special section to take out single species, single component,
2231 // moles. These are species which have non-zero entries in the
2232 // formula matrix, and no other species have zero values either.
2233 bool changed = false;
2234 for (size_t i = 0; i < m_nelem; ++i) {
2235 int numNonZero = 0;
2236 bool multisign = false;
2237 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
2239 double eval = m_formulaMatrix(kspec,i);
2240 if (eval < 0.0) {
2241 multisign = true;
2242 }
2243 if (eval != 0.0) {
2244 numNonZero++;
2245 }
2246 }
2247 }
2248 if (!multisign) {
2249 if (numNonZero < 2) {
2250 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
2252 double eval = m_formulaMatrix(kspec,i);
2253 if (eval > 0.0) {
2254 m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval;
2255 changed = true;
2256 }
2257 }
2258 }
2259 } else {
2260 int numCompNonZero = 0;
2261 size_t compID = npos;
2262 for (size_t kspec = 0; kspec < m_numComponents; kspec++) {
2264 double eval = m_formulaMatrix(kspec,i);
2265 if (eval > 0.0) {
2266 compID = kspec;
2267 numCompNonZero++;
2268 }
2269 }
2270 }
2271 if (numCompNonZero == 1) {
2272 double diff = m_elemAbundancesGoal[i];
2273 for (size_t kspec = m_numComponents; kspec < m_nsp; kspec++) {
2275 double eval = m_formulaMatrix(kspec,i);
2276 diff -= eval * m_molNumSpecies_old[kspec];
2277 }
2278 m_molNumSpecies_old[compID] = std::max(0.0,diff/m_formulaMatrix(compID,i));
2279 changed = true;
2280 }
2281 }
2282 }
2283 }
2284 }
2285 if (changed) {
2286 vcs_elab();
2287 }
2288
2289 // Section to check for maximum bounds errors on all species due to
2290 // elements. This may only be tried on element types which are
2291 // VCS_ELEM_TYPE_ABSPOS. This is because no other species may have a
2292 // negative number of these.
2293 //
2294 // Note, also we can do this over ne, the number of elements, not just the
2295 // number of components.
2296 changed = false;
2297 for (size_t i = 0; i < m_nelem; ++i) {
2298 int elType = m_elType[i];
2299 if (elType == VCS_ELEM_TYPE_ABSPOS) {
2300 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
2302 double atomComp = m_formulaMatrix(kspec,i);
2303 if (atomComp > 0.0) {
2304 double maxPermissible = m_elemAbundancesGoal[i] / atomComp;
2305 if (m_molNumSpecies_old[kspec] > maxPermissible) {
2306 if (m_debug_print_lvl >= 3) {
2307 plogf(" --- vcs_elcorr: Reduced species %s from %g to %g "
2308 "due to %s max bounds constraint\n",
2309 m_speciesName[kspec], m_molNumSpecies_old[kspec],
2310 maxPermissible, m_elementName[i]);
2311 }
2312 m_molNumSpecies_old[kspec] = maxPermissible;
2313 changed = true;
2315 m_molNumSpecies_old[kspec] = 0.0;
2316 if (m_SSPhase[kspec]) {
2318 } else {
2320 }
2321 if (m_debug_print_lvl >= 2) {
2322 plogf(" --- vcs_elcorr: Zeroed species %s and changed "
2323 "status to %d due to max bounds constraint\n",
2324 m_speciesName[kspec], m_speciesStatus[kspec]);
2325 }
2326 }
2327 }
2328 }
2329 }
2330 }
2331 }
2332 }
2333
2334 // Recalculate the element abundances if something has changed.
2335 if (changed) {
2336 vcs_elab();
2337 }
2338
2339 // Ok, do the general case. Linear algebra problem is of length nc, not ne,
2340 // as there may be degenerate rows when nc .ne. ne.
2342 for (size_t i = 0; i < m_numComponents; ++i) {
2344 if (fabs(x[i]) > 1.0E-13) {
2345 retn = 1;
2346 }
2347 for (size_t j = 0; j < m_numComponents; ++j) {
2348 A(j, i) = - m_formulaMatrix(i,j);
2349 }
2350 }
2351
2352 solve(A, x, 1, m_nelem);
2353
2354 // Now apply the new direction without creating negative species.
2355 double par = 0.5;
2356 for (size_t i = 0; i < m_numComponents; ++i) {
2357 if (m_molNumSpecies_old[i] > 0.0) {
2358 par = std::max(par, -x[i] / m_molNumSpecies_old[i]);
2359 }
2360 }
2361 par = std::min(par, 100.0);
2362 par = 1.0 / par;
2363 if (par < 1.0 && par > 0.0) {
2364 retn = 2;
2365 par *= 0.9999;
2366 for (size_t i = 0; i < m_numComponents; ++i) {
2367 double tmp = m_molNumSpecies_old[i] + par * x[i];
2368 if (tmp > 0.0) {
2369 m_molNumSpecies_old[i] = tmp;
2370 } else {
2371 if (m_SSPhase[i]) {
2372 m_molNumSpecies_old[i] = 0.0;
2373 } else {
2375 }
2376 }
2377 }
2378 } else {
2379 for (size_t i = 0; i < m_numComponents; ++i) {
2380 double tmp = m_molNumSpecies_old[i] + x[i];
2381 if (tmp > 0.0) {
2382 m_molNumSpecies_old[i] = tmp;
2383 } else {
2384 if (m_SSPhase[i]) {
2385 m_molNumSpecies_old[i] = 0.0;
2386 } else {
2388 }
2389 }
2390 }
2391 }
2392
2393 // We have changed the element abundances. Calculate them again
2394 vcs_elab();
2395
2396 // We have changed the total moles in each phase. Calculate them again
2397 vcs_tmoles();
2398
2399 // Try some ad hoc procedures for fixing the problem
2400 if (retn >= 2) {
2401 // First find a species whose adjustment is a win-win situation.
2402 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
2404 continue;
2405 }
2406 double saveDir = 0.0;
2407 bool goodSpec = true;
2408 for (size_t i = 0; i < m_numComponents; ++i) {
2409 double dir = m_formulaMatrix(kspec,i) * (m_elemAbundancesGoal[i] - m_elemAbundances[i]);
2410 if (fabs(dir) > 1.0E-10) {
2411 if (dir > 0.0) {
2412 if (saveDir < 0.0) {
2413 goodSpec = false;
2414 break;
2415 }
2416 } else {
2417 if (saveDir > 0.0) {
2418 goodSpec = false;
2419 break;
2420 }
2421 }
2422 saveDir = dir;
2423 } else {
2424 if (m_formulaMatrix(kspec,i) != 0.) {
2425 goodSpec = false;
2426 break;
2427 }
2428 }
2429 }
2430 if (goodSpec) {
2431 int its = 0;
2432 double xx = 0.0;
2433 for (size_t i = 0; i < m_numComponents; ++i) {
2434 if (m_formulaMatrix(kspec,i) != 0.0) {
2435 xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix(kspec,i);
2436 its++;
2437 }
2438 }
2439 if (its > 0) {
2440 xx /= its;
2441 }
2442 m_molNumSpecies_old[kspec] += xx;
2443 m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 1.0E-10);
2444
2445 // If we are dealing with a deleted species, then we need to
2446 // reinsert it into the active list.
2447 if (kspec >= m_numSpeciesRdc) {
2448 vcs_reinsert_deleted(kspec);
2450 vcs_elab();
2451 goto L_CLEANUP;
2452 }
2453 vcs_elab();
2454 }
2455 }
2456 }
2457 if (vcs_elabcheck(0)) {
2458 retn = 1;
2459 goto L_CLEANUP;
2460 }
2461
2462 for (size_t i = 0; i < m_nelem; ++i) {
2464 (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) {
2465 for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
2466 if (m_elemAbundances[i] > 0.0 && m_formulaMatrix(kspec,i) < 0.0) {
2468 m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
2469 vcs_elab();
2470 break;
2471 }
2472 if (m_elemAbundances[i] < 0.0 && m_formulaMatrix(kspec,i) > 0.0) {
2474 m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
2475 vcs_elab();
2476 break;
2477 }
2478 }
2479 }
2480 }
2481 if (vcs_elabcheck(1)) {
2482 retn = 1;
2483 goto L_CLEANUP;
2484 }
2485
2486 // For electron charges element types, we try positive deltas in the species
2487 // concentrations to match the desired electron charge exactly.
2488 for (size_t i = 0; i < m_nelem; ++i) {
2489 double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i];
2490 if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) {
2491 bool useZeroed = true;
2492 for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
2493 if (dev < 0.0) {
2494 if (m_formulaMatrix(kspec,i) < 0.0 && m_molNumSpecies_old[kspec] > 0.0) {
2495 useZeroed = false;
2496 }
2497 } else {
2498 if (m_formulaMatrix(kspec,i) > 0.0 && m_molNumSpecies_old[kspec] > 0.0) {
2499 useZeroed = false;
2500 }
2501 }
2502 }
2503 for (size_t kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
2504 if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) {
2505 if (dev < 0.0 && m_formulaMatrix(kspec,i) < 0.0) {
2506 double delta = dev / m_formulaMatrix(kspec,i);
2507 m_molNumSpecies_old[kspec] += delta;
2508 m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
2509 vcs_elab();
2510 break;
2511 }
2512 if (dev > 0.0 && m_formulaMatrix(kspec,i) > 0.0) {
2513 double delta = dev / m_formulaMatrix(kspec,i);
2514 m_molNumSpecies_old[kspec] += delta;
2515 m_molNumSpecies_old[kspec] = std::max(m_molNumSpecies_old[kspec], 0.0);
2516 vcs_elab();
2517 break;
2518 }
2519 }
2520 }
2521 }
2522 }
2523 if (vcs_elabcheck(1)) {
2524 retn = 1;
2525 goto L_CLEANUP;
2526 }
2527
2528L_CLEANUP:
2529 ;
2530 vcs_tmoles();
2531 double l2after = 0.0;
2532 for (size_t i = 0; i < m_nelem; ++i) {
2533 l2after += pow(m_elemAbundances[i] - m_elemAbundancesGoal[i], 2);
2534 }
2535 l2after = sqrt(l2after/m_nelem);
2536 if (m_debug_print_lvl >= 2) {
2537 plogf(" --- Elem_Abund: Correct Initial "
2538 " Final\n");
2539 for (size_t i = 0; i < m_nelem; ++i) {
2540 plogf(" --- ");
2541 plogf("%-2.2s", m_elementName[i]);
2542 plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]);
2543 }
2544 plogf(" --- Diff_Norm: %20.12E %20.12E\n",
2545 l2before, l2after);
2546 }
2547 return retn;
2548}
2549
2551{
2552 const char* pprefix = " --- vcs_inest: ";
2553 int retn = 0;
2554 clockWC tickTock;
2555 if (m_doEstimateEquil > 0) {
2556 // Calculate the elemental abundances
2557 vcs_elab();
2558 if (vcs_elabcheck(0)) {
2559 if (m_debug_print_lvl >= 2) {
2560 plogf("%s Initial guess passed element abundances on input\n", pprefix);
2561 plogf("%s m_doEstimateEquil = 1 so will use the input mole "
2562 "numbers as estimates\n", pprefix);
2563 }
2564 return retn;
2565 } else if (m_debug_print_lvl >= 2) {
2566 plogf("%s Initial guess failed element abundances on input\n", pprefix);
2567 plogf("%s m_doEstimateEquil = 1 so will discard input "
2568 "mole numbers and find our own estimate\n", pprefix);
2569 }
2570 }
2571
2572 // temporary space for usage in this routine and in subroutines
2573 vector<double> sm(m_nelem*m_nelem, 0.0);
2574 vector<double> ss(m_nelem, 0.0);
2575 vector<double> sa(m_nelem, 0.0);
2576 vector<double> aw(m_nsp + m_nelem, 0.0);
2577
2578 // Go get the estimate of the solution
2579 if (m_debug_print_lvl >= 2) {
2580 plogf("%sGo find an initial estimate for the equilibrium problem\n",
2581 pprefix);
2582 }
2583 double test = -1.0E20;
2584 vcs_inest(&aw[0], &sa[0], &sm[0], &ss[0], test);
2585
2586 // Calculate the elemental abundances
2587 vcs_elab();
2588
2589 // If we still fail to achieve the correct elemental abundances, try to fix
2590 // the problem again by calling the main elemental abundances fixer routine,
2591 // used in the main program. This attempts to tweak the mole numbers of the
2592 // component species to satisfy the element abundance constraints.
2593 //
2594 // Note: We won't do this unless we have to since it involves inverting a
2595 // matrix.
2596 bool rangeCheck = vcs_elabcheck(1);
2597 if (!vcs_elabcheck(0)) {
2598 if (m_debug_print_lvl >= 2) {
2599 plogf("%sInitial guess failed element abundances\n", pprefix);
2600 plogf("%sCall vcs_elcorr to attempt fix\n", pprefix);
2601 }
2602 vcs_elcorr(&sm[0], &aw[0]);
2603 rangeCheck = vcs_elabcheck(1);
2604 if (!vcs_elabcheck(0)) {
2605 plogf("%sInitial guess still fails element abundance equations\n",
2606 pprefix);
2607 plogf("%s - Inability to ever satisfy element abundance "
2608 "constraints is probable\n", pprefix);
2609 retn = -1;
2610 } else {
2611 if (m_debug_print_lvl >= 2) {
2612 if (rangeCheck) {
2613 plogf("%sInitial guess now satisfies element abundances\n", pprefix);
2614 } else {
2615 plogf("%sElement Abundances RANGE ERROR\n", pprefix);
2616 plogf("%s - Initial guess satisfies NC=%d element abundances, "
2617 "BUT not NE=%d element abundances\n", pprefix,
2619 }
2620 }
2621 }
2622 } else {
2623 if (m_debug_print_lvl >= 2) {
2624 if (rangeCheck) {
2625 plogf("%sInitial guess satisfies element abundances\n", pprefix);
2626 } else {
2627 plogf("%sElement Abundances RANGE ERROR\n", pprefix);
2628 plogf("%s - Initial guess satisfies NC=%d element abundances, "
2629 "BUT not NE=%d element abundances\n", pprefix,
2631 }
2632 }
2633 }
2634
2635 if (m_debug_print_lvl >= 2) {
2636 plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E\n", pprefix,
2638 &m_tPhaseMoles_old[0]));
2639 }
2640
2641 // Record time
2642 m_VCount->T_Time_inest += tickTock.secondsWC();
2644 return retn;
2645}
2646
2648{
2649 double test = -1.0E-10;
2650
2651 if (m_debug_print_lvl >= 2) {
2652 plogf(" --- call setInitialMoles\n");
2653 }
2654
2655 double dxi_min = 1.0e10;
2656 int retn;
2657 int iter = 0;
2658 bool abundancesOK = true;
2659 bool usedZeroedSpecies;
2660 vector<double> sm(m_nelem * m_nelem, 0.0);
2661 vector<double> ss(m_nelem, 0.0);
2662 vector<double> sa(m_nelem, 0.0);
2663 vector<double> wx(m_nelem, 0.0);
2664 vector<double> aw(m_nsp, 0.0);
2665
2666 for (size_t ik = 0; ik < m_nsp; ik++) {
2668 m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]);
2669 }
2670 }
2671
2672 if (m_debug_print_lvl >= 2) {
2674 }
2675
2676 bool redo = true;
2677 while (redo) {
2678 if (!vcs_elabcheck(0)) {
2679 if (m_debug_print_lvl >= 2) {
2680 plogf(" --- seMolesLinProg Mole numbers failing element abundances\n");
2681 plogf(" --- seMolesLinProg Call vcs_elcorr to attempt fix\n");
2682 }
2683 retn = vcs_elcorr(&sm[0], &wx[0]);
2684 if (retn >= 2) {
2685 abundancesOK = false;
2686 } else {
2687 abundancesOK = true;
2688 }
2689 } else {
2690 abundancesOK = true;
2691 }
2692
2693 // Now find the optimized basis that spans the stoichiometric
2694 // coefficient matrix, based on the current composition,
2695 // m_molNumSpecies_old[] We also calculate sc[][], the reaction matrix.
2696 retn = vcs_basopt(false, &aw[0], &sa[0], &sm[0], &ss[0],
2697 test, &usedZeroedSpecies);
2698 if (retn != VCS_SUCCESS) {
2699 return retn;
2700 }
2701
2702 if (m_debug_print_lvl >= 2) {
2703 plogf("iteration %d\n", iter);
2704 }
2705 redo = false;
2706 iter++;
2707 if (iter > 15) {
2708 break;
2709 }
2710
2711 // loop over all reactions
2712 for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
2713 // dg_rt is the Delta_G / RT value for the reaction
2714 size_t ik = m_numComponents + irxn;
2715 double dg_rt = m_SSfeSpecies[ik];
2716 dxi_min = 1.0e10;
2717 const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
2718 for (size_t jcomp = 0; jcomp < m_nelem; jcomp++) {
2719 dg_rt += m_SSfeSpecies[jcomp] * sc_irxn[jcomp];
2720 }
2721 // fwd or rev direction.
2722 // idir > 0 implies increasing the current species
2723 // idir < 0 implies decreasing the current species
2724 int idir = (dg_rt < 0.0 ? 1 : -1);
2725 if (idir < 0) {
2726 dxi_min = m_molNumSpecies_old[ik];
2727 }
2728
2729 for (size_t jcomp = 0; jcomp < m_numComponents; jcomp++) {
2730 double nu = sc_irxn[jcomp];
2731 // set max change in progress variable by
2732 // non-negativity requirement
2733 if (nu*idir < 0) {
2734 double delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu);
2735 // if a component has nearly zero moles, redo
2736 // with a new set of components
2737 if (!redo && delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) {
2738 if (m_debug_print_lvl >= 2) {
2739 plogf(" --- Component too small: %s\n", m_speciesName[jcomp]);
2740 }
2741 redo = true;
2742 }
2743 dxi_min = std::min(dxi_min, delta_xi);
2744 }
2745 }
2746
2747 // step the composition by dxi_min, check against zero, since
2748 // we are zeroing components and species on every step.
2749 // Redo the iteration, if a component went from positive to zero on this step.
2750 double dsLocal = idir*dxi_min;
2751 m_molNumSpecies_old[ik] += dsLocal;
2752 m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]);
2753 for (size_t jcomp = 0; jcomp < m_numComponents; jcomp++) {
2754 bool full = false;
2755 if (m_molNumSpecies_old[jcomp] > 1.0E-15) {
2756 full = true;
2757 }
2758 m_molNumSpecies_old[jcomp] += sc_irxn[jcomp] * dsLocal;
2759 m_molNumSpecies_old[jcomp] = max(0.0, m_molNumSpecies_old[jcomp]);
2760 if (full && m_molNumSpecies_old[jcomp] < 1.0E-60) {
2761 redo = true;
2762 }
2763 }
2764 }
2765
2766 if (m_debug_print_lvl >= 2) {
2768 }
2769 }
2770
2771 if (m_debug_print_lvl == 1) {
2773 plogf(" --- setInitialMoles end\n");
2774 }
2775 retn = 0;
2776 if (!abundancesOK) {
2777 retn = -1;
2778 } else if (iter > 15) {
2779 retn = 1;
2780 }
2781 return retn;
2782}
2783
2784double VCS_SOLVE::vcs_Total_Gibbs(double* molesSp, double* chemPot,
2785 double* tPhMoles)
2786{
2787 double g = 0.0;
2788
2789 for (size_t iph = 0; iph < m_numPhases; iph++) {
2790 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
2791 if ((TPhInertMoles[iph] > 0.0) && (tPhMoles[iph] > 0.0)) {
2792 g += TPhInertMoles[iph] *
2793 log(TPhInertMoles[iph] / tPhMoles[iph]);
2794 if (Vphase->m_gasPhase) {
2795 g += TPhInertMoles[iph] * log(m_pressurePA/(1.01325E5));
2796 }
2797 }
2798 }
2799
2800 for (size_t kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
2802 g += molesSp[kspec] * chemPot[kspec];
2803 }
2804 }
2805
2806 return g;
2807}
2808
2809double VCS_SOLVE::vcs_GibbsPhase(size_t iphase, const double* const w,
2810 const double* const fe)
2811{
2812 double g = 0.0;
2813 double phaseMols = 0.0;
2814 for (size_t kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
2815 if (m_phaseID[kspec] == iphase && m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
2816 g += w[kspec] * fe[kspec];
2817 phaseMols += w[kspec];
2818 }
2819 }
2820
2821 if (TPhInertMoles[iphase] > 0.0) {
2822 phaseMols += TPhInertMoles[iphase];
2823 g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols);
2824 vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
2825 if (Vphase->m_gasPhase) {
2826 g += TPhInertMoles[iphase] * log(m_pressurePA/1.01325E5);
2827 }
2828 }
2829
2830 return g;
2831}
2832
2834{
2835 // Transfer the information back to the MultiPhase object. Note we don't
2836 // just call setMoles, because some multispecies solution phases may be
2837 // zeroed out, and that would cause a problem for that routine. Also, the
2838 // mole fractions of such zeroed out phases actually contain information
2839 // about likely reemergent states.
2841 for (size_t ip = 0; ip < m_numPhases; ip++) {
2842 m_mix->setPhaseMoles(ip, m_VolPhaseList[ip]->totalMoles());
2843 }
2844}
2845
2847{
2848 // Whether we have an estimate or not gets overwritten on
2849 // the call to the equilibrium solver.
2850 m_temperature = m_mix->temperature();
2851 m_pressurePA = m_mix->pressure();
2853 m_totalVol = m_mix->volume();
2854
2855 vector<size_t> invSpecies(m_nsp);
2856 for (size_t k = 0; k < m_nsp; k++) {
2857 invSpecies[m_speciesMapIndex[k]] = k;
2858 }
2859
2860 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
2861 ThermoPhase* tPhase = &m_mix->phase(iphase);
2862 vcs_VolPhase* volPhase = m_VolPhaseList[iphase].get();
2863
2865
2866 // Loop through each species in the current phase
2867 size_t nSpPhase = tPhase->nSpecies();
2868 if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) {
2870 } else if (volPhase->totalMoles() > 0.0) {
2872 } else {
2874 }
2875 }
2876
2877 // Printout the species information: PhaseID's and mole nums
2878 if (m_printLvl > 1) {
2879 writeline('=', 80, true, true);
2880 writeline('=', 20, false);
2881 plogf(" Cantera_to_vprob: START OF PROBLEM STATEMENT ");
2882 writeline('=', 20);
2883 writeline('=', 80);
2884 plogf("\n");
2885 plogf(" Phase IDs of species\n");
2886 plogf(" species phaseID phaseName ");
2887 plogf(" Initial_Estimated_kMols\n");
2888 for (size_t i = 0; i < m_nsp; i++) {
2889 size_t iphase = m_phaseID[i];
2890 vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
2891 plogf("%16s %5d %16s", m_speciesName[i].c_str(), iphase,
2892 VolPhase->PhaseName.c_str());
2894 plogf(" Volts = %-10.5g\n", m_molNumSpecies_old[i]);
2895 } else {
2896 plogf(" %-10.5g\n", m_molNumSpecies_old[i]);
2897 }
2898 }
2899
2900 // Printout of the Phase structure information
2901 writeline('-', 80, true, true);
2902 plogf(" Information about phases\n");
2903 plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
2904 plogf(" TMolesInert Tmoles(kmol)\n");
2905
2906 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
2907 vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
2908 plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
2909 VolPhase->VP_ID_, VolPhase->m_singleSpecies,
2910 VolPhase->m_gasPhase, VolPhase->eos_name(),
2911 VolPhase->nSpecies(), VolPhase->totalMolesInert());
2912 plogf("%16e\n", VolPhase->totalMoles());
2913 }
2914
2915 writeline('=', 80, true, true);
2916 writeline('=', 20, false);
2917 plogf(" Cantera_to_vprob: END OF PROBLEM STATEMENT ");
2918 writeline('=', 20);
2919 writeline('=', 80);
2920 plogf("\n");
2921 }
2922
2923 // m_numRxnTot = number of noncomponents, also equal to the number of
2924 // reactions. Note, it's possible that the number of elements is greater
2925 // than the number of species. In that case set the number of reactions to
2926 // zero.
2927 if (m_nelem > m_nsp) {
2928 m_numRxnTot = 0;
2929 } else {
2931 }
2933}
2934
2935void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
2936 double* const ss, double test)
2937{
2938 const char* pprefix = " --- vcs_inest: ";
2939 size_t nrxn = m_numRxnTot;
2940
2941 // CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB AND
2942 // molNum(I) .GE. 0.0. Note, both of these programs do this.
2944
2945 if (m_debug_print_lvl >= 2) {
2946 plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
2947 pprefix);
2948 plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
2949 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
2950 plogf("%s ", pprefix);
2951 plogf("%-12.12s", m_speciesName[kspec]);
2952 plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]);
2953 }
2954 plogf("%s Element Abundance Agreement returned from linear "
2955 "programming (vcs_inest initial guess):\n", pprefix);
2956 plogf("%s Element Goal Actual\n", pprefix);
2957 for (size_t j = 0; j < m_nelem; j++) {
2958 if (m_elementActive[j]) {
2959 double tmp = 0.0;
2960 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
2961 tmp += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec];
2962 }
2963 plogf("%s ", pprefix);
2964 plogf(" %-9.9s", m_elementName[j]);
2965 plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp);
2966 }
2967 }
2968 writelogendl();
2969 }
2970
2971 // Make sure all species have positive definite mole numbers Set voltages to
2972 // zero for now, until we figure out what to do
2973 m_deltaMolNumSpecies.assign(m_deltaMolNumSpecies.size(), 0.0);
2974 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
2976 if (m_molNumSpecies_old[kspec] <= 0.0) {
2977 // HKM Should eventually include logic here for non SS phases
2978 if (!m_SSPhase[kspec]) {
2979 m_molNumSpecies_old[kspec] = 1.0e-30;
2980 }
2981 }
2982 } else {
2983 m_molNumSpecies_old[kspec] = 0.0;
2984 }
2985 }
2986
2987 // Now find the optimized basis that spans the stoichiometric coefficient
2988 // matrix
2989 bool conv;
2990 vcs_basopt(false, aw, sa, sm, ss, test, &conv);
2991
2992 // CALCULATE TOTAL MOLES, CHEMICAL POTENTIALS OF BASIS
2993
2994 // Calculate TMoles and m_tPhaseMoles_old[]
2995 vcs_tmoles();
2996
2997 // m_tPhaseMoles_new[] will consist of just the component moles
2998 for (size_t iph = 0; iph < m_numPhases; iph++) {
2999 m_tPhaseMoles_new[iph] = TPhInertMoles[iph] + 1.0E-20;
3000 }
3001 for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3004 }
3005 }
3006 double TMolesMultiphase = 0.0;
3007 for (size_t iph = 0; iph < m_numPhases; iph++) {
3008 if (! m_VolPhaseList[iph]->m_singleSpecies) {
3009 TMolesMultiphase += m_tPhaseMoles_new[iph];
3010 }
3011 }
3013 for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3015 m_molNumSpecies_new[kspec] = 0.0;
3016 }
3017 }
3019
3020 for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3022 if (! m_SSPhase[kspec]) {
3023 size_t iph = m_phaseID[kspec];
3024 m_feSpecies_new[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
3025 }
3026 } else {
3027 m_molNumSpecies_new[kspec] = 0.0;
3028 }
3029 }
3030 vcs_deltag(0, true, VCS_STATECALC_NEW);
3031 if (m_debug_print_lvl >= 2) {
3032 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
3033 plogf("%s", pprefix);
3034 plogf("%-12.12s", m_speciesName[kspec]);
3035 if (kspec < m_numComponents) {
3036 plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_new[kspec],
3037 m_SSfeSpecies[kspec]);
3038 } else {
3039 plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n",
3041 }
3042 }
3043 }
3044
3045 // ESTIMATE REACTION ADJUSTMENTS
3046 vector<double>& xtphMax = m_TmpPhase;
3047 vector<double>& xtphMin = m_TmpPhase2;
3048 m_deltaPhaseMoles.assign(m_deltaPhaseMoles.size(), 0.0);
3049 for (size_t iph = 0; iph < m_numPhases; iph++) {
3050 xtphMax[iph] = log(m_tPhaseMoles_new[iph] * 1.0E32);
3051 xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32);
3052 }
3053 for (size_t irxn = 0; irxn < nrxn; ++irxn) {
3054 size_t kspec = m_indexRxnToSpecies[irxn];
3055
3056 // For single species phases, we will not estimate the mole numbers. If
3057 // the phase exists, it stays. If it doesn't exist in the estimate, it
3058 // doesn't come into existence here.
3059 if (! m_SSPhase[kspec]) {
3060 size_t iph = m_phaseID[kspec];
3061 if (m_deltaGRxn_new[irxn] > xtphMax[iph]) {
3062 m_deltaGRxn_new[irxn] = 0.8 * xtphMax[iph];
3063 }
3064 if (m_deltaGRxn_new[irxn] < xtphMin[iph]) {
3065 m_deltaGRxn_new[irxn] = 0.8 * xtphMin[iph];
3066 }
3067
3068 // HKM -> The TMolesMultiphase is a change of mine. It more evenly
3069 // distributes the initial moles amongst multiple multispecies
3070 // phases according to the relative values of the standard state
3071 // free energies. There is no change for problems with one
3072 // multispecies phase. It cut diamond4.vin iterations down from 62
3073 // to 14.
3074 m_deltaMolNumSpecies[kspec] = 0.5 * (m_tPhaseMoles_new[iph] + TMolesMultiphase)
3075 * exp(-m_deltaGRxn_new[irxn]);
3076
3077 for (size_t k = 0; k < m_numComponents; ++k) {
3079 }
3080
3081 for (iph = 0; iph < m_numPhases; iph++) {
3083 }
3084 }
3085 }
3086 if (m_debug_print_lvl >= 2) {
3087 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
3089 plogf("%sdirection (", pprefix);
3090 plogf("%-12.12s", m_speciesName[kspec]);
3091 plogf(") = %g", m_deltaMolNumSpecies[kspec]);
3092 if (m_SSPhase[kspec]) {
3093 if (m_molNumSpecies_old[kspec] > 0.0) {
3094 plogf(" (ssPhase exists at w = %g moles)", m_molNumSpecies_old[kspec]);
3095 } else {
3096 plogf(" (ssPhase doesn't exist -> stability not checked)");
3097 }
3098 }
3099 writelogendl();
3100 }
3101 }
3102 }
3103
3104 // KEEP COMPONENT SPECIES POSITIVE
3105 double par = 0.5;
3106 for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3108 par < -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]) {
3109 par = -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec];
3110 }
3111 }
3112 par = 1. / par;
3113 if (par <= 1.0 && par > 0.0) {
3114 par *= 0.8;
3115 } else {
3116 par = 1.0;
3117 }
3118
3119 // CALCULATE NEW MOLE NUMBERS
3120 size_t lt = 0;
3121 size_t ikl = 0;
3122 double s1 = 0.0;
3123 while (true) {
3124 for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
3126 m_molNumSpecies_old[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec];
3127 } else {
3128 m_deltaMolNumSpecies[kspec] = 0.0;
3129 }
3130 }
3131 for (size_t kspec = m_numComponents; kspec < m_nsp; ++kspec) {
3133 m_deltaMolNumSpecies[kspec] != 0.0) {
3134 m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par;
3135 }
3136 }
3137
3138 // We have a new w[] estimate, go get the TMoles and m_tPhaseMoles_old[]
3139 // values
3140 vcs_tmoles();
3141 if (lt > 0) {
3142 break;
3143 }
3144
3145 // CONVERGENCE FORCING SECTION
3146 vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
3148 double s = 0.0;
3149 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
3150 s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec];
3151 }
3152 if (s == 0.0) {
3153 break;
3154 }
3155 if (s < 0.0 && ikl == 0) {
3156 break;
3157 }
3158
3159 // TRY HALF STEP SIZE
3160 if (ikl == 0) {
3161 s1 = s;
3162 par *= 0.5;
3163 ikl = 1;
3164 continue;
3165 }
3166
3167 // FIT PARABOLA THROUGH HALF AND FULL STEPS
3168 double xl = (1.0 - s / (s1 - s)) * 0.5;
3169 if (xl < 0.0) {
3170 // POOR DIRECTION, REDUCE STEP SIZE TO 0.2
3171 par *= 0.2;
3172 } else {
3173 if (xl > 1.0) {
3174 // TOO BIG A STEP, TAKE ORIGINAL FULL STEP
3175 par *= 2.0;
3176 } else {
3177 // ACCEPT RESULTS OF FORCER
3178 par = par * 2.0 * xl;
3179 }
3180 }
3181 lt = 1;
3182 }
3183
3184 if (m_debug_print_lvl >= 2) {
3185 plogf("%s Final Mole Numbers produced by inest:\n",
3186 pprefix);
3187 plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
3188 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
3189 plogf("%s %-12.12s %g\n",
3190 pprefix, m_speciesName[kspec], m_molNumSpecies_old[kspec]);
3191 }
3192 }
3193}
3194
3196{
3197 vector<int> numPhSpecies(m_numPhases, 0);
3198 for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
3199 numPhSpecies[m_phaseID[kspec]]++;
3200 }
3201
3202 // Handle the special case of a single species in a phase that has been
3203 // earmarked as a multispecies phase. Treat that species as a single-species
3204 // phase
3205 for (size_t iph = 0; iph < m_numPhases; iph++) {
3206 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
3207 Vphase->m_singleSpecies = false;
3208 if (TPhInertMoles[iph] > 0.0) {
3209 Vphase->setExistence(2);
3210 }
3211 if (numPhSpecies[iph] <= 1 && TPhInertMoles[iph] == 0.0) {
3212 Vphase->m_singleSpecies = true;
3213 }
3214 }
3215
3216 // Fill in some useful arrays here that have to do with the static
3217 // information concerning the phase ID of species. SSPhase = Boolean
3218 // indicating whether a species is in a single species phase or not.
3219 for (size_t kspec = 0; kspec < m_nsp; kspec++) {
3220 size_t iph = m_phaseID[kspec];
3221 vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
3222 if (Vphase->m_singleSpecies) {
3223 m_SSPhase[kspec] = true;
3224 } else {
3225 m_SSPhase[kspec] = false;
3226 }
3227 }
3228}
3229
3231{
3232 delete m_VCount;
3233 m_VCount = 0;
3234
3235 m_nsp = 0;
3236 m_nelem = 0;
3237 m_numComponents = 0;
3238
3239}
3240
3242{
3243 m_VCount->Its = 0;
3244 m_VCount->Basis_Opts = 0;
3245 m_VCount->Time_vcs_TP = 0.0;
3246 m_VCount->Time_basopt = 0.0;
3247 if (ifunc) {
3248 m_VCount->T_Its = 0;
3252 m_VCount->T_Time_vcs_TP = 0.0;
3253 m_VCount->T_Time_basopt = 0.0;
3254 m_VCount->T_Time_inest = 0.0;
3255 m_VCount->T_Time_vcs = 0.0;
3256 }
3257}
3258
3259void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl)
3260{
3261 plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n");
3262 if (timing_print_lvl > 0) {
3263 plogf(" vcs_basopt: %5d %5d %11.5E\n",
3266 plogf(" vcs_TP: %5d %5d %11.5E\n",
3269 plogf(" vcs_inest: %5d %11.5E\n",
3271 plogf(" vcs_TotalTime: %11.5E\n",
3273 } else {
3274 plogf(" vcs_basopt: %5d %5d %11s\n",
3276 plogf(" vcs_TP: %5d %5d %11s\n",
3278 plogf(" vcs_inest: %5d %11s\n",
3279 m_VCount->T_Calls_Inest, " NA ");
3280 plogf(" vcs_TotalTime: %11s\n",
3281 " NA ");
3282 }
3283}
3284
3285void VCS_SOLVE::prob_report(int print_lvl)
3286{
3287 m_printLvl = print_lvl;
3288
3289 // Printout the species information: PhaseID's and mole nums
3290 if (m_printLvl > 0) {
3291 writeline('=', 80, true, true);
3292 writeline('=', 20, false);
3293 plogf(" VCS_PROB: PROBLEM STATEMENT ");
3294 writeline('=', 31);
3295 writeline('=', 80);
3296 plogf("\n");
3297 plogf("\tSolve a constant T, P problem:\n");
3298 plogf("\t\tT = %g K\n", m_temperature);
3299 double pres_atm = m_pressurePA / 1.01325E5;
3300
3301 plogf("\t\tPres = %g atm\n", pres_atm);
3302 plogf("\n");
3303 plogf(" Phase IDs of species\n");
3304 plogf(" species phaseID phaseName ");
3305 plogf(" Initial_Estimated_Moles Species_Type\n");
3306 for (size_t i = 0; i < m_nsp; i++) {
3307 vcs_VolPhase* Vphase = m_VolPhaseList[m_phaseID[i]].get();
3308 plogf("%16s %5d %16s", m_mix->speciesName(i), m_phaseID[i],
3309 Vphase->PhaseName);
3310 if (m_doEstimateEquil >= 0) {
3311 plogf(" %-10.5g", m_molNumSpecies_old[i]);
3312 } else {
3313 plogf(" N/A");
3314 }
3316 plogf(" Mol_Num");
3318 plogf(" Voltage");
3319 } else {
3320 plogf(" ");
3321 }
3322 plogf("\n");
3323 }
3324
3325 // Printout of the Phase structure information
3326 writeline('-', 80, true, true);
3327 plogf(" Information about phases\n");
3328 plogf(" PhaseName PhaseNum SingSpec GasPhase "
3329 " EqnState NumSpec");
3330 plogf(" TMolesInert TKmoles\n");
3331
3332 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
3333 vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
3334 plogf("%16s %5d %5d %8d ", Vphase->PhaseName,
3335 Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase);
3336 plogf("%16s %8d %16e ", Vphase->eos_name(),
3337 Vphase->nSpecies(), Vphase->totalMolesInert());
3338 if (m_doEstimateEquil >= 0) {
3339 plogf("%16e\n", Vphase->totalMoles());
3340 } else {
3341 plogf(" N/A\n");
3342 }
3343 }
3344
3345 plogf("\nElemental Abundances: ");
3346 plogf(" Target_kmol ElemType ElActive\n");
3347 for (size_t i = 0; i < m_nelem; ++i) {
3348 writeline(' ', 26, false);
3349 plogf("%-2.2s", m_elementName[i]);
3350 plogf("%20.12E ", m_elemAbundancesGoal[i]);
3351 plogf("%3d %3d\n", m_elType[i], m_elementActive[i]);
3352 }
3353
3354 plogf("\nChemical Potentials: (J/kmol)\n");
3355 plogf(" Species (phase) "
3356 " SS0ChemPot StarChemPot\n");
3357 for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
3358 vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
3360 for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) {
3361 size_t kglob = Vphase->spGlobalIndexVCS(kindex);
3362 plogf("%16s ", m_mix->speciesName(kglob));
3363 if (kindex == 0) {
3364 plogf("%16s", Vphase->PhaseName);
3365 } else {
3366 plogf(" ");
3367 }
3368
3369 plogf("%16g %16g\n", Vphase->G0_calc_one(kindex),
3370 Vphase->GStar_calc_one(kindex));
3371 }
3372 }
3373 writeline('=', 80, true, true);
3374 writeline('=', 20, false);
3375 plogf(" VCS_PROB: END OF PROBLEM STATEMENT ");
3376 writeline('=', 24);
3377 writeline('=', 80);
3378 plogf("\n");
3379 }
3380}
3381
3383{
3384 size_t neVP = volPhase->nElemConstraints();
3385 // Loop through the elements in the vol phase object
3386 for (size_t eVP = 0; eVP < neVP; eVP++) {
3387 size_t foundPos = npos;
3388 string enVP = volPhase->elementName(eVP);
3389
3390 // Search for matches with the existing elements. If found, then fill in
3391 // the entry in the global mapping array.
3392 for (size_t e = 0; e < m_nelem; e++) {
3393 string en = m_elementName[e];
3394 if (!strcmp(enVP.c_str(), en.c_str())) {
3395 volPhase->setElemGlobalIndex(eVP, e);
3396 foundPos = e;
3397 }
3398 }
3399 if (foundPos == npos) {
3400 int elType = volPhase->elementType(eVP);
3401 int elactive = volPhase->elementActive(eVP);
3402 size_t e = addElement(enVP.c_str(), elType, elactive);
3403 volPhase->setElemGlobalIndex(eVP, e);
3404 }
3405 }
3406}
3407
3408size_t VCS_SOLVE::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
3409{
3410 if (kT > m_nsp) {
3411 // Need to expand the number of species here
3412 throw CanteraError("VCS_SOLVE::addOnePhaseSpecies", "Shouldn't be here");
3413 }
3414 const Array2D& fm = volPhase->getFormulaMatrix();
3415 for (size_t eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) {
3416 size_t e = volPhase->elemGlobalIndex(eVP);
3417 AssertThrowMsg(e != npos, "VCS_PROB::addOnePhaseSpecies",
3418 "element not found");
3419 m_formulaMatrix(kT,e) = fm(k,eVP);
3420 }
3421
3422 // Tell the phase object about the current position of the species within
3423 // the global species vector
3424 volPhase->setSpGlobalIndexVCS(k, kT);
3425 return kT;
3426}
3427
3428size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive)
3429{
3430 if (!elNameNew) {
3431 throw CanteraError("VCS_SOLVE::addElement",
3432 "error: element must have a name");
3433 }
3434 m_nelem++;
3436
3439 m_elType.push_back(elType);
3440 m_elementActive.push_back(elactive);
3441 m_elemAbundances.push_back(0.0);
3442 m_elemAbundancesGoal.push_back(0.0);
3443 m_elementMapIndex.push_back(0);
3444 m_elementName.push_back(elNameNew);
3445 return m_nelem - 1;
3446}
3447
3450}
3451
3452}
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Chemica...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
void zero()
Set all of the entries to zero.
Definition Array.h:143
size_t nColumns() const
Number of columns.
Definition Array.h:197
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition Array.h:235
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition Array.cpp:47
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
A class for multiphase mixtures.
Definition MultiPhase.h:57
double speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
double pressure() const
Pressure [Pa].
Definition MultiPhase.h:388
double temperature() const
Temperature [K].
Definition MultiPhase.h:328
string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
double volume() const
The total mixture volume [m^3].
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
size_t nElements() const
Number of elements.
Definition MultiPhase.h:97
ThermoPhase & phase(size_t n)
Return a reference to phase n.
void setPhaseMoles(const size_t n, const double moles)
Set the number of moles of phase with index n.
double elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
A species thermodynamic property manager for a phase.
virtual int reportType(size_t index) const
This utility function reports the type of parameterization used for the species with index number ind...
virtual void reportParams(size_t index, int &type, double *const c, double &minTemp, double &maxTemp, double &refPressure) const
This utility function reports back the type of parameterization and all of the parameters for the spe...
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:245
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition Phase.h:646
size_t nElements() const
Number of elements.
Definition Phase.cpp:30
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition Phase.cpp:482
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition Phase.h:638
string name() const
Return the name of the phase.
Definition Phase.cpp:20
Base class for a phase with thermodynamic properties.
double electricPotential() const
Returns the electric potential of this phase (V).
string type() const override
String indicating the thermodynamic model implemented.
virtual int activityConvention() const
This method returns the convention used in specification of the activities, of which there are curren...
virtual MultiSpeciesThermo & speciesThermo(int k=-1)
Return a changeable reference to the calculation manager for species reference-state thermodynamic pr...
Class to keep track of time and iterations.
int Its
Current number of iterations in the main loop of vcs_TP() to solve for thermo equilibrium.
int T_Calls_Inest
Current number of times the initial thermo equilibrium estimator has been called.
double T_Time_basopt
Total Time spent in basopt.
double T_Time_inest
Time spent in initial estimator.
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
double T_Time_vcs
Time spent in the vcs suite of programs.
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
VCS_SOLVE(MultiPhase *mphase, int printLvl=0)
Initialize the sizes within the VCS_SOLVE object.
Definition vcs_solve.cpp:44
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
int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
Evaluate the standard state free energies at the current temperature and pressure.
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
int vcs_inest_TP()
Create an initial estimate of the solution to the thermodynamic equilibrium problem.
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
double vcs_VolTotal(const double tkelvin, const double pres, const double w[], double volPM[])
Calculation of the total volume and the partial molar volumes.
int vcs_elem_rearrange(double *const aw, double *const sa, double *const sm, double *const ss)
Rearrange the constraint equations represented by the Formula Matrix so that the operational ones are...
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
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
double m_totalVol
Total volume of all phases. Units are m^3.
Definition vcs_solve.h:1472
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.
int vcs_prep(int printLvl)
This routine is mostly concerned with changing the private data to be consistent with what's needed f...
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
vector< size_t > m_speciesMapIndex
Index vector that keeps track of the species vector rearrangement.
Definition vcs_solve.h:1307
void vcs_TCounters_report(int timing_print_lvl=1)
Create a report on the plog file containing timing and its information.
Array2D m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise.
Definition vcs_solve.h:1177
static void disableTiming()
Disable printing of timing information.
vector< int > m_actConventionSpecies
specifies the activity convention of the phase containing the species
Definition vcs_solve.h:1404
double vcs_phaseStabilityTest(const size_t iph)
Main program to test whether a deleted phase should be brought back into existence.
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
void vcs_inest(double *const aw, double *const sa, double *const sm, double *const ss, double test)
Estimate equilibrium compositions.
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
bool vcs_popPhasePossible(const size_t iphasePop) const
Utility function that evaluates whether a phase can be popped into existence.
int vcs_TP(int ipr, int ip1, int maxit, double T, double pres)
Solve an equilibrium problem at a particular fixed temperature and pressure.
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.
void prob_report(int print_lvl)
Print out the problem specification in all generality as it currently exists in the VCS_SOLVE object.
int m_timing_print_lvl
printing level of timing information
Definition vcs_solve.h:1506
int vcs_setMolesLinProg()
Estimate the initial mole numbers by constrained linear programming.
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.
void vcs_CalcLnActCoeffJac(const double *const moleSpeciesVCS)
Recalculate all of the activity coefficients in all of the phases based on input mole numbers.
vector< size_t > m_elementMapIndex
Index vector that keeps track of the rearrangement of the elements.
Definition vcs_solve.h:1334
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Computes the elemental abundances vector for a single phase, elemAbundPhase[], and returns it through...
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition vcs_solve.h:1058
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
void addPhaseElements(vcs_VolPhase *volPhase)
Add elements to the local element list.
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.
size_t addOnePhaseSpecies(vcs_VolPhase *volPhase, size_t k, size_t kT)
This routines adds entries for the formula matrix for one species.
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
int vcs(int ipr, int ip1, int maxit)
Solve an equilibrium problem.
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
double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
double vcs_GibbsPhase(size_t iphase, const double *const w, const double *const fe)
Calculate the total dimensionless Gibbs free energy of a single phase.
size_t m_numComponents
Number of components calculated for the problem.
Definition vcs_solve.h:1051
void vcs_prob_specifyFully()
Fully specify the problem to be solved.
double m_pressurePA
Pressure.
Definition vcs_solve.h:1274
void vcs_SSPhase()
Calculate the status of single species phases.
int vcs_report(int iconv)
Print out a report on the state of the equilibrium problem to standard output.
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
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
void vcs_delete_memory()
Delete memory that isn't just resizable STL containers.
vector< double > m_phasePhi
electric potential of the iph phase
Definition vcs_solve.h:1180
double vcs_Hessian_actCoeff_diag(size_t irxn)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
vector< double > m_elemAbundancesGoal
Element abundances vector Goals.
Definition vcs_solve.h:1232
void vcs_prob_update()
Transfer the results of the equilibrium calculation back from VCS_SOLVE.
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_phaseActConvention
specifies the activity convention of the phase.
Definition vcs_solve.h:1413
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
vector< double > m_TmpPhase2
Temporary vector of length NPhase.
Definition vcs_solve.h:1262
void vcs_switch_elem_pos(size_t ipos, size_t jpos)
Swaps the indices for all of the global data for two elements, ipos and jpos.
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
int m_printLvl
Print level for print routines.
Definition vcs_solve.h:985
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
void vcs_fePrep_TP()
Initialize the chemical potential of single species phases.
double m_temperature
Temperature (Kelvin)
Definition vcs_solve.h:1271
vector< double > m_deltaGRxn_old
Last deltag[irxn] from the previous step.
Definition vcs_solve.h:1197
size_t addElement(const char *elNameNew, int elType, int elactive)
This routine resizes the number of elements in the VCS_SOLVE object by adding a new element to the en...
double m_totalMolNum
Total number of kmoles in all phases.
Definition vcs_solve.h:1239
Identifies the thermo model for the species.
size_t IndexSpeciesPhase
Index of this species in the current phase.
int SSStar_Vol_Model
Models for the standard state volume of each species.
double SS0_Cp0
Base heat capacity used in the VCS_SS0_CONSTANT_CP model.
double SSStar_Vol0
parameter that is used in the VCS_SSVOL_CONSTANT model.
double SS0_H0
Base enthalpy used in the VCS_SS0_CONSTANT_CP model.
double SS0_T0
Base temperature used in the VCS_SS0_CONSTANT_CP model.
double SS0_TSave
Internal storage of the last temperature used in the calculation of the reference naught Gibbs free e...
int SSStar_Model
Integer value representing the star state model.
double SS0_S0
Base entropy used in the VCS_SS0_CONSTANT_CP model.
int SS0_Model
Integer representing the models for the species standard state Naught temperature dependence.
double SS0_feSave
Internal storage of the last calculation of the reference naught Gibbs free energy at SS0_TSave.
size_t IndexPhase
Index of the phase that this species belongs to.
vcs_VolPhase * OwningPhase
Pointer to the owning phase object.
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
Properties of a single species.
double VolPM
Partial molar volume of the species.
int SurfaceSpecies
True if this species belongs to a surface phase.
string SpName
Name of the species.
double WtSpecies
Molecular Weight of the species (gm/mol)
VCS_SPECIES_THERMO * SpeciesThermo
Pointer to the thermo structure for this species.
vector< double > FormulaMatrixCol
Column of the formula matrix, comprising the element composition of the species.
double Charge
Charge state of the species -> This may be duplication of what's in the FormulaMatrixCol entries.
Phase information and Phase calculations for vcs.
void setCreationMoleNumbers(const double *const n_k, const vector< size_t > &creationGlobalRxnNumbers)
Sets the creationMoleNum's within the phase object.
void setElectricPotential(const double phi)
set the electric potential of the phase
int speciesUnknownType(const size_t k) const
Returns the type of the species unknown.
size_t elemGlobalIndex(const size_t e) const
Returns the global index of the local element index for the phase.
size_t nSpecies() const
Return the number of species in the phase.
string PhaseName
String name for the phase.
void setMoleFractionsState(const double molNum, const double *const moleFracVec, const int vcsStateStatus)
Set the moles and/or mole fractions within the phase.
size_t VP_ID_
Original ID of the phase in the problem.
const Array2D & getFormulaMatrix() const
Get a constant form of the Species Formula Matrix.
vcs_SpeciesProperties * speciesProperty(const size_t kindex)
Retrieve the kth Species structure for the species belonging to this phase.
size_t phiVarIndex() const
Return the index of the species that represents the the voltage of the phase.
double totalMolesInert() const
Returns the value of the total kmol of inert in the phase.
bool m_singleSpecies
If true, this phase consists of a single species.
void setMolesFromVCS(const int stateCalc, const double *molesSpeciesVCS=0)
Set the moles within the phase.
void setElemGlobalIndex(const size_t eLocal, const size_t eGlobal)
sets a local phase element to a global index value
void resize(const size_t phaseNum, const size_t numSpecies, const size_t numElem, const char *const phaseName, const double molesInert=0.0)
The resize() function fills in all of the initial information if it is not given in the constructor.
int exists() const
int indicating whether the phase exists or not
size_t nElemConstraints() const
Returns the number of element constraints.
size_t spGlobalIndexVCS(const size_t spIndex) const
Return the Global VCS index of the kth species in the phase.
size_t transferElementsFM(const ThermoPhase *const tPhase)
Transfer all of the element information from the ThermoPhase object to the vcs_VolPhase object.
double GStar_calc_one(size_t kspec) const
Gibbs free energy calculation for standard state of one species.
bool m_gasPhase
If true, this phase is a gas-phase like phase.
void setTotalMoles(const double totalMols)
Sets the total moles in the phase.
int p_activityConvention
Convention for the activity formulation.
void sendToVCS_GStar(double *const gstar) const
Fill in the standard state Gibbs free energy vector for VCS.
void sendToVCS_ActCoeff(const int stateCalc, double *const AC)
Fill in an activity coefficients vector within a VCS_SOLVE object.
double sendToVCS_VolPM(double *const VolPM) const
Fill in the partial molar volume vector for VCS.
void sendToVCS_LnActCoeffJac(Array2D &LnACJac_VCS)
Downloads the ln ActCoeff Jacobian into the VCS version of the ln ActCoeff Jacobian.
int m_eqnState
Type of the equation of state.
double G0_calc_one(size_t kspec) const
Gibbs free energy calculation at a temperature for the reference state of a species,...
void setState_TP(const double temperature_Kelvin, const double pressure_PA)
Sets the temperature and pressure in this object and underlying ThermoPhase objects.
string eos_name() const
Return the name corresponding to the equation of state.
int elementType(const size_t e) const
Type of the element constraint with index e.
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.
const vector< double > & creationMoleNumbers(vector< size_t > &creationGlobalRxnNumbers) const
Return a const reference to the creationMoleNumbers stored in the object.
string elementName(const size_t e) const
Name of the element constraint with index e.
bool isIdealSoln() const
Returns whether the phase is an ideal solution phase.
double totalMoles() const
Return the total moles in the phase.
void setPtrThermoPhase(ThermoPhase *tp_ptr)
Set the pointer for Cantera's ThermoPhase parameter.
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 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
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition utilities.h:104
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition global.h:329
const double Faraday
Faraday constant [C/kmol].
Definition ct_defs.h:131
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
int vcs_timing_print_lvl
Global hook for turning on and off time printing.
Definition vcs_solve.cpp:42
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
double vcs_l2norm(const vector< double > &vec)
determine the l2 norm of a vector of doubles
Definition vcs_util.cpp:19
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Contains const definitions for types of species reference-state thermodynamics managers (see Species ...
#define SIMPLE
Constant Cp thermo.
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_SSVOL_IDEALGAS
Models for the standard state volume of each species.
Definition vcs_defs.h:29
#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_ZEROEDPHASE
Species lies in a multicomponent phase that is zeroed atm.
Definition vcs_defs.h:155
#define VCS_ELEM_TYPE_LATTICERATIO
Constraint associated with maintaining a fixed lattice stoichiometry in the solids.
Definition vcs_defs.h:246
#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_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
Definition vcs_defs.h:232
#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_SMALL_MULTIPHASE_SPECIES
Relative value of multiphase species mole number for a multiphase species which is small.
Definition vcs_defs.h:52
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
Definition vcs_defs.h:238
#define VCS_SPECIES_MINOR
Species is a major species.
Definition vcs_defs.h:111
#define VCS_DELETE_ELEMENTABS_CUTOFF
Cutoff moles below which a phase or species which comprises the bulk of an element's total concentrat...
Definition vcs_defs.h:72
#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_STATECALC_PHASESTABILITY
State Calculation based on tentative mole numbers for a phase which is currently zeroed,...
Definition vcs_defs.h:305
#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...