Cantera  2.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
9  * variable pressure standard state methods for calculating
10  * thermodynamic properties that are further based upon expressions
11  * for the excess gibbs free energy expressed as a function of
12  * the mole fractions.
13  */
14 /*
15  * Copyright (2009) Sandia Corporation. Under the terms of
16  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
17  * U.S. Government retains certain rights in this software.
18  */
22 #include "cantera/thermo/mix_defs.h"
24 
25 #include <cmath>
26 #include <iomanip>
27 #include <fstream>
28 
29 using namespace std;
30 
31 namespace Cantera
32 {
33 
34 static const double xxSmall = 1.0E-150;
35 //====================================================================================================================
36 /*
37  * Default constructor.
38  *
39  */
40 IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
42  ionSolnType_(cIonSolnType_SINGLEANION),
43  numNeutralMoleculeSpecies_(0),
44  indexSpecialSpecies_(npos),
45  indexSecondSpecialSpecies_(npos),
46  numCationSpecies_(0),
47  numAnionSpecies_(0),
48  numPassThroughSpecies_(0),
49  neutralMoleculePhase_(0),
50  IOwnNThermoPhase_(true),
51  moleFractionsTmp_(0),
52  muNeutralMolecule_(0),
53  lnActCoeff_NeutralMolecule_(0)
54 {
55 }
56 
57 //====================================================================================================================
58 // Construct and initialize an IonsFromNeutralVPSSTP object
59 // directly from an ASCII input file
60 /*
61  * Working constructors
62  *
63  * The two constructors below are the normal way
64  * the phase initializes itself. They are shells that call
65  * the routine initThermo(), with a reference to the
66  * XML database to get the info for the phase.
67  *
68  * @param inputFile Name of the input file containing the phase XML data
69  * to set up the object
70  * @param id ID of the phase in the input file. Defaults to the
71  * empty string.
72  * @param neutralPhase The object takes a neutralPhase ThermoPhase
73  * object as input. It can either take a pointer
74  * to an existing object in the parameter list,
75  * in which case it does not own the object, or
76  * it can construct a neutral Phase as a slave
77  * object, in which case, it does own the slave
78  * object, for purposes of who gets to destroy
79  * the object.
80  * If this parameter is zero, then a slave
81  * neutral phase object is created and used.
82  */
83 IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(std::string inputFile, std::string id,
84  ThermoPhase* neutralPhase) :
86  ionSolnType_(cIonSolnType_SINGLEANION),
87  numNeutralMoleculeSpecies_(0),
88  indexSpecialSpecies_(npos),
89  indexSecondSpecialSpecies_(npos),
90  numCationSpecies_(0),
91  numAnionSpecies_(0),
92  numPassThroughSpecies_(0),
93  neutralMoleculePhase_(neutralPhase),
94  IOwnNThermoPhase_(true),
95  moleFractionsTmp_(0),
96  muNeutralMolecule_(0),
97  lnActCoeff_NeutralMolecule_(0)
98 {
99  if (neutralPhase) {
100  IOwnNThermoPhase_ = false;
101  }
102  constructPhaseFile(inputFile, id);
103 }
104 //====================================================================================================================
106  ThermoPhase* neutralPhase) :
108  ionSolnType_(cIonSolnType_SINGLEANION),
109  numNeutralMoleculeSpecies_(0),
110  indexSpecialSpecies_(npos),
111  indexSecondSpecialSpecies_(npos),
112  numCationSpecies_(0),
113  numAnionSpecies_(0),
114  numPassThroughSpecies_(0),
115  neutralMoleculePhase_(neutralPhase),
116  IOwnNThermoPhase_(true),
117  moleFractionsTmp_(0),
118  muNeutralMolecule_(0),
119 
120  lnActCoeff_NeutralMolecule_(0)
121 {
122  if (neutralPhase) {
123  IOwnNThermoPhase_ = false;
124  }
125  constructPhaseXML(phaseRoot, id);
126 }
127 
128 //====================================================================================================================
129 
130 /*
131  * Copy Constructor:
132  *
133  * Note this stuff will not work until the underlying phase
134  * has a working copy constructor
135  */
138  ionSolnType_(cIonSolnType_SINGLEANION),
139  numNeutralMoleculeSpecies_(0),
140  indexSpecialSpecies_(npos),
141  indexSecondSpecialSpecies_(npos),
142  numCationSpecies_(0),
143  numAnionSpecies_(0),
144  numPassThroughSpecies_(0),
145  neutralMoleculePhase_(0),
146  IOwnNThermoPhase_(true),
147  moleFractionsTmp_(0),
148  muNeutralMolecule_(0),
149 
150  lnActCoeff_NeutralMolecule_(0)
151 {
153 }
154 //====================================================================================================================
155 /*
156  * operator=()
157  *
158  * Note this stuff will not work until the underlying phase
159  * has a working assignment operator
160  */
163 {
164  if (&b == this) {
165  return *this;
166  }
167 
168  /*
169  * If we own the underlying neutral molecule phase, then we do a deep
170  * copy. If not, we do a shallow copy. We get a valid pointer for
171  * neutralMoleculePhase_ first, because we need it to assign the pointers
172  * within the PDSS_IonsFromNeutral object. which is done in the
173  * GibbsExcessVPSSTP::operator=(b) step.
174  */
175  if (IOwnNThermoPhase_) {
176  if (b.neutralMoleculePhase_) {
177  if (neutralMoleculePhase_) {
178  delete neutralMoleculePhase_;
179  }
181  } else {
183  }
184  } else {
186  }
187 
189 
203 
207  // gammaNeutralMolecule_ = b.gammaNeutralMolecule_;
213 
214  return *this;
215 }
216 
217 /*
218  *
219  * ~IonsFromNeutralVPSSTP(): (virtual)
220  *
221  * Destructor: does nothing:
222  *
223  */
225 {
226  if (IOwnNThermoPhase_) {
227  delete neutralMoleculePhase_;
229  }
230 }
231 
232 /*
233  * This routine duplicates the current object and returns
234  * a pointer to ThermoPhase.
235  */
238 {
240  return (ThermoPhase*) mtp;
241 }
242 
243 /*
244  * -------------- Utilities -------------------------------
245  */
246 
247 
248 // Equation of state type flag.
249 /*
250  * The ThermoPhase base class returns
251  * zero. Subclasses should define this to return a unique
252  * non-zero value. Known constants defined for this purpose are
253  * listed in mix_defs.h. The IonsFromNeutralVPSSTP class also returns
254  * zero, as it is a non-complete class.
255  */
257 {
258  return cIonsFromNeutral;
259 }
260 
261 
262 
263 /*
264  * ------------ Molar Thermodynamic Properties ----------------------
265  */
266 /*
267  * Molar enthalpy of the solution. Units: J/kmol.
268  */
270 {
272  return mean_X(DATA_PTR(m_pp));
273 }
274 
275 /**
276  * Molar internal energy of the solution. Units: J/kmol.
277  *
278  * This is calculated from the soln enthalpy and then
279  * subtracting pV.
280  */
282 {
283  double hh = enthalpy_mole();
284  double pres = pressure();
285  double molarV = 1.0/molarDensity();
286  double uu = hh - pres * molarV;
287  return uu;
288 }
289 
290 /**
291  * Molar soln entropy at constant pressure. Units: J/kmol/K.
292  *
293  * This is calculated from the partial molar entropies.
294  */
296 {
298  return mean_X(DATA_PTR(m_pp));
299 }
300 
301 /// Molar Gibbs function. Units: J/kmol.
303 {
305  return mean_X(DATA_PTR(m_pp));
306 }
307 /** Molar heat capacity at constant pressure. Units: J/kmol/K.
308  *
309  * Returns the solution heat capacition at constant pressure.
310  * This is calculated from the partial molar heat capacities.
311  */
313 {
315  double val = mean_X(DATA_PTR(m_pp));
316  return val;
317 }
318 
319 /// Molar heat capacity at constant volume. Units: J/kmol/K.
321 {
322  // Need to revisit this, as it is wrong
324  return mean_X(DATA_PTR(m_pp));
325  //err("not implemented");
326  //return 0.0;
327 }
328 //===========================================================================================================
329 /*
330  * - Activities, Standard States, Activity Concentrations -----------
331  */
332 //===========================================================================================================
334  vector_fp& charges, std::vector<size_t>& neutMolIndex) const
335 {
336  coeffs = fm_neutralMolec_ions_;
337  charges = m_speciesCharge;
338  neutMolIndex = fm_invert_ionForNeutral;
339 }
340 //===========================================================================================================
341 // Get the array of non-dimensional molar-based activity coefficients at
342 // the current solution temperature, pressure, and solution concentration.
343 /*
344  * @param ac Output vector of activity coefficients. Length: m_kk.
345  */
347 {
348 
349  // This stuff has moved to the setState routines
350  // calcNeutralMoleculeMoleFractions();
351  // neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_));
352  // neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_));
353 
354  /*
355  * Update the activity coefficients
356  */
358 
359  /*
360  * take the exp of the internally stored coefficients.
361  */
362  for (size_t k = 0; k < m_kk; k++) {
363  ac[k] = exp(lnActCoeff_Scaled_[k]);
364  }
365 }
366 
367 /*
368  * --------- Partial Molar Properties of the Solution -------------------------------
369  */
370 
371 // Get the species chemical potentials. Units: J/kmol.
372 /*
373  * This function returns a vector of chemical potentials of the
374  * species in solution at the current temperature, pressure
375  * and mole fraction of the solution.
376  *
377  * @param mu Output vector of species chemical
378  * potentials. Length: m_kk. Units: J/kmol
379  */
380 void
382 {
383  size_t icat, jNeut;
384  doublereal xx, fact2;
385  /*
386  * Transfer the mole fractions to the slave neutral molecule
387  * phase
388  * Note we may move this in the future.
389  */
390  //calcNeutralMoleculeMoleFractions();
391  //neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_));
392 
393  /*
394  * Get the standard chemical potentials of netural molecules
395  */
397 
398  doublereal RT_ = GasConstant * temperature();
399 
400  switch (ionSolnType_) {
401  case cIonSolnType_PASSTHROUGH:
403  break;
404  case cIonSolnType_SINGLEANION:
405  // neutralMoleculePhase_->getActivityCoefficients(DATA_PTR(gammaNeutralMolecule_));
407 
408  fact2 = 2.0 * RT_ * log(2.0);
409 
410  // Do the cation list
411  for (size_t k = 0; k < cationList_.size(); k++) {
412  //! Get the id for the next cation
413  icat = cationList_[k];
414  jNeut = fm_invert_ionForNeutral[icat];
415  xx = std::max(SmallNumber, moleFractions_[icat]);
416  mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
417  }
418 
419  // Do the anion list
420  icat = anionList_[0];
421  jNeut = fm_invert_ionForNeutral[icat];
422  xx = std::max(SmallNumber, moleFractions_[icat]);
423  mu[icat] = RT_ * log(xx);
424 
425  // Do the list of neutral molecules
426  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
427  icat = passThroughList_[k];
428  jNeut = fm_invert_ionForNeutral[icat];
429  xx = std::max(SmallNumber, moleFractions_[icat]);
430  mu[icat] = muNeutralMolecule_[jNeut] + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
431  }
432  break;
433 
434  case cIonSolnType_SINGLECATION:
435  throw CanteraError("eosType", "Unknown type");
436  break;
437  case cIonSolnType_MULTICATIONANION:
438  throw CanteraError("eosType", "Unknown type");
439  break;
440  default:
441  throw CanteraError("eosType", "Unknown type");
442  break;
443  }
444 }
445 
446 
447 // Returns an array of partial molar enthalpies for the species
448 // in the mixture.
449 /*
450  * Units (J/kmol)
451  *
452  * For this phase, the partial molar enthalpies are equal to the
453  * standard state enthalpies modified by the derivative of the
454  * molality-based activity coefficient wrt temperature
455  *
456  * \f[
457  * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
458  * \f]
459  *
460  */
462 {
463  /*
464  * Get the nondimensional standard state enthalpies
465  */
466  getEnthalpy_RT(hbar);
467  /*
468  * dimensionalize it.
469  */
470  double T = temperature();
471  double RT = GasConstant * T;
472  for (size_t k = 0; k < m_kk; k++) {
473  hbar[k] *= RT;
474  }
475  /*
476  * Update the activity coefficients, This also update the
477  * internally stored molalities.
478  */
481  double RTT = RT * T;
482  for (size_t k = 0; k < m_kk; k++) {
483  hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
484  }
485 }
486 
487 // Returns an array of partial molar entropies for the species
488 // in the mixture.
489 /*
490  * Units (J/kmol)
491  *
492  * For this phase, the partial molar enthalpies are equal to the
493  * standard state enthalpies modified by the derivative of the
494  * activity coefficient wrt temperature
495  *
496  * \f[
497  * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
498  * \f]
499  *
500  */
502 {
503  double xx;
504  /*
505  * Get the nondimensional standard state entropies
506  */
507  getEntropy_R(sbar);
508  double T = temperature();
509  /*
510  * Update the activity coefficients, This also update the
511  * internally stored molalities.
512  */
515 
516  for (size_t k = 0; k < m_kk; k++) {
517  xx = std::max(moleFractions_[k], xxSmall);
518  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
519  }
520  /*
521  * dimensionalize it.
522  */
523  for (size_t k = 0; k < m_kk; k++) {
524  sbar[k] *= GasConstant;
525  }
526 }
527 
528 
529 // Get the array of log concentration-like derivatives of the
530 // log activity coefficients
531 /*
532  * This function is a virtual method. For ideal mixtures
533  * (unity activity coefficients), this can return zero.
534  * Implementations should take the derivative of the
535  * logarithm of the activity coefficient with respect to the
536  * logarithm of the concentration-like variable (i.e. mole fraction,
537  * molality, etc.) that represents the standard state.
538  * This quantity is to be used in conjunction with derivatives of
539  * that concentration-like variable when the derivative of the chemical
540  * potential is taken.
541  *
542  * units = dimensionless
543  *
544  * @param dlnActCoeffdlnX Output vector of log(mole fraction)
545  * derivatives of the log Activity Coefficients.
546  * length = m_kk
547  */
548 void IonsFromNeutralVPSSTP::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
549 {
552 
553  for (size_t k = 0; k < m_kk; k++) {
554  dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
555  }
556 }
557 //====================================================================================================================
558 // Get the array of log concentration-like derivatives of the
559 // log activity coefficients
560 /*
561  * This function is a virtual method. For ideal mixtures
562  * (unity activity coefficients), this can return zero.
563  * Implementations should take the derivative of the
564  * logarithm of the activity coefficient with respect to the
565  * logarithm of the concentration-like variable (i.e. moles)
566  * that represents the standard state.
567  * This quantity is to be used in conjunction with derivatives of
568  * that concentration-like variable when the derivative of the chemical
569  * potential is taken.
570  *
571  * units = dimensionless
572  *
573  * @param dlnActCoeffdlnN_diag Output vector of log(mole fraction)
574  * derivatives of the log Activity Coefficients.
575  * length = m_kk
576  */
577 void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const
578 {
581 
582  for (size_t k = 0; k < m_kk; k++) {
583  dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
584  }
585 }
586 //====================================================================================================================
587 void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeffdlnN)
588 {
591  double* data = & dlnActCoeffdlnN_(0,0);
592  for (size_t k = 0; k < m_kk; k++) {
593  for (size_t m = 0; m < m_kk; m++) {
594  dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
595  }
596  }
597 }
598 //====================================================================================================================
599 void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp)
600 {
601  double p = pressure();
603 }
604 //====================================================================================================================
606 {
607  double t = temperature();
609 }
610 //====================================================================================================================
611 // Set the temperature (K) and pressure (Pa)
612 /*
613  * Setting the pressure may involve the solution of a nonlinear equation.
614  *
615  * @param t Temperature (K)
616  * @param p Pressure (Pa)
617  */
618 void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p)
619 {
620  /*
621  * This is a two phase process. First, we calculate the standard states
622  * within the neutral molecule phase.
623  */
626 
627  /*
628  * Calculate the partial molar volumes, and then the density of the fluid
629  */
630 
631  //calcDensity();
632  double dd = neutralMoleculePhase_->density();
633  Phase::setDensity(dd);
634 }
635 
636 // Calculate ion mole fractions from neutral molecule
637 // mole fractions.
638 /*
639  * @param mf Dump the mole fractions into this vector.
640  */
641 void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
642 {
643  doublereal fmij;
644  /*
645  * Download the neutral mole fraction vector into the
646  * vector, NeutralMolecMoleFractions_[]
647  */
649 
650  // Zero the mole fractions
651  for (size_t k = 0; k < m_kk; k++) {
652  mf[k] = 0.0;
653  }
654 
655  /*
656  * Use the formula matrix to calculate the relative mole numbers.
657  */
658  for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
659  for (size_t k = 0; k < m_kk; k++) {
660  fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
661  mf[k] += fmij * NeutralMolecMoleFractions_[jNeut];
662  }
663  }
664 
665  /*
666  * Normalize the new mole fractions
667  */
668  doublereal sum = 0.0;
669  for (size_t k = 0; k < m_kk; k++) {
670  sum += mf[k];
671  }
672  for (size_t k = 0; k < m_kk; k++) {
673  mf[k] /= sum;
674  }
675 
676 }
677 //====================================================================================================================
678 // Calculate neutral molecule mole fractions
679 /*
680  * This routine calculates the neutral molecule mole
681  * fraction given the vector of ion mole fractions,
682  * i.e., the mole fractions from this ThermoPhase.
683  * Note, this routine basically assumes that there
684  * is charge neutrality. If there isn't, then it wouldn't
685  * make much sense.
686  *
687  * for the case of cIonSolnType_SINGLEANION, some slough
688  * in the charge neutrality is allowed. The cation number
689  * is followed, while the difference in charge neutrality
690  * is dumped into the anion mole number to fix the imbalance.
691  */
693 {
694  size_t icat, jNeut;
695  doublereal fmij;
696  doublereal sum = 0.0;
697 
698  //! Zero the vector we are trying to find.
699  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
701  }
702 #ifdef DEBUG_MODE
703  sum = -1.0;
704  for (int k = 0; k < m_kk; k++) {
705  sum += moleFractions_[k];
706  }
707  if (fabs(sum) > 1.0E-11) {
708  throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
709  "molefracts don't sum to one: " + fp2str(sum));
710  }
711 #endif
712 
713  // bool fmSimple = true;
714 
715  switch (ionSolnType_) {
716 
717  case cIonSolnType_PASSTHROUGH:
718  for (size_t k = 0; k < m_kk; k++) {
720  }
721  break;
722 
723  case cIonSolnType_SINGLEANION:
724  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
726  }
727 
728  for (size_t k = 0; k < cationList_.size(); k++) {
729  //! Get the id for the next cation
730  icat = cationList_[k];
731  jNeut = fm_invert_ionForNeutral[icat];
732  if (jNeut != npos) {
733  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
734  AssertTrace(fmij != 0.0);
735  NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
736  }
737  }
738 
739  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
740  icat = passThroughList_[k];
741  jNeut = fm_invert_ionForNeutral[icat];
742  fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
743  NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
744  }
745 
746 #ifdef DEBUG_MODE
747  for (size_t k = 0; k < m_kk; k++) {
749  }
750  for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
751  for (size_t k = 0; k < m_kk; k++) {
752  fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
753  moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut];
754  }
755  }
756  for (size_t k = 0; k < m_kk; k++) {
757  if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
758  //! Check to see if we have in fact found the inverse.
759  if (anionList_[0] != k) {
760  throw CanteraError("", "neutral molecule calc error");
761  } else {
762  //! For the single anion case, we will allow some slippage
763  if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
764  throw CanteraError("", "neutral molecule calc error - anion");
765  }
766  }
767  }
768  }
769 #endif
770 
771  // Normalize the Neutral Molecule mole fractions
772  sum = 0.0;
773  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
774  sum += NeutralMolecMoleFractions_[k];
775  }
776  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
777  NeutralMolecMoleFractions_[k] /= sum;
778  }
779 
780  break;
781 
782  case cIonSolnType_SINGLECATION:
783 
784  throw CanteraError("eosType", "Unknown type");
785 
786  break;
787 
788  case cIonSolnType_MULTICATIONANION:
789 
790  throw CanteraError("eosType", "Unknown type");
791  break;
792 
793  default:
794 
795  throw CanteraError("eosType", "Unknown type");
796  break;
797 
798  }
799 }
800 //====================================================================================================================
801 // Calculate neutral molecule mole fractions
802 /*
803  * This routine calculates the neutral molecule mole
804  * fraction given the vector of ion mole fractions,
805  * i.e., the mole fractions from this ThermoPhase.
806  * Note, this routine basically assumes that there
807  * is charge neutrality. If there isn't, then it wouldn't
808  * make much sense.
809  *
810  * for the case of cIonSolnType_SINGLEANION, some slough
811  * in the charge neutrality is allowed. The cation number
812  * is followed, while the difference in charge neutrality
813  * is dumped into the anion mole number to fix the imbalance.
814  */
815 void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
816 {
817  doublereal fmij;
818  vector_fp y;
819  y.resize(numNeutralMoleculeSpecies_,0.0);
820  doublereal sumy, sumdy;
821 
822  //check sum dx = 0
823 
824  //! Zero the vector we are trying to find.
825  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
826  dy[k] = 0.0;
827  }
828 
829 
830  // bool fmSimple = true;
831 
832  switch (ionSolnType_) {
833 
834  case cIonSolnType_PASSTHROUGH:
835  for (size_t k = 0; k < m_kk; k++) {
836  dy[k] = dx[k];
837  }
838  break;
839 
840  case cIonSolnType_SINGLEANION:
841  for (size_t k = 0; k < cationList_.size(); k++) {
842  //! Get the id for the next cation
843  size_t icat = cationList_[k];
844  size_t jNeut = fm_invert_ionForNeutral[icat];
845  if (jNeut != npos) {
846  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
847  AssertTrace(fmij != 0.0);
848  dy[jNeut] += dx[icat] / fmij;
849  y[jNeut] += moleFractions_[icat] / fmij;
850  }
851  }
852 
853  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
854  size_t icat = passThroughList_[k];
855  size_t jNeut = fm_invert_ionForNeutral[icat];
856  fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
857  dy[jNeut] += dx[icat] / fmij;
858  y[jNeut] += moleFractions_[icat] / fmij;
859  }
860 #ifdef DEBUG_MODE_NOT
861  //check dy sum to zero
862  for (size_t k = 0; k < m_kk; k++) {
863  moleFractionsTmp_[k] = dx[k];
864  }
865  for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
866  for (size_t k = 0; k < m_kk; k++) {
867  fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
868  moleFractionsTmp_[k] -= fmij * dy[jNeut];
869  }
870  }
871  for (size_t k = 0; k < m_kk; k++) {
872  if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
873  //! Check to see if we have in fact found the inverse.
874  if (anionList_[0] != k) {
875  throw CanteraError("", "neutral molecule calc error");
876  } else {
877  //! For the single anion case, we will allow some slippage
878  if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
879  throw CanteraError("", "neutral molecule calc error - anion");
880  }
881  }
882  }
883  }
884 #endif
885  // Normalize the Neutral Molecule mole fractions
886  sumy = 0.0;
887  sumdy = 0.0;
888  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
889  sumy += y[k];
890  sumdy += dy[k];
891  }
892  for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
893  dy[k] = dy[k]/sumy - y[k]*sumdy/sumy/sumy;
894  }
895 
896  break;
897 
898  case cIonSolnType_SINGLECATION:
899 
900  throw CanteraError("eosType", "Unknown type");
901 
902  break;
903 
904  case cIonSolnType_MULTICATIONANION:
905 
906  throw CanteraError("eosType", "Unknown type");
907  break;
908 
909  default:
910 
911  throw CanteraError("eosType", "Unknown type");
912  break;
913 
914  }
915 }
916 
917 
918 void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y)
919 {
923 }
924 
925 void IonsFromNeutralVPSSTP::setMassFractions_NoNorm(const doublereal* const y)
926 {
930 }
931 
932 void IonsFromNeutralVPSSTP::setMoleFractions(const doublereal* const x)
933 {
937 }
938 
939 void IonsFromNeutralVPSSTP::setMoleFractions_NoNorm(const doublereal* const x)
940 {
944 }
945 
946 
947 void IonsFromNeutralVPSSTP::setConcentrations(const doublereal* const c)
948 {
952 }
953 
954 /*
955  * ------------ Partial Molar Properties of the Solution ------------
956  */
957 
958 
959 doublereal IonsFromNeutralVPSSTP::err(std::string msg) const
960 {
961  throw CanteraError("IonsFromNeutralVPSSTP","Base class method "
962  +msg+" called. Equation of state type: "+int2str(eosType()));
963  return 0;
964 }
965 /*
966  * Import, construct, and initialize a phase
967  * specification from an XML tree into the current object.
968  *
969  * This routine is a precursor to constructPhaseXML(XML_Node*)
970  * routine, which does most of the work.
971  *
972  * @param infile XML file containing the description of the
973  * phase
974  *
975  * @param id Optional parameter identifying the name of the
976  * phase. If none is given, the first XML
977  * phase element will be used.
978  */
979 void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::string id)
980 {
981 
982  if (inputFile.size() == 0) {
983  throw CanteraError("MargulesVPSSTP:constructPhaseFile",
984  "input file is null");
985  }
986  string path = findInputFile(inputFile);
987  std::ifstream fin(path.c_str());
988  if (!fin) {
989  throw CanteraError("MargulesVPSSTP:constructPhaseFile","could not open "
990  +path+" for reading.");
991  }
992  /*
993  * The phase object automatically constructs an XML object.
994  * Use this object to store information.
995  */
996  XML_Node& phaseNode_XML = xml();
997  XML_Node* fxml = new XML_Node();
998  fxml->build(fin);
999  XML_Node* fxml_phase = findXMLPhase(fxml, id);
1000  if (!fxml_phase) {
1001  throw CanteraError("MargulesVPSSTP:constructPhaseFile",
1002  "ERROR: Can not find phase named " +
1003  id + " in file named " + inputFile);
1004  }
1005  fxml_phase->copy(&phaseNode_XML);
1006  constructPhaseXML(*fxml_phase, id);
1007  delete fxml;
1008 }
1009 
1010 /*
1011  * Import, construct, and initialize a HMWSoln phase
1012  * specification from an XML tree into the current object.
1013  *
1014  * Most of the work is carried out by the cantera base
1015  * routine, importPhase(). That routine imports all of the
1016  * species and element data, including the standard states
1017  * of the species.
1018  *
1019  * Then, In this routine, we read the information
1020  * particular to the specification of the activity
1021  * coefficient model for the Pitzer parameterization.
1022  *
1023  * We also read information about the molar volumes of the
1024  * standard states if present in the XML file.
1025  *
1026  * @param phaseNode This object must be the phase node of a
1027  * complete XML tree
1028  * description of the phase, including all of the
1029  * species data. In other words while "phase" must
1030  * point to an XML phase object, it must have
1031  * sibling nodes "speciesData" that describe
1032  * the species in the phase.
1033  * @param id ID of the phase. If nonnull, a check is done
1034  * to see if phaseNode is pointing to the phase
1035  * with the correct id.
1036  */
1037 void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id)
1038 {
1039  string stemp;
1040  if (id.size() > 0) {
1041  string idp = phaseNode.id();
1042  if (idp != id) {
1043  throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
1044  "phasenode and Id are incompatible");
1045  }
1046  }
1047 
1048  /*
1049  * Find the Thermo XML node
1050  */
1051  if (!phaseNode.hasChild("thermo")) {
1052  throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
1053  "no thermo XML node");
1054  }
1055  XML_Node& thermoNode = phaseNode.child("thermo");
1056 
1057 
1058 
1059  /*
1060  * Make sure that the thermo model is IonsFromNeutralMolecule
1061  */
1062  stemp = thermoNode.attrib("model");
1063  string formString = lowercase(stemp);
1064  if (formString != "ionsfromneutralmolecule") {
1065  throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
1066  "model name isn't IonsFromNeutralMolecule: " + formString);
1067  }
1068 
1069  /*
1070  * Find the Neutral Molecule Phase
1071  */
1072  if (!thermoNode.hasChild("neutralMoleculePhase")) {
1073  throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
1074  "no neutralMoleculePhase XML node");
1075  }
1076  XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");
1077 
1078  string nsource = neutralMoleculeNode["datasrc"];
1079  XML_Node* neut_ptr = get_XML_Node(nsource, 0);
1080  if (!neut_ptr) {
1081  throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
1082  "neut_ptr = 0");
1083  }
1084 
1085  /*
1086  * Create the neutralMolecule ThermoPhase if we haven't already
1087  */
1088  if (!neutralMoleculePhase_) {
1089  neutralMoleculePhase_ = newPhase(*neut_ptr);
1090  }
1091 
1092  /*
1093  * Call the Cantera importPhase() function. This will import
1094  * all of the species into the phase. This will also handle
1095  * all of the solvent and solute standard states
1096  */
1097  bool m_ok = importPhase(phaseNode, this);
1098  if (!m_ok) {
1099  throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
1100  "importPhase failed ");
1101  }
1102 
1103 }
1104 
1105 //====================================================================================================================
1106 /*
1107  * @internal Initialize. This method is provided to allow
1108  * subclasses to perform any initialization required after all
1109  * species have been added. For example, it might be used to
1110  * resize internal work arrays that must have an entry for
1111  * each species. The base class implementation does nothing,
1112  * and subclasses that do not require initialization do not
1113  * need to overload this method. When importing a CTML phase
1114  * description, this method is called just prior to returning
1115  * from function importPhase.
1116  *
1117  * @see importCTML.cpp
1118  */
1120 {
1121  initLengths();
1123 }
1124 //====================================================================================================================
1125 
1126 // Initialize lengths of local variables after all species have
1127 // been identified.
1129 {
1130  m_kk = nSpecies();
1132  moleFractions_.resize(m_kk);
1134  fm_invert_ionForNeutral.resize(m_kk);
1136  cationList_.resize(m_kk);
1137  anionList_.resize(m_kk);
1138  passThroughList_.resize(m_kk);
1139  moleFractionsTmp_.resize(m_kk);
1146 }
1147 //====================================================================================================================
1148 //! Return the factor overlap
1149 /*!
1150  * @param elnamesVN
1151  * @param elemVectorN
1152  * @param nElementsN
1153  * @param elnamesVI
1154  * @param elemVectorI
1155  * @param nElementsI
1156  *
1157  */
1158 static double factorOverlap(const std::vector<std::string>& elnamesVN ,
1159  const std::vector<double>& elemVectorN,
1160  const size_t nElementsN,
1161  const std::vector<std::string>& elnamesVI ,
1162  const std::vector<double>& elemVectorI,
1163  const size_t nElementsI)
1164 {
1165  double fMax = 1.0E100;
1166  for (size_t mi = 0; mi < nElementsI; mi++) {
1167  if (elnamesVI[mi] != "E") {
1168  if (elemVectorI[mi] > 1.0E-13) {
1169  double eiNum = elemVectorI[mi];
1170  for (size_t mn = 0; mn < nElementsN; mn++) {
1171  if (elnamesVI[mi] == elnamesVN[mn]) {
1172  if (elemVectorN[mn] <= 1.0E-13) {
1173  return 0.0;
1174  }
1175  fMax = std::min(fMax, elemVectorN[mn]/eiNum);
1176  }
1177  }
1178  }
1179  }
1180  }
1181  return fMax;
1182 }
1183 //====================================================================================================================
1184 /*
1185  * initThermoXML() (virtual from ThermoPhase)
1186  * Import and initialize a ThermoPhase object
1187  *
1188  * @param phaseNode This object must be the phase node of a
1189  * complete XML tree
1190  * description of the phase, including all of the
1191  * species data. In other words while "phase" must
1192  * point to an XML phase object, it must have
1193  * sibling nodes "speciesData" that describe
1194  * the species in the phase.
1195  * @param id ID of the phase. If nonnull, a check is done
1196  * to see if phaseNode is pointing to the phase
1197  * with the correct id.
1198  */
1199 void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, std::string id)
1200 {
1201  size_t k;
1202  /*
1203  * variables that need to be populated
1204  *
1205  * cationList_
1206  * numCationSpecies_;
1207  */
1208 
1209  numCationSpecies_ = 0;
1210  cationList_.clear();
1211  for (k = 0; k < m_kk; k++) {
1212  if (charge(k) > 0) {
1213  cationList_.push_back(k);
1215  }
1216  }
1217 
1218  numAnionSpecies_ = 0;
1219  anionList_.clear();
1220  for (k = 0; k < m_kk; k++) {
1221  if (charge(k) < 0) {
1222  anionList_.push_back(k);
1223  numAnionSpecies_++;
1224  }
1225  }
1226 
1228  passThroughList_.clear();
1229  for (k = 0; k < m_kk; k++) {
1230  if (charge(k) == 0) {
1231  passThroughList_.push_back(k);
1233  }
1234  }
1235 
1236  PDSS_IonsFromNeutral* speciesSS = 0;
1238  for (k = 0; k < m_kk; k++) {
1239  speciesSS = dynamic_cast<PDSS_IonsFromNeutral*>(providePDSS(k));
1240  if (!speciesSS) {
1241  throw CanteraError("initThermoXML", "Dynamic cast failed");
1242  }
1243  if (speciesSS->specialSpecies_ == 1) {
1245  }
1246  if (speciesSS->specialSpecies_ == 2) {
1248  }
1249  }
1250 
1251 
1252  size_t nElementsN = neutralMoleculePhase_->nElements();
1253  const std::vector<std::string>& elnamesVN = neutralMoleculePhase_->elementNames();
1254  std::vector<double> elemVectorN(nElementsN);
1255  std::vector<double> elemVectorN_orig(nElementsN);
1256 
1257  size_t nElementsI = nElements();
1258  const std::vector<std::string>& elnamesVI = elementNames();
1259  std::vector<double> elemVectorI(nElementsI);
1260 
1261  vector<doublereal> fm_tmp(m_kk);
1262  for (size_t k = 0; k < m_kk; k++) {
1264  }
1265  /* for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
1266  fm_invert_ionForNeutral[jNeut] = -1;
1267  }*/
1268  for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
1269  for (size_t m = 0; m < nElementsN; m++) {
1270  elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m);
1271  }
1272  elemVectorN_orig = elemVectorN;
1273  fm_tmp.assign(m_kk, 0.0);
1274 
1275  for (size_t m = 0; m < nElementsI; m++) {
1276  elemVectorI[m] = nAtoms(indexSpecialSpecies_, m);
1277  }
1278  double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
1279  elnamesVI ,elemVectorI, nElementsI);
1280  if (fac > 0.0) {
1281  for (size_t m = 0; m < nElementsN; m++) {
1282  std::string mName = elnamesVN[m];
1283  for (size_t mi = 0; mi < nElementsI; mi++) {
1284  std::string eName = elnamesVI[mi];
1285  if (mName == eName) {
1286  elemVectorN[m] -= fac * elemVectorI[mi];
1287  }
1288 
1289  }
1290  }
1291  }
1292  fm_neutralMolec_ions_[indexSpecialSpecies_ + jNeut * m_kk ] += fac;
1293 
1294 
1295  for (k = 0; k < m_kk; k++) {
1296  for (size_t m = 0; m < nElementsI; m++) {
1297  elemVectorI[m] = nAtoms(k, m);
1298  }
1299  double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
1300  elnamesVI ,elemVectorI, nElementsI);
1301  if (fac > 0.0) {
1302  for (size_t m = 0; m < nElementsN; m++) {
1303  std::string mName = elnamesVN[m];
1304  for (size_t mi = 0; mi < nElementsI; mi++) {
1305  std::string eName = elnamesVI[mi];
1306  if (mName == eName) {
1307  elemVectorN[m] -= fac * elemVectorI[mi];
1308  }
1309 
1310  }
1311  }
1312  bool notTaken = true;
1313  for (size_t iNeut = 0; iNeut < jNeut; iNeut++) {
1314  if (fm_invert_ionForNeutral[k] == iNeut) {
1315  notTaken = false;
1316  }
1317  }
1318  if (notTaken) {
1319  fm_invert_ionForNeutral[k] = jNeut;
1320  } else {
1321  throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
1322  "Simple formula matrix generation failed, one cation is shared between two salts");
1323  }
1324  }
1325  fm_neutralMolec_ions_[k + jNeut * m_kk] += fac;
1326  }
1327 
1328  // Ok check the work
1329  for (size_t m = 0; m < nElementsN; m++) {
1330  if (fabs(elemVectorN[m]) > 1.0E-13) {
1331  throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
1332  "Simple formula matrix generation failed");
1333  }
1334  }
1335 
1336 
1337  }
1338  /*
1339  * This includes the setStateFromXML calls
1340  */
1341  GibbsExcessVPSSTP::initThermoXML(phaseNode, id);
1342 
1343  /*
1344  * There is one extra step here. We assure ourselves that we
1345  * have charge conservation.
1346  */
1347 }
1348 //====================================================================================================================
1349 // Update the activity coefficients
1350 /*
1351  * This function will be called to update the internally stored
1352  * natural logarithm of the activity coefficients
1353  *
1354  */
1356 {
1357  size_t icat, jNeut;
1358  doublereal fmij;
1359  /*
1360  * Get the activity coefficiens of the neutral molecules
1361  */
1363 
1364  switch (ionSolnType_) {
1365  case cIonSolnType_PASSTHROUGH:
1366  break;
1367  case cIonSolnType_SINGLEANION:
1368 
1369  // Do the cation list
1370  for (size_t k = 0; k < cationList_.size(); k++) {
1371  //! Get the id for the next cation
1372  icat = cationList_[k];
1373  jNeut = fm_invert_ionForNeutral[icat];
1374  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
1375  lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut] / fmij;
1376  }
1377 
1378  // Do the anion list
1379  icat = anionList_[0];
1380  jNeut = fm_invert_ionForNeutral[icat];
1381  lnActCoeff_Scaled_[icat]= 0.0;
1382 
1383  // Do the list of neutral molecules
1384  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
1385  icat = passThroughList_[k];
1386  jNeut = fm_invert_ionForNeutral[icat];
1388  }
1389  break;
1390 
1391  case cIonSolnType_SINGLECATION:
1392  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
1393  break;
1394  case cIonSolnType_MULTICATIONANION:
1395  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
1396  break;
1397  default:
1398  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
1399  break;
1400  }
1401 
1402 }
1403 //====================================================================================================================
1404 // Get the change in activity coefficients w.r.t. change in state (temp, mole fraction, etc.) along
1405 // a line in parameter space or along a line in physical space
1406 /*
1407  *
1408  * @param dTds Input of temperature change along the path
1409  * @param dXds Input vector of changes in mole fraction along the path. length = m_kk
1410  * Along the path length it must be the case that the mole fractions sum to one.
1411  * @param dlnActCoeffds Output vector of the directional derivatives of the
1412  * log Activity Coefficients along the path. length = m_kk
1413  */
1414 void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
1415  doublereal* dlnActCoeffds) const
1416 {
1417  size_t icat, jNeut;
1418  doublereal fmij;
1419  /*
1420  * Get the activity coefficients of the neutral molecules
1421  */
1422  GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
1423  if (!geThermo) {
1424  for (size_t k = 0; k < m_kk; k++) {
1425  dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
1426  }
1427  return;
1428  }
1429 
1430  size_t numNeutMolSpec = geThermo->nSpecies();
1431  vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
1432  vector_fp dX_NeutralMolecule(numNeutMolSpec);
1433 
1434 
1435  getNeutralMoleculeMoleGrads(DATA_PTR(dXds),DATA_PTR(dX_NeutralMolecule));
1436 
1437  // All mole fractions returned to normal
1438 
1439  geThermo->getdlnActCoeffds(dTds, DATA_PTR(dX_NeutralMolecule), DATA_PTR(dlnActCoeff_NeutralMolecule));
1440 
1441  switch (ionSolnType_) {
1442  case cIonSolnType_PASSTHROUGH:
1443  break;
1444  case cIonSolnType_SINGLEANION:
1445 
1446  // Do the cation list
1447  for (size_t k = 0; k < cationList_.size(); k++) {
1448  //! Get the id for the next cation
1449  icat = cationList_[k];
1450  jNeut = fm_invert_ionForNeutral[icat];
1451  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
1452  dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule[jNeut]/fmij;
1453  }
1454 
1455  // Do the anion list
1456  icat = anionList_[0];
1457  jNeut = fm_invert_ionForNeutral[icat];
1458  dlnActCoeffds[icat]= 0.0;
1459 
1460  // Do the list of neutral molecules
1461  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
1462  icat = passThroughList_[k];
1463  jNeut = fm_invert_ionForNeutral[icat];
1464  dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule[jNeut];
1465  }
1466  break;
1467 
1468  case cIonSolnType_SINGLECATION:
1469  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
1470  break;
1471  case cIonSolnType_MULTICATIONANION:
1472  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
1473  break;
1474  default:
1475  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
1476  break;
1477  }
1478 
1479 }
1480 //====================================================================================================================
1481 // Update the temperature derivative of the ln activity coefficients
1482 /*
1483  * This function will be called to update the internally stored
1484  * temperature derivative of the natural logarithm of the activity coefficients
1485  */
1487 {
1488  size_t icat, jNeut;
1489  doublereal fmij;
1490  /*
1491  * Get the activity coefficients of the neutral molecules
1492  */
1493  GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
1494  if (!geThermo) {
1495  dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
1496  return;
1497  }
1498 
1500 
1501  switch (ionSolnType_) {
1502  case cIonSolnType_PASSTHROUGH:
1503  break;
1504  case cIonSolnType_SINGLEANION:
1505 
1506  // Do the cation list
1507  for (size_t k = 0; k < cationList_.size(); k++) {
1508  //! Get the id for the next cation
1509  icat = cationList_[k];
1510  jNeut = fm_invert_ionForNeutral[icat];
1511  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
1513  }
1514 
1515  // Do the anion list
1516  icat = anionList_[0];
1517  jNeut = fm_invert_ionForNeutral[icat];
1518  dlnActCoeffdT_Scaled_[icat]= 0.0;
1519 
1520  // Do the list of neutral molecules
1521  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
1522  icat = passThroughList_[k];
1523  jNeut = fm_invert_ionForNeutral[icat];
1525  }
1526  break;
1527 
1528  case cIonSolnType_SINGLECATION:
1529  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
1530  break;
1531  case cIonSolnType_MULTICATIONANION:
1532  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
1533  break;
1534  default:
1535  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
1536  break;
1537  }
1538 
1539 }
1540 //====================================================================================================================
1541 /*
1542  * This function will be called to update the internally stored
1543  * temperature derivative of the natural logarithm of the activity coefficients
1544  */
1546 {
1547  size_t icat, jNeut;
1548  doublereal fmij;
1549  /*
1550  * Get the activity coefficients of the neutral molecules
1551  */
1552  GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
1553  if (!geThermo) {
1554  dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
1555  return;
1556  }
1557 
1559 
1560  switch (ionSolnType_) {
1561  case cIonSolnType_PASSTHROUGH:
1562  break;
1563  case cIonSolnType_SINGLEANION:
1564 
1565  // Do the cation list
1566  for (size_t k = 0; k < cationList_.size(); k++) {
1567  //! Get the id for the next cation
1568  icat = cationList_[k];
1569  jNeut = fm_invert_ionForNeutral[icat];
1570  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
1572  }
1573 
1574  // Do the anion list
1575  icat = anionList_[0];
1576  jNeut = fm_invert_ionForNeutral[icat];
1577  dlnActCoeffdlnX_diag_[icat]= 0.0;
1578 
1579  // Do the list of neutral molecules
1580  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
1581  icat = passThroughList_[k];
1582  jNeut = fm_invert_ionForNeutral[icat];
1584  }
1585  break;
1586 
1587  case cIonSolnType_SINGLECATION:
1588  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
1589  break;
1590  case cIonSolnType_MULTICATIONANION:
1591  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
1592  break;
1593  default:
1594  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
1595  break;
1596  }
1597 
1598 }
1599 //====================================================================================================================
1600 /*
1601  * This function will be called to update the internally stored
1602  * temperature derivative of the natural logarithm of the activity coefficients
1603  */
1605 {
1606  size_t icat, jNeut;
1607  doublereal fmij;
1608  /*
1609  * Get the activity coefficients of the neutral molecules
1610  */
1611  GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
1612  if (!geThermo) {
1613  dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
1614  return;
1615  }
1616 
1618 
1619  switch (ionSolnType_) {
1620  case cIonSolnType_PASSTHROUGH:
1621  break;
1622  case cIonSolnType_SINGLEANION:
1623 
1624  // Do the cation list
1625  for (size_t k = 0; k < cationList_.size(); k++) {
1626  //! Get the id for the next cation
1627  icat = cationList_[k];
1628  jNeut = fm_invert_ionForNeutral[icat];
1629  fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
1631  }
1632 
1633  // Do the anion list
1634  icat = anionList_[0];
1635  jNeut = fm_invert_ionForNeutral[icat];
1636  dlnActCoeffdlnN_diag_[icat]= 0.0;
1637 
1638  // Do the list of neutral molecules
1639  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
1640  icat = passThroughList_[k];
1641  jNeut = fm_invert_ionForNeutral[icat];
1643  }
1644  break;
1645 
1646  case cIonSolnType_SINGLECATION:
1647  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
1648  break;
1649  case cIonSolnType_MULTICATIONANION:
1650  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
1651  break;
1652  default:
1653  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
1654  break;
1655  }
1656 
1657 }
1658 //====================================================================================================================
1659 // Update the derivative of the log of the activity coefficients
1660 // wrt log(number of moles) - diagonal components
1661 /*
1662  * This function will be called to update the internally stored
1663  * derivative of the natural logarithm of the activity coefficients
1664  * wrt logarithm of the number of moles of given species.
1665  */
1667 {
1668  size_t kcat, kNeut, mcat, mNeut;
1669  doublereal fmij, mfmij;
1671  /*
1672  * Get the activity coefficients of the neutral molecules
1673  */
1674  GibbsExcessVPSSTP* geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
1675  if (!geThermo) {
1676  throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed");
1677  }
1678  size_t nsp_ge = geThermo->nSpecies();
1679  geThermo->getdlnActCoeffdlnN(nsp_ge, &(dlnActCoeffdlnN_NeutralMolecule_(0,0)));
1680 
1681  switch (ionSolnType_) {
1682  case cIonSolnType_PASSTHROUGH:
1683  break;
1684  case cIonSolnType_SINGLEANION:
1685 
1686  // Do the cation list
1687  for (size_t k = 0; k < cationList_.size(); k++) {
1688  for (size_t m = 0; m < cationList_.size(); m++) {
1689  kcat = cationList_[k];
1690 
1691  kNeut = fm_invert_ionForNeutral[kcat];
1692  fmij = fm_neutralMolec_ions_[kcat + kNeut * m_kk];
1694 
1695  mcat = cationList_[m];
1696  mNeut = fm_invert_ionForNeutral[mcat];
1697  mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
1698 
1699  dlnActCoeffdlnN_(kcat,mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij;
1700 
1701  }
1702  for (size_t m = 0; m < numPassThroughSpecies_; m++) {
1703  mcat = passThroughList_[m];
1704  mNeut = fm_invert_ionForNeutral[mcat];
1705  dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut) / fmij;
1706  }
1707  }
1708 
1709  // Do the anion list -> anion activity coefficient is one
1710  kcat = anionList_[0];
1711  kNeut = fm_invert_ionForNeutral[kcat];
1712  for (size_t k = 0; k < m_kk; k++) {
1713  dlnActCoeffdlnN_(kcat, k) = 0.0;
1714  dlnActCoeffdlnN_(k, kcat) = 0.0;
1715  }
1716 
1717  // Do the list of neutral molecules
1718  for (size_t k = 0; k < numPassThroughSpecies_; k++) {
1719  kcat = passThroughList_[k];
1720  kNeut = fm_invert_ionForNeutral[kcat];
1722 
1723  for (size_t m = 0; m < m_kk; m++) {
1724  mcat = passThroughList_[m];
1725  mNeut = fm_invert_ionForNeutral[mcat];
1726  dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut);
1727  }
1728 
1729 
1730  for (size_t m = 0; m < cationList_.size(); m++) {
1731  mcat = cationList_[m];
1732  mNeut = fm_invert_ionForNeutral[mcat];
1733  mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
1734  dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut);
1735  }
1736 
1737  }
1738  break;
1739 
1740  case cIonSolnType_SINGLECATION:
1741  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
1742  break;
1743  case cIonSolnType_MULTICATIONANION:
1744  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
1745  break;
1746  default:
1747  throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
1748  break;
1749  }
1750 }
1751 //====================================================================================================================
1752 }
1753 //======================================================================================================================