CoolProp 6.8.1dev
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
7namespace CoolProp {
8
13class 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};
53static PredefinedMixturesLibrary predefined_mixtures_library;
54
55std::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
64bool 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
74void 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
83class 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
286static MixtureBinaryPairLibrary mixturebinarypairlibrary;
287// A fixed parameter library containing the default values
288static MixtureBinaryPairLibrary mixturebinarypairlibrary_default;
289
291void 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
295std::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
305std::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}
355void 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
382std::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
398void 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
405class 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};
516static MixtureDepartureFunctionsLibrary mixturedeparturefunctionslibrary;
517
518DepartureFunction* 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}
561void 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
651void 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
797void 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 Mon Jun 2 2025 05:21:52 for CoolProp by doxygen 1.9.4