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