Cantera  2.3.0
HMWSoln_input.cpp
Go to the documentation of this file.
1 /**
2  * @file HMWSoln_input.cpp
3  * Definitions for the HMWSoln ThermoPhase object, which models concentrated
4  * electrolyte solutions
5  * (see \ref thermoprops and \link Cantera::HMWSoln HMWSoln \endlink) .
6  *
7  * This file contains definitions for reading in the interaction terms
8  * in the formulation.
9  */
10 
11 // This file is part of Cantera. See License.txt in the top-level directory or
12 // at http://www.cantera.org/license.txt for license and copyright information.
13 
14 #include "cantera/thermo/HMWSoln.h"
15 #include "cantera/thermo/ThermoFactory.h"
16 #include "cantera/thermo/PDSS_Water.h"
17 #include "cantera/thermo/electrolytes.h"
18 #include "cantera/base/stringUtils.h"
19 #include "cantera/base/ctml.h"
20 
21 #include <fstream>
22 
23 using namespace std;
24 
25 namespace Cantera
26 {
27 int HMWSoln::interp_est(const std::string& estString)
28 {
29  if (ba::iequals(estString, "solvent")) {
30  return cEST_solvent;
31  } else if (ba::iequals(estString, "chargedspecies")) {
32  return cEST_chargedSpecies;
33  } else if (ba::iequals(estString, "weakacidassociated")) {
34  return cEST_weakAcidAssociated;
35  } else if (ba::iequals(estString, "strongacidassociated")) {
36  return cEST_strongAcidAssociated;
37  } else if (ba::iequals(estString, "polarneutral")) {
38  return cEST_polarNeutral;
39  } else if (ba::iequals(estString, "nonpolarneutral")) {
40  return cEST_nonpolarNeutral;
41  }
42  int retn, rval;
43  if ((retn = sscanf(estString.c_str(), "%d", &rval)) != 1) {
44  return -1;
45  }
46  return rval;
47 }
48 
49 void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
50 {
51  string xname = BinSalt.name();
52  if (xname != "binarySaltParameters") {
53  throw CanteraError("HMWSoln::readXMLBinarySalt",
54  "Incorrect name for processing this routine: " + xname);
55  }
56 
57  vector_fp vParams;
58  string iName = BinSalt.attrib("cation");
59  if (iName == "") {
60  throw CanteraError("HMWSoln::readXMLBinarySalt", "no cation attrib");
61  }
62  string jName = BinSalt.attrib("anion");
63  if (jName == "") {
64  throw CanteraError("HMWSoln::readXMLBinarySalt", "no anion attrib");
65  }
66 
67  // Find the index of the species in the current phase. It's not an error to
68  // not find the species
69  size_t iSpecies = speciesIndex(iName);
70  if (iSpecies == npos) {
71  return;
72  }
73  string ispName = speciesName(iSpecies);
74  if (charge(iSpecies) <= 0) {
75  throw CanteraError("HMWSoln::readXMLBinarySalt", "cation charge problem");
76  }
77  size_t jSpecies = speciesIndex(jName);
78  if (jSpecies == npos) {
79  return;
80  }
81  string jspName = speciesName(jSpecies);
82  if (charge(jSpecies) >= 0) {
83  throw CanteraError("HMWSoln::readXMLBinarySalt", "anion charge problem");
84  }
85 
86  size_t n = iSpecies * m_kk + jSpecies;
87  int counter = m_CounterIJ[n];
88  for (size_t iChild = 0; iChild < BinSalt.nChildren(); iChild++) {
89  XML_Node& xmlChild = BinSalt.child(iChild);
90  string nodeName = xmlChild.name();
91 
92  // Process the binary salt child elements
93  if (ba::iequals(nodeName, "beta0")) {
94  // Get the string containing all of the values
95  getFloatArray(xmlChild, vParams, false, "", "beta0");
96  size_t nParamsFound = vParams.size();
97  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
98  if (nParamsFound != 1) {
99  throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
100  + "::" + jspName,
101  "wrong number of params found");
102  }
103  m_Beta0MX_ij[counter] = vParams[0];
104  m_Beta0MX_ij_coeff(0,counter) = m_Beta0MX_ij[counter];
105  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
106  if (nParamsFound != 2) {
107  throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
108  + "::" + jspName,
109  "wrong number of params found");
110  }
111  m_Beta0MX_ij_coeff(0,counter) = vParams[0];
112  m_Beta0MX_ij_coeff(1,counter) = vParams[1];
113  m_Beta0MX_ij[counter] = vParams[0];
114  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
115  if (nParamsFound != 5) {
116  throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
117  + "::" + jspName,
118  "wrong number of params found");
119  }
120  for (size_t i = 0; i < nParamsFound; i++) {
121  m_Beta0MX_ij_coeff(i, counter) = vParams[i];
122  }
123  m_Beta0MX_ij[counter] = vParams[0];
124  }
125  }
126  if (ba::iequals(nodeName, "beta1")) {
127  // Get the string containing all of the values
128  getFloatArray(xmlChild, vParams, false, "", "beta1");
129  size_t nParamsFound = vParams.size();
130  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
131  if (nParamsFound != 1) {
132  throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
133  + "::" + jspName,
134  "wrong number of params found");
135  }
136  m_Beta1MX_ij[counter] = vParams[0];
137  m_Beta1MX_ij_coeff(0,counter) = m_Beta1MX_ij[counter];
138  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
139  if (nParamsFound != 2) {
140  throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
141  + "::" + jspName,
142  "wrong number of params found");
143  }
144  m_Beta1MX_ij_coeff(0,counter) = vParams[0];
145  m_Beta1MX_ij_coeff(1,counter) = vParams[1];
146  m_Beta1MX_ij[counter] = vParams[0];
147  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
148  if (nParamsFound != 5) {
149  throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
150  + "::" + jspName,
151  "wrong number of params found");
152  }
153  for (size_t i = 0; i < nParamsFound; i++) {
154  m_Beta1MX_ij_coeff(i, counter) = vParams[i];
155  }
156  m_Beta1MX_ij[counter] = vParams[0];
157  }
158  }
159  if (ba::iequals(nodeName, "beta2")) {
160  getFloatArray(xmlChild, vParams, false, "", "beta2");
161  size_t nParamsFound = vParams.size();
162  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
163  if (nParamsFound != 1) {
164  throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
165  + "::" + jspName,
166  "wrong number of params found");
167  }
168  m_Beta2MX_ij[counter] = vParams[0];
169  m_Beta2MX_ij_coeff(0,counter) = m_Beta2MX_ij[counter];
170  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
171  if (nParamsFound != 2) {
172  throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
173  + "::" + jspName,
174  "wrong number of params found");
175  }
176  m_Beta2MX_ij_coeff(0,counter) = vParams[0];
177  m_Beta2MX_ij_coeff(1,counter) = vParams[1];
178  m_Beta2MX_ij[counter] = vParams[0];
179  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
180  if (nParamsFound != 5) {
181  throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
182  + "::" + jspName,
183  "wrong number of params found");
184  }
185  for (size_t i = 0; i < nParamsFound; i++) {
186  m_Beta2MX_ij_coeff(i, counter) = vParams[i];
187  }
188  m_Beta2MX_ij[counter] = vParams[0];
189  }
190  }
191  if (ba::iequals(nodeName, "cphi")) {
192  // Get the string containing all of the values
193  getFloatArray(xmlChild, vParams, false, "", "Cphi");
194  size_t nParamsFound = vParams.size();
195  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
196  if (nParamsFound != 1) {
197  throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
198  + "::" + jspName,
199  "wrong number of params found");
200  }
201  m_CphiMX_ij[counter] = vParams[0];
202  m_CphiMX_ij_coeff(0,counter) = m_CphiMX_ij[counter];
203  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
204  if (nParamsFound != 2) {
205  throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
206  + "::" + jspName,
207  "wrong number of params found");
208  }
209  m_CphiMX_ij_coeff(0,counter) = vParams[0];
210  m_CphiMX_ij_coeff(1,counter) = vParams[1];
211  m_CphiMX_ij[counter] = vParams[0];
212  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
213  if (nParamsFound != 5) {
214  throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
215  + "::" + jspName,
216  "wrong number of params found");
217  }
218  for (size_t i = 0; i < nParamsFound; i++) {
219  m_CphiMX_ij_coeff(i, counter) = vParams[i];
220  }
221  m_CphiMX_ij[counter] = vParams[0];
222  }
223  }
224 
225  if (ba::iequals(nodeName, "alpha1")) {
226  m_Alpha1MX_ij[counter] = fpValueCheck(xmlChild.value());
227  }
228 
229  if (ba::iequals(nodeName, "alpha2")) {
230  m_Alpha2MX_ij[counter] = fpValueCheck(xmlChild.value());
231  }
232  }
233 }
234 
235 void HMWSoln::readXMLThetaAnion(XML_Node& BinSalt)
236 {
237  string xname = BinSalt.name();
238  vector_fp vParams;
239  if (xname != "thetaAnion") {
240  throw CanteraError("HMWSoln::readXMLThetaAnion",
241  "Incorrect name for processing this routine: " + xname);
242  }
243  string ispName = BinSalt.attrib("anion1");
244  if (ispName == "") {
245  throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion1 attrib");
246  }
247  string jspName = BinSalt.attrib("anion2");
248  if (jspName == "") {
249  throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion2 attrib");
250  }
251 
252  // Find the index of the species in the current phase. It's not an error to
253  // not find the species
254  size_t iSpecies = speciesIndex(ispName);
255  if (iSpecies == npos) {
256  return;
257  }
258  if (charge(iSpecies) >= 0) {
259  throw CanteraError("HMWSoln::readXMLThetaAnion", "anion1 charge problem");
260  }
261  size_t jSpecies = speciesIndex(jspName);
262  if (jSpecies == npos) {
263  return;
264  }
265  if (charge(jSpecies) >= 0) {
266  throw CanteraError("HMWSoln::readXMLThetaAnion", "anion2 charge problem");
267  }
268 
269  size_t n = iSpecies * m_kk + jSpecies;
270  int counter = m_CounterIJ[n];
271  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
272  XML_Node& xmlChild = BinSalt.child(i);
273  if (ba::iequals(xmlChild.name(), "theta")) {
274  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
275  size_t nParamsFound = vParams.size();
276  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
277  if (nParamsFound != 1) {
278  throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
279  + "::" + jspName,
280  "wrong number of params found");
281  }
282  m_Theta_ij_coeff(0,counter) = vParams[0];
283  m_Theta_ij[counter] = vParams[0];
284  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
285  if (nParamsFound != 2) {
286  throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
287  + "::" + jspName,
288  "wrong number of params found");
289  }
290  m_Theta_ij_coeff(0,counter) = vParams[0];
291  m_Theta_ij_coeff(1,counter) = vParams[1];
292  m_Theta_ij[counter] = vParams[0];
293  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
294  if (nParamsFound == 1) {
295  vParams.resize(5, 0.0);
296  nParamsFound = 5;
297  } else if (nParamsFound != 5) {
298  throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
299  + "::" + jspName,
300  "wrong number of params found");
301  }
302  for (size_t j = 0; j < nParamsFound; j++) {
303  m_Theta_ij_coeff(j, counter) = vParams[j];
304  }
305  m_Theta_ij[counter] = vParams[0];
306  }
307  }
308  }
309 }
310 
311 void HMWSoln::readXMLThetaCation(XML_Node& BinSalt)
312 {
313  string xname = BinSalt.name();
314  vector_fp vParams;
315  if (xname != "thetaCation") {
316  throw CanteraError("HMWSoln::readXMLThetaCation",
317  "Incorrect name for processing this routine: " + xname);
318  }
319  string ispName = BinSalt.attrib("cation1");
320  if (ispName == "") {
321  throw CanteraError("HMWSoln::readXMLThetaCation", "no cation1 attrib");
322  }
323  string jspName = BinSalt.attrib("cation2");
324  if (jspName == "") {
325  throw CanteraError("HMWSoln::readXMLThetaCation", "no cation2 attrib");
326  }
327 
328  // Find the index of the species in the current phase. It's not an error to
329  // not find the species
330  size_t iSpecies = speciesIndex(ispName);
331  if (iSpecies == npos) {
332  return;
333  }
334  if (charge(iSpecies) <= 0) {
335  throw CanteraError("HMWSoln::readXMLThetaCation", "cation1 charge problem");
336  }
337  size_t jSpecies = speciesIndex(jspName);
338  if (jSpecies == npos) {
339  return;
340  }
341  if (charge(jSpecies) <= 0) {
342  throw CanteraError("HMWSoln::readXMLThetaCation", "cation2 charge problem");
343  }
344 
345  size_t n = iSpecies * m_kk + jSpecies;
346  int counter = m_CounterIJ[n];
347  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
348  XML_Node& xmlChild = BinSalt.child(i);
349  if (ba::iequals(xmlChild.name(), "theta")) {
350  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
351  size_t nParamsFound = vParams.size();
352  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
353  if (nParamsFound != 1) {
354  throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
355  + "::" + jspName,
356  "wrong number of params found");
357  }
358  m_Theta_ij_coeff(0,counter) = vParams[0];
359  m_Theta_ij[counter] = vParams[0];
360  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
361  if (nParamsFound != 2) {
362  throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
363  + "::" + jspName,
364  "wrong number of params found");
365  }
366  m_Theta_ij_coeff(0,counter) = vParams[0];
367  m_Theta_ij_coeff(1,counter) = vParams[1];
368  m_Theta_ij[counter] = vParams[0];
369  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
370  if (nParamsFound == 1) {
371  vParams.resize(5, 0.0);
372  nParamsFound = 5;
373  } else if (nParamsFound != 5) {
374  throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
375  + "::" + jspName,
376  "wrong number of params found");
377  }
378  for (size_t j = 0; j < nParamsFound; j++) {
379  m_Theta_ij_coeff(j, counter) = vParams[j];
380  }
381  m_Theta_ij[counter] = vParams[0];
382  }
383  }
384  }
385 }
386 
387 void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
388 {
389  string xname = BinSalt.name();
390  if (xname != "psiCommonCation") {
391  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
392  "Incorrect name for processing this routine: " + xname);
393  }
394  vector_fp vParams;
395  string kName = BinSalt.attrib("cation");
396  if (kName == "") {
397  throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no cation attrib");
398  }
399  string iName = BinSalt.attrib("anion1");
400  if (iName == "") {
401  throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion1 attrib");
402  }
403  string jName = BinSalt.attrib("anion2");
404  if (jName == "") {
405  throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion2 attrib");
406  }
407 
408  // Find the index of the species in the current phase. It's not an error to
409  // not find the species
410  size_t kSpecies = speciesIndex(kName);
411  if (kSpecies == npos) {
412  return;
413  }
414  if (charge(kSpecies) <= 0) {
415  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
416  "cation charge problem");
417  }
418  size_t iSpecies = speciesIndex(iName);
419  if (iSpecies == npos) {
420  return;
421  }
422  if (charge(iSpecies) >= 0) {
423  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
424  "anion1 charge problem");
425  }
426  size_t jSpecies = speciesIndex(jName);
427  if (jSpecies == npos) {
428  return;
429  }
430  if (charge(jSpecies) >= 0) {
431  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
432  "anion2 charge problem");
433  }
434 
435  size_t n = iSpecies * m_kk + jSpecies;
436  int counter = m_CounterIJ[n];
437  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
438  XML_Node& xmlChild = BinSalt.child(i);
439  if (ba::iequals(xmlChild.name(), "theta")) {
440  double old = m_Theta_ij[counter];
441  m_Theta_ij[counter] = fpValueCheck(xmlChild.value());
442  if (old != 0.0 && old != m_Theta_ij[counter]) {
443  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
444  "conflicting values");
445  }
446  }
447  if (ba::iequals(xmlChild.name(), "psi")) {
448  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
449  size_t nParamsFound = vParams.size();
450  n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies;
451 
452  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
453  if (nParamsFound != 1) {
454  throw CanteraError("HMWSoln::readXMLPsiCommonCation::Psi for "
455  + kName + "::" + iName + "::" + jName,
456  "wrong number of params found");
457  }
458  m_Psi_ijk_coeff(0,n) = vParams[0];
459  m_Psi_ijk[n] = vParams[0];
460  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
461  if (nParamsFound != 2) {
462  throw CanteraError("HMWSoln::readXMLPsiCation::Psi for "
463  + kName + "::" + iName + "::" + jName,
464  "wrong number of params found");
465  }
466  m_Psi_ijk_coeff(0,n) = vParams[0];
467  m_Psi_ijk_coeff(1,n) = vParams[1];
468  m_Psi_ijk[n] = vParams[0];
469  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
470  if (nParamsFound == 1) {
471  vParams.resize(5, 0.0);
472  nParamsFound = 5;
473  } else if (nParamsFound != 5) {
474  throw CanteraError("HMWSoln::readXMLPsiCation::Psi for "
475  + kName + "::" + iName + "::" + jName,
476  "wrong number of params found");
477  }
478  for (size_t j = 0; j < nParamsFound; j++) {
479  m_Psi_ijk_coeff(j, n) = vParams[j];
480  }
481  m_Psi_ijk[n] = vParams[0];
482  }
483 
484  // fill in the duplicate entries
485  n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies;
486  for (size_t j = 0; j < nParamsFound; j++) {
487  m_Psi_ijk_coeff(j, n) = vParams[j];
488  }
489  m_Psi_ijk[n] = vParams[0];
490 
491  n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies;
492  for (size_t j = 0; j < nParamsFound; j++) {
493  m_Psi_ijk_coeff(j, n) = vParams[j];
494  }
495  m_Psi_ijk[n] = vParams[0];
496 
497  n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies;
498  for (size_t j = 0; j < nParamsFound; j++) {
499  m_Psi_ijk_coeff(j, n) = vParams[j];
500  }
501  m_Psi_ijk[n] = vParams[0];
502 
503  n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies;
504  for (size_t j = 0; j < nParamsFound; j++) {
505  m_Psi_ijk_coeff(j, n) = vParams[j];
506  }
507  m_Psi_ijk[n] = vParams[0];
508 
509  n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies;
510  for (size_t j = 0; j < nParamsFound; j++) {
511  m_Psi_ijk_coeff(j, n) = vParams[j];
512  }
513  m_Psi_ijk[n] = vParams[0];
514  }
515  }
516 }
517 
518 void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
519 {
520  string xname = BinSalt.name();
521  if (xname != "psiCommonAnion") {
522  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
523  "Incorrect name for processing this routine: " + xname);
524  }
525  vector_fp vParams;
526  string kName = BinSalt.attrib("anion");
527  if (kName == "") {
528  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no anion attrib");
529  }
530  string iName = BinSalt.attrib("cation1");
531  if (iName == "") {
532  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation1 attrib");
533  }
534  string jName = BinSalt.attrib("cation2");
535  if (jName == "") {
536  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation2 attrib");
537  }
538 
539  // Find the index of the species in the current phase. It's not an error to
540  // not find the species
541  size_t kSpecies = speciesIndex(kName);
542  if (kSpecies == npos) {
543  return;
544  }
545  if (charge(kSpecies) >= 0) {
546  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "anion charge problem");
547  }
548  size_t iSpecies = speciesIndex(iName);
549  if (iSpecies == npos) {
550  return;
551  }
552  if (charge(iSpecies) <= 0) {
553  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
554  "cation1 charge problem");
555  }
556  size_t jSpecies = speciesIndex(jName);
557  if (jSpecies == npos) {
558  return;
559  }
560  if (charge(jSpecies) <= 0) {
561  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
562  "cation2 charge problem");
563  }
564 
565  size_t n = iSpecies * m_kk + jSpecies;
566  int counter = m_CounterIJ[n];
567  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
568  XML_Node& xmlChild = BinSalt.child(i);
569  if (ba::iequals(xmlChild.name(), "theta")) {
570  double old = m_Theta_ij[counter];
571  m_Theta_ij[counter] = fpValueCheck(xmlChild.value());
572  if (old != 0.0 && old != m_Theta_ij[counter]) {
573  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
574  "conflicting values");
575  }
576  }
577  if (ba::iequals(xmlChild.name(), "psi")) {
578  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
579  size_t nParamsFound = vParams.size();
580  n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies;
581 
582  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
583  if (nParamsFound != 1) {
584  throw CanteraError("HMWSoln::readXMLPsiCommonAnion::Psi for "
585  + kName + "::" + iName + "::" + jName,
586  "wrong number of params found");
587  }
588  m_Psi_ijk_coeff(0,n) = vParams[0];
589  m_Psi_ijk[n] = vParams[0];
590  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
591  if (nParamsFound != 2) {
592  throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for "
593  + kName + "::" + iName + "::" + jName,
594  "wrong number of params found");
595  }
596  m_Psi_ijk_coeff(0,n) = vParams[0];
597  m_Psi_ijk_coeff(1,n) = vParams[1];
598  m_Psi_ijk[n] = vParams[0];
599  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
600  if (nParamsFound == 1) {
601  vParams.resize(5, 0.0);
602  nParamsFound = 5;
603  } else if (nParamsFound != 5) {
604  throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for "
605  + kName + "::" + iName + "::" + jName,
606  "wrong number of params found");
607  }
608  for (size_t j = 0; j < nParamsFound; j++) {
609  m_Psi_ijk_coeff(j, n) = vParams[j];
610  }
611  m_Psi_ijk[n] = vParams[0];
612  }
613 
614  // fill in the duplicate entries
615  n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies;
616  for (size_t j = 0; j < nParamsFound; j++) {
617  m_Psi_ijk_coeff(j, n) = vParams[j];
618  }
619  m_Psi_ijk[n] = vParams[0];
620 
621  n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies;
622  for (size_t j = 0; j < nParamsFound; j++) {
623  m_Psi_ijk_coeff(j, n) = vParams[j];
624  }
625  m_Psi_ijk[n] = vParams[0];
626 
627  n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies;
628  for (size_t j = 0; j < nParamsFound; j++) {
629  m_Psi_ijk_coeff(j, n) = vParams[j];
630  }
631  m_Psi_ijk[n] = vParams[0];
632 
633  n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies;
634  for (size_t j = 0; j < nParamsFound; j++) {
635  m_Psi_ijk_coeff(j, n) = vParams[j];
636  }
637  m_Psi_ijk[n] = vParams[0];
638 
639  n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies;
640  for (size_t j = 0; j < nParamsFound; j++) {
641  m_Psi_ijk_coeff(j, n) = vParams[j];
642  }
643  m_Psi_ijk[n] = vParams[0];
644  }
645  }
646 }
647 
648 void HMWSoln::readXMLLambdaNeutral(XML_Node& BinSalt)
649 {
650  string xname = BinSalt.name();
651  vector_fp vParams;
652  if (xname != "lambdaNeutral") {
653  throw CanteraError("HMWSoln::readXMLLanbdaNeutral",
654  "Incorrect name for processing this routine: " + xname);
655  }
656  string stemp;
657  string iName = BinSalt.attrib("species1");
658  if (iName == "") {
659  throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no species1 attrib");
660  }
661  string jName = BinSalt.attrib("species2");
662  if (jName == "") {
663  throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no species2 attrib");
664  }
665 
666  // Find the index of the species in the current phase. It's not an error to
667  // not find the species
668  size_t iSpecies = speciesIndex(iName);
669  if (iSpecies == npos) {
670  return;
671  }
672  if (charge(iSpecies) != 0) {
673  throw CanteraError("HMWSoln::readXMLLambdaNeutral",
674  "neutral charge problem");
675  }
676  size_t jSpecies = speciesIndex(jName);
677  if (jSpecies == npos) {
678  return;
679  }
680 
681  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
682  XML_Node& xmlChild = BinSalt.child(i);
683  if (ba::iequals(xmlChild.name(), "lambda")) {
684  size_t nCount = iSpecies*m_kk + jSpecies;
685  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
686  size_t nParamsFound = vParams.size();
687  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
688  if (nParamsFound != 1) {
689  throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
690  + "::" + jName,
691  "wrong number of params found");
692  }
693  m_Lambda_nj_coeff(0,nCount) = vParams[0];
694  m_Lambda_nj(iSpecies,jSpecies) = vParams[0];
695  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
696  if (nParamsFound != 2) {
697  throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
698  + "::" + jName,
699  "wrong number of params found");
700  }
701  m_Lambda_nj_coeff(0,nCount) = vParams[0];
702  m_Lambda_nj_coeff(1,nCount) = vParams[1];
703  m_Lambda_nj(iSpecies, jSpecies) = vParams[0];
704  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
705  if (nParamsFound == 1) {
706  vParams.resize(5, 0.0);
707  nParamsFound = 5;
708  } else if (nParamsFound != 5) {
709  throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
710  + "::" + jName,
711  "wrong number of params found");
712  }
713  for (size_t j = 0; j < nParamsFound; j++) {
714  m_Lambda_nj_coeff(j,nCount) = vParams[j];
715  }
716  m_Lambda_nj(iSpecies, jSpecies) = vParams[0];
717  }
718  }
719  }
720 }
721 
722 void HMWSoln::readXMLMunnnNeutral(XML_Node& BinSalt)
723 {
724  string xname = BinSalt.name();
725  vector_fp vParams;
726  if (xname != "MunnnNeutral") {
727  throw CanteraError("HMWSoln::readXMLMunnnNeutral",
728  "Incorrect name for processing this routine: " + xname);
729  }
730  string iName = BinSalt.attrib("species1");
731  if (iName == "") {
732  throw CanteraError("HMWSoln::readXMLMunnnNeutral", "no species1 attrib");
733  }
734 
735  // Find the index of the species in the current phase. It's not an error to
736  // not find the species
737  size_t iSpecies = speciesIndex(iName);
738  if (iSpecies == npos) {
739  return;
740  }
741  if (charge(iSpecies) != 0) {
742  throw CanteraError("HMWSoln::readXMLMunnnNeutral",
743  "neutral charge problem");
744  }
745 
746  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
747  XML_Node& xmlChild = BinSalt.child(i);
748  if (ba::iequals(xmlChild.name(), "Munnn")) {
749  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
750  size_t nParamsFound = vParams.size();
751  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
752  if (nParamsFound != 1) {
753  throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
754  "wrong number of params found");
755  }
756  m_Mu_nnn_coeff(0,iSpecies) = vParams[0];
757  m_Mu_nnn[iSpecies] = vParams[0];
758  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
759  if (nParamsFound != 2) {
760  throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
761  "wrong number of params found");
762  }
763  m_Mu_nnn_coeff(0, iSpecies) = vParams[0];
764  m_Mu_nnn_coeff(1, iSpecies) = vParams[1];
765  m_Mu_nnn[iSpecies] = vParams[0];
766  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
767  if (nParamsFound == 1) {
768  vParams.resize(5, 0.0);
769  nParamsFound = 5;
770  } else if (nParamsFound != 5) {
771  throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
772  "wrong number of params found");
773  }
774  for (size_t j = 0; j < nParamsFound; j++) {
775  m_Mu_nnn_coeff(j, iSpecies) = vParams[j];
776  }
777  m_Mu_nnn[iSpecies] = vParams[0];
778  }
779  }
780  }
781 }
782 
783 void HMWSoln::readXMLZetaCation(const XML_Node& BinSalt)
784 {
785  string xname = BinSalt.name();
786  if (xname != "zetaCation") {
787  throw CanteraError("HMWSoln::readXMLZetaCation",
788  "Incorrect name for processing this routine: " + xname);
789  }
790  vector_fp vParams;
791 
792  string iName = BinSalt.attrib("neutral");
793  if (iName == "") {
794  throw CanteraError("HMWSoln::readXMLZetaCation", "no neutral attrib");
795  }
796 
797  string jName = BinSalt.attrib("cation1");
798  if (jName == "") {
799  throw CanteraError("HMWSoln::readXMLZetaCation", "no cation1 attrib");
800  }
801 
802  string kName = BinSalt.attrib("anion1");
803  if (kName == "") {
804  throw CanteraError("HMWSoln::readXMLZetaCation", "no anion1 attrib");
805  }
806 
807  // Find the index of the species in the current phase. It's not an error to
808  // not find the species
809  size_t iSpecies = speciesIndex(iName);
810  if (iSpecies == npos) {
811  return;
812  }
813  if (charge(iSpecies) != 0.0) {
814  throw CanteraError("HMWSoln::readXMLZetaCation", "neutral charge problem");
815  }
816 
817  size_t jSpecies = speciesIndex(jName);
818  if (jSpecies == npos) {
819  return;
820  }
821  if (charge(jSpecies) <= 0.0) {
822  throw CanteraError("HMWSoln::readXLZetaCation", "cation1 charge problem");
823  }
824 
825  size_t kSpecies = speciesIndex(kName);
826  if (kSpecies == npos) {
827  return;
828  }
829  if (charge(kSpecies) >= 0.0) {
830  throw CanteraError("HMWSoln::readXMLZetaCation", "anion1 charge problem");
831  }
832 
833  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
834  XML_Node& xmlChild = BinSalt.child(i);
835  if (ba::iequals(xmlChild.name(), "zeta")) {
836  getFloatArray(xmlChild, vParams, false, "", xmlChild.name());
837  size_t nParamsFound = vParams.size();
838  size_t n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies;
839 
840  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
841  if (nParamsFound != 1) {
842  throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
843  + iName + "::" + jName + "::" + kName,
844  "wrong number of params found");
845  }
846  m_Psi_ijk_coeff(0,n) = vParams[0];
847  m_Psi_ijk[n] = vParams[0];
848  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
849  if (nParamsFound != 2) {
850  throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
851  + iName + "::" + jName + "::" + kName,
852  "wrong number of params found");
853  }
854  m_Psi_ijk_coeff(0,n) = vParams[0];
855  m_Psi_ijk_coeff(1,n) = vParams[1];
856  m_Psi_ijk[n] = vParams[0];
857  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
858  if (nParamsFound == 1) {
859  vParams.resize(5, 0.0);
860  nParamsFound = 5;
861  } else if (nParamsFound != 5) {
862  throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
863  + iName + "::" + jName + "::" + kName,
864  "wrong number of params found");
865  }
866  for (size_t j = 0; j < nParamsFound; j++) {
867  m_Psi_ijk_coeff(j, n) = vParams[j];
868  }
869  m_Psi_ijk[n] = vParams[0];
870  }
871  // There are no duplicate entries
872  }
873  }
874 }
875 
876 void HMWSoln::readXMLCroppingCoefficients(const XML_Node& acNode)
877 {
878  if (acNode.hasChild("croppingCoefficients")) {
879  XML_Node& cropNode = acNode.child("croppingCoefficients");
880  if (cropNode.hasChild("ln_gamma_k_min")) {
881  XML_Node& gkminNode = cropNode.child("ln_gamma_k_min");
882  getOptionalFloat(gkminNode, "pureSolventValue", CROP_ln_gamma_k_min);
883  }
884  if (cropNode.hasChild("ln_gamma_k_max")) {
885  XML_Node& gkmaxNode = cropNode.child("ln_gamma_k_max");
886  getOptionalFloat(gkmaxNode, "pureSolventValue", CROP_ln_gamma_k_max);
887  }
888 
889  if (cropNode.hasChild("ln_gamma_o_min")) {
890  XML_Node& gominNode = cropNode.child("ln_gamma_o_min");
891  getOptionalFloat(gominNode, "pureSolventValue", CROP_ln_gamma_o_min);
892  }
893 
894  if (cropNode.hasChild("ln_gamma_o_max")) {
895  XML_Node& gomaxNode = cropNode.child("ln_gamma_o_max");
896  getOptionalFloat(gomaxNode, "pureSolventValue", CROP_ln_gamma_o_max);
897  }
898  }
899 }
900 
901 void HMWSoln::initThermo()
902 {
903  MolalityVPSSTP::initThermo();
904  for (int i = 0; i < 17; i++) {
905  elambda[i] = 0.0;
906  elambda1[i] = 0.0;
907  }
908  initLengths();
909 }
910 
911 void HMWSoln::constructPhaseFile(std::string inputFile, std::string id_)
912 {
913  warn_deprecated("HMWSoln::constructPhaseFile",
914  "Use initThermoFile instead. To be removed after Cantera 2.3.");
915 
916  initThermoFile(inputFile, id_);
917 }
918 
919 void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id_)
920 {
921  warn_deprecated("HMWSoln::constructPhaseXML",
922  "Use importPhase instead. To be removed after Cantera 2.3.");
923  importPhase(phaseNode, this);
924 }
925 
926 void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
927 {
928  if (id_.size() > 0) {
929  string idp = phaseNode.id();
930  if (idp != id_) {
931  throw CanteraError("HMWSoln::initThermoXML",
932  "phasenode and Id are incompatible");
933  }
934  }
935 
936  // Find the Thermo XML node
937  if (!phaseNode.hasChild("thermo")) {
938  throw CanteraError("HMWSoln::initThermoXML",
939  "no thermo XML node");
940  }
941  XML_Node& thermoNode = phaseNode.child("thermo");
942 
943  // Possibly change the form of the standard concentrations
944  if (thermoNode.hasChild("standardConc")) {
945  XML_Node& scNode = thermoNode.child("standardConc");
946  m_formGC = 2;
947  string formString = scNode.attrib("model");
948  if (formString != "") {
949  if (ba::iequals(formString, "unity")) {
950  m_formGC = 0;
951  throw CanteraError("HMWSoln::initThermoXML",
952  "standardConc = unity not done");
953  } else if (ba::iequals(formString, "molar_volume")) {
954  m_formGC = 1;
955  throw CanteraError("HMWSoln::initThermoXML",
956  "standardConc = molar_volume not done");
957  } else if (ba::iequals(formString, "solvent_volume")) {
958  m_formGC = 2;
959  } else {
960  throw CanteraError("HMWSoln::initThermoXML",
961  "Unknown standardConc model: " + formString);
962  }
963  }
964  }
965 
966  // Determine the form of the Pitzer model, We will use this information to
967  // size arrays below.
968  if (thermoNode.hasChild("activityCoefficients")) {
969  XML_Node& scNode = thermoNode.child("activityCoefficients");
970  string formString = scNode.attrib("model");
971  if (formString != "") {
972  if (ba::iequals(formString, "pitzer") || ba::iequals(formString, "default")) {
973  m_formPitzer = PITZERFORM_BASE;
974  } else if (ba::iequals(formString, "base")) {
975  m_formPitzer = PITZERFORM_BASE;
976  } else {
977  throw CanteraError("HMWSoln::initThermoXML",
978  "Unknown Pitzer ActivityCoeff model: "
979  + formString);
980  }
981  }
982 
983  // Determine the form of the temperature dependence of the Pitzer
984  // activity coefficient model.
985  formString = scNode.attrib("TempModel");
986  if (formString != "") {
987  if (ba::iequals(formString, "constant") || ba::iequals(formString, "default")) {
988  m_formPitzerTemp = PITZER_TEMP_CONSTANT;
989  } else if (ba::iequals(formString, "linear")) {
990  m_formPitzerTemp = PITZER_TEMP_LINEAR;
991  } else if (ba::iequals(formString, "complex") || ba::iequals(formString, "complex1")) {
992  m_formPitzerTemp = PITZER_TEMP_COMPLEX1;
993  } else {
994  throw CanteraError("HMWSoln::initThermoXML",
995  "Unknown Pitzer ActivityCoeff Temp model: "
996  + formString);
997  }
998  }
999 
1000  // Determine the reference temperature of the Pitzer activity
1001  // coefficient model's temperature dependence formulation: defaults to
1002  // 25C
1003  formString = scNode.attrib("TempReference");
1004  if (formString != "") {
1005  m_TempPitzerRef = fpValueCheck(formString);
1006  } else {
1007  m_TempPitzerRef = 273.15 + 25;
1008  }
1009  }
1010 
1011  // Get the Name of the Solvent:
1012  // <solvent> solventName </solvent>
1013  string solventName = "";
1014  if (thermoNode.hasChild("solvent")) {
1015  XML_Node& scNode = thermoNode.child("solvent");
1016  vector<string> nameSolventa;
1017  getStringArray(scNode, nameSolventa);
1018  if (nameSolventa.size() != 1) {
1019  throw CanteraError("HMWSoln::initThermoXML",
1020  "badly formed solvent XML node");
1021  }
1022  solventName = nameSolventa[0];
1023  }
1024 
1025  // Initialize all of the lengths of arrays in the object
1026  // now that we know what species are in the phase.
1027  initLengths();
1028 
1029  // Reconcile the solvent name and index.
1030  for (size_t k = 0; k < m_kk; k++) {
1031  string sname = speciesName(k);
1032  if (solventName == sname) {
1033  setSolvent(k);
1034  if (k != 0) {
1035  throw CanteraError("HMWSoln::initThermoXML",
1036  "Solvent must be species 0 atm");
1037  }
1038  m_indexSolvent = k;
1039  break;
1040  }
1041  }
1042  if (m_indexSolvent == npos) {
1043  std::cout << "HMWSoln::initThermo: Solvent Name not found"
1044  << std::endl;
1045  throw CanteraError("HMWSoln::initThermoXML",
1046  "Solvent name not found");
1047  }
1048  if (m_indexSolvent != 0) {
1049  throw CanteraError("HMWSoln::initThermoXML",
1050  "Solvent " + solventName +
1051  " should be first species");
1052  }
1053 
1054  // Now go get the specification of the standard states for species in the
1055  // solution. This includes the molar volumes data blocks for incompressible
1056  // species.
1057  XML_Node& speciesList = phaseNode.child("speciesArray");
1058  XML_Node* speciesDB =
1059  get_XML_NameID("speciesData", speciesList["datasrc"],
1060  &phaseNode.root());
1061  const vector<string>&sss = speciesNames();
1062 
1063  for (size_t k = 0; k < m_kk; k++) {
1064  XML_Node* s = speciesDB->findByAttr("name", sss[k]);
1065  if (!s) {
1066  throw CanteraError("HMWSoln::initThermoXML",
1067  "Species Data Base " + sss[k] + " not found");
1068  }
1069  XML_Node* ss = s->findByName("standardState");
1070  if (!ss) {
1071  throw CanteraError("HMWSoln::initThermoXML",
1072  "Species " + sss[k] +
1073  " standardState XML block not found");
1074  }
1075  string modelString = ba::to_lower_copy(ss->attrib("model"));
1076  if (modelString == "") {
1077  throw CanteraError("HMWSoln::initThermoXML",
1078  "Species " + sss[k] +
1079  " standardState XML block model attribute not found");
1080  }
1081  if (k == 0) {
1082  if (modelString == "wateriapws" || modelString == "real_water" ||
1083  modelString == "waterpdss") {
1084 
1085  // Store a local pointer to the water standard state model.
1086  // We've hardcoded it to a PDSS_Water model, so this is ok.
1087  m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0));
1088  if (!m_waterSS) {
1089  throw CanteraError("HMWSoln::initThermoXML",
1090  "Dynamic cast to PDSS_Water failed");
1091  }
1092 
1093  // Fill in the molar volume of water (m3/kmol) at standard
1094  // conditions to fill in the m_speciesSize entry with something
1095  // reasonable.
1096  m_waterSS->setState_TP(300., OneAtm);
1097  double dens = m_waterSS->density();
1098  double mw = m_waterSS->molecularWeight();
1099  m_speciesSize[0] = mw / dens;
1100  } else {
1101  m_waterSS = providePDSS(0);
1102  m_waterSS->setState_TP(300., OneAtm);
1103  double dens = m_waterSS->density();
1104  double mw = m_waterSS->molecularWeight();
1105  m_speciesSize[0] = mw / dens;
1106  }
1107  } else {
1108  if (modelString != "constant_incompressible" && modelString != "hkft") {
1109  throw CanteraError("HMWSoln::initThermoXML",
1110  "Solute SS Model \"" + modelString +
1111  "\" is not known");
1112  }
1113  if (modelString == "constant_incompressible") {
1114  m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSI");
1115  }
1116  // HKM Note, have to fill up m_speciesSize[] for HKFT species
1117  }
1118  }
1119 
1120  // Initialize the water property calculator. It will share the internal eos
1121  // water calculator.
1122  m_waterProps.reset(new WaterProps(&dynamic_cast<PDSS_Water&>(*m_waterSS)));
1123 
1124  // Fill in parameters for the calculation of the stoichiometric Ionic
1125  // Strength. The default is that stoich charge is the same as the regular
1126  // charge.
1127  for (size_t k = 0; k < m_kk; k++) {
1128  m_speciesCharge_Stoich[k] = charge(k);
1129  }
1130 
1131  // Go get all of the coefficients and factors in the activityCoefficients
1132  // XML block
1133  XML_Node* acNodePtr = 0;
1134  if (thermoNode.hasChild("activityCoefficients")) {
1135  XML_Node& acNode = thermoNode.child("activityCoefficients");
1136  acNodePtr = &acNode;
1137 
1138  // Look for parameters for A_Debye
1139  if (acNode.hasChild("A_Debye")) {
1140  XML_Node& ADebye = acNode.child("A_Debye");
1141  m_form_A_Debye = A_DEBYE_CONST;
1142  string stemp = "model";
1143  if (ADebye.hasAttrib("model")) {
1144  if (ba::iequals(ADebye.attrib("model"), "water")) {
1145  m_form_A_Debye = A_DEBYE_WATER;
1146  }
1147  }
1148  if (m_form_A_Debye == A_DEBYE_CONST) {
1149  m_A_Debye = getFloat(acNode, "A_Debye");
1150  }
1151  }
1152 
1153  // Look for Parameters for the Maximum Ionic Strength
1154  if (acNode.hasChild("maxIonicStrength")) {
1155  m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
1156  }
1157 
1158  // Look for parameters for the Ionic radius
1159  if (acNode.hasChild("ionicRadius")) {
1160  XML_Node& irNode = acNode.child("ionicRadius");
1161  double Afactor = 1.0;
1162  if (irNode.hasAttrib("units")) {
1163  string Aunits = irNode.attrib("units");
1164  Afactor = toSI(Aunits);
1165  }
1166 
1167  if (irNode.hasAttrib("default")) {
1168  string ads = irNode.attrib("default");
1169  double ad = fpValue(ads);
1170  for (size_t k = 0; k < m_kk; k++) {
1171  m_Aionic[k] = ad * Afactor;
1172  }
1173  }
1174  }
1175 
1176  // First look at the species database. Look for the subelement
1177  // "stoichIsMods" in each of the species SS databases.
1178  std::vector<const XML_Node*> xspecies = speciesData();
1179  for (size_t k = 0; k < m_kk; k++) {
1180  size_t jmap = npos;
1181  string kname = speciesName(k);
1182  for (size_t j = 0; j < xspecies.size(); j++) {
1183  const XML_Node& sp = *xspecies[j];
1184  string jname = sp["name"];
1185  if (jname == kname) {
1186  jmap = j;
1187  break;
1188  }
1189  }
1190  if (jmap != npos) {
1191  const XML_Node& sp = *xspecies[jmap];
1192  getOptionalFloat(sp, "stoichIsMods", m_speciesCharge_Stoich[k]);
1193  }
1194  }
1195 
1196  // Now look at the activity coefficient database
1197  if (acNodePtr && acNodePtr->hasChild("stoichIsMods")) {
1198  XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
1199  map<string, string> msIs;
1200  getMap(sIsNode, msIs);
1201  for (const auto& b : msIs) {
1202  size_t kk = speciesIndex(b.first);
1203  if (kk != npos) {
1204  double val = fpValue(b.second);
1205  m_speciesCharge_Stoich[kk] = val;
1206  }
1207  }
1208  }
1209 
1210  // Loop through the children getting multiple instances of parameters
1211  if (acNodePtr) {
1212  for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
1213  XML_Node& xmlACChild = acNodePtr->child(i);
1214  string nodeName = xmlACChild.name();
1215 
1216  // Process a binary salt field, or any of the other XML fields
1217  // that make up the Pitzer Database. Entries will be ignored
1218  // if any of the species in the entry isn't in the solution.
1219  if (ba::iequals(nodeName, "binarysaltparameters")) {
1220  readXMLBinarySalt(xmlACChild);
1221  } else if (ba::iequals(nodeName, "thetaanion")) {
1222  readXMLThetaAnion(xmlACChild);
1223  } else if (ba::iequals(nodeName, "thetacation")) {
1224  readXMLThetaCation(xmlACChild);
1225  } else if (ba::iequals(nodeName, "psicommonanion")) {
1226  readXMLPsiCommonAnion(xmlACChild);
1227  } else if (ba::iequals(nodeName, "psicommoncation")) {
1228  readXMLPsiCommonCation(xmlACChild);
1229  } else if (ba::iequals(nodeName, "lambdaneutral")) {
1230  readXMLLambdaNeutral(xmlACChild);
1231  } else if (ba::iequals(nodeName, "zetacation")) {
1232  readXMLZetaCation(xmlACChild);
1233  }
1234  }
1235  }
1236 
1237  // Go look up the optional Cropping parameters
1238  readXMLCroppingCoefficients(acNode);
1239  }
1240 
1241  // Fill in the vector specifying the electrolyte species type
1242  //
1243  // First fill in default values. Everything is either a charge species, a
1244  // nonpolar neutral, or the solvent.
1245  for (size_t k = 0; k < m_kk; k++) {
1246  if (fabs(charge(k)) > 0.0001) {
1247  m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
1248  if (fabs(m_speciesCharge_Stoich[k] - charge(k)) > 0.0001) {
1249  m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
1250  }
1251  } else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
1252  m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
1253  } else {
1254  m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
1255  }
1256  }
1257  m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
1258 
1259  // First look at the species database. Look for the subelement
1260  // "stoichIsMods" in each of the species SS databases.
1261  std::vector<const XML_Node*> xspecies = speciesData();
1262  for (size_t k = 0; k < m_kk; k++) {
1263  const XML_Node* spPtr = xspecies[k];
1264  if (spPtr && spPtr->hasChild("electrolyteSpeciesType")) {
1265  string est = getChildValue(*spPtr, "electrolyteSpeciesType");
1266  if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
1267  throw CanteraError("HMWSoln::initThermoXML",
1268  "Bad electrolyte type: " + est);
1269  }
1270  }
1271  }
1272 
1273  // Then look at the phase thermo specification
1274  if (acNodePtr && acNodePtr->hasChild("electrolyteSpeciesType")) {
1275  XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
1276  map<string, string> msEST;
1277  getMap(ESTNode, msEST);
1278  for (const auto& b : msEST) {
1279  size_t kk = speciesIndex(b.first);
1280  if (kk != npos) {
1281  string est = b.second;
1282  if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
1283  throw CanteraError("HMWSoln::initThermoXML",
1284  "Bad electrolyte type: " + est);
1285  }
1286  }
1287  }
1288  }
1289 
1290  IMS_typeCutoff_ = 2;
1291  if (IMS_typeCutoff_ == 2) {
1292  calcIMSCutoffParams_();
1293  }
1294  calcMCCutoffParams_();
1295  setMoleFSolventMin(1.0E-5);
1296 
1297  MolalityVPSSTP::initThermoXML(phaseNode, id_);
1298 
1299  // Lastly calculate the charge balance and then add stuff until the charges
1300  // compensate
1301  vector_fp mf(m_kk, 0.0);
1302  getMoleFractions(mf.data());
1303  bool notDone = true;
1304 
1305  while (notDone) {
1306  double sum = 0.0;
1307  size_t kMaxC = npos;
1308  double MaxC = 0.0;
1309  for (size_t k = 0; k < m_kk; k++) {
1310  sum += mf[k] * charge(k);
1311  if (fabs(mf[k] * charge(k)) > MaxC) {
1312  kMaxC = k;
1313  }
1314  }
1315  size_t kHp = speciesIndex("H+");
1316  size_t kOHm = speciesIndex("OH-");
1317 
1318  if (fabs(sum) > 1.0E-30) {
1319  if (kHp != npos) {
1320  if (mf[kHp] > sum * 1.1) {
1321  mf[kHp] -= sum;
1322  mf[0] += sum;
1323  notDone = false;
1324  } else {
1325  if (sum > 0.0) {
1326  mf[kHp] *= 0.5;
1327  mf[0] += mf[kHp];
1328  sum -= mf[kHp];
1329  }
1330  }
1331  }
1332  if (notDone) {
1333  if (kOHm != npos) {
1334  if (mf[kOHm] > -sum * 1.1) {
1335  mf[kOHm] += sum;
1336  mf[0] -= sum;
1337  notDone = false;
1338  } else {
1339  if (sum < 0.0) {
1340  mf[kOHm] *= 0.5;
1341  mf[0] += mf[kOHm];
1342  sum += mf[kOHm];
1343  }
1344  }
1345  }
1346  if (notDone && kMaxC != npos) {
1347  if (mf[kMaxC] > (1.1 * sum / charge(kMaxC))) {
1348  mf[kMaxC] -= sum / charge(kMaxC);
1349  mf[0] += sum / charge(kMaxC);
1350  } else {
1351  mf[kMaxC] *= 0.5;
1352  mf[0] += mf[kMaxC];
1353  notDone = true;
1354  }
1355  }
1356  }
1357  setMoleFractions(mf.data());
1358  } else {
1359  notDone = false;
1360  }
1361  }
1362 }
1363 
1364 void HMWSoln::calcIMSCutoffParams_()
1365 {
1366  IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
1367  IMS_efCut_ = 0.0;
1368  bool converged = false;
1369  double oldV = 0.0;
1370  for (int its = 0; its < 100 && !converged; its++) {
1371  oldV = IMS_efCut_;
1372  IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
1373  IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
1374  IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
1375  /
1376  (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
1377  double tmp = IMS_afCut_ + IMS_X_o_cutoff_*(IMS_bfCut_ + IMS_dfCut_ *IMS_X_o_cutoff_);
1378  double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1379  IMS_efCut_ = - eterm * tmp;
1380  if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
1381  converged = true;
1382  }
1383  }
1384  if (!converged) {
1385  throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1386  " failed to converge on the f polynomial");
1387  }
1388  converged = false;
1389  double f_0 = IMS_afCut_ + IMS_efCut_;
1390  double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
1391  IMS_egCut_ = 0.0;
1392  for (int its = 0; its < 100 && !converged; its++) {
1393  oldV = IMS_egCut_;
1394  double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
1395  IMS_agCut_ = exp(lng_0) - IMS_egCut_;
1396  IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
1397  IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
1398  /
1399  (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
1400  double tmp = IMS_agCut_ + IMS_X_o_cutoff_*(IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
1401  double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1402  IMS_egCut_ = - eterm * tmp;
1403  if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
1404  converged = true;
1405  }
1406  }
1407  if (!converged) {
1408  throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1409  " failed to converge on the g polynomial");
1410  }
1411 }
1412 
1413 void HMWSoln::calcMCCutoffParams_()
1414 {
1415  MC_X_o_min_ = 0.35;
1416  MC_X_o_cutoff_ = 0.6;
1417  MC_slopepCut_ = 0.02;
1418  MC_cpCut_ = 0.25;
1419 
1420  // Initial starting values
1421  MC_apCut_ = MC_X_o_min_;
1422  MC_epCut_ = 0.0;
1423  bool converged = false;
1424  double oldV = 0.0;
1425  double damp = 0.5;
1426  for (int its = 0; its < 500 && !converged; its++) {
1427  oldV = MC_epCut_;
1428  MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
1429  double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;
1430  MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
1431  double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ * MC_X_o_cutoff_/MC_cpCut_)
1432  /
1433  (MC_X_o_cutoff_ * MC_X_o_cutoff_/MC_cpCut_ - 2.0 * MC_X_o_cutoff_));
1434  MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
1435  double tmp = MC_apCut_ + MC_X_o_cutoff_*(MC_bpCut_ + MC_dpCut_ * MC_X_o_cutoff_);
1436  double eterm = std::exp(- MC_X_o_cutoff_ / MC_cpCut_);
1437  MC_epCut_ = - eterm * tmp;
1438  double diff = MC_epCut_ - oldV;
1439  if (fabs(diff) < 1.0E-14) {
1440  converged = true;
1441  }
1442  }
1443  if (!converged) {
1444  throw CanteraError("HMWSoln::calcMCCutoffParams_()",
1445  " failed to converge on the p polynomial");
1446  }
1447 }
1448 
1449 }
Cantera::getMap
void getMap(const XML_Node &node, std::map< std::string, std::string > &m)
This routine is used to interpret the value portions of XML elements that contain colon separated pai...
Definition: ctml.cpp:373
Cantera::cEST_solvent
const int cEST_solvent
Electrolyte species type.
Definition: electrolytes.h:18
Cantera::fpValue
doublereal fpValue(const std::string &val)
Translate a string into one doublereal value.
Definition: stringUtils.cpp:182
Cantera::getFloatArray
size_t getFloatArray(const XML_Node &node, vector_fp &v, const bool convert, const std::string &unitsString, const std::string &nodeName)
This function reads the current node or a child node of the current node with the default name...
Definition: ctml.cpp:299
ctml.h
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
Cantera::XML_Node::name
std::string name() const
Returns the name of the XML node.
Definition: xml.h:370
Cantera::getChildValue
std::string getChildValue(const XML_Node &parent, const std::string &nameString)
This function reads a child node with the name, nameString, and returns its XML value as the return s...
Definition: ctml.cpp:145
Cantera::interp_est
static int interp_est(const std::string &estString)
Utility function to assign an integer value from a string for the ElectrolyteSpeciesType field...
Definition: DebyeHuckel.cpp:398
Cantera::XML_Node::findByName
const XML_Node * findByName(const std::string &nm, int depth=100000) const
This routine carries out a recursive search for an XML node based on the name of the node...
Definition: xml.cpp:695
Cantera::OneAtm
const doublereal OneAtm
One atmosphere [Pa].
Definition: ct_defs.h:69
Cantera::toSI
doublereal toSI(const std::string &unit)
Return the conversion factor to convert unit std::string &#39;unit&#39; to SI units.
Definition: global.cpp:160
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
electrolytes.h
Header file for a common definitions used in electrolytes thermodynamics.
ThermoFactory.h
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
Cantera::XML_Node
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
PDSS_Water.h
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
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:54
std
STL namespace.
Cantera::XML_Node::hasAttrib
bool hasAttrib(const std::string &a) const
Tests whether the current node has an attribute with a particular name.
Definition: xml.cpp:541
HMWSoln.h
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
Cantera::WaterProps
The WaterProps class is used to house several approximation routines for properties of water...
Definition: WaterProps.h:93
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:37
Cantera::PDSS_Water
Class for the liquid water pressure dependent standard state.
Definition: PDSS_Water.h:49
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
Cantera::XML_Node::value
std::string value() const
Return the value of an XML node as a string.
Definition: xml.cpp:449
Cantera::importPhase
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
Definition: ThermoFactory.cpp:240
Cantera::XML_Node::hasChild
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:536
Cantera::XML_Node::child
XML_Node & child(const size_t n) const
Return a changeable reference to the n&#39;th child of the current node.
Definition: xml.cpp:546
Cantera::XML_Node::root
XML_Node & root() const
Return the root of the current XML_Node tree.
Definition: xml.cpp:1025
Cantera::XML_Node::attrib
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:500
Cantera::vector_fp
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
Cantera::getStringArray
void getStringArray(const XML_Node &node, std::vector< std::string > &v)
This function interprets the value portion of an XML element as a string.
Definition: ctml.cpp:470
Cantera::XML_Node::id
std::string id() const
Return the id attribute, if present.
Definition: xml.cpp:428
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
Cantera::getFloat
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:178
Cantera::XML_Node::findByAttr
XML_Node * findByAttr(const std::string &attr, const std::string &val, int depth=100000) const
This routine carries out a recursive search for an XML node based on an attribute of each XML node...
Definition: xml.cpp:661
Cantera
Namespace for the Cantera kernel.
Definition: application.cpp:29
Cantera::fpValueCheck
doublereal fpValueCheck(const std::string &val)
Translate a string into one doublereal value, with error checking.
Definition: stringUtils.cpp:191
Cantera::XML_Node::nChildren
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:556
Cantera::get_XML_NameID
XML_Node * get_XML_NameID(const std::string &nameTarget, const std::string &file_ID, XML_Node *root)
This routine will locate an XML node in either the input XML tree or in another input file specified ...
Definition: global.cpp:252
Cantera::getOptionalFloat
bool getOptionalFloat(const XML_Node &parent, const std::string &name, doublereal &fltRtn, const std::string &type)
Get an optional floating-point value from a child element.
Definition: ctml.cpp:226
Generated by   doxygen 1.8.14