Cantera  2.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  * Copyright (2006) Sandia Corporation. Under the terms of
12  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
13  * U.S. Government retains certain rights in this software.
14  */
15 #include "cantera/thermo/HMWSoln.h"
16 #include "cantera/thermo/ThermoFactory.h"
17 #include "cantera/thermo/WaterProps.h"
18 #include "cantera/thermo/PDSS_Water.h"
19 #include "cantera/base/stringUtils.h"
20 
21 #include <cstring>
22 #include <cstdlib>
23 #include <cstdio>
24 #include <fstream>
25 
26 using namespace std;
27 using namespace ctml;
28 
29 namespace Cantera
30 {
31 
32 
33 //! utility function to assign an integer value from a string
34 //! for the ElectrolyteSpeciesType field.
35 /*!
36  * @param estString string name of the electrolyte species type
37  */
38 int HMWSoln::interp_est(std::string estString)
39 {
40  const char* cc = estString.c_str();
41  string lcs = lowercase(estString);
42  const char* ccl = lcs.c_str();
43  if (!strcmp(ccl, "solvent")) {
44  return cEST_solvent;
45  } else if (!strcmp(ccl, "chargedspecies")) {
46  return cEST_chargedSpecies;
47  } else if (!strcmp(ccl, "weakacidassociated")) {
48  return cEST_weakAcidAssociated;
49  } else if (!strcmp(ccl, "strongacidassociated")) {
50  return cEST_strongAcidAssociated;
51  } else if (!strcmp(ccl, "polarneutral")) {
52  return cEST_polarNeutral;
53  } else if (!strcmp(ccl, "nonpolarneutral")) {
54  return cEST_nonpolarNeutral;
55  }
56  int retn, rval;
57  if ((retn = sscanf(cc, "%d", &rval)) != 1) {
58  return -1;
59  }
60  return rval;
61 }
62 
63 /*
64  * Process an XML node called "SimpleSaltParameters.
65  * This node contains all of the parameters necessary to describe
66  * the Pitzer model for that particular binary salt.
67  * This function reads the XML file and writes the coefficients
68  * it finds to an internal data structures.
69  */
70 void HMWSoln::readXMLBinarySalt(XML_Node& BinSalt)
71 {
72  string xname = BinSalt.name();
73  if (xname != "binarySaltParameters") {
74  throw CanteraError("HMWSoln::readXMLBinarySalt",
75  "Incorrect name for processing this routine: " + xname);
76  }
77  double* charge = DATA_PTR(m_speciesCharge);
78  string stemp;
79  size_t nParamsFound, i;
80  vector_fp vParams;
81  string iName = BinSalt.attrib("cation");
82  if (iName == "") {
83  throw CanteraError("HMWSoln::readXMLBinarySalt", "no cation attrib");
84  }
85  string jName = BinSalt.attrib("anion");
86  if (jName == "") {
87  throw CanteraError("HMWSoln::readXMLBinarySalt", "no anion attrib");
88  }
89  /*
90  * Find the index of the species in the current phase. It's not
91  * an error to not find the species
92  */
93  size_t iSpecies = speciesIndex(iName);
94  if (iSpecies == npos) {
95  return;
96  }
97  string ispName = speciesName(iSpecies);
98  if (charge[iSpecies] <= 0) {
99  throw CanteraError("HMWSoln::readXMLBinarySalt", "cation charge problem");
100  }
101  size_t jSpecies = speciesIndex(jName);
102  if (jSpecies == npos) {
103  return;
104  }
105  string jspName = speciesName(jSpecies);
106  if (charge[jSpecies] >= 0) {
107  throw CanteraError("HMWSoln::readXMLBinarySalt", "anion charge problem");
108  }
109 
110  size_t n = iSpecies * m_kk + jSpecies;
111  int counter = m_CounterIJ[n];
112  for (size_t iChild = 0; iChild < BinSalt.nChildren(); iChild++) {
113  XML_Node& xmlChild = BinSalt.child(iChild);
114  stemp = xmlChild.name();
115  string nodeName = lowercase(stemp);
116  /*
117  * Process the binary salt child elements
118  */
119  if (nodeName == "beta0") {
120  /*
121  * Get the string containing all of the values
122  */
123  getFloatArray(xmlChild, vParams, false, "", "beta0");
124  nParamsFound = vParams.size();
125  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
126  if (nParamsFound != 1) {
127  throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
128  + "::" + jspName,
129  "wrong number of params found");
130  }
131  m_Beta0MX_ij[counter] = vParams[0];
132  m_Beta0MX_ij_coeff(0,counter) = m_Beta0MX_ij[counter];
133  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
134  if (nParamsFound != 2) {
135  throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
136  + "::" + jspName,
137  "wrong number of params found");
138  }
139  m_Beta0MX_ij_coeff(0,counter) = vParams[0];
140  m_Beta0MX_ij_coeff(1,counter) = vParams[1];
141  m_Beta0MX_ij[counter] = vParams[0];
142  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
143  if (nParamsFound != 5) {
144  throw CanteraError("HMWSoln::readXMLBinarySalt::beta0 for " + ispName
145  + "::" + jspName,
146  "wrong number of params found");
147  }
148  for (i = 0; i < nParamsFound; i++) {
149  m_Beta0MX_ij_coeff(i, counter) = vParams[i];
150  }
151  m_Beta0MX_ij[counter] = vParams[0];
152  }
153  }
154  if (nodeName == "beta1") {
155 
156  /*
157  * Get the string containing all of the values
158  */
159  getFloatArray(xmlChild, vParams, false, "", "beta1");
160  nParamsFound = vParams.size();
161  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
162  if (nParamsFound != 1) {
163  throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
164  + "::" + jspName,
165  "wrong number of params found");
166  }
167  m_Beta1MX_ij[counter] = vParams[0];
168  m_Beta1MX_ij_coeff(0,counter) = m_Beta1MX_ij[counter];
169  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
170  if (nParamsFound != 2) {
171  throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
172  + "::" + jspName,
173  "wrong number of params found");
174  }
175  m_Beta1MX_ij_coeff(0,counter) = vParams[0];
176  m_Beta1MX_ij_coeff(1,counter) = vParams[1];
177  m_Beta1MX_ij[counter] = vParams[0];
178  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
179  if (nParamsFound != 5) {
180  throw CanteraError("HMWSoln::readXMLBinarySalt::beta1 for " + ispName
181  + "::" + jspName,
182  "wrong number of params found");
183  }
184  for (i = 0; i < nParamsFound; i++) {
185  m_Beta1MX_ij_coeff(i, counter) = vParams[i];
186  }
187  m_Beta1MX_ij[counter] = vParams[0];
188  }
189  }
190  if (nodeName == "beta2") {
191  getFloatArray(xmlChild, vParams, false, "", "beta2");
192  nParamsFound = vParams.size();
193  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
194  if (nParamsFound != 1) {
195  throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
196  + "::" + jspName,
197  "wrong number of params found");
198  }
199  m_Beta2MX_ij[counter] = vParams[0];
200  m_Beta2MX_ij_coeff(0,counter) = m_Beta2MX_ij[counter];
201  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
202  if (nParamsFound != 2) {
203  throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
204  + "::" + jspName,
205  "wrong number of params found");
206  }
207  m_Beta2MX_ij_coeff(0,counter) = vParams[0];
208  m_Beta2MX_ij_coeff(1,counter) = vParams[1];
209  m_Beta2MX_ij[counter] = vParams[0];
210  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
211  if (nParamsFound != 5) {
212  throw CanteraError("HMWSoln::readXMLBinarySalt::beta2 for " + ispName
213  + "::" + jspName,
214  "wrong number of params found");
215  }
216  for (i = 0; i < nParamsFound; i++) {
217  m_Beta2MX_ij_coeff(i, counter) = vParams[i];
218  }
219  m_Beta2MX_ij[counter] = vParams[0];
220  }
221 
222  }
223  if (nodeName == "cphi") {
224  /*
225  * Get the string containing all of the values
226  */
227  getFloatArray(xmlChild, vParams, false, "", "Cphi");
228  nParamsFound = vParams.size();
229  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
230  if (nParamsFound != 1) {
231  throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
232  + "::" + jspName,
233  "wrong number of params found");
234  }
235  m_CphiMX_ij[counter] = vParams[0];
236  m_CphiMX_ij_coeff(0,counter) = m_CphiMX_ij[counter];
237  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
238  if (nParamsFound != 2) {
239  throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
240  + "::" + jspName,
241  "wrong number of params found");
242  }
243  m_CphiMX_ij_coeff(0,counter) = vParams[0];
244  m_CphiMX_ij_coeff(1,counter) = vParams[1];
245  m_CphiMX_ij[counter] = vParams[0];
246  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
247  if (nParamsFound != 5) {
248  throw CanteraError("HMWSoln::readXMLBinarySalt::Cphi for " + ispName
249  + "::" + jspName,
250  "wrong number of params found");
251  }
252  for (i = 0; i < nParamsFound; i++) {
253  m_CphiMX_ij_coeff(i, counter) = vParams[i];
254  }
255  m_CphiMX_ij[counter] = vParams[0];
256  }
257  }
258 
259  if (nodeName == "alpha1") {
260  stemp = xmlChild.value();
261  m_Alpha1MX_ij[counter] = atofCheck(stemp.c_str());
262  }
263 
264  if (nodeName == "alpha2") {
265  stemp = xmlChild.value();
266  m_Alpha2MX_ij[counter] = atofCheck(stemp.c_str());
267  }
268  }
269 }
270 
271 /**
272  * Process an XML node called "thetaAnion".
273  * This node contains all of the parameters necessary to describe
274  * the binary interactions between two anions.
275  */
276 void HMWSoln::readXMLThetaAnion(XML_Node& BinSalt)
277 {
278  string xname = BinSalt.name();
279  vector_fp vParams;
280  size_t nParamsFound = 0;
281  if (xname != "thetaAnion") {
282  throw CanteraError("HMWSoln::readXMLThetaAnion",
283  "Incorrect name for processing this routine: " + xname);
284  }
285  double* charge = DATA_PTR(m_speciesCharge);
286  string stemp;
287  string ispName = BinSalt.attrib("anion1");
288  if (ispName == "") {
289  throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion1 attrib");
290  }
291  string jspName = BinSalt.attrib("anion2");
292  if (jspName == "") {
293  throw CanteraError("HMWSoln::readXMLThetaAnion", "no anion2 attrib");
294  }
295  /*
296  * Find the index of the species in the current phase. It's not
297  * an error to not find the species
298  */
299  size_t iSpecies = speciesIndex(ispName);
300  if (iSpecies == npos) {
301  return;
302  }
303  if (charge[iSpecies] >= 0) {
304  throw CanteraError("HMWSoln::readXMLThetaAnion", "anion1 charge problem");
305  }
306  size_t jSpecies = speciesIndex(jspName);
307  if (jSpecies == npos) {
308  return;
309  }
310  if (charge[jSpecies] >= 0) {
311  throw CanteraError("HMWSoln::readXMLThetaAnion", "anion2 charge problem");
312  }
313 
314  size_t n = iSpecies * m_kk + jSpecies;
315  int counter = m_CounterIJ[n];
316  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
317  XML_Node& xmlChild = BinSalt.child(i);
318  stemp = xmlChild.name();
319  string nodeName = lowercase(stemp);
320  if (nodeName == "theta") {
321  getFloatArray(xmlChild, vParams, false, "", stemp);
322  nParamsFound = vParams.size();
323  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
324  if (nParamsFound != 1) {
325  throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
326  + "::" + jspName,
327  "wrong number of params found");
328  }
329  m_Theta_ij_coeff(0,counter) = vParams[0];
330  m_Theta_ij[counter] = vParams[0];
331  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
332  if (nParamsFound != 2) {
333  throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
334  + "::" + jspName,
335  "wrong number of params found");
336  }
337  m_Theta_ij_coeff(0,counter) = vParams[0];
338  m_Theta_ij_coeff(1,counter) = vParams[1];
339  m_Theta_ij[counter] = vParams[0];
340  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
341  if (nParamsFound == 1) {
342  vParams.resize(5, 0.0);
343  nParamsFound = 5;
344  } else if (nParamsFound != 5) {
345  throw CanteraError("HMWSoln::readXMLThetaAnion::Theta for " + ispName
346  + "::" + jspName,
347  "wrong number of params found");
348  }
349  for (size_t j = 0; j < nParamsFound; j++) {
350  m_Theta_ij_coeff(j, counter) = vParams[j];
351  }
352  m_Theta_ij[counter] = vParams[0];
353  }
354  }
355  }
356 }
357 
358 /**
359  * Process an XML node called "thetaCation".
360  * This node contains all of the parameters necessary to describe
361  * the binary interactions between two cation.
362  */
363 void HMWSoln::readXMLThetaCation(XML_Node& BinSalt)
364 {
365  string xname = BinSalt.name();
366  vector_fp vParams;
367  size_t nParamsFound = 0;
368  if (xname != "thetaCation") {
369  throw CanteraError("HMWSoln::readXMLThetaCation",
370  "Incorrect name for processing this routine: " + xname);
371  }
372  double* charge = DATA_PTR(m_speciesCharge);
373  string stemp;
374  string ispName = BinSalt.attrib("cation1");
375  if (ispName == "") {
376  throw CanteraError("HMWSoln::readXMLThetaCation", "no cation1 attrib");
377  }
378  string jspName = BinSalt.attrib("cation2");
379  if (jspName == "") {
380  throw CanteraError("HMWSoln::readXMLThetaCation", "no cation2 attrib");
381  }
382  /*
383  * Find the index of the species in the current phase. It's not
384  * an error to not find the species
385  */
386  size_t iSpecies = speciesIndex(ispName);
387  if (iSpecies == npos) {
388  return;
389  }
390  if (charge[iSpecies] <= 0) {
391  throw CanteraError("HMWSoln::readXMLThetaCation", "cation1 charge problem");
392  }
393  size_t jSpecies = speciesIndex(jspName);
394  if (jSpecies == npos) {
395  return;
396  }
397  if (charge[jSpecies] <= 0) {
398  throw CanteraError("HMWSoln::readXMLThetaCation", "cation2 charge problem");
399  }
400 
401  size_t n = iSpecies * m_kk + jSpecies;
402  int counter = m_CounterIJ[n];
403  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
404  XML_Node& xmlChild = BinSalt.child(i);
405  stemp = xmlChild.name();
406  string nodeName = lowercase(stemp);
407  if (nodeName == "theta") {
408  getFloatArray(xmlChild, vParams, false, "", stemp);
409  nParamsFound = vParams.size();
410  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
411  if (nParamsFound != 1) {
412  throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
413  + "::" + jspName,
414  "wrong number of params found");
415  }
416  m_Theta_ij_coeff(0,counter) = vParams[0];
417  m_Theta_ij[counter] = vParams[0];
418  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
419  if (nParamsFound != 2) {
420  throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
421  + "::" + jspName,
422  "wrong number of params found");
423  }
424  m_Theta_ij_coeff(0,counter) = vParams[0];
425  m_Theta_ij_coeff(1,counter) = vParams[1];
426  m_Theta_ij[counter] = vParams[0];
427  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
428  if (nParamsFound == 1) {
429  vParams.resize(5, 0.0);
430  nParamsFound = 5;
431  } else if (nParamsFound != 5) {
432  throw CanteraError("HMWSoln::readXMLThetaCation::Theta for " + ispName
433  + "::" + jspName,
434  "wrong number of params found");
435  }
436  for (size_t j = 0; j < nParamsFound; j++) {
437  m_Theta_ij_coeff(j, counter) = vParams[j];
438  }
439  m_Theta_ij[counter] = vParams[0];
440  }
441  }
442  }
443 }
444 
445 /*
446  * Process an XML node called "readXMLPsiCommonCation".
447  * This node contains all of the parameters necessary to describe
448  * the binary interactions between two anions and one common cation.
449  */
450 void HMWSoln::readXMLPsiCommonCation(XML_Node& BinSalt)
451 {
452  string xname = BinSalt.name();
453  if (xname != "psiCommonCation") {
454  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
455  "Incorrect name for processing this routine: " + xname);
456  }
457  double* charge = DATA_PTR(m_speciesCharge);
458  string stemp;
459  vector_fp vParams;
460  size_t nParamsFound = 0;
461  string kName = BinSalt.attrib("cation");
462  if (kName == "") {
463  throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no cation attrib");
464  }
465  string iName = BinSalt.attrib("anion1");
466  if (iName == "") {
467  throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion1 attrib");
468  }
469  string jName = BinSalt.attrib("anion2");
470  if (jName == "") {
471  throw CanteraError("HMWSoln::readXMLPsiCommonCation", "no anion2 attrib");
472  }
473  /*
474  * Find the index of the species in the current phase. It's not
475  * an error to not find the species
476  */
477  size_t kSpecies = speciesIndex(kName);
478  if (kSpecies == npos) {
479  return;
480  }
481  if (charge[kSpecies] <= 0) {
482  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
483  "cation charge problem");
484  }
485  size_t iSpecies = speciesIndex(iName);
486  if (iSpecies == npos) {
487  return;
488  }
489  if (charge[iSpecies] >= 0) {
490  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
491  "anion1 charge problem");
492  }
493  size_t jSpecies = speciesIndex(jName);
494  if (jSpecies == npos) {
495  return;
496  }
497  if (charge[jSpecies] >= 0) {
498  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
499  "anion2 charge problem");
500  }
501 
502  size_t n = iSpecies * m_kk + jSpecies;
503  int counter = m_CounterIJ[n];
504  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
505  XML_Node& xmlChild = BinSalt.child(i);
506  stemp = xmlChild.name();
507  string nodeName = lowercase(stemp);
508  if (nodeName == "theta") {
509  stemp = xmlChild.value();
510  double old = m_Theta_ij[counter];
511  m_Theta_ij[counter] = atofCheck(stemp.c_str());
512  if (old != 0.0) {
513  if (old != m_Theta_ij[counter]) {
514  throw CanteraError("HMWSoln::readXMLPsiCommonCation",
515  "conflicting values");
516  }
517  }
518  }
519  if (nodeName == "psi") {
520  getFloatArray(xmlChild, vParams, false, "", stemp);
521  nParamsFound = vParams.size();
522  n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
523 
524  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
525  if (nParamsFound != 1) {
526  throw CanteraError("HMWSoln::readXMLPsiCommonCation::Psi for "
527  + kName + "::" + iName + "::" + jName,
528  "wrong number of params found");
529  }
530  m_Psi_ijk_coeff(0,n) = vParams[0];
531  m_Psi_ijk[n] = vParams[0];
532  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
533  if (nParamsFound != 2) {
534  throw CanteraError("HMWSoln::readXMLPsiCation::Psi for "
535  + kName + "::" + iName + "::" + jName,
536  "wrong number of params found");
537  }
538  m_Psi_ijk_coeff(0,n) = vParams[0];
539  m_Psi_ijk_coeff(1,n) = vParams[1];
540  m_Psi_ijk[n] = vParams[0];
541  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
542  if (nParamsFound == 1) {
543  vParams.resize(5, 0.0);
544  nParamsFound = 5;
545  } else if (nParamsFound != 5) {
546  throw CanteraError("HMWSoln::readXMLPsiCation::Psi for "
547  + kName + "::" + iName + "::" + jName,
548  "wrong number of params found");
549  }
550  for (size_t j = 0; j < nParamsFound; j++) {
551  m_Psi_ijk_coeff(j, n) = vParams[j];
552  }
553  m_Psi_ijk[n] = vParams[0];
554  }
555 
556 
557  // fill in the duplicate entries
558  n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ;
559  for (size_t j = 0; j < nParamsFound; j++) {
560  m_Psi_ijk_coeff(j, n) = vParams[j];
561  }
562  m_Psi_ijk[n] = vParams[0];
563 
564  n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ;
565  for (size_t j = 0; j < nParamsFound; j++) {
566  m_Psi_ijk_coeff(j, n) = vParams[j];
567  }
568  m_Psi_ijk[n] = vParams[0];
569 
570  n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ;
571  for (size_t j = 0; j < nParamsFound; j++) {
572  m_Psi_ijk_coeff(j, n) = vParams[j];
573  }
574  m_Psi_ijk[n] = vParams[0];
575 
576  n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ;
577  for (size_t j = 0; j < nParamsFound; j++) {
578  m_Psi_ijk_coeff(j, n) = vParams[j];
579  }
580  m_Psi_ijk[n] = vParams[0];
581 
582  n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ;
583  for (size_t j = 0; j < nParamsFound; j++) {
584  m_Psi_ijk_coeff(j, n) = vParams[j];
585  }
586  m_Psi_ijk[n] = vParams[0];
587  }
588  }
589 }
590 
591 /**
592  * Process an XML node called "PsiCommonAnion".
593  * This node contains all of the parameters necessary to describe
594  * the binary interactions between two cations and one common anion.
595  */
596 void HMWSoln::readXMLPsiCommonAnion(XML_Node& BinSalt)
597 {
598  string xname = BinSalt.name();
599  if (xname != "psiCommonAnion") {
600  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
601  "Incorrect name for processing this routine: " + xname);
602  }
603  double* charge = DATA_PTR(m_speciesCharge);
604  string stemp;
605  vector_fp vParams;
606  size_t nParamsFound = 0;
607  string kName = BinSalt.attrib("anion");
608  if (kName == "") {
609  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no anion attrib");
610  }
611  string iName = BinSalt.attrib("cation1");
612  if (iName == "") {
613  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation1 attrib");
614  }
615  string jName = BinSalt.attrib("cation2");
616  if (jName == "") {
617  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "no cation2 attrib");
618  }
619  /*
620  * Find the index of the species in the current phase. It's not
621  * an error to not find the species
622  */
623  size_t kSpecies = speciesIndex(kName);
624  if (kSpecies == npos) {
625  return;
626  }
627  if (charge[kSpecies] >= 0) {
628  throw CanteraError("HMWSoln::readXMLPsiCommonAnion", "anion charge problem");
629  }
630  size_t iSpecies = speciesIndex(iName);
631  if (iSpecies == npos) {
632  return;
633  }
634  if (charge[iSpecies] <= 0) {
635  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
636  "cation1 charge problem");
637  }
638  size_t jSpecies = speciesIndex(jName);
639  if (jSpecies == npos) {
640  return;
641  }
642  if (charge[jSpecies] <= 0) {
643  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
644  "cation2 charge problem");
645  }
646 
647  size_t n = iSpecies * m_kk + jSpecies;
648  int counter = m_CounterIJ[n];
649  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
650  XML_Node& xmlChild = BinSalt.child(i);
651  stemp = xmlChild.name();
652  string nodeName = lowercase(stemp);
653  if (nodeName == "theta") {
654  stemp = xmlChild.value();
655  double old = m_Theta_ij[counter];
656  m_Theta_ij[counter] = atofCheck(stemp.c_str());
657  if (old != 0.0) {
658  if (old != m_Theta_ij[counter]) {
659  throw CanteraError("HMWSoln::readXMLPsiCommonAnion",
660  "conflicting values");
661  }
662  }
663  }
664  if (nodeName == "psi") {
665 
666  getFloatArray(xmlChild, vParams, false, "", stemp);
667  nParamsFound = vParams.size();
668  n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
669 
670  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
671  if (nParamsFound != 1) {
672  throw CanteraError("HMWSoln::readXMLPsiCommonAnion::Psi for "
673  + kName + "::" + iName + "::" + jName,
674  "wrong number of params found");
675  }
676  m_Psi_ijk_coeff(0,n) = vParams[0];
677  m_Psi_ijk[n] = vParams[0];
678  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
679  if (nParamsFound != 2) {
680  throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for "
681  + kName + "::" + iName + "::" + jName,
682  "wrong number of params found");
683  }
684  m_Psi_ijk_coeff(0,n) = vParams[0];
685  m_Psi_ijk_coeff(1,n) = vParams[1];
686  m_Psi_ijk[n] = vParams[0];
687  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
688  if (nParamsFound == 1) {
689  vParams.resize(5, 0.0);
690  nParamsFound = 5;
691  } else if (nParamsFound != 5) {
692  throw CanteraError("HMWSoln::readXMLPsiAnion::Psi for "
693  + kName + "::" + iName + "::" + jName,
694  "wrong number of params found");
695  }
696  for (size_t j = 0; j < nParamsFound; j++) {
697  m_Psi_ijk_coeff(j, n) = vParams[j];
698  }
699  m_Psi_ijk[n] = vParams[0];
700  }
701 
702 
703  // fill in the duplicate entries
704  n = iSpecies * m_kk *m_kk + kSpecies * m_kk + jSpecies ;
705  for (size_t j = 0; j < nParamsFound; j++) {
706  m_Psi_ijk_coeff(j, n) = vParams[j];
707  }
708  m_Psi_ijk[n] = vParams[0];
709 
710  n = jSpecies * m_kk *m_kk + iSpecies * m_kk + kSpecies ;
711  for (size_t j = 0; j < nParamsFound; j++) {
712  m_Psi_ijk_coeff(j, n) = vParams[j];
713  }
714  m_Psi_ijk[n] = vParams[0];
715 
716  n = jSpecies * m_kk *m_kk + kSpecies * m_kk + iSpecies ;
717  for (size_t j = 0; j < nParamsFound; j++) {
718  m_Psi_ijk_coeff(j, n) = vParams[j];
719  }
720  m_Psi_ijk[n] = vParams[0];
721 
722  n = kSpecies * m_kk *m_kk + jSpecies * m_kk + iSpecies ;
723  for (size_t j = 0; j < nParamsFound; j++) {
724  m_Psi_ijk_coeff(j, n) = vParams[j];
725  }
726  m_Psi_ijk[n] = vParams[0];
727 
728  n = kSpecies * m_kk *m_kk + iSpecies * m_kk + jSpecies ;
729  for (size_t j = 0; j < nParamsFound; j++) {
730  m_Psi_ijk_coeff(j, n) = vParams[j];
731  }
732  m_Psi_ijk[n] = vParams[0];
733 
734  }
735  }
736 }
737 
738 
739 
740 /**
741  * Process an XML node called "LambdaNeutral".
742  * This node contains all of the parameters necessary to describe
743  * the binary interactions between one neutral species and
744  * any other species (neutral or otherwise) in the mechanism.
745  */
746 void HMWSoln::readXMLLambdaNeutral(XML_Node& BinSalt)
747 {
748  string xname = BinSalt.name();
749  vector_fp vParams;
750  size_t nParamsFound;
751  if (xname != "lambdaNeutral") {
752  throw CanteraError("HMWSoln::readXMLLanbdaNeutral",
753  "Incorrect name for processing this routine: " + xname);
754  }
755  double* charge = DATA_PTR(m_speciesCharge);
756  string stemp;
757  string iName = BinSalt.attrib("species1");
758  if (iName == "") {
759  throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no species1 attrib");
760  }
761  string jName = BinSalt.attrib("species2");
762  if (jName == "") {
763  throw CanteraError("HMWSoln::readXMLLambdaNeutral", "no species2 attrib");
764  }
765  /*
766  * Find the index of the species in the current phase. It's not
767  * an error to not find the species
768  */
769  size_t iSpecies = speciesIndex(iName);
770  if (iSpecies == npos) {
771  return;
772  }
773  if (charge[iSpecies] != 0) {
774  throw CanteraError("HMWSoln::readXMLLambdaNeutral",
775  "neutral charge problem");
776  }
777  size_t jSpecies = speciesIndex(jName);
778  if (jSpecies == npos) {
779  return;
780  }
781 
782  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
783  XML_Node& xmlChild = BinSalt.child(i);
784  stemp = xmlChild.name();
785  string nodeName = lowercase(stemp);
786  if (nodeName == "lambda") {
787  size_t nCount = iSpecies*m_kk + jSpecies;
788  getFloatArray(xmlChild, vParams, false, "", stemp);
789  nParamsFound = vParams.size();
790  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
791  if (nParamsFound != 1) {
792  throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
793  + "::" + jName,
794  "wrong number of params found");
795  }
796  m_Lambda_nj_coeff(0,nCount) = vParams[0];
797  m_Lambda_nj(iSpecies,jSpecies) = vParams[0];
798 
799  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
800  if (nParamsFound != 2) {
801  throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
802  + "::" + jName,
803  "wrong number of params found");
804  }
805  m_Lambda_nj_coeff(0,nCount) = vParams[0];
806  m_Lambda_nj_coeff(1,nCount) = vParams[1];
807  m_Lambda_nj(iSpecies, jSpecies) = vParams[0];
808 
809  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
810  if (nParamsFound == 1) {
811  vParams.resize(5, 0.0);
812  nParamsFound = 5;
813  } else if (nParamsFound != 5) {
814  throw CanteraError("HMWSoln::readXMLLambdaNeutral::Lambda for " + iName
815  + "::" + jName,
816  "wrong number of params found");
817  }
818  for (size_t j = 0; j < nParamsFound; j++) {
819  m_Lambda_nj_coeff(j,nCount) = vParams[j];
820  }
821  m_Lambda_nj(iSpecies, jSpecies) = vParams[0];
822  }
823  }
824  }
825 }
826 
827 /**
828  * Process an XML node called "MunnnNeutral".
829  * This node contains all of the parameters necessary to describe
830  * the self-ternary interactions for one neutral species.
831  */
832 void HMWSoln::readXMLMunnnNeutral(XML_Node& BinSalt)
833 {
834  string xname = BinSalt.name();
835  vector_fp vParams;
836  size_t nParamsFound;
837  if (xname != "MunnnNeutral") {
838  throw CanteraError("HMWSoln::readXMLMunnnNeutral",
839  "Incorrect name for processing this routine: " + xname);
840  }
841  double* charge = DATA_PTR(m_speciesCharge);
842  string stemp;
843  string iName = BinSalt.attrib("species1");
844  if (iName == "") {
845  throw CanteraError("HMWSoln::readXMLMunnnNeutral", "no species1 attrib");
846  }
847 
848  /*
849  * Find the index of the species in the current phase. It's not
850  * an error to not find the species
851  */
852  size_t iSpecies = speciesIndex(iName);
853  if (iSpecies == npos) {
854  return;
855  }
856  if (charge[iSpecies] != 0) {
857  throw CanteraError("HMWSoln::readXMLMunnnNeutral",
858  "neutral charge problem");
859  }
860 
861  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
862  XML_Node& xmlChild = BinSalt.child(i);
863  stemp = xmlChild.name();
864  string nodeName = lowercase(stemp);
865  if (nodeName == "munnn") {
866  getFloatArray(xmlChild, vParams, false, "", "Munnn");
867  nParamsFound = vParams.size();
868  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
869  if (nParamsFound != 1) {
870  throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
871  "wrong number of params found");
872  }
873  m_Mu_nnn_coeff(0,iSpecies) = vParams[0];
874  m_Mu_nnn[iSpecies] = vParams[0];
875 
876  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
877  if (nParamsFound != 2) {
878  throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
879  "wrong number of params found");
880  }
881  m_Mu_nnn_coeff(0, iSpecies) = vParams[0];
882  m_Mu_nnn_coeff(1, iSpecies) = vParams[1];
883  m_Mu_nnn[iSpecies] = vParams[0];
884 
885  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
886  if (nParamsFound == 1) {
887  vParams.resize(5, 0.0);
888  nParamsFound = 5;
889  } else if (nParamsFound != 5) {
890  throw CanteraError("HMWSoln::readXMLMunnnNeutral::Munnn for " + iName,
891  "wrong number of params found");
892  }
893  for (size_t j = 0; j < nParamsFound; j++) {
894  m_Mu_nnn_coeff(j, iSpecies) = vParams[j];
895  }
896  m_Mu_nnn[iSpecies] = vParams[0];
897  }
898  }
899  }
900 }
901 
902 /*
903  * Process an XML node called "readXMLZetaCation".
904  * This node contains all of the parameters necessary to describe
905  * the ternary interactions between a neutral, a cation and an anion
906  */
907 void HMWSoln::readXMLZetaCation(const XML_Node& BinSalt)
908 {
909  string xname = BinSalt.name();
910  if (xname != "zetaCation") {
911  throw CanteraError("HMWSoln::readXMLZetaCation",
912  "Incorrect name for processing this routine: " + xname);
913  }
914  double* charge = DATA_PTR(m_speciesCharge);
915  string stemp;
916  vector_fp vParams;
917  size_t nParamsFound = 0;
918 
919  string iName = BinSalt.attrib("neutral");
920  if (iName == "") {
921  throw CanteraError("HMWSoln::readXMLZetaCation", "no neutral attrib");
922  }
923 
924  string jName = BinSalt.attrib("cation1");
925  if (jName == "") {
926  throw CanteraError("HMWSoln::readXMLZetaCation", "no cation1 attrib");
927  }
928 
929  string kName = BinSalt.attrib("anion1");
930  if (kName == "") {
931  throw CanteraError("HMWSoln::readXMLZetaCation", "no anion1 attrib");
932  }
933  /*
934  * Find the index of the species in the current phase. It's not
935  * an error to not find the species
936  */
937  size_t iSpecies = speciesIndex(iName);
938  if (iSpecies == npos) {
939  return;
940  }
941  if (charge[iSpecies] != 0.0) {
942  throw CanteraError("HMWSoln::readXMLZetaCation", "neutral charge problem");
943  }
944 
945  size_t jSpecies = speciesIndex(jName);
946  if (jSpecies == npos) {
947  return;
948  }
949  if (charge[jSpecies] <= 0.0) {
950  throw CanteraError("HMWSoln::readXLZetaCation", "cation1 charge problem");
951  }
952 
953  size_t kSpecies = speciesIndex(kName);
954  if (kSpecies == npos) {
955  return;
956  }
957  if (charge[kSpecies] >= 0.0) {
958  throw CanteraError("HMWSoln::readXMLZetaCation", "anion1 charge problem");
959  }
960 
961  for (size_t i = 0; i < BinSalt.nChildren(); i++) {
962  XML_Node& xmlChild = BinSalt.child(i);
963  stemp = xmlChild.name();
964  string nodeName = lowercase(stemp);
965  if (nodeName == "zeta") {
966  getFloatArray(xmlChild, vParams, false, "", "zeta");
967  nParamsFound = vParams.size();
968  size_t n = iSpecies * m_kk *m_kk + jSpecies * m_kk + kSpecies ;
969 
970  if (m_formPitzerTemp == PITZER_TEMP_CONSTANT) {
971  if (nParamsFound != 1) {
972  throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
973  + iName + "::" + jName + "::" + kName,
974  "wrong number of params found");
975  }
976  m_Psi_ijk_coeff(0,n) = vParams[0];
977  m_Psi_ijk[n] = vParams[0];
978  } else if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
979  if (nParamsFound != 2) {
980  throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
981  + iName + "::" + jName + "::" + kName,
982  "wrong number of params found");
983  }
984  m_Psi_ijk_coeff(0,n) = vParams[0];
985  m_Psi_ijk_coeff(1,n) = vParams[1];
986  m_Psi_ijk[n] = vParams[0];
987  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
988  if (nParamsFound == 1) {
989  vParams.resize(5, 0.0);
990  nParamsFound = 5;
991  } else if (nParamsFound != 5) {
992  throw CanteraError("HMWSoln::readXMLZetaCation::Zeta for "
993  + iName + "::" + jName + "::" + kName,
994  "wrong number of params found");
995  }
996  for (size_t j = 0; j < nParamsFound; j++) {
997  m_Psi_ijk_coeff(j, n) = vParams[j];
998  }
999  m_Psi_ijk[n] = vParams[0];
1000  }
1001 
1002  // There are no duplicate entries
1003  }
1004  }
1005 }
1006 
1007 // Process an XML node called "croppingCoefficients"
1008 // for the cropping coefficients values
1009 /*
1010  * @param acNode Activity Coefficient XML Node
1011  */
1012 void HMWSoln::readXMLCroppingCoefficients(const XML_Node& acNode)
1013 {
1014 
1015  if (acNode.hasChild("croppingCoefficients")) {
1016  XML_Node& cropNode = acNode.child("croppingCoefficients");
1017  if (cropNode.hasChild("ln_gamma_k_min")) {
1018  XML_Node& gkminNode = cropNode.child("ln_gamma_k_min");
1019  getOptionalFloat(gkminNode, "pureSolventValue", CROP_ln_gamma_k_min);
1020  }
1021  if (cropNode.hasChild("ln_gamma_k_max")) {
1022  XML_Node& gkmaxNode = cropNode.child("ln_gamma_k_max");
1023  getOptionalFloat(gkmaxNode, "pureSolventValue", CROP_ln_gamma_k_max);
1024  }
1025 
1026  if (cropNode.hasChild("ln_gamma_o_min")) {
1027  XML_Node& gominNode = cropNode.child("ln_gamma_o_min");
1028  getOptionalFloat(gominNode, "pureSolventValue", CROP_ln_gamma_o_min);
1029  }
1030 
1031  if (cropNode.hasChild("ln_gamma_o_max")) {
1032  XML_Node& gomaxNode = cropNode.child("ln_gamma_o_max");
1033  getOptionalFloat(gomaxNode, "pureSolventValue", CROP_ln_gamma_o_max);
1034  }
1035  }
1036 }
1037 
1038 /*
1039  * Initialization routine for a HMWSoln phase.
1040  *
1041  * This is a virtual routine. This routine will call initThermo()
1042  * for the parent class as well.
1043  */
1044 void HMWSoln::initThermo()
1045 {
1046  MolalityVPSSTP::initThermo();
1047  initLengths();
1048 }
1049 
1050 /*
1051  * Import, construct, and initialize a HMWSoln phase
1052  * specification from an XML tree into the current object.
1053  *
1054  * This routine is a precursor to constructPhaseXML(XML_Node*)
1055  * routine, which does most of the work.
1056  *
1057  * @param infile XML file containing the description of the
1058  * phase
1059  *
1060  * @param id Optional parameter identifying the name of the
1061  * phase. If none is given, the first XML
1062  * phase element will be used.
1063  */
1064 void HMWSoln::constructPhaseFile(std::string inputFile, std::string id)
1065 {
1066 
1067  if (inputFile.size() == 0) {
1068  throw CanteraError("HMWSoln:constructPhaseFile",
1069  "input file is null");
1070  }
1071  string path = findInputFile(inputFile);
1072  std::ifstream fin(path.c_str());
1073  if (!fin) {
1074  throw CanteraError("HMWSoln:constructPhaseFile","could not open "
1075  +path+" for reading.");
1076  }
1077  /*
1078  * The phase object automatically constructs an XML object.
1079  * Use this object to store information.
1080  */
1081  XML_Node& phaseNode_XML = xml();
1082  XML_Node* fxml = new XML_Node();
1083  fxml->build(fin);
1084  XML_Node* fxml_phase = findXMLPhase(fxml, id);
1085  if (!fxml_phase) {
1086  throw CanteraError("HMWSoln:constructPhaseFile",
1087  "ERROR: Can not find phase named " +
1088  id + " in file named " + inputFile);
1089  }
1090  fxml_phase->copy(&phaseNode_XML);
1091  constructPhaseXML(*fxml_phase, id);
1092  delete fxml;
1093 }
1094 
1095 /*
1096  * Import, construct, and initialize a HMWSoln phase
1097  * specification from an XML tree into the current object.
1098  *
1099  * Most of the work is carried out by the cantera base
1100  * routine, importPhase(). That routine imports all of the
1101  * species and element data, including the standard states
1102  * of the species.
1103  *
1104  * Then, In this routine, we read the information
1105  * particular to the specification of the activity
1106  * coefficient model for the Pitzer parameterization.
1107  *
1108  * We also read information about the molar volumes of the
1109  * standard states if present in the XML file.
1110  *
1111  * @param phaseNode This object must be the phase node of a
1112  * complete XML tree
1113  * description of the phase, including all of the
1114  * species data. In other words while "phase" must
1115  * point to an XML phase object, it must have
1116  * sibling nodes "speciesData" that describe
1117  * the species in the phase.
1118  * @param id ID of the phase. If nonnull, a check is done
1119  * to see if phaseNode is pointing to the phase
1120  * with the correct id.
1121  */
1122 void HMWSoln::constructPhaseXML(XML_Node& phaseNode, std::string id)
1123 {
1124  string stemp;
1125  if (id.size() > 0) {
1126  string idp = phaseNode.id();
1127  if (idp != id) {
1128  throw CanteraError("HMWSoln::constructPhaseXML",
1129  "phasenode and Id are incompatible");
1130  }
1131  }
1132 
1133  /*
1134  * Find the Thermo XML node
1135  */
1136  if (!phaseNode.hasChild("thermo")) {
1137  throw CanteraError("HMWSoln::constructPhaseXML",
1138  "no thermo XML node");
1139  }
1140  XML_Node& thermoNode = phaseNode.child("thermo");
1141 
1142  /*
1143  * Possibly change the form of the standard concentrations
1144  */
1145  if (thermoNode.hasChild("standardConc")) {
1146  XML_Node& scNode = thermoNode.child("standardConc");
1147  m_formGC = 2;
1148  stemp = scNode.attrib("model");
1149  string formString = lowercase(stemp);
1150  if (formString != "") {
1151  if (formString == "unity") {
1152  m_formGC = 0;
1153  printf("exit standardConc = unity not done\n");
1154  exit(EXIT_FAILURE);
1155  } else if (formString == "molar_volume") {
1156  m_formGC = 1;
1157  printf("exit standardConc = molar_volume not done\n");
1158  exit(EXIT_FAILURE);
1159  } else if (formString == "solvent_volume") {
1160  m_formGC = 2;
1161  } else {
1162  throw CanteraError("HMWSoln::constructPhaseXML",
1163  "Unknown standardConc model: " + formString);
1164  }
1165  }
1166  }
1167  /*
1168  * Get the Name of the Solvent:
1169  * <solvent> solventName </solvent>
1170  */
1171  string solventName = "";
1172  if (thermoNode.hasChild("solvent")) {
1173  XML_Node& scNode = thermoNode.child("solvent");
1174  vector<string> nameSolventa;
1175  getStringArray(scNode, nameSolventa);
1176  int nsp = static_cast<int>(nameSolventa.size());
1177  if (nsp != 1) {
1178  throw CanteraError("HMWSoln::constructPhaseXML",
1179  "badly formed solvent XML node");
1180  }
1181  solventName = nameSolventa[0];
1182  }
1183 
1184  /*
1185  * Determine the form of the Pitzer model,
1186  * We will use this information to size arrays below.
1187  */
1188  if (thermoNode.hasChild("activityCoefficients")) {
1189  XML_Node& scNode = thermoNode.child("activityCoefficients");
1190  stemp = scNode.attrib("model");
1191  string formString = lowercase(stemp);
1192  if (formString != "") {
1193  if (formString == "pitzer" || formString == "default") {
1194  m_formPitzer = PITZERFORM_BASE;
1195  } else if (formString == "base") {
1196  m_formPitzer = PITZERFORM_BASE;
1197  } else {
1198  throw CanteraError("HMWSoln::constructPhaseXML",
1199  "Unknown Pitzer ActivityCoeff model: "
1200  + formString);
1201  }
1202  }
1203 
1204  /*
1205  * Determine the form of the temperature dependence
1206  * of the Pitzer activity coefficient model.
1207  */
1208  stemp = scNode.attrib("TempModel");
1209  formString = lowercase(stemp);
1210  if (formString != "") {
1211  if (formString == "constant" || formString == "default") {
1212  m_formPitzerTemp = PITZER_TEMP_CONSTANT;
1213  } else if (formString == "linear") {
1214  m_formPitzerTemp = PITZER_TEMP_LINEAR;
1215  } else if (formString == "complex" || formString == "complex1") {
1216  m_formPitzerTemp = PITZER_TEMP_COMPLEX1;
1217  } else {
1218  throw CanteraError("HMWSoln::constructPhaseXML",
1219  "Unknown Pitzer ActivityCoeff Temp model: "
1220  + formString);
1221  }
1222  }
1223 
1224  /*
1225  * Determine the reference temperature
1226  * of the Pitzer activity coefficient model's temperature
1227  * dependence formulation: defaults to 25C
1228  */
1229  stemp = scNode.attrib("TempReference");
1230  formString = lowercase(stemp);
1231  if (formString != "") {
1232  m_TempPitzerRef = atofCheck(formString.c_str());
1233  } else {
1234  m_TempPitzerRef = 273.15 + 25;
1235  }
1236 
1237  }
1238 
1239  /*
1240  * Call the Cantera importPhase() function. This will import
1241  * all of the species into the phase. This will also handle
1242  * all of the solvent and solute standard states
1243  */
1244  bool m_ok = importPhase(phaseNode, this);
1245  if (!m_ok) {
1246  throw CanteraError("HMWSoln::constructPhaseXML","importPhase failed ");
1247  }
1248 
1249 }
1250 
1251 /**
1252  * Process the XML file after species are set up.
1253  *
1254  * This gets called from importPhase(). It processes the XML file
1255  * after the species are set up. This is the main routine for
1256  * reading in activity coefficient parameters.
1257  *
1258  * @param phaseNode This object must be the phase node of a
1259  * complete XML tree
1260  * description of the phase, including all of the
1261  * species data. In other words while "phase" must
1262  * point to an XML phase object, it must have
1263  * sibling nodes "speciesData" that describe
1264  * the species in the phase.
1265  * @param id ID of the phase. If nonnull, a check is done
1266  * to see if phaseNode is pointing to the phase
1267  * with the correct id.
1268  */
1269 void HMWSoln::
1270 initThermoXML(XML_Node& phaseNode, std::string id)
1271 {
1272  string stemp;
1273  /*
1274  * Find the Thermo XML node
1275  */
1276  if (!phaseNode.hasChild("thermo")) {
1277  throw CanteraError("HMWSoln::initThermoXML",
1278  "no thermo XML node");
1279  }
1280  XML_Node& thermoNode = phaseNode.child("thermo");
1281 
1282  /*
1283  * Get the Name of the Solvent:
1284  * <solvent> solventName </solvent>
1285  */
1286  string solventName = "";
1287  if (thermoNode.hasChild("solvent")) {
1288  XML_Node& scNode = thermoNode.child("solvent");
1289  vector<string> nameSolventa;
1290  getStringArray(scNode, nameSolventa);
1291  int nsp = static_cast<int>(nameSolventa.size());
1292  if (nsp != 1) {
1293  throw CanteraError("HMWSoln::initThermoXML",
1294  "badly formed solvent XML node");
1295  }
1296  solventName = nameSolventa[0];
1297  }
1298 
1299  /*
1300  * Initialize all of the lengths of arrays in the object
1301  * now that we know what species are in the phase.
1302  */
1303  initLengths();
1304 
1305  /*
1306  * Reconcile the solvent name and index.
1307  */
1308  for (size_t k = 0; k < m_kk; k++) {
1309  string sname = speciesName(k);
1310  if (solventName == sname) {
1311  setSolvent(k);
1312  if (k != 0) {
1313  throw CanteraError("HMWSoln::initThermoXML",
1314  "Solvent must be species 0 atm");
1315  }
1316  m_indexSolvent = k;
1317  break;
1318  }
1319  }
1320  if (m_indexSolvent == npos) {
1321  std::cout << "HMWSoln::initThermo: Solvent Name not found"
1322  << std::endl;
1323  throw CanteraError("HMWSoln::initThermoXML",
1324  "Solvent name not found");
1325  }
1326  if (m_indexSolvent != 0) {
1327  throw CanteraError("HMWSoln::initThermoXML",
1328  "Solvent " + solventName +
1329  " should be first species");
1330  }
1331 
1332  /*
1333  * Now go get the specification of the standard states for
1334  * species in the solution. This includes the molar volumes
1335  * data blocks for incompressible species.
1336  */
1337  XML_Node& speciesList = phaseNode.child("speciesArray");
1338  XML_Node* speciesDB =
1339  get_XML_NameID("speciesData", speciesList["datasrc"],
1340  &phaseNode.root());
1341  const vector<string>&sss = speciesNames();
1342 
1343  for (size_t k = 0; k < m_kk; k++) {
1344  XML_Node* s = speciesDB->findByAttr("name", sss[k]);
1345  if (!s) {
1346  throw CanteraError("HMWSoln::initThermoXML",
1347  "Species Data Base " + sss[k] + " not found");
1348  }
1349  XML_Node* ss = s->findByName("standardState");
1350  if (!ss) {
1351  throw CanteraError("HMWSoln::initThermoXML",
1352  "Species " + sss[k] +
1353  " standardState XML block not found");
1354  }
1355  string modelStringa = ss->attrib("model");
1356  if (modelStringa == "") {
1357  throw CanteraError("HMWSoln::initThermoXML",
1358  "Species " + sss[k] +
1359  " standardState XML block model attribute not found");
1360  }
1361  string modelString = lowercase(modelStringa);
1362  if (k == 0) {
1363  if (modelString == "wateriapws" || modelString == "real_water" ||
1364  modelString == "waterpdss") {
1365  /*
1366  * Store a local pointer to the water standard state model.
1367  * -> We've hardcoded it to a PDSS_Water model, so this is ok.
1368  */
1369  m_waterSS = dynamic_cast<PDSS_Water*>(providePDSS(0)) ;
1370  if (!m_waterSS) {
1371  throw CanteraError("HMWSoln::initThermoXML",
1372  "Dynamic cast to PDSS_Water failed");
1373  }
1374  /*
1375  * Fill in the molar volume of water (m3/kmol)
1376  * at standard conditions to fill in the m_speciesSize entry
1377  * with something reasonable.
1378  */
1379  m_waterSS->setState_TP(300., OneAtm);
1380  double dens = m_waterSS->density();
1381  double mw = m_waterSS->molecularWeight();
1382  m_speciesSize[0] = mw / dens;
1383 #ifdef DEBUG_HKM_NOT
1384  cout << "Solvent species " << sss[k] << " has volume " <<
1385  m_speciesSize[k] << endl;
1386 #endif
1387  } else {
1388  // throw CanteraError("HMWSoln::initThermoXML",
1389  // "Solvent SS Model \"" + modelStringa +
1390  // "\" is not allowed, name = " + sss[0]);
1391  m_waterSS = providePDSS(0);
1392  m_waterSS->setState_TP(300., OneAtm);
1393  double dens = m_waterSS->density();
1394  double mw = m_waterSS->molecularWeight();
1395  m_speciesSize[0] = mw / dens;
1396  }
1397  } else {
1398  if (modelString != "constant_incompressible" && modelString != "hkft") {
1399  throw CanteraError("HMWSoln::initThermoXML",
1400  "Solute SS Model \"" + modelStringa +
1401  "\" is not known");
1402  }
1403  if (modelString == "constant_incompressible") {
1404  m_speciesSize[k] = getFloat(*ss, "molarVolume", "toSI");
1405 #ifdef DEBUG_HKM_NOT
1406  cout << "species " << sss[k] << " has volume " <<
1407  m_speciesSize[k] << endl;
1408 #endif
1409  }
1410  // HKM Note, have to fill up m_speciesSize[] for HKFT species
1411  }
1412  }
1413 
1414  /*
1415  * Initialize the water property calculator. It will share
1416  * the internal eos water calculator.
1417  */
1418  m_waterProps = new WaterProps(dynamic_cast<PDSS_Water*>(m_waterSS));
1419 
1420  /*
1421  * Fill in parameters for the calculation of the
1422  * stoichiometric Ionic Strength
1423  *
1424  * The default is that stoich charge is the same as the
1425  * regular charge.
1426  */
1427  for (size_t k = 0; k < m_kk; k++) {
1428  m_speciesCharge_Stoich[k] = m_speciesCharge[k];
1429  }
1430 
1431  /*
1432  * Go get all of the coefficients and factors in the
1433  * activityCoefficients XML block
1434  */
1435  XML_Node* acNodePtr = 0;
1436  if (thermoNode.hasChild("activityCoefficients")) {
1437  XML_Node& acNode = thermoNode.child("activityCoefficients");
1438  acNodePtr = &acNode;
1439  /*
1440  * Look for parameters for A_Debye
1441  */
1442  if (acNode.hasChild("A_Debye")) {
1443  XML_Node& ADebye = acNode.child("A_Debye");
1444  m_form_A_Debye = A_DEBYE_CONST;
1445  stemp = "model";
1446  if (ADebye.hasAttrib(stemp)) {
1447  string atemp = ADebye.attrib(stemp);
1448  stemp = lowercase(atemp);
1449  if (stemp == "water") {
1450  m_form_A_Debye = A_DEBYE_WATER;
1451  }
1452  }
1453  if (m_form_A_Debye == A_DEBYE_CONST) {
1454  m_A_Debye = getFloat(acNode, "A_Debye");
1455  }
1456 #ifdef DEBUG_HKM_NOT
1457  cout << "A_Debye = " << m_A_Debye << endl;
1458 #endif
1459  }
1460 
1461  /*
1462  * Look for Parameters for the Maximum Ionic Strength
1463  */
1464  if (acNode.hasChild("maxIonicStrength")) {
1465  m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
1466 #ifdef DEBUG_HKM_NOT
1467  cout << "m_maxIionicStrength = "
1468  <<m_maxIionicStrength << endl;
1469 #endif
1470  }
1471 
1472 
1473  /*
1474  * Look for parameters for the Ionic radius
1475  */
1476  if (acNode.hasChild("ionicRadius")) {
1477  XML_Node& irNode = acNode.child("ionicRadius");
1478 
1479  string Aunits = "";
1480  double Afactor = 1.0;
1481  if (irNode.hasAttrib("units")) {
1482  string Aunits = irNode.attrib("units");
1483  Afactor = toSI(Aunits);
1484  }
1485 
1486  if (irNode.hasAttrib("default")) {
1487  string ads = irNode.attrib("default");
1488  double ad = fpValue(ads);
1489  for (size_t k = 0; k < m_kk; k++) {
1490  m_Aionic[k] = ad * Afactor;
1491  }
1492  }
1493 
1494  }
1495 
1496  /*
1497  * First look at the species database.
1498  * -> Look for the subelement "stoichIsMods"
1499  * in each of the species SS databases.
1500  */
1501  std::vector<const XML_Node*> xspecies = speciesData();
1502 
1503  string kname, jname;
1504  size_t jj = xspecies.size();
1505  for (size_t k = 0; k < m_kk; k++) {
1506  size_t jmap = npos;
1507  kname = speciesName(k);
1508  for (size_t j = 0; j < jj; j++) {
1509  const XML_Node& sp = *xspecies[j];
1510  jname = sp["name"];
1511  if (jname == kname) {
1512  jmap = j;
1513  break;
1514  }
1515  }
1516  if (jmap != npos) {
1517  const XML_Node& sp = *xspecies[jmap];
1518  getOptionalFloat(sp, "stoichIsMods", m_speciesCharge_Stoich[k]);
1519  // if (sp.hasChild("stoichIsMods")) {
1520  // double val = getFloat(sp, "stoichIsMods");
1521  //m_speciesCharge_Stoich[k] = val;
1522  //}
1523  }
1524  }
1525 
1526  /*
1527  * Now look at the activity coefficient database
1528  */
1529  if (acNodePtr) {
1530  if (acNodePtr->hasChild("stoichIsMods")) {
1531  XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
1532 
1533  map<string, string> msIs;
1534  getMap(sIsNode, msIs);
1535  map<string,string>::const_iterator _b = msIs.begin();
1536  for (; _b != msIs.end(); ++_b) {
1537  size_t kk = speciesIndex(_b->first);
1538  if (kk != npos) {
1539  double val = fpValue(_b->second);
1540  m_speciesCharge_Stoich[kk] = val;
1541  }
1542  }
1543  }
1544  }
1545 
1546 
1547  /*
1548  * Loop through the children getting multiple instances of
1549  * parameters
1550  */
1551  if (acNodePtr) {
1552  for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
1553  XML_Node& xmlACChild = acNodePtr->child(i);
1554  stemp = xmlACChild.name();
1555  string nodeName = lowercase(stemp);
1556  /*
1557  * Process a binary salt field, or any of the other XML fields
1558  * that make up the Pitzer Database. Entries will be ignored
1559  * if any of the species in the entry isn't in the solution.
1560  */
1561  if (nodeName == "binarysaltparameters") {
1562  readXMLBinarySalt(xmlACChild);
1563  } else if (nodeName == "thetaanion") {
1564  readXMLThetaAnion(xmlACChild);
1565  } else if (nodeName == "thetacation") {
1566  readXMLThetaCation(xmlACChild);
1567  } else if (nodeName == "psicommonanion") {
1568  readXMLPsiCommonAnion(xmlACChild);
1569  } else if (nodeName == "psicommoncation") {
1570  readXMLPsiCommonCation(xmlACChild);
1571  } else if (nodeName == "lambdaneutral") {
1572  readXMLLambdaNeutral(xmlACChild);
1573  } else if (nodeName == "zetacation") {
1574  readXMLZetaCation(xmlACChild);
1575  }
1576  }
1577  }
1578 
1579  // Go look up the optional Cropping parameters
1580  readXMLCroppingCoefficients(acNode);
1581 
1582  }
1583 
1584  /*
1585  * Fill in the vector specifying the electrolyte species
1586  * type
1587  *
1588  * First fill in default values. Everthing is either
1589  * a charge species, a nonpolar neutral, or the solvent.
1590  */
1591  for (size_t k = 0; k < m_kk; k++) {
1592  if (fabs(m_speciesCharge[k]) > 0.0001) {
1593  m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
1594  if (fabs(m_speciesCharge_Stoich[k] - m_speciesCharge[k])
1595  > 0.0001) {
1596  m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
1597  }
1598  } else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
1599  m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
1600  } else {
1601  m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
1602  }
1603  }
1604  m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
1605  /*
1606  * First look at the species database.
1607  * -> Look for the subelement "stoichIsMods"
1608  * in each of the species SS databases.
1609  */
1610  std::vector<const XML_Node*> xspecies = speciesData();
1611  const XML_Node* spPtr = 0;
1612  string kname;
1613  for (size_t k = 0; k < m_kk; k++) {
1614  kname = speciesName(k);
1615  spPtr = xspecies[k];
1616  if (!spPtr) {
1617  if (spPtr->hasChild("electrolyteSpeciesType")) {
1618  string est = getChildValue(*spPtr, "electrolyteSpeciesType");
1619  if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
1620  throw CanteraError("HMWSoln::initThermoXML",
1621  "Bad electrolyte type: " + est);
1622  }
1623  }
1624  }
1625  }
1626  /*
1627  * Then look at the phase thermo specification
1628  */
1629  if (acNodePtr) {
1630  if (acNodePtr->hasChild("electrolyteSpeciesType")) {
1631  XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
1632  map<string, string> msEST;
1633  getMap(ESTNode, msEST);
1634  map<string,string>::const_iterator _b = msEST.begin();
1635  for (; _b != msEST.end(); ++_b) {
1636  size_t kk = speciesIndex(_b->first);
1637  if (kk != npos) {
1638  string est = _b->second;
1639  if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
1640  throw CanteraError("HMWSoln::initThermoXML",
1641  "Bad electrolyte type: " + est);
1642  }
1643  }
1644  }
1645  }
1646  }
1647 
1648  IMS_typeCutoff_ = 2;
1649  if (IMS_typeCutoff_ == 2) {
1650  calcIMSCutoffParams_();
1651  }
1652  calcMCCutoffParams_();
1653  setMoleFSolventMin(1.0E-5);
1654 
1655  MolalityVPSSTP::initThermoXML(phaseNode, id);
1656  /*
1657  * Lastly calculate the charge balance and then add stuff until the charges compensate
1658  */
1659 
1660  vector_fp mf(m_kk, 0.0);
1661  getMoleFractions(DATA_PTR(mf));
1662  bool notDone = true;
1663 
1664  do {
1665  double sum = 0.0;
1666  size_t kMaxC = npos;
1667  double MaxC = 0.0;
1668  for (size_t k = 0; k < m_kk; k++) {
1669  sum += mf[k] * m_speciesCharge[k];
1670  if (fabs(mf[k] * m_speciesCharge[k]) > MaxC) {
1671  kMaxC = k;
1672  }
1673  }
1674  size_t kHp = speciesIndex("H+");
1675  size_t kOHm = speciesIndex("OH-");
1676 
1677 
1678  if (fabs(sum) > 1.0E-30) {
1679  if (kHp != npos) {
1680  if (mf[kHp] > sum * 1.1) {
1681  mf[kHp] -= sum;
1682  mf[0] += sum;
1683  notDone = false;
1684  } else {
1685  if (sum > 0.0) {
1686  mf[kHp] *= 0.5;
1687  mf[0] += mf[kHp];
1688  sum -= mf[kHp];
1689  }
1690  }
1691  }
1692  if (notDone) {
1693  if (kOHm != npos) {
1694  if (mf[kOHm] > -sum * 1.1) {
1695  mf[kOHm] += sum;
1696  mf[0] -= sum;
1697  notDone = false;
1698  } else {
1699  if (sum < 0.0) {
1700  mf[kOHm] *= 0.5;
1701  mf[0] += mf[kOHm];
1702  sum += mf[kOHm];
1703  }
1704  }
1705  }
1706  if (notDone) {
1707  if (kMaxC != npos) {
1708  if (mf[kMaxC] > (1.1 * sum / m_speciesCharge[kMaxC])) {
1709  mf[kMaxC] -= sum / m_speciesCharge[kMaxC];
1710  mf[0] += sum / m_speciesCharge[kMaxC];
1711  } else {
1712  mf[kMaxC] *= 0.5;
1713  mf[0] += mf[kMaxC];
1714  notDone = true;
1715  }
1716  }
1717  }
1718  }
1719  setMoleFractions(DATA_PTR(mf));
1720  } else {
1721  notDone = false;
1722  }
1723  } while (notDone);
1724 
1725 
1726 
1727 
1728 
1729 
1730 
1731 
1732  // if (phaseNode.hasChild("state")) {
1733  // XML_Node& stateNode = phaseNode.child("state");
1734  // setStateFromXML(stateNode);
1735  //}
1736 
1737 }
1738 //====================================================================================================================
1739 // Precalculate the IMS Cutoff parameters for typeCutoff = 2
1740 void HMWSoln::calcIMSCutoffParams_()
1741 {
1742  IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
1743  IMS_efCut_ = 0.0;
1744  bool converged = false;
1745  double oldV = 0.0;
1746  int its;
1747  for (its = 0; its < 100 && !converged; its++) {
1748  oldV = IMS_efCut_;
1749  IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
1750  IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
1751  IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
1752  /
1753  (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
1754  double tmp = IMS_afCut_ + IMS_X_o_cutoff_*(IMS_bfCut_ + IMS_dfCut_ *IMS_X_o_cutoff_);
1755  double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1756  IMS_efCut_ = - eterm * (tmp);
1757  if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
1758  converged = true;
1759  }
1760  }
1761  if (!converged) {
1762  throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1763  " failed to converge on the f polynomial");
1764  }
1765  converged = false;
1766  double f_0 = IMS_afCut_ + IMS_efCut_;
1767  double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
1768  IMS_egCut_ = 0.0;
1769  for (its = 0; its < 100 && !converged; its++) {
1770  oldV = IMS_egCut_;
1771  double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
1772  IMS_agCut_ = exp(lng_0) - IMS_egCut_;
1773  IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
1774  IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
1775  /
1776  (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
1777  double tmp = IMS_agCut_ + IMS_X_o_cutoff_*(IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
1778  double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
1779  IMS_egCut_ = - eterm * (tmp);
1780  if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
1781  converged = true;
1782  }
1783  }
1784  if (!converged) {
1785  throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
1786  " failed to converge on the g polynomial");
1787  }
1788 }
1789 
1790 // Precalculate the MC Cutoff parameters
1791 void HMWSoln::calcMCCutoffParams_()
1792 {
1793  MC_X_o_min_ = 0.35;
1794  MC_X_o_cutoff_ = 0.6;
1795  MC_slopepCut_ = 0.02;
1796  MC_cpCut_ = 0.25;
1797 
1798  // Initial starting values
1799  MC_apCut_ = MC_X_o_min_;
1800  MC_epCut_ = 0.0;
1801  bool converged = false;
1802  double oldV = 0.0;
1803  int its;
1804  double damp = 0.5;
1805  for (its = 0; its < 500 && !converged; its++) {
1806  oldV = MC_epCut_;
1807  MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
1808  double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;
1809  MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
1810  double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ * MC_X_o_cutoff_/MC_cpCut_)
1811  /
1812  (MC_X_o_cutoff_ * MC_X_o_cutoff_/MC_cpCut_ - 2.0 * MC_X_o_cutoff_));
1813  MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
1814  double tmp = MC_apCut_ + MC_X_o_cutoff_*(MC_bpCut_ + MC_dpCut_ * MC_X_o_cutoff_);
1815  double eterm = std::exp(- MC_X_o_cutoff_ / MC_cpCut_);
1816  MC_epCut_ = - eterm * (tmp);
1817  double diff = MC_epCut_ - oldV;
1818  if (fabs(diff) < 1.0E-14) {
1819  converged = true;
1820  }
1821  }
1822  if (!converged) {
1823  throw CanteraError("HMWSoln::calcMCCutoffParams_()",
1824  " failed to converge on the p polynomial");
1825  }
1826 }
1827 
1828 }
Generated by   doxygen 1.8.2