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