1 /**
2 * @file vcs_VolPhase.cpp
3 */
4
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7
8 #include "cantera/equil/vcs_VolPhase.h"
9 #include "cantera/equil/vcs_species_thermo.h"
10 #include "cantera/equil/vcs_solve.h"
11
12 #include "cantera/thermo/ThermoPhase.h"
13 #include "cantera/base/stringUtils.h"
14
15 #include <cstdio>
16
17 namespace Cantera
18 {
19
vcs_VolPhase(VCS_SOLVE * owningSolverObject)20 vcs_VolPhase::vcs_VolPhase(VCS_SOLVE* owningSolverObject) :
21 m_owningSolverObject(0),
22 VP_ID_(npos),
23 m_singleSpecies(true),
24 m_gasPhase(false),
25 m_eqnState(VCS_EOS_CONSTANT),
26 ChargeNeutralityElement(npos),
27 p_activityConvention(0),
28 m_numElemConstraints(0),
29 m_elemGlobalIndex(0),
30 m_numSpecies(0),
31 m_totalMolesInert(0.0),
32 m_isIdealSoln(false),
33 m_existence(VCS_PHASE_EXIST_NO),
34 m_MFStartIndex(0),
35 IndSpecies(0),
36 TP_ptr(0),
37 v_totalMoles(0.0),
38 m_phiVarIndex(npos),
39 m_totalVol(0.0),
40 m_vcsStateStatus(VCS_STATECALC_OLD),
41 m_phi(0.0),
42 m_UpToDate(false),
43 m_UpToDate_AC(false),
44 m_UpToDate_VolStar(false),
45 m_UpToDate_VolPM(false),
46 m_UpToDate_GStar(false),
47 m_UpToDate_G0(false),
48 Temp_(273.15),
49 Pres_(1.01325E5)
50 {
51 m_owningSolverObject = owningSolverObject;
52 }
53
~vcs_VolPhase()54 vcs_VolPhase::~vcs_VolPhase()
55 {
56 for (size_t k = 0; k < m_numSpecies; k++) {
57 delete ListSpeciesPtr[k];
58 }
59 }
60
resize(const size_t phaseNum,const size_t nspecies,const size_t numElem,const char * const phaseName,const double molesInert)61 void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies,
62 const size_t numElem, const char* const phaseName,
63 const double molesInert)
64 {
65 AssertThrowMsg(nspecies > 0, "vcs_VolPhase::resize", "nspecies Error");
66 setTotalMolesInert(molesInert);
67 m_phi = 0.0;
68 m_phiVarIndex = npos;
69
70 if (phaseNum == VP_ID_) {
71 if (strcmp(PhaseName.c_str(), phaseName)) {
72 throw CanteraError("vcs_VolPhase::resize",
73 "Strings are different: " + PhaseName + " " +
74 phaseName + " :unknown situation");
75 }
76 } else {
77 VP_ID_ = phaseNum;
78 if (!phaseName) {
79 PhaseName = fmt::format("Phase_{}", VP_ID_);
80 } else {
81 PhaseName = phaseName;
82 }
83 }
84 if (nspecies > 1) {
85 m_singleSpecies = false;
86 } else {
87 m_singleSpecies = true;
88 }
89
90 if (m_numSpecies == nspecies && numElem == m_numElemConstraints) {
91 return;
92 }
93
94 m_numSpecies = nspecies;
95 if (nspecies > 1) {
96 m_singleSpecies = false;
97 }
98
99 IndSpecies.resize(nspecies, npos);
100
101 if (ListSpeciesPtr.size() >= m_numSpecies) {
102 for (size_t i = 0; i < m_numSpecies; i++) {
103 if (ListSpeciesPtr[i]) {
104 delete ListSpeciesPtr[i];
105 ListSpeciesPtr[i] = 0;
106 }
107 }
108 }
109 ListSpeciesPtr.resize(nspecies, 0);
110 for (size_t i = 0; i < nspecies; i++) {
111 ListSpeciesPtr[i] = new vcs_SpeciesProperties(phaseNum, i, this);
112 }
113
114 Xmol_.resize(nspecies, 0.0);
115 creationMoleNumbers_.resize(nspecies, 0.0);
116 creationGlobalRxnNumbers_.resize(nspecies, npos);
117 for (size_t i = 0; i < nspecies; i++) {
118 Xmol_[i] = 1.0/nspecies;
119 creationMoleNumbers_[i] = 1.0/nspecies;
120 if (IndSpecies[i] >= m_numElemConstraints) {
121 creationGlobalRxnNumbers_[i] = IndSpecies[i] - m_numElemConstraints;
122 } else {
123 creationGlobalRxnNumbers_[i] = npos;
124 }
125 }
126
127 SS0ChemicalPotential.resize(nspecies, -1.0);
128 StarChemicalPotential.resize(nspecies, -1.0);
129 StarMolarVol.resize(nspecies, -1.0);
130 PartialMolarVol.resize(nspecies, -1.0);
131 ActCoeff.resize(nspecies, 1.0);
132 np_dLnActCoeffdMolNumber.resize(nspecies, nspecies, 0.0);
133
134 m_speciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM);
135 m_UpToDate = false;
136 m_vcsStateStatus = VCS_STATECALC_OLD;
137 m_UpToDate_AC = false;
138 m_UpToDate_VolStar = false;
139 m_UpToDate_VolPM = false;
140 m_UpToDate_GStar = false;
141 m_UpToDate_G0 = false;
142
143 elemResize(numElem);
144 }
145
elemResize(const size_t numElemConstraints)146 void vcs_VolPhase::elemResize(const size_t numElemConstraints)
147 {
148 m_elementNames.resize(numElemConstraints);
149 m_elementActive.resize(numElemConstraints+1, 1);
150 m_elementType.resize(numElemConstraints, VCS_ELEM_TYPE_ABSPOS);
151 m_formulaMatrix.resize(m_numSpecies, numElemConstraints, 0.0);
152 m_elementNames.resize(numElemConstraints, "");
153 m_elemGlobalIndex.resize(numElemConstraints, npos);
154 m_numElemConstraints = numElemConstraints;
155 }
156
_updateActCoeff() const157 void vcs_VolPhase::_updateActCoeff() const
158 {
159 if (m_isIdealSoln) {
160 m_UpToDate_AC = true;
161 return;
162 }
163 TP_ptr->getActivityCoefficients(&ActCoeff[0]);
164 m_UpToDate_AC = true;
165 }
166
_updateG0() const167 void vcs_VolPhase::_updateG0() const
168 {
169 TP_ptr->getGibbs_ref(&SS0ChemicalPotential[0]);
170 m_UpToDate_G0 = true;
171 }
172
G0_calc_one(size_t kspec) const173 double vcs_VolPhase::G0_calc_one(size_t kspec) const
174 {
175 if (!m_UpToDate_G0) {
176 _updateG0();
177 }
178 return SS0ChemicalPotential[kspec];
179 }
180
_updateGStar() const181 void vcs_VolPhase::_updateGStar() const
182 {
183 TP_ptr->getStandardChemPotentials(&StarChemicalPotential[0]);
184 m_UpToDate_GStar = true;
185 }
186
GStar_calc_one(size_t kspec) const187 double vcs_VolPhase::GStar_calc_one(size_t kspec) const
188 {
189 if (!m_UpToDate_GStar) {
190 _updateGStar();
191 }
192 return StarChemicalPotential[kspec];
193 }
194
setMoleFractions(const double * const xmol)195 void vcs_VolPhase::setMoleFractions(const double* const xmol)
196 {
197 double sum = -1.0;
198 for (size_t k = 0; k < m_numSpecies; k++) {
199 Xmol_[k] = xmol[k];
200 sum+= xmol[k];
201 }
202 if (std::fabs(sum) > 1.0E-13) {
203 for (size_t k = 0; k < m_numSpecies; k++) {
204 Xmol_[k] /= sum;
205 }
206 }
207 _updateMoleFractionDependencies();
208 m_UpToDate = false;
209 m_vcsStateStatus = VCS_STATECALC_TMP;
210 }
211
_updateMoleFractionDependencies()212 void vcs_VolPhase::_updateMoleFractionDependencies()
213 {
214 if (TP_ptr) {
215 TP_ptr->setState_PX(Pres_, &Xmol_[m_MFStartIndex]);
216 }
217 if (!m_isIdealSoln) {
218 m_UpToDate_AC = false;
219 m_UpToDate_VolPM = false;
220 }
221 }
222
moleFractions() const223 const vector_fp & vcs_VolPhase::moleFractions() const
224 {
225 return Xmol_;
226 }
227
moleFraction(size_t k) const228 double vcs_VolPhase::moleFraction(size_t k) const
229 {
230 return Xmol_[k];
231 }
232
setMoleFractionsState(const double totalMoles,const double * const moleFractions,const int vcsStateStatus)233 void vcs_VolPhase::setMoleFractionsState(const double totalMoles,
234 const double* const moleFractions,
235 const int vcsStateStatus)
236 {
237 if (totalMoles != 0.0) {
238 // There are other ways to set the mole fractions when VCS_STATECALC
239 // is set to a normal settting.
240 if (vcsStateStatus != VCS_STATECALC_TMP) {
241 throw CanteraError("vcs_VolPhase::setMolesFractionsState",
242 "inappropriate usage");
243 }
244 m_UpToDate = false;
245 m_vcsStateStatus = VCS_STATECALC_TMP;
246 if (m_existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
247 throw CanteraError("vcs_VolPhase::setMolesFractionsState",
248 "inappropriate usage");
249 }
250 m_existence = VCS_PHASE_EXIST_YES;
251 } else {
252 m_UpToDate = true;
253 m_vcsStateStatus = vcsStateStatus;
254 m_existence = std::min(m_existence, VCS_PHASE_EXIST_NO);
255 }
256 double fractotal = 1.0;
257 v_totalMoles = totalMoles;
258 if (m_totalMolesInert > 0.0) {
259 if (m_totalMolesInert > v_totalMoles) {
260 throw CanteraError("vcs_VolPhase::setMolesFractionsState",
261 "inerts greater than total: {} {}",
262 v_totalMoles, m_totalMolesInert);
263 }
264 fractotal = 1.0 - m_totalMolesInert/v_totalMoles;
265 }
266 double sum = 0.0;
267 for (size_t k = 0; k < m_numSpecies; k++) {
268 Xmol_[k] = moleFractions[k];
269 sum += moleFractions[k];
270 }
271 if (sum == 0.0) {
272 throw CanteraError("vcs_VolPhase::setMolesFractionsState",
273 "inappropriate usage");
274 }
275 if (sum != fractotal) {
276 for (size_t k = 0; k < m_numSpecies; k++) {
277 Xmol_[k] *= (fractotal /sum);
278 }
279 }
280 _updateMoleFractionDependencies();
281 }
282
setMolesFromVCS(const int stateCalc,const double * molesSpeciesVCS)283 void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
284 const double* molesSpeciesVCS)
285 {
286 v_totalMoles = m_totalMolesInert;
287
288 if (molesSpeciesVCS == 0) {
289 AssertThrowMsg(m_owningSolverObject, "vcs_VolPhase::setMolesFromVCS",
290 "shouldn't be here");
291 if (stateCalc == VCS_STATECALC_OLD) {
292 molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_old[0];
293 } else if (stateCalc == VCS_STATECALC_NEW) {
294 molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_new[0];
295 } else {
296 throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
297 }
298 } else if (m_owningSolverObject) {
299 // if (stateCalc == VCS_STATECALC_OLD) {
300 // if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_old[0]) {
301 // throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
302 // }
303 // } else if (stateCalc == VCS_STATECALC_NEW) {
304 // if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_new[0]) {
305 // throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
306 // }
307 // }
308 }
309
310 for (size_t k = 0; k < m_numSpecies; k++) {
311 if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
312 size_t kglob = IndSpecies[k];
313 v_totalMoles += std::max(0.0, molesSpeciesVCS[kglob]);
314 }
315 }
316 if (v_totalMoles > 0.0) {
317 for (size_t k = 0; k < m_numSpecies; k++) {
318 if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
319 size_t kglob = IndSpecies[k];
320 double tmp = std::max(0.0, molesSpeciesVCS[kglob]);
321 Xmol_[k] = tmp / v_totalMoles;
322 }
323 }
324 m_existence = VCS_PHASE_EXIST_YES;
325 } else {
326 // This is where we will start to store a better approximation
327 // for the mole fractions, when the phase doesn't exist.
328 // This is currently unimplemented.
329 m_existence = VCS_PHASE_EXIST_NO;
330 }
331
332 // Update the electric potential if it is a solution variable in the
333 // equation system
334 if (m_phiVarIndex != npos) {
335 size_t kglob = IndSpecies[m_phiVarIndex];
336 if (m_numSpecies == 1) {
337 Xmol_[m_phiVarIndex] = 1.0;
338 } else {
339 Xmol_[m_phiVarIndex] = 0.0;
340 }
341 double phi = molesSpeciesVCS[kglob];
342 setElectricPotential(phi);
343 if (m_numSpecies == 1) {
344 m_existence = VCS_PHASE_EXIST_YES;
345 }
346 }
347 _updateMoleFractionDependencies();
348 if (m_totalMolesInert > 0.0) {
349 m_existence = VCS_PHASE_EXIST_ALWAYS;
350 }
351
352 // If stateCalc is old and the total moles is positive, then we have a valid
353 // state. If the phase went away, it would be a valid starting point for
354 // F_k's. So, save the state.
355 if (stateCalc == VCS_STATECALC_OLD && v_totalMoles > 0.0) {
356 creationMoleNumbers_ = Xmol_;
357 }
358
359 // Set flags indicating we are up to date with the VCS state vector.
360 m_UpToDate = true;
361 m_vcsStateStatus = stateCalc;
362 }
363
setMolesFromVCSCheck(const int vcsStateStatus,const double * molesSpeciesVCS,const double * const TPhMoles)364 void vcs_VolPhase::setMolesFromVCSCheck(const int vcsStateStatus,
365 const double* molesSpeciesVCS,
366 const double* const TPhMoles)
367 {
368 setMolesFromVCS(vcsStateStatus, molesSpeciesVCS);
369
370 // Check for consistency with TPhMoles[]
371 double Tcheck = TPhMoles[VP_ID_];
372 if (Tcheck != v_totalMoles) {
373 if (vcs_doubleEqual(Tcheck, v_totalMoles)) {
374 Tcheck = v_totalMoles;
375 } else {
376 throw CanteraError("vcs_VolPhase::setMolesFromVCSCheck",
377 "We have a consistency problem: {} {}", Tcheck, v_totalMoles);
378 }
379 }
380 }
381
updateFromVCS_MoleNumbers(const int vcsStateStatus)382 void vcs_VolPhase::updateFromVCS_MoleNumbers(const int vcsStateStatus)
383 {
384 if ((!m_UpToDate || vcsStateStatus != m_vcsStateStatus) && m_owningSolverObject &&
385 (vcsStateStatus == VCS_STATECALC_OLD || vcsStateStatus == VCS_STATECALC_NEW)) {
386 setMolesFromVCS(vcsStateStatus);
387 }
388 }
389
sendToVCS_ActCoeff(const int vcsStateStatus,double * const AC)390 void vcs_VolPhase::sendToVCS_ActCoeff(const int vcsStateStatus,
391 double* const AC)
392 {
393 updateFromVCS_MoleNumbers(vcsStateStatus);
394 if (!m_UpToDate_AC) {
395 _updateActCoeff();
396 }
397 for (size_t k = 0; k < m_numSpecies; k++) {
398 size_t kglob = IndSpecies[k];
399 AC[kglob] = ActCoeff[k];
400 }
401 }
402
sendToVCS_VolPM(double * const VolPM) const403 double vcs_VolPhase::sendToVCS_VolPM(double* const VolPM) const
404 {
405 if (!m_UpToDate_VolPM) {
406 _updateVolPM();
407 }
408 for (size_t k = 0; k < m_numSpecies; k++) {
409 size_t kglob = IndSpecies[k];
410 VolPM[kglob] = PartialMolarVol[k];
411 }
412 return m_totalVol;
413 }
414
sendToVCS_GStar(double * const gstar) const415 void vcs_VolPhase::sendToVCS_GStar(double* const gstar) const
416 {
417 if (!m_UpToDate_GStar) {
418 _updateGStar();
419 }
420 for (size_t k = 0; k < m_numSpecies; k++) {
421 size_t kglob = IndSpecies[k];
422 gstar[kglob] = StarChemicalPotential[k];
423 }
424 }
425
setElectricPotential(const double phi)426 void vcs_VolPhase::setElectricPotential(const double phi)
427 {
428 m_phi = phi;
429 TP_ptr->setElectricPotential(m_phi);
430 // We have changed the state variable. Set uptodate flags to false
431 m_UpToDate_AC = false;
432 m_UpToDate_VolStar = false;
433 m_UpToDate_VolPM = false;
434 m_UpToDate_GStar = false;
435 }
436
electricPotential() const437 double vcs_VolPhase::electricPotential() const
438 {
439 return m_phi;
440 }
441
setState_TP(const double temp,const double pres)442 void vcs_VolPhase::setState_TP(const double temp, const double pres)
443 {
444 if (Temp_ == temp && Pres_ == pres) {
445 return;
446 }
447 TP_ptr->setElectricPotential(m_phi);
448 TP_ptr->setState_TP(temp, pres);
449 Temp_ = temp;
450 Pres_ = pres;
451 m_UpToDate_AC = false;
452 m_UpToDate_VolStar = false;
453 m_UpToDate_VolPM = false;
454 m_UpToDate_GStar = false;
455 m_UpToDate_G0 = false;
456 }
457
setState_T(const double temp)458 void vcs_VolPhase::setState_T(const double temp)
459 {
460 setState_TP(temp, Pres_);
461 }
462
_updateVolStar() const463 void vcs_VolPhase::_updateVolStar() const
464 {
465 TP_ptr->getStandardVolumes(&StarMolarVol[0]);
466 m_UpToDate_VolStar = true;
467 }
468
VolStar_calc_one(size_t kspec) const469 double vcs_VolPhase::VolStar_calc_one(size_t kspec) const
470 {
471 if (!m_UpToDate_VolStar) {
472 _updateVolStar();
473 }
474 return StarMolarVol[kspec];
475 }
476
_updateVolPM() const477 double vcs_VolPhase::_updateVolPM() const
478 {
479 TP_ptr->getPartialMolarVolumes(&PartialMolarVol[0]);
480 m_totalVol = 0.0;
481 for (size_t k = 0; k < m_numSpecies; k++) {
482 m_totalVol += PartialMolarVol[k] * Xmol_[k];
483 }
484 m_totalVol *= v_totalMoles;
485
486 if (m_totalMolesInert > 0.0) {
487 if (m_gasPhase) {
488 double volI = m_totalMolesInert * GasConstant * Temp_ / Pres_;
489 m_totalVol += volI;
490 } else {
491 throw CanteraError("vcs_VolPhase::_updateVolPM", "unknown situation");
492 }
493 }
494 m_UpToDate_VolPM = true;
495 return m_totalVol;
496 }
497
_updateLnActCoeffJac()498 void vcs_VolPhase::_updateLnActCoeffJac()
499 {
500 double phaseTotalMoles = v_totalMoles;
501 if (phaseTotalMoles < 1.0E-14) {
502 phaseTotalMoles = 1.0;
503 }
504
505 // Evaluate the current base activity coefficients if necessary
506 if (!m_UpToDate_AC) {
507 _updateActCoeff();
508 }
509 if (!TP_ptr) {
510 return;
511 }
512 TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &np_dLnActCoeffdMolNumber(0,0));
513 for (size_t j = 0; j < m_numSpecies; j++) {
514 double moles_j_base = phaseTotalMoles * Xmol_[j];
515 double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber.ptrColumn(j);
516 if (moles_j_base < 1.0E-200) {
517 moles_j_base = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150;
518 }
519 for (size_t k = 0; k < m_numSpecies; k++) {
520 np_lnActCoeffCol[k] = np_lnActCoeffCol[k] * phaseTotalMoles / moles_j_base;
521 }
522 }
523
524 double deltaMoles_j = 0.0;
525 // Make copies of ActCoeff and Xmol_ for use in taking differences
526 vector_fp ActCoeff_Base(ActCoeff);
527 vector_fp Xmol_Base(Xmol_);
528 double TMoles_base = phaseTotalMoles;
529
530 // Loop over the columns species to be deltad
531 for (size_t j = 0; j < m_numSpecies; j++) {
532 // Calculate a value for the delta moles of species j. Note Xmol_[] and
533 // Tmoles are always positive or zero quantities.
534 double moles_j_base = phaseTotalMoles * Xmol_Base[j];
535 deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150;
536
537 // Now, update the total moles in the phase and all of the mole
538 // fractions based on this.
539 phaseTotalMoles = TMoles_base + deltaMoles_j;
540 for (size_t k = 0; k < m_numSpecies; k++) {
541 Xmol_[k] = Xmol_Base[k] * TMoles_base / phaseTotalMoles;
542 }
543 Xmol_[j] = (moles_j_base + deltaMoles_j) / phaseTotalMoles;
544
545 // Go get new values for the activity coefficients. Note this calls
546 // setState_PX();
547 _updateMoleFractionDependencies();
548 _updateActCoeff();
549
550 // Revert to the base case Xmol_, v_totalMoles
551 v_totalMoles = TMoles_base;
552 Xmol_ = Xmol_Base;
553 }
554
555 // Go get base values for the activity coefficients. Note this calls
556 // setState_TPX() again; Just wanted to make sure that cantera is in sync
557 // with VolPhase after this call.
558 setMoleFractions(&Xmol_Base[0]);
559 _updateMoleFractionDependencies();
560 _updateActCoeff();
561 }
562
sendToVCS_LnActCoeffJac(Array2D & np_LnACJac_VCS)563 void vcs_VolPhase::sendToVCS_LnActCoeffJac(Array2D& np_LnACJac_VCS)
564 {
565 // update the Ln Act Coeff Jacobian entries with respect to the mole number
566 // of species in the phase -> we always assume that they are out of date.
567 _updateLnActCoeffJac();
568
569 // Now copy over the values
570 for (size_t j = 0; j < m_numSpecies; j++) {
571 size_t jglob = IndSpecies[j];
572 for (size_t k = 0; k < m_numSpecies; k++) {
573 size_t kglob = IndSpecies[k];
574 np_LnACJac_VCS(kglob,jglob) = np_dLnActCoeffdMolNumber(k,j);
575 }
576 }
577 }
578
setPtrThermoPhase(ThermoPhase * tp_ptr)579 void vcs_VolPhase::setPtrThermoPhase(ThermoPhase* tp_ptr)
580 {
581 TP_ptr = tp_ptr;
582 Temp_ = TP_ptr->temperature();
583 Pres_ = TP_ptr->pressure();
584 setState_TP(Temp_, Pres_);
585 m_phi = TP_ptr->electricPotential();
586 size_t nsp = TP_ptr->nSpecies();
587 size_t nelem = TP_ptr->nElements();
588 if (nsp != m_numSpecies) {
589 if (m_numSpecies != 0) {
590 warn_user("vcs_VolPhase::setPtrThermoPhase",
591 "Nsp != NVolSpeces: {} {}", nsp, m_numSpecies);
592 }
593 resize(VP_ID_, nsp, nelem, PhaseName.c_str());
594 }
595 TP_ptr->getMoleFractions(&Xmol_[0]);
596 creationMoleNumbers_ = Xmol_;
597 _updateMoleFractionDependencies();
598
599 // figure out ideal solution tag
600 if (nsp == 1) {
601 m_isIdealSoln = true;
602 } else {
603 std::string eos = TP_ptr->type();
604 if (eos == "IdealGas" || eos == "ConstDensity" || eos == "Surf"
605 || eos == "Metal" || eos == "StoichSubstance"
606 || eos == "LatticeSolid"
607 || eos == "Lattice" || eos == "Edge" || eos == "IdealSolidSoln") {
608 m_isIdealSoln = true;
609 } else {
610 m_isIdealSoln = false;
611 };
612 }
613 }
614
ptrThermoPhase() const615 const ThermoPhase* vcs_VolPhase::ptrThermoPhase() const
616 {
617 return TP_ptr;
618 }
619
totalMoles() const620 double vcs_VolPhase::totalMoles() const
621 {
622 return v_totalMoles;
623 }
624
molefraction(size_t k) const625 double vcs_VolPhase::molefraction(size_t k) const
626 {
627 return Xmol_[k];
628 }
629
setCreationMoleNumbers(const double * const n_k,const std::vector<size_t> & creationGlobalRxnNumbers)630 void vcs_VolPhase::setCreationMoleNumbers(const double* const n_k,
631 const std::vector<size_t> &creationGlobalRxnNumbers)
632 {
633 creationMoleNumbers_.assign(n_k, n_k+m_numSpecies);
634 for (size_t k = 0; k < m_numSpecies; k++) {
635 creationGlobalRxnNumbers_[k] = creationGlobalRxnNumbers[k];
636 }
637 }
638
creationMoleNumbers(std::vector<size_t> & creationGlobalRxnNumbers) const639 const vector_fp& vcs_VolPhase::creationMoleNumbers(std::vector<size_t> &creationGlobalRxnNumbers) const
640 {
641 creationGlobalRxnNumbers = creationGlobalRxnNumbers_;
642 return creationMoleNumbers_;
643 }
644
setTotalMoles(const double totalMols)645 void vcs_VolPhase::setTotalMoles(const double totalMols)
646 {
647 v_totalMoles = totalMols;
648 if (m_totalMolesInert > 0.0) {
649 m_existence = VCS_PHASE_EXIST_ALWAYS;
650 AssertThrowMsg(totalMols >= m_totalMolesInert,
651 "vcs_VolPhase::setTotalMoles",
652 "totalMoles less than inert moles: {} {}",
653 totalMols, m_totalMolesInert);
654 } else {
655 if (m_singleSpecies && (m_phiVarIndex == 0)) {
656 m_existence = VCS_PHASE_EXIST_ALWAYS;
657 } else {
658 if (totalMols > 0.0) {
659 m_existence = VCS_PHASE_EXIST_YES;
660 } else {
661 m_existence = VCS_PHASE_EXIST_NO;
662 }
663 }
664 }
665 }
666
setMolesOutOfDate(int stateCalc)667 void vcs_VolPhase::setMolesOutOfDate(int stateCalc)
668 {
669 m_UpToDate = false;
670 if (stateCalc != -1) {
671 m_vcsStateStatus = stateCalc;
672 }
673 }
674
setMolesCurrent(int stateCalc)675 void vcs_VolPhase::setMolesCurrent(int stateCalc)
676 {
677 m_UpToDate = true;
678 m_vcsStateStatus = stateCalc;
679 }
680
isIdealSoln() const681 bool vcs_VolPhase::isIdealSoln() const
682 {
683 return m_isIdealSoln;
684 }
685
phiVarIndex() const686 size_t vcs_VolPhase::phiVarIndex() const
687 {
688 return m_phiVarIndex;
689 }
690
setPhiVarIndex(size_t phiVarIndex)691 void vcs_VolPhase::setPhiVarIndex(size_t phiVarIndex)
692 {
693 m_phiVarIndex = phiVarIndex;
694 m_speciesUnknownType[m_phiVarIndex] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE;
695 if (m_singleSpecies && m_phiVarIndex == 0) {
696 m_existence = VCS_PHASE_EXIST_ALWAYS;
697 }
698 }
699
speciesProperty(const size_t kindex)700 vcs_SpeciesProperties* vcs_VolPhase::speciesProperty(const size_t kindex)
701 {
702 return ListSpeciesPtr[kindex];
703 }
704
exists() const705 int vcs_VolPhase::exists() const
706 {
707 return m_existence;
708 }
709
setExistence(const int existence)710 void vcs_VolPhase::setExistence(const int existence)
711 {
712 if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
713 if (v_totalMoles != 0.0) {
714 throw CanteraError("vcs_VolPhase::setExistence",
715 "setting false existence for phase with moles");
716 }
717 } else if (m_totalMolesInert == 0.0) {
718 if (v_totalMoles == 0.0 && (!m_singleSpecies || m_phiVarIndex != 0)) {
719 throw CanteraError("vcs_VolPhase::setExistence",
720 "setting true existence for phase with no moles");
721 }
722 }
723 if (m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) {
724 throw CanteraError("vcs_VolPhase::setExistence",
725 "Trying to set existence of an electron phase to false");
726 }
727 m_existence = existence;
728 }
729
spGlobalIndexVCS(const size_t spIndex) const730 size_t vcs_VolPhase::spGlobalIndexVCS(const size_t spIndex) const
731 {
732 return IndSpecies[spIndex];
733 }
734
setSpGlobalIndexVCS(const size_t spIndex,const size_t spGlobalIndex)735 void vcs_VolPhase::setSpGlobalIndexVCS(const size_t spIndex,
736 const size_t spGlobalIndex)
737 {
738 IndSpecies[spIndex] = spGlobalIndex;
739 if (spGlobalIndex >= m_numElemConstraints) {
740 creationGlobalRxnNumbers_[spIndex] = spGlobalIndex - m_numElemConstraints;
741 }
742 }
743
setTotalMolesInert(const double tMolesInert)744 void vcs_VolPhase::setTotalMolesInert(const double tMolesInert)
745 {
746 if (m_totalMolesInert != tMolesInert) {
747 m_UpToDate = false;
748 m_UpToDate_AC = false;
749 m_UpToDate_VolStar = false;
750 m_UpToDate_VolPM = false;
751 m_UpToDate_GStar = false;
752 m_UpToDate_G0 = false;
753 v_totalMoles += (tMolesInert - m_totalMolesInert);
754 m_totalMolesInert = tMolesInert;
755 }
756 if (m_totalMolesInert > 0.0) {
757 m_existence = VCS_PHASE_EXIST_ALWAYS;
758 } else if (m_singleSpecies && (m_phiVarIndex == 0)) {
759 m_existence = VCS_PHASE_EXIST_ALWAYS;
760 } else {
761 if (v_totalMoles > 0.0) {
762 m_existence = VCS_PHASE_EXIST_YES;
763 } else {
764 m_existence = VCS_PHASE_EXIST_NO;
765 }
766 }
767 }
768
totalMolesInert() const769 double vcs_VolPhase::totalMolesInert() const
770 {
771 return m_totalMolesInert;
772 }
773
elemGlobalIndex(const size_t e) const774 size_t vcs_VolPhase::elemGlobalIndex(const size_t e) const
775 {
776 AssertThrow(e < m_numElemConstraints, " vcs_VolPhase::elemGlobalIndex");
777 return m_elemGlobalIndex[e];
778 }
779
setElemGlobalIndex(const size_t eLocal,const size_t eGlobal)780 void vcs_VolPhase::setElemGlobalIndex(const size_t eLocal, const size_t eGlobal)
781 {
782 AssertThrow(eLocal < m_numElemConstraints,
783 "vcs_VolPhase::setElemGlobalIndex");
784 m_elemGlobalIndex[eLocal] = eGlobal;
785 }
786
nElemConstraints() const787 size_t vcs_VolPhase::nElemConstraints() const
788 {
789 return m_numElemConstraints;
790 }
791
elementName(const size_t e) const792 std::string vcs_VolPhase::elementName(const size_t e) const
793 {
794 return m_elementNames[e];
795 }
796
797 //! This function decides whether a phase has charged species or not.
hasChargedSpecies(const ThermoPhase * const tPhase)798 static bool hasChargedSpecies(const ThermoPhase* const tPhase)
799 {
800 for (size_t k = 0; k < tPhase->nSpecies(); k++) {
801 if (tPhase->charge(k) != 0.0) {
802 return true;
803 }
804 }
805 return false;
806 }
807
808 //! This utility routine decides whether a Cantera ThermoPhase needs
809 //! a constraint equation representing the charge neutrality of the
810 //! phase. It does this by searching for charged species. If it
811 //! finds one, and if the phase needs one, then it returns true.
chargeNeutralityElement(const ThermoPhase * const tPhase)812 static bool chargeNeutralityElement(const ThermoPhase* const tPhase)
813 {
814 int hasCharge = hasChargedSpecies(tPhase);
815 if (tPhase->chargeNeutralityNecessary() && hasCharge) {
816 return true;
817 }
818 return false;
819 }
820
transferElementsFM(const ThermoPhase * const tPhase)821 size_t vcs_VolPhase::transferElementsFM(const ThermoPhase* const tPhase)
822 {
823 size_t nebase = tPhase->nElements();
824 size_t ne = nebase;
825 size_t ns = tPhase->nSpecies();
826
827 // Decide whether we need an extra element constraint for charge
828 // neutrality of the phase
829 bool cne = chargeNeutralityElement(tPhase);
830 if (cne) {
831 ChargeNeutralityElement = ne;
832 ne++;
833 }
834
835 // Assign and resize structures
836 elemResize(ne);
837
838 if (ChargeNeutralityElement != npos) {
839 m_elementType[ChargeNeutralityElement] = VCS_ELEM_TYPE_CHARGENEUTRALITY;
840 }
841
842 size_t eFound = npos;
843 if (hasChargedSpecies(tPhase)) {
844 if (cne) {
845 // We need a charge neutrality constraint. We also have an Electron
846 // Element. These are duplicates of each other. To avoid trouble
847 // with possible range error conflicts, sometimes we eliminate the
848 // Electron condition. Flag that condition for elimination by
849 // toggling the ElActive variable. If we find we need it later, we
850 // will retoggle ElActive to true.
851 for (size_t eT = 0; eT < nebase; eT++) {
852 if (tPhase->elementName(eT) == "E") {
853 eFound = eT;
854 m_elementActive[eT] = 0;
855 m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
856 }
857 }
858 } else {
859 for (size_t eT = 0; eT < nebase; eT++) {
860 if (tPhase->elementName(eT) == "E") {
861 eFound = eT;
862 m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
863 }
864 }
865 }
866 if (eFound == npos) {
867 eFound = ne;
868 m_elementType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE;
869 m_elementActive[ne] = 0;
870 std::string ename = "E";
871 m_elementNames[ne] = ename;
872 ne++;
873 elemResize(ne);
874 }
875 }
876
877 m_formulaMatrix.resize(ns, ne, 0.0);
878 m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM);
879 elemResize(ne);
880
881 size_t e = 0;
882 for (size_t eT = 0; eT < nebase; eT++) {
883 m_elementNames[e] = tPhase->elementName(eT);
884 m_elementType[e] = tPhase->elementType(eT);
885 e++;
886 }
887
888 if (cne) {
889 std::string pname = tPhase->name();
890 if (pname == "") {
891 pname = fmt::format("phase{}", VP_ID_);
892 }
893 e = ChargeNeutralityElement;
894 m_elementNames[e] = "cn_" + pname;
895 }
896
897 for (size_t k = 0; k < ns; k++) {
898 e = 0;
899 for (size_t eT = 0; eT < nebase; eT++) {
900 m_formulaMatrix(k,e) = tPhase->nAtoms(k, eT);
901 e++;
902 }
903 if (eFound != npos) {
904 m_formulaMatrix(k,eFound) = - tPhase->charge(k);
905 }
906 }
907
908 if (cne) {
909 for (size_t k = 0; k < ns; k++) {
910 m_formulaMatrix(k,ChargeNeutralityElement) = tPhase->charge(k);
911 }
912 }
913
914 // Here, we figure out what is the species types are The logic isn't set in
915 // stone, and is just for a particular type of problem that I'm solving
916 // first.
917 if (ns == 1 && tPhase->charge(0) != 0.0) {
918 m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE;
919 setPhiVarIndex(0);
920 }
921
922 return ne;
923 }
924
elementType(const size_t e) const925 int vcs_VolPhase::elementType(const size_t e) const
926 {
927 return m_elementType[e];
928 }
929
getFormulaMatrix() const930 const Array2D& vcs_VolPhase::getFormulaMatrix() const
931 {
932 return m_formulaMatrix;
933 }
934
speciesUnknownType(const size_t k) const935 int vcs_VolPhase::speciesUnknownType(const size_t k) const
936 {
937 return m_speciesUnknownType[k];
938 }
939
elementActive(const size_t e) const940 int vcs_VolPhase::elementActive(const size_t e) const
941 {
942 return m_elementActive[e];
943 }
944
nSpecies() const945 size_t vcs_VolPhase::nSpecies() const
946 {
947 return m_numSpecies;
948 }
949
eos_name() const950 std::string vcs_VolPhase::eos_name() const
951 {
952 switch (m_eqnState) {
953 case VCS_EOS_CONSTANT:
954 return "Constant";
955 case VCS_EOS_IDEAL_GAS:
956 return "Ideal Gas";
957 case VCS_EOS_STOICH_SUB:
958 return "Stoich Sub";
959 case VCS_EOS_IDEAL_SOLN:
960 return "Ideal Soln";
961 case VCS_EOS_DEBEYE_HUCKEL:
962 return "Debeye Huckel";
963 case VCS_EOS_REDLICH_KWONG:
964 return "Redlich_Kwong";
965 case VCS_EOS_REGULAR_SOLN:
966 return "Regular Soln";
967 default:
968 return fmt::format("UnkType: {:7d}", m_eqnState);
969 break;
970 }
971 }
972
973 }
974