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