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