Cantera  2.4.0
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 http://www.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  plogf("Warning Nsp != NVolSpeces: %d %d \n", nsp, m_numSpecies);
591  }
592  resize(VP_ID_, nsp, nelem, PhaseName.c_str());
593  }
594  TP_ptr->getMoleFractions(&Xmol_[0]);
595  creationMoleNumbers_ = Xmol_;
596  _updateMoleFractionDependencies();
597 
598  // figure out ideal solution tag
599  if (nsp == 1) {
600  m_isIdealSoln = true;
601  } else {
602  std::string eos = TP_ptr->type();
603  if (eos == "IdealGas" || eos == "ConstDensity" || eos == "Surf"
604  || eos == "Metal" || eos == "StoichSubstance"
605  || eos == "Semiconductor" || eos == "LatticeSolid"
606  || eos == "Lattice" || eos == "Edge" || eos == "IdealSolidSoln") {
607  m_isIdealSoln = true;
608  } else {
609  m_isIdealSoln = false;
610  };
611  }
612 }
613 
614 const ThermoPhase* vcs_VolPhase::ptrThermoPhase() const
615 {
616  return TP_ptr;
617 }
618 
619 double vcs_VolPhase::totalMoles() const
620 {
621  return v_totalMoles;
622 }
623 
624 double vcs_VolPhase::molefraction(size_t k) const
625 {
626  return Xmol_[k];
627 }
628 
629 void vcs_VolPhase::setCreationMoleNumbers(const double* const n_k,
630  const std::vector<size_t> &creationGlobalRxnNumbers)
631 {
632  creationMoleNumbers_.assign(n_k, n_k+m_numSpecies);
633  for (size_t k = 0; k < m_numSpecies; k++) {
634  creationGlobalRxnNumbers_[k] = creationGlobalRxnNumbers[k];
635  }
636 }
637 
638 const vector_fp& vcs_VolPhase::creationMoleNumbers(std::vector<size_t> &creationGlobalRxnNumbers) const
639 {
640  creationGlobalRxnNumbers = creationGlobalRxnNumbers_;
641  return creationMoleNumbers_;
642 }
643 
644 void vcs_VolPhase::setTotalMoles(const double totalMols)
645 {
646  v_totalMoles = totalMols;
647  if (m_totalMolesInert > 0.0) {
648  m_existence = VCS_PHASE_EXIST_ALWAYS;
649  AssertThrowMsg(totalMols >= m_totalMolesInert,
650  "vcs_VolPhase::setTotalMoles",
651  "totalMoles less than inert moles: {} {}",
652  totalMols, m_totalMolesInert);
653  } else {
654  if (m_singleSpecies && (m_phiVarIndex == 0)) {
655  m_existence = VCS_PHASE_EXIST_ALWAYS;
656  } else {
657  if (totalMols > 0.0) {
658  m_existence = VCS_PHASE_EXIST_YES;
659  } else {
660  m_existence = VCS_PHASE_EXIST_NO;
661  }
662  }
663  }
664 }
665 
666 void vcs_VolPhase::setMolesOutOfDate(int stateCalc)
667 {
668  m_UpToDate = false;
669  if (stateCalc != -1) {
670  m_vcsStateStatus = stateCalc;
671  }
672 }
673 
674 void vcs_VolPhase::setMolesCurrent(int stateCalc)
675 {
676  m_UpToDate = true;
677  m_vcsStateStatus = stateCalc;
678 }
679 
680 bool vcs_VolPhase::isIdealSoln() const
681 {
682  return m_isIdealSoln;
683 }
684 
685 size_t vcs_VolPhase::phiVarIndex() const
686 {
687  return m_phiVarIndex;
688 }
689 
690 void vcs_VolPhase::setPhiVarIndex(size_t phiVarIndex)
691 {
692  m_phiVarIndex = phiVarIndex;
693  m_speciesUnknownType[m_phiVarIndex] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE;
694  if (m_singleSpecies && m_phiVarIndex == 0) {
695  m_existence = VCS_PHASE_EXIST_ALWAYS;
696  }
697 }
698 
699 vcs_SpeciesProperties* vcs_VolPhase::speciesProperty(const size_t kindex)
700 {
701  return ListSpeciesPtr[kindex];
702 }
703 
704 int vcs_VolPhase::exists() const
705 {
706  return m_existence;
707 }
708 
709 void vcs_VolPhase::setExistence(const int existence)
710 {
711  if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
712  if (v_totalMoles != 0.0) {
713  throw CanteraError("vcs_VolPhase::setExistence",
714  "setting false existence for phase with moles");
715  }
716  } else if (m_totalMolesInert == 0.0) {
717  if (v_totalMoles == 0.0 && (!m_singleSpecies || m_phiVarIndex != 0)) {
718  throw CanteraError("vcs_VolPhase::setExistence",
719  "setting true existence for phase with no moles");
720  }
721  }
722  if (m_singleSpecies && m_phiVarIndex == 0 && (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE)) {
723  throw CanteraError("vcs_VolPhase::setExistence",
724  "Trying to set existence of an electron phase to false");
725  }
726  m_existence = existence;
727 }
728 
729 size_t vcs_VolPhase::spGlobalIndexVCS(const size_t spIndex) const
730 {
731  return IndSpecies[spIndex];
732 }
733 
734 void vcs_VolPhase::setSpGlobalIndexVCS(const size_t spIndex,
735  const size_t spGlobalIndex)
736 {
737  IndSpecies[spIndex] = spGlobalIndex;
738  if (spGlobalIndex >= m_numElemConstraints) {
739  creationGlobalRxnNumbers_[spIndex] = spGlobalIndex - m_numElemConstraints;
740  }
741 }
742 
743 void vcs_VolPhase::setTotalMolesInert(const double tMolesInert)
744 {
745  if (m_totalMolesInert != tMolesInert) {
746  m_UpToDate = false;
747  m_UpToDate_AC = false;
748  m_UpToDate_VolStar = false;
749  m_UpToDate_VolPM = false;
750  m_UpToDate_GStar = false;
751  m_UpToDate_G0 = false;
752  v_totalMoles += (tMolesInert - m_totalMolesInert);
753  m_totalMolesInert = tMolesInert;
754  }
755  if (m_totalMolesInert > 0.0) {
756  m_existence = VCS_PHASE_EXIST_ALWAYS;
757  } else if (m_singleSpecies && (m_phiVarIndex == 0)) {
758  m_existence = VCS_PHASE_EXIST_ALWAYS;
759  } else {
760  if (v_totalMoles > 0.0) {
761  m_existence = VCS_PHASE_EXIST_YES;
762  } else {
763  m_existence = VCS_PHASE_EXIST_NO;
764  }
765  }
766 }
767 
768 double vcs_VolPhase::totalMolesInert() const
769 {
770  return m_totalMolesInert;
771 }
772 
773 size_t vcs_VolPhase::elemGlobalIndex(const size_t e) const
774 {
775  AssertThrow(e < m_numElemConstraints, " vcs_VolPhase::elemGlobalIndex");
776  return m_elemGlobalIndex[e];
777 }
778 
779 void vcs_VolPhase::setElemGlobalIndex(const size_t eLocal, const size_t eGlobal)
780 {
781  AssertThrow(eLocal < m_numElemConstraints,
782  "vcs_VolPhase::setElemGlobalIndex");
783  m_elemGlobalIndex[eLocal] = eGlobal;
784 }
785 
786 size_t vcs_VolPhase::nElemConstraints() const
787 {
788  return m_numElemConstraints;
789 }
790 
791 std::string vcs_VolPhase::elementName(const size_t e) const
792 {
793  return m_elementNames[e];
794 }
795 
796 //! This function decides whether a phase has charged species or not.
797 static bool hasChargedSpecies(const ThermoPhase* const tPhase)
798 {
799  for (size_t k = 0; k < tPhase->nSpecies(); k++) {
800  if (tPhase->charge(k) != 0.0) {
801  return true;
802  }
803  }
804  return false;
805 }
806 
807 //! This utility routine decides whether a Cantera ThermoPhase needs
808 //! a constraint equation representing the charge neutrality of the
809 //! phase. It does this by searching for charged species. If it
810 //! finds one, and if the phase needs one, then it returns true.
811 static bool chargeNeutralityElement(const ThermoPhase* const tPhase)
812 {
813  int hasCharge = hasChargedSpecies(tPhase);
814  if (tPhase->chargeNeutralityNecessary() && hasCharge) {
815  return true;
816  }
817  return false;
818 }
819 
820 size_t vcs_VolPhase::transferElementsFM(const ThermoPhase* const tPhase)
821 {
822  size_t nebase = tPhase->nElements();
823  size_t ne = nebase;
824  size_t ns = tPhase->nSpecies();
825 
826  // Decide whether we need an extra element constraint for charge
827  // neutrality of the phase
828  bool cne = chargeNeutralityElement(tPhase);
829  if (cne) {
830  ChargeNeutralityElement = ne;
831  ne++;
832  }
833 
834  // Assign and resize structures
835  elemResize(ne);
836 
837  if (ChargeNeutralityElement != npos) {
838  m_elementType[ChargeNeutralityElement] = VCS_ELEM_TYPE_CHARGENEUTRALITY;
839  }
840 
841  size_t eFound = npos;
842  if (hasChargedSpecies(tPhase)) {
843  if (cne) {
844  // We need a charge neutrality constraint. We also have an Electron
845  // Element. These are duplicates of each other. To avoid trouble
846  // with possible range error conflicts, sometimes we eliminate the
847  // Electron condition. Flag that condition for elimination by
848  // toggling the ElActive variable. If we find we need it later, we
849  // will retoggle ElActive to true.
850  for (size_t eT = 0; eT < nebase; eT++) {
851  if (tPhase->elementName(eT) == "E") {
852  eFound = eT;
853  m_elementActive[eT] = 0;
854  m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
855  }
856  }
857  } else {
858  for (size_t eT = 0; eT < nebase; eT++) {
859  if (tPhase->elementName(eT) == "E") {
860  eFound = eT;
861  m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
862  }
863  }
864  }
865  if (eFound == npos) {
866  eFound = ne;
867  m_elementType[ne] = VCS_ELEM_TYPE_ELECTRONCHARGE;
868  m_elementActive[ne] = 0;
869  std::string ename = "E";
870  m_elementNames[ne] = ename;
871  ne++;
872  elemResize(ne);
873  }
874  }
875 
876  m_formulaMatrix.resize(ns, ne, 0.0);
877  m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM);
878  elemResize(ne);
879 
880  size_t e = 0;
881  for (size_t eT = 0; eT < nebase; eT++) {
882  m_elementNames[e] = tPhase->elementName(eT);
883  m_elementType[e] = tPhase->elementType(eT);
884  e++;
885  }
886 
887  if (cne) {
888  std::string pname = tPhase->id();
889  if (pname == "") {
890  pname = fmt::format("phase{}", VP_ID_);
891  }
892  e = ChargeNeutralityElement;
893  m_elementNames[e] = "cn_" + pname;
894  }
895 
896  for (size_t k = 0; k < ns; k++) {
897  e = 0;
898  for (size_t eT = 0; eT < nebase; eT++) {
899  m_formulaMatrix(k,e) = tPhase->nAtoms(k, eT);
900  e++;
901  }
902  if (eFound != npos) {
903  m_formulaMatrix(k,eFound) = - tPhase->charge(k);
904  }
905  }
906 
907  if (cne) {
908  for (size_t k = 0; k < ns; k++) {
909  m_formulaMatrix(k,ChargeNeutralityElement) = tPhase->charge(k);
910  }
911  }
912 
913  // Here, we figure out what is the species types are The logic isn't set in
914  // stone, and is just for a particular type of problem that I'm solving
915  // first.
916  if (ns == 1 && tPhase->charge(0) != 0.0) {
917  m_speciesUnknownType[0] = VCS_SPECIES_TYPE_INTERFACIALVOLTAGE;
918  setPhiVarIndex(0);
919  }
920 
921  return ne;
922 }
923 
924 int vcs_VolPhase::elementType(const size_t e) const
925 {
926  return m_elementType[e];
927 }
928 
929 const Array2D& vcs_VolPhase::getFormulaMatrix() const
930 {
931  return m_formulaMatrix;
932 }
933 
934 int vcs_VolPhase::speciesUnknownType(const size_t k) const
935 {
936  return m_speciesUnknownType[k];
937 }
938 
939 int vcs_VolPhase::elementActive(const size_t e) const
940 {
941  return m_elementActive[e];
942 }
943 
944 size_t vcs_VolPhase::nSpecies() const
945 {
946  return m_numSpecies;
947 }
948 
949 std::string vcs_VolPhase::eos_name() const
950 {
951  switch (m_eqnState) {
952  case VCS_EOS_CONSTANT:
953  return "Constant";
954  case VCS_EOS_IDEAL_GAS:
955  return "Ideal Gas";
956  case VCS_EOS_STOICH_SUB:
957  return "Stoich Sub";
958  case VCS_EOS_IDEAL_SOLN:
959  return "Ideal Soln";
960  case VCS_EOS_DEBEYE_HUCKEL:
961  return "Debeye Huckel";
962  case VCS_EOS_REDLICH_KWONG:
963  return "Redlich_Kwong";
964  case VCS_EOS_REGULAR_SOLN:
965  return "Regular Soln";
966  default:
967  return fmt::format("UnkType: {:7d}", m_eqnState);
968  break;
969  }
970 }
971 
972 }
#define VCS_ELEM_TYPE_ELECTRONCHARGE
This refers to conservation of electrons.
Definition: vcs_defs.h:234
#define VCS_PHASE_EXIST_NO
Phase doesn&#39;t currently exist in the mixture.
Definition: vcs_defs.h:198
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
size_t nElements() const
Number of elements.
Definition: Phase.cpp:88
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
#define VCS_PHASE_EXIST_ALWAYS
Always exists because it contains inerts which can&#39;t exist in any other phase.
Definition: vcs_defs.h:185
#define VCS_PHASE_EXIST_ZEROEDPHASE
Phase currently is zeroed due to a programmatic issue.
Definition: vcs_defs.h:204
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
static bool hasChargedSpecies(const ThermoPhase *const tPhase)
This function decides whether a phase has charged species or not.
static bool chargeNeutralityElement(const ThermoPhase *const tPhase)
This utility routine decides whether a Cantera ThermoPhase needs a constraint equation representing t...
#define VCS_PHASE_EXIST_YES
Phase is a normal phase that currently exists.
Definition: vcs_defs.h:188
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:31
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
Header file for the internal object that holds the vcs equilibrium problem (see Class VCS_SOLVE and E...
Header for the object representing each phase within vcs.
#define VCS_SPECIES_TYPE_MOLNUM
Unknown refers to mole number of a single species.
Definition: vcs_defs.h:281
#define VCS_ELEM_TYPE_ABSPOS
Normal element constraint consisting of positive coefficients for the formula matrix.
Definition: vcs_defs.h:228
Properties of a single species.
#define VCS_ELEM_TYPE_CHARGENEUTRALITY
This refers to a charge neutrality of a single phase.
Definition: vcs_defs.h:240
#define AssertThrow(expr, procedure)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:247
#define VCS_STATECALC_NEW
State Calculation based on the new or tentative mole numbers.
Definition: vcs_defs.h:303
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:264
#define VCS_STATECALC_OLD
State Calculation based on the old or base mole numbers.
Definition: vcs_defs.h:300
std::string id() const
Return the string id for the phase.
Definition: Phase.cpp:68
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:157
int elementType(size_t m) const
Return the element constraint type Possible types include:
Definition: Phase.cpp:149
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
Contains declarations for string manipulation functions within Cantera.
#define plogf
define this Cantera function to replace printf
Definition: vcs_internal.h:18
bool chargeNeutralityNecessary() const
Returns the chargeNeutralityNecessity boolean.
Definition: ThermoPhase.h:196
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE
Unknown refers to the voltage level of a phase.
Definition: vcs_defs.h:289
#define VCS_STATECALC_TMP
State Calculation based on a temporary set of mole numbers.
Definition: vcs_defs.h:311
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:161
Header file for class ThermoPhase, the base class for phases with thermodynamic properties, and the text for the Module thermoprops (see Thermodynamic Properties and class ThermoPhase).
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:577
std::string elementName(size_t m) const
Name of the element with index m.
Definition: Phase.cpp:107