Cantera  3.1.0a1
HMWSoln.cpp
Go to the documentation of this file.
1 /**
2  * @file HMWSoln.cpp
3  * Definitions for the HMWSoln ThermoPhase object, which
4  * models concentrated electrolyte solutions
5  * (see @ref thermoprops and @link Cantera::HMWSoln HMWSoln @endlink) .
6  *
7  * Class HMWSoln represents a concentrated liquid electrolyte phase which obeys
8  * the Pitzer formulation for nonideality using molality-based standard states.
9  *
10  * This version of the code was modified to have the binary Beta2 Pitzer
11  * parameter consistent with the temperature expansions used for Beta0,
12  * Beta1, and Cphi.(CFJC, SNL)
13  */
14 
15 // This file is part of Cantera. See License.txt in the top-level directory or
16 // at https://cantera.org/license.txt for license and copyright information.
17 
18 #include "cantera/thermo/HMWSoln.h"
23 
24 namespace Cantera
25 {
26 
27 namespace {
28 double A_Debye_default = 1.172576; // units = sqrt(kg/gmol)
29 double maxIionicStrength_default = 100.0;
30 double crop_ln_gamma_o_min_default = -6.0;
31 double crop_ln_gamma_o_max_default = 3.0;
32 double crop_ln_gamma_k_min_default = -5.0;
33 double crop_ln_gamma_k_max_default = 15.0;
34 }
35 
36 HMWSoln::~HMWSoln()
37 {
38  // Defined in .cpp to limit dependence on WaterProps.h
39 }
40 
41 HMWSoln::HMWSoln(const string& inputFile, const string& id_) :
42  m_maxIionicStrength(maxIionicStrength_default),
43  CROP_ln_gamma_o_min(crop_ln_gamma_o_min_default),
44  CROP_ln_gamma_o_max(crop_ln_gamma_o_max_default),
45  CROP_ln_gamma_k_min(crop_ln_gamma_k_min_default),
46  CROP_ln_gamma_k_max(crop_ln_gamma_k_max_default),
47  m_last_is(-1.0)
48 {
49  initThermoFile(inputFile, id_);
50 }
51 
52 // -------- Molar Thermodynamic Properties of the Solution ---------------
53 
54 double HMWSoln::enthalpy_mole() const
55 {
57  return mean_X(m_tmpV);
58 }
59 
61 {
63  double hbar = mean_X(m_tmpV);
65  for (size_t k = 0; k < m_kk; k++) {
66  m_gamma_tmp[k] *= RT();
67  }
68  double h0bar = mean_X(m_gamma_tmp);
69  return hbar - h0bar;
70 }
71 
73 {
74  double L = relative_enthalpy();
75  getMoleFractions(m_tmpV.data());
76  double xanion = 0.0;
77  size_t kcation = npos;
78  double xcation = 0.0;
79  size_t kanion = npos;
80  for (size_t k = 0; k < m_kk; k++) {
81  if (charge(k) > 0.0) {
82  if (m_tmpV[k] > xanion) {
83  xanion = m_tmpV[k];
84  kanion = k;
85  }
86  } else if (charge(k) < 0.0) {
87  if (m_tmpV[k] > xcation) {
88  xcation = m_tmpV[k];
89  kcation = k;
90  }
91  }
92  }
93  if (kcation == npos || kanion == npos) {
94  return L;
95  }
96  double xuse = xcation;
97  double factor = 1;
98  if (xanion < xcation) {
99  xuse = xanion;
100  if (charge(kcation) != 1.0) {
101  factor = charge(kcation);
102  }
103  } else {
104  if (charge(kanion) != 1.0) {
105  factor = charge(kanion);
106  }
107  }
108  xuse = xuse / factor;
109  return L / xuse;
110 }
111 
112 double HMWSoln::entropy_mole() const
113 {
115  return mean_X(m_tmpV);
116 }
117 
118 double HMWSoln::gibbs_mole() const
119 {
120  getChemPotentials(m_tmpV.data());
121  return mean_X(m_tmpV);
122 }
123 
124 double HMWSoln::cp_mole() const
125 {
126  getPartialMolarCp(m_tmpV.data());
127  return mean_X(m_tmpV);
128 }
129 
130 double HMWSoln::cv_mole() const
131 {
132  double kappa_t = isothermalCompressibility();
133  double beta = thermalExpansionCoeff();
134  double cp = cp_mole();
135  double tt = temperature();
136  double molarV = molarVolume();
137  return cp - beta * beta * tt * molarV / kappa_t;
138 }
139 
140 // ------- Mechanical Equation of State Properties ------------------------
141 
143 {
144  static const int cacheId = m_cache.getId();
145  CachedScalar cached = m_cache.getScalar(cacheId);
146  if(cached.validate(temperature(), pressure(), stateMFNumber())) {
147  return;
148  }
149 
150  // Calculate all of the other standard volumes. Note these are constant for
151  // now
153  double dd = meanMolecularWeight() / mean_X(m_tmpV);
155 }
156 
157 // ------- Activities and Activity Concentrations
158 
160 {
161  double cs_solvent = standardConcentration();
162  getActivities(c);
163  c[0] *= cs_solvent;
164  if (m_kk > 1) {
165  double cs_solute = standardConcentration(1);
166  for (size_t k = 1; k < m_kk; k++) {
167  c[k] *= cs_solute;
168  }
169  }
170 }
171 
172 double HMWSoln::standardConcentration(size_t k) const
173 {
174  getStandardVolumes(m_tmpV.data());
175  double mvSolvent = m_tmpV[0];
176  if (k > 0) {
177  return m_Mnaught / mvSolvent;
178  }
179  return 1.0 / mvSolvent;
180 }
181 
182 void HMWSoln::getActivities(double* ac) const
183 {
185 
186  // Update the molality array, m_molalities(). This requires an update due to
187  // mole fractions
189 
190  // Now calculate the array of activities.
191  for (size_t k = 1; k < m_kk; k++) {
192  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
193  }
194  double xmolSolvent = moleFraction(0);
195  ac[0] = exp(m_lnActCoeffMolal_Scaled[0]) * xmolSolvent;
196 }
197 
199 {
201  A_Debye_TP(-1.0, -1.0);
203  std::copy(m_lnActCoeffMolal_Unscaled.begin(), m_lnActCoeffMolal_Unscaled.end(), acMolality);
204  for (size_t k = 0; k < m_kk; k++) {
205  acMolality[k] = exp(acMolality[k]);
206  }
207 }
208 
209 // ------ Partial Molar Properties of the Solution -----------------
210 
211 void HMWSoln::getChemPotentials(double* mu) const
212 {
213  double xx;
214 
215  // First get the standard chemical potentials in molar form. This requires
216  // updates of standard state as a function of T and P
218 
219  // Update the activity coefficients. This also updates the internal molality
220  // array.
222  double xmolSolvent = moleFraction(0);
223  for (size_t k = 1; k < m_kk; k++) {
224  xx = std::max(m_molalities[k], SmallNumber);
225  mu[k] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
226  }
227  xx = std::max(xmolSolvent, SmallNumber);
228  mu[0] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[0]);
229 }
230 
231 void HMWSoln::getPartialMolarEnthalpies(double* hbar) const
232 {
233  // Get the nondimensional standard state enthalpies
234  getEnthalpy_RT(hbar);
235 
236  // dimensionalize it.
237  for (size_t k = 0; k < m_kk; k++) {
238  hbar[k] *= RT();
239  }
240 
241  // Update the activity coefficients, This also update the internally stored
242  // molalities.
245  for (size_t k = 0; k < m_kk; k++) {
246  hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT_Scaled[k];
247  }
248 }
249 
250 void HMWSoln::getPartialMolarEntropies(double* sbar) const
251 {
252  // Get the standard state entropies at the temperature and pressure of the
253  // solution.
254  getEntropy_R(sbar);
255 
256  // Dimensionalize the entropies
257  for (size_t k = 0; k < m_kk; k++) {
258  sbar[k] *= GasConstant;
259  }
260 
261  // Update the activity coefficients, This also update the internally stored
262  // molalities.
264 
265  // First we will add in the obvious dependence on the T term out front of
266  // the log activity term
267  double mm;
268  for (size_t k = 1; k < m_kk; k++) {
269  mm = std::max(SmallNumber, m_molalities[k]);
270  sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
271  }
272  double xmolSolvent = moleFraction(0);
273  mm = std::max(SmallNumber, xmolSolvent);
274  sbar[0] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[0]);
275 
276  // Check to see whether activity coefficients are temperature dependent. If
277  // they are, then calculate the their temperature derivatives and add them
278  // into the result.
280  for (size_t k = 0; k < m_kk; k++) {
281  sbar[k] -= RT() * m_dlnActCoeffMolaldT_Scaled[k];
282  }
283 }
284 
285 void HMWSoln::getPartialMolarVolumes(double* vbar) const
286 {
287  // Get the standard state values in m^3 kmol-1
288  getStandardVolumes(vbar);
289 
290  // Update the derivatives wrt the activity coefficients.
293  for (size_t k = 0; k < m_kk; k++) {
294  vbar[k] += RT() * m_dlnActCoeffMolaldP_Scaled[k];
295  }
296 }
297 
298 void HMWSoln::getPartialMolarCp(double* cpbar) const
299 {
300  getCp_R(cpbar);
301  for (size_t k = 0; k < m_kk; k++) {
302  cpbar[k] *= GasConstant;
303  }
304 
305  // Update the activity coefficients, This also update the internally stored
306  // molalities.
310  for (size_t k = 0; k < m_kk; k++) {
311  cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT_Scaled[k] +
313  }
314 }
315 
316 // -------------- Utilities -------------------------------
317 
318 double HMWSoln::satPressure(double t) {
319  double p_old = pressure();
320  double t_old = temperature();
321  double pres = m_waterSS->satPressure(t);
322 
323  // Set the underlying object back to its original state.
324  m_waterSS->setState_TP(t_old, p_old);
325  return pres;
326 }
327 
328 static void check_nParams(const string& method, size_t nParams, size_t m_formPitzerTemp)
329 {
330  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT && nParams != 1) {
331  throw CanteraError(method, "'constant' temperature model requires one"
332  " coefficient for each of parameter, but {} were given", nParams);
333  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR && nParams != 2) {
334  throw CanteraError(method, "'linear' temperature model requires two"
335  " coefficients for each parameter, but {} were given", nParams);
336  }
337  if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1 && nParams != 5) {
338  throw CanteraError(method, "'complex' temperature model requires five"
339  " coefficients for each parameter, but {} were given", nParams);
340  }
341 }
342 
343 void HMWSoln::setBinarySalt(const string& sp1, const string& sp2,
344  size_t nParams, double* beta0, double* beta1, double* beta2,
345  double* Cphi, double alpha1, double alpha2)
346 {
347  size_t k1 = speciesIndex(sp1);
348  size_t k2 = speciesIndex(sp2);
349  if (k1 == npos) {
350  throw CanteraError("HMWSoln::setBinarySalt", "Species '{}' not found", sp1);
351  } else if (k2 == npos) {
352  throw CanteraError("HMWSoln::setBinarySalt", "Species '{}' not found", sp2);
353  }
354  if (charge(k1) < 0 && charge(k2) > 0) {
355  std::swap(k1, k2);
356  } else if (charge(k1) * charge(k2) >= 0) {
357  throw CanteraError("HMWSoln::setBinarySalt", "Species '{}' and '{}' "
358  "do not have opposite charges ({}, {})", sp1, sp2,
359  charge(k1), charge(k2));
360  }
361  check_nParams("HMWSoln::setBinarySalt", nParams, m_formPitzerTemp);
362 
363  size_t c = m_CounterIJ[k1 * m_kk + k2];
364  m_Beta0MX_ij[c] = beta0[0];
365  m_Beta1MX_ij[c] = beta1[0];
366  m_Beta2MX_ij[c] = beta2[0];
367  m_CphiMX_ij[c] = Cphi[0];
368  for (size_t n = 0; n < nParams; n++) {
369  m_Beta0MX_ij_coeff(n, c) = beta0[n];
370  m_Beta1MX_ij_coeff(n, c) = beta1[n];
371  m_Beta2MX_ij_coeff(n, c) = beta2[n];
372  m_CphiMX_ij_coeff(n, c) = Cphi[n];
373  }
374  m_Alpha1MX_ij[c] = alpha1;
375  m_Alpha2MX_ij[c] = alpha2;
376 }
377 
378 void HMWSoln::setTheta(const string& sp1, const string& sp2,
379  size_t nParams, double* theta)
380 {
381  size_t k1 = speciesIndex(sp1);
382  size_t k2 = speciesIndex(sp2);
383  if (k1 == npos) {
384  throw CanteraError("HMWSoln::setTheta", "Species '{}' not found", sp1);
385  } else if (k2 == npos) {
386  throw CanteraError("HMWSoln::setTheta", "Species '{}' not found", sp2);
387  }
388  if (charge(k1) * charge(k2) <= 0) {
389  throw CanteraError("HMWSoln::setTheta", "Species '{}' and '{}' "
390  "should both have the same (non-zero) charge ({}, {})", sp1, sp2,
391  charge(k1), charge(k2));
392  }
393  check_nParams("HMWSoln::setTheta", nParams, m_formPitzerTemp);
394  size_t c = m_CounterIJ[k1 * m_kk + k2];
395  m_Theta_ij[c] = theta[0];
396  for (size_t n = 0; n < nParams; n++) {
397  m_Theta_ij_coeff(n, c) = theta[n];
398  }
399 }
400 
401 void HMWSoln::setPsi(const string& sp1, const string& sp2,
402  const string& sp3, size_t nParams, double* psi)
403 {
404  size_t k1 = speciesIndex(sp1);
405  size_t k2 = speciesIndex(sp2);
406  size_t k3 = speciesIndex(sp3);
407  if (k1 == npos) {
408  throw CanteraError("HMWSoln::setPsi", "Species '{}' not found", sp1);
409  } else if (k2 == npos) {
410  throw CanteraError("HMWSoln::setPsi", "Species '{}' not found", sp2);
411  } else if (k3 == npos) {
412  throw CanteraError("HMWSoln::setPsi", "Species '{}' not found", sp3);
413  }
414 
415  if (!charge(k1) || !charge(k2) || !charge(k3) ||
416  std::abs(sign(charge(k1) + sign(charge(k2)) + sign(charge(k3)))) != 1) {
417  throw CanteraError("HMWSoln::setPsi", "All species must be ions and"
418  " must include at least one cation and one anion, but given species"
419  " (charges) were: {} ({}), {} ({}), and {} ({}).",
420  sp1, charge(k1), sp2, charge(k2), sp3, charge(k3));
421  }
422  check_nParams("HMWSoln::setPsi", nParams, m_formPitzerTemp);
423  auto cc = {k1*m_kk*m_kk + k2*m_kk + k3,
424  k1*m_kk*m_kk + k3*m_kk + k2,
425  k2*m_kk*m_kk + k1*m_kk + k3,
426  k2*m_kk*m_kk + k3*m_kk + k1,
427  k3*m_kk*m_kk + k2*m_kk + k1,
428  k3*m_kk*m_kk + k1*m_kk + k2};
429  for (auto c : cc) {
430  for (size_t n = 0; n < nParams; n++) {
431  m_Psi_ijk_coeff(n, c) = psi[n];
432  }
433  m_Psi_ijk[c] = psi[0];
434  }
435 }
436 
437 void HMWSoln::setLambda(const string& sp1, const string& sp2,
438  size_t nParams, double* lambda)
439 {
440  size_t k1 = speciesIndex(sp1);
441  size_t k2 = speciesIndex(sp2);
442  if (k1 == npos) {
443  throw CanteraError("HMWSoln::setLambda", "Species '{}' not found", sp1);
444  } else if (k2 == npos) {
445  throw CanteraError("HMWSoln::setLambda", "Species '{}' not found", sp2);
446  }
447 
448  if (charge(k1) != 0 && charge(k2) != 0) {
449  throw CanteraError("HMWSoln::setLambda", "Expected at least one neutral"
450  " species, but given species (charges) were: {} ({}) and {} ({}).",
451  sp1, charge(k1), sp2, charge(k2));
452  }
453  if (charge(k1) != 0) {
454  std::swap(k1, k2);
455  }
456  check_nParams("HMWSoln::setLambda", nParams, m_formPitzerTemp);
457  size_t c = k1*m_kk + k2;
458  for (size_t n = 0; n < nParams; n++) {
459  m_Lambda_nj_coeff(n, c) = lambda[n];
460  }
461  m_Lambda_nj(k1, k2) = lambda[0];
462 }
463 
464 void HMWSoln::setMunnn(const string& sp, size_t nParams, double* munnn)
465 {
466  size_t k = speciesIndex(sp);
467  if (k == npos) {
468  throw CanteraError("HMWSoln::setMunnn", "Species '{}' not found", sp);
469  }
470 
471  if (charge(k) != 0) {
472  throw CanteraError("HMWSoln::setMunnn", "Expected a neutral species,"
473  " got {} ({}).", sp, charge(k));
474  }
475  check_nParams("HMWSoln::setMunnn", nParams, m_formPitzerTemp);
476  for (size_t n = 0; n < nParams; n++) {
477  m_Mu_nnn_coeff(n, k) = munnn[n];
478  }
479  m_Mu_nnn[k] = munnn[0];
480 }
481 
482 void HMWSoln::setZeta(const string& sp1, const string& sp2,
483  const string& sp3, size_t nParams, double* psi)
484 {
485  size_t k1 = speciesIndex(sp1);
486  size_t k2 = speciesIndex(sp2);
487  size_t k3 = speciesIndex(sp3);
488  if (k1 == npos) {
489  throw CanteraError("HMWSoln::setZeta", "Species '{}' not found", sp1);
490  } else if (k2 == npos) {
491  throw CanteraError("HMWSoln::setZeta", "Species '{}' not found", sp2);
492  } else if (k3 == npos) {
493  throw CanteraError("HMWSoln::setZeta", "Species '{}' not found", sp3);
494  }
495 
496  if (charge(k1)*charge(k2)*charge(k3) != 0 ||
497  sign(charge(k1)) + sign(charge(k2)) + sign(charge(k3)) != 0) {
498  throw CanteraError("HMWSoln::setZeta", "Requires one neutral species, "
499  "one cation, and one anion, but given species (charges) were: "
500  "{} ({}), {} ({}), and {} ({}).",
501  sp1, charge(k1), sp2, charge(k2), sp3, charge(k3));
502  }
503 
504  //! Make k1 the neutral species
505  if (charge(k2) == 0) {
506  std::swap(k1, k2);
507  } else if (charge(k3) == 0) {
508  std::swap(k1, k3);
509  }
510 
511  // Make k2 the cation
512  if (charge(k3) > 0) {
513  std::swap(k2, k3);
514  }
515 
516  check_nParams("HMWSoln::setZeta", nParams, m_formPitzerTemp);
517  // In contrast to setPsi, there are no duplicate entries
518  size_t c = k1 * m_kk *m_kk + k2 * m_kk + k3;
519  for (size_t n = 0; n < nParams; n++) {
520  m_Psi_ijk_coeff(n, c) = psi[n];
521  }
522  m_Psi_ijk[c] = psi[0];
523 }
524 
525 void HMWSoln::setPitzerTempModel(const string& model)
526 {
527  if (caseInsensitiveEquals(model, "constant") || caseInsensitiveEquals(model, "default")) {
528  m_formPitzerTemp = PITZER_TEMP_CONSTANT;
529  } else if (caseInsensitiveEquals(model, "linear")) {
530  m_formPitzerTemp = PITZER_TEMP_LINEAR;
531  } else if (caseInsensitiveEquals(model, "complex") || caseInsensitiveEquals(model, "complex1")) {
532  m_formPitzerTemp = PITZER_TEMP_COMPLEX1;
533  } else {
534  throw CanteraError("HMWSoln::setPitzerTempModel",
535  "Unknown Pitzer ActivityCoeff Temp model: {}", model);
536  }
537 }
538 
539 void HMWSoln::setA_Debye(double A)
540 {
541  if (A < 0) {
542  m_form_A_Debye = A_DEBYE_WATER;
543  } else {
544  m_form_A_Debye = A_DEBYE_CONST;
545  m_A_Debye = A;
546  }
547 }
548 
549 void HMWSoln::setCroppingCoefficients(double ln_gamma_k_min,
550  double ln_gamma_k_max, double ln_gamma_o_min, double ln_gamma_o_max)
551 {
552  CROP_ln_gamma_k_min = ln_gamma_k_min;
553  CROP_ln_gamma_k_max = ln_gamma_k_max;
554  CROP_ln_gamma_o_min = ln_gamma_o_min;
555  CROP_ln_gamma_o_max = ln_gamma_o_max;
556 }
557 
558 vector<double> getSizedVector(const AnyMap& item, const string& key, size_t nCoeffs)
559 {
560  vector<double> v;
561  if (item[key].is<double>()) {
562  // Allow a single value to be given directly, rather than as a list of
563  // one item
564  v.push_back(item[key].asDouble());
565  } else {
566  v = item[key].asVector<double>(1, nCoeffs);
567  }
568  if (v.size() == 1 && nCoeffs == 5) {
569  // Adapt constant-temperature data to be compatible with the "complex"
570  // temperature model
571  v.resize(5, 0.0);
572  }
573  return v;
574 }
575 
577 {
579  if (m_input.hasKey("activity-data")) {
580  auto& actData = m_input["activity-data"].as<AnyMap>();
581  setPitzerTempModel(actData["temperature-model"].asString());
582  initLengths();
583  size_t nCoeffs = 1;
584  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
585  nCoeffs = 2;
586  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
587  nCoeffs = 5;
588  }
589  if (actData.hasKey("A_Debye")) {
590  if (actData["A_Debye"] == "variable") {
591  setA_Debye(-1);
592  } else {
593  setA_Debye(actData.convert("A_Debye", "kg^0.5/gmol^0.5"));
594  }
595  }
596  if (actData.hasKey("max-ionic-strength")) {
597  setMaxIonicStrength(actData["max-ionic-strength"].asDouble());
598  }
599  if (actData.hasKey("interactions")) {
600  for (auto& item : actData["interactions"].asVector<AnyMap>()) {
601  auto& species = item["species"].asVector<string>(1, 3);
602  size_t nsp = species.size();
603  double q0 = charge(speciesIndex(species[0]));
604  double q1 = (nsp > 1) ? charge(speciesIndex(species[1])) : 0;
605  double q2 = (nsp == 3) ? charge(speciesIndex(species[2])) : 0;
606  if (nsp == 2 && q0 * q1 < 0) {
607  // Two species with opposite charges - binary salt
608  vector<double> beta0 = getSizedVector(item, "beta0", nCoeffs);
609  vector<double> beta1 = getSizedVector(item, "beta1", nCoeffs);
610  vector<double> beta2 = getSizedVector(item, "beta2", nCoeffs);
611  vector<double> Cphi = getSizedVector(item, "Cphi", nCoeffs);
612  if (beta0.size() != beta1.size() || beta0.size() != beta2.size()
613  || beta0.size() != Cphi.size()) {
614  throw InputFileError("HMWSoln::initThermo", item,
615  "Inconsistent binary salt array sizes ({}, {}, {}, {})",
616  beta0.size(), beta1.size(), beta2.size(), Cphi.size());
617  }
618  double alpha1 = item["alpha1"].asDouble();
619  double alpha2 = item.getDouble("alpha2", 0.0);
620  setBinarySalt(species[0], species[1], beta0.size(),
621  beta0.data(), beta1.data(), beta2.data(), Cphi.data(),
622  alpha1, alpha2);
623  } else if (nsp == 2 && q0 * q1 > 0) {
624  // Two species with like charges - "theta" interaction
625  vector<double> theta = getSizedVector(item, "theta", nCoeffs);
626  setTheta(species[0], species[1], theta.size(), theta.data());
627  } else if (nsp == 2 && q0 * q1 == 0) {
628  // Two species, including at least one neutral
629  vector<double> lambda = getSizedVector(item, "lambda", nCoeffs);
630  setLambda(species[0], species[1], lambda.size(), lambda.data());
631  } else if (nsp == 3 && q0 * q1 * q2 != 0) {
632  // Three charged species - "psi" interaction
633  vector<double> psi = getSizedVector(item, "psi", nCoeffs);
634  setPsi(species[0], species[1], species[2],
635  psi.size(), psi.data());
636  } else if (nsp == 3 && q0 * q1 * q2 == 0) {
637  // Three species, including one neutral
638  vector<double> zeta = getSizedVector(item, "zeta", nCoeffs);
639  setZeta(species[0], species[1], species[2],
640  zeta.size(), zeta.data());
641  } else if (nsp == 1) {
642  // single species (should be neutral)
643  vector<double> mu = getSizedVector(item, "mu", nCoeffs);
644  setMunnn(species[0], mu.size(), mu.data());
645  }
646  }
647  }
648  if (actData.hasKey("cropping-coefficients")) {
649  auto& crop = actData["cropping-coefficients"].as<AnyMap>();
650  setCroppingCoefficients(
651  crop.getDouble("ln_gamma_k_min", crop_ln_gamma_k_min_default),
652  crop.getDouble("ln_gamma_k_max", crop_ln_gamma_k_max_default),
653  crop.getDouble("ln_gamma_o_min", crop_ln_gamma_o_min_default),
654  crop.getDouble("ln_gamma_o_max", crop_ln_gamma_o_max_default));
655  }
656  } else {
657  initLengths();
658  }
659 
660  for (int i = 0; i < 17; i++) {
661  elambda[i] = 0.0;
662  elambda1[i] = 0.0;
663  }
664 
665  // Store a local pointer to the water standard state model.
666  m_waterSS = providePDSS(0);
667 
668  // Initialize the water property calculator. It will share the internal eos
669  // water calculator.
670  m_waterProps = make_unique<WaterProps>(dynamic_cast<PDSS_Water*>(m_waterSS));
671 
672  // Lastly calculate the charge balance and then add stuff until the charges
673  // compensate
674  vector<double> mf(m_kk, 0.0);
675  getMoleFractions(mf.data());
676  bool notDone = true;
677 
678  while (notDone) {
679  double sum = 0.0;
680  size_t kMaxC = npos;
681  double MaxC = 0.0;
682  for (size_t k = 0; k < m_kk; k++) {
683  sum += mf[k] * charge(k);
684  if (fabs(mf[k] * charge(k)) > MaxC) {
685  kMaxC = k;
686  }
687  }
688  size_t kHp = speciesIndex("H+");
689  size_t kOHm = speciesIndex("OH-");
690 
691  if (fabs(sum) > 1.0E-30) {
692  if (kHp != npos) {
693  if (mf[kHp] > sum * 1.1) {
694  mf[kHp] -= sum;
695  mf[0] += sum;
696  notDone = false;
697  } else {
698  if (sum > 0.0) {
699  mf[kHp] *= 0.5;
700  mf[0] += mf[kHp];
701  sum -= mf[kHp];
702  }
703  }
704  }
705  if (notDone) {
706  if (kOHm != npos) {
707  if (mf[kOHm] > -sum * 1.1) {
708  mf[kOHm] += sum;
709  mf[0] -= sum;
710  notDone = false;
711  } else {
712  if (sum < 0.0) {
713  mf[kOHm] *= 0.5;
714  mf[0] += mf[kOHm];
715  sum += mf[kOHm];
716  }
717  }
718  }
719  if (notDone && kMaxC != npos) {
720  if (mf[kMaxC] > (1.1 * sum / charge(kMaxC))) {
721  mf[kMaxC] -= sum / charge(kMaxC);
722  mf[0] += sum / charge(kMaxC);
723  } else {
724  mf[kMaxC] *= 0.5;
725  mf[0] += mf[kMaxC];
726  notDone = true;
727  }
728  }
729  }
730  setMoleFractions(mf.data());
731  } else {
732  notDone = false;
733  }
734  }
735 
738  setMoleFSolventMin(1.0E-5);
739 }
740 
741 void assignTrimmed(AnyMap& interaction, const string& key, vector<double>& values) {
742  while (values.size() > 1 && values.back() == 0) {
743  values.pop_back();
744  }
745  if (values.size() == 1) {
746  interaction[key] = values[0];
747  } else {
748  interaction[key] = values;
749  }
750 }
751 
752 void HMWSoln::getParameters(AnyMap& phaseNode) const
753 {
755  AnyMap activityNode;
756  size_t nParams = 1;
757  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
758  activityNode["temperature-model"] = "linear";
759  nParams = 2;
760  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
761  activityNode["temperature-model"] = "complex";
762  nParams = 5;
763  }
764 
765  if (m_form_A_Debye == A_DEBYE_WATER) {
766  activityNode["A_Debye"] = "variable";
767  } else if (m_A_Debye != A_Debye_default) {
768  activityNode["A_Debye"].setQuantity(m_A_Debye, "kg^0.5/gmol^0.5");
769  }
770  if (m_maxIionicStrength != maxIionicStrength_default) {
771  activityNode["max-ionic-strength"] = m_maxIionicStrength;
772  }
773 
774  vector<AnyMap> interactions;
775 
776  // Binary interactions
777  for (size_t i = 1; i < m_kk; i++) {
778  for (size_t j = 1; j < m_kk; j++) {
779  size_t c = i*m_kk + j;
780  // lambda: neutral-charged / neutral-neutral interactions
781  bool lambda_found = false;
782  for (size_t n = 0; n < nParams; n++) {
783  if (m_Lambda_nj_coeff(n, c)) {
784  lambda_found = true;
785  break;
786  }
787  }
788  if (lambda_found) {
789  AnyMap interaction;
790  interaction["species"] = vector<string>{
791  speciesName(i), speciesName(j)};
792  vector<double> lambda(nParams);
793  for (size_t n = 0; n < nParams; n++) {
794  lambda[n] = m_Lambda_nj_coeff(n, c);
795  }
796  assignTrimmed(interaction, "lambda", lambda);
797  interactions.push_back(std::move(interaction));
798  continue;
799  }
800 
801  c = static_cast<size_t>(m_CounterIJ[m_kk * i + j]);
802  if (c == 0 || i > j) {
803  continue;
804  }
805 
806  // beta: opposite charged interactions
807  bool salt_found = false;
808  for (size_t n = 0; n < nParams; n++) {
809  if (m_Beta0MX_ij_coeff(n, c) || m_Beta1MX_ij_coeff(n, c) ||
810  m_Beta2MX_ij_coeff(n, c) || m_CphiMX_ij_coeff(n, c))
811  {
812  salt_found = true;
813  break;
814  }
815  }
816  if (salt_found) {
817  AnyMap interaction;
818  interaction["species"] = vector<string>{
819  speciesName(i), speciesName(j)};
820  vector<double> beta0(nParams), beta1(nParams), beta2(nParams), Cphi(nParams);
821  size_t last_nonzero = 0;
822  for (size_t n = 0; n < nParams; n++) {
823  beta0[n] = m_Beta0MX_ij_coeff(n, c);
824  beta1[n] = m_Beta1MX_ij_coeff(n, c);
825  beta2[n] = m_Beta2MX_ij_coeff(n, c);
826  Cphi[n] = m_CphiMX_ij_coeff(n, c);
827  if (beta0[n] || beta1[n] || beta2[n] || Cphi[n]) {
828  last_nonzero = n;
829  }
830  }
831  if (last_nonzero == 0) {
832  interaction["beta0"] = beta0[0];
833  interaction["beta1"] = beta1[0];
834  interaction["beta2"] = beta2[0];
835  interaction["Cphi"] = Cphi[0];
836  } else {
837  beta0.resize(1 + last_nonzero);
838  beta1.resize(1 + last_nonzero);
839  beta2.resize(1 + last_nonzero);
840  Cphi.resize(1 + last_nonzero);
841  interaction["beta0"] = beta0;
842  interaction["beta1"] = beta1;
843  interaction["beta2"] = beta2;
844  interaction["Cphi"] = Cphi;
845  }
846  interaction["alpha1"] = m_Alpha1MX_ij[c];
847  if (m_Alpha2MX_ij[c]) {
848  interaction["alpha2"] = m_Alpha2MX_ij[c];
849  }
850  interactions.push_back(std::move(interaction));
851  continue;
852  }
853 
854  // theta: like-charge interactions
855  bool theta_found = false;
856  for (size_t n = 0; n < nParams; n++) {
857  if (m_Theta_ij_coeff(n, c)) {
858  theta_found = true;
859  break;
860  }
861  }
862  if (theta_found) {
863  AnyMap interaction;
864  interaction["species"] = vector<string>{
865  speciesName(i), speciesName(j)};
866  vector<double> theta(nParams);
867  for (size_t n = 0; n < nParams; n++) {
868  theta[n] = m_Theta_ij_coeff(n, c);
869  }
870  assignTrimmed(interaction, "theta", theta);
871  interactions.push_back(std::move(interaction));
872  continue;
873  }
874  }
875  }
876 
877  // psi: ternary charged species interactions
878  // Need to check species charges because both psi and zeta parameters
879  // are stored in m_Psi_ijk_coeff
880  for (size_t i = 1; i < m_kk; i++) {
881  if (charge(i) == 0) {
882  continue;
883  }
884  for (size_t j = i + 1; j < m_kk; j++) {
885  if (charge(j) == 0) {
886  continue;
887  }
888  for (size_t k = j + 1; k < m_kk; k++) {
889  if (charge(k) == 0) {
890  continue;
891  }
892  size_t c = i*m_kk*m_kk + j*m_kk + k;
893  for (size_t n = 0; n < nParams; n++) {
894  if (m_Psi_ijk_coeff(n, c) != 0) {
895  AnyMap interaction;
896  interaction["species"] = vector<string>{
897  speciesName(i), speciesName(j), speciesName(k)};
898  vector<double> psi(nParams);
899  for (size_t m = 0; m < nParams; m++) {
900  psi[m] = m_Psi_ijk_coeff(m, c);
901  }
902  assignTrimmed(interaction, "psi", psi);
903  interactions.push_back(std::move(interaction));
904  break;
905  }
906  }
907  }
908  }
909  }
910 
911  // zeta: neutral-cation-anion interactions
912  for (size_t i = 1; i < m_kk; i++) {
913  if (charge(i) != 0) {
914  continue; // first species must be neutral
915  }
916  for (size_t j = 1; j < m_kk; j++) {
917  if (charge(j) <= 0) {
918  continue; // second species must be cation
919  }
920  for (size_t k = 1; k < m_kk; k++) {
921  if (charge(k) >= 0) {
922  continue; // third species must be anion
923  }
924  size_t c = i*m_kk*m_kk + j*m_kk + k;
925  for (size_t n = 0; n < nParams; n++) {
926  if (m_Psi_ijk_coeff(n, c) != 0) {
927  AnyMap interaction;
928  interaction["species"] = vector<string>{
929  speciesName(i), speciesName(j), speciesName(k)};
930  vector<double> zeta(nParams);
931  for (size_t m = 0; m < nParams; m++) {
932  zeta[m] = m_Psi_ijk_coeff(m, c);
933  }
934  assignTrimmed(interaction, "zeta", zeta);
935  interactions.push_back(std::move(interaction));
936  break;
937  }
938  }
939  }
940  }
941  }
942 
943  // mu: neutral self-interaction
944  for (size_t i = 1; i < m_kk; i++) {
945  for (size_t n = 0; n < nParams; n++) {
946  if (m_Mu_nnn_coeff(n, i) != 0) {
947  AnyMap interaction;
948  interaction["species"] = vector<string>{speciesName(i)};
949  vector<double> mu(nParams);
950  for (size_t m = 0; m < nParams; m++) {
951  mu[m] = m_Mu_nnn_coeff(m, i);
952  }
953  assignTrimmed(interaction, "mu", mu);
954  interactions.push_back(std::move(interaction));
955  break;
956  }
957  }
958  }
959 
960  activityNode["interactions"] = std::move(interactions);
961 
962  AnyMap croppingCoeffs;
963  if (CROP_ln_gamma_k_min != crop_ln_gamma_k_min_default) {
964  croppingCoeffs["ln_gamma_k_min"] = CROP_ln_gamma_k_min;
965  }
966  if (CROP_ln_gamma_k_max != crop_ln_gamma_k_max_default) {
967  croppingCoeffs["ln_gamma_k_max"] = CROP_ln_gamma_k_max;
968  }
969  if (CROP_ln_gamma_o_min != crop_ln_gamma_o_min_default) {
970  croppingCoeffs["ln_gamma_o_min"] = CROP_ln_gamma_o_min;
971  }
972  if (CROP_ln_gamma_o_max != crop_ln_gamma_o_max_default) {
973  croppingCoeffs["ln_gamma_o_max"] = CROP_ln_gamma_o_max;
974  }
975  if (croppingCoeffs.size()) {
976  activityNode["cropping-coefficients"] = std::move(croppingCoeffs);
977  }
978 
979  phaseNode["activity-data"] = std::move(activityNode);
980 }
981 
982 double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
983 {
984  double T = temperature();
985  double A;
986  if (tempArg != -1.0) {
987  T = tempArg;
988  }
989  double P = pressure();
990  if (presArg != -1.0) {
991  P = presArg;
992  }
993 
994  static const int cacheId = m_cache.getId();
995  CachedScalar cached = m_cache.getScalar(cacheId);
996  if(cached.validate(T, P)) {
997  return m_A_Debye;
998  }
999 
1000  switch (m_form_A_Debye) {
1001  case A_DEBYE_CONST:
1002  A = m_A_Debye;
1003  break;
1004  case A_DEBYE_WATER:
1005  A = m_waterProps->ADebye(T, P, 0);
1006  m_A_Debye = A;
1007  break;
1008  default:
1009  throw CanteraError("HMWSoln::A_Debye_TP", "shouldn't be here");
1010  }
1011  return A;
1012 }
1013 
1014 double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
1015 {
1016  double T = temperature();
1017  if (tempArg != -1.0) {
1018  T = tempArg;
1019  }
1020  double P = pressure();
1021  if (presArg != -1.0) {
1022  P = presArg;
1023  }
1024  double dAdT;
1025  switch (m_form_A_Debye) {
1026  case A_DEBYE_CONST:
1027  dAdT = 0.0;
1028  break;
1029  case A_DEBYE_WATER:
1030  dAdT = m_waterProps->ADebye(T, P, 1);
1031  break;
1032  default:
1033  throw CanteraError("HMWSoln::dA_DebyedT_TP", "shouldn't be here");
1034  }
1035  return dAdT;
1036 }
1037 
1038 double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
1039 {
1040  double T = temperature();
1041  if (tempArg != -1.0) {
1042  T = tempArg;
1043  }
1044  double P = pressure();
1045  if (presArg != -1.0) {
1046  P = presArg;
1047  }
1048 
1049  double dAdP;
1050  static const int cacheId = m_cache.getId();
1051  CachedScalar cached = m_cache.getScalar(cacheId);
1052  switch (m_form_A_Debye) {
1053  case A_DEBYE_CONST:
1054  dAdP = 0.0;
1055  break;
1056  case A_DEBYE_WATER:
1057  if(cached.validate(T, P)) {
1058  dAdP = cached.value;
1059  } else {
1060  dAdP = m_waterProps->ADebye(T, P, 3);
1061  cached.value = dAdP;
1062  }
1063  break;
1064  default:
1065  throw CanteraError("HMWSoln::dA_DebyedP_TP", "shouldn't be here");
1066  }
1067  return dAdP;
1068 }
1069 
1070 double HMWSoln::ADebye_L(double tempArg, double presArg) const
1071 {
1072  double dAdT = dA_DebyedT_TP();
1073  double dAphidT = dAdT /3.0;
1074  double T = temperature();
1075  if (tempArg != -1.0) {
1076  T = tempArg;
1077  }
1078  return dAphidT * (4.0 * GasConstant * T * T);
1079 }
1080 
1081 double HMWSoln::ADebye_V(double tempArg, double presArg) const
1082 {
1083  double dAdP = dA_DebyedP_TP();
1084  double dAphidP = dAdP /3.0;
1085  double T = temperature();
1086  if (tempArg != -1.0) {
1087  T = tempArg;
1088  }
1089  return - dAphidP * (4.0 * GasConstant * T);
1090 }
1091 
1092 double HMWSoln::ADebye_J(double tempArg, double presArg) const
1093 {
1094  double T = temperature();
1095  if (tempArg != -1.0) {
1096  T = tempArg;
1097  }
1098  double A_L = ADebye_L(T, presArg);
1099  double d2 = d2A_DebyedT2_TP(T, presArg);
1100  double d2Aphi = d2 / 3.0;
1101  return 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi;
1102 }
1103 
1104 double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
1105 {
1106  double T = temperature();
1107  if (tempArg != -1.0) {
1108  T = tempArg;
1109  }
1110  double P = pressure();
1111  if (presArg != -1.0) {
1112  P = presArg;
1113  }
1114  double d2AdT2;
1115  switch (m_form_A_Debye) {
1116  case A_DEBYE_CONST:
1117  d2AdT2 = 0.0;
1118  break;
1119  case A_DEBYE_WATER:
1120  d2AdT2 = m_waterProps->ADebye(T, P, 2);
1121  break;
1122  default:
1123  throw CanteraError("HMWSoln::d2A_DebyedT2_TP", "shouldn't be here");
1124  }
1125  return d2AdT2;
1126 }
1127 
1128 // ---------- Other Property Functions
1129 
1130 // ------------ Private and Restricted Functions ------------------
1131 
1133 {
1134  m_tmpV.resize(m_kk, 0.0);
1135  m_molalitiesCropped.resize(m_kk, 0.0);
1136 
1137  size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
1138 
1139  // Figure out the size of the temperature coefficient arrays
1140  int TCoeffLength = 1;
1141  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1142  TCoeffLength = 2;
1143  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1144  TCoeffLength = 5;
1145  }
1146 
1147  m_Beta0MX_ij.resize(maxCounterIJlen, 0.0);
1148  m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0);
1149  m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0);
1150  m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0);
1151  m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1152 
1153  m_Beta1MX_ij.resize(maxCounterIJlen, 0.0);
1154  m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0);
1155  m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0);
1156  m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0);
1157  m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1158 
1159  m_Beta2MX_ij.resize(maxCounterIJlen, 0.0);
1160  m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0);
1161  m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0);
1162  m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0);
1163  m_Beta2MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1164 
1165  m_CphiMX_ij.resize(maxCounterIJlen, 0.0);
1166  m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0);
1167  m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0);
1168  m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0);
1169  m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1170 
1171  m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1172  m_Alpha2MX_ij.resize(maxCounterIJlen, 12.0);
1173  m_Theta_ij.resize(maxCounterIJlen, 0.0);
1174  m_Theta_ij_L.resize(maxCounterIJlen, 0.0);
1175  m_Theta_ij_LL.resize(maxCounterIJlen, 0.0);
1176  m_Theta_ij_P.resize(maxCounterIJlen, 0.0);
1177  m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1178 
1179  size_t n = m_kk*m_kk*m_kk;
1180  m_Psi_ijk.resize(n, 0.0);
1181  m_Psi_ijk_L.resize(n, 0.0);
1182  m_Psi_ijk_LL.resize(n, 0.0);
1183  m_Psi_ijk_P.resize(n, 0.0);
1184  m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0);
1185 
1186  m_Lambda_nj.resize(m_kk, m_kk, 0.0);
1187  m_Lambda_nj_L.resize(m_kk, m_kk, 0.0);
1188  m_Lambda_nj_LL.resize(m_kk, m_kk, 0.0);
1189  m_Lambda_nj_P.resize(m_kk, m_kk, 0.0);
1190  m_Lambda_nj_coeff.resize(TCoeffLength, m_kk * m_kk, 0.0);
1191 
1192  m_Mu_nnn.resize(m_kk, 0.0);
1193  m_Mu_nnn_L.resize(m_kk, 0.0);
1194  m_Mu_nnn_LL.resize(m_kk, 0.0);
1195  m_Mu_nnn_P.resize(m_kk, 0.0);
1196  m_Mu_nnn_coeff.resize(TCoeffLength, m_kk, 0.0);
1197 
1198  m_lnActCoeffMolal_Scaled.resize(m_kk, 0.0);
1199  m_dlnActCoeffMolaldT_Scaled.resize(m_kk, 0.0);
1200  m_d2lnActCoeffMolaldT2_Scaled.resize(m_kk, 0.0);
1201  m_dlnActCoeffMolaldP_Scaled.resize(m_kk, 0.0);
1202  m_lnActCoeffMolal_Unscaled.resize(m_kk, 0.0);
1203  m_dlnActCoeffMolaldT_Unscaled.resize(m_kk, 0.0);
1205  m_dlnActCoeffMolaldP_Unscaled.resize(m_kk, 0.0);
1206 
1207  m_CounterIJ.resize(m_kk*m_kk, 0);
1208  m_gfunc_IJ.resize(maxCounterIJlen, 0.0);
1209  m_g2func_IJ.resize(maxCounterIJlen, 0.0);
1210  m_hfunc_IJ.resize(maxCounterIJlen, 0.0);
1211  m_h2func_IJ.resize(maxCounterIJlen, 0.0);
1212  m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1213  m_BMX_IJ_L.resize(maxCounterIJlen, 0.0);
1214  m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1215  m_BMX_IJ_P.resize(maxCounterIJlen, 0.0);
1216  m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0);
1217  m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0);
1218  m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1219  m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0);
1220  m_BphiMX_IJ.resize(maxCounterIJlen, 0.0);
1221  m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0);
1222  m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1223  m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0);
1224  m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1225  m_Phi_IJ_L.resize(maxCounterIJlen, 0.0);
1226  m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0);
1227  m_Phi_IJ_P.resize(maxCounterIJlen, 0.0);
1228  m_Phiprime_IJ.resize(maxCounterIJlen, 0.0);
1229  m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0);
1230  m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0);
1231  m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0);
1232  m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0);
1233  m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1234  m_CMX_IJ_L.resize(maxCounterIJlen, 0.0);
1235  m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1236  m_CMX_IJ_P.resize(maxCounterIJlen, 0.0);
1237 
1238  m_gamma_tmp.resize(m_kk, 0.0);
1239  IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
1240  CROP_speciesCropped_.resize(m_kk, 0);
1241 
1242  counterIJ_setup();
1243 }
1244 
1246 {
1247  static const int cacheId = m_cache.getId();
1248  CachedScalar cached = m_cache.getScalar(cacheId);
1249  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
1250  return;
1251  }
1252 
1253  // Calculate the molalities. Currently, the molalities may not be current
1254  // with respect to the contents of the State objects' data.
1255  calcMolalities();
1256 
1257  // Calculate a cropped set of molalities that will be used in all activity
1258  // coefficient calculations.
1260 
1261  // Update the temperature dependence of the Pitzer coefficients and their
1262  // derivatives
1264 
1265  // Calculate the IMS cutoff factors
1267 
1268  // Now do the main calculation.
1270 
1271  double xmolSolvent = moleFraction(0);
1272  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
1273  double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1274  double lnxs = log(xx);
1275 
1276  for (size_t k = 1; k < m_kk; k++) {
1277  CROP_speciesCropped_[k] = 0;
1279  if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
1280  CROP_speciesCropped_[k] = 2;
1281  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
1282  }
1283  if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
1284  // -1.0 and -1.5 caused multiple solutions
1285  CROP_speciesCropped_[k] = 2;
1286  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
1287  }
1288  }
1289  CROP_speciesCropped_[0] = 0;
1290  m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
1291  if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
1292  CROP_speciesCropped_[0] = 2;
1293  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
1294  }
1295  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
1296  CROP_speciesCropped_[0] = 2;
1297  // -0.5 caused multiple solutions
1298  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
1299  }
1300  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
1301  CROP_speciesCropped_[0] = 2;
1302  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
1303  }
1304 
1305  // Now do the pH Scaling
1307 }
1308 
1310 {
1311  double Imax = 0.0;
1312  m_molalitiesAreCropped = false;
1313 
1314  for (size_t k = 0; k < m_kk; k++) {
1316  Imax = std::max(m_molalities[k] * charge(k) * charge(k), Imax);
1317  }
1318 
1319  int cropMethod = 1;
1320  if (cropMethod == 0) {
1321  // Quick return
1322  if (Imax < m_maxIionicStrength) {
1323  return;
1324  }
1325 
1326  m_molalitiesAreCropped = true;
1327  for (size_t i = 1; i < (m_kk - 1); i++) {
1328  double charge_i = charge(i);
1329  double abs_charge_i = fabs(charge_i);
1330  if (charge_i == 0.0) {
1331  continue;
1332  }
1333  for (size_t j = (i+1); j < m_kk; j++) {
1334  double charge_j = charge(j);
1335  double abs_charge_j = fabs(charge_j);
1336 
1337  // Only loop over oppositely charge species
1338  if (charge_i * charge_j < 0) {
1339  double Iac_max = m_maxIionicStrength;
1340 
1342  Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
1343  if (Imax > Iac_max) {
1344  m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
1345  }
1346  Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
1347  if (Imax > Iac_max) {
1348  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
1349  }
1350  } else {
1351  Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
1352  if (Imax > Iac_max) {
1353  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
1354  }
1355  Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
1356  if (Imax > Iac_max) {
1357  m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
1358  }
1359  }
1360  }
1361  }
1362  }
1363 
1364  // Do this loop 10 times until we have achieved charge neutrality in the
1365  // cropped molalities
1366  for (int times = 0; times< 10; times++) {
1367  double anion_charge = 0.0;
1368  double cation_charge = 0.0;
1369  size_t anion_contrib_max_i = npos;
1370  double anion_contrib_max = -1.0;
1371  size_t cation_contrib_max_i = npos;
1372  double cation_contrib_max = -1.0;
1373  for (size_t i = 0; i < m_kk; i++) {
1374  double charge_i = charge(i);
1375  if (charge_i < 0.0) {
1376  double anion_contrib = - m_molalitiesCropped[i] * charge_i;
1377  anion_charge += anion_contrib;
1378  if (anion_contrib > anion_contrib_max) {
1379  anion_contrib_max = anion_contrib;
1380  anion_contrib_max_i = i;
1381  }
1382  } else if (charge_i > 0.0) {
1383  double cation_contrib = m_molalitiesCropped[i] * charge_i;
1384  cation_charge += cation_contrib;
1385  if (cation_contrib > cation_contrib_max) {
1386  cation_contrib_max = cation_contrib;
1387  cation_contrib_max_i = i;
1388  }
1389  }
1390  }
1391  double total_charge = cation_charge - anion_charge;
1392  if (total_charge > 1.0E-8) {
1393  double desiredCrop = total_charge/charge(cation_contrib_max_i);
1394  double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
1395  if (desiredCrop < maxCrop) {
1396  m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
1397  break;
1398  } else {
1399  m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
1400  }
1401  } else if (total_charge < -1.0E-8) {
1402  double desiredCrop = total_charge/charge(anion_contrib_max_i);
1403  double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
1404  if (desiredCrop < maxCrop) {
1405  m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
1406  break;
1407  } else {
1408  m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
1409  }
1410  } else {
1411  break;
1412  }
1413  }
1414  }
1415 
1416  if (cropMethod == 1) {
1417  double* molF = m_gamma_tmp.data();
1418  getMoleFractions(molF);
1419  double xmolSolvent = molF[0];
1420  if (xmolSolvent >= MC_X_o_cutoff_) {
1421  return;
1422  }
1423 
1424  m_molalitiesAreCropped = true;
1425  double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1426  double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1427  double denomInv = 1.0/ (m_Mnaught * p);
1428  for (size_t k = 0; k < m_kk; k++) {
1429  m_molalitiesCropped[k] = molF[k] * denomInv;
1430  }
1431 
1432  // Do a further check to see if the Ionic strength is below a max value
1433  // Reduce the molalities to enforce this. Note, this algorithm preserves
1434  // the charge neutrality of the solution after cropping.
1435  double Itmp = 0.0;
1436  for (size_t k = 0; k < m_kk; k++) {
1437  Itmp += m_molalitiesCropped[k] * charge(k) * charge(k);
1438  }
1439  if (Itmp > m_maxIionicStrength) {
1440  double ratio = Itmp / m_maxIionicStrength;
1441  for (size_t k = 0; k < m_kk; k++) {
1442  if (charge(k) != 0.0) {
1443  m_molalitiesCropped[k] *= ratio;
1444  }
1445  }
1446  }
1447  }
1448 }
1449 
1451 {
1452  m_CounterIJ.resize(m_kk * m_kk);
1453  int counter = 0;
1454  for (size_t i = 0; i < m_kk; i++) {
1455  size_t n = i;
1456  size_t nc = m_kk * i;
1457  m_CounterIJ[n] = 0;
1458  m_CounterIJ[nc] = 0;
1459  }
1460  for (size_t i = 1; i < (m_kk - 1); i++) {
1461  size_t n = m_kk * i + i;
1462  m_CounterIJ[n] = 0;
1463  for (size_t j = (i+1); j < m_kk; j++) {
1464  n = m_kk * j + i;
1465  size_t nc = m_kk * i + j;
1466  counter++;
1467  m_CounterIJ[n] = counter;
1468  m_CounterIJ[nc] = counter;
1469  }
1470  }
1471 }
1472 
1474 {
1475  double IMS_gamma_o_min_ = 1.0E-5; // value at the zero solvent point
1476  double IMS_gamma_k_min_ = 10.0; // minimum at the zero solvent point
1477  double IMS_slopefCut_ = 0.6; // slope of the f function at the zero solvent point
1478 
1479  IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
1480  IMS_efCut_ = 0.0;
1481  bool converged = false;
1482  double oldV = 0.0;
1483  for (int its = 0; its < 100 && !converged; its++) {
1484  oldV = IMS_efCut_;
1485  IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
1486  IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
1487  IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
1488  /
1490  double tmp = IMS_afCut_ + IMS_X_o_cutoff_*(IMS_bfCut_ + IMS_dfCut_ *IMS_X_o_cutoff_);
1491  double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1492  IMS_efCut_ = - eterm * tmp;
1493  if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
1494  converged = true;
1495  }
1496  }
1497  if (!converged) {
1498  throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1499  " failed to converge on the f polynomial");
1500  }
1501  converged = false;
1502  double f_0 = IMS_afCut_ + IMS_efCut_;
1503  double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
1504  IMS_egCut_ = 0.0;
1505  for (int its = 0; its < 100 && !converged; its++) {
1506  oldV = IMS_egCut_;
1507  double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
1508  IMS_agCut_ = exp(lng_0) - IMS_egCut_;
1509  IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
1510  IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
1511  /
1513  double tmp = IMS_agCut_ + IMS_X_o_cutoff_*(IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
1514  double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1515  IMS_egCut_ = - eterm * tmp;
1516  if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
1517  converged = true;
1518  }
1519  }
1520  if (!converged) {
1521  throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1522  " failed to converge on the g polynomial");
1523  }
1524 }
1525 
1527 {
1528  double MC_X_o_min_ = 0.35; // value at the zero solvent point
1529  MC_X_o_cutoff_ = 0.6;
1530  double MC_slopepCut_ = 0.02; // slope of the p function at the zero solvent point
1531  MC_cpCut_ = 0.25;
1532 
1533  // Initial starting values
1534  MC_apCut_ = MC_X_o_min_;
1535  MC_epCut_ = 0.0;
1536  bool converged = false;
1537  double oldV = 0.0;
1538  double damp = 0.5;
1539  for (int its = 0; its < 500 && !converged; its++) {
1540  oldV = MC_epCut_;
1541  MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
1542  double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;
1543  MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
1544  double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ * MC_X_o_cutoff_/MC_cpCut_)
1545  /
1546  (MC_X_o_cutoff_ * MC_X_o_cutoff_/MC_cpCut_ - 2.0 * MC_X_o_cutoff_));
1547  MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
1548  double tmp = MC_apCut_ + MC_X_o_cutoff_*(MC_bpCut_ + MC_dpCut_ * MC_X_o_cutoff_);
1549  double eterm = std::exp(- MC_X_o_cutoff_ / MC_cpCut_);
1550  MC_epCut_ = - eterm * tmp;
1551  double diff = MC_epCut_ - oldV;
1552  if (fabs(diff) < 1.0E-14) {
1553  converged = true;
1554  }
1555  }
1556  if (!converged) {
1557  throw CanteraError("HMWSoln::calcMCCutoffParams_()",
1558  " failed to converge on the p polynomial");
1559  }
1560 }
1561 
1562 void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
1563 {
1564  double T = temperature();
1565  const double twoT = 2.0 * T;
1566  const double invT = 1.0 / T;
1567  const double invT2 = invT * invT;
1568  const double twoinvT3 = 2.0 * invT * invT2;
1569  double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1570  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1571  tlin = T - m_TempPitzerRef;
1572  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1573  tlin = T - m_TempPitzerRef;
1574  tquad = T * T - m_TempPitzerRef * m_TempPitzerRef;
1575  tln = log(T/ m_TempPitzerRef);
1576  tinv = 1.0/T - 1.0/m_TempPitzerRef;
1577  }
1578 
1579  for (size_t i = 1; i < (m_kk - 1); i++) {
1580  for (size_t j = (i+1); j < m_kk; j++) {
1581  // Find the counterIJ for the symmetric binary interaction
1582  size_t n = m_kk*i + j;
1583  size_t counterIJ = m_CounterIJ[n];
1584 
1585  const double* beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ);
1586  const double* beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ);
1587  const double* beta2MX_coeff = m_Beta2MX_ij_coeff.ptrColumn(counterIJ);
1588  const double* CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ);
1589  const double* Theta_coeff = m_Theta_ij_coeff.ptrColumn(counterIJ);
1590 
1591  switch (m_formPitzerTemp) {
1592  case PITZER_TEMP_CONSTANT:
1593  break;
1594  case PITZER_TEMP_LINEAR:
1595 
1596  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1597  + beta0MX_coeff[1]*tlin;
1598  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1];
1599  m_Beta0MX_ij_LL[counterIJ] = 0.0;
1600  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1601  + beta1MX_coeff[1]*tlin;
1602  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1];
1603  m_Beta1MX_ij_LL[counterIJ] = 0.0;
1604  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1605  + beta2MX_coeff[1]*tlin;
1606  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1];
1607  m_Beta2MX_ij_LL[counterIJ] = 0.0;
1608  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1609  + CphiMX_coeff[1]*tlin;
1610  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1];
1611  m_CphiMX_ij_LL[counterIJ] = 0.0;
1612  m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1613  m_Theta_ij_L[counterIJ] = Theta_coeff[1];
1614  m_Theta_ij_LL[counterIJ] = 0.0;
1615  break;
1616 
1617  case PITZER_TEMP_COMPLEX1:
1618  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1619  + beta0MX_coeff[1]*tlin
1620  + beta0MX_coeff[2]*tquad
1621  + beta0MX_coeff[3]*tinv
1622  + beta0MX_coeff[4]*tln;
1623  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1624  + beta1MX_coeff[1]*tlin
1625  + beta1MX_coeff[2]*tquad
1626  + beta1MX_coeff[3]*tinv
1627  + beta1MX_coeff[4]*tln;
1628  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1629  + beta2MX_coeff[1]*tlin
1630  + beta2MX_coeff[2]*tquad
1631  + beta2MX_coeff[3]*tinv
1632  + beta2MX_coeff[4]*tln;
1633  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1634  + CphiMX_coeff[1]*tlin
1635  + CphiMX_coeff[2]*tquad
1636  + CphiMX_coeff[3]*tinv
1637  + CphiMX_coeff[4]*tln;
1638  m_Theta_ij[counterIJ] = Theta_coeff[0]
1639  + Theta_coeff[1]*tlin
1640  + Theta_coeff[2]*tquad
1641  + Theta_coeff[3]*tinv
1642  + Theta_coeff[4]*tln;
1643  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]
1644  + beta0MX_coeff[2]*twoT
1645  - beta0MX_coeff[3]*invT2
1646  + beta0MX_coeff[4]*invT;
1647  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]
1648  + beta1MX_coeff[2]*twoT
1649  - beta1MX_coeff[3]*invT2
1650  + beta1MX_coeff[4]*invT;
1651  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1]
1652  + beta2MX_coeff[2]*twoT
1653  - beta2MX_coeff[3]*invT2
1654  + beta2MX_coeff[4]*invT;
1655  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]
1656  + CphiMX_coeff[2]*twoT
1657  - CphiMX_coeff[3]*invT2
1658  + CphiMX_coeff[4]*invT;
1659  m_Theta_ij_L[counterIJ] = Theta_coeff[1]
1660  + Theta_coeff[2]*twoT
1661  - Theta_coeff[3]*invT2
1662  + Theta_coeff[4]*invT;
1663  doDerivs = 2;
1664  if (doDerivs > 1) {
1665  m_Beta0MX_ij_LL[counterIJ] =
1666  + beta0MX_coeff[2]*2.0
1667  + beta0MX_coeff[3]*twoinvT3
1668  - beta0MX_coeff[4]*invT2;
1669  m_Beta1MX_ij_LL[counterIJ] =
1670  + beta1MX_coeff[2]*2.0
1671  + beta1MX_coeff[3]*twoinvT3
1672  - beta1MX_coeff[4]*invT2;
1673  m_Beta2MX_ij_LL[counterIJ] =
1674  + beta2MX_coeff[2]*2.0
1675  + beta2MX_coeff[3]*twoinvT3
1676  - beta2MX_coeff[4]*invT2;
1677  m_CphiMX_ij_LL[counterIJ] =
1678  + CphiMX_coeff[2]*2.0
1679  + CphiMX_coeff[3]*twoinvT3
1680  - CphiMX_coeff[4]*invT2;
1681  m_Theta_ij_LL[counterIJ] =
1682  + Theta_coeff[2]*2.0
1683  + Theta_coeff[3]*twoinvT3
1684  - Theta_coeff[4]*invT2;
1685  }
1686  break;
1687  }
1688  }
1689  }
1690 
1691  // Lambda interactions and Mu_nnn
1692  // i must be neutral for this term to be nonzero. We take advantage of this
1693  // here to lower the operation count.
1694  for (size_t i = 1; i < m_kk; i++) {
1695  if (charge(i) == 0.0) {
1696  for (size_t j = 1; j < m_kk; j++) {
1697  size_t n = i * m_kk + j;
1698  const double* Lambda_coeff = m_Lambda_nj_coeff.ptrColumn(n);
1699  switch (m_formPitzerTemp) {
1700  case PITZER_TEMP_CONSTANT:
1701  m_Lambda_nj(i,j) = Lambda_coeff[0];
1702  break;
1703  case PITZER_TEMP_LINEAR:
1704  m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1705  m_Lambda_nj_L(i,j) = Lambda_coeff[1];
1706  m_Lambda_nj_LL(i,j) = 0.0;
1707  break;
1708  case PITZER_TEMP_COMPLEX1:
1709  m_Lambda_nj(i,j) = Lambda_coeff[0]
1710  + Lambda_coeff[1]*tlin
1711  + Lambda_coeff[2]*tquad
1712  + Lambda_coeff[3]*tinv
1713  + Lambda_coeff[4]*tln;
1714 
1715  m_Lambda_nj_L(i,j) = Lambda_coeff[1]
1716  + Lambda_coeff[2]*twoT
1717  - Lambda_coeff[3]*invT2
1718  + Lambda_coeff[4]*invT;
1719 
1720  m_Lambda_nj_LL(i,j) =
1721  Lambda_coeff[2]*2.0
1722  + Lambda_coeff[3]*twoinvT3
1723  - Lambda_coeff[4]*invT2;
1724  }
1725 
1726  if (j == i) {
1727  const double* Mu_coeff = m_Mu_nnn_coeff.ptrColumn(i);
1728  switch (m_formPitzerTemp) {
1729  case PITZER_TEMP_CONSTANT:
1730  m_Mu_nnn[i] = Mu_coeff[0];
1731  break;
1732  case PITZER_TEMP_LINEAR:
1733  m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1734  m_Mu_nnn_L[i] = Mu_coeff[1];
1735  m_Mu_nnn_LL[i] = 0.0;
1736  break;
1737  case PITZER_TEMP_COMPLEX1:
1738  m_Mu_nnn[i] = Mu_coeff[0]
1739  + Mu_coeff[1]*tlin
1740  + Mu_coeff[2]*tquad
1741  + Mu_coeff[3]*tinv
1742  + Mu_coeff[4]*tln;
1743  m_Mu_nnn_L[i] = Mu_coeff[1]
1744  + Mu_coeff[2]*twoT
1745  - Mu_coeff[3]*invT2
1746  + Mu_coeff[4]*invT;
1747  m_Mu_nnn_LL[i] =
1748  Mu_coeff[2]*2.0
1749  + Mu_coeff[3]*twoinvT3
1750  - Mu_coeff[4]*invT2;
1751  }
1752  }
1753  }
1754  }
1755  }
1756 
1757  switch(m_formPitzerTemp) {
1758  case PITZER_TEMP_CONSTANT:
1759  for (size_t i = 1; i < m_kk; i++) {
1760  for (size_t j = 1; j < m_kk; j++) {
1761  for (size_t k = 1; k < m_kk; k++) {
1762  size_t n = i * m_kk *m_kk + j * m_kk + k;
1763  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1764  m_Psi_ijk[n] = Psi_coeff[0];
1765  }
1766  }
1767  }
1768  break;
1769  case PITZER_TEMP_LINEAR:
1770  for (size_t i = 1; i < m_kk; i++) {
1771  for (size_t j = 1; j < m_kk; j++) {
1772  for (size_t k = 1; k < m_kk; k++) {
1773  size_t n = i * m_kk *m_kk + j * m_kk + k;
1774  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1775  m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1776  m_Psi_ijk_L[n] = Psi_coeff[1];
1777  m_Psi_ijk_LL[n] = 0.0;
1778  }
1779  }
1780  }
1781  break;
1782  case PITZER_TEMP_COMPLEX1:
1783  for (size_t i = 1; i < m_kk; i++) {
1784  for (size_t j = 1; j < m_kk; j++) {
1785  for (size_t k = 1; k < m_kk; k++) {
1786  size_t n = i * m_kk *m_kk + j * m_kk + k;
1787  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1788  m_Psi_ijk[n] = Psi_coeff[0]
1789  + Psi_coeff[1]*tlin
1790  + Psi_coeff[2]*tquad
1791  + Psi_coeff[3]*tinv
1792  + Psi_coeff[4]*tln;
1793  m_Psi_ijk_L[n] = Psi_coeff[1]
1794  + Psi_coeff[2]*twoT
1795  - Psi_coeff[3]*invT2
1796  + Psi_coeff[4]*invT;
1797  m_Psi_ijk_LL[n] =
1798  Psi_coeff[2]*2.0
1799  + Psi_coeff[3]*twoinvT3
1800  - Psi_coeff[4]*invT2;
1801  }
1802  }
1803  }
1804  break;
1805  }
1806 }
1807 
1809 {
1810  // Use the CROPPED molality of the species in solution.
1811  const vector<double>& molality = m_molalitiesCropped;
1812 
1813  // These are data inputs about the Pitzer correlation. They come from the
1814  // input file for the Pitzer model.
1815  vector<double>& gamma_Unscaled = m_gamma_tmp;
1816 
1817  // Local variables defined by Coltrin
1818  double etheta[5][5], etheta_prime[5][5], sqrtIs;
1819 
1820  // Molality based ionic strength of the solution
1821  double Is = 0.0;
1822 
1823  // Molar charge of the solution: In Pitzer's notation, this is his variable
1824  // called "Z".
1825  double molarcharge = 0.0;
1826 
1827  // molalitysum is the sum of the molalities over all solutes, even those
1828  // with zero charge.
1829  double molalitysumUncropped = 0.0;
1830 
1831  // Make sure the counter variables are setup
1832  counterIJ_setup();
1833 
1834  // ---------- Calculate common sums over solutes ---------------------
1835  for (size_t n = 1; n < m_kk; n++) {
1836  // ionic strength
1837  Is += charge(n) * charge(n) * molality[n];
1838  // total molar charge
1839  molarcharge += fabs(charge(n)) * molality[n];
1840  molalitysumUncropped += m_molalities[n];
1841  }
1842  Is *= 0.5;
1843 
1844  // Store the ionic molality in the object for reference.
1845  m_IionicMolality = Is;
1846  sqrtIs = sqrt(Is);
1847 
1848  // The following call to calc_lambdas() calculates all 16 elements of the
1849  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
1850  calc_lambdas(Is);
1851 
1852  // Step 2: Find the coefficients E-theta and E-thetaprime for all
1853  // combinations of positive unlike charges up to 4
1854  for (int z1 = 1; z1 <=4; z1++) {
1855  for (int z2 =1; z2 <=4; z2++) {
1856  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
1857  }
1858  }
1859 
1860  // calculate g(x) and hfunc(x) for each cation-anion pair MX. In the
1861  // original literature, hfunc, was called gprime. However, it's not the
1862  // derivative of g(x), so I renamed it.
1863  for (size_t i = 1; i < (m_kk - 1); i++) {
1864  for (size_t j = (i+1); j < m_kk; j++) {
1865  // Find the counterIJ for the symmetric binary interaction
1866  size_t n = m_kk*i + j;
1867  size_t counterIJ = m_CounterIJ[n];
1868 
1869  // Only loop over oppositely charge species
1870  if (charge(i)*charge(j) < 0) {
1871  // x is a reduced function variable
1872  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
1873  if (x1 > 1.0E-100) {
1874  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
1875  m_hfunc_IJ[counterIJ] = -2.0 *
1876  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
1877  } else {
1878  m_gfunc_IJ[counterIJ] = 0.0;
1879  m_hfunc_IJ[counterIJ] = 0.0;
1880  }
1881 
1882  if (m_Beta2MX_ij[counterIJ] != 0.0) {
1883  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
1884  if (x2 > 1.0E-100) {
1885  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
1886  m_h2func_IJ[counterIJ] = -2.0 *
1887  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
1888  } else {
1889  m_g2func_IJ[counterIJ] = 0.0;
1890  m_h2func_IJ[counterIJ] = 0.0;
1891  }
1892  }
1893  } else {
1894  m_gfunc_IJ[counterIJ] = 0.0;
1895  m_hfunc_IJ[counterIJ] = 0.0;
1896  }
1897  }
1898  }
1899 
1900  // SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX
1901  // Agrees with Pitzer, Eq. (49), (51), (55)
1902  for (size_t i = 1; i < m_kk - 1; i++) {
1903  for (size_t j = i+1; j < m_kk; j++) {
1904  // Find the counterIJ for the symmetric binary interaction
1905  size_t n = m_kk*i + j;
1906  size_t counterIJ = m_CounterIJ[n];
1907 
1908  // both species have a non-zero charge, and one is positive and the
1909  // other is negative
1910  if (charge(i)*charge(j) < 0.0) {
1911  m_BMX_IJ[counterIJ] = m_Beta0MX_ij[counterIJ]
1912  + m_Beta1MX_ij[counterIJ] * m_gfunc_IJ[counterIJ]
1913  + m_Beta2MX_ij[counterIJ] * m_g2func_IJ[counterIJ];
1914 
1915  if (Is > 1.0E-150) {
1916  m_BprimeMX_IJ[counterIJ] = (m_Beta1MX_ij[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
1917  m_Beta2MX_ij[counterIJ] * m_h2func_IJ[counterIJ]/Is);
1918  } else {
1919  m_BprimeMX_IJ[counterIJ] = 0.0;
1920  }
1921  m_BphiMX_IJ[counterIJ] = m_BMX_IJ[counterIJ] + Is*m_BprimeMX_IJ[counterIJ];
1922  } else {
1923  m_BMX_IJ[counterIJ] = 0.0;
1924  m_BprimeMX_IJ[counterIJ] = 0.0;
1925  m_BphiMX_IJ[counterIJ] = 0.0;
1926  }
1927  }
1928  }
1929 
1930  // SUBSECTION TO CALCULATE CMX
1931  // Agrees with Pitzer, Eq. (53).
1932  for (size_t i = 1; i < m_kk-1; i++) {
1933  for (size_t j = i+1; j < m_kk; j++) {
1934  // Find the counterIJ for the symmetric binary interaction
1935  size_t n = m_kk*i + j;
1936  size_t counterIJ = m_CounterIJ[n];
1937 
1938  // both species have a non-zero charge, and one is positive
1939  // and the other is negative
1940  if (charge(i)*charge(j) < 0.0) {
1941  m_CMX_IJ[counterIJ] = m_CphiMX_ij[counterIJ]/
1942  (2.0* sqrt(fabs(charge(i)*charge(j))));
1943  } else {
1944  m_CMX_IJ[counterIJ] = 0.0;
1945  }
1946  }
1947  }
1948 
1949  // SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi
1950  // Agrees with Pitzer, Eq. 72, 73, 74
1951  for (size_t i = 1; i < m_kk-1; i++) {
1952  for (size_t j = i+1; j < m_kk; j++) {
1953  // Find the counterIJ for the symmetric binary interaction
1954  size_t n = m_kk*i + j;
1955  size_t counterIJ = m_CounterIJ[n];
1956 
1957  // both species have a non-zero charge, and one is positive and the
1958  // other is negative
1959  if (charge(i)*charge(j) > 0) {
1960  int z1 = (int) fabs(charge(i));
1961  int z2 = (int) fabs(charge(j));
1962  m_Phi_IJ[counterIJ] = m_Theta_ij[counterIJ] + etheta[z1][z2];
1963  m_Phiprime_IJ[counterIJ] = etheta_prime[z1][z2];
1964  m_PhiPhi_IJ[counterIJ] = m_Phi_IJ[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
1965  } else {
1966  m_Phi_IJ[counterIJ] = 0.0;
1967  m_Phiprime_IJ[counterIJ] = 0.0;
1968  m_PhiPhi_IJ[counterIJ] = 0.0;
1969  }
1970  }
1971  }
1972 
1973  // SUBSECTION FOR CALCULATION OF F
1974  // Agrees with Pitzer Eqn. (65)
1975  double Aphi = A_Debye_TP() / 3.0;
1976  double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
1977  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
1978  for (size_t i = 1; i < m_kk-1; i++) {
1979  for (size_t j = i+1; j < m_kk; j++) {
1980  // Find the counterIJ for the symmetric binary interaction
1981  size_t n = m_kk*i + j;
1982  size_t counterIJ = m_CounterIJ[n];
1983 
1984  // both species have a non-zero charge, and one is positive and the
1985  // other is negative
1986  if (charge(i)*charge(j) < 0) {
1987  F += molality[i]*molality[j] * m_BprimeMX_IJ[counterIJ];
1988  }
1989 
1990  // Both species have a non-zero charge, and they
1991  // have the same sign
1992  if (charge(i)*charge(j) > 0) {
1993  F += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
1994  }
1995  }
1996  }
1997 
1998  for (size_t i = 1; i < m_kk; i++) {
1999 
2000  // SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS
2001  // equations agree with my notes, Eqn. (118).
2002  // Equations agree with Pitzer, eqn.(63)
2003  if (charge(i) > 0.0) {
2004  // species i is the cation (positive) to calc the actcoeff
2005  double zsqF = charge(i)*charge(i)*F;
2006  double sum1 = 0.0;
2007  double sum2 = 0.0;
2008  double sum3 = 0.0;
2009  double sum4 = 0.0;
2010  double sum5 = 0.0;
2011  for (size_t j = 1; j < m_kk; j++) {
2012  // Find the counterIJ for the symmetric binary interaction
2013  size_t n = m_kk*i + j;
2014  size_t counterIJ = m_CounterIJ[n];
2015 
2016  if (charge(j) < 0.0) {
2017  // sum over all anions
2018  sum1 += molality[j] *
2019  (2.0*m_BMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
2020  if (j < m_kk-1) {
2021  // This term is the ternary interaction involving the
2022  // non-duplicate sum over double anions, j, k, with
2023  // respect to the cation, i.
2024  for (size_t k = j+1; k < m_kk; k++) {
2025  // an inner sum over all anions
2026  if (charge(k) < 0.0) {
2027  n = k + j * m_kk + i * m_kk * m_kk;
2028  sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
2029  }
2030  }
2031  }
2032  }
2033 
2034  if (charge(j) > 0.0) {
2035  // sum over all cations
2036  if (j != i) {
2037  sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
2038  }
2039  for (size_t k = 1; k < m_kk; k++) {
2040  if (charge(k) < 0.0) {
2041  // two inner sums over anions
2042  n = k + j * m_kk + i * m_kk * m_kk;
2043  sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
2044 
2045  // Find the counterIJ for the j,k interaction
2046  n = m_kk*j + k;
2047  size_t counterIJ2 = m_CounterIJ[n];
2048  sum4 += (fabs(charge(i))*
2049  molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
2050  }
2051  }
2052  }
2053 
2054  // Handle neutral j species
2055  if (charge(j) == 0) {
2056  sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
2057 
2058  // Zeta interaction term
2059  for (size_t k = 1; k < m_kk; k++) {
2060  if (charge(k) < 0.0) {
2061  size_t izeta = j;
2062  size_t jzeta = i;
2063  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2064  double zeta = m_Psi_ijk[n];
2065  if (zeta != 0.0) {
2066  sum5 += molality[j]*molality[k]*zeta;
2067  }
2068  }
2069  }
2070  }
2071  }
2072 
2073  // Add all of the contributions up to yield the log of the solute
2074  // activity coefficients (molality scale)
2075  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2076  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2077  }
2078 
2079  // SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS
2080  // equations agree with my notes, Eqn. (119).
2081  // Equations agree with Pitzer, eqn.(64)
2082  if (charge(i) < 0) {
2083  // species i is an anion (negative)
2084  double zsqF = charge(i)*charge(i)*F;
2085  double sum1 = 0.0;
2086  double sum2 = 0.0;
2087  double sum3 = 0.0;
2088  double sum4 = 0.0;
2089  double sum5 = 0.0;
2090  for (size_t j = 1; j < m_kk; j++) {
2091  // Find the counterIJ for the symmetric binary interaction
2092  size_t n = m_kk*i + j;
2093  size_t counterIJ = m_CounterIJ[n];
2094 
2095  // For Anions, do the cation interactions.
2096  if (charge(j) > 0) {
2097  sum1 += molality[j]*
2098  (2.0*m_BMX_IJ[counterIJ]+molarcharge*m_CMX_IJ[counterIJ]);
2099  if (j < m_kk-1) {
2100  for (size_t k = j+1; k < m_kk; k++) {
2101  // an inner sum over all cations
2102  if (charge(k) > 0) {
2103  n = k + j * m_kk + i * m_kk * m_kk;
2104  sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
2105  }
2106  }
2107  }
2108  }
2109 
2110  // For Anions, do the other anion interactions.
2111  if (charge(j) < 0.0) {
2112  // sum over all anions
2113  if (j != i) {
2114  sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
2115  }
2116  for (size_t k = 1; k < m_kk; k++) {
2117  if (charge(k) > 0.0) {
2118  // two inner sums over cations
2119  n = k + j * m_kk + i * m_kk * m_kk;
2120  sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
2121  // Find the counterIJ for the symmetric binary interaction
2122  n = m_kk*j + k;
2123  size_t counterIJ2 = m_CounterIJ[n];
2124  sum4 += fabs(charge(i))*
2125  molality[j]*molality[k]*m_CMX_IJ[counterIJ2];
2126  }
2127  }
2128  }
2129 
2130  // for Anions, do the neutral species interaction
2131  if (charge(j) == 0.0) {
2132  sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
2133  // Zeta interaction term
2134  for (size_t k = 1; k < m_kk; k++) {
2135  if (charge(k) > 0.0) {
2136  size_t izeta = j;
2137  size_t jzeta = k;
2138  size_t kzeta = i;
2139  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2140  double zeta = m_Psi_ijk[n];
2141  if (zeta != 0.0) {
2142  sum5 += molality[j]*molality[k]*zeta;
2143  }
2144  }
2145  }
2146  }
2147  }
2148  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2149  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2150  }
2151 
2152  // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
2153  // equations agree with my notes,
2154  // Equations agree with Pitzer,
2155  if (charge(i) == 0.0) {
2156  double sum1 = 0.0;
2157  double sum3 = 0.0;
2158  for (size_t j = 1; j < m_kk; j++) {
2159  sum1 += molality[j]*2.0*m_Lambda_nj(i,j);
2160  // Zeta term -> we piggyback on the psi term
2161  if (charge(j) > 0.0) {
2162  for (size_t k = 1; k < m_kk; k++) {
2163  if (charge(k) < 0.0) {
2164  size_t n = k + j * m_kk + i * m_kk * m_kk;
2165  sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
2166  }
2167  }
2168  }
2169  }
2170  double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
2171  m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
2172  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2173  }
2174  }
2175 
2176  // SUBSECTION FOR CALCULATING THE OSMOTIC COEFF
2177  // equations agree with my notes, Eqn. (117).
2178  // Equations agree with Pitzer, eqn.(62)
2179  double sum1 = 0.0;
2180  double sum2 = 0.0;
2181  double sum3 = 0.0;
2182  double sum4 = 0.0;
2183  double sum5 = 0.0;
2184  double sum6 = 0.0;
2185  double sum7 = 0.0;
2186 
2187  // term1 is the DH term in the osmotic coefficient expression
2188  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
2189  // implementations.
2190  // Is = Ionic strength on the molality scale (units of (gmol/kg))
2191  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2192  double term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2193 
2194  for (size_t j = 1; j < m_kk; j++) {
2195  // Loop Over Cations
2196  if (charge(j) > 0.0) {
2197  for (size_t k = 1; k < m_kk; k++) {
2198  if (charge(k) < 0.0) {
2199  // Find the counterIJ for the symmetric j,k binary interaction
2200  size_t n = m_kk*j + k;
2201  size_t counterIJ = m_CounterIJ[n];
2202 
2203  sum1 += molality[j]*molality[k]*
2204  (m_BphiMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
2205  }
2206  }
2207 
2208  for (size_t k = j+1; k < m_kk; k++) {
2209  if (j == (m_kk-1)) {
2210  // we should never reach this step
2211  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2212  "logic error 1 in Step 9 of hmw_act");
2213  }
2214  if (charge(k) > 0.0) {
2215  // Find the counterIJ for the symmetric j,k binary interaction
2216  // between 2 cations.
2217  size_t n = m_kk*j + k;
2218  size_t counterIJ = m_CounterIJ[n];
2219  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ[counterIJ];
2220  for (size_t m = 1; m < m_kk; m++) {
2221  if (charge(m) < 0.0) {
2222  // species m is an anion
2223  n = m + k * m_kk + j * m_kk * m_kk;
2224  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk[n];
2225  }
2226  }
2227  }
2228  }
2229  }
2230 
2231  // Loop Over Anions
2232  if (charge(j) < 0) {
2233  for (size_t k = j+1; k < m_kk; k++) {
2234  if (j == m_kk-1) {
2235  // we should never reach this step
2236  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2237  "logic error 2 in Step 9 of hmw_act");
2238  }
2239  if (charge(k) < 0) {
2240  // Find the counterIJ for the symmetric j,k binary interaction
2241  // between two anions
2242  size_t n = m_kk*j + k;
2243  size_t counterIJ = m_CounterIJ[n];
2244  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ[counterIJ];
2245  for (size_t m = 1; m < m_kk; m++) {
2246  if (charge(m) > 0.0) {
2247  n = m + k * m_kk + j * m_kk * m_kk;
2248  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk[n];
2249  }
2250  }
2251  }
2252  }
2253  }
2254 
2255  // Loop Over Neutral Species
2256  if (charge(j) == 0) {
2257  for (size_t k = 1; k < m_kk; k++) {
2258  if (charge(k) < 0.0) {
2259  sum4 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2260  }
2261  if (charge(k) > 0.0) {
2262  sum5 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2263  }
2264  if (charge(k) == 0.0) {
2265  if (k > j) {
2266  sum6 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2267  } else if (k == j) {
2268  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj(j,k);
2269  }
2270  }
2271  if (charge(k) < 0.0) {
2272  size_t izeta = j;
2273  for (size_t m = 1; m < m_kk; m++) {
2274  if (charge(m) > 0.0) {
2275  size_t jzeta = m;
2276  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2277  double zeta = m_Psi_ijk[n];
2278  if (zeta != 0.0) {
2279  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2280  }
2281  }
2282  }
2283  }
2284  }
2285  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j];
2286  }
2287  }
2288  double sum_m_phi_minus_1 = 2.0 *
2289  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2290  // Calculate the osmotic coefficient from
2291  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2292  double osmotic_coef;
2293  if (molalitysumUncropped > 1.0E-150) {
2294  osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2295  } else {
2296  osmotic_coef = 1.0;
2297  }
2298  double lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2299 
2300  // In Cantera, we define the activity coefficient of the solvent as
2301  //
2302  // act_0 = actcoeff_0 * Xmol_0
2303  //
2304  // We have just computed act_0. However, this routine returns
2305  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2306  double xmolSolvent = moleFraction(0);
2307  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
2308  m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
2309 }
2310 
2312 {
2313  static const int cacheId = m_cache.getId();
2314  CachedScalar cached = m_cache.getScalar(cacheId);
2315  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2316  return;
2317  }
2318 
2319  // Zero the unscaled 2nd derivatives
2320  m_dlnActCoeffMolaldT_Unscaled.assign(m_kk, 0.0);
2321 
2322  // Do the actual calculation of the unscaled temperature derivatives
2324 
2325  for (size_t k = 1; k < m_kk; k++) {
2326  if (CROP_speciesCropped_[k] == 2) {
2328  }
2329  }
2330 
2331  if (CROP_speciesCropped_[0]) {
2333  }
2334 
2335  // Do the pH scaling to the derivatives
2337 }
2338 
2340 {
2341  // It may be assumed that the Pitzer activity coefficient routine is called
2342  // immediately preceding the calling of this routine. Therefore, some
2343  // quantities do not need to be recalculated in this routine.
2344 
2345  const vector<double>& molality = m_molalitiesCropped;
2346  double* d_gamma_dT_Unscaled = m_gamma_tmp.data();
2347 
2348  // Local variables defined by Coltrin
2349  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2350 
2351  // Molality based ionic strength of the solution
2352  double Is = 0.0;
2353 
2354  // Molar charge of the solution: In Pitzer's notation, this is his variable
2355  // called "Z".
2356  double molarcharge = 0.0;
2357 
2358  // molalitysum is the sum of the molalities over all solutes, even those
2359  // with zero charge.
2360  double molalitysum = 0.0;
2361 
2362  // Make sure the counter variables are setup
2363  counterIJ_setup();
2364 
2365  // ---------- Calculate common sums over solutes ---------------------
2366  for (size_t n = 1; n < m_kk; n++) {
2367  // ionic strength
2368  Is += charge(n) * charge(n) * molality[n];
2369  // total molar charge
2370  molarcharge += fabs(charge(n)) * molality[n];
2371  molalitysum += molality[n];
2372  }
2373  Is *= 0.5;
2374 
2375  // Store the ionic molality in the object for reference.
2376  m_IionicMolality = Is;
2377  sqrtIs = sqrt(Is);
2378 
2379  // The following call to calc_lambdas() calculates all 16 elements of the
2380  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
2381  calc_lambdas(Is);
2382 
2383  // Step 2: Find the coefficients E-theta and E-thetaprime for all
2384  // combinations of positive unlike charges up to 4
2385  for (int z1 = 1; z1 <=4; z1++) {
2386  for (int z2 =1; z2 <=4; z2++) {
2387  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2388  }
2389  }
2390 
2391  // calculate g(x) and hfunc(x) for each cation-anion pair MX
2392  // In the original literature, hfunc, was called gprime. However,
2393  // it's not the derivative of g(x), so I renamed it.
2394  for (size_t i = 1; i < (m_kk - 1); i++) {
2395  for (size_t j = (i+1); j < m_kk; j++) {
2396  // Find the counterIJ for the symmetric binary interaction
2397  size_t n = m_kk*i + j;
2398  size_t counterIJ = m_CounterIJ[n];
2399 
2400  // Only loop over oppositely charge species
2401  if (charge(i)*charge(j) < 0) {
2402  // x is a reduced function variable
2403  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2404  if (x1 > 1.0E-100) {
2405  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2406  m_hfunc_IJ[counterIJ] = -2.0 *
2407  (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
2408  } else {
2409  m_gfunc_IJ[counterIJ] = 0.0;
2410  m_hfunc_IJ[counterIJ] = 0.0;
2411  }
2412 
2413  if (m_Beta2MX_ij_L[counterIJ] != 0.0) {
2414  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
2415  if (x2 > 1.0E-100) {
2416  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2417  m_h2func_IJ[counterIJ] = -2.0 *
2418  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2419  } else {
2420  m_g2func_IJ[counterIJ] = 0.0;
2421  m_h2func_IJ[counterIJ] = 0.0;
2422  }
2423  }
2424  } else {
2425  m_gfunc_IJ[counterIJ] = 0.0;
2426  m_hfunc_IJ[counterIJ] = 0.0;
2427  }
2428  }
2429  }
2430 
2431  // SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L
2432  // These are now temperature derivatives of the previously calculated
2433  // quantities.
2434  for (size_t i = 1; i < m_kk - 1; i++) {
2435  for (size_t j = i+1; j < m_kk; j++) {
2436  // Find the counterIJ for the symmetric binary interaction
2437  size_t n = m_kk*i + j;
2438  size_t counterIJ = m_CounterIJ[n];
2439 
2440  // both species have a non-zero charge, and one is positive
2441  // and the other is negative
2442  if (charge(i)*charge(j) < 0.0) {
2443  m_BMX_IJ_L[counterIJ] = m_Beta0MX_ij_L[counterIJ]
2444  + m_Beta1MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ]
2445  + m_Beta2MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ];
2446  if (Is > 1.0E-150) {
2447  m_BprimeMX_IJ_L[counterIJ] = (m_Beta1MX_ij_L[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
2448  m_Beta2MX_ij_L[counterIJ] * m_h2func_IJ[counterIJ]/Is);
2449  } else {
2450  m_BprimeMX_IJ_L[counterIJ] = 0.0;
2451  }
2452  m_BphiMX_IJ_L[counterIJ] = m_BMX_IJ_L[counterIJ] + Is*m_BprimeMX_IJ_L[counterIJ];
2453  } else {
2454  m_BMX_IJ_L[counterIJ] = 0.0;
2455  m_BprimeMX_IJ_L[counterIJ] = 0.0;
2456  m_BphiMX_IJ_L[counterIJ] = 0.0;
2457  }
2458  }
2459  }
2460 
2461  // --------- SUBSECTION TO CALCULATE CMX_L ----------
2462  for (size_t i = 1; i < m_kk-1; i++) {
2463  for (size_t j = i+1; j < m_kk; j++) {
2464  // Find the counterIJ for the symmetric binary interaction
2465  size_t n = m_kk*i + j;
2466  size_t counterIJ = m_CounterIJ[n];
2467 
2468  // both species have a non-zero charge, and one is positive
2469  // and the other is negative
2470  if (charge(i)*charge(j) < 0.0) {
2471  m_CMX_IJ_L[counterIJ] = m_CphiMX_ij_L[counterIJ]/
2472  (2.0* sqrt(fabs(charge(i)*charge(j))));
2473  } else {
2474  m_CMX_IJ_L[counterIJ] = 0.0;
2475  }
2476  }
2477  }
2478 
2479  // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2480  for (size_t i = 1; i < m_kk-1; i++) {
2481  for (size_t j = i+1; j < m_kk; j++) {
2482  // Find the counterIJ for the symmetric binary interaction
2483  size_t n = m_kk*i + j;
2484  size_t counterIJ = m_CounterIJ[n];
2485 
2486  // both species have a non-zero charge, and one is positive
2487  // and the other is negative
2488  if (charge(i)*charge(j) > 0) {
2489  m_Phi_IJ_L[counterIJ] = m_Theta_ij_L[counterIJ];
2490  m_Phiprime_IJ[counterIJ] = 0.0;
2491  m_PhiPhi_IJ_L[counterIJ] = m_Phi_IJ_L[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
2492  } else {
2493  m_Phi_IJ_L[counterIJ] = 0.0;
2494  m_Phiprime_IJ[counterIJ] = 0.0;
2495  m_PhiPhi_IJ_L[counterIJ] = 0.0;
2496  }
2497  }
2498  }
2499 
2500  // ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
2501  double dA_DebyedT = dA_DebyedT_TP();
2502  double dAphidT = dA_DebyedT /3.0;
2503  double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2504  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2505  for (size_t i = 1; i < m_kk-1; i++) {
2506  for (size_t j = i+1; j < m_kk; j++) {
2507  // Find the counterIJ for the symmetric binary interaction
2508  size_t n = m_kk*i + j;
2509  size_t counterIJ = m_CounterIJ[n];
2510 
2511  // both species have a non-zero charge, and one is positive
2512  // and the other is negative
2513  if (charge(i)*charge(j) < 0) {
2514  dFdT += molality[i]*molality[j] * m_BprimeMX_IJ_L[counterIJ];
2515  }
2516 
2517  // Both species have a non-zero charge, and they
2518  // have the same sign, that is, both positive or both negative.
2519  if (charge(i)*charge(j) > 0) {
2520  dFdT += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
2521  }
2522  }
2523  }
2524 
2525  for (size_t i = 1; i < m_kk; i++) {
2526  // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
2527  if (charge(i) > 0) {
2528  // species i is the cation (positive) to calc the actcoeff
2529  double zsqdFdT = charge(i)*charge(i)*dFdT;
2530  double sum1 = 0.0;
2531  double sum2 = 0.0;
2532  double sum3 = 0.0;
2533  double sum4 = 0.0;
2534  double sum5 = 0.0;
2535  for (size_t j = 1; j < m_kk; j++) {
2536  // Find the counterIJ for the symmetric binary interaction
2537  size_t n = m_kk*i + j;
2538  size_t counterIJ = m_CounterIJ[n];
2539 
2540  if (charge(j) < 0.0) {
2541  // sum over all anions
2542  sum1 += molality[j]*
2543  (2.0*m_BMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2544  if (j < m_kk-1) {
2545  // This term is the ternary interaction involving the
2546  // non-duplicate sum over double anions, j, k, with
2547  // respect to the cation, i.
2548  for (size_t k = j+1; k < m_kk; k++) {
2549  // an inner sum over all anions
2550  if (charge(k) < 0.0) {
2551  n = k + j * m_kk + i * m_kk * m_kk;
2552  sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2553  }
2554  }
2555  }
2556  }
2557 
2558  if (charge(j) > 0.0) {
2559  // sum over all cations
2560  if (j != i) {
2561  sum2 += molality[j]*(2.0*m_Phi_IJ_L[counterIJ]);
2562  }
2563  for (size_t k = 1; k < m_kk; k++) {
2564  if (charge(k) < 0.0) {
2565  // two inner sums over anions
2566  n = k + j * m_kk + i * m_kk * m_kk;
2567  sum2 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2568 
2569  // Find the counterIJ for the j,k interaction
2570  n = m_kk*j + k;
2571  size_t counterIJ2 = m_CounterIJ[n];
2572  sum4 += fabs(charge(i))*
2573  molality[j]*molality[k]*m_CMX_IJ_L[counterIJ2];
2574  }
2575  }
2576  }
2577 
2578  // Handle neutral j species
2579  if (charge(j) == 0) {
2580  sum5 += molality[j]*2.0*m_Lambda_nj_L(j,i);
2581  }
2582 
2583  // Zeta interaction term
2584  for (size_t k = 1; k < m_kk; k++) {
2585  if (charge(k) < 0.0) {
2586  size_t izeta = j;
2587  size_t jzeta = i;
2588  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2589  double zeta_L = m_Psi_ijk_L[n];
2590  if (zeta_L != 0.0) {
2591  sum5 += molality[j]*molality[k]*zeta_L;
2592  }
2593  }
2594  }
2595  }
2596 
2597  // Add all of the contributions up to yield the log of the
2598  // solute activity coefficients (molality scale)
2600  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2601  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2602  }
2603 
2604  // ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------
2605  if (charge(i) < 0) {
2606  // species i is an anion (negative)
2607  double zsqdFdT = charge(i)*charge(i)*dFdT;
2608  double sum1 = 0.0;
2609  double sum2 = 0.0;
2610  double sum3 = 0.0;
2611  double sum4 = 0.0;
2612  double sum5 = 0.0;
2613  for (size_t j = 1; j < m_kk; j++) {
2614  // Find the counterIJ for the symmetric binary interaction
2615  size_t n = m_kk*i + j;
2616  size_t counterIJ = m_CounterIJ[n];
2617 
2618  // For Anions, do the cation interactions.
2619  if (charge(j) > 0) {
2620  sum1 += molality[j]*
2621  (2.0*m_BMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2622  if (j < m_kk-1) {
2623  for (size_t k = j+1; k < m_kk; k++) {
2624  // an inner sum over all cations
2625  if (charge(k) > 0) {
2626  n = k + j * m_kk + i * m_kk * m_kk;
2627  sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2628  }
2629  }
2630  }
2631  }
2632 
2633  // For Anions, do the other anion interactions.
2634  if (charge(j) < 0.0) {
2635  // sum over all anions
2636  if (j != i) {
2637  sum2 += molality[j]*(2.0*m_Phi_IJ_L[counterIJ]);
2638  }
2639  for (size_t k = 1; k < m_kk; k++) {
2640  if (charge(k) > 0.0) {
2641  // two inner sums over cations
2642  n = k + j * m_kk + i * m_kk * m_kk;
2643  sum2 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2644  // Find the counterIJ for the symmetric binary interaction
2645  n = m_kk*j + k;
2646  size_t counterIJ2 = m_CounterIJ[n];
2647  sum4 += fabs(charge(i)) *
2648  molality[j]*molality[k]*m_CMX_IJ_L[counterIJ2];
2649  }
2650  }
2651  }
2652 
2653  // for Anions, do the neutral species interaction
2654  if (charge(j) == 0.0) {
2655  sum5 += molality[j]*2.0*m_Lambda_nj_L(j,i);
2656  for (size_t k = 1; k < m_kk; k++) {
2657  if (charge(k) > 0.0) {
2658  size_t izeta = j;
2659  size_t jzeta = k;
2660  size_t kzeta = i;
2661  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2662  double zeta_L = m_Psi_ijk_L[n];
2663  if (zeta_L != 0.0) {
2664  sum5 += molality[j]*molality[k]*zeta_L;
2665  }
2666  }
2667  }
2668  }
2669  }
2671  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2672  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2673  }
2674 
2675  // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
2676  // equations agree with my notes,
2677  // Equations agree with Pitzer,
2678  if (charge(i) == 0.0) {
2679  double sum1 = 0.0;
2680  double sum3 = 0.0;
2681  for (size_t j = 1; j < m_kk; j++) {
2682  sum1 += molality[j]*2.0*m_Lambda_nj_L(i,j);
2683  // Zeta term -> we piggyback on the psi term
2684  if (charge(j) > 0.0) {
2685  for (size_t k = 1; k < m_kk; k++) {
2686  if (charge(k) < 0.0) {
2687  size_t n = k + j * m_kk + i * m_kk * m_kk;
2688  sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2689  }
2690  }
2691  }
2692  }
2693  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
2694  m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
2695  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2696  }
2697  }
2698 
2699  // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
2700  double sum1 = 0.0;
2701  double sum2 = 0.0;
2702  double sum3 = 0.0;
2703  double sum4 = 0.0;
2704  double sum5 = 0.0;
2705  double sum6 = 0.0;
2706  double sum7 = 0.0;
2707 
2708  // term1 is the temperature derivative of the DH term in the osmotic
2709  // coefficient expression
2710  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
2711  // Is = Ionic strength on the molality scale (units of (gmol/kg))
2712  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2713  double term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
2714 
2715  for (size_t j = 1; j < m_kk; j++) {
2716  // Loop Over Cations
2717  if (charge(j) > 0.0) {
2718  for (size_t k = 1; k < m_kk; k++) {
2719  if (charge(k) < 0.0) {
2720  // Find the counterIJ for the symmetric j,k binary interaction
2721  size_t n = m_kk*j + k;
2722  size_t counterIJ = m_CounterIJ[n];
2723  sum1 += molality[j]*molality[k]*
2724  (m_BphiMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2725  }
2726  }
2727 
2728  for (size_t k = j+1; k < m_kk; k++) {
2729  if (j == (m_kk-1)) {
2730  // we should never reach this step
2731  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2732  "logic error 1 in Step 9 of hmw_act");
2733  }
2734  if (charge(k) > 0.0) {
2735  // Find the counterIJ for the symmetric j,k binary interaction
2736  // between 2 cations.
2737  size_t n = m_kk*j + k;
2738  size_t counterIJ = m_CounterIJ[n];
2739  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_L[counterIJ];
2740  for (size_t m = 1; m < m_kk; m++) {
2741  if (charge(m) < 0.0) {
2742  // species m is an anion
2743  n = m + k * m_kk + j * m_kk * m_kk;
2744  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_L[n];
2745  }
2746  }
2747  }
2748  }
2749  }
2750 
2751  // Loop Over Anions
2752  if (charge(j) < 0) {
2753  for (size_t k = j+1; k < m_kk; k++) {
2754  if (j == m_kk-1) {
2755  // we should never reach this step
2756  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2757  "logic error 2 in Step 9 of hmw_act");
2758  }
2759  if (charge(k) < 0) {
2760  // Find the counterIJ for the symmetric j,k binary interaction
2761  // between two anions
2762  size_t n = m_kk*j + k;
2763  size_t counterIJ = m_CounterIJ[n];
2764  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_L[counterIJ];
2765  for (size_t m = 1; m < m_kk; m++) {
2766  if (charge(m) > 0.0) {
2767  n = m + k * m_kk + j * m_kk * m_kk;
2768  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_L[n];
2769  }
2770  }
2771  }
2772  }
2773  }
2774 
2775  // Loop Over Neutral Species
2776  if (charge(j) == 0) {
2777  for (size_t k = 1; k < m_kk; k++) {
2778  if (charge(k) < 0.0) {
2779  sum4 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2780  }
2781  if (charge(k) > 0.0) {
2782  sum5 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2783  }
2784  if (charge(k) == 0.0) {
2785  if (k > j) {
2786  sum6 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2787  } else if (k == j) {
2788  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2789  }
2790  }
2791  if (charge(k) < 0.0) {
2792  size_t izeta = j;
2793  for (size_t m = 1; m < m_kk; m++) {
2794  if (charge(m) > 0.0) {
2795  size_t jzeta = m;
2796  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2797  double zeta_L = m_Psi_ijk_L[n];
2798  if (zeta_L != 0.0) {
2799  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
2800  }
2801  }
2802  }
2803  }
2804  }
2805  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j];
2806  }
2807  }
2808  double sum_m_phi_minus_1 = 2.0 *
2809  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2810  // Calculate the osmotic coefficient from
2811  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2812  double d_osmotic_coef_dT;
2813  if (molalitysum > 1.0E-150) {
2814  d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
2815  } else {
2816  d_osmotic_coef_dT = 0.0;
2817  }
2818 
2819  double d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
2820 
2821  // In Cantera, we define the activity coefficient of the solvent as
2822  //
2823  // act_0 = actcoeff_0 * Xmol_0
2824  //
2825  // We have just computed act_0. However, this routine returns
2826  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2827  m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
2828 }
2829 
2831 {
2832  static const int cacheId = m_cache.getId();
2833  CachedScalar cached = m_cache.getScalar(cacheId);
2834  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2835  return;
2836  }
2837 
2838  // Zero the unscaled 2nd derivatives
2840 
2841  //! Calculate the unscaled 2nd derivatives
2843 
2844  for (size_t k = 1; k < m_kk; k++) {
2845  if (CROP_speciesCropped_[k] == 2) {
2847  }
2848  }
2849 
2850  if (CROP_speciesCropped_[0]) {
2852  }
2853 
2854  // Scale the 2nd derivatives
2856 }
2857 
2859 {
2860  const double* molality = m_molalitiesCropped.data();
2861 
2862  // Local variables defined by Coltrin
2863  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2864 
2865  // Molality based ionic strength of the solution
2866  double Is = 0.0;
2867 
2868  // Molar charge of the solution: In Pitzer's notation, this is his variable
2869  // called "Z".
2870  double molarcharge = 0.0;
2871 
2872  // molalitysum is the sum of the molalities over all solutes, even those
2873  // with zero charge.
2874  double molalitysum = 0.0;
2875 
2876  // Make sure the counter variables are setup
2877  counterIJ_setup();
2878 
2879  // ---------- Calculate common sums over solutes ---------------------
2880  for (size_t n = 1; n < m_kk; n++) {
2881  // ionic strength
2882  Is += charge(n) * charge(n) * molality[n];
2883  // total molar charge
2884  molarcharge += fabs(charge(n)) * molality[n];
2885  molalitysum += molality[n];
2886  }
2887  Is *= 0.5;
2888 
2889  // Store the ionic molality in the object for reference.
2890  m_IionicMolality = Is;
2891  sqrtIs = sqrt(Is);
2892 
2893  // The following call to calc_lambdas() calculates all 16 elements of the
2894  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
2895  calc_lambdas(Is);
2896 
2897  // Step 2: Find the coefficients E-theta and E-thetaprime for all
2898  // combinations of positive unlike charges up to 4
2899  for (int z1 = 1; z1 <=4; z1++) {
2900  for (int z2 =1; z2 <=4; z2++) {
2901  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2902  }
2903  }
2904 
2905  // calculate gfunc(x) and hfunc(x) for each cation-anion pair MX. In the
2906  // original literature, hfunc, was called gprime. However, it's not the
2907  // derivative of gfunc(x), so I renamed it.
2908  for (size_t i = 1; i < (m_kk - 1); i++) {
2909  for (size_t j = (i+1); j < m_kk; j++) {
2910  // Find the counterIJ for the symmetric binary interaction
2911  size_t n = m_kk*i + j;
2912  size_t counterIJ = m_CounterIJ[n];
2913 
2914  // Only loop over oppositely charge species
2915  if (charge(i)*charge(j) < 0) {
2916  // x is a reduced function variable
2917  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2918  if (x1 > 1.0E-100) {
2919  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
2920  m_hfunc_IJ[counterIJ] = -2.0*
2921  (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
2922  } else {
2923  m_gfunc_IJ[counterIJ] = 0.0;
2924  m_hfunc_IJ[counterIJ] = 0.0;
2925  }
2926 
2927  if (m_Beta2MX_ij_LL[counterIJ] != 0.0) {
2928  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
2929  if (x2 > 1.0E-100) {
2930  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2931  m_h2func_IJ[counterIJ] = -2.0 *
2932  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2933  } else {
2934  m_g2func_IJ[counterIJ] = 0.0;
2935  m_h2func_IJ[counterIJ] = 0.0;
2936  }
2937  }
2938  } else {
2939  m_gfunc_IJ[counterIJ] = 0.0;
2940  m_hfunc_IJ[counterIJ] = 0.0;
2941  }
2942  }
2943  }
2944 
2945  // SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L
2946  // These are now temperature derivatives of the previously calculated
2947  // quantities.
2948  for (size_t i = 1; i < m_kk - 1; i++) {
2949  for (size_t j = i+1; j < m_kk; j++) {
2950  // Find the counterIJ for the symmetric binary interaction
2951  size_t n = m_kk*i + j;
2952  size_t counterIJ = m_CounterIJ[n];
2953 
2954  // both species have a non-zero charge, and one is positive
2955  // and the other is negative
2956  if (charge(i)*charge(j) < 0.0) {
2957  m_BMX_IJ_LL[counterIJ] = m_Beta0MX_ij_LL[counterIJ]
2958  + m_Beta1MX_ij_LL[counterIJ] * m_gfunc_IJ[counterIJ]
2959  + m_Beta2MX_ij_LL[counterIJ] * m_g2func_IJ[counterIJ];
2960  if (Is > 1.0E-150) {
2961  m_BprimeMX_IJ_LL[counterIJ] = (m_Beta1MX_ij_LL[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
2962  m_Beta2MX_ij_LL[counterIJ] * m_h2func_IJ[counterIJ]/Is);
2963  } else {
2964  m_BprimeMX_IJ_LL[counterIJ] = 0.0;
2965  }
2966  m_BphiMX_IJ_LL[counterIJ] = m_BMX_IJ_LL[counterIJ] + Is*m_BprimeMX_IJ_LL[counterIJ];
2967  } else {
2968  m_BMX_IJ_LL[counterIJ] = 0.0;
2969  m_BprimeMX_IJ_LL[counterIJ] = 0.0;
2970  m_BphiMX_IJ_LL[counterIJ] = 0.0;
2971  }
2972  }
2973  }
2974 
2975  // --------- SUBSECTION TO CALCULATE CMX_LL ----------
2976  for (size_t i = 1; i < m_kk-1; i++) {
2977  for (size_t j = i+1; j < m_kk; j++) {
2978  // Find the counterIJ for the symmetric binary interaction
2979  size_t n = m_kk*i + j;
2980  size_t counterIJ = m_CounterIJ[n];
2981 
2982  // both species have a non-zero charge, and one is positive
2983  // and the other is negative
2984  if (charge(i)*charge(j) < 0.0) {
2985  m_CMX_IJ_LL[counterIJ] = m_CphiMX_ij_LL[counterIJ]/
2986  (2.0* sqrt(fabs(charge(i)*charge(j))));
2987  } else {
2988  m_CMX_IJ_LL[counterIJ] = 0.0;
2989  }
2990  }
2991  }
2992 
2993  // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2994  for (size_t i = 1; i < m_kk-1; i++) {
2995  for (size_t j = i+1; j < m_kk; j++) {
2996  // Find the counterIJ for the symmetric binary interaction
2997  size_t n = m_kk*i + j;
2998  size_t counterIJ = m_CounterIJ[n];
2999 
3000  // both species have a non-zero charge, and one is positive
3001  // and the other is negative
3002  if (charge(i)*charge(j) > 0) {
3003  m_Phi_IJ_LL[counterIJ] = m_Theta_ij_LL[counterIJ];
3004  m_Phiprime_IJ[counterIJ] = 0.0;
3005  m_PhiPhi_IJ_LL[counterIJ] = m_Phi_IJ_LL[counterIJ];
3006  } else {
3007  m_Phi_IJ_LL[counterIJ] = 0.0;
3008  m_Phiprime_IJ[counterIJ] = 0.0;
3009  m_PhiPhi_IJ_LL[counterIJ] = 0.0;
3010  }
3011  }
3012  }
3013 
3014  // ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
3015  double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
3016  double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3017  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3018  for (size_t i = 1; i < m_kk-1; i++) {
3019  for (size_t j = i+1; j < m_kk; j++) {
3020  // Find the counterIJ for the symmetric binary interaction
3021  size_t n = m_kk*i + j;
3022  size_t counterIJ = m_CounterIJ[n];
3023 
3024  // both species have a non-zero charge, and one is positive
3025  // and the other is negative
3026  if (charge(i)*charge(j) < 0) {
3027  d2FdT2 += molality[i]*molality[j] * m_BprimeMX_IJ_LL[counterIJ];
3028  }
3029 
3030  // Both species have a non-zero charge, and they
3031  // have the same sign, that is, both positive or both negative.
3032  if (charge(i)*charge(j) > 0) {
3033  d2FdT2 += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
3034  }
3035  }
3036  }
3037 
3038  for (size_t i = 1; i < m_kk; i++) {
3039  // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
3040  if (charge(i) > 0) {
3041  // species i is the cation (positive) to calc the actcoeff
3042  double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3043  double sum1 = 0.0;
3044  double sum2 = 0.0;
3045  double sum3 = 0.0;
3046  double sum4 = 0.0;
3047  double sum5 = 0.0;
3048  for (size_t j = 1; j < m_kk; j++) {
3049  // Find the counterIJ for the symmetric binary interaction
3050  size_t n = m_kk*i + j;
3051  size_t counterIJ = m_CounterIJ[n];
3052 
3053  if (charge(j) < 0.0) {
3054  // sum over all anions
3055  sum1 += molality[j]*
3056  (2.0*m_BMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3057  if (j < m_kk-1) {
3058  // This term is the ternary interaction involving the
3059  // non-duplicate sum over double anions, j, k, with
3060  // respect to the cation, i.
3061  for (size_t k = j+1; k < m_kk; k++) {
3062  // an inner sum over all anions
3063  if (charge(k) < 0.0) {
3064  n = k + j * m_kk + i * m_kk * m_kk;
3065  sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3066  }
3067  }
3068  }
3069  }
3070 
3071  if (charge(j) > 0.0) {
3072  // sum over all cations
3073  if (j != i) {
3074  sum2 += molality[j]*(2.0*m_Phi_IJ_LL[counterIJ]);
3075  }
3076  for (size_t k = 1; k < m_kk; k++) {
3077  if (charge(k) < 0.0) {
3078  // two inner sums over anions
3079  n = k + j * m_kk + i * m_kk * m_kk;
3080  sum2 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3081 
3082  // Find the counterIJ for the j,k interaction
3083  n = m_kk*j + k;
3084  size_t counterIJ2 = m_CounterIJ[n];
3085  sum4 += fabs(charge(i)) *
3086  molality[j]*molality[k]*m_CMX_IJ_LL[counterIJ2];
3087  }
3088  }
3089  }
3090 
3091  // Handle neutral j species
3092  if (charge(j) == 0) {
3093  sum5 += molality[j]*2.0*m_Lambda_nj_LL(j,i);
3094  // Zeta interaction term
3095  for (size_t k = 1; k < m_kk; k++) {
3096  if (charge(k) < 0.0) {
3097  size_t izeta = j;
3098  size_t jzeta = i;
3099  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3100  double zeta_LL = m_Psi_ijk_LL[n];
3101  if (zeta_LL != 0.0) {
3102  sum5 += molality[j]*molality[k]*zeta_LL;
3103  }
3104  }
3105  }
3106  }
3107  }
3108  // Add all of the contributions up to yield the log of the
3109  // solute activity coefficients (molality scale)
3111  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3112  }
3113 
3114  // ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------
3115  if (charge(i) < 0) {
3116  // species i is an anion (negative)
3117  double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3118  double sum1 = 0.0;
3119  double sum2 = 0.0;
3120  double sum3 = 0.0;
3121  double sum4 = 0.0;
3122  double sum5 = 0.0;
3123  for (size_t j = 1; j < m_kk; j++) {
3124  // Find the counterIJ for the symmetric binary interaction
3125  size_t n = m_kk*i + j;
3126  size_t counterIJ = m_CounterIJ[n];
3127 
3128  // For Anions, do the cation interactions.
3129  if (charge(j) > 0) {
3130  sum1 += molality[j]*
3131  (2.0*m_BMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3132  if (j < m_kk-1) {
3133  for (size_t k = j+1; k < m_kk; k++) {
3134  // an inner sum over all cations
3135  if (charge(k) > 0) {
3136  n = k + j * m_kk + i * m_kk * m_kk;
3137  sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3138  }
3139  }
3140  }
3141  }
3142 
3143  // For Anions, do the other anion interactions.
3144  if (charge(j) < 0.0) {
3145  // sum over all anions
3146  if (j != i) {
3147  sum2 += molality[j]*(2.0*m_Phi_IJ_LL[counterIJ]);
3148  }
3149  for (size_t k = 1; k < m_kk; k++) {
3150  if (charge(k) > 0.0) {
3151  // two inner sums over cations
3152  n = k + j * m_kk + i * m_kk * m_kk;
3153  sum2 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3154  // Find the counterIJ for the symmetric binary interaction
3155  n = m_kk*j + k;
3156  size_t counterIJ2 = m_CounterIJ[n];
3157  sum4 += fabs(charge(i)) *
3158  molality[j]*molality[k]*m_CMX_IJ_LL[counterIJ2];
3159  }
3160  }
3161  }
3162 
3163  // for Anions, do the neutral species interaction
3164  if (charge(j) == 0.0) {
3165  sum5 += molality[j]*2.0*m_Lambda_nj_LL(j,i);
3166  // Zeta interaction term
3167  for (size_t k = 1; k < m_kk; k++) {
3168  if (charge(k) > 0.0) {
3169  size_t izeta = j;
3170  size_t jzeta = k;
3171  size_t kzeta = i;
3172  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3173  double zeta_LL = m_Psi_ijk_LL[n];
3174  if (zeta_LL != 0.0) {
3175  sum5 += molality[j]*molality[k]*zeta_LL;
3176  }
3177  }
3178  }
3179  }
3180  }
3182  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3183  }
3184 
3185  // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
3186  // equations agree with my notes,
3187  // Equations agree with Pitzer,
3188  if (charge(i) == 0.0) {
3189  double sum1 = 0.0;
3190  double sum3 = 0.0;
3191  for (size_t j = 1; j < m_kk; j++) {
3192  sum1 += molality[j]*2.0*m_Lambda_nj_LL(i,j);
3193  // Zeta term -> we piggyback on the psi term
3194  if (charge(j) > 0.0) {
3195  for (size_t k = 1; k < m_kk; k++) {
3196  if (charge(k) < 0.0) {
3197  size_t n = k + j * m_kk + i * m_kk * m_kk;
3198  sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3199  }
3200  }
3201  }
3202  }
3203  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
3204  m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
3205  }
3206  }
3207 
3208  // ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
3209  double sum1 = 0.0;
3210  double sum2 = 0.0;
3211  double sum3 = 0.0;
3212  double sum4 = 0.0;
3213  double sum5 = 0.0;
3214  double sum6 = 0.0;
3215  double sum7 = 0.0;
3216 
3217  // term1 is the temperature derivative of the DH term in the osmotic
3218  // coefficient expression
3219  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
3220  // Is = Ionic strength on the molality scale (units of (gmol/kg))
3221  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3222  double term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3223 
3224  for (size_t j = 1; j < m_kk; j++) {
3225  // Loop Over Cations
3226  if (charge(j) > 0.0) {
3227  for (size_t k = 1; k < m_kk; k++) {
3228  if (charge(k) < 0.0) {
3229  // Find the counterIJ for the symmetric j,k binary interaction
3230  size_t n = m_kk*j + k;
3231  size_t counterIJ = m_CounterIJ[n];
3232 
3233  sum1 += molality[j]*molality[k] *
3234  (m_BphiMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3235  }
3236  }
3237 
3238  for (size_t k = j+1; k < m_kk; k++) {
3239  if (j == (m_kk-1)) {
3240  // we should never reach this step
3241  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3242  "logic error 1 in Step 9 of hmw_act");
3243  }
3244  if (charge(k) > 0.0) {
3245  // Find the counterIJ for the symmetric j,k binary interaction
3246  // between 2 cations.
3247  size_t n = m_kk*j + k;
3248  size_t counterIJ = m_CounterIJ[n];
3249  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_LL[counterIJ];
3250  for (size_t m = 1; m < m_kk; m++) {
3251  if (charge(m) < 0.0) {
3252  // species m is an anion
3253  n = m + k * m_kk + j * m_kk * m_kk;
3254  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_LL[n];
3255  }
3256  }
3257  }
3258  }
3259  }
3260 
3261  // Loop Over Anions
3262  if (charge(j) < 0) {
3263  for (size_t k = j+1; k < m_kk; k++) {
3264  if (j == m_kk-1) {
3265  // we should never reach this step
3266  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3267  "logic error 2 in Step 9 of hmw_act");
3268  }
3269  if (charge(k) < 0) {
3270  // Find the counterIJ for the symmetric j,k binary interaction
3271  // between two anions
3272  size_t n = m_kk*j + k;
3273  size_t counterIJ = m_CounterIJ[n];
3274 
3275  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_LL[counterIJ];
3276  for (size_t m = 1; m < m_kk; m++) {
3277  if (charge(m) > 0.0) {
3278  n = m + k * m_kk + j * m_kk * m_kk;
3279  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_LL[n];
3280  }
3281  }
3282  }
3283  }
3284  }
3285 
3286  // Loop Over Neutral Species
3287  if (charge(j) == 0) {
3288  for (size_t k = 1; k < m_kk; k++) {
3289  if (charge(k) < 0.0) {
3290  sum4 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3291  }
3292  if (charge(k) > 0.0) {
3293  sum5 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3294  }
3295  if (charge(k) == 0.0) {
3296  if (k > j) {
3297  sum6 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3298  } else if (k == j) {
3299  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3300  }
3301  }
3302  if (charge(k) < 0.0) {
3303  size_t izeta = j;
3304  for (size_t m = 1; m < m_kk; m++) {
3305  if (charge(m) > 0.0) {
3306  size_t jzeta = m;
3307  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3308  double zeta_LL = m_Psi_ijk_LL[n];
3309  if (zeta_LL != 0.0) {
3310  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
3311  }
3312  }
3313  }
3314  }
3315  }
3316 
3317  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j];
3318  }
3319  }
3320  double sum_m_phi_minus_1 = 2.0 *
3321  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3322  // Calculate the osmotic coefficient from
3323  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3324  double d2_osmotic_coef_dT2;
3325  if (molalitysum > 1.0E-150) {
3326  d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3327  } else {
3328  d2_osmotic_coef_dT2 = 0.0;
3329  }
3330  double d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
3331 
3332  // In Cantera, we define the activity coefficient of the solvent as
3333  //
3334  // act_0 = actcoeff_0 * Xmol_0
3335  //
3336  // We have just computed act_0. However, this routine returns
3337  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3338  m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
3339 }
3340 
3342 {
3343  static const int cacheId = m_cache.getId();
3344  CachedScalar cached = m_cache.getScalar(cacheId);
3345  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
3346  return;
3347  }
3348 
3349  m_dlnActCoeffMolaldP_Unscaled.assign(m_kk, 0.0);
3351 
3352  for (size_t k = 1; k < m_kk; k++) {
3353  if (CROP_speciesCropped_[k] == 2) {
3355  }
3356  }
3357 
3358  if (CROP_speciesCropped_[0]) {
3360  }
3361 
3363 }
3364 
3366 {
3367  const double* molality = m_molalitiesCropped.data();
3368 
3369  // Local variables defined by Coltrin
3370  double etheta[5][5], etheta_prime[5][5], sqrtIs;
3371 
3372  // Molality based ionic strength of the solution
3373  double Is = 0.0;
3374 
3375  // Molar charge of the solution: In Pitzer's notation, this is his variable
3376  // called "Z".
3377  double molarcharge = 0.0;
3378 
3379  // molalitysum is the sum of the molalities over all solutes, even those
3380  // with zero charge.
3381  double molalitysum = 0.0;
3382  double currTemp = temperature();
3383  double currPres = pressure();
3384 
3385  // Make sure the counter variables are setup
3386  counterIJ_setup();
3387 
3388  // ---------- Calculate common sums over solutes ---------------------
3389  for (size_t n = 1; n < m_kk; n++) {
3390  // ionic strength
3391  Is += charge(n) * charge(n) * molality[n];
3392  // total molar charge
3393  molarcharge += fabs(charge(n)) * molality[n];
3394  molalitysum += molality[n];
3395  }
3396  Is *= 0.5;
3397 
3398  // Store the ionic molality in the object for reference.
3399  m_IionicMolality = Is;
3400  sqrtIs = sqrt(Is);
3401 
3402  // The following call to calc_lambdas() calculates all 16 elements of the
3403  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
3404  calc_lambdas(Is);
3405 
3406  // Step 2: Find the coefficients E-theta and E-thetaprime for all
3407  // combinations of positive unlike charges up to 4
3408  for (int z1 = 1; z1 <=4; z1++) {
3409  for (int z2 =1; z2 <=4; z2++) {
3410  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3411  }
3412  }
3413 
3414  // calculate g(x) and hfunc(x) for each cation-anion pair MX
3415  // In the original literature, hfunc, was called gprime. However,
3416  // it's not the derivative of g(x), so I renamed it.
3417  for (size_t i = 1; i < (m_kk - 1); i++) {
3418  for (size_t j = (i+1); j < m_kk; j++) {
3419  // Find the counterIJ for the symmetric binary interaction
3420  size_t n = m_kk*i + j;
3421  size_t counterIJ = m_CounterIJ[n];
3422 
3423  // Only loop over oppositely charge species
3424  if (charge(i)*charge(j) < 0) {
3425  // x is a reduced function variable
3426  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
3427  if (x1 > 1.0E-100) {
3428  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3429  m_hfunc_IJ[counterIJ] = -2.0*
3430  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
3431  } else {
3432  m_gfunc_IJ[counterIJ] = 0.0;
3433  m_hfunc_IJ[counterIJ] = 0.0;
3434  }
3435 
3436  if (m_Beta2MX_ij_P[counterIJ] != 0.0) {
3437  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
3438  if (x2 > 1.0E-100) {
3439  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3440  m_h2func_IJ[counterIJ] = -2.0 *
3441  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3442  } else {
3443  m_g2func_IJ[counterIJ] = 0.0;
3444  m_h2func_IJ[counterIJ] = 0.0;
3445  }
3446  }
3447  } else {
3448  m_gfunc_IJ[counterIJ] = 0.0;
3449  m_hfunc_IJ[counterIJ] = 0.0;
3450  }
3451  }
3452  }
3453 
3454  // SUBSECTION TO CALCULATE BMX_P, BprimeMX_P, BphiMX_P
3455  // These are now temperature derivatives of the previously calculated
3456  // quantities.
3457  for (size_t i = 1; i < m_kk - 1; i++) {
3458  for (size_t j = i+1; j < m_kk; j++) {
3459  // Find the counterIJ for the symmetric binary interaction
3460  size_t n = m_kk*i + j;
3461  size_t counterIJ = m_CounterIJ[n];
3462 
3463  // both species have a non-zero charge, and one is positive
3464  // and the other is negative
3465  if (charge(i)*charge(j) < 0.0) {
3466  m_BMX_IJ_P[counterIJ] = m_Beta0MX_ij_P[counterIJ]
3467  + m_Beta1MX_ij_P[counterIJ] * m_gfunc_IJ[counterIJ]
3468  + m_Beta2MX_ij_P[counterIJ] * m_g2func_IJ[counterIJ];
3469  if (Is > 1.0E-150) {
3470  m_BprimeMX_IJ_P[counterIJ] = (m_Beta1MX_ij_P[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
3471  m_Beta2MX_ij_P[counterIJ] * m_h2func_IJ[counterIJ]/Is);
3472  } else {
3473  m_BprimeMX_IJ_P[counterIJ] = 0.0;
3474  }
3475  m_BphiMX_IJ_P[counterIJ] = m_BMX_IJ_P[counterIJ] + Is*m_BprimeMX_IJ_P[counterIJ];
3476  } else {
3477  m_BMX_IJ_P[counterIJ] = 0.0;
3478  m_BprimeMX_IJ_P[counterIJ] = 0.0;
3479  m_BphiMX_IJ_P[counterIJ] = 0.0;
3480  }
3481  }
3482  }
3483 
3484  // --------- SUBSECTION TO CALCULATE CMX_P ----------
3485  for (size_t i = 1; i < m_kk-1; i++) {
3486  for (size_t j = i+1; j < m_kk; j++) {
3487  // Find the counterIJ for the symmetric binary interaction
3488  size_t n = m_kk*i + j;
3489  size_t counterIJ = m_CounterIJ[n];
3490 
3491  // both species have a non-zero charge, and one is positive
3492  // and the other is negative
3493  if (charge(i)*charge(j) < 0.0) {
3494  m_CMX_IJ_P[counterIJ] = m_CphiMX_ij_P[counterIJ]/
3495  (2.0* sqrt(fabs(charge(i)*charge(j))));
3496  } else {
3497  m_CMX_IJ_P[counterIJ] = 0.0;
3498  }
3499  }
3500  }
3501 
3502  // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
3503  for (size_t i = 1; i < m_kk-1; i++) {
3504  for (size_t j = i+1; j < m_kk; j++) {
3505  // Find the counterIJ for the symmetric binary interaction
3506  size_t n = m_kk*i + j;
3507  size_t counterIJ = m_CounterIJ[n];
3508 
3509  // both species have a non-zero charge, and one is positive
3510  // and the other is negative
3511  if (charge(i)*charge(j) > 0) {
3512  m_Phi_IJ_P[counterIJ] = m_Theta_ij_P[counterIJ];
3513  m_Phiprime_IJ[counterIJ] = 0.0;
3514  m_PhiPhi_IJ_P[counterIJ] = m_Phi_IJ_P[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
3515  } else {
3516  m_Phi_IJ_P[counterIJ] = 0.0;
3517  m_Phiprime_IJ[counterIJ] = 0.0;
3518  m_PhiPhi_IJ_P[counterIJ] = 0.0;
3519  }
3520  }
3521  }
3522 
3523  // ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
3524  double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
3525  double dAphidP = dA_DebyedP /3.0;
3526  double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3527  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3528  for (size_t i = 1; i < m_kk-1; i++) {
3529  for (size_t j = i+1; j < m_kk; j++) {
3530  // Find the counterIJ for the symmetric binary interaction
3531  size_t n = m_kk*i + j;
3532  size_t counterIJ = m_CounterIJ[n];
3533 
3534  // both species have a non-zero charge, and one is positive
3535  // and the other is negative
3536  if (charge(i)*charge(j) < 0) {
3537  dFdP += molality[i]*molality[j] * m_BprimeMX_IJ_P[counterIJ];
3538  }
3539 
3540  // Both species have a non-zero charge, and they
3541  // have the same sign, that is, both positive or both negative.
3542  if (charge(i)*charge(j) > 0) {
3543  dFdP += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
3544  }
3545  }
3546  }
3547 
3548  for (size_t i = 1; i < m_kk; i++) {
3549  // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
3550  if (charge(i) > 0) {
3551  // species i is the cation (positive) to calc the actcoeff
3552  double zsqdFdP = charge(i)*charge(i)*dFdP;
3553  double sum1 = 0.0;
3554  double sum2 = 0.0;
3555  double sum3 = 0.0;
3556  double sum4 = 0.0;
3557  double sum5 = 0.0;
3558  for (size_t j = 1; j < m_kk; j++) {
3559  // Find the counterIJ for the symmetric binary interaction
3560  size_t n = m_kk*i + j;
3561  size_t counterIJ = m_CounterIJ[n];
3562 
3563  if (charge(j) < 0.0) {
3564  // sum over all anions
3565  sum1 += molality[j]*
3566  (2.0*m_BMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3567  if (j < m_kk-1) {
3568  // This term is the ternary interaction involving the
3569  // non-duplicate sum over double anions, j, k, with
3570  // respect to the cation, i.
3571  for (size_t k = j+1; k < m_kk; k++) {
3572  // an inner sum over all anions
3573  if (charge(k) < 0.0) {
3574  n = k + j * m_kk + i * m_kk * m_kk;
3575  sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3576  }
3577  }
3578  }
3579  }
3580 
3581  if (charge(j) > 0.0) {
3582  // sum over all cations
3583  if (j != i) {
3584  sum2 += molality[j]*(2.0*m_Phi_IJ_P[counterIJ]);
3585  }
3586  for (size_t k = 1; k < m_kk; k++) {
3587  if (charge(k) < 0.0) {
3588  // two inner sums over anions
3589  n = k + j * m_kk + i * m_kk * m_kk;
3590  sum2 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3591 
3592  // Find the counterIJ for the j,k interaction
3593  n = m_kk*j + k;
3594  size_t counterIJ2 = m_CounterIJ[n];
3595  sum4 += fabs(charge(i)) *
3596  molality[j]*molality[k]*m_CMX_IJ_P[counterIJ2];
3597  }
3598  }
3599  }
3600 
3601  // for Anions, do the neutral species interaction
3602  if (charge(j) == 0) {
3603  sum5 += molality[j]*2.0*m_Lambda_nj_P(j,i);
3604  // Zeta interaction term
3605  for (size_t k = 1; k < m_kk; k++) {
3606  if (charge(k) < 0.0) {
3607  size_t izeta = j;
3608  size_t jzeta = i;
3609  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3610  double zeta_P = m_Psi_ijk_P[n];
3611  if (zeta_P != 0.0) {
3612  sum5 += molality[j]*molality[k]*zeta_P;
3613  }
3614  }
3615  }
3616  }
3617  }
3618 
3619  // Add all of the contributions up to yield the log of the
3620  // solute activity coefficients (molality scale)
3622  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3623  }
3624 
3625  // ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR ANIONS ------
3626  if (charge(i) < 0) {
3627  // species i is an anion (negative)
3628  double zsqdFdP = charge(i)*charge(i)*dFdP;
3629  double sum1 = 0.0;
3630  double sum2 = 0.0;
3631  double sum3 = 0.0;
3632  double sum4 = 0.0;
3633  double sum5 = 0.0;
3634  for (size_t j = 1; j < m_kk; j++) {
3635  // Find the counterIJ for the symmetric binary interaction
3636  size_t n = m_kk*i + j;
3637  size_t counterIJ = m_CounterIJ[n];
3638 
3639  // For Anions, do the cation interactions.
3640  if (charge(j) > 0) {
3641  sum1 += molality[j] *
3642  (2.0*m_BMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3643  if (j < m_kk-1) {
3644  for (size_t k = j+1; k < m_kk; k++) {
3645  // an inner sum over all cations
3646  if (charge(k) > 0) {
3647  n = k + j * m_kk + i * m_kk * m_kk;
3648  sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3649  }
3650  }
3651  }
3652  }
3653 
3654  // For Anions, do the other anion interactions.
3655  if (charge(j) < 0.0) {
3656  // sum over all anions
3657  if (j != i) {
3658  sum2 += molality[j]*(2.0*m_Phi_IJ_P[counterIJ]);
3659  }
3660  for (size_t k = 1; k < m_kk; k++) {
3661  if (charge(k) > 0.0) {
3662  // two inner sums over cations
3663  n = k + j * m_kk + i * m_kk * m_kk;
3664  sum2 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3665  // Find the counterIJ for the symmetric binary interaction
3666  n = m_kk*j + k;
3667  size_t counterIJ2 = m_CounterIJ[n];
3668  sum4 += fabs(charge(i))*
3669  molality[j]*molality[k]*m_CMX_IJ_P[counterIJ2];
3670  }
3671  }
3672  }
3673 
3674  // for Anions, do the neutral species interaction
3675  if (charge(j) == 0.0) {
3676  sum5 += molality[j]*2.0*m_Lambda_nj_P(j,i);
3677  // Zeta interaction term
3678  for (size_t k = 1; k < m_kk; k++) {
3679  if (charge(k) > 0.0) {
3680  size_t izeta = j;
3681  size_t jzeta = k;
3682  size_t kzeta = i;
3683  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3684  double zeta_P = m_Psi_ijk_P[n];
3685  if (zeta_P != 0.0) {
3686  sum5 += molality[j]*molality[k]*zeta_P;
3687  }
3688  }
3689  }
3690  }
3691  }
3693  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3694  }
3695 
3696  // ------ SUBSECTION FOR CALCULATING d NEUTRAL SOLUTE ACT COEFF dP -----
3697  if (charge(i) == 0.0) {
3698  double sum1 = 0.0;
3699  double sum3 = 0.0;
3700  for (size_t j = 1; j < m_kk; j++) {
3701  sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j);
3702  // Zeta term -> we piggyback on the psi term
3703  if (charge(j) > 0.0) {
3704  for (size_t k = 1; k < m_kk; k++) {
3705  if (charge(k) < 0.0) {
3706  size_t n = k + j * m_kk + i * m_kk * m_kk;
3707  sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3708  }
3709  }
3710  }
3711  }
3712  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
3713  m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
3714  }
3715  }
3716 
3717  // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
3718  double sum1 = 0.0;
3719  double sum2 = 0.0;
3720  double sum3 = 0.0;
3721  double sum4 = 0.0;
3722  double sum5 = 0.0;
3723  double sum6 = 0.0;
3724  double sum7 = 0.0;
3725 
3726  // term1 is the temperature derivative of the DH term in the osmotic
3727  // coefficient expression
3728  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
3729  // Is = Ionic strength on the molality scale (units of (gmol/kg))
3730  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3731  double term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3732 
3733  for (size_t j = 1; j < m_kk; j++) {
3734  // Loop Over Cations
3735  if (charge(j) > 0.0) {
3736  for (size_t k = 1; k < m_kk; k++) {
3737  if (charge(k) < 0.0) {
3738  // Find the counterIJ for the symmetric j,k binary interaction
3739  size_t n = m_kk*j + k;
3740  size_t counterIJ = m_CounterIJ[n];
3741  sum1 += molality[j]*molality[k]*
3742  (m_BphiMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3743  }
3744  }
3745 
3746  for (size_t k = j+1; k < m_kk; k++) {
3747  if (j == (m_kk-1)) {
3748  // we should never reach this step
3749  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3750  "logic error 1 in Step 9 of hmw_act");
3751  }
3752  if (charge(k) > 0.0) {
3753  // Find the counterIJ for the symmetric j,k binary interaction
3754  // between 2 cations.
3755  size_t n = m_kk*j + k;
3756  size_t counterIJ = m_CounterIJ[n];
3757  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_P[counterIJ];
3758  for (size_t m = 1; m < m_kk; m++) {
3759  if (charge(m) < 0.0) {
3760  // species m is an anion
3761  n = m + k * m_kk + j * m_kk * m_kk;
3762  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_P[n];
3763  }
3764  }
3765  }
3766  }
3767  }
3768 
3769  // Loop Over Anions
3770  if (charge(j) < 0) {
3771  for (size_t k = j+1; k < m_kk; k++) {
3772  if (j == m_kk-1) {
3773  // we should never reach this step
3774  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3775  "logic error 2 in Step 9 of hmw_act");
3776  }
3777  if (charge(k) < 0) {
3778  // Find the counterIJ for the symmetric j,k binary interaction
3779  // between two anions
3780  size_t n = m_kk*j + k;
3781  size_t counterIJ = m_CounterIJ[n];
3782 
3783  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_P[counterIJ];
3784  for (size_t m = 1; m < m_kk; m++) {
3785  if (charge(m) > 0.0) {
3786  n = m + k * m_kk + j * m_kk * m_kk;
3787  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_P[n];
3788  }
3789  }
3790  }
3791  }
3792  }
3793 
3794  // Loop Over Neutral Species
3795  if (charge(j) == 0) {
3796  for (size_t k = 1; k < m_kk; k++) {
3797  if (charge(k) < 0.0) {
3798  sum4 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3799  }
3800  if (charge(k) > 0.0) {
3801  sum5 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3802  }
3803  if (charge(k) == 0.0) {
3804  if (k > j) {
3805  sum6 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3806  } else if (k == j) {
3807  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3808  }
3809  }
3810  if (charge(k) < 0.0) {
3811  size_t izeta = j;
3812  for (size_t m = 1; m < m_kk; m++) {
3813  if (charge(m) > 0.0) {
3814  size_t jzeta = m;
3815  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3816  double zeta_P = m_Psi_ijk_P[n];
3817  if (zeta_P != 0.0) {
3818  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
3819  }
3820  }
3821  }
3822  }
3823  }
3824 
3825  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j];
3826  }
3827  }
3828  double sum_m_phi_minus_1 = 2.0 *
3829  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3830 
3831  // Calculate the osmotic coefficient from
3832  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3833  double d_osmotic_coef_dP;
3834  if (molalitysum > 1.0E-150) {
3835  d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3836  } else {
3837  d_osmotic_coef_dP = 0.0;
3838  }
3839  double d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
3840 
3841  // In Cantera, we define the activity coefficient of the solvent as
3842  //
3843  // act_0 = actcoeff_0 * Xmol_0
3844  //
3845  // We have just computed act_0. However, this routine returns
3846  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3847  m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
3848 }
3849 
3850 void HMWSoln::calc_lambdas(double is) const
3851 {
3852  if( m_last_is == is ) {
3853  return;
3854  }
3855  m_last_is = is;
3856 
3857  // Coefficients c1-c4 are used to approximate the integral function "J";
3858  // aphi is the Debye-Huckel constant at 25 C
3859  double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
3860  double aphi = 0.392; /* Value at 25 C */
3861  if (is < 1.0E-150) {
3862  for (int i = 0; i < 17; i++) {
3863  elambda[i] = 0.0;
3864  elambda1[i] = 0.0;
3865  }
3866  return;
3867  }
3868 
3869  // Calculate E-lambda terms for charge combinations of like sign,
3870  // using method of Pitzer (1975). Charges up to 4 are calculated.
3871  for (int i=1; i<=4; i++) {
3872  for (int j=i; j<=4; j++) {
3873  int ij = i*j;
3874 
3875  // calculate the product of the charges
3876  double zprod = (double)ij;
3877 
3878  // calculate Xmn (A1) from Harvie, Weare (1980).
3879  double x = 6.0* zprod * aphi * sqrt(is); // eqn 23
3880 
3881  double jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); // eqn 47
3882 
3883  double t = c3 * c4 * pow(x,c4);
3884  double dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
3885  double jprime = (jfunc/x)*(1.0 + jfunc*dj);
3886 
3887  elambda[ij] = zprod*jfunc / (4.0*is); // eqn 14
3888  elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
3889  - elambda[ij])/is;
3890  }
3891  }
3892 }
3893 
3894 void HMWSoln::calc_thetas(int z1, int z2,
3895  double* etheta, double* etheta_prime) const
3896 {
3897  // Calculate E-theta(i) and E-theta'(I) using method of Pitzer (1987)
3898  int i = abs(z1);
3899  int j = abs(z2);
3900 
3901  AssertThrowMsg(i <= 4 && j <= 4, "HMWSoln::calc_thetas",
3902  "we shouldn't be here");
3903  AssertThrowMsg(i != 0 && j != 0, "HMWSoln::calc_thetas",
3904  "called with one species being neutral");
3905 
3906  // Check to see if the charges are of opposite sign. If they are of opposite
3907  // sign then their etheta interaction is zero.
3908  if (z1*z2 < 0) {
3909  *etheta = 0.0;
3910  *etheta_prime = 0.0;
3911  } else {
3912  // Actually calculate the interaction.
3913  double f1 = (double)i / (2.0 * j);
3914  double f2 = (double)j / (2.0 * i);
3915  *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i];
3916  *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
3917  }
3918 }
3919 
3921 {
3922  // Calculate the molalities. Currently, the molalities may not be current
3923  // with respect to the contents of the State objects' data.
3924  calcMolalities();
3925  double xmolSolvent = moleFraction(0);
3926  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
3927  // Exponentials - trial 2
3928  if (xmolSolvent > IMS_X_o_cutoff_) {
3929  for (size_t k = 1; k < m_kk; k++) {
3930  IMS_lnActCoeffMolal_[k]= 0.0;
3931  }
3932  IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
3933  return;
3934  } else {
3935  double xoverc = xmolSolvent/IMS_cCut_;
3936  double eterm = std::exp(-xoverc);
3937 
3938  double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
3939  + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
3940  double f_prime = 1.0 + eterm*fptmp;
3941  double f = xmolSolvent + IMS_efCut_
3942  + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
3943 
3944  double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
3945  + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
3946  double g_prime = 1.0 + eterm*gptmp;
3947  double g = xmolSolvent + IMS_egCut_
3948  + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
3949 
3950  double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
3951  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
3952  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
3953 
3954  tmp = log(xx) + lngammak;
3955  for (size_t k = 1; k < m_kk; k++) {
3956  IMS_lnActCoeffMolal_[k]= tmp;
3957  }
3958  IMS_lnActCoeffMolal_[0] = lngammao;
3959  }
3960  return;
3961 }
3962 
3964 {
3965  calcMolalities();
3966  vector<double>& moleF = m_tmpV;
3967 
3968  // Update the coefficients wrt Temperature. Calculate the derivatives as well
3970  getMoleFractions(moleF.data());
3971 
3972  writelog("Index Name MoleF MolalityCropped Charge\n");
3973  for (size_t k = 0; k < m_kk; k++) {
3974  writelogf("%2d %-16s %14.7le %14.7le %5.1f \n",
3975  k, speciesName(k), moleF[k], m_molalitiesCropped[k], charge(k));
3976  }
3977 
3978  writelog("\n Species Species beta0MX "
3979  "beta1MX beta2MX CphiMX alphaMX thetaij\n");
3980  for (size_t i = 1; i < m_kk - 1; i++) {
3981  for (size_t j = i+1; j < m_kk; j++) {
3982  size_t n = i * m_kk + j;
3983  size_t ct = m_CounterIJ[n];
3984  writelogf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
3985  speciesName(i), speciesName(j),
3986  m_Beta0MX_ij[ct], m_Beta1MX_ij[ct],
3987  m_Beta2MX_ij[ct], m_CphiMX_ij[ct],
3988  m_Alpha1MX_ij[ct], m_Theta_ij[ct]);
3989  }
3990  }
3991 
3992  writelog("\n Species Species Species psi \n");
3993  for (size_t i = 1; i < m_kk; i++) {
3994  for (size_t j = 1; j < m_kk; j++) {
3995  for (size_t k = 1; k < m_kk; k++) {
3996  size_t n = k + j * m_kk + i * m_kk * m_kk;
3997  if (m_Psi_ijk[n] != 0.0) {
3998  writelogf(" %-16s %-16s %-16s %9.5f \n",
3999  speciesName(i), speciesName(j),
4000  speciesName(k), m_Psi_ijk[n]);
4001  }
4002  }
4003  }
4004  }
4005 }
4006 
4007 void HMWSoln::applyphScale(double* acMolality) const
4008 {
4010  return;
4011  }
4013  double lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
4014  double lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
4015  double afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4016  for (size_t k = 0; k < m_kk; k++) {
4017  acMolality[k] *= exp(charge(k) * afac);
4018  }
4019 }
4020 
4022 {
4025  return;
4026  }
4028  double lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
4029  double lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
4030  double afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4031  for (size_t k = 0; k < m_kk; k++) {
4033  }
4034 }
4035 
4037 {
4040  return;
4041  }
4043  double dlnGammaClM_dT_s2 = s_NBS_CLM_dlnMolalityActCoeff_dT();
4044  double dlnGammaCLM_dT_s1 = m_dlnActCoeffMolaldT_Unscaled[m_indexCLM];
4045  double afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
4046  for (size_t k = 0; k < m_kk; k++) {
4048  }
4049 }
4050 
4052 {
4055  return;
4056  }
4058  double d2lnGammaClM_dT2_s2 = s_NBS_CLM_d2lnMolalityActCoeff_dT2();
4059  double d2lnGammaCLM_dT2_s1 = m_d2lnActCoeffMolaldT2_Unscaled[m_indexCLM];
4060  double afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
4061  for (size_t k = 0; k < m_kk; k++) {
4063  }
4064 }
4065 
4067 {
4070  return;
4071  }
4073  double dlnGammaClM_dP_s2 = s_NBS_CLM_dlnMolalityActCoeff_dP();
4074  double dlnGammaCLM_dP_s1 = m_dlnActCoeffMolaldP_Unscaled[m_indexCLM];
4075  double afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
4076  for (size_t k = 0; k < m_kk; k++) {
4078  }
4079 }
4080 
4082 {
4083  double sqrtIs = sqrt(m_IionicMolality);
4084  double A = A_Debye_TP();
4085  double lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
4086  return lnGammaClMs2;
4087 }
4088 
4090 {
4091  double sqrtIs = sqrt(m_IionicMolality);
4092  double dAdT = dA_DebyedT_TP();
4093  return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
4094 }
4095 
4097 {
4098  double sqrtIs = sqrt(m_IionicMolality);
4099  double d2AdT2 = d2A_DebyedT2_TP();
4100  return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
4101 }
4102 
4104 {
4105  double sqrtIs = sqrt(m_IionicMolality);
4106  double dAdP = dA_DebyedP_TP();
4107  return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
4108 }
4109 
4110 }
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
Definition: AnyMap.h:427
size_t size() const
Returns the number of elements in this map.
Definition: AnyMap.h:622
bool hasKey(const string &key) const
Returns true if the map contains an item named key.
Definition: AnyMap.cpp:1423
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.cpp:47
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:203
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
void calcIMSCutoffParams_()
Precalculate the IMS Cutoff parameters for typeCutoff = 2.
Definition: HMWSoln.cpp:1473
vector< double > m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT.
Definition: HMWSoln.h:1721
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1575
vector< double > m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:1668
vector< double > m_Beta1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1459
double enthalpy_mole() const override
Molar enthalpy. Units: J/kmol.
Definition: HMWSoln.cpp:54
void applyphScale(double *acMolality) const override
Apply the current phScale to a set of activity Coefficients or activities.
Definition: HMWSoln.cpp:4007
unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition: HMWSoln.h:1423
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
Definition: HMWSoln.h:1452
double m_A_Debye
A_Debye: this expression appears on the top of the ln actCoeff term in the general Debye-Huckel expre...
Definition: HMWSoln.h:1414
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
Definition: HMWSoln.cpp:1104
vector< double > m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
Definition: HMWSoln.h:1772
vector< double > m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1821
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
Definition: HMWSoln.cpp:231
vector< double > m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1435
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
Definition: HMWSoln.cpp:211
double IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
Definition: HMWSoln.h:1876
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
Definition: HMWSoln.h:1635
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature.
Definition: HMWSoln.cpp:4036
vector< double > m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
Definition: HMWSoln.h:1761
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
Definition: HMWSoln.h:1358
void calcMCCutoffParams_()
Calculate molality cut-off parameters.
Definition: HMWSoln.cpp:1526
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer's papers.
Definition: HMWSoln.h:1502
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:2830
vector< double > m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1492
vector< double > m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1441
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
Definition: HMWSoln.h:1382
vector< double > m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
Definition: HMWSoln.h:1735
HMWSoln(const string &inputFile="", const string &id="")
Construct and initialize an HMWSoln ThermoPhase object directly from an input file.
Definition: HMWSoln.cpp:41
vector< double > m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:1698
vector< double > m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1533
vector< double > m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1835
vector< double > m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1560
vector< double > m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1557
vector< double > m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1530
vector< double > m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1838
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
Definition: HMWSoln.h:1629
double s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:4089
vector< double > m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Definition: HMWSoln.h:1592
void getParameters(AnyMap &phaseNode) const override
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
Definition: HMWSoln.cpp:752
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
Definition: HMWSoln.cpp:576
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure.
Definition: HMWSoln.cpp:4066
vector< double > m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition: HMWSoln.h:1426
void printCoeffs() const
Print out all of the input Pitzer coefficients.
Definition: HMWSoln.cpp:3963
void getActivityConcentrations(double *c) const override
This method returns an array of generalized activity concentrations.
Definition: HMWSoln.cpp:159
vector< double > m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1588
vector< double > m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
Definition: HMWSoln.h:1768
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
Definition: HMWSoln.cpp:2311
void s_update_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition: HMWSoln.cpp:1245
vector< double > m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1468
vector< int > CROP_speciesCropped_
This is a boolean-type vector indicating whether a species's activity coefficient is in the cropped r...
Definition: HMWSoln.h:1915
vector< double > m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1517
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_V.
Definition: HMWSoln.cpp:1081
vector< double > m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1805
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
Definition: HMWSoln.cpp:285
vector< double > m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1824
vector< double > m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1811
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
Definition: HMWSoln.cpp:1450
vector< double > m_BMX_IJ
Intermediate variable called BMX in Pitzer's paper.
Definition: HMWSoln.h:1776
vector< double > m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P.
Definition: HMWSoln.h:1729
vector< double > m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1851
vector< double > m_BMX_IJ_P
Derivative of BMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1785
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition: HMWSoln.cpp:130
vector< double > m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:1688
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1616
void getUnscaledMolalityActivityCoefficients(double *acMolality) const override
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
Definition: HMWSoln.cpp:198
void setA_Debye(double A)
Set the A_Debye parameter.
Definition: HMWSoln.cpp:539
vector< double > IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
Definition: HMWSoln.h:1862
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model.
Definition: HMWSoln.h:1353
vector< double > m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T.
Definition: HMWSoln.h:1713
vector< double > m_BMX_IJ_L
Derivative of BMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1779
virtual double relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
Definition: HMWSoln.cpp:60
vector< double > m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1438
vector< double > m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1444
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients' dependence on the temperature.
Definition: HMWSoln.cpp:1562
double s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:4081
vector< double > m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1554
vector< double > m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1486
vector< double > m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:1678
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients.
Definition: HMWSoln.cpp:3365
vector< double > m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
Definition: HMWSoln.h:1802
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
Definition: HMWSoln.h:1626
vector< double > m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer's paper.
Definition: HMWSoln.h:1789
vector< double > m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1563
vector< double > m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
Definition: HMWSoln.h:1832
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
Definition: HMWSoln.cpp:3894
vector< int > m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species.
Definition: HMWSoln.h:1746
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
Definition: HMWSoln.h:1691
virtual double relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
Definition: HMWSoln.cpp:72
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_L.
Definition: HMWSoln.cpp:1070
double entropy_mole() const override
Molar entropy. Units: J/kmol/K.
Definition: HMWSoln.cpp:112
PDSS * m_waterSS
Water standard state calculator.
Definition: HMWSoln.h:1420
void calcDensity() override
Calculate the density of the mixture using the partial molar volumes and mole fractions as input.
Definition: HMWSoln.cpp:142
void calc_lambdas(double is) const
Calculate the lambda interactions.
Definition: HMWSoln.cpp:3850
double elambda1[17]
This is elambda1, MEC.
Definition: HMWSoln.h:1752
vector< double > m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1465
vector< double > m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:1705
double s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
Definition: HMWSoln.cpp:4103
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: HMWSoln.h:1362
vector< double > m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T.
Definition: HMWSoln.h:1709
vector< double > m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1483
vector< double > m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1795
vector< double > m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
Definition: HMWSoln.h:1828
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients.
Definition: HMWSoln.cpp:2339
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
Definition: HMWSoln.cpp:1808
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition: HMWSoln.cpp:124
vector< double > m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1841
void getActivities(double *ac) const override
Get the array of non-dimensional activities at the current solution temperature, pressure,...
Definition: HMWSoln.cpp:182
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature as a function of tem...
Definition: HMWSoln.cpp:1014
vector< double > m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1792
vector< double > m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
Definition: HMWSoln.h:1815
vector< double > m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1782
vector< double > m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT.
Definition: HMWSoln.h:1717
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
Definition: HMWSoln.h:1632
double IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process;
Definition: HMWSoln.h:1866
double IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
Definition: HMWSoln.h:1869
vector< double > m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1808
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:2858
vector< double > m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1848
void getPartialMolarCp(double *cpbar) const override
Return an array of partial molar heat capacities for the species in the mixture.
Definition: HMWSoln.cpp:298
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1650
vector< double > m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
Definition: HMWSoln.h:1596
void initLengths()
Initialize all of the species-dependent lengths in the object.
Definition: HMWSoln.cpp:1132
vector< double > m_g2func_IJ
This is the value of g2(x2) in Pitzer's papers. Vector index is counterIJ.
Definition: HMWSoln.h:1764
double gibbs_mole() const override
Molar Gibbs function. Units: J/kmol.
Definition: HMWSoln.cpp:118
double standardConcentration(size_t k=0) const override
Return the standard concentration for the kth species.
Definition: HMWSoln.cpp:172
vector< double > m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:1524
vector< double > m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1798
vector< double > m_CMX_IJ
Intermediate variable called CMX in Pitzer's paper.
Definition: HMWSoln.h:1845
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Definition: HMWSoln.cpp:3341
double MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process;
Definition: HMWSoln.h:1892
double satPressure(double T) override
Get the saturation pressure for a given temperature.
Definition: HMWSoln.cpp:318
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
Definition: HMWSoln.h:1365
vector< double > m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:1489
double elambda[17]
This is elambda, MEC.
Definition: HMWSoln.h:1749
vector< double > m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P.
Definition: HMWSoln.h:1725
vector< double > m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
Definition: HMWSoln.h:1658
vector< double > m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
Definition: HMWSoln.h:1600
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition: HMWSoln.cpp:3920
void calcMolalitiesCropped() const
Calculate the cropped molalities.
Definition: HMWSoln.cpp:1309
vector< double > m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1527
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_J.
Definition: HMWSoln.cpp:1092
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature.
Definition: HMWSoln.cpp:4051
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
Definition: HMWSoln.cpp:250
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer's papers.
Definition: HMWSoln.h:1476
vector< double > m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:1854
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
Definition: HMWSoln.cpp:4021
vector< double > m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1462
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
Definition: HMWSoln.h:1542
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the Debye Huckel constant as a function of temperature and pressure.
Definition: HMWSoln.cpp:982
virtual double dA_DebyedP_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to pressure, as a function of tempe...
Definition: HMWSoln.cpp:1038
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
Definition: HMWSoln.h:1738
vector< double > m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:1818
vector< double > m_gamma_tmp
Intermediate storage of the activity coefficient itself.
Definition: HMWSoln.h:1858
double s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:4096
Error thrown for problems processing information contained in an AnyMap or AnyValue.
Definition: AnyMap.h:738
double m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
size_t m_indexCLM
Index of the phScale species.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
double m_xmolSolventMIN
In any molality implementation, it makes sense to have a minimum solvent mole fraction requirement,...
vector< double > m_molalities
Current value of the molalities of the species in the phase.
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
void setMoleFSolventMin(double xmolSolventMIN)
Sets the minimum mole fraction in the molality formulation.
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
double m_weightSolvent
Molecular weight of the Solvent.
Class for the liquid water pressure dependent standard state.
Definition: PDSS_Water.h:50
virtual double satPressure(double T)
saturation pressure
Definition: PDSS.cpp:157
virtual void setState_TP(double temp, double pres)
Set the internal temperature and pressure.
Definition: PDSS.cpp:152
void assignDensity(const double density_)
Set the internally stored constant density (kg/m^3) of the phase.
Definition: Phase.cpp:597
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
Definition: Phase.cpp:289
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition: Phase.h:822
size_t m_kk
Number of species in the phase.
Definition: Phase.h:842
double temperature() const
Temperature (K).
Definition: Phase.h:562
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:655
string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:142
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:434
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:129
double moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:439
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition: Phase.h:761
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:616
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition: Phase.cpp:856
virtual double molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:581
double 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:538
virtual double thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
Definition: ThermoPhase.h:580
virtual void getParameters(AnyMap &phaseNode) const
Store the parameters of a ThermoPhase object such that an identical one could be reconstructed using ...
double RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:1062
virtual double isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
Definition: ThermoPhase.h:569
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
AnyMap m_input
Data supplied via setParameters.
Definition: ThermoPhase.h:1966
double pressure() const override
Returns the current pressure of the phase.
void getEntropy_R(double *sr) const override
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
void getStandardChemPotentials(double *mu) const override
Get the array of chemical potentials at unit activity for the species at their standard states at the...
void getCp_R(double *cpr) const override
Get the nondimensional Heat Capacities at constant pressure for the species standard states at the cu...
void getEnthalpy_RT(double *hrt) const override
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
void getStandardVolumes(double *vol) const override
Get the molar volumes of the species standard states at the current T and P of the solution.
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (double) with the given id.
Definition: ValueCache.h:161
int getId()
Get a unique id for a cached value.
Definition: ValueCache.cpp:21
Header file for a common definitions used in electrolytes thermodynamics.
bool caseInsensitiveEquals(const string &input, const string &test)
Case insensitive equality predicate.
#define AssertTrace(expr)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:247
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:278
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:175
int sign(T x)
Sign of a number.
Definition: global.h:336
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer.
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:158
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
Contains declarations for string manipulation functions within Cantera.
A cached property value and the state at which it was evaluated.
Definition: ValueCache.h:33
T value
The value of the cached property.
Definition: ValueCache.h:113
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
Definition: ValueCache.h:39