Cantera  2.5.1
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 
17 namespace Cantera
18 {
19 
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 
54 vcs_VolPhase::~vcs_VolPhase()
55 {
56  for (size_t k = 0; k < m_numSpecies; k++) {
57  delete ListSpeciesPtr[k];
58  }
59 }
60 
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 
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 
157 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 
167 void vcs_VolPhase::_updateG0() const
168 {
169  TP_ptr->getGibbs_ref(&SS0ChemicalPotential[0]);
170  m_UpToDate_G0 = true;
171 }
172 
173 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 
181 void vcs_VolPhase::_updateGStar() const
182 {
183  TP_ptr->getStandardChemPotentials(&StarChemicalPotential[0]);
184  m_UpToDate_GStar = true;
185 }
186 
187 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 
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 
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 
223 const vector_fp & vcs_VolPhase::moleFractions() const
224 {
225  return Xmol_;
226 }
227 
228 double vcs_VolPhase::moleFraction(size_t k) const
229 {
230  return Xmol_[k];
231 }
232 
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 
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 
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 
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 
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 
403 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 
415 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 
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 
437 double vcs_VolPhase::electricPotential() const
438 {
439  return m_phi;
440 }
441 
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 
458 void vcs_VolPhase::setState_T(const double temp)
459 {
460  setState_TP(temp, Pres_);
461 }
462 
463 void vcs_VolPhase::_updateVolStar() const
464 {
465  TP_ptr->getStandardVolumes(&StarMolarVol[0]);
466  m_UpToDate_VolStar = true;
467 }
468 
469 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 
477 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 
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 
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 
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 
615 const ThermoPhase* vcs_VolPhase::ptrThermoPhase() const
616 {
617  return TP_ptr;
618 }
619 
620 double vcs_VolPhase::totalMoles() const
621 {
622  return v_totalMoles;
623 }
624 
625 double vcs_VolPhase::molefraction(size_t k) const
626 {
627  return Xmol_[k];
628 }
629 
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 
639 const vector_fp& vcs_VolPhase::creationMoleNumbers(std::vector<size_t> &creationGlobalRxnNumbers) const
640 {
641  creationGlobalRxnNumbers = creationGlobalRxnNumbers_;
642  return creationMoleNumbers_;
643 }
644 
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 
667 void vcs_VolPhase::setMolesOutOfDate(int stateCalc)
668 {
669  m_UpToDate = false;
670  if (stateCalc != -1) {
671  m_vcsStateStatus = stateCalc;
672  }
673 }
674 
675 void vcs_VolPhase::setMolesCurrent(int stateCalc)
676 {
677  m_UpToDate = true;
678  m_vcsStateStatus = stateCalc;
679 }
680 
681 bool vcs_VolPhase::isIdealSoln() const
682 {
683  return m_isIdealSoln;
684 }
685 
686 size_t vcs_VolPhase::phiVarIndex() const
687 {
688  return m_phiVarIndex;
689 }
690 
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 
700 vcs_SpeciesProperties* vcs_VolPhase::speciesProperty(const size_t kindex)
701 {
702  return ListSpeciesPtr[kindex];
703 }
704 
705 int vcs_VolPhase::exists() const
706 {
707  return m_existence;
708 }
709 
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 
730 size_t vcs_VolPhase::spGlobalIndexVCS(const size_t spIndex) const
731 {
732  return IndSpecies[spIndex];
733 }
734 
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 
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 
769 double vcs_VolPhase::totalMolesInert() const
770 {
771  return m_totalMolesInert;
772 }
773 
774 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 
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 
787 size_t vcs_VolPhase::nElemConstraints() const
788 {
789  return m_numElemConstraints;
790 }
791 
792 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.
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.
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 
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->id();
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 
925 int vcs_VolPhase::elementType(const size_t e) const
926 {
927  return m_elementType[e];
928 }
929 
930 const Array2D& vcs_VolPhase::getFormulaMatrix() const
931 {
932  return m_formulaMatrix;
933 }
934 
935 int vcs_VolPhase::speciesUnknownType(const size_t k) const
936 {
937  return m_speciesUnknownType[k];
938 }
939 
940 int vcs_VolPhase::elementActive(const size_t e) const
941 {
942  return m_elementActive[e];
943 }
944 
945 size_t vcs_VolPhase::nSpecies() const
946 {
947  return m_numSpecies;
948 }
949 
950 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 }
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.
Definition: ctexceptions.h:61
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:285
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:643
std::string id() const
Return the string id for the phase.
Definition: Phase.cpp:69
int elementType(size_t m) const
Return the element constraint type Possible types include:
Definition: Phase.cpp:156
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:168
doublereal temperature() const
Temperature (K).
Definition: Phase.h:667
size_t nElements() const
Number of elements.
Definition: Phase.cpp:95
std::string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:114
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:102
bool chargeNeutralityNecessary() const
Returns the chargeNeutralityNecessity boolean.
Definition: ThermoPhase.h:225
Properties of a single species.
#define AssertThrow(expr, procedure)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:248
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:265
void warn_user(const std::string &method, const std::string &msg, const Args &... args)
Definition: global.h:206
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:188
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:180
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:109
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
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:91
static bool chargeNeutralityElement(const ThermoPhase *const tPhase)
This utility routine decides whether a Cantera ThermoPhase needs a constraint equation representing t...
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:289
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:300
#define VCS_STATECALC_TMP
State Calculation based on a temporary set of mole numbers.
Definition: vcs_defs.h:311
#define VCS_PHASE_EXIST_NO
Phase doesn't currently exist in the mixture.
Definition: vcs_defs.h:198
#define VCS_PHASE_EXIST_ALWAYS
Always exists because it contains inerts which can't exist in any other phase.
Definition: vcs_defs.h:185
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:228
#define VCS_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
Definition: vcs_defs.h:303
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:281
#define VCS_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
Definition: vcs_defs.h:234
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
Definition: vcs_defs.h:240
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
Definition: vcs_defs.h:188
#define VCS_PHASE_EXIST_ZEROEDPHASE
Phase currently is zeroed due to a programmatic issue.
Definition: vcs_defs.h:204
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...