Cantera  2.4.0
IonsFromNeutralVPSSTP.cpp
Go to the documentation of this file.
1 /**
2  * @file IonsFromNeutralVPSSTP.cpp
3  * Definitions for the object which treats ionic liquids as made of ions as species
4  * even though the thermodynamics is obtained from the neutral molecule representation.
5  * (see \ref thermoprops
6  * and class \link Cantera::IonsFromNeutralVPSSTP IonsFromNeutralVPSSTP\endlink).
7  *
8  * Header file for a derived class of ThermoPhase that handles variable pressure
9  * standard state methods for calculating thermodynamic properties that are
10  * further based upon expressions for the excess Gibbs free energy expressed as
11  * a function of the mole fractions.
12  */
13 
14 // This file is part of Cantera. See License.txt in the top-level directory or
15 // at http://www.cantera.org/license.txt for license and copyright information.
16 
21 
22 #include <fstream>
23 
24 using namespace std;
25 
26 namespace Cantera
27 {
28 
29 IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
30  ionSolnType_(cIonSolnType_SINGLEANION),
31  numNeutralMoleculeSpecies_(0),
32  indexSpecialSpecies_(npos),
33  neutralMoleculePhase_(0),
34  geThermo(0)
35 {
36 }
37 
38 IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile,
39  const std::string& id_) :
40  ionSolnType_(cIonSolnType_SINGLEANION),
41  numNeutralMoleculeSpecies_(0),
42  indexSpecialSpecies_(npos)
43 {
44  initThermoFile(inputFile, id_);
45 }
46 
48  const std::string& id_) :
49  ionSolnType_(cIonSolnType_SINGLEANION),
50  numNeutralMoleculeSpecies_(0),
51  indexSpecialSpecies_(npos)
52 {
53  importPhase(phaseRoot, this);
54 }
55 
56 // ------------ Molar Thermodynamic Properties ----------------------
57 
59 {
60  getPartialMolarEnthalpies(m_work.data());
61  return mean_X(m_work);
62 }
63 
65 {
66  getPartialMolarEntropies(m_work.data());
67  return mean_X(m_work);
68 }
69 
71 {
72  getChemPotentials(m_work.data());
73  return mean_X(m_work);
74 }
75 
77 {
78  getPartialMolarCp(m_work.data());
79  return mean_X(m_work);
80 }
81 
83 {
84  // Need to revisit this, as it is wrong
85  getPartialMolarCp(m_work.data());
86  return mean_X(m_work);
87 }
88 
89 // -- Activities, Standard States, Activity Concentrations -----------
90 
92  vector_fp& charges, std::vector<size_t>& neutMolIndex) const
93 {
94  coeffs = fm_neutralMolec_ions_;
95  charges = m_speciesCharge;
96  neutMolIndex = fm_invert_ionForNeutral;
97 }
98 
100 {
101  // Update the activity coefficients
103 
104  // take the exp of the internally stored coefficients.
105  for (size_t k = 0; k < m_kk; k++) {
106  ac[k] = exp(lnActCoeff_Scaled_[k]);
107  }
108 }
109 
110 // --------- Partial Molar Properties of the Solution -------------
111 
113 {
114  size_t icat, jNeut;
115  doublereal xx, fact2;
116 
117  // Get the standard chemical potentials of netural molecules
118  neutralMoleculePhase_->getStandardChemPotentials(muNeutralMolecule_.data());
119 
120  switch (ionSolnType_) {
121  case cIonSolnType_PASSTHROUGH:
122  neutralMoleculePhase_->getChemPotentials(mu);
123  break;
124  case cIonSolnType_SINGLEANION:
125  neutralMoleculePhase_->getLnActivityCoefficients(lnActCoeff_NeutralMolecule_.data());
126  fact2 = 2.0 * RT() * log(2.0);
127 
128  // Do the cation list
129  for (size_t k = 0; k < cationList_.size(); k++) {
130  // Get the id for the next cation
131  icat = cationList_[k];
132  jNeut = fm_invert_ionForNeutral[icat];
133  xx = std::max(SmallNumber, moleFractions_[icat]);
134  mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT() * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
135  }
136 
137  // Do the anion list
138  icat = anionList_[0];
139  jNeut = fm_invert_ionForNeutral[icat];
140  xx = std::max(SmallNumber, moleFractions_[icat]);
141  mu[icat] = RT() * log(xx);
142 
143  // Do the list of neutral molecules
144  for (size_t k = 0; k < passThroughList_.size(); k++) {
145  icat = passThroughList_[k];
146  jNeut = fm_invert_ionForNeutral[icat];
147  xx = std::max(SmallNumber, moleFractions_[icat]);
148  mu[icat] = muNeutralMolecule_[jNeut] + RT() * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
149  }
150  break;
151 
152  case cIonSolnType_SINGLECATION:
153  throw CanteraError("eosType", "Unknown type");
154  break;
155  case cIonSolnType_MULTICATIONANION:
156  throw CanteraError("eosType", "Unknown type");
157  break;
158  default:
159  throw CanteraError("eosType", "Unknown type");
160  break;
161  }
162 }
163 
165 {
166  // Get the nondimensional standard state enthalpies
167  getEnthalpy_RT(hbar);
168 
169  // dimensionalize it.
170  for (size_t k = 0; k < m_kk; k++) {
171  hbar[k] *= RT();
172  }
173 
174  // Update the activity coefficients, This also update the internally stored
175  // molalities.
178  for (size_t k = 0; k < m_kk; k++) {
179  hbar[k] -= RT() * temperature() * dlnActCoeffdT_Scaled_[k];
180  }
181 }
182 
184 {
185  // Get the nondimensional standard state entropies
186  getEntropy_R(sbar);
187 
188  // Update the activity coefficients, This also update the internally stored
189  // molalities.
192 
193  for (size_t k = 0; k < m_kk; k++) {
194  double xx = std::max(moleFractions_[k], SmallNumber);
195  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - temperature() * dlnActCoeffdT_Scaled_[k];
196  }
197 
198  // dimensionalize it.
199  for (size_t k = 0; k < m_kk; k++) {
200  sbar[k] *= GasConstant;
201  }
202 }
203 
204 void IonsFromNeutralVPSSTP::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
205 {
208 
209  for (size_t k = 0; k < m_kk; k++) {
210  dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
211  }
212 }
213 
214 void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const
215 {
218 
219  for (size_t k = 0; k < m_kk; k++) {
220  dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
221  }
222 }
223 
224 void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeffdlnN)
225 {
228  double* data = & dlnActCoeffdlnN_(0,0);
229  for (size_t k = 0; k < m_kk; k++) {
230  for (size_t m = 0; m < m_kk; m++) {
231  dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
232  }
233  }
234 }
235 
237 {
238  // This is a two phase process. First, we calculate the standard states
239  // within the neutral molecule phase.
240  neutralMoleculePhase_->setState_TP(temperature(), pressure());
241 
242  // Calculate the partial molar volumes, and then the density of the fluid
244 }
245 
246 void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
247 {
248  // Download the neutral mole fraction vector into the vector,
249  // NeutralMolecMoleFractions_[]
250  neutralMoleculePhase_->getMoleFractions(NeutralMolecMoleFractions_.data());
251 
252  // Zero the mole fractions
253  for (size_t k = 0; k < m_kk; k++) {
254  mf[k] = 0.0;
255  }
256 
257  // Use the formula matrix to calculate the relative mole numbers.
258  for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
259  for (size_t k = 0; k < m_kk; k++) {
260  double fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
261  mf[k] += fmij * NeutralMolecMoleFractions_[jNeut];
262  }
263  }
264 
265  // Normalize the new mole fractions
266  doublereal sum = 0.0;
267  for (size_t k = 0; k < m_kk; k++) {
268  sum += mf[k];
269  }
270  for (size_t k = 0; k < m_kk; k++) {
271  mf[k] /= sum;
272  }
273 }
274 
276 {
277  size_t icat, jNeut;
278  doublereal fmij;
279  doublereal sum = 0.0;
280 
281  // Zero the vector we are trying to find.
282  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
284  }
285  sum = -1.0;
286  for (size_t k = 0; k < m_kk; k++) {
287  sum += moleFractions_[k];
288  }
289  if (fabs(sum) > 1.0E-11) {
290  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
291  "molefracts don't sum to one: {}", sum);
292  }
293 
294  switch (ionSolnType_) {
295  case cIonSolnType_PASSTHROUGH:
296  for (size_t k = 0; k < m_kk; k++) {
298  }
299  break;
300 
301  case cIonSolnType_SINGLEANION:
302  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
304  }
305 
306  for (size_t k = 0; k < cationList_.size(); k++) {
307  // Get the id for the next cation
308  icat = cationList_[k];
309  jNeut = fm_invert_ionForNeutral[icat];
310  if (jNeut != npos) {
311  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
312  AssertTrace(fmij != 0.0);
313  NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
314  }
315  }
316 
317  for (size_t k = 0; k < passThroughList_.size(); k++) {
318  icat = passThroughList_[k];
319  jNeut = fm_invert_ionForNeutral[icat];
320  fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
321  NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
322  }
323 
324  for (size_t k = 0; k < m_kk; k++) {
326  }
327  for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
328  for (size_t k = 0; k < m_kk; k++) {
329  fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
330  moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut];
331  }
332  }
333  for (size_t k = 0; k < m_kk; k++) {
334  if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
335  // Check to see if we have in fact found the inverse.
336  if (anionList_[0] != k) {
337  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
338  "neutral molecule calc error");
339  } else {
340  // For the single anion case, we will allow some slippage
341  if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
342  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
343  "neutral molecule calc error - anion");
344  }
345  }
346  }
347  }
348 
349  // Normalize the Neutral Molecule mole fractions
350  sum = 0.0;
351  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
352  sum += NeutralMolecMoleFractions_[k];
353  }
354  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
355  NeutralMolecMoleFractions_[k] /= sum;
356  }
357  break;
358 
359  case cIonSolnType_SINGLECATION:
360  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", "Unknown type");
361  break;
362  case cIonSolnType_MULTICATIONANION:
363  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", "Unknown type");
364  break;
365  default:
366  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions", "Unknown type");
367  break;
368  }
369 }
370 
371 void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
372 {
373  doublereal sumy, sumdy;
374 
375  // check sum dx = 0
376  // Zero the vector we are trying to find.
377  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
378  y_[k] = 0.0;
379  dy[k] = 0.0;
380  }
381 
382  switch (ionSolnType_) {
383 
384  case cIonSolnType_PASSTHROUGH:
385  for (size_t k = 0; k < m_kk; k++) {
386  dy[k] = dx[k];
387  }
388  break;
389 
390  case cIonSolnType_SINGLEANION:
391  for (size_t k = 0; k < cationList_.size(); k++) {
392  // Get the id for the next cation
393  size_t icat = cationList_[k];
394  size_t jNeut = fm_invert_ionForNeutral[icat];
395  if (jNeut != npos) {
396  double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
397  AssertTrace(fmij != 0.0);
398  const doublereal temp = 1.0/fmij;
399  dy[jNeut] += dx[icat] * temp;
400  y_[jNeut] += moleFractions_[icat] * temp;
401  }
402  }
403 
404  for (size_t k = 0; k < passThroughList_.size(); k++) {
405  size_t icat = passThroughList_[k];
406  size_t jNeut = fm_invert_ionForNeutral[icat];
407  double fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
408  const doublereal temp = 1.0/fmij;
409  dy[jNeut] += dx[icat] * temp;
410  y_[jNeut] += moleFractions_[icat] * temp;
411  }
412  // Normalize the Neutral Molecule mole fractions
413  sumy = 0.0;
414  sumdy = 0.0;
415  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
416  sumy += y_[k];
417  sumdy += dy[k];
418  }
419  sumy = 1.0 / sumy;
420  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
421  dy[k] = dy[k] * sumy - y_[k]*sumdy*sumy*sumy;
422  }
423 
424  break;
425 
426  case cIonSolnType_SINGLECATION:
427  throw CanteraError("IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads",
428  "Unknown type");
429  break;
430  case cIonSolnType_MULTICATIONANION:
431  throw CanteraError("IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads",
432  "Unknown type");
433  break;
434  default:
435  throw CanteraError("IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads",
436  "Unknown type");
437  break;
438  }
439 }
440 
442 {
445  neutralMoleculePhase_->setMoleFractions(NeutralMolecMoleFractions_.data());
446 }
447 
448 // ------------ Partial Molar Properties of the Solution ------------
449 
450 //! Return the factor overlap
451 /*!
452  * @param elnamesVN
453  * @param elemVectorN
454  * @param nElementsN
455  * @param elnamesVI
456  * @param elemVectorI
457  * @param nElementsI
458  */
459 static double factorOverlap(const std::vector<std::string>& elnamesVN ,
460  const vector_fp& elemVectorN,
461  const size_t nElementsN,
462  const std::vector<std::string>& elnamesVI ,
463  const vector_fp& elemVectorI,
464  const size_t nElementsI)
465 {
466  double fMax = 1.0E100;
467  for (size_t mi = 0; mi < nElementsI; mi++) {
468  if (elnamesVI[mi] != "E" && elemVectorI[mi] > 1.0E-13) {
469  double eiNum = elemVectorI[mi];
470  for (size_t mn = 0; mn < nElementsN; mn++) {
471  if (elnamesVI[mi] == elnamesVN[mn]) {
472  if (elemVectorN[mn] <= 1.0E-13) {
473  return 0.0;
474  }
475  fMax = std::min(fMax, elemVectorN[mn]/eiNum);
476  }
477  }
478  }
479  }
480  return fMax;
481 }
482 
484 {
485  size_t nElementsN = neutralMoleculePhase_->nElements();
486  const std::vector<std::string>& elnamesVN = neutralMoleculePhase_->elementNames();
487  vector_fp elemVectorN(nElementsN);
488 
489  size_t nElementsI = nElements();
490  const std::vector<std::string>& elnamesVI = elementNames();
491  vector_fp elemVectorI(nElementsI);
492 
493  for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
494  for (size_t m = 0; m < nElementsN; m++) {
495  elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m);
496  }
497 
498  for (size_t m = 0; m < nElementsI; m++) {
499  elemVectorI[m] = nAtoms(indexSpecialSpecies_, m);
500  }
501  double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
502  elnamesVI ,elemVectorI, nElementsI);
503  if (fac > 0.0) {
504  for (size_t m = 0; m < nElementsN; m++) {
505  for (size_t mi = 0; mi < nElementsI; mi++) {
506  if (elnamesVN[m] == elnamesVI[mi]) {
507  elemVectorN[m] -= fac * elemVectorI[mi];
508  }
509 
510  }
511  }
512  }
514 
515  for (size_t k = 0; k < m_kk; k++) {
516  for (size_t m = 0; m < nElementsI; m++) {
517  elemVectorI[m] = nAtoms(k, m);
518  }
519  fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
520  elnamesVI ,elemVectorI, nElementsI);
521  if (fac > 0.0) {
522  for (size_t m = 0; m < nElementsN; m++) {
523  for (size_t mi = 0; mi < nElementsI; mi++) {
524  if (elnamesVN[m] == elnamesVI[mi]) {
525  elemVectorN[m] -= fac * elemVectorI[mi];
526  }
527  }
528  }
529  bool notTaken = true;
530  for (size_t iNeut = 0; iNeut < jNeut; iNeut++) {
531  if (fm_invert_ionForNeutral[k] == iNeut) {
532  notTaken = false;
533  }
534  }
535  if (notTaken) {
536  fm_invert_ionForNeutral[k] = jNeut;
537  } else {
538  throw CanteraError("IonsFromNeutralVPSSTP::initThermo",
539  "Simple formula matrix generation failed, one cation is shared between two salts");
540  }
541  }
542  fm_neutralMolec_ions_[k + jNeut * m_kk] += fac;
543  }
544 
545  // Ok check the work
546  for (size_t m = 0; m < nElementsN; m++) {
547  if (fabs(elemVectorN[m]) > 1.0E-13) {
548  throw CanteraError("IonsFromNeutralVPSSTP::initThermo",
549  "Simple formula matrix generation failed");
550  }
551  }
552  }
553 
555 }
556 
557 void IonsFromNeutralVPSSTP::setNeutralMoleculePhase(shared_ptr<ThermoPhase> neutral)
558 {
559  neutralMoleculePhase_ = neutral;
560  geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_.get());
570  y_.resize(numNeutralMoleculeSpecies_, 0.0);
571  dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
572  dX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
573 }
574 
575 shared_ptr<ThermoPhase> IonsFromNeutralVPSSTP::getNeutralMoleculePhase()
576 {
577  return neutralMoleculePhase_;
578 }
579 
580 bool IonsFromNeutralVPSSTP::addSpecies(shared_ptr<Species> spec)
581 {
582  bool added = GibbsExcessVPSSTP::addSpecies(spec);
583  if (added) {
584  moleFractions_.push_back(0.0);
585  moleFractionsTmp_.push_back(0.0);
586  m_work.push_back(0.0);
587  fm_invert_ionForNeutral.push_back(npos);
589 
590  if (spec->charge > 0) {
591  cationList_.push_back(m_kk-1);
592  } else if (spec->charge < 0) {
593  anionList_.push_back(m_kk-1);
594  } else {
595  passThroughList_.push_back(m_kk-1);
596  }
597 
598  if (spec->extra.hasKey("special_species")
599  && spec->extra["special_species"].asBool()) {
601  }
602  }
603  return added;
604 }
605 
607 {
609  // Find the Neutral Molecule Phase
610  if (!thermoNode.hasChild("neutralMoleculePhase")) {
611  throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
612  "no neutralMoleculePhase XML node");
613  }
614  XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");
615 
616  XML_Node* neut_ptr = get_XML_Node(neutralMoleculeNode["datasrc"], 0);
617  if (!neut_ptr) {
618  throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
619  "neut_ptr = 0");
620  }
621 
622  setNeutralMoleculePhase(shared_ptr<ThermoPhase>(newPhase(*neut_ptr)));
623 }
624 
626 {
627  size_t icat, jNeut;
628  // Get the activity coefficiens of the neutral molecules
629  neutralMoleculePhase_->getLnActivityCoefficients(lnActCoeff_NeutralMolecule_.data());
630 
631  switch (ionSolnType_) {
632  case cIonSolnType_PASSTHROUGH:
633  break;
634  case cIonSolnType_SINGLEANION:
635  // Do the cation list
636  for (size_t k = 0; k < cationList_.size(); k++) {
637  // Get the id for the next cation
638  icat = cationList_[k];
639  jNeut = fm_invert_ionForNeutral[icat];
640  double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
641  lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut] / fmij;
642  }
643 
644  // Do the anion list
645  icat = anionList_[0];
646  jNeut = fm_invert_ionForNeutral[icat];
647  lnActCoeff_Scaled_[icat]= 0.0;
648 
649  // Do the list of neutral molecules
650  for (size_t k = 0; k < passThroughList_.size(); k++) {
651  icat = passThroughList_[k];
652  jNeut = fm_invert_ionForNeutral[icat];
654  }
655  break;
656 
657  case cIonSolnType_SINGLECATION:
658  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
659  break;
660  case cIonSolnType_MULTICATIONANION:
661  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
662  break;
663  default:
664  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
665  break;
666  }
667 }
668 
669 void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
670  doublereal* dlnActCoeffds) const
671 {
672  size_t icat, jNeut;
673  // Get the activity coefficients of the neutral molecules
674  if (!geThermo) {
675  for (size_t k = 0; k < m_kk; k++) {
676  dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
677  }
678  return;
679  }
680 
681  getNeutralMoleculeMoleGrads(dXds, dX_NeutralMolecule_.data());
682 
683  // All mole fractions returned to normal
684  geThermo->getdlnActCoeffds(dTds, dX_NeutralMolecule_.data(), dlnActCoeff_NeutralMolecule_.data());
685 
686  switch (ionSolnType_) {
687  case cIonSolnType_PASSTHROUGH:
688  break;
689  case cIonSolnType_SINGLEANION:
690  // Do the cation list
691  for (size_t k = 0; k < cationList_.size(); k++) {
692  // Get the id for the next cation
693  icat = cationList_[k];
694  jNeut = fm_invert_ionForNeutral[icat];
695  double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
696  dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule_[jNeut]/fmij;
697  }
698 
699  // Do the anion list
700  icat = anionList_[0];
701  jNeut = fm_invert_ionForNeutral[icat];
702  dlnActCoeffds[icat]= 0.0;
703 
704  // Do the list of neutral molecules
705  for (size_t k = 0; k < passThroughList_.size(); k++) {
706  icat = passThroughList_[k];
707  jNeut = fm_invert_ionForNeutral[icat];
708  dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule_[jNeut];
709  }
710  break;
711 
712  case cIonSolnType_SINGLECATION:
713  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
714  break;
715  case cIonSolnType_MULTICATIONANION:
716  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
717  break;
718  default:
719  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
720  break;
721  }
722 }
723 
725 {
726  size_t icat, jNeut;
727 
728  // Get the activity coefficients of the neutral molecules
729  if (!geThermo) {
730  dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
731  return;
732  }
733 
735 
736  switch (ionSolnType_) {
737  case cIonSolnType_PASSTHROUGH:
738  break;
739  case cIonSolnType_SINGLEANION:
740  // Do the cation list
741  for (size_t k = 0; k < cationList_.size(); k++) {
742  //! Get the id for the next cation
743  icat = cationList_[k];
744  jNeut = fm_invert_ionForNeutral[icat];
745  double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
747  }
748 
749  // Do the anion list
750  icat = anionList_[0];
751  jNeut = fm_invert_ionForNeutral[icat];
752  dlnActCoeffdT_Scaled_[icat]= 0.0;
753 
754  // Do the list of neutral molecules
755  for (size_t k = 0; k < passThroughList_.size(); k++) {
756  icat = passThroughList_[k];
757  jNeut = fm_invert_ionForNeutral[icat];
759  }
760  break;
761 
762  case cIonSolnType_SINGLECATION:
763  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
764  break;
765  case cIonSolnType_MULTICATIONANION:
766  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
767  break;
768  default:
769  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
770  break;
771  }
772 }
773 
775 {
776  size_t icat, jNeut;
777 
778  // Get the activity coefficients of the neutral molecules
779  if (!geThermo) {
780  dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
781  return;
782  }
783 
785 
786  switch (ionSolnType_) {
787  case cIonSolnType_PASSTHROUGH:
788  break;
789  case cIonSolnType_SINGLEANION:
790  // Do the cation list
791  for (size_t k = 0; k < cationList_.size(); k++) {
792  // Get the id for the next cation
793  icat = cationList_[k];
794  jNeut = fm_invert_ionForNeutral[icat];
795  double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
797  }
798 
799  // Do the anion list
800  icat = anionList_[0];
801  jNeut = fm_invert_ionForNeutral[icat];
802  dlnActCoeffdlnX_diag_[icat]= 0.0;
803 
804  // Do the list of neutral molecules
805  for (size_t k = 0; k < passThroughList_.size(); k++) {
806  icat = passThroughList_[k];
807  jNeut = fm_invert_ionForNeutral[icat];
809  }
810  break;
811 
812  case cIonSolnType_SINGLECATION:
813  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
814  break;
815  case cIonSolnType_MULTICATIONANION:
816  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
817  break;
818  default:
819  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
820  break;
821  }
822 }
823 
825 {
826  size_t icat, jNeut;
827 
828  // Get the activity coefficients of the neutral molecules
829  if (!geThermo) {
830  dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
831  return;
832  }
833 
835 
836  switch (ionSolnType_) {
837  case cIonSolnType_PASSTHROUGH:
838  break;
839  case cIonSolnType_SINGLEANION:
840  // Do the cation list
841  for (size_t k = 0; k < cationList_.size(); k++) {
842  // Get the id for the next cation
843  icat = cationList_[k];
844  jNeut = fm_invert_ionForNeutral[icat];
845  double fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
847  }
848 
849  // Do the anion list
850  icat = anionList_[0];
851  jNeut = fm_invert_ionForNeutral[icat];
852  dlnActCoeffdlnN_diag_[icat]= 0.0;
853 
854  // Do the list of neutral molecules
855  for (size_t k = 0; k < passThroughList_.size(); k++) {
856  icat = passThroughList_[k];
857  jNeut = fm_invert_ionForNeutral[icat];
859  }
860  break;
861 
862  case cIonSolnType_SINGLECATION:
863  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
864  break;
865  case cIonSolnType_MULTICATIONANION:
866  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
867  break;
868  default:
869  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
870  break;
871  }
872 }
873 
875 {
876  size_t kcat = 0, kNeut = 0, mcat = 0, mNeut = 0;
877  doublereal fmij = 0.0;
879  // Get the activity coefficients of the neutral molecules
880  if (!geThermo) {
881  throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed");
882  }
883  size_t nsp_ge = geThermo->nSpecies();
884  geThermo->getdlnActCoeffdlnN(nsp_ge, &dlnActCoeffdlnN_NeutralMolecule_(0,0));
885 
886  switch (ionSolnType_) {
887  case cIonSolnType_PASSTHROUGH:
888  break;
889  case cIonSolnType_SINGLEANION:
890  // Do the cation list
891  for (size_t k = 0; k < cationList_.size(); k++) {
892  for (size_t m = 0; m < cationList_.size(); m++) {
893  kcat = cationList_[k];
894 
895  kNeut = fm_invert_ionForNeutral[kcat];
896  fmij = fm_neutralMolec_ions_[kcat + kNeut * m_kk];
898 
899  mcat = cationList_[m];
900  mNeut = fm_invert_ionForNeutral[mcat];
901  double mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
902 
903  dlnActCoeffdlnN_(kcat,mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij;
904 
905  }
906  for (size_t m = 0; m < passThroughList_.size(); m++) {
907  mcat = passThroughList_[m];
908  mNeut = fm_invert_ionForNeutral[mcat];
909  dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut) / fmij;
910  }
911  }
912 
913  // Do the anion list -> anion activity coefficient is one
914  kcat = anionList_[0];
915  kNeut = fm_invert_ionForNeutral[kcat];
916  for (size_t k = 0; k < m_kk; k++) {
917  dlnActCoeffdlnN_(kcat, k) = 0.0;
918  dlnActCoeffdlnN_(k, kcat) = 0.0;
919  }
920 
921  // Do the list of neutral molecules
922  for (size_t k = 0; k < passThroughList_.size(); k++) {
923  kcat = passThroughList_[k];
924  kNeut = fm_invert_ionForNeutral[kcat];
926 
927  for (size_t m = 0; m < m_kk; m++) {
928  mcat = passThroughList_[m];
929  mNeut = fm_invert_ionForNeutral[mcat];
930  dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut);
931  }
932 
933  for (size_t m = 0; m < cationList_.size(); m++) {
934  mcat = cationList_[m];
935  mNeut = fm_invert_ionForNeutral[mcat];
936  dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut);
937  }
938  }
939  break;
940 
941  case cIonSolnType_SINGLECATION:
942  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
943  break;
944  case cIonSolnType_MULTICATIONANION:
945  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
946  break;
947  default:
948  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
949  break;
950  }
951 }
952 
953 }
virtual void getdlnActCoeffdlnN(const size_t ld, doublereal *const dlnActCoeffdlnN)
Get the array of derivatives of the log activity coefficients with respect to the log of the species ...
XML_Node * get_XML_Node(const std::string &file_ID, XML_Node *root)
This routine will locate an XML node in either the input XML tree or in another input file specified ...
Definition: global.cpp:189
vector_fp muNeutralMolecule_
Storage vector for the neutral molecule chemical potentials.
ThermoPhase * newPhase(XML_Node &xmlphase)
Create a new ThermoPhase object and initializes it according to the XML tree.
size_t nElements() const
Number of elements.
Definition: Phase.cpp:88
vector_fp NeutralMolecMoleFractions_
Mole fractions using the Neutral Molecule Mole fraction basis.
vector_fp dlnActCoeffdlnX_diag_
Storage for the current derivative values of the gradients with respect to logarithm of the mole frac...
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with &#39;v&#39;.
Definition: Array.h:112
Header for intermediate ThermoPhase object for phases which consist of ions whose thermodynamics is c...
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
vector_fp m_speciesCharge
Vector of species charges. length m_kk.
Definition: Phase.h:799
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
size_t indexSpecialSpecies_
Index of special species.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:614
vector_fp lnActCoeff_NeutralMolecule_
Storage vector for the neutral molecule ln activity coefficients.
vector_fp dlnActCoeffdT_NeutralMolecule_
Storage vector for the neutral molecule d ln activity coefficients dT.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:748
void s_update_dlnActCoeff_dlnX_diag() const
Update the derivative of the log of the activity coefficients wrt log(mole fraction) ...
vector_fp dlnActCoeffdlnN_diag_NeutralMolecule_
Storage vector for the neutral molecule d ln activity coefficients dlnN.
virtual void calcNeutralMoleculeMoleFractions() const
Calculate neutral molecule mole fractions.
vector_fp dlnActCoeffdlnN_diag_
Storage for the current derivative values of the gradients with respect to logarithm of the mole frac...
shared_ptr< ThermoPhase > neutralMoleculePhase_
This is a pointer to the neutral Molecule Phase.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual void getdlnActCoeffds(const doublereal dTds, const doublereal *const dXds, doublereal *dlnActCoeffds) const
Get the change in activity coefficients wrt changes in state (temp, mole fraction, etc) along a line in parameter space or along a line in physical space.
void s_update_dlnActCoeffdT() const
Update the temperature derivative of the ln activity coefficients.
IonSolnType_enumType ionSolnType_
Ion solution type.
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const
Get the array of log species mole number derivatives of the log activity coefficients.
virtual void getdlnActCoeffdlnN(const size_t ld, doublereal *const dlnActCoeffdlnN)
Get the array of derivatives of the log activity coefficients with respect to the log of the species ...
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
void s_update_lnActCoeff() const
Update the activity coefficients.
Array2D dlnActCoeffdlnN_NeutralMolecule_
Storage vector for the neutral molecule d ln activity coefficients dlnN.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
virtual bool addSpecies(shared_ptr< Species > spec)
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
virtual bool addSpecies(shared_ptr< Species > spec)
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies for the species in the mixture.
void s_update_dlnActCoeff_dlnN() const
Update the derivative of the log of the activity coefficients wrt log(number of moles) - diagonal com...
std::vector< size_t > anionList_
List of the species in this ThermoPhase which are anion species.
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
vector_fp lnActCoeff_Scaled_
Storage for the current values of the activity coefficients of the species.
virtual doublereal pressure() const
Returns the current pressure of the phase.
std::vector< size_t > passThroughList_
List of the species in this ThermoPhase which are passed through to the neutralMoleculePhase ThermoPh...
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:536
XML_Node & child(const size_t n) const
Return a changeable reference to the n&#39;th child of the current node.
Definition: xml.cpp:546
#define AssertTrace(expr)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:233
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
virtual void calcIonMoleFractions(doublereal *const mf) const
Calculate ion mole fractions from neutral molecule mole fractions.
virtual void getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const
Get the array of ln mole fraction derivatives of the log activity coefficients - diagonal component o...
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
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
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const
Get the array of log species mole number derivatives of the log activity coefficients.
std::vector< size_t > fm_invert_ionForNeutral
Mapping between ion species and neutral molecule for quick invert.
virtual doublereal enthalpy_mole() const
Return the Molar enthalpy. Units: J/kmol.
virtual void setParametersFromXML(const XML_Node &eosdata)
Set equation of state parameter values from XML entries.
Definition: ThermoPhase.h:1468
virtual void setParametersFromXML(const XML_Node &thermoNode)
Set equation of state parameter values from XML entries.
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
Contains declarations for string manipulation functions within Cantera.
std::vector< size_t > cationList_
List of the species in this ThermoPhase which are cation species.
vector_fp moleFractions_
Storage for the current values of the mole fractions of the species.
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional molar-based activity coefficients at the current solution temperatur...
vector_fp fm_neutralMolec_ions_
Formula Matrix for composition of neutral molecules in terms of the molecules in this ThermoPhase...
size_t m_kk
Number of species in the phase.
Definition: Phase.h:788
Array2D dlnActCoeffdlnN_
Storage for the current derivative values of the gradients with respect to logarithm of the species m...
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
void zero()
Set all of the entries to zero.
Definition: Array.h:198
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
void getDissociationCoeffs(vector_fp &fm_neutralMolec_ions, vector_fp &charges, std::vector< size_t > &neutMolIndex) const
Get the Salt Dissociation Coefficients.
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
Definition: ThermoPhase.h:518
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
static double factorOverlap(const std::vector< std::string > &elnamesVN, const vector_fp &elemVectorN, const size_t nElementsN, const std::vector< std::string > &elnamesVI, const vector_fp &elemVectorI, const size_t nElementsI)
Return the factor overlap.
Declarations for the class PDSS_IonsFromNeutral ( which handles calculations for a single ion in a fl...
doublereal nAtoms(size_t k, size_t m) const
Number of atoms of element m in species k.
Definition: Phase.cpp:161
size_t numNeutralMoleculeSpecies_
Number of neutral molecule species.
virtual void getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const
Get the array of ln mole fraction derivatives of the log activity coefficients - diagonal component o...
Definition: ThermoPhase.h:1520
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual void getdlnActCoeffds(const doublereal dTds, const doublereal *const dXds, doublereal *dlnActCoeffds) const
Get the change in activity coefficients wrt changes in state (temp, mole fraction, etc) along a line in parameter space or along a line in physical space.
Definition: ThermoPhase.h:1500
const std::vector< std::string > & elementNames() const
Return a read-only reference to the vector of element names.
Definition: Phase.cpp:123
vector_fp moleFractionsTmp_
Temporary mole fraction vector.
void getNeutralMoleculeMoleGrads(const doublereal *const dx, doublereal *const dy) const
Calculate neutral molecule mole fractions.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const
Get the array of temperature derivatives of the log activity coefficients.
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
Definition: Phase.h:622
vector_fp dlnActCoeffdT_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
void s_update_dlnActCoeff_dlnN_diag() const
Update the derivative of the log of the activity coefficients wrt log(number of moles) - diagonal com...
vector_fp dlnActCoeffdlnX_diag_NeutralMolecule_
Storage vector for the neutral molecule d ln activity coefficients dX - diagonal component.