Cantera  2.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
8  * obeys the Pitzer formulation for nonideality using molality-based
9  * standard states.
10  *
11  * This version of the code was modified to have the binary Beta2 Pitzer
12  * parameter consistent with the temperature expansions used for Beta0,
13  * Beta1, and Cphi.(CFJC, SNL)
14  */
15 /*
16  * Copyright (2006) Sandia Corporation. Under the terms of
17  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
18  * U.S. Government retains certain rights in this software.
19  */
20 
21 #include "cantera/thermo/HMWSoln.h"
22 #include "cantera/thermo/ThermoFactory.h"
23 #include "cantera/thermo/PDSS_Water.h"
24 #include "cantera/thermo/electrolytes.h"
25 #include "cantera/base/stringUtils.h"
26 #include <cstdio>
27 
28 namespace Cantera
29 {
30 
31 HMWSoln::HMWSoln() :
32  m_formPitzer(PITZERFORM_BASE),
33  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
34  m_formGC(2),
35  m_IionicMolality(0.0),
36  m_maxIionicStrength(100.0),
37  m_TempPitzerRef(298.15),
38  m_IionicMolalityStoich(0.0),
39  m_form_A_Debye(A_DEBYE_WATER),
40  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
41  m_waterSS(0),
42  m_densWaterSS(1000.),
43  m_waterProps(0),
44  m_molalitiesAreCropped(false),
45  IMS_typeCutoff_(0),
46  IMS_X_o_cutoff_(0.2),
47  IMS_gamma_o_min_(1.0E-5),
48  IMS_gamma_k_min_(10.0),
49  IMS_cCut_(0.05),
50  IMS_slopefCut_(0.6),
51  IMS_slopegCut_(0.0),
52  IMS_dfCut_(0.0),
53  IMS_efCut_(0.0),
54  IMS_afCut_(0.0),
55  IMS_bfCut_(0.0),
56  IMS_dgCut_(0.0),
57  IMS_egCut_(0.0),
58  IMS_agCut_(0.0),
59  IMS_bgCut_(0.0),
60  MC_X_o_cutoff_(0.0),
61  MC_X_o_min_(0.0),
62  MC_slopepCut_(0.0),
63  MC_dpCut_(0.0),
64  MC_epCut_(0.0),
65  MC_apCut_(0.0),
66  MC_bpCut_(0.0),
67  MC_cpCut_(0.0),
68  CROP_ln_gamma_o_min(-6.0),
69  CROP_ln_gamma_o_max(3.0),
70  CROP_ln_gamma_k_min(-5.0),
71  CROP_ln_gamma_k_max(15.0),
72  m_last_is(-1.0),
73  m_debugCalc(0)
74 {
75  for (size_t i = 0; i < 17; i++) {
76  elambda[i] = 0.0;
77  elambda1[i] = 0.0;
78  }
79 }
80 
81 HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
82  m_formPitzer(PITZERFORM_BASE),
83  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
84  m_formGC(2),
85  m_IionicMolality(0.0),
86  m_maxIionicStrength(100.0),
87  m_TempPitzerRef(298.15),
88  m_IionicMolalityStoich(0.0),
89  m_form_A_Debye(A_DEBYE_WATER),
90  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
91  m_waterSS(0),
92  m_densWaterSS(1000.),
93  m_waterProps(0),
94  m_molalitiesAreCropped(false),
95  IMS_typeCutoff_(0),
96  IMS_X_o_cutoff_(0.2),
97  IMS_gamma_o_min_(1.0E-5),
98  IMS_gamma_k_min_(10.0),
99  IMS_cCut_(0.05),
100  IMS_slopefCut_(0.6),
101  IMS_slopegCut_(0.0),
102  IMS_dfCut_(0.0),
103  IMS_efCut_(0.0),
104  IMS_afCut_(0.0),
105  IMS_bfCut_(0.0),
106  IMS_dgCut_(0.0),
107  IMS_egCut_(0.0),
108  IMS_agCut_(0.0),
109  IMS_bgCut_(0.0),
110  MC_X_o_cutoff_(0.0),
111  MC_X_o_min_(0.0),
112  MC_slopepCut_(0.0),
113  MC_dpCut_(0.0),
114  MC_epCut_(0.0),
115  MC_apCut_(0.0),
116  MC_bpCut_(0.0),
117  MC_cpCut_(0.0),
118  CROP_ln_gamma_o_min(-6.0),
119  CROP_ln_gamma_o_max(3.0),
120  CROP_ln_gamma_k_min(-5.0),
121  CROP_ln_gamma_k_max(15.0),
122  m_last_is(-1.0),
123  m_debugCalc(0)
124 {
125  for (int i = 0; i < 17; i++) {
126  elambda[i] = 0.0;
127  elambda1[i] = 0.0;
128  }
129  initThermoFile(inputFile, id_);
130 }
131 
132 HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
133  m_formPitzer(PITZERFORM_BASE),
134  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
135  m_formGC(2),
136  m_IionicMolality(0.0),
137  m_maxIionicStrength(100.0),
138  m_TempPitzerRef(298.15),
139  m_IionicMolalityStoich(0.0),
140  m_form_A_Debye(A_DEBYE_WATER),
141  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
142  m_waterSS(0),
143  m_densWaterSS(1000.),
144  m_waterProps(0),
145  m_molalitiesAreCropped(false),
146  IMS_typeCutoff_(0),
147  IMS_X_o_cutoff_(0.2),
148  IMS_gamma_o_min_(1.0E-5),
149  IMS_gamma_k_min_(10.0),
150  IMS_cCut_(0.05),
151  IMS_slopefCut_(0.6),
152  IMS_slopegCut_(0.0),
153  IMS_dfCut_(0.0),
154  IMS_efCut_(0.0),
155  IMS_afCut_(0.0),
156  IMS_bfCut_(0.0),
157  IMS_dgCut_(0.0),
158  IMS_egCut_(0.0),
159  IMS_agCut_(0.0),
160  IMS_bgCut_(0.0),
161  MC_X_o_cutoff_(0.0),
162  MC_X_o_min_(0.0),
163  MC_slopepCut_(0.0),
164  MC_dpCut_(0.0),
165  MC_epCut_(0.0),
166  MC_apCut_(0.0),
167  MC_bpCut_(0.0),
168  MC_cpCut_(0.0),
169  CROP_ln_gamma_o_min(-6.0),
170  CROP_ln_gamma_o_max(3.0),
171  CROP_ln_gamma_k_min(-5.0),
172  CROP_ln_gamma_k_max(15.0),
173  m_last_is(-1.0),
174  m_debugCalc(0)
175 {
176  for (int i = 0; i < 17; i++) {
177  elambda[i] = 0.0;
178  elambda1[i] = 0.0;
179  }
180  importPhase(*findXMLPhase(&phaseRoot, id_), this);
181 }
182 
183 HMWSoln::HMWSoln(const HMWSoln& b) :
184  m_formPitzer(PITZERFORM_BASE),
185  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
186  m_formGC(2),
187  m_IionicMolality(0.0),
188  m_maxIionicStrength(100.0),
189  m_TempPitzerRef(298.15),
190  m_IionicMolalityStoich(0.0),
191  m_form_A_Debye(A_DEBYE_WATER),
192  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
193  m_waterSS(0),
194  m_densWaterSS(1000.),
195  m_waterProps(0),
196  m_molalitiesAreCropped(false),
197  IMS_typeCutoff_(0),
198  IMS_X_o_cutoff_(0.2),
199  IMS_gamma_o_min_(1.0E-5),
200  IMS_gamma_k_min_(10.0),
201  IMS_cCut_(0.05),
202  IMS_slopefCut_(0.6),
203  IMS_slopegCut_(0.0),
204  IMS_dfCut_(0.0),
205  IMS_efCut_(0.0),
206  IMS_afCut_(0.0),
207  IMS_bfCut_(0.0),
208  IMS_dgCut_(0.0),
209  IMS_egCut_(0.0),
210  IMS_agCut_(0.0),
211  IMS_bgCut_(0.0),
212  MC_X_o_cutoff_(0.0),
213  MC_X_o_min_(0.0),
214  MC_slopepCut_(0.0),
215  MC_dpCut_(0.0),
216  MC_epCut_(0.0),
217  MC_apCut_(0.0),
218  MC_bpCut_(0.0),
219  MC_cpCut_(0.0),
220  CROP_ln_gamma_o_min(-6.0),
221  CROP_ln_gamma_o_max(3.0),
222  CROP_ln_gamma_k_min(-5.0),
223  CROP_ln_gamma_k_max(15.0),
224  m_last_is(-1.0),
225  m_debugCalc(0)
226 {
227  /*
228  * Use the assignment operator to do the brunt
229  * of the work for the copy constructor.
230  */
231  *this = b;
232 }
233 
234 HMWSoln& HMWSoln::operator=(const HMWSoln& b)
235 {
236  if (&b != this) {
237  MolalityVPSSTP::operator=(b);
238  m_formPitzer = b.m_formPitzer;
239  m_formPitzerTemp = b.m_formPitzerTemp;
240  m_formGC = b.m_formGC;
241  m_Aionic = b.m_Aionic;
242  m_IionicMolality = b.m_IionicMolality;
243  m_maxIionicStrength = b.m_maxIionicStrength;
244  m_TempPitzerRef = b.m_TempPitzerRef;
245  m_IionicMolalityStoich= b.m_IionicMolalityStoich;
246  m_form_A_Debye = b.m_form_A_Debye;
247  m_A_Debye = b.m_A_Debye;
248 
249  // This is an internal shallow copy of the PDSS_Water pointer
250  m_waterSS = providePDSS(0);
251  if (!m_waterSS) {
252  throw CanteraError("HMWSoln::operator=()", "Dynamic cast to PDSS_Water failed");
253  }
254 
255  m_densWaterSS = b.m_densWaterSS;
256 
257  delete m_waterProps;
258  m_waterProps = 0;
259  if (b.m_waterProps) {
260  m_waterProps = new WaterProps(dynamic_cast<PDSS_Water*>(m_waterSS));
261  }
262 
263  m_pp = b.m_pp;
264  m_tmpV = b.m_tmpV;
265  m_speciesCharge_Stoich= b.m_speciesCharge_Stoich;
266  m_Beta0MX_ij = b.m_Beta0MX_ij;
267  m_Beta0MX_ij_L = b.m_Beta0MX_ij_L;
268  m_Beta0MX_ij_LL = b.m_Beta0MX_ij_LL;
269  m_Beta0MX_ij_P = b.m_Beta0MX_ij_P;
270  m_Beta0MX_ij_coeff = b.m_Beta0MX_ij_coeff;
271  m_Beta1MX_ij = b.m_Beta1MX_ij;
272  m_Beta1MX_ij_L = b.m_Beta1MX_ij_L;
273  m_Beta1MX_ij_LL = b.m_Beta1MX_ij_LL;
274  m_Beta1MX_ij_P = b.m_Beta1MX_ij_P;
275  m_Beta1MX_ij_coeff = b.m_Beta1MX_ij_coeff;
276  m_Beta2MX_ij = b.m_Beta2MX_ij;
277  m_Beta2MX_ij_L = b.m_Beta2MX_ij_L;
278  m_Beta2MX_ij_LL = b.m_Beta2MX_ij_LL;
279  m_Beta2MX_ij_P = b.m_Beta2MX_ij_P;
280  m_Beta2MX_ij_coeff = b.m_Beta2MX_ij_coeff;
281  m_Alpha1MX_ij = b.m_Alpha1MX_ij;
282  m_Alpha2MX_ij = b.m_Alpha2MX_ij;
283  m_CphiMX_ij = b.m_CphiMX_ij;
284  m_CphiMX_ij_L = b.m_CphiMX_ij_L;
285  m_CphiMX_ij_LL = b.m_CphiMX_ij_LL;
286  m_CphiMX_ij_P = b.m_CphiMX_ij_P;
287  m_CphiMX_ij_coeff = b.m_CphiMX_ij_coeff;
288  m_Theta_ij = b.m_Theta_ij;
289  m_Theta_ij_L = b.m_Theta_ij_L;
290  m_Theta_ij_LL = b.m_Theta_ij_LL;
291  m_Theta_ij_P = b.m_Theta_ij_P;
292  m_Theta_ij_coeff = b.m_Theta_ij_coeff;
293  m_Psi_ijk = b.m_Psi_ijk;
294  m_Psi_ijk_L = b.m_Psi_ijk_L;
295  m_Psi_ijk_LL = b.m_Psi_ijk_LL;
296  m_Psi_ijk_P = b.m_Psi_ijk_P;
297  m_Psi_ijk_coeff = b.m_Psi_ijk_coeff;
298  m_Lambda_nj = b.m_Lambda_nj;
299  m_Lambda_nj_L = b.m_Lambda_nj_L;
300  m_Lambda_nj_LL = b.m_Lambda_nj_LL;
301  m_Lambda_nj_P = b.m_Lambda_nj_P;
302  m_Lambda_nj_coeff = b.m_Lambda_nj_coeff;
303 
304  m_Mu_nnn = b.m_Mu_nnn;
305  m_Mu_nnn_L = b.m_Mu_nnn_L;
306  m_Mu_nnn_LL = b.m_Mu_nnn_LL;
307  m_Mu_nnn_P = b.m_Mu_nnn_P;
308  m_Mu_nnn_coeff = b.m_Mu_nnn_coeff;
309 
310  m_lnActCoeffMolal_Scaled = b.m_lnActCoeffMolal_Scaled;
311  m_lnActCoeffMolal_Unscaled = b.m_lnActCoeffMolal_Unscaled;
312 
313  m_dlnActCoeffMolaldT_Scaled = b.m_dlnActCoeffMolaldT_Scaled;
314  m_dlnActCoeffMolaldT_Unscaled = b.m_dlnActCoeffMolaldT_Unscaled;
315 
316  m_d2lnActCoeffMolaldT2_Scaled = b.m_d2lnActCoeffMolaldT2_Scaled;
317  m_d2lnActCoeffMolaldT2_Unscaled= b.m_d2lnActCoeffMolaldT2_Unscaled;
318 
319  m_dlnActCoeffMolaldP_Scaled = b.m_dlnActCoeffMolaldP_Scaled;
320  m_dlnActCoeffMolaldP_Unscaled = b.m_dlnActCoeffMolaldP_Unscaled;
321 
322  m_molalitiesCropped = b.m_molalitiesCropped;
323  m_molalitiesAreCropped = b.m_molalitiesAreCropped;
324  m_CounterIJ = b.m_CounterIJ;
325 
326  m_gfunc_IJ = b.m_gfunc_IJ;
327  m_g2func_IJ = b.m_g2func_IJ;
328  m_hfunc_IJ = b.m_hfunc_IJ;
329  m_h2func_IJ = b.m_h2func_IJ;
330  m_BMX_IJ = b.m_BMX_IJ;
331  m_BMX_IJ_L = b.m_BMX_IJ_L;
332  m_BMX_IJ_LL = b.m_BMX_IJ_LL;
333  m_BMX_IJ_P = b.m_BMX_IJ_P;
334  m_BprimeMX_IJ = b.m_BprimeMX_IJ;
335  m_BprimeMX_IJ_L = b.m_BprimeMX_IJ_L;
336  m_BprimeMX_IJ_LL = b.m_BprimeMX_IJ_LL;
337  m_BprimeMX_IJ_P = b.m_BprimeMX_IJ_P;
338  m_BphiMX_IJ = b.m_BphiMX_IJ;
339  m_BphiMX_IJ_L = b.m_BphiMX_IJ_L;
340  m_BphiMX_IJ_LL = b.m_BphiMX_IJ_LL;
341  m_BphiMX_IJ_P = b.m_BphiMX_IJ_P;
342  m_Phi_IJ = b.m_Phi_IJ;
343  m_Phi_IJ_L = b.m_Phi_IJ_L;
344  m_Phi_IJ_LL = b.m_Phi_IJ_LL;
345  m_Phi_IJ_P = b.m_Phi_IJ_P;
346  m_Phiprime_IJ = b.m_Phiprime_IJ;
347  m_PhiPhi_IJ = b.m_PhiPhi_IJ;
348  m_PhiPhi_IJ_L = b.m_PhiPhi_IJ_L;
349  m_PhiPhi_IJ_LL = b.m_PhiPhi_IJ_LL;
350  m_PhiPhi_IJ_P = b.m_PhiPhi_IJ_P;
351  m_CMX_IJ = b.m_CMX_IJ;
352  m_CMX_IJ_L = b.m_CMX_IJ_L;
353  m_CMX_IJ_LL = b.m_CMX_IJ_LL;
354  m_CMX_IJ_P = b.m_CMX_IJ_P;
355  m_gamma_tmp = b.m_gamma_tmp;
356 
357  IMS_lnActCoeffMolal_ = b.IMS_lnActCoeffMolal_;
358  IMS_typeCutoff_ = b.IMS_typeCutoff_;
359  IMS_X_o_cutoff_ = b.IMS_X_o_cutoff_;
360  IMS_gamma_o_min_ = b.IMS_gamma_o_min_;
361  IMS_gamma_k_min_ = b.IMS_gamma_k_min_;
362  IMS_cCut_ = b.IMS_cCut_;
363  IMS_slopefCut_ = b.IMS_slopefCut_;
364  IMS_dfCut_ = b.IMS_dfCut_;
365  IMS_efCut_ = b.IMS_efCut_;
366  IMS_afCut_ = b.IMS_afCut_;
367  IMS_bfCut_ = b.IMS_bfCut_;
368  IMS_slopegCut_ = b.IMS_slopegCut_;
369  IMS_dgCut_ = b.IMS_dgCut_;
370  IMS_egCut_ = b.IMS_egCut_;
371  IMS_agCut_ = b.IMS_agCut_;
372  IMS_bgCut_ = b.IMS_bgCut_;
373  MC_X_o_cutoff_ = b.MC_X_o_cutoff_;
374  MC_X_o_min_ = b.MC_X_o_min_;
375  MC_slopepCut_ = b.MC_slopepCut_;
376  MC_dpCut_ = b.MC_dpCut_;
377  MC_epCut_ = b.MC_epCut_;
378  MC_apCut_ = b.MC_apCut_;
379  MC_bpCut_ = b.MC_bpCut_;
380  MC_cpCut_ = b.MC_cpCut_;
381  CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min;
382  CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max;
383  CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min;
384  CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max;
385  CROP_speciesCropped_ = b.CROP_speciesCropped_;
386 
387  m_debugCalc = b.m_debugCalc;
388  }
389  return *this;
390 }
391 
392 HMWSoln::~HMWSoln()
393 {
394  delete m_waterProps;
395  m_waterProps = 0;
396 }
397 
398 ThermoPhase* HMWSoln::duplMyselfAsThermoPhase() const
399 {
400  return new HMWSoln(*this);
401 }
402 
403 int HMWSoln::eosType() const
404 {
405  int res;
406  switch (m_formGC) {
407  case 0:
408  res = cHMWSoln0;
409  break;
410  case 1:
411  res = cHMWSoln1;
412  break;
413  case 2:
414  res = cHMWSoln2;
415  break;
416  default:
417  throw CanteraError("eosType", "Unknown type");
418  }
419  return res;
420 }
421 
422 //
423 // -------- Molar Thermodynamic Properties of the Solution ---------------
424 //
425 doublereal HMWSoln::enthalpy_mole() const
426 {
427  getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
428  getMoleFractions(DATA_PTR(m_pp));
429  return mean_X(m_tmpV);
430 }
431 
432 doublereal HMWSoln::relative_enthalpy() const
433 {
434  getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
435  double hbar = mean_X(m_tmpV);
436  getEnthalpy_RT(DATA_PTR(m_gamma_tmp));
437  double RT = GasConstant * temperature();
438  for (size_t k = 0; k < m_kk; k++) {
439  m_gamma_tmp[k] *= RT;
440  }
441  double h0bar = mean_X(m_gamma_tmp);
442  return hbar - h0bar;
443 }
444 
445 doublereal HMWSoln::relative_molal_enthalpy() const
446 {
447  double L = relative_enthalpy();
448  getMoleFractions(DATA_PTR(m_tmpV));
449  double xanion = 0.0;
450  size_t kcation = npos;
451  double xcation = 0.0;
452  size_t kanion = npos;
453  for (size_t k = 0; k < m_kk; k++) {
454  if (charge(k) > 0.0) {
455  if (m_tmpV[k] > xanion) {
456  xanion = m_tmpV[k];
457  kanion = k;
458  }
459  } else if (charge(k) < 0.0) {
460  if (m_tmpV[k] > xcation) {
461  xcation = m_tmpV[k];
462  kcation = k;
463  }
464  }
465  }
466  if (kcation == npos || kanion == npos) {
467  return L;
468  }
469  double xuse = xcation;
470  double factor = 1;
471  if (xanion < xcation) {
472  xuse = xanion;
473  if (charge(kcation) != 1.0) {
474  factor = charge(kcation);
475  }
476  } else {
477  if (charge(kanion) != 1.0) {
478  factor = charge(kanion);
479  }
480  }
481  xuse = xuse / factor;
482  return L / xuse;
483 }
484 
485 doublereal HMWSoln::entropy_mole() const
486 {
487  getPartialMolarEntropies(DATA_PTR(m_tmpV));
488  return mean_X(m_tmpV);
489 }
490 
491 doublereal HMWSoln::gibbs_mole() const
492 {
493  getChemPotentials(DATA_PTR(m_tmpV));
494  return mean_X(m_tmpV);
495 }
496 
497 doublereal HMWSoln::cp_mole() const
498 {
499  getPartialMolarCp(DATA_PTR(m_tmpV));
500  return mean_X(m_tmpV);
501 }
502 
503 doublereal HMWSoln::cv_mole() const
504 {
505  double kappa_t = isothermalCompressibility();
506  double beta = thermalExpansionCoeff();
507  double cp = cp_mole();
508  double tt = temperature();
509  double molarV = molarVolume();
510  return cp - beta * beta * tt * molarV / kappa_t;
511 }
512 
513 //
514 // ------- Mechanical Equation of State Properties ------------------------
515 //
516 
517 doublereal HMWSoln::pressure() const
518 {
519  return m_Pcurrent;
520 }
521 
522 void HMWSoln::setPressure(doublereal p)
523 {
524  setState_TP(temperature(), p);
525 }
526 
527 void HMWSoln::calcDensity()
528 {
529  static const int cacheId = m_cache.getId();
530  CachedScalar cached = m_cache.getScalar(cacheId);
531  if(cached.validate(temperature(), pressure(), stateMFNumber())) {
532  return;
533  }
534 
535  double* vbar = &m_pp[0];
536  getPartialMolarVolumes(vbar);
537  double* x = &m_tmpV[0];
538  getMoleFractions(x);
539  doublereal vtotal = 0.0;
540  for (size_t i = 0; i < m_kk; i++) {
541  vtotal += vbar[i] * x[i];
542  }
543  doublereal dd = meanMolecularWeight() / vtotal;
544  Phase::setDensity(dd);
545 }
546 
547 double HMWSoln::density() const
548 {
549  return Phase::density();
550 }
551 
552 void HMWSoln::setDensity(const doublereal rho)
553 {
554  double dens_old = density();
555 
556  if (rho != dens_old) {
557  throw CanteraError("HMWSoln::setDensity",
558  "Density is not an independent variable");
559  }
560 }
561 
562 void HMWSoln::setMolarDensity(const doublereal rho)
563 {
564  throw CanteraError("HMWSoln::setMolarDensity",
565  "Density is not an independent variable");
566 }
567 
568 void HMWSoln::setTemperature(const doublereal temp)
569 {
570  setState_TP(temp, m_Pcurrent);
571 }
572 
573 void HMWSoln::setState_TP(doublereal temp, doublereal pres)
574 {
575  Phase::setTemperature(temp);
576  /*
577  * Store the current pressure
578  */
579  m_Pcurrent = pres;
580 
581  /*
582  * update the standard state thermo
583  * -> This involves calling the water function and setting the pressure
584  */
585  updateStandardStateThermo();
586 
587  /*
588  * Store the internal density of the water SS.
589  * Note, we would have to do this for all other
590  * species if they had pressure dependent properties.
591  */
592  m_densWaterSS = m_waterSS->density();
593  /*
594  * Calculate all of the other standard volumes
595  * -> note these are constant for now
596  */
597  calcDensity();
598 }
599 
600 //
601 // ------- Activities and Activity Concentrations
602 //
603 
604 void HMWSoln::getActivityConcentrations(doublereal* c) const
605 {
606  double cs_solvent = standardConcentration();
607  getActivities(c);
608  c[0] *= cs_solvent;
609  if (m_kk > 1) {
610  double cs_solute = standardConcentration(1);
611  for (size_t k = 1; k < m_kk; k++) {
612  c[k] *= cs_solute;
613  }
614  }
615 }
616 
617 doublereal HMWSoln::standardConcentration(size_t k) const
618 {
619  getStandardVolumes(DATA_PTR(m_tmpV));
620  double mvSolvent = m_tmpV[m_indexSolvent];
621  if (k > 0) {
622  return m_Mnaught / mvSolvent;
623  }
624  return 1.0 / mvSolvent;
625 }
626 
627 void HMWSoln::getUnitsStandardConc(double* uA, int k, int sizeUA) const
628 {
629  warn_deprecated("HMWSoln::getUnitsStandardConc",
630  "To be removed after Cantera 2.2.");
631  for (int i = 0; i < sizeUA; i++) {
632  if (i == 0) {
633  uA[0] = 1.0;
634  }
635  if (i == 1) {
636  uA[1] = -int(nDim());
637  }
638  if (i == 2) {
639  uA[2] = 0.0;
640  }
641  if (i == 3) {
642  uA[3] = 0.0;
643  }
644  if (i == 4) {
645  uA[4] = 0.0;
646  }
647  if (i == 5) {
648  uA[5] = 0.0;
649  }
650  }
651 }
652 
653 void HMWSoln::getActivities(doublereal* ac) const
654 {
655  updateStandardStateThermo();
656  /*
657  * Update the molality array, m_molalities()
658  * This requires an update due to mole fractions
659  */
660  s_update_lnMolalityActCoeff();
661  /*
662  * Now calculate the array of activities.
663  */
664  for (size_t k = 0; k < m_kk; k++) {
665  if (k != m_indexSolvent) {
666  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
667  }
668  }
669  double xmolSolvent = moleFraction(m_indexSolvent);
670  ac[m_indexSolvent] =
671  exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent;
672 }
673 
674 void HMWSoln::getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
675 {
676  updateStandardStateThermo();
677  A_Debye_TP(-1.0, -1.0);
678  s_update_lnMolalityActCoeff();
679  std::copy(m_lnActCoeffMolal_Unscaled.begin(), m_lnActCoeffMolal_Unscaled.end(), acMolality);
680  for (size_t k = 0; k < m_kk; k++) {
681  acMolality[k] = exp(acMolality[k]);
682  }
683 }
684 
685 //
686 // ------ Partial Molar Properties of the Solution -----------------
687 //
688 
689 void HMWSoln::getChemPotentials(doublereal* mu) const
690 {
691  double xx;
692  /*
693  * First get the standard chemical potentials in
694  * molar form.
695  * -> this requires updates of standard state as a function
696  * of T and P
697  */
698  getStandardChemPotentials(mu);
699  /*
700  * Update the activity coefficients
701  * This also updates the internal molality array.
702  */
703  s_update_lnMolalityActCoeff();
704  doublereal RT = GasConstant * temperature();
705  double xmolSolvent = moleFraction(m_indexSolvent);
706  for (size_t k = 0; k < m_kk; k++) {
707  if (m_indexSolvent != k) {
708  xx = std::max(m_molalities[k], SmallNumber);
709  mu[k] += RT * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
710  }
711  }
712  xx = std::max(xmolSolvent, SmallNumber);
713  mu[m_indexSolvent] +=
714  RT * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
715 }
716 
717 void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const
718 {
719  /*
720  * Get the nondimensional standard state enthalpies
721  */
722  getEnthalpy_RT(hbar);
723  /*
724  * dimensionalize it.
725  */
726  double T = temperature();
727  double RT = GasConstant * T;
728  for (size_t k = 0; k < m_kk; k++) {
729  hbar[k] *= RT;
730  }
731  /*
732  * Update the activity coefficients, This also update the
733  * internally stored molalities.
734  */
735  s_update_lnMolalityActCoeff();
736  s_update_dlnMolalityActCoeff_dT();
737  double RTT = RT * T;
738  for (size_t k = 0; k < m_kk; k++) {
739  hbar[k] -= RTT * m_dlnActCoeffMolaldT_Scaled[k];
740  }
741 }
742 
743 void HMWSoln::getPartialMolarEntropies(doublereal* sbar) const
744 {
745  /*
746  * Get the standard state entropies at the temperature
747  * and pressure of the solution.
748  */
749  getEntropy_R(sbar);
750  /*
751  * Dimensionalize the entropies
752  */
753  doublereal R = GasConstant;
754  for (size_t k = 0; k < m_kk; k++) {
755  sbar[k] *= R;
756  }
757  /*
758  * Update the activity coefficients, This also update the
759  * internally stored molalities.
760  */
761  s_update_lnMolalityActCoeff();
762  /*
763  * First we will add in the obvious dependence on the T
764  * term out front of the log activity term
765  */
766  doublereal mm;
767  for (size_t k = 0; k < m_kk; k++) {
768  if (k != m_indexSolvent) {
769  mm = std::max(SmallNumber, m_molalities[k]);
770  sbar[k] -= R * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
771  }
772  }
773  double xmolSolvent = moleFraction(m_indexSolvent);
774  mm = std::max(SmallNumber, xmolSolvent);
775  sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
776  /*
777  * Check to see whether activity coefficients are temperature
778  * dependent. If they are, then calculate the their temperature
779  * derivatives and add them into the result.
780  */
781  s_update_dlnMolalityActCoeff_dT();
782  double RT = R * temperature();
783  for (size_t k = 0; k < m_kk; k++) {
784  sbar[k] -= RT * m_dlnActCoeffMolaldT_Scaled[k];
785  }
786 }
787 
788 void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const
789 {
790  /*
791  * Get the standard state values in m^3 kmol-1
792  */
793  getStandardVolumes(vbar);
794  /*
795  * Update the derivatives wrt the activity coefficients.
796  */
797  s_update_lnMolalityActCoeff();
798  s_update_dlnMolalityActCoeff_dP();
799  double T = temperature();
800  double RT = GasConstant * T;
801  for (size_t k = 0; k < m_kk; k++) {
802  vbar[k] += RT * m_dlnActCoeffMolaldP_Scaled[k];
803  }
804 }
805 
806 void HMWSoln::getPartialMolarCp(doublereal* cpbar) const
807 {
808  /*
809  * Get the nondimensional Gibbs standard state of the
810  * species at the T and P of the solution.
811  */
812  getCp_R(cpbar);
813 
814  for (size_t k = 0; k < m_kk; k++) {
815  cpbar[k] *= GasConstant;
816  }
817  /*
818  * Update the activity coefficients, This also update the
819  * internally stored molalities.
820  */
821  s_update_lnMolalityActCoeff();
822  s_update_dlnMolalityActCoeff_dT();
823  s_update_d2lnMolalityActCoeff_dT2();
824  double T = temperature();
825  double RT = GasConstant * T;
826  double RTT = RT * T;
827  for (size_t k = 0; k < m_kk; k++) {
828  cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT_Scaled[k] +
829  RTT * m_d2lnActCoeffMolaldT2_Scaled[k]);
830  }
831 }
832 
833 /*
834  * -------------- Utilities -------------------------------
835  */
836 
837 doublereal HMWSoln::satPressure(doublereal t) {
838  double p_old = pressure();
839  double t_old = temperature();
840  double pres = m_waterSS->satPressure(t);
841  /*
842  * Set the underlying object back to its original state.
843  */
844  m_waterSS->setState_TP(t_old, p_old);
845  return pres;
846 }
847 
848 double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
849 {
850  double T = temperature();
851  double A;
852  if (tempArg != -1.0) {
853  T = tempArg;
854  }
855  double P = pressure();
856  if (presArg != -1.0) {
857  P = presArg;
858  }
859 
860  static const int cacheId = m_cache.getId();
861  CachedScalar cached = m_cache.getScalar(cacheId);
862  if(cached.validate(T, P)) {
863  return m_A_Debye;
864  }
865 
866  switch (m_form_A_Debye) {
867  case A_DEBYE_CONST:
868  A = m_A_Debye;
869  break;
870  case A_DEBYE_WATER:
871  A = m_waterProps->ADebye(T, P, 0);
872  m_A_Debye = A;
873  break;
874  default:
875  throw CanteraError("HMWSoln::A_Debye_TP", "shouldn't be here");
876  }
877  return A;
878 }
879 
880 double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
881 {
882  doublereal T = temperature();
883  if (tempArg != -1.0) {
884  T = tempArg;
885  }
886  doublereal P = pressure();
887  if (presArg != -1.0) {
888  P = presArg;
889  }
890  doublereal dAdT;
891  switch (m_form_A_Debye) {
892  case A_DEBYE_CONST:
893  dAdT = 0.0;
894  break;
895  case A_DEBYE_WATER:
896  dAdT = m_waterProps->ADebye(T, P, 1);
897  break;
898  default:
899  throw CanteraError("HMWSoln::dA_DebyedT_TP", "shouldn't be here");
900  }
901  return dAdT;
902 }
903 
904 double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
905 {
906  double T = temperature();
907  if (tempArg != -1.0) {
908  T = tempArg;
909  }
910  double P = pressure();
911  if (presArg != -1.0) {
912  P = presArg;
913  }
914 
915  double dAdP;
916  static const int cacheId = m_cache.getId();
917  CachedScalar cached = m_cache.getScalar(cacheId);
918  switch (m_form_A_Debye) {
919  case A_DEBYE_CONST:
920  dAdP = 0.0;
921  break;
922  case A_DEBYE_WATER:
923  if(cached.validate(T, P)) {
924  dAdP = cached.value;
925  } else {
926  dAdP = m_waterProps->ADebye(T, P, 3);
927  cached.value = dAdP;
928  }
929  break;
930  default:
931  throw CanteraError("HMWSoln::dA_DebyedP_TP", "shouldn't be here");
932  }
933  return dAdP;
934 }
935 
936 double HMWSoln::ADebye_L(double tempArg, double presArg) const
937 {
938  double dAdT = dA_DebyedT_TP();
939  double dAphidT = dAdT /3.0;
940  double T = temperature();
941  if (tempArg != -1.0) {
942  T = tempArg;
943  }
944  return dAphidT * (4.0 * GasConstant * T * T);
945 }
946 
947 double HMWSoln::ADebye_V(double tempArg, double presArg) const
948 {
949  double dAdP = dA_DebyedP_TP();
950  double dAphidP = dAdP /3.0;
951  double T = temperature();
952  if (tempArg != -1.0) {
953  T = tempArg;
954  }
955  return - dAphidP * (4.0 * GasConstant * T);
956 }
957 
958 double HMWSoln::ADebye_J(double tempArg, double presArg) const
959 {
960  double T = temperature();
961  if (tempArg != -1.0) {
962  T = tempArg;
963  }
964  double A_L = ADebye_L(T, presArg);
965  double d2 = d2A_DebyedT2_TP(T, presArg);
966  double d2Aphi = d2 / 3.0;
967  return 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi;
968 }
969 
970 double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
971 {
972  double T = temperature();
973  if (tempArg != -1.0) {
974  T = tempArg;
975  }
976  double P = pressure();
977  if (presArg != -1.0) {
978  P = presArg;
979  }
980  double d2AdT2;
981  switch (m_form_A_Debye) {
982  case A_DEBYE_CONST:
983  d2AdT2 = 0.0;
984  break;
985  case A_DEBYE_WATER:
986  d2AdT2 = m_waterProps->ADebye(T, P, 2);
987  break;
988  default:
989  throw CanteraError("HMWSoln::d2A_DebyedT2_TP", "shouldn't be here");
990  }
991  return d2AdT2;
992 }
993 
994 /*
995  * ---------- Other Property Functions
996  */
997 double HMWSoln::AionicRadius(int k) const
998 {
999  return m_Aionic[k];
1000 }
1001 
1002 /*
1003  * ------------ Private and Restricted Functions ------------------
1004  */
1005 
1006 void HMWSoln::initLengths()
1007 {
1008  /*
1009  * Resize lengths equal to the number of species in
1010  * the phase.
1011  */
1012  m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
1013  m_speciesSize.resize(m_kk);
1014  m_speciesCharge_Stoich.resize(m_kk, 0.0);
1015  m_Aionic.resize(m_kk, 0.0);
1016 
1017  m_pp.resize(m_kk, 0.0);
1018  m_tmpV.resize(m_kk, 0.0);
1019  m_molalitiesCropped.resize(m_kk, 0.0);
1020 
1021  size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
1022 
1023  /*
1024  * Figure out the size of the temperature coefficient
1025  * arrays
1026  */
1027  int TCoeffLength = 1;
1028  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1029  TCoeffLength = 2;
1030  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1031  TCoeffLength = 5;
1032  }
1033 
1034  m_Beta0MX_ij.resize(maxCounterIJlen, 0.0);
1035  m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0);
1036  m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0);
1037  m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0);
1038  m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1039 
1040  m_Beta1MX_ij.resize(maxCounterIJlen, 0.0);
1041  m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0);
1042  m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0);
1043  m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0);
1044  m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1045 
1046  m_Beta2MX_ij.resize(maxCounterIJlen, 0.0);
1047  m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0);
1048  m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0);
1049  m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0);
1050  m_Beta2MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1051 
1052  m_CphiMX_ij.resize(maxCounterIJlen, 0.0);
1053  m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0);
1054  m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0);
1055  m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0);
1056  m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1057 
1058  m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
1059  m_Alpha2MX_ij.resize(maxCounterIJlen, 12.0);
1060  m_Theta_ij.resize(maxCounterIJlen, 0.0);
1061  m_Theta_ij_L.resize(maxCounterIJlen, 0.0);
1062  m_Theta_ij_LL.resize(maxCounterIJlen, 0.0);
1063  m_Theta_ij_P.resize(maxCounterIJlen, 0.0);
1064  m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
1065 
1066  size_t n = m_kk*m_kk*m_kk;
1067  m_Psi_ijk.resize(n, 0.0);
1068  m_Psi_ijk_L.resize(n, 0.0);
1069  m_Psi_ijk_LL.resize(n, 0.0);
1070  m_Psi_ijk_P.resize(n, 0.0);
1071  m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0);
1072 
1073  m_Lambda_nj.resize(m_kk, m_kk, 0.0);
1074  m_Lambda_nj_L.resize(m_kk, m_kk, 0.0);
1075  m_Lambda_nj_LL.resize(m_kk, m_kk, 0.0);
1076  m_Lambda_nj_P.resize(m_kk, m_kk, 0.0);
1077  m_Lambda_nj_coeff.resize(TCoeffLength, m_kk * m_kk, 0.0);
1078 
1079  m_Mu_nnn.resize(m_kk, 0.0);
1080  m_Mu_nnn_L.resize(m_kk, 0.0);
1081  m_Mu_nnn_LL.resize(m_kk, 0.0);
1082  m_Mu_nnn_P.resize(m_kk, 0.0);
1083  m_Mu_nnn_coeff.resize(TCoeffLength, m_kk, 0.0);
1084 
1085  m_lnActCoeffMolal_Scaled.resize(m_kk, 0.0);
1086  m_dlnActCoeffMolaldT_Scaled.resize(m_kk, 0.0);
1087  m_d2lnActCoeffMolaldT2_Scaled.resize(m_kk, 0.0);
1088  m_dlnActCoeffMolaldP_Scaled.resize(m_kk, 0.0);
1089 
1090  m_lnActCoeffMolal_Unscaled.resize(m_kk, 0.0);
1091  m_dlnActCoeffMolaldT_Unscaled.resize(m_kk, 0.0);
1092  m_d2lnActCoeffMolaldT2_Unscaled.resize(m_kk, 0.0);
1093  m_dlnActCoeffMolaldP_Unscaled.resize(m_kk, 0.0);
1094 
1095  m_CounterIJ.resize(m_kk*m_kk, 0);
1096 
1097  m_gfunc_IJ.resize(maxCounterIJlen, 0.0);
1098  m_g2func_IJ.resize(maxCounterIJlen, 0.0);
1099  m_hfunc_IJ.resize(maxCounterIJlen, 0.0);
1100  m_h2func_IJ.resize(maxCounterIJlen, 0.0);
1101  m_BMX_IJ.resize(maxCounterIJlen, 0.0);
1102  m_BMX_IJ_L.resize(maxCounterIJlen, 0.0);
1103  m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1104  m_BMX_IJ_P.resize(maxCounterIJlen, 0.0);
1105  m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0);
1106  m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0);
1107  m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1108  m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0);
1109  m_BphiMX_IJ.resize(maxCounterIJlen, 0.0);
1110  m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0);
1111  m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1112  m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0);
1113  m_Phi_IJ.resize(maxCounterIJlen, 0.0);
1114  m_Phi_IJ_L.resize(maxCounterIJlen, 0.0);
1115  m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0);
1116  m_Phi_IJ_P.resize(maxCounterIJlen, 0.0);
1117  m_Phiprime_IJ.resize(maxCounterIJlen, 0.0);
1118  m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0);
1119  m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0);
1120  m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0);
1121  m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0);
1122  m_CMX_IJ.resize(maxCounterIJlen, 0.0);
1123  m_CMX_IJ_L.resize(maxCounterIJlen, 0.0);
1124  m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0);
1125  m_CMX_IJ_P.resize(maxCounterIJlen, 0.0);
1126 
1127  m_gamma_tmp.resize(m_kk, 0.0);
1128 
1129  IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
1130  CROP_speciesCropped_.resize(m_kk, 0);
1131 
1132  counterIJ_setup();
1133 }
1134 
1135 void HMWSoln::s_update_lnMolalityActCoeff() const
1136 {
1137  static const int cacheId = m_cache.getId();
1138  CachedScalar cached = m_cache.getScalar(cacheId);
1139  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
1140  return;
1141  }
1142 
1143  /*
1144  * Calculate the molalities. Currently, the molalities
1145  * may not be current with respect to the contents of the
1146  * State objects' data.
1147  */
1148  calcMolalities();
1149  /*
1150  * Calculate a cropped set of molalities that will be used
1151  * in all activity coefficient calculations.
1152  */
1153  calcMolalitiesCropped();
1154 
1155  /*
1156  * Calculate the stoichiometric ionic charge. This isn't used in the
1157  * Pitzer formulation.
1158  */
1159  m_IionicMolalityStoich = 0.0;
1160  for (size_t k = 0; k < m_kk; k++) {
1161  double z_k = charge(k);
1162  double zs_k1 = m_speciesCharge_Stoich[k];
1163  if (z_k == zs_k1) {
1164  m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
1165  } else {
1166  double zs_k2 = z_k - zs_k1;
1167  m_IionicMolalityStoich
1168  += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
1169  }
1170  }
1171 
1172  /*
1173  * Update the temperature dependence of the pitzer coefficients
1174  * and their derivatives
1175  */
1176  s_updatePitzer_CoeffWRTemp();
1177 
1178  /*
1179  * Calculate the IMS cutoff factors
1180  */
1181  s_updateIMS_lnMolalityActCoeff();
1182 
1183  /*
1184  * Now do the main calculation.
1185  */
1186  s_updatePitzer_lnMolalityActCoeff();
1187 
1188  double xmolSolvent = moleFraction(m_indexSolvent);
1189  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
1190  double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1191  double lnxs = log(xx);
1192 
1193  for (size_t k = 1; k < m_kk; k++) {
1194  CROP_speciesCropped_[k] = 0;
1195  m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
1196  if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
1197  CROP_speciesCropped_[k] = 2;
1198  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
1199  }
1200  if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
1201  // -1.0 and -1.5 caused multiple solutions
1202  CROP_speciesCropped_[k] = 2;
1203  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
1204  }
1205  }
1206  CROP_speciesCropped_[0] = 0;
1207  m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
1208  if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
1209  CROP_speciesCropped_[0] = 2;
1210  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
1211  }
1212  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
1213  CROP_speciesCropped_[0] = 2;
1214  // -0.5 caused multiple solutions
1215  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
1216  }
1217  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
1218  CROP_speciesCropped_[0] = 2;
1219  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
1220  }
1221 
1222  /*
1223  * Now do the pH Scaling
1224  */
1225  s_updateScaling_pHScaling();
1226 }
1227 
1228 void HMWSoln::calcMolalitiesCropped() const
1229 {
1230  doublereal Imax = 0.0;
1231  m_molalitiesAreCropped = false;
1232 
1233  for (size_t k = 0; k < m_kk; k++) {
1234  m_molalitiesCropped[k] = m_molalities[k];
1235  Imax = std::max(m_molalities[k] * charge(k) * charge(k), Imax);
1236  }
1237 
1238  int cropMethod = 1;
1239 
1240 
1241  if (cropMethod == 0) {
1242 
1243  /*
1244  * Quick return
1245  */
1246  if (Imax < m_maxIionicStrength) {
1247  return;
1248  }
1249 
1250  m_molalitiesAreCropped = true;
1251 
1252  for (size_t i = 1; i < (m_kk - 1); i++) {
1253  double charge_i = charge(i);
1254  double abs_charge_i = fabs(charge_i);
1255  if (charge_i == 0.0) {
1256  continue;
1257  }
1258  for (size_t j = (i+1); j < m_kk; j++) {
1259  double charge_j = charge(j);
1260  double abs_charge_j = fabs(charge_j);
1261  /*
1262  * Only loop over oppositely charge species
1263  */
1264  if (charge_i * charge_j < 0) {
1265  double Iac_max = m_maxIionicStrength;
1266 
1267  if (m_molalitiesCropped[i] > m_molalitiesCropped[j]) {
1268  Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
1269  if (Imax > Iac_max) {
1270  m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
1271  }
1272  Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
1273  if (Imax > Iac_max) {
1274  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
1275  }
1276  } else {
1277  Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
1278  if (Imax > Iac_max) {
1279  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
1280  }
1281  Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
1282  if (Imax > Iac_max) {
1283  m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
1284  }
1285  }
1286  }
1287  }
1288  }
1289 
1290  /*
1291  * Do this loop 10 times until we have achieved charge neutrality
1292  * in the cropped molalities
1293  */
1294  for (int times = 0; times< 10; times++) {
1295  double anion_charge = 0.0;
1296  double cation_charge = 0.0;
1297  size_t anion_contrib_max_i = npos;
1298  double anion_contrib_max = -1.0;
1299  size_t cation_contrib_max_i = npos;
1300  double cation_contrib_max = -1.0;
1301  for (size_t i = 0; i < m_kk; i++) {
1302  double charge_i = charge(i);
1303  if (charge_i < 0.0) {
1304  double anion_contrib = - m_molalitiesCropped[i] * charge_i;
1305  anion_charge += anion_contrib ;
1306  if (anion_contrib > anion_contrib_max) {
1307  anion_contrib_max = anion_contrib;
1308  anion_contrib_max_i = i;
1309  }
1310  } else if (charge_i > 0.0) {
1311  double cation_contrib = m_molalitiesCropped[i] * charge_i;
1312  cation_charge += cation_contrib ;
1313  if (cation_contrib > cation_contrib_max) {
1314  cation_contrib_max = cation_contrib;
1315  cation_contrib_max_i = i;
1316  }
1317  }
1318  }
1319  double total_charge = cation_charge - anion_charge;
1320  if (total_charge > 1.0E-8) {
1321  double desiredCrop = total_charge/charge(cation_contrib_max_i);
1322  double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
1323  if (desiredCrop < maxCrop) {
1324  m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
1325  break;
1326  } else {
1327  m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
1328  }
1329  } else if (total_charge < -1.0E-8) {
1330  double desiredCrop = total_charge/charge(anion_contrib_max_i);
1331  double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
1332  if (desiredCrop < maxCrop) {
1333  m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
1334  break;
1335  } else {
1336  m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
1337  }
1338  } else {
1339  break;
1340  }
1341  }
1342  }
1343 
1344  if (cropMethod == 1) {
1345  double* molF = DATA_PTR(m_gamma_tmp);
1346  getMoleFractions(molF);
1347  double xmolSolvent = molF[m_indexSolvent];
1348  if (xmolSolvent >= MC_X_o_cutoff_) {
1349  return;
1350  }
1351 
1352  m_molalitiesAreCropped = true;
1353 
1354  double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1355  double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1356  double denomInv = 1.0/ (m_Mnaught * p);
1357 
1358  for (size_t k = 0; k < m_kk; k++) {
1359  m_molalitiesCropped[k] = molF[k] * denomInv ;
1360  }
1361 
1362  // Do a further check to see if the Ionic strength is below a max value
1363  // Reduce the molalities to enforce this. Note, this algorithm preserves
1364  // the charge neutrality of the solution after cropping.
1365  double Itmp = 0.0;
1366  for (size_t k = 0; k < m_kk; k++) {
1367  Itmp += m_molalitiesCropped[k] * charge(k) * charge(k);
1368  }
1369  if (Itmp > m_maxIionicStrength) {
1370  double ratio = Itmp / m_maxIionicStrength;
1371  for (size_t k = 0; k < m_kk; k++) {
1372  if (charge(k) != 0.0) {
1373  m_molalitiesCropped[k] *= ratio;
1374  }
1375  }
1376  }
1377  }
1378 
1379 }
1380 
1381 void HMWSoln::counterIJ_setup(void) const
1382 {
1383  m_CounterIJ.resize(m_kk * m_kk);
1384  int counter = 0;
1385  for (size_t i = 0; i < m_kk; i++) {
1386  size_t n = i;
1387  size_t nc = m_kk * i;
1388  m_CounterIJ[n] = 0;
1389  m_CounterIJ[nc] = 0;
1390  }
1391  for (size_t i = 1; i < (m_kk - 1); i++) {
1392  size_t n = m_kk * i + i;
1393  m_CounterIJ[n] = 0;
1394  for (size_t j = (i+1); j < m_kk; j++) {
1395  n = m_kk * j + i;
1396  size_t nc = m_kk * i + j;
1397  counter++;
1398  m_CounterIJ[n] = counter;
1399  m_CounterIJ[nc] = counter;
1400  }
1401  }
1402 }
1403 
1404 void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
1405 {
1406  double T = temperature();
1407  const double twoT = 2.0 * T;
1408  const double invT = 1.0 / T;
1409  const double invT2 = invT * invT;
1410  const double twoinvT3 = 2.0 * invT * invT2;
1411  double Tr = m_TempPitzerRef;
1412  double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1413  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1414  tlin = T - Tr;
1415  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1416  tlin = T - Tr;
1417  tquad = T * T - Tr * Tr;
1418  tln = log(T/ Tr);
1419  tinv = 1.0/T - 1.0/Tr;
1420  }
1421 
1422  for (size_t i = 1; i < (m_kk - 1); i++) {
1423  for (size_t j = (i+1); j < m_kk; j++) {
1424 
1425  /*
1426  * Find the counterIJ for the symmetric binary interaction
1427  */
1428  size_t n = m_kk*i + j;
1429  size_t counterIJ = m_CounterIJ[n];
1430 
1431  const double* beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ);
1432  const double* beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ);
1433  const double* beta2MX_coeff = m_Beta2MX_ij_coeff.ptrColumn(counterIJ);
1434  const double* CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ);
1435  const double* Theta_coeff = m_Theta_ij_coeff.ptrColumn(counterIJ);
1436 
1437  switch (m_formPitzerTemp) {
1438  case PITZER_TEMP_CONSTANT:
1439  break;
1440  case PITZER_TEMP_LINEAR:
1441 
1442  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1443  + beta0MX_coeff[1]*tlin;
1444  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1];
1445  m_Beta0MX_ij_LL[counterIJ] = 0.0;
1446 
1447  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1448  + beta1MX_coeff[1]*tlin;
1449  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1];
1450  m_Beta1MX_ij_LL[counterIJ] = 0.0;
1451 
1452  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1453  + beta2MX_coeff[1]*tlin;
1454  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1];
1455  m_Beta2MX_ij_LL[counterIJ] = 0.0;
1456 
1457  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1458  + CphiMX_coeff[1]*tlin;
1459  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1];
1460  m_CphiMX_ij_LL[counterIJ] = 0.0;
1461 
1462  m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1463  m_Theta_ij_L[counterIJ] = Theta_coeff[1];
1464  m_Theta_ij_LL[counterIJ] = 0.0;
1465 
1466  break;
1467 
1468  case PITZER_TEMP_COMPLEX1:
1469  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1470  + beta0MX_coeff[1]*tlin
1471  + beta0MX_coeff[2]*tquad
1472  + beta0MX_coeff[3]*tinv
1473  + beta0MX_coeff[4]*tln;
1474 
1475  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1476  + beta1MX_coeff[1]*tlin
1477  + beta1MX_coeff[2]*tquad
1478  + beta1MX_coeff[3]*tinv
1479  + beta1MX_coeff[4]*tln;
1480 
1481  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1482  + beta2MX_coeff[1]*tlin
1483  + beta2MX_coeff[2]*tquad
1484  + beta2MX_coeff[3]*tinv
1485  + beta2MX_coeff[4]*tln;
1486 
1487  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1488  + CphiMX_coeff[1]*tlin
1489  + CphiMX_coeff[2]*tquad
1490  + CphiMX_coeff[3]*tinv
1491  + CphiMX_coeff[4]*tln;
1492 
1493  m_Theta_ij[counterIJ] = Theta_coeff[0]
1494  + Theta_coeff[1]*tlin
1495  + Theta_coeff[2]*tquad
1496  + Theta_coeff[3]*tinv
1497  + Theta_coeff[4]*tln;
1498 
1499  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]
1500  + beta0MX_coeff[2]*twoT
1501  - beta0MX_coeff[3]*invT2
1502  + beta0MX_coeff[4]*invT;
1503 
1504  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]
1505  + beta1MX_coeff[2]*twoT
1506  - beta1MX_coeff[3]*invT2
1507  + beta1MX_coeff[4]*invT;
1508 
1509  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1]
1510  + beta2MX_coeff[2]*twoT
1511  - beta2MX_coeff[3]*invT2
1512  + beta2MX_coeff[4]*invT;
1513 
1514  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]
1515  + CphiMX_coeff[2]*twoT
1516  - CphiMX_coeff[3]*invT2
1517  + CphiMX_coeff[4]*invT;
1518 
1519  m_Theta_ij_L[counterIJ] = Theta_coeff[1]
1520  + Theta_coeff[2]*twoT
1521  - Theta_coeff[3]*invT2
1522  + Theta_coeff[4]*invT;
1523 
1524  doDerivs = 2;
1525  if (doDerivs > 1) {
1526  m_Beta0MX_ij_LL[counterIJ] =
1527  + beta0MX_coeff[2]*2.0
1528  + beta0MX_coeff[3]*twoinvT3
1529  - beta0MX_coeff[4]*invT2;
1530 
1531  m_Beta1MX_ij_LL[counterIJ] =
1532  + beta1MX_coeff[2]*2.0
1533  + beta1MX_coeff[3]*twoinvT3
1534  - beta1MX_coeff[4]*invT2;
1535 
1536  m_Beta2MX_ij_LL[counterIJ] =
1537  + beta2MX_coeff[2]*2.0
1538  + beta2MX_coeff[3]*twoinvT3
1539  - beta2MX_coeff[4]*invT2;
1540 
1541  m_CphiMX_ij_LL[counterIJ] =
1542  + CphiMX_coeff[2]*2.0
1543  + CphiMX_coeff[3]*twoinvT3
1544  - CphiMX_coeff[4]*invT2;
1545 
1546  m_Theta_ij_LL[counterIJ] =
1547  + Theta_coeff[2]*2.0
1548  + Theta_coeff[3]*twoinvT3
1549  - Theta_coeff[4]*invT2;
1550  }
1551  break;
1552  }
1553  }
1554  }
1555 
1556  // Lambda interactions and Mu_nnn
1557  // i must be neutral for this term to be nonzero. We take advantage of this
1558  // here to lower the operation count.
1559  for (size_t i = 1; i < m_kk; i++) {
1560  if (charge(i) == 0.0) {
1561  for (size_t j = 1; j < m_kk; j++) {
1562  size_t n = i * m_kk + j;
1563  const double* Lambda_coeff = m_Lambda_nj_coeff.ptrColumn(n);
1564  switch (m_formPitzerTemp) {
1565  case PITZER_TEMP_CONSTANT:
1566  m_Lambda_nj(i,j) = Lambda_coeff[0];
1567  break;
1568  case PITZER_TEMP_LINEAR:
1569  m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1570  m_Lambda_nj_L(i,j) = Lambda_coeff[1];
1571  m_Lambda_nj_LL(i,j) = 0.0;
1572  break;
1573  case PITZER_TEMP_COMPLEX1:
1574  m_Lambda_nj(i,j) = Lambda_coeff[0]
1575  + Lambda_coeff[1]*tlin
1576  + Lambda_coeff[2]*tquad
1577  + Lambda_coeff[3]*tinv
1578  + Lambda_coeff[4]*tln;
1579 
1580  m_Lambda_nj_L(i,j) = Lambda_coeff[1]
1581  + Lambda_coeff[2]*twoT
1582  - Lambda_coeff[3]*invT2
1583  + Lambda_coeff[4]*invT;
1584 
1585  m_Lambda_nj_LL(i,j) =
1586  Lambda_coeff[2]*2.0
1587  + Lambda_coeff[3]*twoinvT3
1588  - Lambda_coeff[4]*invT2;
1589  }
1590 
1591  if (j == i) {
1592  const double* Mu_coeff = m_Mu_nnn_coeff.ptrColumn(i);
1593  switch (m_formPitzerTemp) {
1594  case PITZER_TEMP_CONSTANT:
1595  m_Mu_nnn[i] = Mu_coeff[0];
1596  break;
1597  case PITZER_TEMP_LINEAR:
1598  m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1599  m_Mu_nnn_L[i] = Mu_coeff[1];
1600  m_Mu_nnn_LL[i] = 0.0;
1601  break;
1602  case PITZER_TEMP_COMPLEX1:
1603  m_Mu_nnn[i] = Mu_coeff[0]
1604  + Mu_coeff[1]*tlin
1605  + Mu_coeff[2]*tquad
1606  + Mu_coeff[3]*tinv
1607  + Mu_coeff[4]*tln;
1608 
1609  m_Mu_nnn_L[i] = Mu_coeff[1]
1610  + Mu_coeff[2]*twoT
1611  - Mu_coeff[3]*invT2
1612  + Mu_coeff[4]*invT;
1613 
1614  m_Mu_nnn_LL[i] =
1615  Mu_coeff[2]*2.0
1616  + Mu_coeff[3]*twoinvT3
1617  - Mu_coeff[4]*invT2;
1618  }
1619  }
1620  }
1621  }
1622  }
1623 
1624  switch(m_formPitzerTemp) {
1625  case PITZER_TEMP_CONSTANT:
1626  for (size_t i = 1; i < m_kk; i++) {
1627  for (size_t j = 1; j < m_kk; j++) {
1628  for (size_t k = 1; k < m_kk; k++) {
1629  size_t n = i * m_kk *m_kk + j * m_kk + k ;
1630  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1631  m_Psi_ijk[n] = Psi_coeff[0];
1632  }
1633  }
1634  }
1635  break;
1636  case PITZER_TEMP_LINEAR:
1637  for (size_t i = 1; i < m_kk; i++) {
1638  for (size_t j = 1; j < m_kk; j++) {
1639  for (size_t k = 1; k < m_kk; k++) {
1640  size_t n = i * m_kk *m_kk + j * m_kk + k ;
1641  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1642  m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1643  m_Psi_ijk_L[n] = Psi_coeff[1];
1644  m_Psi_ijk_LL[n] = 0.0;
1645  }
1646  }
1647  }
1648  break;
1649  case PITZER_TEMP_COMPLEX1:
1650  for (size_t i = 1; i < m_kk; i++) {
1651  for (size_t j = 1; j < m_kk; j++) {
1652  for (size_t k = 1; k < m_kk; k++) {
1653  size_t n = i * m_kk *m_kk + j * m_kk + k ;
1654  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1655  m_Psi_ijk[n] = Psi_coeff[0]
1656  + Psi_coeff[1]*tlin
1657  + Psi_coeff[2]*tquad
1658  + Psi_coeff[3]*tinv
1659  + Psi_coeff[4]*tln;
1660 
1661  m_Psi_ijk_L[n] = Psi_coeff[1]
1662  + Psi_coeff[2]*twoT
1663  - Psi_coeff[3]*invT2
1664  + Psi_coeff[4]*invT;
1665 
1666  m_Psi_ijk_LL[n] =
1667  Psi_coeff[2]*2.0
1668  + Psi_coeff[3]*twoinvT3
1669  - Psi_coeff[4]*invT2;
1670  }
1671  }
1672  }
1673  break;
1674  }
1675 
1676 }
1677 
1678 void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
1679 {
1680  /*
1681  * HKM -> Assumption is made that the solvent is
1682  * species 0.
1683  */
1684  if (m_indexSolvent != 0) {
1685  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
1686  "Wrong index solvent value!");
1687  }
1688 
1689  /*
1690  * Use the CROPPED molality of the species in solution.
1691  */
1692  const double* molality = DATA_PTR(m_molalitiesCropped);
1693 
1694  /*
1695  * These are data inputs about the Pitzer correlation. They come
1696  * from the input file for the Pitzer model.
1697  */
1698  const double* beta0MX = DATA_PTR(m_Beta0MX_ij);
1699  const double* beta1MX = DATA_PTR(m_Beta1MX_ij);
1700  const double* beta2MX = DATA_PTR(m_Beta2MX_ij);
1701  const double* CphiMX = DATA_PTR(m_CphiMX_ij);
1702  const double* thetaij = DATA_PTR(m_Theta_ij);
1703  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
1704  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
1705 
1706  const double* psi_ijk = DATA_PTR(m_Psi_ijk);
1707  //n = k + j * m_kk + i * m_kk * m_kk;
1708 
1709 
1710  double* gamma_Unscaled = DATA_PTR(m_gamma_tmp);
1711  /*
1712  * Local variables defined by Coltrin
1713  */
1714  double etheta[5][5], etheta_prime[5][5], sqrtIs;
1715  /*
1716  * Molality based ionic strength of the solution
1717  */
1718  double Is = 0.0;
1719  /*
1720  * Molarcharge of the solution: In Pitzer's notation,
1721  * this is his variable called "Z".
1722  */
1723  double molarcharge = 0.0;
1724  /*
1725  * molalitysum is the sum of the molalities over all solutes,
1726  * even those with zero charge.
1727  */
1728  double molalitysumUncropped = 0.0;
1729 
1730  double* gfunc = DATA_PTR(m_gfunc_IJ);
1731  double* g2func = DATA_PTR(m_g2func_IJ);
1732  double* hfunc = DATA_PTR(m_hfunc_IJ);
1733  double* h2func = DATA_PTR(m_h2func_IJ);
1734  double* BMX = DATA_PTR(m_BMX_IJ);
1735  double* BprimeMX = DATA_PTR(m_BprimeMX_IJ);
1736  double* BphiMX = DATA_PTR(m_BphiMX_IJ);
1737  double* Phi = DATA_PTR(m_Phi_IJ);
1738  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
1739  double* Phiphi = DATA_PTR(m_PhiPhi_IJ);
1740  double* CMX = DATA_PTR(m_CMX_IJ);
1741 
1742  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1743  printf("\n Debugging information from hmw_act \n");
1744  }
1745  /*
1746  * Make sure the counter variables are setup
1747  */
1748  counterIJ_setup();
1749 
1750  /*
1751  * ---------- Calculate common sums over solutes ---------------------
1752  */
1753  for (size_t n = 1; n < m_kk; n++) {
1754  // ionic strength
1755  Is += charge(n) * charge(n) * molality[n];
1756  // total molar charge
1757  molarcharge += fabs(charge(n)) * molality[n];
1758  molalitysumUncropped += m_molalities[n];
1759  }
1760  Is *= 0.5;
1761 
1762  /*
1763  * Store the ionic molality in the object for reference.
1764  */
1765  m_IionicMolality = Is;
1766  sqrtIs = sqrt(Is);
1767  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1768  printf(" Step 1: \n");
1769  printf(" ionic strenth = %14.7le \n total molar "
1770  "charge = %14.7le \n", Is, molarcharge);
1771  }
1772 
1773  /*
1774  * The following call to calc_lambdas() calculates all 16 elements
1775  * of the elambda and elambda1 arrays, given the value of the
1776  * ionic strength (Is)
1777  */
1778  calc_lambdas(Is);
1779 
1780  /*
1781  * ----- Step 2: Find the coefficients E-theta and -------------------
1782  * E-thetaprime for all combinations of positive
1783  * unlike charges up to 4
1784  */
1785  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1786  printf(" Step 2: \n");
1787  }
1788  for (int z1 = 1; z1 <=4; z1++) {
1789  for (int z2 =1; z2 <=4; z2++) {
1790  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
1791  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1792  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
1793  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
1794  }
1795  }
1796  }
1797 
1798  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1799  printf(" Step 3: \n");
1800  printf(" Species Species g(x) "
1801  " hfunc(x) \n");
1802  }
1803 
1804  /*
1805  * calculate g(x) and hfunc(x) for each cation-anion pair MX
1806  * In the original literature, hfunc, was called gprime. However,
1807  * it's not the derivative of g(x), so I renamed it.
1808  */
1809  for (size_t i = 1; i < (m_kk - 1); i++) {
1810  for (size_t j = (i+1); j < m_kk; j++) {
1811  /*
1812  * Find the counterIJ for the symmetric binary interaction
1813  */
1814  size_t n = m_kk*i + j;
1815  size_t counterIJ = m_CounterIJ[n];
1816 
1817  /*
1818  * Only loop over oppositely charge species
1819  */
1820  if (charge(i)*charge(j) < 0) {
1821  /*
1822  * x is a reduced function variable
1823  */
1824  double x1 = sqrtIs * alpha1MX[counterIJ];
1825  if (x1 > 1.0E-100) {
1826  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
1827  hfunc[counterIJ] = -2.0 *
1828  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
1829  } else {
1830  gfunc[counterIJ] = 0.0;
1831  hfunc[counterIJ] = 0.0;
1832  }
1833 
1834  if (beta2MX[counterIJ] != 0.0) {
1835  double x2 = sqrtIs * alpha2MX[counterIJ];
1836  if (x2 > 1.0E-100) {
1837  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
1838  h2func[counterIJ] = -2.0 *
1839  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
1840  } else {
1841  g2func[counterIJ] = 0.0;
1842  h2func[counterIJ] = 0.0;
1843  }
1844  }
1845  } else {
1846  gfunc[counterIJ] = 0.0;
1847  hfunc[counterIJ] = 0.0;
1848  }
1849  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1850  std::string sni = speciesName(i);
1851  std::string snj = speciesName(j);
1852  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
1853  gfunc[counterIJ], hfunc[counterIJ]);
1854  }
1855  }
1856  }
1857 
1858  /*
1859  * --------- SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX ----------
1860  * --------- Agrees with Pitzer, Eq. (49), (51), (55)
1861  */
1862  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1863  printf(" Step 4: \n");
1864  printf(" Species Species BMX "
1865  "BprimeMX BphiMX \n");
1866  }
1867 
1868  for (size_t i = 1; i < m_kk - 1; i++) {
1869  for (size_t j = i+1; j < m_kk; j++) {
1870  /*
1871  * Find the counterIJ for the symmetric binary interaction
1872  */
1873  size_t n = m_kk*i + j;
1874  size_t counterIJ = m_CounterIJ[n];
1875 
1876  /*
1877  * both species have a non-zero charge, and one is positive
1878  * and the other is negative
1879  */
1880  if (charge(i)*charge(j) < 0.0) {
1881  BMX[counterIJ] = beta0MX[counterIJ]
1882  + beta1MX[counterIJ] * gfunc[counterIJ]
1883  + beta2MX[counterIJ] * g2func[counterIJ];
1884 
1885  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1886  printf("%d %g: %g %g %g %g\n",
1887  (int) counterIJ, BMX[counterIJ], beta0MX[counterIJ],
1888  beta1MX[counterIJ], beta2MX[counterIJ], gfunc[counterIJ]);
1889  }
1890  if (Is > 1.0E-150) {
1891  BprimeMX[counterIJ] = (beta1MX[counterIJ] * hfunc[counterIJ]/Is +
1892  beta2MX[counterIJ] * h2func[counterIJ]/Is);
1893  } else {
1894  BprimeMX[counterIJ] = 0.0;
1895  }
1896  BphiMX[counterIJ] = BMX[counterIJ] + Is*BprimeMX[counterIJ];
1897  } else {
1898  BMX[counterIJ] = 0.0;
1899  BprimeMX[counterIJ] = 0.0;
1900  BphiMX[counterIJ] = 0.0;
1901  }
1902  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1903  std::string sni = speciesName(i);
1904  std::string snj = speciesName(j);
1905  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
1906  sni.c_str(), snj.c_str(),
1907  BMX[counterIJ], BprimeMX[counterIJ], BphiMX[counterIJ]);
1908  }
1909  }
1910  }
1911 
1912  /*
1913  * --------- SUBSECTION TO CALCULATE CMX ----------
1914  * --------- Agrees with Pitzer, Eq. (53).
1915  */
1916  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1917  printf(" Step 5: \n");
1918  printf(" Species Species CMX \n");
1919  }
1920  for (size_t i = 1; i < m_kk-1; i++) {
1921  for (size_t j = i+1; j < m_kk; j++) {
1922  /*
1923  * Find the counterIJ for the symmetric binary interaction
1924  */
1925  size_t n = m_kk*i + j;
1926  size_t counterIJ = m_CounterIJ[n];
1927  /*
1928  * both species have a non-zero charge, and one is positive
1929  * and the other is negative
1930  */
1931  if (charge(i)*charge(j) < 0.0) {
1932  CMX[counterIJ] = CphiMX[counterIJ]/
1933  (2.0* sqrt(fabs(charge(i)*charge(j))));
1934  } else {
1935  CMX[counterIJ] = 0.0;
1936  }
1937  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1938  std::string sni = speciesName(i);
1939  std::string snj = speciesName(j);
1940  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
1941  CMX[counterIJ]);
1942  }
1943  }
1944  }
1945 
1946  /*
1947  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
1948  * --------- Agrees with Pitzer, Eq. 72, 73, 74
1949  */
1950  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1951  printf(" Step 6: \n");
1952  printf(" Species Species Phi_ij "
1953  " Phiprime_ij Phi^phi_ij \n");
1954  }
1955  for (size_t i = 1; i < m_kk-1; i++) {
1956  for (size_t j = i+1; j < m_kk; j++) {
1957  /*
1958  * Find the counterIJ for the symmetric binary interaction
1959  */
1960  size_t n = m_kk*i + j;
1961  size_t counterIJ = m_CounterIJ[n];
1962  /*
1963  * both species have a non-zero charge, and one is positive
1964  * and the other is negative
1965  */
1966  if (charge(i)*charge(j) > 0) {
1967  int z1 = (int) fabs(charge(i));
1968  int z2 = (int) fabs(charge(j));
1969  Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2];
1970  Phiprime[counterIJ] = etheta_prime[z1][z2];
1971  Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ];
1972  } else {
1973  Phi[counterIJ] = 0.0;
1974  Phiprime[counterIJ] = 0.0;
1975  Phiphi[counterIJ] = 0.0;
1976  }
1977  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1978  std::string sni = speciesName(i);
1979  std::string snj = speciesName(j);
1980  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
1981  sni.c_str(), snj.c_str(),
1982  Phi[counterIJ], Phiprime[counterIJ], Phiphi[counterIJ]);
1983  }
1984  }
1985  }
1986 
1987  /*
1988  * ------------- SUBSECTION FOR CALCULATION OF F ----------------------
1989  * ------------ Agrees with Pitzer Eqn. (65) --------------------------
1990  */
1991  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1992  printf(" Step 7: \n");
1993  }
1994  double Aphi = A_Debye_TP() / 3.0;
1995  double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
1996  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
1997  if (DEBUG_MODE_ENABLED && m_debugCalc) {
1998  printf(" initial value of F = %10.6f \n", F);
1999  }
2000  for (size_t i = 1; i < m_kk-1; i++) {
2001  for (size_t j = i+1; j < m_kk; j++) {
2002  /*
2003  * Find the counterIJ for the symmetric binary interaction
2004  */
2005  size_t n = m_kk*i + j;
2006  size_t counterIJ = m_CounterIJ[n];
2007  /*
2008  * both species have a non-zero charge, and one is positive
2009  * and the other is negative
2010  */
2011  if (charge(i)*charge(j) < 0) {
2012  F = F + molality[i]*molality[j] * BprimeMX[counterIJ];
2013  }
2014  /*
2015  * Both species have a non-zero charge, and they
2016  * have the same sign
2017  */
2018  if (charge(i)*charge(j) > 0) {
2019  F = F + molality[i]*molality[j] * Phiprime[counterIJ];
2020  }
2021  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2022  printf(" F = %10.6f \n", F);
2023  }
2024  }
2025  }
2026  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2027  printf(" Step 8: Summing in All Contributions to Activity Coefficients \n");
2028  }
2029 
2030  for (size_t i = 1; i < m_kk; i++) {
2031 
2032  /*
2033  * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS -----
2034  * -------- -> equations agree with my notes, Eqn. (118).
2035  * -> Equations agree with Pitzer, eqn.(63)
2036  */
2037  if (charge(i) > 0.0) {
2038  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2039  std::string sni = speciesName(i);
2040  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2041  }
2042  // species i is the cation (positive) to calc the actcoeff
2043  double zsqF = charge(i)*charge(i)*F;
2044  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2045  printf(" Unary term: z*z*F = %10.5f\n", zsqF);
2046  }
2047  double sum1 = 0.0;
2048  double sum2 = 0.0;
2049  double sum3 = 0.0;
2050  double sum4 = 0.0;
2051  double sum5 = 0.0;
2052  for (size_t j = 1; j < m_kk; j++) {
2053  /*
2054  * Find the counterIJ for the symmetric binary interaction
2055  */
2056  size_t n = m_kk*i + j;
2057  size_t counterIJ = m_CounterIJ[n];
2058 
2059  if (charge(j) < 0.0) {
2060  // sum over all anions
2061  sum1 = sum1 + molality[j]*
2062  (2.0*BMX[counterIJ] + molarcharge*CMX[counterIJ]);
2063  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2064  std::string snj = speciesName(j) + ":";
2065  printf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2066  molality[j]*2.0*BMX[counterIJ]);
2067  printf(" m_j Z CMX = %10.5f\n",
2068  molality[j]* molarcharge*CMX[counterIJ]);
2069  }
2070  if (j < m_kk-1) {
2071  /*
2072  * This term is the ternary interaction involving the
2073  * non-duplicate sum over double anions, j, k, with
2074  * respect to the cation, i.
2075  */
2076  for (size_t k = j+1; k < m_kk; k++) {
2077  // an inner sum over all anions
2078  if (charge(k) < 0.0) {
2079  n = k + j * m_kk + i * m_kk * m_kk;
2080  sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2081  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2082  if (psi_ijk[n] != 0.0) {
2083  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2084  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2085  molality[j]*molality[k]*psi_ijk[n]);
2086  }
2087  }
2088  }
2089  }
2090  }
2091  }
2092 
2093 
2094  if (charge(j) > 0.0) {
2095  // sum over all cations
2096  if (j != i) {
2097  sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2098  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2099  if ((molality[j] * Phi[counterIJ])!= 0.0) {
2100  std::string snj = speciesName(j) + ":";
2101  printf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj.c_str(),
2102  molality[j]*(2.0*Phi[counterIJ]));
2103  }
2104  }
2105  }
2106  for (size_t k = 1; k < m_kk; k++) {
2107  if (charge(k) < 0.0) {
2108  // two inner sums over anions
2109 
2110  n = k + j * m_kk + i * m_kk * m_kk;
2111  sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2112  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2113  if (psi_ijk[n] != 0.0) {
2114  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2115  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2116  molality[j]*molality[k]*psi_ijk[n]);
2117  }
2118  }
2119  /*
2120  * Find the counterIJ for the j,k interaction
2121  */
2122  n = m_kk*j + k;
2123  size_t counterIJ2 = m_CounterIJ[n];
2124  sum4 = sum4 + (fabs(charge(i))*
2125  molality[j]*molality[k]*CMX[counterIJ2]);
2126  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2127  if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2128  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2129  printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2130  fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]);
2131  }
2132  }
2133  }
2134  }
2135  }
2136 
2137  /*
2138  * Handle neutral j species
2139  */
2140  if (charge(j) == 0) {
2141  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i);
2142  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2143  if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
2144  std::string snj = speciesName(j) + ":";
2145  printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2146  molality[j]*2.0*m_Lambda_nj(j,i));
2147  }
2148  }
2149  /*
2150  * Zeta interaction term
2151  */
2152  for (size_t k = 1; k < m_kk; k++) {
2153  if (charge(k) < 0.0) {
2154  size_t izeta = j;
2155  size_t jzeta = i;
2156  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2157  double zeta = psi_ijk[n];
2158  if (zeta != 0.0) {
2159  sum5 = sum5 + molality[j]*molality[k]*zeta;
2160  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2161  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2162  printf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj.c_str(),
2163  molality[j]*molality[k]*psi_ijk[n]);
2164  }
2165  }
2166  }
2167  }
2168  }
2169  }
2170  /*
2171  * Add all of the contributions up to yield the log of the
2172  * solute activity coefficients (molality scale)
2173  */
2174  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2175  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2176  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2177  std::string sni = speciesName(i);
2178  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
2179  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
2180  }
2181  }
2182 
2183  /*
2184  * -------- SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS ------
2185  * -------- -> equations agree with my notes, Eqn. (119).
2186  * -> Equations agree with Pitzer, eqn.(64)
2187  */
2188  if (charge(i) < 0) {
2189  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2190  std::string sni = speciesName(i);
2191  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2192  }
2193  // species i is an anion (negative)
2194  double zsqF = charge(i)*charge(i)*F;
2195  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2196  printf(" Unary term: z*z*F = %10.5f\n", zsqF);
2197  }
2198  double sum1 = 0.0;
2199  double sum2 = 0.0;
2200  double sum3 = 0.0;
2201  double sum4 = 0.0;
2202  double sum5 = 0.0;
2203  for (size_t j = 1; j < m_kk; j++) {
2204  /*
2205  * Find the counterIJ for the symmetric binary interaction
2206  */
2207  size_t n = m_kk*i + j;
2208  size_t counterIJ = m_CounterIJ[n];
2209 
2210  /*
2211  * For Anions, do the cation interactions.
2212  */
2213  if (charge(j) > 0) {
2214  sum1 = sum1 + molality[j]*
2215  (2.0*BMX[counterIJ]+molarcharge*CMX[counterIJ]);
2216  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2217  std::string snj = speciesName(j) + ":";
2218  printf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj.c_str(),
2219  molality[j]*2.0*BMX[counterIJ]);
2220  printf(" m_j Z CMX = %10.5f\n",
2221  molality[j]* molarcharge*CMX[counterIJ]);
2222  }
2223  if (j < m_kk-1) {
2224  for (size_t k = j+1; k < m_kk; k++) {
2225  // an inner sum over all cations
2226  if (charge(k) > 0) {
2227  n = k + j * m_kk + i * m_kk * m_kk;
2228  sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2229  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2230  if (psi_ijk[n] != 0.0) {
2231  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2232  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2233  molality[j]*molality[k]*psi_ijk[n]);
2234  }
2235  }
2236  }
2237  }
2238  }
2239  }
2240 
2241  /*
2242  * For Anions, do the other anion interactions.
2243  */
2244  if (charge(j) < 0.0) {
2245  // sum over all anions
2246  if (j != i) {
2247  sum2 = sum2 + molality[j]*(2.0*Phi[counterIJ]);
2248  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2249  if ((molality[j] * Phi[counterIJ])!= 0.0) {
2250  std::string snj = speciesName(j) + ":";
2251  printf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj.c_str(),
2252  molality[j]*(2.0*Phi[counterIJ]));
2253  }
2254  }
2255  }
2256  for (size_t k = 1; k < m_kk; k++) {
2257  if (charge(k) > 0.0) {
2258  // two inner sums over cations
2259  n = k + j * m_kk + i * m_kk * m_kk;
2260  sum2 = sum2 + molality[j]*molality[k]*psi_ijk[n];
2261  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2262  if (psi_ijk[n] != 0.0) {
2263  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2264  printf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2265  molality[j]*molality[k]*psi_ijk[n]);
2266  }
2267  }
2268  /*
2269  * Find the counterIJ for the symmetric binary interaction
2270  */
2271  n = m_kk*j + k;
2272  size_t counterIJ2 = m_CounterIJ[n];
2273  sum4 = sum4 +
2274  (fabs(charge(i))*
2275  molality[j]*molality[k]*CMX[counterIJ2]);
2276  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2277  if ((molality[j]*molality[k]*CMX[counterIJ2]) != 0.0) {
2278  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2279  printf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj.c_str(),
2280  fabs(charge(i))* molality[j]*molality[k]*CMX[counterIJ2]);
2281  }
2282  }
2283  }
2284  }
2285  }
2286 
2287  /*
2288  * for Anions, do the neutral species interaction
2289  */
2290  if (charge(j) == 0.0) {
2291  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj(j,i);
2292  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2293  if ((molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
2294  std::string snj = speciesName(j) + ":";
2295  printf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj.c_str(),
2296  molality[j]*2.0*m_Lambda_nj(j,i));
2297  }
2298  }
2299  /*
2300  * Zeta interaction term
2301  */
2302  for (size_t k = 1; k < m_kk; k++) {
2303  if (charge(k) > 0.0) {
2304  size_t izeta = j;
2305  size_t jzeta = k;
2306  size_t kzeta = i;
2307  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2308  double zeta = psi_ijk[n];
2309  if (zeta != 0.0) {
2310  sum5 = sum5 + molality[j]*molality[k]*zeta;
2311  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2312  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2313  printf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj.c_str(),
2314  molality[j]*molality[k]*psi_ijk[n]);
2315  }
2316  }
2317  }
2318  }
2319  }
2320  }
2321  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
2322  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2323  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2324  std::string sni = speciesName(i);
2325  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
2326  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
2327  }
2328  }
2329  /*
2330  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
2331  * ------ -> equations agree with my notes,
2332  * -> Equations agree with Pitzer,
2333  */
2334  if (charge(i) == 0.0) {
2335  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2336  std::string sni = speciesName(i);
2337  printf(" Contributions to ln(ActCoeff_%s):\n", sni.c_str());
2338  }
2339  double sum1 = 0.0;
2340  double sum3 = 0.0;
2341  for (size_t j = 1; j < m_kk; j++) {
2342  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj(i,j);
2343  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2344  if (m_Lambda_nj(i,j) != 0.0) {
2345  std::string snj = speciesName(j) + ":";
2346  printf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj.c_str(),
2347  molality[j]*2.0*m_Lambda_nj(i,j));
2348  }
2349  }
2350  /*
2351  * Zeta term -> we piggyback on the psi term
2352  */
2353  if (charge(j) > 0.0) {
2354  for (size_t k = 1; k < m_kk; k++) {
2355  if (charge(k) < 0.0) {
2356  size_t n = k + j * m_kk + i * m_kk * m_kk;
2357  sum3 = sum3 + molality[j]*molality[k]*psi_ijk[n];
2358  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2359  if (psi_ijk[n] != 0.0) {
2360  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
2361  printf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj.c_str(),
2362  molality[j]*molality[k]*psi_ijk[n]);
2363  }
2364  }
2365  }
2366  }
2367  }
2368  }
2369  double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
2370  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2371  if (m_Mu_nnn[i] != 0.0) {
2372  printf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n",
2373  3.0 * molality[i]* molality[i] * m_Mu_nnn[i]);
2374  }
2375  }
2376  m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
2377  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
2378  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2379  std::string sni = speciesName(i);
2380  printf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
2381  sni.c_str(), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
2382  }
2383  }
2384 
2385  }
2386  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2387  printf(" Step 9: \n");
2388  }
2389  /*
2390  * -------- SUBSECTION FOR CALCULATING THE OSMOTIC COEFF ---------
2391  * -------- -> equations agree with my notes, Eqn. (117).
2392  * -> Equations agree with Pitzer, eqn.(62)
2393  */
2394  double sum1 = 0.0;
2395  double sum2 = 0.0;
2396  double sum3 = 0.0;
2397  double sum4 = 0.0;
2398  double sum5 = 0.0;
2399  double sum6 = 0.0;
2400  double sum7 = 0.0;
2401  /*
2402  * term1 is the DH term in the osmotic coefficient expression
2403  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
2404  * implementations.
2405  * Is = Ionic strength on the molality scale (units of (gmol/kg))
2406  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2407  */
2408  double term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2409 
2410  for (size_t j = 1; j < m_kk; j++) {
2411  /*
2412  * Loop Over Cations
2413  */
2414  if (charge(j) > 0.0) {
2415  for (size_t k = 1; k < m_kk; k++) {
2416  if (charge(k) < 0.0) {
2417  /*
2418  * Find the counterIJ for the symmetric j,k binary interaction
2419  */
2420  size_t n = m_kk*j + k;
2421  size_t counterIJ = m_CounterIJ[n];
2422 
2423  sum1 = sum1 + molality[j]*molality[k]*
2424  (BphiMX[counterIJ] + molarcharge*CMX[counterIJ]);
2425  }
2426  }
2427 
2428  for (size_t k = j+1; k < m_kk; k++) {
2429  if (j == (m_kk-1)) {
2430  // we should never reach this step
2431  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2432  "logic error 1 in Step 9 of hmw_act");
2433  }
2434  if (charge(k) > 0.0) {
2435  /*
2436  * Find the counterIJ for the symmetric j,k binary interaction
2437  * between 2 cations.
2438  */
2439  size_t n = m_kk*j + k;
2440  size_t counterIJ = m_CounterIJ[n];
2441  sum2 = sum2 + molality[j]*molality[k]*Phiphi[counterIJ];
2442  for (size_t m = 1; m < m_kk; m++) {
2443  if (charge(m) < 0.0) {
2444  // species m is an anion
2445  n = m + k * m_kk + j * m_kk * m_kk;
2446  sum2 = sum2 +
2447  molality[j]*molality[k]*molality[m]*psi_ijk[n];
2448  }
2449  }
2450  }
2451  }
2452  }
2453 
2454  /*
2455  * Loop Over Anions
2456  */
2457  if (charge(j) < 0) {
2458  for (size_t k = j+1; k < m_kk; k++) {
2459  if (j == m_kk-1) {
2460  // we should never reach this step
2461  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2462  "logic error 2 in Step 9 of hmw_act");
2463  }
2464  if (charge(k) < 0) {
2465  /*
2466  * Find the counterIJ for the symmetric j,k binary interaction
2467  * between two anions
2468  */
2469  size_t n = m_kk*j + k;
2470  size_t counterIJ = m_CounterIJ[n];
2471 
2472  sum3 = sum3 + molality[j]*molality[k]*Phiphi[counterIJ];
2473  for (size_t m = 1; m < m_kk; m++) {
2474  if (charge(m) > 0.0) {
2475  n = m + k * m_kk + j * m_kk * m_kk;
2476  sum3 = sum3 +
2477  molality[j]*molality[k]*molality[m]*psi_ijk[n];
2478  }
2479  }
2480  }
2481  }
2482  }
2483 
2484  /*
2485  * Loop Over Neutral Species
2486  */
2487  if (charge(j) == 0) {
2488  for (size_t k = 1; k < m_kk; k++) {
2489  if (charge(k) < 0.0) {
2490  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj(j,k);
2491  }
2492  if (charge(k) > 0.0) {
2493  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj(j,k);
2494  }
2495  if (charge(k) == 0.0) {
2496  if (k > j) {
2497  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj(j,k);
2498  } else if (k == j) {
2499  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj(j,k);
2500  }
2501  }
2502  if (charge(k) < 0.0) {
2503  size_t izeta = j;
2504  for (size_t m = 1; m < m_kk; m++) {
2505  if (charge(m) > 0.0) {
2506  size_t jzeta = m;
2507  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2508  double zeta = psi_ijk[n];
2509  if (zeta != 0.0) {
2510  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2511  }
2512  }
2513  }
2514  }
2515  }
2516  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j];
2517  }
2518  }
2519  double sum_m_phi_minus_1 = 2.0 *
2520  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2521  /*
2522  * Calculate the osmotic coefficient from
2523  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2524  */
2525  double osmotic_coef;
2526  if (molalitysumUncropped > 1.0E-150) {
2527  osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2528  } else {
2529  osmotic_coef = 1.0;
2530  }
2531  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2532  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
2533  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
2534  term1, sum1, sum2, sum3, sum4, sum5);
2535  printf(" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n",
2536  sum_m_phi_minus_1, osmotic_coef);
2537  printf(" Step 10: \n");
2538  }
2539  double lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2540 
2541  /*
2542  * In Cantera, we define the activity coefficient of the solvent as
2543  *
2544  * act_0 = actcoeff_0 * Xmol_0
2545  *
2546  * We have just computed act_0. However, this routine returns
2547  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2548  */
2549  double xmolSolvent = moleFraction(m_indexSolvent);
2550  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
2551  m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
2552  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2553  double wateract = exp(lnwateract);
2554  printf(" Weight of Solvent = %16.7g\n", m_weightSolvent);
2555  printf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
2556  printf(" ln_a_water=%10.6f a_water=%10.6f\n\n",
2557  lnwateract, wateract);
2558  }
2559 }
2560 
2561 void HMWSoln::s_update_dlnMolalityActCoeff_dT() const
2562 {
2563  static const int cacheId = m_cache.getId();
2564  CachedScalar cached = m_cache.getScalar(cacheId);
2565  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2566  return;
2567  }
2568 
2569  /*
2570  * Zero the unscaled 2nd derivatives
2571  */
2572  m_dlnActCoeffMolaldT_Unscaled.assign(m_kk, 0.0);
2573  /*
2574  * Do the actual calculation of the unscaled temperature derivatives
2575  */
2576  s_updatePitzer_dlnMolalityActCoeff_dT();
2577 
2578  for (size_t k = 1; k < m_kk; k++) {
2579  if (CROP_speciesCropped_[k] == 2) {
2580  m_dlnActCoeffMolaldT_Unscaled[k] = 0.0;
2581  }
2582  }
2583 
2584  if (CROP_speciesCropped_[0]) {
2585  m_dlnActCoeffMolaldT_Unscaled[0] = 0.0;
2586  }
2587 
2588  /*
2589  * Do the pH scaling to the derivatives
2590  */
2591  s_updateScaling_pHScaling_dT();
2592 
2593 
2594 }
2595 
2596 void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
2597 {
2598  /*
2599  * It may be assumed that the Pitzer activity coefficient routine is
2600  * called immediately preceding the calling of this routine. Therefore,
2601  * some quantities do not need to be recalculated in this routine.
2602  */
2603 
2604  /*
2605  * HKM -> Assumption is made that the solvent is
2606  * species 0.
2607  */
2608 #ifdef DEBUG_MODE
2609  m_debugCalc = 0;
2610 #endif
2611  if (m_indexSolvent != 0) {
2612  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2613  "Wrong index solvent value!");
2614  }
2615 
2616  const double* molality = DATA_PTR(m_molalitiesCropped);
2617  const double* beta0MX_L = DATA_PTR(m_Beta0MX_ij_L);
2618  const double* beta1MX_L = DATA_PTR(m_Beta1MX_ij_L);
2619  const double* beta2MX_L = DATA_PTR(m_Beta2MX_ij_L);
2620  const double* CphiMX_L = DATA_PTR(m_CphiMX_ij_L);
2621  const double* thetaij_L = DATA_PTR(m_Theta_ij_L);
2622  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
2623  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
2624  const double* psi_ijk_L = DATA_PTR(m_Psi_ijk_L);
2625  double* d_gamma_dT_Unscaled = DATA_PTR(m_gamma_tmp);
2626  /*
2627  * Local variables defined by Coltrin
2628  */
2629  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2630  /*
2631  * Molality based ionic strength of the solution
2632  */
2633  double Is = 0.0;
2634  /*
2635  * Molarcharge of the solution: In Pitzer's notation,
2636  * this is his variable called "Z".
2637  */
2638  double molarcharge = 0.0;
2639  /*
2640  * molalitysum is the sum of the molalities over all solutes,
2641  * even those with zero charge.
2642  */
2643  double molalitysum = 0.0;
2644 
2645  double* gfunc = DATA_PTR(m_gfunc_IJ);
2646  double* g2func = DATA_PTR(m_g2func_IJ);
2647  double* hfunc = DATA_PTR(m_hfunc_IJ);
2648  double* h2func = DATA_PTR(m_h2func_IJ);
2649  double* BMX_L = DATA_PTR(m_BMX_IJ_L);
2650  double* BprimeMX_L= DATA_PTR(m_BprimeMX_IJ_L);
2651  double* BphiMX_L = DATA_PTR(m_BphiMX_IJ_L);
2652  double* Phi_L = DATA_PTR(m_Phi_IJ_L);
2653  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
2654  double* Phiphi_L = DATA_PTR(m_PhiPhi_IJ_L);
2655  double* CMX_L = DATA_PTR(m_CMX_IJ_L);
2656 
2657  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2658  printf("\n Debugging information from "
2659  "s_Pitzer_dlnMolalityActCoeff_dT()\n");
2660  }
2661  /*
2662  * Make sure the counter variables are setup
2663  */
2664  counterIJ_setup();
2665 
2666  /*
2667  * ---------- Calculate common sums over solutes ---------------------
2668  */
2669  for (size_t n = 1; n < m_kk; n++) {
2670  // ionic strength
2671  Is += charge(n) * charge(n) * molality[n];
2672  // total molar charge
2673  molarcharge += fabs(charge(n)) * molality[n];
2674  molalitysum += molality[n];
2675  }
2676  Is *= 0.5;
2677 
2678  /*
2679  * Store the ionic molality in the object for reference.
2680  */
2681  m_IionicMolality = Is;
2682  sqrtIs = sqrt(Is);
2683  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2684  printf(" Step 1: \n");
2685  printf(" ionic strenth = %14.7le \n total molar "
2686  "charge = %14.7le \n", Is, molarcharge);
2687  }
2688 
2689  /*
2690  * The following call to calc_lambdas() calculates all 16 elements
2691  * of the elambda and elambda1 arrays, given the value of the
2692  * ionic strength (Is)
2693  */
2694  calc_lambdas(Is);
2695 
2696  /*
2697  * ----- Step 2: Find the coefficients E-theta and -------------------
2698  * E-thetaprime for all combinations of positive
2699  * unlike charges up to 4
2700  */
2701  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2702  printf(" Step 2: \n");
2703  }
2704  for (int z1 = 1; z1 <=4; z1++) {
2705  for (int z2 =1; z2 <=4; z2++) {
2706  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2707  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2708  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
2709  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
2710  }
2711  }
2712  }
2713 
2714  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2715  printf(" Step 3: \n");
2716  printf(" Species Species g(x) "
2717  " hfunc(x) \n");
2718  }
2719 
2720  /*
2721  * calculate g(x) and hfunc(x) for each cation-anion pair MX
2722  * In the original literature, hfunc, was called gprime. However,
2723  * it's not the derivative of g(x), so I renamed it.
2724  */
2725  for (size_t i = 1; i < (m_kk - 1); i++) {
2726  for (size_t j = (i+1); j < m_kk; j++) {
2727  /*
2728  * Find the counterIJ for the symmetric binary interaction
2729  */
2730  size_t n = m_kk*i + j;
2731  size_t counterIJ = m_CounterIJ[n];
2732  /*
2733  * Only loop over oppositely charge species
2734  */
2735  if (charge(i)*charge(j) < 0) {
2736  /*
2737  * x is a reduced function variable
2738  */
2739  double x1 = sqrtIs * alpha1MX[counterIJ];
2740  if (x1 > 1.0E-100) {
2741  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2742  hfunc[counterIJ] = -2.0 *
2743  (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
2744  } else {
2745  gfunc[counterIJ] = 0.0;
2746  hfunc[counterIJ] = 0.0;
2747  }
2748 
2749  if (beta2MX_L[counterIJ] != 0.0) {
2750  double x2 = sqrtIs * alpha2MX[counterIJ];
2751  if (x2 > 1.0E-100) {
2752  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2753  h2func[counterIJ] = -2.0 *
2754  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2755  } else {
2756  g2func[counterIJ] = 0.0;
2757  h2func[counterIJ] = 0.0;
2758  }
2759  }
2760  } else {
2761  gfunc[counterIJ] = 0.0;
2762  hfunc[counterIJ] = 0.0;
2763  }
2764  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2765  std::string sni = speciesName(i);
2766  std::string snj = speciesName(j);
2767  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
2768  gfunc[counterIJ], hfunc[counterIJ]);
2769  }
2770  }
2771  }
2772 
2773  /*
2774  * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L ----------
2775  * ------- These are now temperature derivatives of the
2776  * previously calculated quantities.
2777  */
2778  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2779  printf(" Step 4: \n");
2780  printf(" Species Species BMX "
2781  "BprimeMX BphiMX \n");
2782  }
2783 
2784  for (size_t i = 1; i < m_kk - 1; i++) {
2785  for (size_t j = i+1; j < m_kk; j++) {
2786  /*
2787  * Find the counterIJ for the symmetric binary interaction
2788  */
2789  size_t n = m_kk*i + j;
2790  size_t counterIJ = m_CounterIJ[n];
2791  /*
2792  * both species have a non-zero charge, and one is positive
2793  * and the other is negative
2794  */
2795  if (charge(i)*charge(j) < 0.0) {
2796  BMX_L[counterIJ] = beta0MX_L[counterIJ]
2797  + beta1MX_L[counterIJ] * gfunc[counterIJ]
2798  + beta2MX_L[counterIJ] * gfunc[counterIJ];
2799  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2800  printf("%d %g: %g %g %g %g\n",
2801  (int) counterIJ, BMX_L[counterIJ], beta0MX_L[counterIJ],
2802  beta1MX_L[counterIJ], beta2MX_L[counterIJ], gfunc[counterIJ]);
2803  }
2804  if (Is > 1.0E-150) {
2805  BprimeMX_L[counterIJ] = (beta1MX_L[counterIJ] * hfunc[counterIJ]/Is +
2806  beta2MX_L[counterIJ] * h2func[counterIJ]/Is);
2807  } else {
2808  BprimeMX_L[counterIJ] = 0.0;
2809  }
2810  BphiMX_L[counterIJ] = BMX_L[counterIJ] + Is*BprimeMX_L[counterIJ];
2811  } else {
2812  BMX_L[counterIJ] = 0.0;
2813  BprimeMX_L[counterIJ] = 0.0;
2814  BphiMX_L[counterIJ] = 0.0;
2815  }
2816  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2817  std::string sni = speciesName(i);
2818  std::string snj = speciesName(j);
2819  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
2820  sni.c_str(), snj.c_str(),
2821  BMX_L[counterIJ], BprimeMX_L[counterIJ], BphiMX_L[counterIJ]);
2822  }
2823  }
2824  }
2825 
2826  /*
2827  * --------- SUBSECTION TO CALCULATE CMX_L ----------
2828  * ---------
2829  */
2830  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2831  printf(" Step 5: \n");
2832  printf(" Species Species CMX \n");
2833  }
2834  for (size_t i = 1; i < m_kk-1; i++) {
2835  for (size_t j = i+1; j < m_kk; j++) {
2836  /*
2837  * Find the counterIJ for the symmetric binary interaction
2838  */
2839  size_t n = m_kk*i + j;
2840  size_t counterIJ = m_CounterIJ[n];
2841  /*
2842  * both species have a non-zero charge, and one is positive
2843  * and the other is negative
2844  */
2845  if (charge(i)*charge(j) < 0.0) {
2846  CMX_L[counterIJ] = CphiMX_L[counterIJ]/
2847  (2.0* sqrt(fabs(charge(i)*charge(j))));
2848  } else {
2849  CMX_L[counterIJ] = 0.0;
2850  }
2851  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2852  std::string sni = speciesName(i);
2853  std::string snj = speciesName(j);
2854  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
2855  CMX_L[counterIJ]);
2856  }
2857  }
2858  }
2859 
2860  /*
2861  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2862  * --------
2863  */
2864  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2865  printf(" Step 6: \n");
2866  printf(" Species Species Phi_ij "
2867  " Phiprime_ij Phi^phi_ij \n");
2868  }
2869  for (size_t i = 1; i < m_kk-1; i++) {
2870  for (size_t j = i+1; j < m_kk; j++) {
2871  /*
2872  * Find the counterIJ for the symmetric binary interaction
2873  */
2874  size_t n = m_kk*i + j;
2875  size_t counterIJ = m_CounterIJ[n];
2876  /*
2877  * both species have a non-zero charge, and one is positive
2878  * and the other is negative
2879  */
2880  if (charge(i)*charge(j) > 0) {
2881  Phi_L[counterIJ] = thetaij_L[counterIJ];
2882  Phiprime[counterIJ] = 0.0;
2883  Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ];
2884  } else {
2885  Phi_L[counterIJ] = 0.0;
2886  Phiprime[counterIJ] = 0.0;
2887  Phiphi_L[counterIJ] = 0.0;
2888  }
2889  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2890  std::string sni = speciesName(i);
2891  std::string snj = speciesName(j);
2892  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
2893  sni.c_str(), snj.c_str(),
2894  Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ]);
2895  }
2896  }
2897  }
2898 
2899  /*
2900  * ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
2901  */
2902  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2903  printf(" Step 7: \n");
2904  }
2905  double dA_DebyedT = dA_DebyedT_TP();
2906  double dAphidT = dA_DebyedT /3.0;
2907  double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2908  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2909  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2910  printf(" initial value of dFdT = %10.6f \n", dFdT);
2911  }
2912  for (size_t i = 1; i < m_kk-1; i++) {
2913  for (size_t j = i+1; j < m_kk; j++) {
2914  /*
2915  * Find the counterIJ for the symmetric binary interaction
2916  */
2917  size_t n = m_kk*i + j;
2918  size_t counterIJ = m_CounterIJ[n];
2919  /*
2920  * both species have a non-zero charge, and one is positive
2921  * and the other is negative
2922  */
2923  if (charge(i)*charge(j) < 0) {
2924  dFdT = dFdT + molality[i]*molality[j] * BprimeMX_L[counterIJ];
2925  }
2926  /*
2927  * Both species have a non-zero charge, and they
2928  * have the same sign, e.g., both positive or both negative.
2929  */
2930  if (charge(i)*charge(j) > 0) {
2931  dFdT = dFdT + molality[i]*molality[j] * Phiprime[counterIJ];
2932  }
2933  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2934  printf(" dFdT = %10.6f \n", dFdT);
2935  }
2936  }
2937  }
2938  if (DEBUG_MODE_ENABLED && m_debugCalc) {
2939  printf(" Step 8: \n");
2940  }
2941 
2942  for (size_t i = 1; i < m_kk; i++) {
2943  /*
2944  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
2945  * --
2946  */
2947  if (charge(i) > 0) {
2948  // species i is the cation (positive) to calc the actcoeff
2949  double zsqdFdT = charge(i)*charge(i)*dFdT;
2950  double sum1 = 0.0;
2951  double sum2 = 0.0;
2952  double sum3 = 0.0;
2953  double sum4 = 0.0;
2954  double sum5 = 0.0;
2955  for (size_t j = 1; j < m_kk; j++) {
2956  /*
2957  * Find the counterIJ for the symmetric binary interaction
2958  */
2959  size_t n = m_kk*i + j;
2960  size_t counterIJ = m_CounterIJ[n];
2961 
2962  if (charge(j) < 0.0) {
2963  // sum over all anions
2964  sum1 = sum1 + molality[j]*
2965  (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
2966  if (j < m_kk-1) {
2967  /*
2968  * This term is the ternary interaction involving the
2969  * non-duplicate sum over double anions, j, k, with
2970  * respect to the cation, i.
2971  */
2972  for (size_t k = j+1; k < m_kk; k++) {
2973  // an inner sum over all anions
2974  if (charge(k) < 0.0) {
2975  n = k + j * m_kk + i * m_kk * m_kk;
2976  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
2977  }
2978  }
2979  }
2980  }
2981 
2982 
2983  if (charge(j) > 0.0) {
2984  // sum over all cations
2985  if (j != i) {
2986  sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
2987  }
2988  for (size_t k = 1; k < m_kk; k++) {
2989  if (charge(k) < 0.0) {
2990  // two inner sums over anions
2991 
2992  n = k + j * m_kk + i * m_kk * m_kk;
2993  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
2994  /*
2995  * Find the counterIJ for the j,k interaction
2996  */
2997  n = m_kk*j + k;
2998  size_t counterIJ2 = m_CounterIJ[n];
2999  sum4 = sum4 + (fabs(charge(i))*
3000  molality[j]*molality[k]*CMX_L[counterIJ2]);
3001  }
3002  }
3003  }
3004 
3005  /*
3006  * Handle neutral j species
3007  */
3008  if (charge(j) == 0) {
3009  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i);
3010  }
3011  /*
3012  * Zeta interaction term
3013  */
3014  for (size_t k = 1; k < m_kk; k++) {
3015  if (charge(k) < 0.0) {
3016  size_t izeta = j;
3017  size_t jzeta = i;
3018  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3019  double zeta_L = psi_ijk_L[n];
3020  if (zeta_L != 0.0) {
3021  sum5 = sum5 + molality[j]*molality[k]*zeta_L;
3022  }
3023  }
3024  }
3025  }
3026  /*
3027  * Add all of the contributions up to yield the log of the
3028  * solute activity coefficients (molality scale)
3029  */
3030  m_dlnActCoeffMolaldT_Unscaled[i] =
3031  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
3032  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
3033  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3034  std::string sni = speciesName(i);
3035  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
3036  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
3037  printf(" %12g %12g %12g %12g %12g %12g\n",
3038  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
3039  }
3040  }
3041 
3042  /*
3043  * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------
3044  *
3045  */
3046  if (charge(i) < 0) {
3047  // species i is an anion (negative)
3048  double zsqdFdT = charge(i)*charge(i)*dFdT;
3049  double sum1 = 0.0;
3050  double sum2 = 0.0;
3051  double sum3 = 0.0;
3052  double sum4 = 0.0;
3053  double sum5 = 0.0;
3054  for (size_t j = 1; j < m_kk; j++) {
3055  /*
3056  * Find the counterIJ for the symmetric binary interaction
3057  */
3058  size_t n = m_kk*i + j;
3059  size_t counterIJ = m_CounterIJ[n];
3060 
3061  /*
3062  * For Anions, do the cation interactions.
3063  */
3064  if (charge(j) > 0) {
3065  sum1 = sum1 + molality[j]*
3066  (2.0*BMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3067  if (j < m_kk-1) {
3068  for (size_t k = j+1; k < m_kk; k++) {
3069  // an inner sum over all cations
3070  if (charge(k) > 0) {
3071  n = k + j * m_kk + i * m_kk * m_kk;
3072  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3073  }
3074  }
3075  }
3076  }
3077 
3078  /*
3079  * For Anions, do the other anion interactions.
3080  */
3081  if (charge(j) < 0.0) {
3082  // sum over all anions
3083  if (j != i) {
3084  sum2 = sum2 + molality[j]*(2.0*Phi_L[counterIJ]);
3085  }
3086  for (size_t k = 1; k < m_kk; k++) {
3087  if (charge(k) > 0.0) {
3088  // two inner sums over cations
3089  n = k + j * m_kk + i * m_kk * m_kk;
3090  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_L[n];
3091  /*
3092  * Find the counterIJ for the symmetric binary interaction
3093  */
3094  n = m_kk*j + k;
3095  size_t counterIJ2 = m_CounterIJ[n];
3096  sum4 = sum4 +
3097  (fabs(charge(i))*
3098  molality[j]*molality[k]*CMX_L[counterIJ2]);
3099  }
3100  }
3101  }
3102 
3103  /*
3104  * for Anions, do the neutral species interaction
3105  */
3106  if (charge(j) == 0.0) {
3107  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_L(j,i);
3108  for (size_t k = 1; k < m_kk; k++) {
3109  if (charge(k) > 0.0) {
3110  size_t izeta = j;
3111  size_t jzeta = k;
3112  size_t kzeta = i;
3113  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3114  double zeta_L = psi_ijk_L[n];
3115  if (zeta_L != 0.0) {
3116  sum5 = sum5 + molality[j]*molality[k]*zeta_L;
3117  }
3118  }
3119  }
3120  }
3121  }
3122  m_dlnActCoeffMolaldT_Unscaled[i] =
3123  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
3124  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
3125  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3126  std::string sni = speciesName(i);
3127  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
3128  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
3129  printf(" %12g %12g %12g %12g %12g %12g\n",
3130  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
3131  }
3132  }
3133  /*
3134  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
3135  * ------ -> equations agree with my notes,
3136  * -> Equations agree with Pitzer,
3137  */
3138  if (charge(i) == 0.0) {
3139  double sum1 = 0.0;
3140  double sum3 = 0.0;
3141  for (size_t j = 1; j < m_kk; j++) {
3142  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_L(i,j);
3143  /*
3144  * Zeta term -> we piggyback on the psi term
3145  */
3146  if (charge(j) > 0.0) {
3147  for (size_t k = 1; k < m_kk; k++) {
3148  if (charge(k) < 0.0) {
3149  size_t n = k + j * m_kk + i * m_kk * m_kk;
3150  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_L[n];
3151  }
3152  }
3153  }
3154  }
3155  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
3156  m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
3157  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
3158  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3159  std::string sni = speciesName(i);
3160  printf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
3161  sni.c_str(), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
3162  }
3163  }
3164 
3165  }
3166  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3167  printf(" Step 9: \n");
3168  }
3169  /*
3170  * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
3171  *
3172  */
3173  double sum1 = 0.0;
3174  double sum2 = 0.0;
3175  double sum3 = 0.0;
3176  double sum4 = 0.0;
3177  double sum5 = 0.0;
3178  double sum6 = 0.0;
3179  double sum7 = 0.0;
3180  /*
3181  * term1 is the temperature derivative of the
3182  * DH term in the osmotic coefficient expression
3183  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
3184  * implementations.
3185  * Is = Ionic strength on the molality scale (units of (gmol/kg))
3186  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3187  */
3188  double term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3189 
3190  for (size_t j = 1; j < m_kk; j++) {
3191  /*
3192  * Loop Over Cations
3193  */
3194  if (charge(j) > 0.0) {
3195  for (size_t k = 1; k < m_kk; k++) {
3196  if (charge(k) < 0.0) {
3197  /*
3198  * Find the counterIJ for the symmetric j,k binary interaction
3199  */
3200  size_t n = m_kk*j + k;
3201  size_t counterIJ = m_CounterIJ[n];
3202 
3203  sum1 = sum1 + molality[j]*molality[k]*
3204  (BphiMX_L[counterIJ] + molarcharge*CMX_L[counterIJ]);
3205  }
3206  }
3207 
3208  for (size_t k = j+1; k < m_kk; k++) {
3209  if (j == (m_kk-1)) {
3210  // we should never reach this step
3211  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
3212  "logic error 1 in Step 9 of hmw_act");
3213  }
3214  if (charge(k) > 0.0) {
3215  /*
3216  * Find the counterIJ for the symmetric j,k binary interaction
3217  * between 2 cations.
3218  */
3219  size_t n = m_kk*j + k;
3220  size_t counterIJ = m_CounterIJ[n];
3221  sum2 = sum2 + molality[j]*molality[k]*Phiphi_L[counterIJ];
3222  for (size_t m = 1; m < m_kk; m++) {
3223  if (charge(m) < 0.0) {
3224  // species m is an anion
3225  n = m + k * m_kk + j * m_kk * m_kk;
3226  sum2 = sum2 +
3227  molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
3228  }
3229  }
3230  }
3231  }
3232  }
3233 
3234  /*
3235  * Loop Over Anions
3236  */
3237  if (charge(j) < 0) {
3238  for (size_t k = j+1; k < m_kk; k++) {
3239  if (j == m_kk-1) {
3240  // we should never reach this step
3241  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
3242  "logic error 2 in Step 9 of hmw_act");
3243  }
3244  if (charge(k) < 0) {
3245  /*
3246  * Find the counterIJ for the symmetric j,k binary interaction
3247  * between two anions
3248  */
3249  size_t n = m_kk*j + k;
3250  size_t counterIJ = m_CounterIJ[n];
3251 
3252  sum3 = sum3 + molality[j]*molality[k]*Phiphi_L[counterIJ];
3253  for (size_t m = 1; m < m_kk; m++) {
3254  if (charge(m) > 0.0) {
3255  n = m + k * m_kk + j * m_kk * m_kk;
3256  sum3 = sum3 +
3257  molality[j]*molality[k]*molality[m]*psi_ijk_L[n];
3258  }
3259  }
3260  }
3261  }
3262  }
3263 
3264  /*
3265  * Loop Over Neutral Species
3266  */
3267  if (charge(j) == 0) {
3268  for (size_t k = 1; k < m_kk; k++) {
3269  if (charge(k) < 0.0) {
3270  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3271  }
3272  if (charge(k) > 0.0) {
3273  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3274  }
3275  if (charge(k) == 0.0) {
3276  if (k > j) {
3277  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3278  } else if (k == j) {
3279  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k);
3280  }
3281  }
3282  if (charge(k) < 0.0) {
3283  size_t izeta = j;
3284  for (size_t m = 1; m < m_kk; m++) {
3285  if (charge(m) > 0.0) {
3286  size_t jzeta = m;
3287  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3288  double zeta_L = psi_ijk_L[n];
3289  if (zeta_L != 0.0) {
3290  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
3291  }
3292  }
3293  }
3294  }
3295  }
3296  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j];
3297  }
3298  }
3299  double sum_m_phi_minus_1 = 2.0 *
3300  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3301  /*
3302  * Calculate the osmotic coefficient from
3303  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3304  */
3305  double d_osmotic_coef_dT;
3306  if (molalitysum > 1.0E-150) {
3307  d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3308  } else {
3309  d_osmotic_coef_dT = 0.0;
3310  }
3311 
3312  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3313  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
3314  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
3315  term1, sum1, sum2, sum3, sum4, sum5);
3316  printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n",
3317  sum_m_phi_minus_1, d_osmotic_coef_dT);
3318  printf(" Step 10: \n");
3319  }
3320  double d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
3321 
3322  /*
3323  * In Cantera, we define the activity coefficient of the solvent as
3324  *
3325  * act_0 = actcoeff_0 * Xmol_0
3326  *
3327  * We have just computed act_0. However, this routine returns
3328  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3329  */
3330  m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
3331  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3332  double d_wateract_dT = exp(d_lnwateract_dT);
3333  printf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n",
3334  d_lnwateract_dT, d_wateract_dT);
3335  }
3336 }
3337 
3338 void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
3339 {
3340  static const int cacheId = m_cache.getId();
3341  CachedScalar cached = m_cache.getScalar(cacheId);
3342  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
3343  return;
3344  }
3345 
3346  /*
3347  * Zero the unscaled 2nd derivatives
3348  */
3349  m_d2lnActCoeffMolaldT2_Unscaled.assign(m_kk, 0.0);
3350  /*
3351  * Calculate the unscaled 2nd derivatives
3352  */
3353  s_updatePitzer_d2lnMolalityActCoeff_dT2();
3354 
3355  for (size_t k = 1; k < m_kk; k++) {
3356  if (CROP_speciesCropped_[k] == 2) {
3357  m_d2lnActCoeffMolaldT2_Unscaled[k] = 0.0;
3358  }
3359  }
3360 
3361  if (CROP_speciesCropped_[0]) {
3362  m_d2lnActCoeffMolaldT2_Unscaled[0] = 0.0;
3363  }
3364 
3365  /*
3366  * Scale the 2nd derivatives
3367  */
3368  s_updateScaling_pHScaling_dT2();
3369 }
3370 
3371 void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
3372 {
3373  /*
3374  * HKM -> Assumption is made that the solvent is species 0.
3375  */
3376 #ifdef DEBUG_MODE
3377  m_debugCalc = 0;
3378 #endif
3379  if (m_indexSolvent != 0) {
3380  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3381  "Wrong index solvent value!");
3382  }
3383 
3384  const double* molality = DATA_PTR(m_molalitiesCropped);
3385  const double* beta0MX_LL= DATA_PTR(m_Beta0MX_ij_LL);
3386  const double* beta1MX_LL= DATA_PTR(m_Beta1MX_ij_LL);
3387  const double* beta2MX_LL= DATA_PTR(m_Beta2MX_ij_LL);
3388  const double* CphiMX_LL = DATA_PTR(m_CphiMX_ij_LL);
3389  const double* thetaij_LL= DATA_PTR(m_Theta_ij_LL);
3390  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
3391  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
3392  const double* psi_ijk_LL= DATA_PTR(m_Psi_ijk_LL);
3393 
3394  /*
3395  * Local variables defined by Coltrin
3396  */
3397  double etheta[5][5], etheta_prime[5][5], sqrtIs;
3398  /*
3399  * Molality based ionic strength of the solution
3400  */
3401  double Is = 0.0;
3402  /*
3403  * Molarcharge of the solution: In Pitzer's notation,
3404  * this is his variable called "Z".
3405  */
3406  double molarcharge = 0.0;
3407  /*
3408  * molalitysum is the sum of the molalities over all solutes,
3409  * even those with zero charge.
3410  */
3411  double molalitysum = 0.0;
3412 
3413  double* gfunc = DATA_PTR(m_gfunc_IJ);
3414  double* g2func = DATA_PTR(m_g2func_IJ);
3415  double* hfunc = DATA_PTR(m_hfunc_IJ);
3416  double* h2func = DATA_PTR(m_h2func_IJ);
3417  double* BMX_LL = DATA_PTR(m_BMX_IJ_LL);
3418  double* BprimeMX_LL=DATA_PTR(m_BprimeMX_IJ_LL);
3419  double* BphiMX_LL= DATA_PTR(m_BphiMX_IJ_LL);
3420  double* Phi_LL = DATA_PTR(m_Phi_IJ_LL);
3421  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
3422  double* Phiphi_LL= DATA_PTR(m_PhiPhi_IJ_LL);
3423  double* CMX_LL = DATA_PTR(m_CMX_IJ_LL);
3424 
3425  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3426  printf("\n Debugging information from "
3427  "s_Pitzer_d2lnMolalityActCoeff_dT2()\n");
3428  }
3429  /*
3430  * Make sure the counter variables are setup
3431  */
3432  counterIJ_setup();
3433 
3434 
3435  /*
3436  * ---------- Calculate common sums over solutes ---------------------
3437  */
3438  for (size_t n = 1; n < m_kk; n++) {
3439  // ionic strength
3440  Is += charge(n) * charge(n) * molality[n];
3441  // total molar charge
3442  molarcharge += fabs(charge(n)) * molality[n];
3443  molalitysum += molality[n];
3444  }
3445  Is *= 0.5;
3446 
3447  /*
3448  * Store the ionic molality in the object for reference.
3449  */
3450  m_IionicMolality = Is;
3451  sqrtIs = sqrt(Is);
3452  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3453  printf(" Step 1: \n");
3454  printf(" ionic strenth = %14.7le \n total molar "
3455  "charge = %14.7le \n", Is, molarcharge);
3456  }
3457 
3458  /*
3459  * The following call to calc_lambdas() calculates all 16 elements
3460  * of the elambda and elambda1 arrays, given the value of the
3461  * ionic strength (Is)
3462  */
3463  calc_lambdas(Is);
3464 
3465  /*
3466  * ----- Step 2: Find the coefficients E-theta and -------------------
3467  * E-thetaprime for all combinations of positive
3468  * unlike charges up to 4
3469  */
3470  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3471  printf(" Step 2: \n");
3472  }
3473  for (int z1 = 1; z1 <=4; z1++) {
3474  for (int z2 =1; z2 <=4; z2++) {
3475  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3476  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3477  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3478  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3479  }
3480  }
3481  }
3482 
3483  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3484  printf(" Step 3: \n");
3485  printf(" Species Species g(x) "
3486  " hfunc(x) \n");
3487  }
3488 
3489  /*
3490  *
3491  * calculate gfunc(x) and hfunc(x) for each cation-anion pair MX
3492  * In the original literature, hfunc, was called gprime. However,
3493  * it's not the derivative of gfunc(x), so I renamed it.
3494  */
3495  for (size_t i = 1; i < (m_kk - 1); i++) {
3496  for (size_t j = (i+1); j < m_kk; j++) {
3497  /*
3498  * Find the counterIJ for the symmetric binary interaction
3499  */
3500  size_t n = m_kk*i + j;
3501  size_t counterIJ = m_CounterIJ[n];
3502  /*
3503  * Only loop over oppositely charge species
3504  */
3505  if (charge(i)*charge(j) < 0) {
3506  /*
3507  * x is a reduced function variable
3508  */
3509  double x1 = sqrtIs * alpha1MX[counterIJ];
3510  if (x1 > 1.0E-100) {
3511  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
3512  hfunc[counterIJ] = -2.0*
3513  (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
3514  } else {
3515  gfunc[counterIJ] = 0.0;
3516  hfunc[counterIJ] = 0.0;
3517  }
3518 
3519  if (beta2MX_LL[counterIJ] != 0.0) {
3520  double x2 = sqrtIs * alpha2MX[counterIJ];
3521  if (x2 > 1.0E-100) {
3522  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3523  h2func[counterIJ] = -2.0 *
3524  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3525  } else {
3526  g2func[counterIJ] = 0.0;
3527  h2func[counterIJ] = 0.0;
3528  }
3529  }
3530  } else {
3531  gfunc[counterIJ] = 0.0;
3532  hfunc[counterIJ] = 0.0;
3533  }
3534  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3535  std::string sni = speciesName(i);
3536  std::string snj = speciesName(j);
3537  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
3538  gfunc[counterIJ], hfunc[counterIJ]);
3539  }
3540  }
3541  }
3542  /*
3543  * ------- SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L ----------
3544  * ------- These are now temperature derivatives of the
3545  * previously calculated quantities.
3546  */
3547  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3548  printf(" Step 4: \n");
3549  printf(" Species Species BMX "
3550  "BprimeMX BphiMX \n");
3551  }
3552 
3553  for (size_t i = 1; i < m_kk - 1; i++) {
3554  for (size_t j = i+1; j < m_kk; j++) {
3555  /*
3556  * Find the counterIJ for the symmetric binary interaction
3557  */
3558  size_t n = m_kk*i + j;
3559  size_t counterIJ = m_CounterIJ[n];
3560  /*
3561  * both species have a non-zero charge, and one is positive
3562  * and the other is negative
3563  */
3564  if (charge(i)*charge(j) < 0.0) {
3565  BMX_LL[counterIJ] = beta0MX_LL[counterIJ]
3566  + beta1MX_LL[counterIJ] * gfunc[counterIJ]
3567  + beta2MX_LL[counterIJ] * g2func[counterIJ];
3568  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3569  printf("%d %g: %g %g %g %g\n",
3570  (int) counterIJ, BMX_LL[counterIJ], beta0MX_LL[counterIJ],
3571  beta1MX_LL[counterIJ], beta2MX_LL[counterIJ], gfunc[counterIJ]);
3572  }
3573  if (Is > 1.0E-150) {
3574  BprimeMX_LL[counterIJ] = (beta1MX_LL[counterIJ] * hfunc[counterIJ]/Is +
3575  beta2MX_LL[counterIJ] * h2func[counterIJ]/Is);
3576  } else {
3577  BprimeMX_LL[counterIJ] = 0.0;
3578  }
3579  BphiMX_LL[counterIJ] = BMX_LL[counterIJ] + Is*BprimeMX_LL[counterIJ];
3580  } else {
3581  BMX_LL[counterIJ] = 0.0;
3582  BprimeMX_LL[counterIJ] = 0.0;
3583  BphiMX_LL[counterIJ] = 0.0;
3584  }
3585  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3586  std::string sni = speciesName(i);
3587  std::string snj = speciesName(j);
3588  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
3589  sni.c_str(), snj.c_str(),
3590  BMX_LL[counterIJ], BprimeMX_LL[counterIJ], BphiMX_LL[counterIJ]);
3591  }
3592  }
3593  }
3594 
3595  /*
3596  * --------- SUBSECTION TO CALCULATE CMX_LL ----------
3597  * ---------
3598  */
3599  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3600  printf(" Step 5: \n");
3601  printf(" Species Species CMX \n");
3602  }
3603  for (size_t i = 1; i < m_kk-1; i++) {
3604  for (size_t j = i+1; j < m_kk; j++) {
3605  /*
3606  * Find the counterIJ for the symmetric binary interaction
3607  */
3608  size_t n = m_kk*i + j;
3609  size_t counterIJ = m_CounterIJ[n];
3610  /*
3611  * both species have a non-zero charge, and one is positive
3612  * and the other is negative
3613  */
3614  if (charge(i)*charge(j) < 0.0) {
3615  CMX_LL[counterIJ] = CphiMX_LL[counterIJ]/
3616  (2.0* sqrt(fabs(charge(i)*charge(j))));
3617  } else {
3618  CMX_LL[counterIJ] = 0.0;
3619  }
3620  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3621  std::string sni = speciesName(i);
3622  std::string snj = speciesName(j);
3623  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
3624  CMX_LL[counterIJ]);
3625  }
3626  }
3627  }
3628 
3629  /*
3630  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
3631  * --------
3632  */
3633  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3634  printf(" Step 6: \n");
3635  printf(" Species Species Phi_ij "
3636  " Phiprime_ij Phi^phi_ij \n");
3637  }
3638  for (size_t i = 1; i < m_kk-1; i++) {
3639  for (size_t j = i+1; j < m_kk; j++) {
3640  /*
3641  * Find the counterIJ for the symmetric binary interaction
3642  */
3643  size_t n = m_kk*i + j;
3644  size_t counterIJ = m_CounterIJ[n];
3645  /*
3646  * both species have a non-zero charge, and one is positive
3647  * and the other is negative
3648  */
3649  if (charge(i)*charge(j) > 0) {
3650  Phi_LL[counterIJ] = thetaij_LL[counterIJ];
3651  Phiprime[counterIJ] = 0.0;
3652  Phiphi_LL[counterIJ] = Phi_LL[counterIJ];
3653  } else {
3654  Phi_LL[counterIJ] = 0.0;
3655  Phiprime[counterIJ] = 0.0;
3656  Phiphi_LL[counterIJ] = 0.0;
3657  }
3658  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3659  std::string sni = speciesName(i);
3660  std::string snj = speciesName(j);
3661  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
3662  sni.c_str(), snj.c_str(),
3663  Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ]);
3664  }
3665  }
3666  }
3667 
3668  /*
3669  * ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
3670  */
3671  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3672  printf(" Step 7: \n");
3673  }
3674  double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
3675  double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3676  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3677  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3678  printf(" initial value of d2FdT2 = %10.6f \n", d2FdT2);
3679  }
3680  for (size_t i = 1; i < m_kk-1; i++) {
3681  for (size_t j = i+1; j < m_kk; j++) {
3682  /*
3683  * Find the counterIJ for the symmetric binary interaction
3684  */
3685  size_t n = m_kk*i + j;
3686  size_t counterIJ = m_CounterIJ[n];
3687  /*
3688  * both species have a non-zero charge, and one is positive
3689  * and the other is negative
3690  */
3691  if (charge(i)*charge(j) < 0) {
3692  d2FdT2 = d2FdT2 + molality[i]*molality[j] * BprimeMX_LL[counterIJ];
3693  }
3694  /*
3695  * Both species have a non-zero charge, and they
3696  * have the same sign, e.g., both positive or both negative.
3697  */
3698  if (charge(i)*charge(j) > 0) {
3699  d2FdT2 = d2FdT2 + molality[i]*molality[j] * Phiprime[counterIJ];
3700  }
3701  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3702  printf(" d2FdT2 = %10.6f \n", d2FdT2);
3703  }
3704  }
3705  }
3706  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3707  printf(" Step 8: \n");
3708  }
3709 
3710  for (size_t i = 1; i < m_kk; i++) {
3711 
3712  /*
3713  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
3714  * --
3715  */
3716  if (charge(i) > 0) {
3717  // species i is the cation (positive) to calc the actcoeff
3718  double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3719  double sum1 = 0.0;
3720  double sum2 = 0.0;
3721  double sum3 = 0.0;
3722  double sum4 = 0.0;
3723  double sum5 = 0.0;
3724  for (size_t j = 1; j < m_kk; j++) {
3725  /*
3726  * Find the counterIJ for the symmetric binary interaction
3727  */
3728  size_t n = m_kk*i + j;
3729  size_t counterIJ = m_CounterIJ[n];
3730 
3731  if (charge(j) < 0.0) {
3732  // sum over all anions
3733  sum1 = sum1 + molality[j]*
3734  (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
3735  if (j < m_kk-1) {
3736  /*
3737  * This term is the ternary interaction involving the
3738  * non-duplicate sum over double anions, j, k, with
3739  * respect to the cation, i.
3740  */
3741  for (size_t k = j+1; k < m_kk; k++) {
3742  // an inner sum over all anions
3743  if (charge(k) < 0.0) {
3744  n = k + j * m_kk + i * m_kk * m_kk;
3745  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
3746  }
3747  }
3748  }
3749  }
3750 
3751 
3752  if (charge(j) > 0.0) {
3753  // sum over all cations
3754  if (j != i) {
3755  sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
3756  }
3757  for (size_t k = 1; k < m_kk; k++) {
3758  if (charge(k) < 0.0) {
3759  // two inner sums over anions
3760 
3761  n = k + j * m_kk + i * m_kk * m_kk;
3762  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
3763  /*
3764  * Find the counterIJ for the j,k interaction
3765  */
3766  n = m_kk*j + k;
3767  size_t counterIJ2 = m_CounterIJ[n];
3768  sum4 = sum4 + (fabs(charge(i))*
3769  molality[j]*molality[k]*CMX_LL[counterIJ2]);
3770  }
3771  }
3772  }
3773 
3774  /*
3775  * Handle neutral j species
3776  */
3777  if (charge(j) == 0) {
3778  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i);
3779  /*
3780  * Zeta interaction term
3781  */
3782  for (size_t k = 1; k < m_kk; k++) {
3783  if (charge(k) < 0.0) {
3784  size_t izeta = j;
3785  size_t jzeta = i;
3786  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3787  double zeta_LL = psi_ijk_LL[n];
3788  if (zeta_LL != 0.0) {
3789  sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
3790  }
3791  }
3792  }
3793  }
3794  }
3795  /*
3796  * Add all of the contributions up to yield the log of the
3797  * solute activity coefficients (molality scale)
3798  */
3799  m_d2lnActCoeffMolaldT2_Unscaled[i] =
3800  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3801  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3802  std::string sni = speciesName(i);
3803  printf(" %-16s d2lngammadT2[i]=%10.6f \n",
3804  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
3805  printf(" %12g %12g %12g %12g %12g %12g\n",
3806  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
3807  }
3808  }
3809 
3810 
3811  /*
3812  * ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------
3813  *
3814  */
3815  if (charge(i) < 0) {
3816  // species i is an anion (negative)
3817  double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3818  double sum1 = 0.0;
3819  double sum2 = 0.0;
3820  double sum3 = 0.0;
3821  double sum4 = 0.0;
3822  double sum5 = 0.0;
3823  for (size_t j = 1; j < m_kk; j++) {
3824  /*
3825  * Find the counterIJ for the symmetric binary interaction
3826  */
3827  size_t n = m_kk*i + j;
3828  size_t counterIJ = m_CounterIJ[n];
3829 
3830  /*
3831  * For Anions, do the cation interactions.
3832  */
3833  if (charge(j) > 0) {
3834  sum1 = sum1 + molality[j]*
3835  (2.0*BMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
3836  if (j < m_kk-1) {
3837  for (size_t k = j+1; k < m_kk; k++) {
3838  // an inner sum over all cations
3839  if (charge(k) > 0) {
3840  n = k + j * m_kk + i * m_kk * m_kk;
3841  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
3842  }
3843  }
3844  }
3845  }
3846 
3847  /*
3848  * For Anions, do the other anion interactions.
3849  */
3850  if (charge(j) < 0.0) {
3851  // sum over all anions
3852  if (j != i) {
3853  sum2 = sum2 + molality[j]*(2.0*Phi_LL[counterIJ]);
3854  }
3855  for (size_t k = 1; k < m_kk; k++) {
3856  if (charge(k) > 0.0) {
3857  // two inner sums over cations
3858  n = k + j * m_kk + i * m_kk * m_kk;
3859  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_LL[n];
3860  /*
3861  * Find the counterIJ for the symmetric binary interaction
3862  */
3863  n = m_kk*j + k;
3864  size_t counterIJ2 = m_CounterIJ[n];
3865  sum4 = sum4 +
3866  (fabs(charge(i))*
3867  molality[j]*molality[k]*CMX_LL[counterIJ2]);
3868  }
3869  }
3870  }
3871 
3872  /*
3873  * for Anions, do the neutral species interaction
3874  */
3875  if (charge(j) == 0.0) {
3876  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_LL(j,i);
3877  /*
3878  * Zeta interaction term
3879  */
3880  for (size_t k = 1; k < m_kk; k++) {
3881  if (charge(k) > 0.0) {
3882  size_t izeta = j;
3883  size_t jzeta = k;
3884  size_t kzeta = i;
3885  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3886  double zeta_LL = psi_ijk_LL[n];
3887  if (zeta_LL != 0.0) {
3888  sum5 = sum5 + molality[j]*molality[k]*zeta_LL;
3889  }
3890  }
3891  }
3892  }
3893  }
3894  m_d2lnActCoeffMolaldT2_Unscaled[i] =
3895  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3896  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3897  std::string sni = speciesName(i);
3898  printf(" %-16s d2lngammadT2[i]=%10.6f\n",
3899  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
3900  printf(" %12g %12g %12g %12g %12g %12g\n",
3901  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
3902  }
3903  }
3904  /*
3905  * ------ SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF -------
3906  * ------ -> equations agree with my notes,
3907  * -> Equations agree with Pitzer,
3908  */
3909  if (charge(i) == 0.0) {
3910  double sum1 = 0.0;
3911  double sum3 = 0.0;
3912  for (size_t j = 1; j < m_kk; j++) {
3913  sum1 = sum1 + molality[j]*2.0*m_Lambda_nj_LL(i,j);
3914  /*
3915  * Zeta term -> we piggyback on the psi term
3916  */
3917  if (charge(j) > 0.0) {
3918  for (size_t k = 1; k < m_kk; k++) {
3919  if (charge(k) < 0.0) {
3920  size_t n = k + j * m_kk + i * m_kk * m_kk;
3921  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_LL[n];
3922  }
3923  }
3924  }
3925  }
3926  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
3927  m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
3928  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3929  std::string sni = speciesName(i);
3930  printf(" %-16s d2lngammadT2[i]=%10.6f \n",
3931  sni.c_str(), m_d2lnActCoeffMolaldT2_Unscaled[i]);
3932  }
3933  }
3934 
3935  }
3936  if (DEBUG_MODE_ENABLED && m_debugCalc) {
3937  printf(" Step 9: \n");
3938  }
3939 
3940  /*
3941  * ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
3942  */
3943  double sum1 = 0.0;
3944  double sum2 = 0.0;
3945  double sum3 = 0.0;
3946  double sum4 = 0.0;
3947  double sum5 = 0.0;
3948  double sum6 = 0.0;
3949  double sum7 = 0.0;
3950  /*
3951  * term1 is the temperature derivative of the
3952  * DH term in the osmotic coefficient expression
3953  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
3954  * implementations.
3955  * Is = Ionic strength on the molality scale (units of (gmol/kg))
3956  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3957  */
3958  double term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3959 
3960  for (size_t j = 1; j < m_kk; j++) {
3961  /*
3962  * Loop Over Cations
3963  */
3964  if (charge(j) > 0.0) {
3965  for (size_t k = 1; k < m_kk; k++) {
3966  if (charge(k) < 0.0) {
3967  /*
3968  * Find the counterIJ for the symmetric j,k binary interaction
3969  */
3970  size_t n = m_kk*j + k;
3971  size_t counterIJ = m_CounterIJ[n];
3972 
3973  sum1 = sum1 + molality[j]*molality[k]*
3974  (BphiMX_LL[counterIJ] + molarcharge*CMX_LL[counterIJ]);
3975  }
3976  }
3977 
3978  for (size_t k = j+1; k < m_kk; k++) {
3979  if (j == (m_kk-1)) {
3980  // we should never reach this step
3981  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3982  "logic error 1 in Step 9 of hmw_act");
3983  }
3984  if (charge(k) > 0.0) {
3985  /*
3986  * Find the counterIJ for the symmetric j,k binary interaction
3987  * between 2 cations.
3988  */
3989  size_t n = m_kk*j + k;
3990  size_t counterIJ = m_CounterIJ[n];
3991  sum2 = sum2 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
3992  for (size_t m = 1; m < m_kk; m++) {
3993  if (charge(m) < 0.0) {
3994  // species m is an anion
3995  n = m + k * m_kk + j * m_kk * m_kk;
3996  sum2 = sum2 +
3997  molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
3998  }
3999  }
4000  }
4001  }
4002  }
4003 
4004  /*
4005  * Loop Over Anions
4006  */
4007  if (charge(j) < 0) {
4008  for (size_t k = j+1; k < m_kk; k++) {
4009  if (j == m_kk-1) {
4010  // we should never reach this step
4011  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
4012  "logic error 2 in Step 9 of hmw_act");
4013  }
4014  if (charge(k) < 0) {
4015  /*
4016  * Find the counterIJ for the symmetric j,k binary interaction
4017  * between two anions
4018  */
4019  size_t n = m_kk*j + k;
4020  size_t counterIJ = m_CounterIJ[n];
4021 
4022  sum3 = sum3 + molality[j]*molality[k]*Phiphi_LL[counterIJ];
4023  for (size_t m = 1; m < m_kk; m++) {
4024  if (charge(m) > 0.0) {
4025  n = m + k * m_kk + j * m_kk * m_kk;
4026  sum3 = sum3 +
4027  molality[j]*molality[k]*molality[m]*psi_ijk_LL[n];
4028  }
4029  }
4030  }
4031  }
4032  }
4033 
4034  /*
4035  * Loop Over Neutral Species
4036  */
4037  if (charge(j) == 0) {
4038  for (size_t k = 1; k < m_kk; k++) {
4039  if (charge(k) < 0.0) {
4040  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4041  }
4042  if (charge(k) > 0.0) {
4043  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4044  }
4045  if (charge(k) == 0.0) {
4046  if (k > j) {
4047  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4048  } else if (k == j) {
4049  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
4050  }
4051  }
4052  if (charge(k) < 0.0) {
4053  size_t izeta = j;
4054  for (size_t m = 1; m < m_kk; m++) {
4055  if (charge(m) > 0.0) {
4056  size_t jzeta = m;
4057  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
4058  double zeta_LL = psi_ijk_LL[n];
4059  if (zeta_LL != 0.0) {
4060  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
4061  }
4062  }
4063  }
4064  }
4065  }
4066 
4067  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j];
4068  }
4069  }
4070  double sum_m_phi_minus_1 = 2.0 *
4071  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
4072  /*
4073  * Calculate the osmotic coefficient from
4074  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
4075  */
4076  double d2_osmotic_coef_dT2;
4077  if (molalitysum > 1.0E-150) {
4078  d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
4079  } else {
4080  d2_osmotic_coef_dT2 = 0.0;
4081  }
4082  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4083  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
4084  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
4085  term1, sum1, sum2, sum3, sum4, sum5);
4086  printf(" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n",
4087  sum_m_phi_minus_1, d2_osmotic_coef_dT2);
4088  printf(" Step 10: \n");
4089  }
4090  double d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
4091 
4092  /*
4093  * In Cantera, we define the activity coefficient of the solvent as
4094  *
4095  * act_0 = actcoeff_0 * Xmol_0
4096  *
4097  * We have just computed act_0. However, this routine returns
4098  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
4099  */
4100  m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
4101 
4102  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4103  double d2_wateract_dT2 = exp(d2_lnwateract_dT2);
4104  printf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n",
4105  d2_lnwateract_dT2, d2_wateract_dT2);
4106  }
4107 }
4108 
4109 void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
4110 {
4111  static const int cacheId = m_cache.getId();
4112  CachedScalar cached = m_cache.getScalar(cacheId);
4113  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
4114  return;
4115  }
4116 
4117  m_dlnActCoeffMolaldP_Unscaled.assign(m_kk, 0.0);
4118  s_updatePitzer_dlnMolalityActCoeff_dP();
4119 
4120  for (size_t k = 1; k < m_kk; k++) {
4121  if (CROP_speciesCropped_[k] == 2) {
4122  m_dlnActCoeffMolaldP_Unscaled[k] = 0.0;
4123  }
4124  }
4125 
4126  if (CROP_speciesCropped_[0]) {
4127  m_dlnActCoeffMolaldP_Unscaled[0] = 0.0;
4128  }
4129 
4130  s_updateScaling_pHScaling_dP();
4131 }
4132 
4133 void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
4134 {
4135  /*
4136  * HKM -> Assumption is made that the solvent is species 0.
4137  */
4138 #ifdef DEBUG_MODE
4139  m_debugCalc = 0;
4140 #endif
4141  if (m_indexSolvent != 0) {
4142  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
4143  "Wrong index solvent value!");
4144  }
4145 
4146  const double* molality = DATA_PTR(m_molalitiesCropped);
4147  const double* beta0MX_P = DATA_PTR(m_Beta0MX_ij_P);
4148  const double* beta1MX_P = DATA_PTR(m_Beta1MX_ij_P);
4149  const double* beta2MX_P = DATA_PTR(m_Beta2MX_ij_P);
4150  const double* CphiMX_P = DATA_PTR(m_CphiMX_ij_P);
4151  const double* thetaij_P = DATA_PTR(m_Theta_ij_P);
4152  const double* alpha1MX = DATA_PTR(m_Alpha1MX_ij);
4153  const double* alpha2MX = DATA_PTR(m_Alpha2MX_ij);
4154  const double* psi_ijk_P = DATA_PTR(m_Psi_ijk_P);
4155 
4156  /*
4157  * Local variables defined by Coltrin
4158  */
4159  double etheta[5][5], etheta_prime[5][5], sqrtIs;
4160  /*
4161  * Molality based ionic strength of the solution
4162  */
4163  double Is = 0.0;
4164  /*
4165  * Molarcharge of the solution: In Pitzer's notation,
4166  * this is his variable called "Z".
4167  */
4168  double molarcharge = 0.0;
4169  /*
4170  * molalitysum is the sum of the molalities over all solutes,
4171  * even those with zero charge.
4172  */
4173  double molalitysum = 0.0;
4174 
4175  double* gfunc = DATA_PTR(m_gfunc_IJ);
4176  double* g2func = DATA_PTR(m_g2func_IJ);
4177  double* hfunc = DATA_PTR(m_hfunc_IJ);
4178  double* h2func = DATA_PTR(m_h2func_IJ);
4179  double* BMX_P = DATA_PTR(m_BMX_IJ_P);
4180  double* BprimeMX_P= DATA_PTR(m_BprimeMX_IJ_P);
4181  double* BphiMX_P = DATA_PTR(m_BphiMX_IJ_P);
4182  double* Phi_P = DATA_PTR(m_Phi_IJ_P);
4183  double* Phiprime = DATA_PTR(m_Phiprime_IJ);
4184  double* Phiphi_P = DATA_PTR(m_PhiPhi_IJ_P);
4185  double* CMX_P = DATA_PTR(m_CMX_IJ_P);
4186 
4187  double currTemp = temperature();
4188  double currPres = pressure();
4189 
4190  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4191  printf("\n Debugging information from "
4192  "s_Pitzer_dlnMolalityActCoeff_dP()\n");
4193  }
4194  /*
4195  * Make sure the counter variables are setup
4196  */
4197  counterIJ_setup();
4198 
4199  /*
4200  * ---------- Calculate common sums over solutes ---------------------
4201  */
4202  for (size_t n = 1; n < m_kk; n++) {
4203  // ionic strength
4204  Is += charge(n) * charge(n) * molality[n];
4205  // total molar charge
4206  molarcharge += fabs(charge(n)) * molality[n];
4207  molalitysum += molality[n];
4208  }
4209  Is *= 0.5;
4210 
4211  /*
4212  * Store the ionic molality in the object for reference.
4213  */
4214  m_IionicMolality = Is;
4215  sqrtIs = sqrt(Is);
4216  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4217  printf(" Step 1: \n");
4218  printf(" ionic strenth = %14.7le \n total molar "
4219  "charge = %14.7le \n", Is, molarcharge);
4220  }
4221 
4222  /*
4223  * The following call to calc_lambdas() calculates all 16 elements
4224  * of the elambda and elambda1 arrays, given the value of the
4225  * ionic strength (Is)
4226  */
4227  calc_lambdas(Is);
4228 
4229 
4230  /*
4231  * ----- Step 2: Find the coefficients E-theta and -------------------
4232  * E-thetaprime for all combinations of positive
4233  * unlike charges up to 4
4234  */
4235  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4236  printf(" Step 2: \n");
4237  }
4238  for (int z1 = 1; z1 <=4; z1++) {
4239  for (int z2 =1; z2 <=4; z2++) {
4240  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
4241  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4242  printf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
4243  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
4244  }
4245  }
4246  }
4247 
4248  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4249  printf(" Step 3: \n");
4250  printf(" Species Species g(x) "
4251  " hfunc(x) \n");
4252  }
4253 
4254  /*
4255  *
4256  * calculate g(x) and hfunc(x) for each cation-anion pair MX
4257  * In the original literature, hfunc, was called gprime. However,
4258  * it's not the derivative of g(x), so I renamed it.
4259  */
4260  for (size_t i = 1; i < (m_kk - 1); i++) {
4261  for (size_t j = (i+1); j < m_kk; j++) {
4262  /*
4263  * Find the counterIJ for the symmetric binary interaction
4264  */
4265  size_t n = m_kk*i + j;
4266  size_t counterIJ = m_CounterIJ[n];
4267  /*
4268  * Only loop over oppositely charge species
4269  */
4270  if (charge(i)*charge(j) < 0) {
4271  /*
4272  * x is a reduced function variable
4273  */
4274  double x1 = sqrtIs * alpha1MX[counterIJ];
4275  if (x1 > 1.0E-100) {
4276  gfunc[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
4277  hfunc[counterIJ] = -2.0*
4278  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
4279  } else {
4280  gfunc[counterIJ] = 0.0;
4281  hfunc[counterIJ] = 0.0;
4282  }
4283 
4284  if (beta2MX_P[counterIJ] != 0.0) {
4285  double x2 = sqrtIs * alpha2MX[counterIJ];
4286  if (x2 > 1.0E-100) {
4287  g2func[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
4288  h2func[counterIJ] = -2.0 *
4289  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
4290  } else {
4291  g2func[counterIJ] = 0.0;
4292  h2func[counterIJ] = 0.0;
4293  }
4294  }
4295  } else {
4296  gfunc[counterIJ] = 0.0;
4297  hfunc[counterIJ] = 0.0;
4298  }
4299  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4300  std::string sni = speciesName(i);
4301  std::string snj = speciesName(j);
4302  printf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
4303  gfunc[counterIJ], hfunc[counterIJ]);
4304  }
4305  }
4306  }
4307 
4308  /*
4309  * ------- SUBSECTION TO CALCULATE BMX_P, BprimeMX_P, BphiMX_P ----------
4310  * ------- These are now temperature derivatives of the
4311  * previously calculated quantities.
4312  */
4313  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4314  printf(" Step 4: \n");
4315  printf(" Species Species BMX "
4316  "BprimeMX BphiMX \n");
4317  }
4318 
4319  for (size_t i = 1; i < m_kk - 1; i++) {
4320  for (size_t j = i+1; j < m_kk; j++) {
4321  /*
4322  * Find the counterIJ for the symmetric binary interaction
4323  */
4324  size_t n = m_kk*i + j;
4325  size_t counterIJ = m_CounterIJ[n];
4326  /*
4327  * both species have a non-zero charge, and one is positive
4328  * and the other is negative
4329  */
4330  if (charge(i)*charge(j) < 0.0) {
4331  BMX_P[counterIJ] = beta0MX_P[counterIJ]
4332  + beta1MX_P[counterIJ] * gfunc[counterIJ]
4333  + beta2MX_P[counterIJ] * g2func[counterIJ];
4334  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4335  printf("%d %g: %g %g %g %g\n",
4336  (int) counterIJ, BMX_P[counterIJ], beta0MX_P[counterIJ],
4337  beta1MX_P[counterIJ], beta2MX_P[counterIJ], gfunc[counterIJ]);
4338  }
4339  if (Is > 1.0E-150) {
4340  BprimeMX_P[counterIJ] = (beta1MX_P[counterIJ] * hfunc[counterIJ]/Is +
4341  beta2MX_P[counterIJ] * h2func[counterIJ]/Is);
4342  } else {
4343  BprimeMX_P[counterIJ] = 0.0;
4344  }
4345  BphiMX_P[counterIJ] = BMX_P[counterIJ] + Is*BprimeMX_P[counterIJ];
4346  } else {
4347  BMX_P[counterIJ] = 0.0;
4348  BprimeMX_P[counterIJ] = 0.0;
4349  BphiMX_P[counterIJ] = 0.0;
4350  }
4351  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4352  std::string sni = speciesName(i);
4353  std::string snj = speciesName(j);
4354  printf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
4355  sni.c_str(), snj.c_str(),
4356  BMX_P[counterIJ], BprimeMX_P[counterIJ], BphiMX_P[counterIJ]);
4357  }
4358  }
4359  }
4360 
4361  /*
4362  * --------- SUBSECTION TO CALCULATE CMX_P ----------
4363  */
4364  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4365  printf(" Step 5: \n");
4366  printf(" Species Species CMX \n");
4367  }
4368  for (size_t i = 1; i < m_kk-1; i++) {
4369  for (size_t j = i+1; j < m_kk; j++) {
4370  /*
4371  * Find the counterIJ for the symmetric binary interaction
4372  */
4373  size_t n = m_kk*i + j;
4374  size_t counterIJ = m_CounterIJ[n];
4375  /*
4376  * both species have a non-zero charge, and one is positive
4377  * and the other is negative
4378  */
4379  if (charge(i)*charge(j) < 0.0) {
4380  CMX_P[counterIJ] = CphiMX_P[counterIJ]/
4381  (2.0* sqrt(fabs(charge(i)*charge(j))));
4382  } else {
4383  CMX_P[counterIJ] = 0.0;
4384  }
4385  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4386  std::string sni = speciesName(i);
4387  std::string snj = speciesName(j);
4388  printf(" %-16s %-16s %11.7f \n", sni.c_str(), snj.c_str(),
4389  CMX_P[counterIJ]);
4390  }
4391  }
4392  }
4393 
4394  /*
4395  * ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
4396  * --------
4397  */
4398  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4399  printf(" Step 6: \n");
4400  printf(" Species Species Phi_ij "
4401  " Phiprime_ij Phi^phi_ij \n");
4402  }
4403  for (size_t i = 1; i < m_kk-1; i++) {
4404  for (size_t j = i+1; j < m_kk; j++) {
4405  /*
4406  * Find the counterIJ for the symmetric binary interaction
4407  */
4408  size_t n = m_kk*i + j;
4409  size_t counterIJ = m_CounterIJ[n];
4410  /*
4411  * both species have a non-zero charge, and one is positive
4412  * and the other is negative
4413  */
4414  if (charge(i)*charge(j) > 0) {
4415  Phi_P[counterIJ] = thetaij_P[counterIJ];
4416  Phiprime[counterIJ] = 0.0;
4417  Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ];
4418  } else {
4419  Phi_P[counterIJ] = 0.0;
4420  Phiprime[counterIJ] = 0.0;
4421  Phiphi_P[counterIJ] = 0.0;
4422  }
4423  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4424  std::string sni = speciesName(i);
4425  std::string snj = speciesName(j);
4426  printf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
4427  sni.c_str(), snj.c_str(),
4428  Phi_P[counterIJ], Phiprime[counterIJ], Phiphi_P[counterIJ]);
4429  }
4430  }
4431  }
4432 
4433  /*
4434  * ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
4435  */
4436  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4437  printf(" Step 7: \n");
4438  }
4439  double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
4440  double dAphidP = dA_DebyedP /3.0;
4441  double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
4442  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
4443  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4444  printf(" initial value of dFdP = %10.6f \n", dFdP);
4445  }
4446  for (size_t i = 1; i < m_kk-1; i++) {
4447  for (size_t j = i+1; j < m_kk; j++) {
4448  /*
4449  * Find the counterIJ for the symmetric binary interaction
4450  */
4451  size_t n = m_kk*i + j;
4452  size_t counterIJ = m_CounterIJ[n];
4453  /*
4454  * both species have a non-zero charge, and one is positive
4455  * and the other is negative
4456  */
4457  if (charge(i)*charge(j) < 0) {
4458  dFdP = dFdP + molality[i]*molality[j] * BprimeMX_P[counterIJ];
4459  }
4460  /*
4461  * Both species have a non-zero charge, and they
4462  * have the same sign, e.g., both positive or both negative.
4463  */
4464  if (charge(i)*charge(j) > 0) {
4465  dFdP = dFdP + molality[i]*molality[j] * Phiprime[counterIJ];
4466  }
4467  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4468  printf(" dFdP = %10.6f \n", dFdP);
4469  }
4470  }
4471  }
4472  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4473  printf(" Step 8: \n");
4474  }
4475 
4476  for (size_t i = 1; i < m_kk; i++) {
4477 
4478  /*
4479  * -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
4480  */
4481  if (charge(i) > 0) {
4482  // species i is the cation (positive) to calc the actcoeff
4483  double zsqdFdP = charge(i)*charge(i)*dFdP;
4484  double sum1 = 0.0;
4485  double sum2 = 0.0;
4486  double sum3 = 0.0;
4487  double sum4 = 0.0;
4488  double sum5 = 0.0;
4489  for (size_t j = 1; j < m_kk; j++) {
4490  /*
4491  * Find the counterIJ for the symmetric binary interaction
4492  */
4493  size_t n = m_kk*i + j;
4494  size_t counterIJ = m_CounterIJ[n];
4495 
4496  if (charge(j) < 0.0) {
4497  // sum over all anions
4498  sum1 = sum1 + molality[j]*
4499  (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
4500  if (j < m_kk-1) {
4501  /*
4502  * This term is the ternary interaction involving the
4503  * non-duplicate sum over double anions, j, k, with
4504  * respect to the cation, i.
4505  */
4506  for (size_t k = j+1; k < m_kk; k++) {
4507  // an inner sum over all anions
4508  if (charge(k) < 0.0) {
4509  n = k + j * m_kk + i * m_kk * m_kk;
4510  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
4511  }
4512  }
4513  }
4514  }
4515 
4516 
4517  if (charge(j) > 0.0) {
4518  // sum over all cations
4519  if (j != i) {
4520  sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
4521  }
4522  for (size_t k = 1; k < m_kk; k++) {
4523  if (charge(k) < 0.0) {
4524  // two inner sums over anions
4525 
4526  n = k + j * m_kk + i * m_kk * m_kk;
4527  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
4528  /*
4529  * Find the counterIJ for the j,k interaction
4530  */
4531  n = m_kk*j + k;
4532  size_t counterIJ2 = m_CounterIJ[n];
4533  sum4 = sum4 + (fabs(charge(i))*
4534  molality[j]*molality[k]*CMX_P[counterIJ2]);
4535  }
4536  }
4537  }
4538 
4539  /*
4540  * for Anions, do the neutral species interaction
4541  */
4542  if (charge(j) == 0) {
4543  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i);
4544  /*
4545  * Zeta interaction term
4546  */
4547  for (size_t k = 1; k < m_kk; k++) {
4548  if (charge(k) < 0.0) {
4549  size_t izeta = j;
4550  size_t jzeta = i;
4551  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
4552  double zeta_P = psi_ijk_P[n];
4553  if (zeta_P != 0.0) {
4554  sum5 = sum5 + molality[j]*molality[k]*zeta_P;
4555  }
4556  }
4557  }
4558  }
4559  }
4560 
4561  /*
4562  * Add all of the contributions up to yield the log of the
4563  * solute activity coefficients (molality scale)
4564  */
4565  m_dlnActCoeffMolaldP_Unscaled[i] =
4566  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
4567 
4568  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4569  std::string sni = speciesName(i);
4570  printf(" %-16s lngamma[i]=%10.6f \n",
4571  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
4572  printf(" %12g %12g %12g %12g %12g %12g\n",
4573  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
4574  }
4575  }
4576 
4577 
4578  /*
4579  * ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR ANIONS ------
4580  */
4581  if (charge(i) < 0) {
4582  // species i is an anion (negative)
4583  double zsqdFdP = charge(i)*charge(i)*dFdP;
4584  double sum1 = 0.0;
4585  double sum2 = 0.0;
4586  double sum3 = 0.0;
4587  double sum4 = 0.0;
4588  double sum5 = 0.0;
4589  for (size_t j = 1; j < m_kk; j++) {
4590  /*
4591  * Find the counterIJ for the symmetric binary interaction
4592  */
4593  size_t n = m_kk*i + j;
4594  size_t counterIJ = m_CounterIJ[n];
4595 
4596  /*
4597  * For Anions, do the cation interactions.
4598  */
4599  if (charge(j) > 0) {
4600  sum1 = sum1 + molality[j]*
4601  (2.0*BMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
4602  if (j < m_kk-1) {
4603  for (size_t k = j+1; k < m_kk; k++) {
4604  // an inner sum over all cations
4605  if (charge(k) > 0) {
4606  n = k + j * m_kk + i * m_kk * m_kk;
4607  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
4608  }
4609  }
4610  }
4611  }
4612 
4613  /*
4614  * For Anions, do the other anion interactions.
4615  */
4616  if (charge(j) < 0.0) {
4617  // sum over all anions
4618  if (j != i) {
4619  sum2 = sum2 + molality[j]*(2.0*Phi_P[counterIJ]);
4620  }
4621  for (size_t k = 1; k < m_kk; k++) {
4622  if (charge(k) > 0.0) {
4623  // two inner sums over cations
4624  n = k + j * m_kk + i * m_kk * m_kk;
4625  sum2 = sum2 + molality[j]*molality[k]*psi_ijk_P[n];
4626  /*
4627  * Find the counterIJ for the symmetric binary interaction
4628  */
4629  n = m_kk*j + k;
4630  size_t counterIJ2 = m_CounterIJ[n];
4631  sum4 = sum4 +
4632  (fabs(charge(i))*
4633  molality[j]*molality[k]*CMX_P[counterIJ2]);
4634  }
4635  }
4636  }
4637 
4638  /*
4639  * for Anions, do the neutral species interaction
4640  */
4641  if (charge(j) == 0.0) {
4642  sum5 = sum5 + molality[j]*2.0*m_Lambda_nj_P(j,i);
4643  /*
4644  * Zeta interaction term
4645  */
4646  for (size_t k = 1; k < m_kk; k++) {
4647  if (charge(k) > 0.0) {
4648  size_t izeta = j;
4649  size_t jzeta = k;
4650  size_t kzeta = i;
4651  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
4652  double zeta_P = psi_ijk_P[n];
4653  if (zeta_P != 0.0) {
4654  sum5 = sum5 + molality[j]*molality[k]*zeta_P;
4655  }
4656  }
4657  }
4658  }
4659  }
4660  m_dlnActCoeffMolaldP_Unscaled[i] =
4661  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
4662  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4663  std::string sni = speciesName(i);
4664  printf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n",
4665  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
4666  printf(" %12g %12g %12g %12g %12g %12g\n",
4667  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
4668  }
4669  }
4670 
4671  /*
4672  * ------ SUBSECTION FOR CALCULATING d NEUTRAL SOLUTE ACT COEFF dP -------
4673  */
4674  if (charge(i) == 0.0) {
4675  double sum1 = 0.0;
4676  double sum3 = 0.0;
4677  for (size_t j = 1; j < m_kk; j++) {
4678  sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j);
4679  /*
4680  * Zeta term -> we piggyback on the psi term
4681  */
4682  if (charge(j) > 0.0) {
4683  for (size_t k = 1; k < m_kk; k++) {
4684  if (charge(k) < 0.0) {
4685  size_t n = k + j * m_kk + i * m_kk * m_kk;
4686  sum3 = sum3 + molality[j]*molality[k]*psi_ijk_P[n];
4687  }
4688  }
4689  }
4690  }
4691  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
4692  m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
4693  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4694  std::string sni = speciesName(i);
4695  printf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n",
4696  sni.c_str(), m_dlnActCoeffMolaldP_Unscaled[i]);
4697  }
4698  }
4699 
4700  }
4701  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4702  printf(" Step 9: \n");
4703  }
4704 
4705  /*
4706  * ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
4707  */
4708  double sum1 = 0.0;
4709  double sum2 = 0.0;
4710  double sum3 = 0.0;
4711  double sum4 = 0.0;
4712  double sum5 = 0.0;
4713  double sum6 = 0.0;
4714  double sum7 = 0.0;
4715  /*
4716  * term1 is the temperature derivative of the
4717  * DH term in the osmotic coefficient expression
4718  * b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
4719  * implementations.
4720  * Is = Ionic strength on the molality scale (units of (gmol/kg))
4721  * Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
4722  */
4723  double term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
4724 
4725  for (size_t j = 1; j < m_kk; j++) {
4726  /*
4727  * Loop Over Cations
4728  */
4729  if (charge(j) > 0.0) {
4730  for (size_t k = 1; k < m_kk; k++) {
4731  if (charge(k) < 0.0) {
4732  /*
4733  * Find the counterIJ for the symmetric j,k binary interaction
4734  */
4735  size_t n = m_kk*j + k;
4736  size_t counterIJ = m_CounterIJ[n];
4737 
4738  sum1 = sum1 + molality[j]*molality[k]*
4739  (BphiMX_P[counterIJ] + molarcharge*CMX_P[counterIJ]);
4740  }
4741  }
4742 
4743  for (size_t k = j+1; k < m_kk; k++) {
4744  if (j == (m_kk-1)) {
4745  // we should never reach this step
4746  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
4747  "logic error 1 in Step 9 of hmw_act");
4748  }
4749  if (charge(k) > 0.0) {
4750  /*
4751  * Find the counterIJ for the symmetric j,k binary interaction
4752  * between 2 cations.
4753  */
4754  size_t n = m_kk*j + k;
4755  size_t counterIJ = m_CounterIJ[n];
4756  sum2 = sum2 + molality[j]*molality[k]*Phiphi_P[counterIJ];
4757  for (size_t m = 1; m < m_kk; m++) {
4758  if (charge(m) < 0.0) {
4759  // species m is an anion
4760  n = m + k * m_kk + j * m_kk * m_kk;
4761  sum2 = sum2 +
4762  molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
4763  }
4764  }
4765  }
4766  }
4767  }
4768 
4769 
4770  /*
4771  * Loop Over Anions
4772  */
4773  if (charge(j) < 0) {
4774  for (size_t k = j+1; k < m_kk; k++) {
4775  if (j == m_kk-1) {
4776  // we should never reach this step
4777  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
4778  "logic error 2 in Step 9 of hmw_act");
4779  }
4780  if (charge(k) < 0) {
4781  /*
4782  * Find the counterIJ for the symmetric j,k binary interaction
4783  * between two anions
4784  */
4785  size_t n = m_kk*j + k;
4786  size_t counterIJ = m_CounterIJ[n];
4787 
4788  sum3 = sum3 + molality[j]*molality[k]*Phiphi_P[counterIJ];
4789  for (size_t m = 1; m < m_kk; m++) {
4790  if (charge(m) > 0.0) {
4791  n = m + k * m_kk + j * m_kk * m_kk;
4792  sum3 = sum3 +
4793  molality[j]*molality[k]*molality[m]*psi_ijk_P[n];
4794  }
4795  }
4796  }
4797  }
4798  }
4799 
4800  /*
4801  * Loop Over Neutral Species
4802  */
4803  if (charge(j) == 0) {
4804  for (size_t k = 1; k < m_kk; k++) {
4805  if (charge(k) < 0.0) {
4806  sum4 = sum4 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
4807  }
4808  if (charge(k) > 0.0) {
4809  sum5 = sum5 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
4810  }
4811  if (charge(k) == 0.0) {
4812  if (k > j) {
4813  sum6 = sum6 + molality[j]*molality[k]*m_Lambda_nj_P(j,k);
4814  } else if (k == j) {
4815  sum6 = sum6 + 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k);
4816  }
4817  }
4818  if (charge(k) < 0.0) {
4819  size_t izeta = j;
4820  for (size_t m = 1; m < m_kk; m++) {
4821  if (charge(m) > 0.0) {
4822  size_t jzeta = m;
4823  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
4824  double zeta_P = psi_ijk_P[n];
4825  if (zeta_P != 0.0) {
4826  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
4827  }
4828  }
4829  }
4830  }
4831  }
4832 
4833  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j];
4834  }
4835  }
4836  double sum_m_phi_minus_1 = 2.0 *
4837  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
4838 
4839  /*
4840  * Calculate the osmotic coefficient from
4841  * osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
4842  */
4843  double d_osmotic_coef_dP;
4844  if (molalitysum > 1.0E-150) {
4845  d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
4846  } else {
4847  d_osmotic_coef_dP = 0.0;
4848  }
4849  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4850  printf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
4851  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
4852  term1, sum1, sum2, sum3, sum4, sum5);
4853  printf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n",
4854  sum_m_phi_minus_1, d_osmotic_coef_dP);
4855  printf(" Step 10: \n");
4856  }
4857  double d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
4858 
4859 
4860  /*
4861  * In Cantera, we define the activity coefficient of the solvent as
4862  *
4863  * act_0 = actcoeff_0 * Xmol_0
4864  *
4865  * We have just computed act_0. However, this routine returns
4866  * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
4867  */
4868  m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
4869  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4870  double d_wateract_dP = exp(d_lnwateract_dP);
4871  printf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n",
4872  d_lnwateract_dP, d_wateract_dP);
4873  }
4874 }
4875 
4876 void HMWSoln::calc_lambdas(double is) const
4877 {
4878  if( m_last_is == is ) {
4879  return;
4880  }
4881  m_last_is = is;
4882 
4883  /*
4884  * Coefficients c1-c4 are used to approximate
4885  * the integral function "J";
4886  * aphi is the Debye-Huckel constant at 25 C
4887  */
4888  double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
4889  double aphi = 0.392; /* Value at 25 C */
4890  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4891  printf(" Is = %g\n", is);
4892  }
4893  if (is < 1.0E-150) {
4894  for (int i = 0; i < 17; i++) {
4895  elambda[i] = 0.0;
4896  elambda1[i] = 0.0;
4897  }
4898  return;
4899  }
4900  /*
4901  * Calculate E-lambda terms for charge combinations of like sign,
4902  * using method of Pitzer (1975). Charges up to 4 are calculated.
4903  */
4904 
4905  for (int i=1; i<=4; i++) {
4906  for (int j=i; j<=4; j++) {
4907  int ij = i*j;
4908  /*
4909  * calculate the product of the charges
4910  */
4911  double zprod = (double)ij;
4912  /*
4913  * calculate Xmn (A1) from Harvie, Weare (1980).
4914  */
4915  double x = 6.0* zprod * aphi * sqrt(is); /* eqn 23 */
4916 
4917  double jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); /* eqn 47 */
4918 
4919  double t = c3 * c4 * pow(x,c4);
4920  double dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
4921  double jprime = (jfunc/x)*(1.0 + jfunc*dj);
4922 
4923  elambda[ij] = zprod*jfunc / (4.0*is); /* eqn 14 */
4924  elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
4925  - elambda[ij])/is;
4926  if (DEBUG_MODE_ENABLED && m_debugCalc) {
4927  printf(" ij = %d, elambda = %g, elambda1 = %g\n",
4928  ij, elambda[ij], elambda1[ij]);
4929  }
4930  }
4931  }
4932 }
4933 
4934 void HMWSoln::calc_thetas(int z1, int z2,
4935  double* etheta, double* etheta_prime) const
4936 {
4937  /*
4938  * Calculate E-theta(i) and E-theta'(I) using method of
4939  * Pitzer (1987)
4940  */
4941  int i = abs(z1);
4942  int j = abs(z2);
4943 
4944  AssertThrowMsg(i <= 4 && j <= 4, "HMWSoln::calc_thetas",
4945  "we shouldn't be here");
4946  AssertThrowMsg(i != 0 && j != 0, "HMWSoln::calc_thetas",
4947  "called with one species being neutral");
4948 
4949  /*
4950  * Check to see if the charges are of opposite sign. If they are of
4951  * opposite sign then their etheta interaction is zero.
4952  */
4953  if (z1*z2 < 0) {
4954  *etheta = 0.0;
4955  *etheta_prime = 0.0;
4956  }
4957  /*
4958  * Actually calculate the interaction.
4959  */
4960  else {
4961  double f1 = (double)i / (2.0 * j);
4962  double f2 = (double)j / (2.0 * i);
4963  *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i];
4964  *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
4965  }
4966 }
4967 
4968 void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
4969 {
4970  /*
4971  * Calculate the molalities. Currently, the molalities
4972  * may not be current with respect to the contents of the
4973  * State objects' data.
4974  */
4975  calcMolalities();
4976 
4977  double xmolSolvent = moleFraction(m_indexSolvent);
4978  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
4979  if (IMS_typeCutoff_ == 0) {
4980  for (size_t k = 1; k < m_kk; k++) {
4981  IMS_lnActCoeffMolal_[k]= 0.0;
4982  }
4983  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
4984  return;
4985  } else if (IMS_typeCutoff_ == 1) {
4986  if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
4987  for (size_t k = 1; k < m_kk; k++) {
4988  IMS_lnActCoeffMolal_[k]= 0.0;
4989  }
4990  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
4991  return;
4992  } else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
4993  double tmp = log(xx * IMS_gamma_k_min_);
4994  for (size_t k = 1; k < m_kk; k++) {
4995  IMS_lnActCoeffMolal_[k]= tmp;
4996  }
4997  IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
4998  return;
4999  } else {
5000  /*
5001  * If we are in the middle region, calculate the connecting polynomials
5002  */
5003  double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
5004  double xminus2 = xminus * xminus;
5005  double xminus3 = xminus2 * xminus;
5006  double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
5007  double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
5008 
5009  double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
5010  double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
5011 
5012  double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
5013  double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
5014 
5015  double h1_g = h1 / IMS_gamma_o_min_;
5016  double h1_g_prime = h1_prime / IMS_gamma_o_min_;
5017 
5018  double alpha = 1.0 / (exp(1.0) * IMS_gamma_k_min_);
5019  double h1_f = h1 * alpha;
5020  double h1_f_prime = h1_prime * alpha;
5021 
5022  double f = h2 + h1_f;
5023  double f_prime = h2_prime + h1_f_prime;
5024 
5025  double g = h2 + h1_g;
5026  double g_prime = h2_prime + h1_g_prime;
5027 
5028  double tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
5029  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
5030  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
5031 
5032  tmp = log(xmolSolvent) + lngammak;
5033  for (size_t k = 1; k < m_kk; k++) {
5034  IMS_lnActCoeffMolal_[k]= tmp;
5035  }
5036  IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
5037  }
5038  }
5039  // Exponentials - trial 2
5040  else if (IMS_typeCutoff_ == 2) {
5041  if (xmolSolvent > IMS_X_o_cutoff_) {
5042  for (size_t k = 1; k < m_kk; k++) {
5043  IMS_lnActCoeffMolal_[k]= 0.0;
5044  }
5045  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
5046  return;
5047  } else {
5048 
5049  double xoverc = xmolSolvent/IMS_cCut_;
5050  double eterm = std::exp(-xoverc);
5051 
5052  double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
5053  + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
5054  double f_prime = 1.0 + eterm*fptmp;
5055  double f = xmolSolvent + IMS_efCut_
5056  + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
5057 
5058  double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
5059  + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
5060  double g_prime = 1.0 + eterm*gptmp;
5061  double g = xmolSolvent + IMS_egCut_
5062  + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
5063 
5064  double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
5065  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
5066  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
5067 
5068  tmp = log(xx) + lngammak;
5069  for (size_t k = 1; k < m_kk; k++) {
5070  IMS_lnActCoeffMolal_[k]= tmp;
5071  }
5072  IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
5073  }
5074  }
5075  return;
5076 }
5077 
5078 void HMWSoln::printCoeffs() const
5079 {
5080  calcMolalities();
5081  double* molality = DATA_PTR(m_molalitiesCropped);
5082  double* moleF = DATA_PTR(m_tmpV);
5083  /*
5084  * Update the coefficients wrt Temperature
5085  * Calculate the derivatives as well
5086  */
5087  s_updatePitzer_CoeffWRTemp(2);
5088  getMoleFractions(moleF);
5089 
5090  printf("Index Name MoleF MolalityCropped Charge\n");
5091  for (size_t k = 0; k < m_kk; k++) {
5092  std::string sni = speciesName(k);
5093  printf("%2s %-16s %14.7le %14.7le %5.1f \n",
5094  int2str(k).c_str(), sni.c_str(), moleF[k], molality[k], charge(k));
5095  }
5096 
5097  printf("\n Species Species beta0MX "
5098  "beta1MX beta2MX CphiMX alphaMX thetaij \n");
5099  for (size_t i = 1; i < m_kk - 1; i++) {
5100  std::string sni = speciesName(i);
5101  for (size_t j = i+1; j < m_kk; j++) {
5102  std::string snj = speciesName(j);
5103  size_t n = i * m_kk + j;
5104  size_t ct = m_CounterIJ[n];
5105  printf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
5106  sni.c_str(), snj.c_str(),
5107  m_Beta0MX_ij[ct], m_Beta1MX_ij[ct],
5108  m_Beta2MX_ij[ct], m_CphiMX_ij[ct],
5109  m_Alpha1MX_ij[ct], m_Theta_ij[ct]);
5110 
5111 
5112  }
5113  }
5114 
5115  printf("\n Species Species Species "
5116  "psi \n");
5117  for (size_t i = 1; i < m_kk; i++) {
5118  std::string sni = speciesName(i);
5119  for (size_t j = 1; j < m_kk; j++) {
5120  std::string snj = speciesName(j);
5121  for (size_t k = 1; k < m_kk; k++) {
5122  std::string snk = speciesName(k);
5123  size_t n = k + j * m_kk + i * m_kk * m_kk;
5124  if (m_Psi_ijk[n] != 0.0) {
5125  printf(" %-16s %-16s %-16s %9.5f \n",
5126  sni.c_str(), snj.c_str(),
5127  snk.c_str(), m_Psi_ijk[n]);
5128  }
5129  }
5130  }
5131  }
5132 }
5133 
5134 void HMWSoln::applyphScale(doublereal* acMolality) const
5135 {
5136  if (m_pHScalingType == PHSCALE_PITZER) {
5137  return;
5138  }
5139  AssertTrace(m_pHScalingType == PHSCALE_NBS);
5140  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
5141  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
5142  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
5143  for (size_t k = 0; k < m_kk; k++) {
5144  acMolality[k] *= exp(charge(k) * afac);
5145  }
5146 }
5147 
5148 void HMWSoln::s_updateScaling_pHScaling() const
5149 {
5150  if (m_pHScalingType == PHSCALE_PITZER) {
5151  m_lnActCoeffMolal_Scaled = m_lnActCoeffMolal_Unscaled;
5152  return;
5153  }
5154  AssertTrace(m_pHScalingType == PHSCALE_NBS);
5155  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
5156  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
5157  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
5158  for (size_t k = 0; k < m_kk; k++) {
5159  m_lnActCoeffMolal_Scaled[k] = m_lnActCoeffMolal_Unscaled[k] + charge(k) * afac;
5160  }
5161 }
5162 
5163 void HMWSoln::s_updateScaling_pHScaling_dT() const
5164 {
5165  if (m_pHScalingType == PHSCALE_PITZER) {
5166  m_dlnActCoeffMolaldT_Scaled = m_dlnActCoeffMolaldT_Unscaled;
5167  return;
5168  }
5169  AssertTrace(m_pHScalingType == PHSCALE_NBS);
5170  doublereal dlnGammaClM_dT_s2 = s_NBS_CLM_dlnMolalityActCoeff_dT();
5171  doublereal dlnGammaCLM_dT_s1 = m_dlnActCoeffMolaldT_Unscaled[m_indexCLM];
5172  doublereal afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
5173  for (size_t k = 0; k < m_kk; k++) {
5174  m_dlnActCoeffMolaldT_Scaled[k] = m_dlnActCoeffMolaldT_Unscaled[k] + charge(k) * afac;
5175  }
5176 }
5177 
5178 void HMWSoln::s_updateScaling_pHScaling_dT2() const
5179 {
5180  if (m_pHScalingType == PHSCALE_PITZER) {
5181  m_d2lnActCoeffMolaldT2_Scaled = m_d2lnActCoeffMolaldT2_Unscaled;
5182  return;
5183  }
5184  AssertTrace(m_pHScalingType == PHSCALE_NBS);
5185  doublereal d2lnGammaClM_dT2_s2 = s_NBS_CLM_d2lnMolalityActCoeff_dT2();
5186  doublereal d2lnGammaCLM_dT2_s1 = m_d2lnActCoeffMolaldT2_Unscaled[m_indexCLM];
5187  doublereal afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
5188  for (size_t k = 0; k < m_kk; k++) {
5189  m_d2lnActCoeffMolaldT2_Scaled[k] = m_d2lnActCoeffMolaldT2_Unscaled[k] + charge(k) * afac;
5190  }
5191 }
5192 
5193 void HMWSoln::s_updateScaling_pHScaling_dP() const
5194 {
5195  if (m_pHScalingType == PHSCALE_PITZER) {
5196  m_dlnActCoeffMolaldP_Scaled = m_dlnActCoeffMolaldP_Unscaled;
5197  return;
5198  }
5199  AssertTrace(m_pHScalingType == PHSCALE_NBS);
5200  doublereal dlnGammaClM_dP_s2 = s_NBS_CLM_dlnMolalityActCoeff_dP();
5201  doublereal dlnGammaCLM_dP_s1 = m_dlnActCoeffMolaldP_Unscaled[m_indexCLM];
5202  doublereal afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
5203  for (size_t k = 0; k < m_kk; k++) {
5204  m_dlnActCoeffMolaldP_Scaled[k] = m_dlnActCoeffMolaldP_Unscaled[k] + charge(k) * afac;
5205  }
5206 }
5207 
5208 doublereal HMWSoln::s_NBS_CLM_lnMolalityActCoeff() const
5209 {
5210  doublereal sqrtIs = sqrt(m_IionicMolality);
5211  doublereal A = A_Debye_TP();
5212  doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
5213  return lnGammaClMs2;
5214 }
5215 
5216 doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dT() const
5217 {
5218  doublereal sqrtIs = sqrt(m_IionicMolality);
5219  doublereal dAdT = dA_DebyedT_TP();
5220  return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
5221 }
5222 
5223 doublereal HMWSoln::s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
5224 {
5225  doublereal sqrtIs = sqrt(m_IionicMolality);
5226  doublereal d2AdT2 = d2A_DebyedT2_TP();
5227  return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
5228 }
5229 
5230 doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP() const
5231 {
5232  doublereal sqrtIs = sqrt(m_IionicMolality);
5233  doublereal dAdP = dA_DebyedP_TP();
5234  return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
5235 }
5236 
5237 int HMWSoln::debugPrinting()
5238 {
5239 #ifdef DEBUG_MODE
5240  return m_debugCalc;
5241 #else
5242  return 0;
5243 #endif
5244 }
5245 
5246 }
Cantera::HMWSoln::m_Theta_ij_P
vector_fp m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P.
Definition: HMWSoln.h:2491
Cantera::HMWSoln::m_BMX_IJ_L
vector_fp m_BMX_IJ_L
Derivative of BMX_IJ wrt T.
Definition: HMWSoln.h:2763
Cantera::HMWSoln::m_PhiPhi_IJ_P
vector_fp m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P.
Definition: HMWSoln.h:2877
Cantera::HMWSoln::s_updateIMS_lnMolalityActCoeff
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition: HMWSoln.cpp:4968
Cantera::HMWSoln::m_CphiMX_ij_LL
vector_fp m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT.
Definition: HMWSoln.h:2444
Cantera::HMWSoln::s_update_dlnMolalityActCoeff_dP
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Definition: HMWSoln.cpp:4109
Cantera::HMWSoln::pressure
virtual doublereal pressure() const
Pressure.
Definition: HMWSoln.cpp:517
Cantera::Phase::density
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:608
Cantera::HMWSoln::m_Theta_ij_L
vector_fp m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T.
Definition: HMWSoln.h:2479
Cantera::int2str
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
Cantera::HMWSoln::m_speciesCharge_Stoich
vector_fp m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
Definition: HMWSoln.h:2283
Cantera::HMWSoln::m_waterProps
WaterProps * m_waterProps
Pointer to the water property calculator.
Definition: HMWSoln.h:2263
Cantera::HMWSoln::m_CMX_IJ_P
vector_fp m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P.
Definition: HMWSoln.h:2901
Cantera::HMWSoln::m_gfunc_IJ
vector_fp m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
Definition: HMWSoln.h:2730
Cantera::ValueCache::getId
int getId()
Get a unique id for a cached value.
Definition: ValueCache.cpp:18
Cantera::HMWSoln::m_Beta2MX_ij_L
vector_fp m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T.
Definition: HMWSoln.h:2375
Cantera::findXMLPhase
XML_Node * findXMLPhase(XML_Node *root, const std::string &idtarget)
Search an XML_Node tree for a named phase XML_Node.
Definition: xml.cpp:1108
Cantera::HMWSoln::m_d2lnActCoeffMolaldT2_Unscaled
vector_fp m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2678
Cantera::HMWSoln::s_NBS_CLM_lnMolalityActCoeff
doublereal s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:5208
Cantera::HMWSoln::MC_X_o_cutoff_
doublereal MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:2959
Cantera::HMWSoln::IMS_lnActCoeffMolal_
vector_fp IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
Definition: HMWSoln.h:2913
Cantera::Array2D::resize
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.h:121
Cantera::HMWSoln::m_Lambda_nj_coeff
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2591
Cantera::HMWSoln::m_maxIionicStrength
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: HMWSoln.h:2184
Cantera::HMWSoln::m_Beta1MX_ij_LL
vector_fp m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT.
Definition: HMWSoln.h:2343
Cantera::MolalityVPSSTP::m_Mnaught
doublereal m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
Definition: MolalityVPSSTP.h:798
Cantera::HMWSoln::m_Mu_nnn_L
vector_fp m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2609
Cantera::HMWSoln::cv_mole
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition: HMWSoln.cpp:503
Cantera::HMWSoln::m_dlnActCoeffMolaldP_Unscaled
vector_fp m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2692
Cantera::HMWSoln::m_Beta0MX_ij_coeff
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer's papers.
Definition: HMWSoln.h:2320
Cantera::HMWSoln::ADebye_J
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_J.
Definition: HMWSoln.cpp:958
Cantera::VPStandardStateTP::getStandardVolumes
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
Definition: VPStandardStateTP.cpp:189
Cantera::HMWSoln::m_Beta2MX_ij_P
vector_fp m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P.
Definition: HMWSoln.h:2387
Cantera::HMWSoln::getPartialMolarCp
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
Definition: HMWSoln.cpp:806
Cantera::HMWSoln::m_Beta2MX_ij
vector_fp m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2369
Cantera::HMWSoln::m_Lambda_nj_P
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
Definition: HMWSoln.h:2573
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
Cantera::HMWSoln::m_PhiPhi_IJ
vector_fp m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer's paper.
Definition: HMWSoln.h:2859
Cantera::HMWSoln::m_Phi_IJ_L
vector_fp m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T.
Definition: HMWSoln.h:2835
electrolytes.h
Header file for a common definitions used in electrolytes thermodynamics.
Cantera::VPStandardStateTP::getCp_R
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
Definition: VPStandardStateTP.cpp:183
Cantera::HMWSoln::IMS_typeCutoff_
int IMS_typeCutoff_
IMS Cutoff type.
Definition: HMWSoln.h:2916
ThermoFactory.h
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
Cantera::HMWSoln::m_dlnActCoeffMolaldT_Scaled
vector_fp m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2657
Cantera::PDSS::satPressure
virtual doublereal satPressure(doublereal T)
saturation pressure
Definition: PDSS.cpp:372
Cantera::HMWSoln::setMolarDensity
void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) for the phase.
Definition: HMWSoln.cpp:562
Cantera::MolalityVPSSTP::m_indexCLM
size_t m_indexCLM
Index of the phScale species.
Definition: MolalityVPSSTP.h:777
Cantera::HMWSoln::m_BMX_IJ_P
vector_fp m_BMX_IJ_P
Derivative of BMX_IJ wrt P.
Definition: HMWSoln.h:2775
Cantera::HMWSoln::m_Beta0MX_ij_L
vector_fp m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T.
Definition: HMWSoln.h:2300
Cantera::HMWSoln::elambda
double elambda[17]
This is elambda, MEC.
Definition: HMWSoln.h:2714
Cantera::XML_Node
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
PDSS_Water.h
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Cantera::HMWSoln::m_dlnActCoeffMolaldP_Scaled
vector_fp m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2685
Cantera::HMWSoln::debugPrinting
int debugPrinting()
Return int specifying the amount of debug printing.
Definition: HMWSoln.cpp:5237
Cantera::HMWSoln::m_Mu_nnn_P
vector_fp m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2627
Cantera::HMWSoln::s_updateScaling_pHScaling
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
Definition: HMWSoln.cpp:5148
Cantera::VPStandardStateTP::updateStandardStateThermo
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
Definition: VPStandardStateTP.cpp:370
Cantera::HMWSoln::m_Lambda_nj_L
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
Definition: HMWSoln.h:2567
Cantera::HMWSoln::setDensity
void setDensity(const doublereal rho)
Set the internally stored density (kg/m^3) of the phase.
Definition: HMWSoln.cpp:552
Cantera::HMWSoln::m_BprimeMX_IJ_P
vector_fp m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P.
Definition: HMWSoln.h:2799
Cantera::HMWSoln::AionicRadius
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
Definition: HMWSoln.cpp:997
Cantera::HMWSoln::m_CphiMX_ij_P
vector_fp m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P.
Definition: HMWSoln.h:2450
Cantera::warn_deprecated
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
Cantera::Phase::m_speciesSize
vector_fp m_speciesSize
Vector of species sizes.
Definition: Phase.h:856
Cantera::HMWSoln::s_NBS_CLM_d2lnMolalityActCoeff_dT2
doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:5223
Cantera::HMWSoln::relative_enthalpy
virtual doublereal relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
Definition: HMWSoln.cpp:432
Cantera::HMWSoln::m_Psi_ijk_L
vector_fp m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Definition: HMWSoln.h:2524
Cantera::Phase::stateMFNumber
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition: Phase.h:823
Cantera::HMWSoln::s_update_d2lnMolalityActCoeff_dT2
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:3338
Cantera::HMWSoln::m_lnActCoeffMolal_Unscaled
vector_fp m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2650
Cantera::Array2D::ptrColumn
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:358
Cantera::HMWSoln::m_Theta_ij
vector_fp m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2473
Cantera::HMWSoln::m_hfunc_IJ
vector_fp m_hfunc_IJ
hfunc, was called gprime in Pitzer's paper.
Definition: HMWSoln.h:2743
Cantera::HMWSoln::elambda1
double elambda1[17]
This is elambda1, MEC.
Definition: HMWSoln.h:2717
Cantera::HMWSoln::relative_molal_enthalpy
virtual doublereal relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
Definition: HMWSoln.cpp:445
Cantera::HMWSoln::m_Beta1MX_ij_P
vector_fp m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P.
Definition: HMWSoln.h:2349
Cantera::HMWSoln::m_BprimeMX_IJ_LL
vector_fp m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT.
Definition: HMWSoln.h:2793
Cantera::HMWSoln::m_PhiPhi_IJ_LL
vector_fp m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT.
Definition: HMWSoln.h:2871
Cantera::HMWSoln::duplMyselfAsThermoPhase
ThermoPhase * duplMyselfAsThermoPhase() const
Duplicator from the ThermoPhase parent class.
Definition: HMWSoln.cpp:398
Cantera::HMWSoln::IMS_cCut_
doublereal IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
Definition: HMWSoln.h:2929
Cantera::HMWSoln::setTemperature
virtual void setTemperature(const doublereal temp)
Set the temperature (K)
Definition: HMWSoln.cpp:568
Cantera::Phase::getMoleFractions
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:556
Cantera::HMWSoln::m_Psi_ijk_LL
vector_fp m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
Definition: HMWSoln.h:2530
Cantera::HMWSoln::m_IionicMolality
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
Definition: HMWSoln.h:2178
Cantera::HMWSoln::m_Phi_IJ
vector_fp m_Phi_IJ
Intermediate variable called Phi in Pitzer's paper.
Definition: HMWSoln.h:2829
Cantera::HMWSoln::m_BphiMX_IJ
vector_fp m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer's paper.
Definition: HMWSoln.h:2805
Cantera::HMWSoln::m_CphiMX_ij_coeff
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
Definition: HMWSoln.h:2459
Cantera::HMWSoln::m_Beta1MX_ij
vector_fp m_Beta1MX_ij
Definition: HMWSoln.h:2331
Cantera::HMWSoln::s_updateScaling_pHScaling_dP
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure...
Definition: HMWSoln.cpp:5193
AssertThrowMsg
#define AssertThrowMsg(expr, procedure, message)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:301
Cantera::HMWSoln::m_tmpV
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition: HMWSoln.h:2269
Cantera::HMWSoln::standardConcentration
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Definition: HMWSoln.cpp:617
Cantera::ThermoPhase::thermalExpansionCoeff
virtual doublereal thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
Definition: ThermoPhase.h:321
Cantera::ThermoPhase
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:97
Cantera::HMWSoln::getUnscaledMolalityActivityCoefficients
void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
Definition: HMWSoln.cpp:674
Cantera::HMWSoln::HMWSoln
HMWSoln()
Default Constructor.
Definition: HMWSoln.cpp:31
Cantera::HMWSoln::m_lnActCoeffMolal_Scaled
vector_fp m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2640
Cantera::HMWSoln::getActivities
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:653
Cantera::Phase::mean_X
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:687
Cantera::HMWSoln::MC_X_o_min_
doublereal MC_X_o_min_
gamma_o value for the cutoff process at the zero solvent point
Definition: HMWSoln.h:2962
Cantera::HMWSoln::m_Phi_IJ_P
vector_fp m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P.
Definition: HMWSoln.h:2847
Cantera::HMWSoln::m_waterSS
PDSS * m_waterSS
Water standard state calculator.
Definition: HMWSoln.h:2254
Cantera::HMWSoln::m_Alpha1MX_ij
vector_fp m_Alpha1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2408
Cantera::importPhase
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty ThermoPhase object.
Definition: ThermoFactory.cpp:316
Cantera::HMWSoln::setPressure
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
Definition: HMWSoln.cpp:522
Cantera::HMWSoln::m_Phiprime_IJ
vector_fp m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer's paper.
Definition: HMWSoln.h:2853
Cantera::HMWSoln::getPartialMolarEntropies
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: HMWSoln.cpp:743
Cantera::HMWSoln::m_h2func_IJ
vector_fp m_h2func_IJ
hfunc2, was called gprime in Pitzer's paper.
Definition: HMWSoln.h:2750
Cantera::HMWSoln::m_Psi_ijk_P
vector_fp m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
Definition: HMWSoln.h:2536
Cantera::HMWSoln::counterIJ_setup
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
Definition: HMWSoln.cpp:1381
Cantera::HMWSoln::m_CMX_IJ_L
vector_fp m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T.
Definition: HMWSoln.h:2889
Cantera::HMWSoln::getPartialMolarVolumes
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
Definition: HMWSoln.cpp:788
Cantera::HMWSoln::m_molalitiesAreCropped
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
Definition: HMWSoln.h:2703
Cantera::HMWSoln::m_dlnActCoeffMolaldT_Unscaled
vector_fp m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2664
Cantera::HMWSoln::m_Psi_ijk_coeff
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2553
Cantera::HMWSoln::s_updateScaling_pHScaling_dT
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature...
Definition: HMWSoln.cpp:5163
HMWSoln.h
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
Cantera::HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients...
Definition: HMWSoln.cpp:2596
Cantera::HMWSoln::m_Mu_nnn
vector_fp m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2600
Cantera::HMWSoln::m_molalitiesCropped
vector_fp m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
Definition: HMWSoln.h:2700
Cantera::Phase::molarVolume
doublereal molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:673
Cantera::WaterProps
The WaterProps class is used to house several approximation routines for properties of water...
Definition: WaterProps.h:96
PITZERFORM_BASE
#define PITZERFORM_BASE
Major Parameters: The form of the Pitzer expression refers to the form of the Gibbs free energy expre...
Definition: HMWSoln.h:39
Cantera::HMWSoln::m_Beta2MX_ij_LL
vector_fp m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT.
Definition: HMWSoln.h:2381
Cantera::ValueCache::getScalar
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (doublereal) with the given id...
Definition: ValueCache.h:162
Cantera::HMWSoln::m_CMX_IJ_LL
vector_fp m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT.
Definition: HMWSoln.h:2895
Cantera::HMWSoln::getUnitsStandardConc
virtual void getUnitsStandardConc(double *uA, int k=0, int sizeUA=6) const
Returns the units of the standard and generalized concentrations.
Definition: HMWSoln.cpp:627
Cantera::HMWSoln::m_Beta2MX_ij_coeff
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer's papers.
Definition: HMWSoln.h:2397
Cantera::MolalityVPSSTP::operator=
MolalityVPSSTP & operator=(const MolalityVPSSTP &b)
Assignment operator.
Definition: MolalityVPSSTP.cpp:58
Cantera::HMWSoln::m_Lambda_nj_LL
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
Definition: HMWSoln.h:2570
Cantera::HMWSoln::s_updatePitzer_lnMolalityActCoeff
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
Definition: HMWSoln.cpp:1678
Cantera::HMWSoln::ADebye_V
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_V.
Definition: HMWSoln.cpp:947
Cantera::HMWSoln::IMS_slopegCut_
doublereal IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
Definition: HMWSoln.h:2943
Cantera::HMWSoln::IMS_slopefCut_
doublereal IMS_slopefCut_
Parameter in the polyExp cutoff treatment.
Definition: HMWSoln.h:2936
Cantera::HMWSoln::m_Beta0MX_ij_LL
vector_fp m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT.
Definition: HMWSoln.h:2306
Cantera::PHSCALE_NBS
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
Definition: MolalityVPSSTP.h:852
Cantera::HMWSoln::m_form_A_Debye
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
Definition: HMWSoln.h:2214
Cantera::Phase::m_cache
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition: Phase.h:833
Cantera::HMWSoln::entropy_mole
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
Definition: HMWSoln.cpp:485
Cantera::HMWSoln::s_update_dlnMolalityActCoeff_dT
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
Definition: HMWSoln.cpp:2561
Cantera::HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients...
Definition: HMWSoln.cpp:4133
Cantera::HMWSoln::m_CounterIJ
vector_int m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species...
Definition: HMWSoln.h:2711
Cantera::HMWSoln::MC_slopepCut_
doublereal MC_slopepCut_
Parameter in the Molality Exp cutoff treatment.
Definition: HMWSoln.h:2969
Cantera::HMWSoln::s_updateScaling_pHScaling_dT2
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:5178
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
Cantera::HMWSoln::m_CMX_IJ
vector_fp m_CMX_IJ
Intermediate variable called CMX in Pitzer's paper.
Definition: HMWSoln.h:2883
Cantera::HMWSoln::enthalpy_mole
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
Definition: HMWSoln.cpp:425
Cantera::HMWSoln::m_TempPitzerRef
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
Definition: HMWSoln.h:2187
Cantera::HMWSoln::m_CphiMX_ij_L
vector_fp m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T.
Definition: HMWSoln.h:2438
Cantera::HMWSoln::calc_lambdas
void calc_lambdas(double is) const
Calculate the lambda interactions.
Definition: HMWSoln.cpp:4876
Cantera::PHSCALE_PITZER
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer...
Definition: MolalityVPSSTP.h:827
Cantera::HMWSoln::calc_thetas
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
Definition: HMWSoln.cpp:4934
Cantera::HMWSoln::m_BphiMX_IJ_L
vector_fp m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T.
Definition: HMWSoln.h:2811
Cantera::HMWSoln::m_g2func_IJ
vector_fp m_g2func_IJ
This is the value of g2(x2) in Pitzer's papers.
Definition: HMWSoln.h:2736
Cantera::MolalityVPSSTP::m_pHScalingType
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
Definition: MolalityVPSSTP.h:769
Cantera::PDSS::density
virtual doublereal density() const
Return the standard state density at standard state.
Definition: PDSS.cpp:257
Cantera::HMWSoln::m_densWaterSS
double m_densWaterSS
density of standard-state water
Definition: HMWSoln.h:2260
Cantera::HMWSoln::m_Lambda_nj
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
Definition: HMWSoln.h:2564
Cantera::HMWSoln::~HMWSoln
virtual ~HMWSoln()
Destructor.
Definition: HMWSoln.cpp:392
Cantera::HMWSoln::m_Phi_IJ_LL
vector_fp m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT.
Definition: HMWSoln.h:2841
Cantera::HMWSoln::m_formPitzer
int m_formPitzer
This is the form of the Pitzer parameterization used in this model.
Definition: HMWSoln.h:2111
Cantera::HMWSoln::IMS_gamma_k_min_
doublereal IMS_gamma_k_min_
gamma_k minimum for the cutoff process at the zero solvent point
Definition: HMWSoln.h:2926
Cantera::HMWSoln::m_Alpha2MX_ij
vector_fp m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2420
Cantera::HMWSoln::A_Debye_TP
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:848
Cantera::HMWSoln::m_electrolyteSpeciesType
vector_int m_electrolyteSpeciesType
Vector containing the electrolyte species type.
Definition: HMWSoln.h:2164
Cantera::HMWSoln::m_Beta0MX_ij
vector_fp m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2294
Cantera::HMWSoln::m_BprimeMX_IJ_L
vector_fp m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T.
Definition: HMWSoln.h:2787
Cantera::cHMWSoln0
const int cHMWSoln0
eosTypes returned for this ThermoPhase Object
Definition: electrolytes.h:41
Cantera::ThermoPhase::isothermalCompressibility
virtual doublereal isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
Definition: ThermoPhase.h:310
Cantera::HMWSoln::IMS_gamma_o_min_
doublereal IMS_gamma_o_min_
gamma_o value for the cutoff process at the zero solvent point
Definition: HMWSoln.h:2923
Cantera::Phase::moleFraction
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:561
Cantera::HMWSoln::dA_DebyedT_TP
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:880
AssertTrace
#define AssertTrace(expr)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:270
Cantera::HMWSoln::m_BMX_IJ_LL
vector_fp m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT.
Definition: HMWSoln.h:2769
Cantera::VPStandardStateTP::getEntropy_R
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
Definition: VPStandardStateTP.cpp:157
Cantera::Phase::temperature
doublereal temperature() const
Temperature (K).
Definition: Phase.h:602
Cantera::Phase::nDim
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition: Phase.h:586
Cantera::HMWSoln::m_formGC
int m_formGC
Format for the generalized concentration:
Definition: HMWSoln.h:2152
Cantera::MolalityVPSSTP::m_molalities
vector_fp m_molalities
Current value of the molalities of the species in the phase.
Definition: MolalityVPSSTP.h:805
Cantera::HMWSoln::m_A_Debye
double m_A_Debye
A_Debye -> this expression appears on the top of the ln actCoeff term in the general Debye-Huckel exp...
Definition: HMWSoln.h:2248
Cantera::HMWSoln::CROP_speciesCropped_
std::vector< int > CROP_speciesCropped_
This is a boolean-type vector indicating whether a species's activity coefficient is in the cropped r...
Definition: HMWSoln.h:2992
Cantera::HMWSoln::initLengths
void initLengths()
Initialize all of the species-dependent lengths in the object.
Definition: HMWSoln.cpp:1006
Cantera::HMWSoln::getActivityConcentrations
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized activity concentrations.
Definition: HMWSoln.cpp:604
Cantera::SmallNumber
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
Cantera::HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dT
doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:5216
Cantera::HMWSoln::ADebye_L
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer's definition of A_L.
Definition: HMWSoln.cpp:936
Cantera::HMWSoln::gibbs_mole
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: HMWSoln.cpp:491
Cantera::CachedValue::validate
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
Definition: ValueCache.h:41
Cantera::Phase::setTemperature
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:638
Cantera::HMWSoln::dA_DebyedP_TP
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:904
Cantera::HMWSoln::applyphScale
virtual void applyphScale(doublereal *acMolality) const
Apply the current phScale to a set of activity Coefficients or activities.
Definition: HMWSoln.cpp:5134
Cantera::VPStandardStateTP::getEnthalpy_RT
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
Definition: VPStandardStateTP.cpp:145
Cantera::HMWSoln::calcDensity
void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
Definition: HMWSoln.cpp:527
Cantera::HMWSoln::m_gamma_tmp
vector_fp m_gamma_tmp
Intermediate storage of the activity coefficient itself.
Definition: HMWSoln.h:2907
Cantera::Phase::meanMolecularWeight
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:669
Cantera::GasConstant
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
Cantera::HMWSoln::eosType
virtual int eosType() const
Equation of state type flag.
Definition: HMWSoln.cpp:403
Cantera::HMWSoln::getChemPotentials
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: HMWSoln.cpp:689
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
DATA_PTR
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
Cantera::MolalityVPSSTP::m_weightSolvent
doublereal m_weightSolvent
Molecular weight of the Solvent.
Definition: MolalityVPSSTP.h:780
Cantera::HMWSoln::IMS_X_o_cutoff_
doublereal IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:2920
Cantera::HMWSoln::calcMolalitiesCropped
void calcMolalitiesCropped() const
Calculate the cropped molalities.
Definition: HMWSoln.cpp:1228
Cantera::HMWSoln::m_Psi_ijk
vector_fp m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2518
Cantera::HMWSoln::m_pp
vector_fp m_pp
Temporary array used in equilibrium calculations.
Definition: HMWSoln.h:2266
Cantera::HMWSoln
Class HMWSoln represents a dilute or concentrated liquid electrolyte phase which obeys the Pitzer for...
Definition: HMWSoln.h:1225
Cantera::WaterProps::ADebye
doublereal ADebye(doublereal T, doublereal P, int ifunc)
ADebye calculates the value of A_Debye as a function of temperature and pressure according to relatio...
Definition: WaterProps.cpp:199
Cantera::HMWSoln::m_BMX_IJ
vector_fp m_BMX_IJ
Intermediate variable called BMX in Pitzer's paper This is the basic cation - anion interaction...
Definition: HMWSoln.h:2757
Cantera::HMWSoln::d2A_DebyedT2_TP
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:970
Cantera::HMWSoln::setState_TP
virtual void setState_TP(doublereal t, doublereal p)
Set the temperature (K) and pressure (Pa)
Definition: HMWSoln.cpp:573
Cantera::Phase::m_kk
size_t m_kk
Number of species in the phase.
Definition: Phase.h:843
Cantera::HMWSoln::m_IionicMolalityStoich
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
Definition: HMWSoln.h:2195
Cantera::ThermoPhase::initThermoFile
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
Definition: ThermoPhase.cpp:637
Cantera::HMWSoln::m_BprimeMX_IJ
vector_fp m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer's paper.
Definition: HMWSoln.h:2781
Cantera::VPStandardStateTP::m_Pcurrent
doublereal m_Pcurrent
Current value of the pressure - state variable.
Definition: VPStandardStateTP.h:523
Cantera::HMWSoln::s_updatePitzer_CoeffWRTemp
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients' dependence on the temperature.
Definition: HMWSoln.cpp:1404
Cantera::CachedValue::value
T value
The value of the cached property.
Definition: ValueCache.h:115
Cantera::HMWSoln::m_d2lnActCoeffMolaldT2_Scaled
vector_fp m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2671
Cantera::VPStandardStateTP::getStandardChemPotentials
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
Definition: VPStandardStateTP.cpp:135
Cantera::HMWSoln::m_CphiMX_ij
vector_fp m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2432
Cantera::HMWSoln::density
virtual doublereal density() const
Returns the current value of the density.
Definition: HMWSoln.cpp:547
Cantera::MolalityVPSSTP::m_indexSolvent
size_t m_indexSolvent
Index of the solvent.
Definition: MolalityVPSSTP.h:759
Cantera::HMWSoln::m_formPitzerTemp
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model...
Definition: HMWSoln.h:2121
Cantera::HMWSoln::printCoeffs
void printCoeffs() const
Print out all of the input Pitzer coefficients.
Definition: HMWSoln.cpp:5078
Cantera::HMWSoln::operator=
HMWSoln & operator=(const HMWSoln &right)
Assignment operator.
Definition: HMWSoln.cpp:234
Cantera::HMWSoln::m_Beta1MX_ij_L
vector_fp m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T.
Definition: HMWSoln.h:2337
Cantera::HMWSoln::m_Theta_ij_LL
vector_fp m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT.
Definition: HMWSoln.h:2485
Cantera::HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:3371
Cantera::PDSS::setState_TP
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
Definition: PDSS.cpp:362
Cantera::HMWSoln::m_Theta_ij_coeff
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2505
Cantera::HMWSoln::cp_mole
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition: HMWSoln.cpp:497
Cantera::HMWSoln::satPressure
virtual doublereal satPressure(doublereal T)
Get the saturation pressure for a given temperature.
Definition: HMWSoln.cpp:837
Cantera::HMWSoln::m_Mu_nnn_coeff
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
Definition: HMWSoln.h:2630
Cantera::HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP
doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
Definition: HMWSoln.cpp:5230
Cantera::HMWSoln::m_Mu_nnn_LL
vector_fp m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2618
Cantera::MolalityVPSSTP::calcMolalities
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
Definition: MolalityVPSSTP.cpp:129
Cantera::Phase::speciesName
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:272
Cantera::HMWSoln::m_Aionic
vector_fp m_Aionic
a_k = Size of the ionic species in the DH formulation units = meters
Definition: HMWSoln.h:2170
Cantera::HMWSoln::m_BphiMX_IJ_P
vector_fp m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P.
Definition: HMWSoln.h:2823
Cantera::HMWSoln::m_Beta1MX_ij_coeff
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer's papers.
Definition: HMWSoln.h:2357
Cantera::HMWSoln::m_PhiPhi_IJ_L
vector_fp m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T.
Definition: HMWSoln.h:2865
Cantera::HMWSoln::getPartialMolarEnthalpies
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
Definition: HMWSoln.cpp:717
Cantera::Phase::setDensity
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase Note the density of a phase is an independent...
Definition: Phase.h:623
Cantera::HMWSoln::m_Beta0MX_ij_P
vector_fp m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P.
Definition: HMWSoln.h:2312
Cantera::HMWSoln::m_BphiMX_IJ_LL
vector_fp m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT.
Definition: HMWSoln.h:2817
Cantera::Phase::charge
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:578
Generated by   doxygen 1.8.6