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