CoolProp  6.6.0
An open-source fluid property and humid air property database
MixtureParameters.cpp
Go to the documentation of this file.
1 #include "MixtureParameters.h"
2 #include "CPstrings.h"
3 #include "mixture_departure_functions_JSON.h" // Creates the variable mixture_departure_functions_JSON
4 #include "mixture_binary_pairs_JSON.h" // Creates the variable mixture_binary_pairs_JSON
5 #include "predefined_mixtures_JSON.h" // Makes a std::string variable called predefined_mixtures_JSON
6 
7 namespace CoolProp {
8 
13 class PredefinedMixturesLibrary
14 {
15  public:
16  std::map<std::string, Dictionary> predefined_mixture_map;
17 
18  PredefinedMixturesLibrary() {
19  load_from_string(predefined_mixtures_JSON);
20  }
21 
22  void load_from_string(const std::string& str) {
23  rapidjson::Document doc;
24  doc.Parse<0>(str.c_str());
25  if (doc.HasParseError()) {
26  std::cout << str << std::endl;
27  throw ValueError("Unable to parse predefined mixture string");
28  }
29  load_from_JSON(doc);
30  }
31 
32  void load_from_JSON(rapidjson::Document & doc){
33  if (!doc.IsArray() || !doc[0].IsObject()){
34  throw ValueError("You must provide an array of objects");
35  }
36  // Iterate over the papers in the listing
37  for (rapidjson::Value::ValueIterator itr = doc.Begin(); itr != doc.End(); ++itr) {
38  // Instantiate the empty dictionary to be filled
39  Dictionary dict;
40  // Get the name
41  std::string name = cpjson::get_string(*itr, "name") + ".mix";
42  // Get the fluid names
43  dict.add_string_vector("fluids", cpjson::get_string_array(*itr, "fluids"));
44  // Get the mole fractions
45  dict.add_double_vector("mole_fractions", cpjson::get_double_array(*itr, "mole_fractions"));
46  // Add to the map
47  predefined_mixture_map.insert(std::pair<std::string, Dictionary>(name, dict));
48  // Also add the uppercase version to the map
49  predefined_mixture_map.insert(std::pair<std::string, Dictionary>(upper(name), dict));
50  }
51  }
52 };
53 static PredefinedMixturesLibrary predefined_mixtures_library;
54 
55 std::string get_csv_predefined_mixtures() {
56  std::vector<std::string> out;
57  for (std::map<std::string, Dictionary>::const_iterator it = predefined_mixtures_library.predefined_mixture_map.begin();
58  it != predefined_mixtures_library.predefined_mixture_map.end(); ++it) {
59  out.push_back(it->first);
60  }
61  return strjoin(out, ",");
62 }
63 
64 bool is_predefined_mixture(const std::string& name, Dictionary& dict) {
65  std::map<std::string, Dictionary>::const_iterator iter = predefined_mixtures_library.predefined_mixture_map.find(name);
66  if (iter != predefined_mixtures_library.predefined_mixture_map.end()) {
67  dict = iter->second;
68  return true;
69  } else {
70  return false;
71  }
72 }
73 
74 void set_predefined_mixtures(const std::string& string_data) {
75  // JSON-encoded string for binary interaction parameters
76  predefined_mixtures_library.load_from_string(string_data);
77 }
78 
83 class MixtureBinaryPairLibrary
84 {
85  private:
87  std::map<std::vector<std::string>, std::vector<Dictionary>> m_binary_pair_map;
88 
89  public:
90  std::map<std::vector<std::string>, std::vector<Dictionary>>& binary_pair_map() {
91  // Set the default departure functions if none have been provided yet
92  if (m_binary_pair_map.size() == 0) {
93  load_defaults();
94  }
95  return m_binary_pair_map;
96  };
97 
98  void load_from_string(const std::string& str) {
99  rapidjson::Document doc;
100  doc.Parse<0>(str.c_str());
101  if (doc.HasParseError()) {
102  std::cout << str << std::endl;
103  throw ValueError("Unable to parse binary interaction function string");
104  }
105  load_from_JSON(doc);
106  }
107 
108  // Load the defaults that come from the JSON-encoded string compiled into library
109  // as the variable mixture_departure_functions_JSON
110  void load_defaults() {
111  load_from_string(mixture_binary_pairs_JSON);
112  }
113 
118  void load_from_JSON(rapidjson::Document& doc) {
119 
120  // Iterate over the papers in the listing
121  for (rapidjson::Value::ValueIterator itr = doc.Begin(); itr != doc.End(); ++itr) {
122  // Get the empty dictionary to be filled by the appropriate reducing parameter filling function
123  Dictionary dict;
124 
125  // Get the vector of CAS numbers
126  std::vector<std::string> CAS;
127  CAS.push_back(cpjson::get_string(*itr, "CAS1"));
128  CAS.push_back(cpjson::get_string(*itr, "CAS2"));
129  std::string name1 = cpjson::get_string(*itr, "Name1");
130  std::string name2 = cpjson::get_string(*itr, "Name2");
131 
132  // Sort the CAS number vector
133  std::sort(CAS.begin(), CAS.end());
134 
135  // A sort was carried out, names/CAS were swapped
136  bool swapped = CAS[0].compare(cpjson::get_string(*itr, "CAS1")) != 0;
137 
138  if (swapped) {
139  std::swap(name1, name2);
140  }
141 
142  // Populate the dictionary with common terms
143  dict.add_string("name1", name1);
144  dict.add_string("name2", name2);
145  dict.add_string("BibTeX", cpjson::get_string(*itr, "BibTeX"));
146  dict.add_number("F", cpjson::get_double(*itr, "F"));
147  if (std::abs(dict.get_number("F")) > DBL_EPSILON) {
148  dict.add_string("function", cpjson::get_string(*itr, "function"));
149  }
150 
151  if (itr->HasMember("xi") && itr->HasMember("zeta")) {
152  dict.add_string("type", "Lemmon-xi-zeta");
153  // Air and HFC mixtures from Lemmon - we could also directly do the conversion
154  dict.add_number("xi", cpjson::get_double(*itr, "xi"));
155  dict.add_number("zeta", cpjson::get_double(*itr, "zeta"));
156  } else if (itr->HasMember("gammaT") && itr->HasMember("gammaV") && itr->HasMember("betaT") && itr->HasMember("betaV")) {
157  dict.add_string("type", "GERG-2008");
158  dict.add_number("gammaV", cpjson::get_double(*itr, "gammaV"));
159  dict.add_number("gammaT", cpjson::get_double(*itr, "gammaT"));
160 
161  double betaV = cpjson::get_double(*itr, "betaV");
162  double betaT = cpjson::get_double(*itr, "betaT");
163  if (swapped) {
164  dict.add_number("betaV", 1 / betaV);
165  dict.add_number("betaT", 1 / betaT);
166  } else {
167  dict.add_number("betaV", betaV);
168  dict.add_number("betaT", betaT);
169  }
170  } else {
171  std::cout << "Loading error: binary pair of " << name1 << " & " << name2
172  << "does not provide either a) xi and zeta b) gammaT, gammaV, betaT, and betaV" << std::endl;
173  continue;
174  }
175 
176  std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator it = m_binary_pair_map.find(CAS);
177  if (it == m_binary_pair_map.end()) {
178  // Add to binary pair map by creating one-element vector
179  m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
180  } else {
181  if (get_config_bool(OVERWRITE_BINARY_INTERACTION)) {
182  // Already there, see http://www.cplusplus.com/reference/map/map/insert/, so we are going to pop it and overwrite it
183  m_binary_pair_map.erase(it);
184  std::pair<std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator, bool> ret;
185  ret =
186  m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
187  assert(ret.second == true);
188  } else {
189  // Error if already in map!
190  throw ValueError(
191  format("CAS pair(%s,%s) already in binary interaction map; considering enabling configuration key OVERWRITE_BINARY_INTERACTION",
192  CAS[0].c_str(), CAS[1].c_str()));
193  }
194  }
195  }
196  }
198  void add_simple_mixing_rule(const std::string& identifier1, const std::string& identifier2, const std::string& rule) {
199  // Get the empty dictionary to be filled by the appropriate reducing parameter filling function
200  Dictionary dict;
201 
202  // Get the names/CAS of the compounds
203  std::string CAS1, CAS2, name1 = identifier1, name2 = identifier2;
204  shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> HEOS1, HEOS2;
205 
206  std::vector<std::string> id1split = strsplit(identifier1, '-');
207  if (id1split.size() == 3) { // Check if identifier is in CAS format
208  CAS1 = identifier1;
209  } else {
210  std::vector<std::string> names1(1, identifier1);
211  HEOS1.reset(new CoolProp::HelmholtzEOSMixtureBackend(names1));
212  CAS1 = HEOS1->fluid_param_string("CAS");
213  }
214 
215  std::vector<std::string> id2split = strsplit(identifier2, '-');
216  if (id2split.size() == 3) { // Check if identifier is in CAS format
217  CAS2 = identifier2;
218  } else {
219  std::vector<std::string> names2(1, identifier2);
220  HEOS2.reset(new CoolProp::HelmholtzEOSMixtureBackend(names2));
221  CAS2 = HEOS2->fluid_param_string("CAS");
222  }
223 
224  // Get the vector of CAS numbers
225  std::vector<std::string> CAS;
226  CAS.push_back(CAS1);
227  CAS.push_back(CAS2);
228 
229  // Sort the CAS number vector
230  std::sort(CAS.begin(), CAS.end());
231 
232  // Swap fluid names if the CAS numbers were swapped
233  if (CAS[0] != CAS1) {
234  std::swap(name1, name2);
235  }
236 
237  // Populate the dictionary with common terms
238  dict.add_string("name1", name1);
239  dict.add_string("name2", name2);
240  dict.add_string("BibTeX", "N/A - " + rule);
241  dict.add_number("F", 0);
242  dict.add_string("type", "GERG-2008");
243 
244  if (rule == "linear") {
245  // Terms for linear mixing
246  HEOS1.reset(new CoolProp::HelmholtzEOSMixtureBackend(std::vector<std::string>(1, name1)));
247  HEOS2.reset(new CoolProp::HelmholtzEOSMixtureBackend(std::vector<std::string>(1, name2)));
248 
249  dict.add_number("gammaT", 0.5 * (HEOS1->T_critical() + HEOS2->T_critical()) / sqrt(HEOS1->T_critical() * HEOS2->T_critical()));
250  double rhoc1 = HEOS1->rhomolar_critical(), rhoc2 = HEOS2->rhomolar_critical();
251  dict.add_number("gammaV", 4 * (1 / rhoc1 + 1 / rhoc2) / pow(1 / pow(rhoc1, 1.0 / 3.0) + 1 / pow(rhoc2, 1.0 / 3.0), 3));
252  dict.add_number("betaV", 1.0);
253  dict.add_number("betaT", 1.0);
254  } else if (rule == "Lorentz-Berthelot") {
255  // Terms for Lorentz-Berthelot quadratic mixing
256 
257  dict.add_number("gammaT", 1.0);
258  dict.add_number("gammaV", 1.0);
259  dict.add_number("betaV", 1.0);
260  dict.add_number("betaT", 1.0);
261  } else {
262  throw ValueError(format("Your simple mixing rule [%s] was not understood", rule.c_str()));
263  }
264 
265  std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator it = m_binary_pair_map.find(CAS);
266  if (it == m_binary_pair_map.end()) {
267  // Add to binary pair map by creating one-element vector
268  m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
269  } else {
270  if (get_config_bool(OVERWRITE_BINARY_INTERACTION)) {
271  // Already there, see http://www.cplusplus.com/reference/map/map/insert/, so we are going to pop it and overwrite it
272  m_binary_pair_map.erase(it);
273  std::pair<std::map<std::vector<std::string>, std::vector<Dictionary>>::iterator, bool> ret;
274  ret = m_binary_pair_map.insert(std::pair<std::vector<std::string>, std::vector<Dictionary>>(CAS, std::vector<Dictionary>(1, dict)));
275  assert(ret.second == true);
276  } else {
277  // Error if already in map!
278  throw ValueError(
279  format("CAS pair(%s,%s) already in binary interaction map; considering enabling configuration key OVERWRITE_BINARY_INTERACTION",
280  CAS[0].c_str(), CAS[1].c_str()));
281  }
282  }
283  }
284 };
285 // The modifiable parameter library
286 static MixtureBinaryPairLibrary mixturebinarypairlibrary;
287 // A fixed parameter library containing the default values
288 static MixtureBinaryPairLibrary mixturebinarypairlibrary_default;
289 
291 void apply_simple_mixing_rule(const std::string& identifier1, const std::string& identifier2, const std::string& rule) {
292  mixturebinarypairlibrary.add_simple_mixing_rule(identifier1, identifier2, rule);
293 }
294 
295 std::string get_csv_mixture_binary_pairs() {
296 
297  std::vector<std::string> out;
298  for (std::map<std::vector<std::string>, std::vector<Dictionary>>::const_iterator it = mixturebinarypairlibrary.binary_pair_map().begin();
299  it != mixturebinarypairlibrary.binary_pair_map().end(); ++it) {
300  out.push_back(strjoin(it->first, "&"));
301  }
302  return strjoin(out, ",");
303 }
304 
305 std::string get_mixture_binary_pair_data(const std::string& CAS1, const std::string& CAS2, const std::string& key) {
306  // Find pair
307  std::vector<std::string> CAS;
308  CAS.push_back(CAS1);
309  CAS.push_back(CAS2);
310 
311  if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
312  std::vector<Dictionary>& v = mixturebinarypairlibrary.binary_pair_map()[CAS];
313  try {
314  if (key == "name1") {
315  return v[0].get_string("name1");
316  } else if (key == "name2") {
317  return v[0].get_string("name2");
318  } else if (key == "BibTeX") {
319  return v[0].get_string("BibTeX");
320  } else if (key == "function") {
321  return v[0].get_string("function");
322  } else if (key == "type") {
323  return v[0].get_string("type");
324  } else if (key == "F") {
325  return format("%0.16g", v[0].get_double("F"));
326  } else if (key == "xi") {
327  return format("%0.16g", v[0].get_double("xi"));
328  } else if (key == "zeta") {
329  return format("%0.16g", v[0].get_double("zeta"));
330  } else if (key == "gammaT") {
331  return format("%0.16g", v[0].get_double("gammaT"));
332  } else if (key == "gammaV") {
333  return format("%0.16g", v[0].get_double("gammaV"));
334  } else if (key == "betaT") {
335  return format("%0.16g", v[0].get_double("betaT"));
336  } else if (key == "betaV") {
337  return format("%0.16g", v[0].get_double("betaV"));
338  } else {
339  }
340  } catch (...) {
341  }
342  throw ValueError(format("Could not match the parameter [%s] for the binary pair [%s,%s] - for now this is an error.", key.c_str(),
343  CAS1.c_str(), CAS2.c_str()));
344  } else {
345  // Sort, see if other order works properly
346  std::sort(CAS.begin(), CAS.end());
347  if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
348  throw ValueError(format("Could not match the binary pair [%s,%s] - order of CAS numbers is backwards; found the swapped CAS numbers.",
349  CAS1.c_str(), CAS2.c_str()));
350  } else {
351  throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS1.c_str(), CAS2.c_str()));
352  }
353  }
354 }
355 void set_mixture_binary_pair_data(const std::string& CAS1, const std::string& CAS2, const std::string& key, const double value) {
356 
357  // Find pair
358  std::vector<std::string> CAS;
359  CAS.push_back(CAS1);
360  CAS.push_back(CAS2);
361 
362  if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
363  std::vector<Dictionary>& v = mixturebinarypairlibrary.binary_pair_map()[CAS];
364  if (v[0].has_number(key)) {
365  v[0].add_number(key, value);
366  } else {
367  throw ValueError(format("Could not set the parameter [%s] for the binary pair [%s,%s] - for now this is an error", key.c_str(),
368  CAS1.c_str(), CAS2.c_str()));
369  }
370  } else {
371  // Sort, see if other order works properly
372  std::sort(CAS.begin(), CAS.end());
373  if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
374  throw ValueError(format("Could not match the binary pair [%s,%s] - order of CAS numbers is backwards; found the swapped CAS numbers.",
375  CAS1.c_str(), CAS2.c_str()));
376  } else {
377  throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS1.c_str(), CAS2.c_str()));
378  }
379  }
380 }
381 
382 std::string get_reducing_function_name(const std::string& CAS1, const std::string& CAS2) {
383 
384  std::vector<std::string> CAS;
385  CAS.push_back(CAS1);
386  CAS.push_back(CAS2);
387 
388  // Sort the CAS number vector - map is based on sorted CAS codes
389  std::sort(CAS.begin(), CAS.end());
390 
391  if (mixturebinarypairlibrary.binary_pair_map().find(CAS) != mixturebinarypairlibrary.binary_pair_map().end()) {
392  return mixturebinarypairlibrary.binary_pair_map()[CAS][0].get_string("function");
393  } else {
394  throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS1.c_str(), CAS2.c_str()));
395  }
396 }
397 
398 void set_interaction_parameters(const std::string& string_data) {
399  // JSON-encoded string for binary interaction parameters
400  mixturebinarypairlibrary.load_from_string(string_data);
401 }
402 
405 class MixtureDepartureFunctionsLibrary
406 {
407  private:
409  std::map<std::string, Dictionary> m_departure_function_map;
410 
411  public:
412  std::map<std::string, Dictionary>& departure_function_map() {
413  // Set the default departure functions if none have been provided yet
414  if (m_departure_function_map.size() == 0) {
415  load_defaults();
416  }
417  return m_departure_function_map;
418  };
419 
420  void load_from_string(const std::string& str) {
421  rapidjson::Document doc;
422  doc.Parse<0>(str.c_str());
423  if (doc.HasParseError()) {
424  std::cout << str << std::endl;
425  throw ValueError("Unable to parse departure function string");
426  }
427  load_from_JSON(doc);
428  }
429  void load_from_JSON(rapidjson::Document& doc) {
430  // Iterate over the departure functions in the listing
431  for (rapidjson::Value::ValueIterator itr = doc.Begin(); itr != doc.End(); ++itr) {
432  // Get the empty dictionary to be filled in
433  Dictionary dict;
434 
435  // Populate the dictionary with common terms
436  std::string Name = cpjson::get_string(*itr, "Name");
437  std::string type = cpjson::get_string(*itr, "type");
438  dict.add_string("Name", Name);
439  dict.add_string("BibTeX", cpjson::get_string(*itr, "BibTeX"));
440  dict.add_string_vector("aliases", cpjson::get_string_array(*itr, "aliases"));
441  dict.add_string("type", type);
442 
443  // Terms for the power (common to both types)
444  dict.add_double_vector("n", cpjson::get_double_array(*itr, "n"));
445  dict.add_double_vector("d", cpjson::get_double_array(*itr, "d"));
446  dict.add_double_vector("t", cpjson::get_double_array(*itr, "t"));
447 
448  // Now we need to load additional terms
449  if (!type.compare("GERG-2008")) {
450  // Number of terms that are power terms
451  dict.add_number("Npower", cpjson::get_double(*itr, "Npower"));
452  // Terms for the gaussian
453  dict.add_double_vector("eta", cpjson::get_double_array(*itr, "eta"));
454  dict.add_double_vector("epsilon", cpjson::get_double_array(*itr, "epsilon"));
455  dict.add_double_vector("beta", cpjson::get_double_array(*itr, "beta"));
456  dict.add_double_vector("gamma", cpjson::get_double_array(*itr, "gamma"));
457  } else if (type == "Gaussian+Exponential") {
458  // Number of terms that are power terms
459  dict.add_number("Npower", cpjson::get_double(*itr, "Npower"));
460  // The decay strength parameters
461  dict.add_double_vector("l", cpjson::get_double_array(*itr, "l"));
462  // Terms for the gaussian part
463  dict.add_double_vector("eta", cpjson::get_double_array(*itr, "eta"));
464  dict.add_double_vector("epsilon", cpjson::get_double_array(*itr, "epsilon"));
465  dict.add_double_vector("beta", cpjson::get_double_array(*itr, "beta"));
466  dict.add_double_vector("gamma", cpjson::get_double_array(*itr, "gamma"));
467  } else if (!type.compare("Exponential")) {
468  dict.add_double_vector("l", cpjson::get_double_array(*itr, "l"));
469  } else {
470  throw ValueError(format("It was not possible to parse departure function with type [%s]", type.c_str()));
471  }
472  // Add the normal name;
473  add_one(Name, dict);
474  std::vector<std::string> aliases = dict.get_string_vector("aliases");
475  // Add the aliases too;
476  for (std::vector<std::string>::const_iterator it = aliases.begin(); it != aliases.end(); ++it) {
477  // Add the alias;
478  add_one(*it, dict);
479  }
480  }
481  }
482  void add_one(const std::string& name, Dictionary& dict) {
483 
484  // Check if this name is already in use
485  std::map<std::string, Dictionary>::iterator it = m_departure_function_map.find(name);
486  if (it == m_departure_function_map.end()) {
487  // Not in map, add new entry to map with dictionary as value and Name as key
488  m_departure_function_map.insert(std::pair<std::string, Dictionary>(name, dict));
489  } else {
490  if (get_config_bool(OVERWRITE_DEPARTURE_FUNCTION)) {
491  // Already there, see http://www.cplusplus.com/reference/map/map/insert/
492  m_departure_function_map.erase(it);
493  std::pair<std::map<std::string, Dictionary>::iterator, bool> ret;
494  ret = m_departure_function_map.insert(std::pair<std::string, Dictionary>(name, dict));
495  assert(ret.second == true);
496  } else {
497  // Error if already in map!
498  //
499  // Collect all the current names for departure functions for a nicer error message
500  std::vector<std::string> names;
501  for (std::map<std::string, Dictionary>::const_iterator it = m_departure_function_map.begin(); it != m_departure_function_map.end();
502  ++it) {
503  names.push_back(it->first);
504  }
505  throw ValueError(format("Name of departure function [%s] is already loaded. Current departure function names are: %s", name.c_str(),
506  strjoin(names, ",").c_str()));
507  }
508  }
509  }
510  // Load the defaults that come from the JSON-encoded string compiled into library
511  // as the variable mixture_departure_functions_JSON
512  void load_defaults() {
513  load_from_string(mixture_departure_functions_JSON);
514  }
515 };
516 static MixtureDepartureFunctionsLibrary mixturedeparturefunctionslibrary;
517 
518 DepartureFunction* get_departure_function(const std::string& Name) {
519  // Get the dictionary itself
520  Dictionary& dict_dep = mixturedeparturefunctionslibrary.departure_function_map()[Name];
521 
522  if (dict_dep.is_empty()) {
523  throw ValueError(format("Departure function name [%s] seems to be invalid", Name.c_str()));
524  }
525 
526  // These terms are common
527  std::vector<double> n = dict_dep.get_double_vector("n");
528  std::vector<double> d = dict_dep.get_double_vector("d");
529  std::vector<double> t = dict_dep.get_double_vector("t");
530 
531  std::string type_dep = dict_dep.get_string("type");
532 
533  if (!type_dep.compare("GERG-2008")) {
534  // Number of power terms needed
535  int Npower = static_cast<int>(dict_dep.get_number("Npower"));
536  // Terms for the gaussian
537  std::vector<double> eta = dict_dep.get_double_vector("eta");
538  std::vector<double> epsilon = dict_dep.get_double_vector("epsilon");
539  std::vector<double> beta = dict_dep.get_double_vector("beta");
540  std::vector<double> gamma = dict_dep.get_double_vector("gamma");
541  return new GERG2008DepartureFunction(n, d, t, eta, epsilon, beta, gamma, Npower);
542  } else if (!type_dep.compare("Exponential")) {
543  // Powers of the exponents inside the exponential term
544  std::vector<double> l = dict_dep.get_double_vector("l");
545  return new ExponentialDepartureFunction(n, d, t, l);
546  } else if (!type_dep.compare("Gaussian+Exponential")) {
547  // Number of power terms needed
548  int Npower = static_cast<int>(dict_dep.get_number("Npower"));
549  // Powers of the exponents inside the exponential term
550  std::vector<double> l = dict_dep.get_double_vector("l");
551  // Terms for the gaussian
552  std::vector<double> eta = dict_dep.get_double_vector("eta");
553  std::vector<double> epsilon = dict_dep.get_double_vector("epsilon");
554  std::vector<double> beta = dict_dep.get_double_vector("beta");
555  std::vector<double> gamma = dict_dep.get_double_vector("gamma");
556  return new GaussianExponentialDepartureFunction(n, d, t, l, eta, epsilon, beta, gamma, Npower);
557  } else {
558  throw ValueError();
559  }
560 }
561 void MixtureParameters::set_mixture_parameters(HelmholtzEOSMixtureBackend& HEOS) {
562 
563  std::vector<CoolPropFluid> components = HEOS.get_components();
564 
565  std::size_t N = components.size();
566 
567  STLMatrix beta_v, gamma_v, beta_T, gamma_T;
568  beta_v.resize(N, std::vector<CoolPropDbl>(N, 0));
569  gamma_v.resize(N, std::vector<CoolPropDbl>(N, 0));
570  beta_T.resize(N, std::vector<CoolPropDbl>(N, 0));
571  gamma_T.resize(N, std::vector<CoolPropDbl>(N, 0));
572 
573  HEOS.residual_helmholtz->Excess.resize(N);
574 
575  for (std::size_t i = 0; i < N; ++i) {
576  for (std::size_t j = 0; j < N; ++j) {
577  if (i == j) {
578  continue;
579  }
580 
581  std::string CAS1 = components[i].CAS;
582  std::vector<std::string> CAS(2, "");
583  CAS[0] = components[i].CAS;
584  CAS[1] = components[j].CAS;
585  std::sort(CAS.begin(), CAS.end());
586 
587  // The variable swapped is true if a swap occurred.
588  bool swapped = (CAS1.compare(CAS[0]) != 0);
589 
590  // ***************************************************
591  // Reducing parameters for binary pair
592  // ***************************************************
593 
594  if (mixturebinarypairlibrary.binary_pair_map().find(CAS) == mixturebinarypairlibrary.binary_pair_map().end()) {
595  throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.", CAS[0].c_str(), CAS[1].c_str()));
596  }
597 
598  // Get a reference to the first matching binary pair in the dictionary
599  Dictionary& dict_red = mixturebinarypairlibrary.binary_pair_map()[CAS][0];
600 
601  // Get the name of the type being used, one of GERG-2008, Lemmon-xi-zeta, etc.
602  std::string type_red = dict_red.get_string("type");
603 
604  if (!type_red.compare("GERG-2008")) {
605  if (swapped) {
606  beta_v[i][j] = 1 / dict_red.get_number("betaV");
607  beta_T[i][j] = 1 / dict_red.get_number("betaT");
608  } else {
609  beta_v[i][j] = dict_red.get_number("betaV");
610  beta_T[i][j] = dict_red.get_number("betaT");
611  }
612  gamma_v[i][j] = dict_red.get_number("gammaV");
613  gamma_T[i][j] = dict_red.get_number("gammaT");
614  } else if (!type_red.compare("Lemmon-xi-zeta")) {
615  LemmonAirHFCReducingFunction::convert_to_GERG(components, i, j, dict_red, beta_T[i][j], beta_v[i][j], gamma_T[i][j], gamma_v[i][j]);
616  } else {
617  throw ValueError(format("type [%s] for reducing function for pair [%s, %s] is invalid", type_red.c_str(),
618  dict_red.get_string("Name1").c_str(), dict_red.get_string("Name2").c_str()));
619  }
620  /*
621  if (i == 0){
622  std::cout << format("betaT %10.9Lg gammaT %10.9Lg betaV %10.9Lg gammaV %10.9Lg %s",
623  beta_T[i][j], gamma_T[i][j], beta_v[i][j], gamma_v[i][j], get_mixture_binary_pair_data(CAS[0],CAS[1],"gammaT").c_str()) << std::endl;
624  }
625  */
626 
627  // ***************************************************
628  // Departure functions used in excess term
629  // ***************************************************
630 
631  // Set the scaling factor F for the excess term
632  HEOS.residual_helmholtz->Excess.F[i][j] = dict_red.get_number("F");
633 
634  if (std::abs(HEOS.residual_helmholtz->Excess.F[i][j]) < DBL_EPSILON) {
635  // Empty departure function that will just return 0
636  std::vector<double> n(1, 0), d(1, 1), t(1, 1), l(1, 0);
637  HEOS.residual_helmholtz->Excess.DepartureFunctionMatrix[i][j].reset(new ExponentialDepartureFunction(n, d, t, l));
638  continue;
639  }
640 
641  // Get the name of the departure function to be used for this binary pair
642  std::string Name = CoolProp::get_reducing_function_name(components[i].CAS, components[j].CAS);
643 
644  HEOS.residual_helmholtz->Excess.DepartureFunctionMatrix[i][j].reset(get_departure_function(Name));
645  }
646  }
647  // We have obtained all the parameters needed for the reducing function, now set the reducing function for the mixture
648  HEOS.Reducing = shared_ptr<ReducingFunction>(new GERG2008ReducingFunction(components, beta_v, gamma_v, beta_T, gamma_T));
649 }
650 
651 void parse_HMX_BNC(const std::string& s, std::vector<REFPROP_binary_element>& BIP, std::vector<REFPROP_departure_function>& functions) {
652  // Capture the betas, gammas, Fij, models
653  bool block_started = false;
654  std::size_t i_started = 0, i_ended = 0, i = 0;
655  std::vector<std::string> lines = strsplit(s, '\n');
656  for (std::vector<std::string>::iterator it = lines.begin(); it != lines.end(); ++it) {
657  if (strstartswith(strstrip(*it), "#BNC")) {
658  block_started = true;
659  i_started = i + 1;
660  }
661  if (block_started && strstrip(*it).empty()) {
662  i_ended = i - 1;
663  break;
664  }
665  i++;
666  }
667  // Find the first line with a !
668  for (i = i_started; i < i_ended; ++i) {
669  if (strstrip(lines[i]) == "!") {
670  i_started = i;
671  break;
672  }
673  }
674  // Find all the lines that '!' - these are delimiters
675  std::vector<std::pair<std::size_t, std::size_t>> bounds; // The left and right indices (inclusive) that form a binary pair
676  std::size_t last_excalamation = i_started;
677  for (i = i_started; i <= i_ended; ++i) {
678  if (strstrip(lines[i]) == "!") {
679  bounds.push_back(std::make_pair<std::size_t, std::size_t>(last_excalamation + 1, i - 1));
680  last_excalamation = i;
681  }
682  }
683  // Parse each chunk
684  std::vector<REFPROP_binary_element> chunks;
685  for (std::vector<std::pair<std::size_t, std::size_t>>::iterator it = bounds.begin(); it != bounds.end(); ++it) {
686  REFPROP_binary_element bnc;
687  for (std::size_t i = it->first; i <= it->second; ++i) {
688  // Store comments
689  if (strstartswith(lines[i], "?")) {
690  bnc.comments.push_back(lines[i]);
691  continue;
692  }
693  // Parse the line with the thermo BIP
694  if (lines[i].find("/") > 0) {
695  // Split at ' '
696  std::vector<std::string> bits = strsplit(strstrip(lines[i]), ' ');
697  // Remove empty elements
698  for (std::size_t j = bits.size() - 1; j > 0; --j) {
699  if (bits[j].empty()) {
700  bits.erase(bits.begin() + j);
701  }
702  }
703  // Get the line that contains the thermo BIP
704  if (bits[0].find("/") > 0 && bits[1].size() == 3) {
705  std::vector<std::string> theCAS = strsplit(bits[0], '/');
706  bnc.CAS1 = theCAS[0];
707  bnc.CAS2 = theCAS[1];
708  bnc.model = bits[1];
709  bnc.betaT = string2double(bits[2]);
710  bnc.gammaT = string2double(bits[3]);
711  bnc.betaV = string2double(bits[4]);
712  bnc.gammaV = string2double(bits[5]);
713  bnc.Fij = string2double(bits[6]);
714  break;
715  } else if (strstrip(bits[0]) == "CAS#") {
716  break;
717  } else {
718  throw CoolProp::ValueError(format("Unable to parse binary interaction line: %s", lines[i]));
719  }
720  }
721  }
722  if (!bnc.CAS1.empty()) {
723  BIP.push_back(bnc);
724  }
725  }
726 
727  // ****************************************
728  // Parse the departure functions
729  // ****************************************
730  for (std::size_t i = i_ended + 1; i < lines.size(); ++i) {
731  std::size_t j_end;
732  // Find the end of this block
733  for (j_end = i + 1; j_end < lines.size(); ++j_end) {
734  if (strstrip(lines[j_end]).empty()) {
735  j_end -= 1;
736  break;
737  }
738  }
739 
740  if (strstartswith(lines[i], "#MXM")) {
741  REFPROP_departure_function dep;
742  dep.Npower = -1;
743  dep.model = std::string(lines[i + 1].begin(), lines[i + 1].begin() + 3);
744  dep.comments.push_back(lines[i + 1]);
745  for (std::size_t j = i + 2; j <= j_end; ++j) {
746  if (strstartswith(strstrip(lines[j]), "?")) {
747  dep.comments.push_back(lines[j]);
748  continue;
749  }
750  if (strstartswith(strstrip(lines[j]), "!")) {
751  j += 2; // Skip the BIP here, not used
752  continue;
753  }
754  std::vector<std::string> bits = strsplit(lines[j], ' ');
755  // Remove empty elements
756  for (std::size_t k = bits.size() - 1; k > 0; --k) {
757  if (bits[k].empty()) {
758  bits.erase(bits.begin() + k);
759  }
760  }
761 
762  if (dep.Npower < 0) { // Not extracted yet, let's do it now
763  // Extract the number of terms
764  dep.Npower = static_cast<short>(strtol(bits[0].c_str(), NULL, 10));
765  dep.Nterms_power = static_cast<short>(strtol(bits[1].c_str(), NULL, 10));
766  dep.Nspecial = static_cast<short>(strtol(bits[3].c_str(), NULL, 10));
767  dep.Nterms_special = static_cast<short>(strtol(bits[4].c_str(), NULL, 10));
768  } else {
769  dep.a.push_back(string2double(bits[0]));
770  dep.t.push_back(string2double(bits[1]));
771  dep.d.push_back(string2double(bits[2]));
772  // Extracting "polynomial" terms
773  if (dep.Nterms_power == 4) {
774  dep.e.push_back(string2double(bits[3]));
775  }
776  if (dep.Nspecial > 0) {
777  if (dep.a.size() - 1 < dep.Npower) {
778  dep.eta.push_back(0);
779  dep.epsilon.push_back(0);
780  dep.beta.push_back(0);
781  dep.gamma.push_back(0);
782  } else {
783  // Extracting "special" terms
784  dep.eta.push_back(string2double(bits[3]));
785  dep.epsilon.push_back(string2double(bits[4]));
786  dep.beta.push_back(string2double(bits[5]));
787  dep.gamma.push_back(string2double(bits[6]));
788  }
789  }
790  }
791  }
792  functions.push_back(dep);
793  }
794  }
795 }
796 
797 void set_departure_functions(const std::string& string_data) {
798  if (string_data.find("#MXM") != std::string::npos) {
799  // REFPROP HMX.BNC file was provided
800  std::vector<REFPROP_binary_element> BIP;
801  std::vector<REFPROP_departure_function> functions;
802  parse_HMX_BNC(string_data, BIP, functions);
803 
804  {
805  rapidjson::Document doc;
806  doc.SetArray();
807  for (std::vector<REFPROP_binary_element>::const_iterator it = BIP.begin(); it < BIP.end(); ++it) {
808  rapidjson::Value el;
809  el.SetObject();
810  el.AddMember("CAS1", rapidjson::Value(it->CAS1.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
811  el.AddMember("CAS2", rapidjson::Value(it->CAS2.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
812  el.AddMember("Name1", "??", doc.GetAllocator());
813  el.AddMember("Name2", "??", doc.GetAllocator());
814  el.AddMember("betaT", it->betaT, doc.GetAllocator());
815  el.AddMember("gammaT", it->gammaT, doc.GetAllocator());
816  el.AddMember("betaV", it->betaV, doc.GetAllocator());
817  el.AddMember("gammaV", it->gammaV, doc.GetAllocator());
818  el.AddMember("F", it->Fij, doc.GetAllocator());
819  el.AddMember("function", rapidjson::Value(it->model.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
820  std::string tex_string = "(from HMX.BNC format)::" + strjoin(it->comments, "\n");
821  el.AddMember("BibTeX", rapidjson::Value(tex_string.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
822  doc.PushBack(el, doc.GetAllocator());
823  }
824  mixturebinarypairlibrary.load_from_JSON(doc);
825  }
826  {
827  rapidjson::Document doc;
828  doc.SetArray();
829  for (std::vector<REFPROP_departure_function>::const_iterator it = functions.begin(); it < functions.end(); ++it) {
830  rapidjson::Value el;
831  el.SetObject();
832  el.AddMember("Name", rapidjson::Value(it->model.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
833  std::vector<std::string> aliases;
834  cpjson::set_string_array("aliases", aliases, el, doc);
835  cpjson::set_double_array("n", it->a, el, doc);
836  cpjson::set_double_array("d", it->d, el, doc);
837  cpjson::set_double_array("t", it->t, el, doc);
838  if (it->Nterms_special > 0 || it->Nterms_power == 3) {
839  el.AddMember("type", "GERG-2008", doc.GetAllocator());
840  el.AddMember("Npower", it->Npower, doc.GetAllocator());
841  if (it->Nterms_power == 3 && it->Nspecial == 0) {
842  std::vector<double> zeros(it->a.size(), 0);
843  cpjson::set_double_array("eta", zeros, el, doc);
844  cpjson::set_double_array("epsilon", zeros, el, doc);
845  cpjson::set_double_array("beta", zeros, el, doc);
846  cpjson::set_double_array("gamma", zeros, el, doc);
847  } else {
848  cpjson::set_double_array("eta", it->eta, el, doc);
849  cpjson::set_double_array("epsilon", it->epsilon, el, doc);
850  cpjson::set_double_array("beta", it->beta, el, doc);
851  cpjson::set_double_array("gamma", it->gamma, el, doc);
852  }
853  } else {
854  el.AddMember("type", "Exponential", doc.GetAllocator());
855  cpjson::set_double_array("l", it->e, el, doc);
856  }
857 
858  std::string tex_string = "(from HMX.BNC format)::" + strjoin(it->comments, "\n");
859  el.AddMember("BibTeX", rapidjson::Value(tex_string.c_str(), doc.GetAllocator()).Move(), doc.GetAllocator());
860  doc.PushBack(el, doc.GetAllocator());
861  }
862  mixturedeparturefunctionslibrary.load_from_JSON(doc);
863  }
864  } else {
865  // JSON-encoded string for departure functions
866  mixturedeparturefunctionslibrary.load_from_string(string_data);
867  }
868 }
869 
870 } /* namespace CoolProp */
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1