CoolProp 6.8.0
An open-source fluid property and humid air property database
AbstractState.cpp
Go to the documentation of this file.
1/*
2 * AbstractState.cpp
3 *
4 * Created on: 21 Dec 2013
5 * Author: jowr
6 */
7
8#ifndef _CRT_SECURE_NO_WARNINGS
9#define _CRT_SECURE_NO_WARNINGS
10#endif
11
12#include <stdlib.h>
13#include "math.h"
14#include "AbstractState.h"
15#include "DataStructures.h"
16#include "Backends/IF97/IF97Backend.h"
17#include "Backends/Cubics/CubicBackend.h"
18#include "Backends/Cubics/VTPRBackend.h"
19#include "Backends/Incompressible/IncompressibleBackend.h"
20#include "Backends/PCSAFT/PCSAFTBackend.h"
21
22#if !defined(NO_TABULAR_BACKENDS)
23#include "Backends/Tabular/TTSEBackend.h"
24#include "Backends/Tabular/BicubicBackend.h"
25#endif
26
27namespace CoolProp {
28
31
32class BackendLibrary
33{
34 private:
35 std::map<backend_families, shared_ptr<AbstractStateGenerator>> backends;
36
37 public:
38 void add_backend(const backend_families& bg, const shared_ptr<AbstractStateGenerator>& asg) {
39 backends[bg] = asg;
40 };
41 void get_generator_iterators(const backend_families& bg,
42 std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& generator,
43 std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator& end) {
44 generator = backends.find(bg);
45 end = backends.end();
46 };
47 std::size_t size() {
48 return backends.size();
49 };
50};
51inline BackendLibrary& get_backend_library() {
52 static BackendLibrary the_library;
53 return the_library;
54}
55
56void register_backend(const backend_families& bf, shared_ptr<AbstractStateGenerator> gen) {
57 get_backend_library().add_backend(bf, gen);
58};
59
60class IF97BackendGenerator : public AbstractStateGenerator
61{
62 public:
63 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
64 if (fluid_names.size() == 1) { // Check that fluid_names[0] has only one component
65 std::string str = fluid_names[0]; // Check that the fluid name is an alias for "Water"
66 if ((upper(str) == "WATER") || (upper(str) == "H2O")) {
67 return new IF97Backend();
68 } else {
69 throw ValueError(format("The IF97 backend returns Water props only; fluid name [%s] not allowed", fluid_names[0].c_str()));
70 }
71 } else {
72 throw ValueError(format("The IF97 backend does not support mixtures, only Water"));
73 };
74 };
75};
76// This static initialization will cause the generator to register
77static GeneratorInitializer<IF97BackendGenerator> if97_gen(IF97_BACKEND_FAMILY);
78class SRKGenerator : public AbstractStateGenerator
79{
80 public:
81 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
82 return new SRKBackend(fluid_names, get_config_double(R_U_CODATA));
83 };
84};
85static GeneratorInitializer<SRKGenerator> srk_gen(CoolProp::SRK_BACKEND_FAMILY);
86class PRGenerator : public AbstractStateGenerator
87{
88 public:
89 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
90 return new PengRobinsonBackend(fluid_names, get_config_double(R_U_CODATA));
91 };
92};
93static GeneratorInitializer<PRGenerator> pr_gen(CoolProp::PR_BACKEND_FAMILY);
94class IncompressibleBackendGenerator : public AbstractStateGenerator
95{
96 public:
97 AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
98 if (fluid_names.size() != 1) {
99 throw ValueError(format("For INCOMP backend, name vector must be one element long"));
100 }
101 return new IncompressibleBackend(fluid_names[0]);
102 };
103};
104// This static initialization will cause the generator to register
105static GeneratorInitializer<IncompressibleBackendGenerator> incomp_gen(INCOMP_BACKEND_FAMILY);
106class VTPRGenerator : public CoolProp::AbstractStateGenerator
107{
108 public:
109 CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
110 return new CoolProp::VTPRBackend(fluid_names, CoolProp::get_config_double(R_U_CODATA));
111 };
112};
113// This static initialization will cause the generator to register
114static CoolProp::GeneratorInitializer<VTPRGenerator> vtpr_gen(CoolProp::VTPR_BACKEND_FAMILY);
115
116class PCSAFTGenerator : public CoolProp::AbstractStateGenerator
117{
118 public:
119 CoolProp::AbstractState* get_AbstractState(const std::vector<std::string>& fluid_names) {
120 return new CoolProp::PCSAFTBackend(fluid_names);
121 };
122};
123// This static initialization will cause the generator to register
124static CoolProp::GeneratorInitializer<PCSAFTGenerator> pcsaft_gen(CoolProp::PCSAFT_BACKEND_FAMILY);
125
126AbstractState* AbstractState::factory(const std::string& backend, const std::vector<std::string>& fluid_names) {
127 if (get_debug_level() > 0) {
128 std::cout << "AbstractState::factory(" << backend << "," << stringvec_to_string(fluid_names) << ")" << std::endl;
129 }
130
131 backend_families f1;
132 std::string f2;
133 extract_backend_families_string(backend, f1, f2);
134
135 std::map<backend_families, shared_ptr<AbstractStateGenerator>>::const_iterator gen, end;
136 get_backend_library().get_generator_iterators(f1, gen, end);
137
138 if (get_debug_level() > 0) {
139 std::cout << "AbstractState::factory backend_library size: " << get_backend_library().size() << std::endl;
140 }
141
142 if (gen != end) {
143 // One of the registered backends was able to match the given backend family
144 return gen->second->get_AbstractState(fluid_names);
145 }
146#if !defined(NO_TABULAR_BACKENDS)
147 else if (f1 == TTSE_BACKEND_FAMILY) {
148 // Will throw if there is a problem with this backend
149 shared_ptr<AbstractState> AS(factory(f2, fluid_names));
150 return new TTSEBackend(AS);
151 } else if (f1 == BICUBIC_BACKEND_FAMILY) {
152 // Will throw if there is a problem with this backend
153 shared_ptr<AbstractState> AS(factory(f2, fluid_names));
154 return new BicubicBackend(AS);
155 }
156#endif
157 else if (!backend.compare("?") || backend.empty()) {
158 std::size_t idel = fluid_names[0].find("::");
159 // Backend has not been specified, and we have to figure out what the backend is by parsing the string
160 if (idel == std::string::npos) // No '::' found, no backend specified, try HEOS, otherwise a failure
161 {
162 // Figure out what backend to use
163 return factory("HEOS", fluid_names);
164 } else {
165 // Split string at the '::' into two std::string, call again
166 return factory(std::string(fluid_names[0].begin(), fluid_names[0].begin() + idel),
167 std::string(fluid_names[0].begin() + (idel + 2), fluid_names[0].end()));
168 }
169 } else {
170 throw ValueError(format("Invalid backend name [%s] to factory function", backend.c_str()));
171 }
172}
173std::vector<std::string> AbstractState::fluid_names(void) {
174 return calc_fluid_names();
175}
176bool AbstractState::clear_comp_change() {
177 // Reset all instances of CachedElement and overwrite
178 // the internal double values with -_HUGE
179 this->_R = _HUGE;
180 this->_gas_constant.clear();
181 this->_molar_mass.clear();
182 this->_critical.fill(_HUGE);
183 this->_reducing.fill(_HUGE);
184 return true;
185}
186bool AbstractState::clear() {
187 // Reset all instances of CachedElement and overwrite
188 // the internal double values with -_HUGE
189 this->_R = _HUGE;
190 this->_gas_constant.clear();
191 this->_molar_mass.clear();
192
194 this->_rhoLanc.clear();
195 this->_rhoVanc.clear();
196 this->_pVanc.clear();
197 this->_pLanc.clear();
198 this->_TVanc.clear();
199 this->_TLanc.clear();
200
201 this->_critical.fill(_HUGE);
202 this->_reducing.fill(_HUGE);
203
205 this->_rhomolar = -_HUGE;
206 this->_T = -_HUGE;
207 this->_p = -_HUGE;
208 this->_Q = -_HUGE;
209 this->_tau.clear();
210 this->_delta.clear();
211
212 this->_umolar.clear();
213 this->_cpmolar.clear();
214 this->_cp0molar.clear();
215 this->_cvmolar.clear();
216 this->_speed_sound.clear();
217 this->_hmolar.clear();
218 this->_smolar.clear();
219 this->_gibbsmolar.clear();
220 this->_helmholtzmolar.clear();
221 this->_logp.clear();
222 this->_logrhomolar.clear();
223
224 this->_hmolar_excess.clear();
225 this->_smolar_excess.clear();
226 this->_gibbsmolar_excess.clear();
227 this->_volumemolar_excess.clear();
228 this->_umolar_excess.clear();
229 this->_helmholtzmolar_excess.clear();
230
231 this->_hmolar_residual.clear();
232 this->_smolar_residual.clear();
233 this->_gibbsmolar_residual.clear();
234
236 this->_rho_spline.clear();
237 this->_drho_spline_dh__constp.clear();
238 this->_drho_spline_dp__consth.clear();
239
241 this->_alpha0.clear();
242 this->_dalpha0_dTau.clear();
243 this->_dalpha0_dDelta.clear();
244 this->_d2alpha0_dTau2.clear();
245 this->_d2alpha0_dDelta_dTau.clear();
246 this->_d2alpha0_dDelta2.clear();
247 this->_d3alpha0_dTau3.clear();
248 this->_d3alpha0_dDelta_dTau2.clear();
249 this->_d3alpha0_dDelta2_dTau.clear();
250 this->_d3alpha0_dDelta3.clear();
251 this->_alphar.clear();
252 this->_dalphar_dTau.clear();
253 this->_dalphar_dDelta.clear();
254 this->_d2alphar_dTau2.clear();
255 this->_d2alphar_dDelta_dTau.clear();
256 this->_d2alphar_dDelta2.clear();
257 this->_d3alphar_dTau3.clear();
258 this->_d3alphar_dDelta_dTau2.clear();
259 this->_d3alphar_dDelta2_dTau.clear();
260 this->_d3alphar_dDelta3.clear();
261
262 this->_dalphar_dDelta_lim.clear();
263 this->_d2alphar_dDelta2_lim.clear();
264 this->_d2alphar_dDelta_dTau_lim.clear();
265 this->_d3alphar_dDelta2_dTau_lim.clear();
266
268 this->_rhoLmolar.clear();
269 this->_rhoVmolar.clear();
270
272 this->_viscosity.clear();
273 this->_conductivity.clear();
274 this->_surface_tension.clear();
275
276 return true;
277}
278void AbstractState::mass_to_molar_inputs(CoolProp::input_pairs& input_pair, CoolPropDbl& value1, CoolPropDbl& value2) {
279 // Check if a mass based input, convert it to molar units
280
281 switch (input_pair) {
282 case DmassT_INPUTS:
283 //case HmassT_INPUTS: ///< Enthalpy in J/kg, Temperature in K (NOT CURRENTLY IMPLEMENTED)
284 case SmassT_INPUTS:
285 //case TUmass_INPUTS: ///< Temperature in K, Internal energy in J/kg (NOT CURRENTLY IMPLEMENTED)
286 case DmassP_INPUTS:
287 case DmassQ_INPUTS:
288 case HmassP_INPUTS:
289 case PSmass_INPUTS:
290 case PUmass_INPUTS:
291 case HmassSmass_INPUTS:
292 case SmassUmass_INPUTS:
293 case DmassHmass_INPUTS:
294 case DmassSmass_INPUTS:
295 case DmassUmass_INPUTS:
296 {
297 // Set the cache value for the molar mass if it hasn't been set yet
298 molar_mass();
299
300 // Molar mass (just for compactness of the following switch)
301 CoolPropDbl mm = static_cast<CoolPropDbl>(_molar_mass);
302
303 switch (input_pair) {
304 case DmassT_INPUTS:
305 input_pair = DmolarT_INPUTS;
306 value1 /= mm;
307 break;
308 //case HmassT_INPUTS: input_pair = HmolarT_INPUTS; value1 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
309 case SmassT_INPUTS:
310 input_pair = SmolarT_INPUTS;
311 value1 *= mm;
312 break;
313 //case TUmass_INPUTS: input_pair = TUmolar_INPUTS; value2 *= mm; break; (NOT CURRENTLY IMPLEMENTED)
314 case DmassP_INPUTS:
315 input_pair = DmolarP_INPUTS;
316 value1 /= mm;
317 break;
318 case DmassQ_INPUTS:
319 input_pair = DmolarQ_INPUTS;
320 value1 /= mm;
321 break;
322 case HmassP_INPUTS:
323 input_pair = HmolarP_INPUTS;
324 value1 *= mm;
325 break;
326 case PSmass_INPUTS:
327 input_pair = PSmolar_INPUTS;
328 value2 *= mm;
329 break;
330 case PUmass_INPUTS:
331 input_pair = PUmolar_INPUTS;
332 value2 *= mm;
333 break;
334 case HmassSmass_INPUTS:
335 input_pair = HmolarSmolar_INPUTS;
336 value1 *= mm;
337 value2 *= mm;
338 break;
339 case SmassUmass_INPUTS:
340 input_pair = SmolarUmolar_INPUTS;
341 value1 *= mm;
342 value2 *= mm;
343 break;
344 case DmassHmass_INPUTS:
345 input_pair = DmolarHmolar_INPUTS;
346 value1 /= mm;
347 value2 *= mm;
348 break;
349 case DmassSmass_INPUTS:
350 input_pair = DmolarSmolar_INPUTS;
351 value1 /= mm;
352 value2 *= mm;
353 break;
354 case DmassUmass_INPUTS:
355 input_pair = DmolarUmolar_INPUTS;
356 value1 /= mm;
357 value2 *= mm;
358 break;
359 default:
360 break;
361 }
362 break;
363 }
364 default:
365 return;
366 }
367}
368double AbstractState::trivial_keyed_output(parameters key) {
369 if (get_debug_level() >= 50)
370 std::cout << format("AbstractState: trivial_keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
371 switch (key) {
372 case imolar_mass:
373 return molar_mass();
374 case iacentric_factor:
375 return acentric_factor();
376 case igas_constant:
377 return gas_constant();
378 case iT_min:
379 return Tmin();
380 case iT_triple:
381 return Ttriple();
382 case iT_max:
383 return Tmax();
384 case iP_max:
385 return pmax();
386 case iP_min:
387 case iP_triple:
388 return this->p_triple();
389 case iT_reducing:
390 return calc_T_reducing();
391 case irhomolar_reducing:
392 return calc_rhomolar_reducing();
393 case iP_reducing:
394 return calc_p_reducing();
395 case iP_critical:
396 return this->p_critical();
397 case iT_critical:
398 return this->T_critical();
399 case irhomolar_critical:
400 return this->rhomolar_critical();
401 case irhomass_critical:
402 return this->rhomass_critical();
403 case iODP:
404 return this->calc_ODP();
405 case iGWP100:
406 return this->calc_GWP100();
407 case iGWP20:
408 return this->calc_GWP20();
409 case iGWP500:
410 return this->calc_GWP500();
411 case ifraction_min:
412 return this->calc_fraction_min();
413 case ifraction_max:
414 return this->calc_fraction_max();
415 case iT_freeze:
416 return this->calc_T_freeze();
417 case iFH:
418 return this->calc_flame_hazard();
419 case iHH:
420 return this->calc_health_hazard();
421 case iPH:
422 return this->calc_physical_hazard();
423 case idipole_moment:
424 return this->calc_dipole_moment();
425 default:
426 throw ValueError(
427 format("This input [%d: \"%s\"] is not valid for trivial_keyed_output", key, get_parameter_information(key, "short").c_str()));
428 }
429}
430double AbstractState::keyed_output(parameters key) {
431 if (get_debug_level() >= 50)
432 std::cout << format("AbstractState: keyed_output called for %s ", get_parameter_information(key, "short").c_str()) << std::endl;
433 // Handle trivial inputs
434 if (is_trivial_parameter(key)) {
435 return trivial_keyed_output(key);
436 }
437 switch (key) {
438 case iQ:
439 return Q();
440 case iT:
441 return T();
442 case iP:
443 return p();
444 case iDmolar:
445 return rhomolar();
446 case iDmass:
447 return rhomass();
448 case iHmolar:
449 return hmolar();
450 case iHmolar_residual:
451 return hmolar_residual();
452 case iHmass:
453 return hmass();
454 case iSmolar:
455 return smolar();
456 case iSmolar_residual:
457 return smolar_residual();
458 case iSmass:
459 return smass();
460 case iUmolar:
461 return umolar();
462 case iUmass:
463 return umass();
464 case iGmolar:
465 return gibbsmolar();
466 case iGmolar_residual:
467 return gibbsmolar_residual();
468 case iGmass:
469 return gibbsmass();
470 case iHelmholtzmolar:
471 return helmholtzmolar();
472 case iHelmholtzmass:
473 return helmholtzmass();
474 case iCvmolar:
475 return cvmolar();
476 case iCvmass:
477 return cvmass();
478 case iCpmolar:
479 return cpmolar();
480 case iCp0molar:
481 return cp0molar();
482 case iCpmass:
483 return cpmass();
484 case iCp0mass:
485 return cp0mass();
486 case imolar_mass:
487 return molar_mass();
488 case iT_reducing:
489 return get_reducing_state().T;
490 case irhomolar_reducing:
491 return get_reducing_state().rhomolar;
492 case ispeed_sound:
493 return speed_sound();
494 case ialphar:
495 return alphar();
496 case ialpha0:
497 return alpha0();
498 case idalpha0_ddelta_consttau:
499 return dalpha0_dDelta();
500 case id2alpha0_ddelta2_consttau:
501 return d2alpha0_dDelta2();
502 case id3alpha0_ddelta3_consttau:
503 return d3alpha0_dDelta3();
504 case idalpha0_dtau_constdelta:
505 return dalpha0_dTau();
506 case idalphar_ddelta_consttau:
507 return dalphar_dDelta();
508 case idalphar_dtau_constdelta:
509 return dalphar_dTau();
510 case iBvirial:
511 return Bvirial();
512 case idBvirial_dT:
513 return dBvirial_dT();
514 case iCvirial:
515 return Cvirial();
516 case idCvirial_dT:
517 return dCvirial_dT();
518 case iisothermal_compressibility:
519 return isothermal_compressibility();
520 case iisobaric_expansion_coefficient:
521 return isobaric_expansion_coefficient();
522 case iisentropic_expansion_coefficient:
523 return isentropic_expansion_coefficient();
524 case iviscosity:
525 return viscosity();
526 case iconductivity:
527 return conductivity();
528 case iPrandtl:
529 return Prandtl();
530 case isurface_tension:
531 return surface_tension();
532 case iPhase:
533 return phase();
534 case iZ:
535 return compressibility_factor();
536 case iPIP:
537 return PIP();
538 case ifundamental_derivative_of_gas_dynamics:
539 return fundamental_derivative_of_gas_dynamics();
540 default:
541 throw ValueError(format("This input [%d: \"%s\"] is not valid for keyed_output", key, get_parameter_information(key, "short").c_str()));
542 }
543}
544
545double AbstractState::tau(void) {
546 if (!_tau) _tau = calc_reciprocal_reduced_temperature();
547 return _tau;
548}
549double AbstractState::delta(void) {
550 if (!_delta) _delta = calc_reduced_density();
551 return _delta;
552}
553double AbstractState::Tmin(void) {
554 return calc_Tmin();
555}
556double AbstractState::Tmax(void) {
557 return calc_Tmax();
558}
559double AbstractState::Ttriple(void) {
560 return calc_Ttriple();
561}
562double AbstractState::pmax(void) {
563 return calc_pmax();
564}
565double AbstractState::T_critical(void) {
566 return calc_T_critical();
567}
568double AbstractState::T_reducing(void) {
569 if (!ValidNumber(_reducing.T)) {
570 calc_reducing_state();
571 }
572 return _reducing.T;
573}
574double AbstractState::p_critical(void) {
575 return calc_p_critical();
576}
577double AbstractState::p_triple(void) {
578 return calc_p_triple();
579}
580double AbstractState::rhomolar_critical(void) {
581 return calc_rhomolar_critical();
582}
583double AbstractState::rhomass_critical(void) {
584 return calc_rhomolar_critical() * molar_mass();
585}
586double AbstractState::rhomolar_reducing(void) {
587 if (!ValidNumber(_reducing.rhomolar)) {
588 calc_reducing_state();
589 }
590 return _reducing.rhomolar;
591}
592double AbstractState::rhomass_reducing(void) {
593 return rhomolar_reducing() * molar_mass();
594}
595double AbstractState::hmolar(void) {
596 if (!_hmolar) _hmolar = calc_hmolar();
597 return _hmolar;
598}
599double AbstractState::hmolar_residual(void) {
600 if (!_hmolar_residual) _hmolar_residual = calc_hmolar_residual();
601 return _hmolar_residual;
602}
603double AbstractState::hmolar_excess(void) {
604 if (!_hmolar_excess) calc_excess_properties();
605 return _hmolar_excess;
606}
607double AbstractState::smolar(void) {
608 if (!_smolar) _smolar = calc_smolar();
609 return _smolar;
610}
611double AbstractState::smolar_residual(void) {
612 if (!_smolar_residual) _smolar_residual = calc_smolar_residual();
613 return _smolar_residual;
614}
615double AbstractState::neff(void) {
616 double tau = calc_T_reducing()/_T;
617 double delta = _rhomolar/calc_rhomolar_reducing();
618 double Ar01 = delta*dalphar_dDelta();
619 double Ar11 = tau*delta*d2alphar_dDelta_dTau();
620 double Ar20 = tau*tau*d2alphar_dTau2();
621 return -3.0*(Ar01-Ar11)/Ar20;
622}
623double AbstractState::smolar_excess(void) {
624 if (!_smolar_excess) calc_excess_properties();
625 return _smolar_excess;
626}
627double AbstractState::umolar(void) {
628 if (!_umolar) _umolar = calc_umolar();
629 return _umolar;
630}
631double AbstractState::umolar_excess(void) {
632 if (!_umolar_excess) calc_excess_properties();
633 return _umolar_excess;
634}
635double AbstractState::gibbsmolar(void) {
636 if (!_gibbsmolar) _gibbsmolar = calc_gibbsmolar();
637 return _gibbsmolar;
638}
639double AbstractState::gibbsmolar_residual(void) {
640 if (!_gibbsmolar_residual) _gibbsmolar_residual = calc_gibbsmolar_residual();
641 return _gibbsmolar_residual;
642}
643double AbstractState::gibbsmolar_excess(void) {
644 if (!_gibbsmolar_excess) calc_excess_properties();
645 return _gibbsmolar_excess;
646}
647double AbstractState::helmholtzmolar(void) {
648 if (!_helmholtzmolar) _helmholtzmolar = calc_helmholtzmolar();
649 return _helmholtzmolar;
650}
651double AbstractState::helmholtzmolar_excess(void) {
652 if (!_helmholtzmolar_excess) calc_excess_properties();
653 return _helmholtzmolar_excess;
654}
655double AbstractState::volumemolar_excess(void) {
656 if (!_volumemolar_excess) calc_excess_properties();
657 return _volumemolar_excess;
658}
659double AbstractState::cpmolar(void) {
660 if (!_cpmolar) _cpmolar = calc_cpmolar();
661 return _cpmolar;
662}
663double AbstractState::cp0molar(void) {
664 return calc_cpmolar_idealgas();
665}
666double AbstractState::cvmolar(void) {
667 if (!_cvmolar) _cvmolar = calc_cvmolar();
668 return _cvmolar;
669}
670double AbstractState::speed_sound(void) {
671 if (!_speed_sound) _speed_sound = calc_speed_sound();
672 return _speed_sound;
673}
674double AbstractState::viscosity(void) {
675 if (!_viscosity) _viscosity = calc_viscosity();
676 return _viscosity;
677}
678double AbstractState::conductivity(void) {
679 if (!_conductivity) _conductivity = calc_conductivity();
680 return _conductivity;
681}
682double AbstractState::melting_line(int param, int given, double value) {
683 return calc_melting_line(param, given, value);
684}
685double AbstractState::acentric_factor() {
686 return calc_acentric_factor();
687}
688double AbstractState::saturation_ancillary(parameters param, int Q, parameters given, double value) {
689 return calc_saturation_ancillary(param, Q, given, value);
690}
691double AbstractState::surface_tension(void) {
692 if (!_surface_tension) _surface_tension = calc_surface_tension();
693 return _surface_tension;
694}
695double AbstractState::molar_mass(void) {
696 if (!_molar_mass) _molar_mass = calc_molar_mass();
697 return _molar_mass;
698}
699double AbstractState::gas_constant(void) {
700 if (!_gas_constant) _gas_constant = calc_gas_constant();
701 return _gas_constant;
702}
703double AbstractState::fugacity_coefficient(std::size_t i) {
704 // TODO: Cache the fug. coeff for each component
705 return calc_fugacity_coefficient(i);
706}
707std::vector<double> AbstractState::fugacity_coefficients() {
708 // TODO: Cache the fug. coeff for each component
709 return calc_fugacity_coefficients();
710}
711double AbstractState::fugacity(std::size_t i) {
712 // TODO: Cache the fug. coeff for each component
713 return calc_fugacity(i);
714}
715double AbstractState::chemical_potential(std::size_t i) {
716 // TODO: Cache the chemical potential for each component
717 return calc_chemical_potential(i);
718}
719void AbstractState::build_phase_envelope(const std::string& type) {
720 calc_phase_envelope(type);
721}
722double AbstractState::isothermal_compressibility(void) {
723 return 1.0 / _rhomolar * first_partial_deriv(iDmolar, iP, iT);
724}
725double AbstractState::isobaric_expansion_coefficient(void) {
726 return -1.0 / _rhomolar * first_partial_deriv(iDmolar, iT, iP);
727}
728double AbstractState::isentropic_expansion_coefficient(void) {
729 return _rhomolar / _p * first_partial_deriv(iP, iDmolar, iSmolar);
730}
731double AbstractState::Bvirial(void) {
732 return calc_Bvirial();
733}
734double AbstractState::Cvirial(void) {
735 return calc_Cvirial();
736}
737double AbstractState::dBvirial_dT(void) {
738 return calc_dBvirial_dT();
739}
740double AbstractState::dCvirial_dT(void) {
741 return calc_dCvirial_dT();
742}
743double AbstractState::compressibility_factor(void) {
744 return calc_compressibility_factor();
745}
746
747double AbstractState::fundamental_derivative_of_gas_dynamics() {
748 // See Colonna, FPE, 2010, Eq. 1
749 return 1 + this->second_partial_deriv(iP, iDmass, iSmolar, iDmass, iSmolar) * this->rhomass() / (2 * powInt(speed_sound(), 2));
750};
751
752// Get the derivatives of the parameters in the partial derivative with respect to T and rho
753void get_dT_drho(AbstractState& AS, parameters index, CoolPropDbl& dT, CoolPropDbl& drho) {
754 CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), dT_dtau = -pow(T, 2) / Tr,
755 R = AS.gas_constant(), delta = rho / rhor, tau = Tr / T;
756
757 switch (index) {
758 case iT:
759 dT = 1;
760 drho = 0;
761 break;
762 case iDmolar:
763 dT = 0;
764 drho = 1;
765 break;
766 case iDmass:
767 dT = 0;
768 drho = AS.molar_mass();
769 break;
770 case iP: {
771 // dp/drho|T
772 drho = R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
773 // dp/dT|rho
774 dT = rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau());
775 break;
776 }
777 case iHmass:
778 case iHmolar: {
779 // dh/dT|rho
780 dT = R
781 * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
782 + (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
783 // dh/drhomolar|T
784 drho = T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2());
785 if (index == iHmass) {
786 // dhmolar/drhomolar|T * dhmass/dhmolar where dhmass/dhmolar = 1/mole mass
787 drho /= AS.molar_mass();
788 dT /= AS.molar_mass();
789 }
790 break;
791 }
792 case iSmass:
793 case iSmolar: {
794 // ds/dT|rho
795 dT = R / T * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
796 // ds/drho|T
797 drho = R / rho * (-(1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()));
798 if (index == iSmass) {
799 // ds/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
800 drho /= AS.molar_mass();
801 dT /= AS.molar_mass();
802 }
803 break;
804 }
805 case iUmass:
806 case iUmolar: {
807 // du/dT|rho
808 dT = R * (-pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
809 // du/drho|T
810 drho = AS.T() * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau());
811 if (index == iUmass) {
812 // du/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
813 drho /= AS.molar_mass();
814 dT /= AS.molar_mass();
815 }
816 break;
817 }
818 case iGmass:
819 case iGmolar: {
820 // dg/dT|rho
821 double dTau_dT = 1 / dT_dtau;
822 dT = R * AS.T() * (AS.dalpha0_dTau() + AS.dalphar_dTau() + AS.delta() * AS.d2alphar_dDelta_dTau()) * dTau_dT
823 + R * (1 + AS.alpha0() + AS.alphar() + AS.delta() * AS.dalphar_dDelta());
824 // dg/drho|T
825 double dDelta_drho = 1 / rhor;
826 drho = AS.T() * R * (AS.dalpha0_dDelta() + AS.dalphar_dDelta() + AS.delta() * AS.d2alphar_dDelta2() + AS.dalphar_dDelta()) * dDelta_drho;
827 if (index == iGmass) {
828 // dg/drho|T / drhomass/drhomolar where drhomass/drhomolar = mole mass
829 drho /= AS.molar_mass();
830 dT /= AS.molar_mass();
831 }
832 break;
833 }
834 case iTau:
835 dT = 1 / dT_dtau;
836 drho = 0;
837 break;
838 case iDelta:
839 dT = 0;
840 drho = 1 / rhor;
841 break;
842 case iCvmolar:
843 case iCvmass: {
844 // use the second order derivative of internal energy
845 // make it cleaner by using the function get_dT_drho_second_derivatives directly?
846 // dcvdT|rho = d2u/dT2|rho
847 dT = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
848 // dcvdrho|T = d2u/dT/drho
849 drho = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
850 if (index == iCvmass) {
851 drho /= AS.molar_mass();
852 dT /= AS.molar_mass();
853 }
854 break;
855 }
856 case iCpmolar:
857 case iCpmass: {
858 // dcp/dT|rho = d2h/dT2 + dh/drho * dP/dT * d2P/drhodT / ( dp/drho )^2 - ( d2h/dTdrho * dP/dT + dh/drho * d2P/dT2 ) / ( dP/drho )
859 dT = R / T * pow(tau, 2)
860 * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
861 + delta * AS.d3alphar_dDelta_dTau2());
862 dT += (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
863 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
864 * (R
865 * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
866 - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau()))
867 / pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
868 dT -= ((R / rho * delta
869 * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
870 - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()))
871 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
872 + (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
873 * (rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2())))
874 / (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
875 // dcpdrho|T = d2h/dTdrho + dh/drho * dP/dT * d2P/drho2 / ( dp/drho )^2 - ( d2h/drho2 * dP/dT + dh/drho * d2P/dTdrho ) / ( dP/drho )
876 drho = R / rho * delta
877 * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
878 - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau()); //d2h/dTdrho
879 drho +=
880 (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
881 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
882 * (T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3()))
883 / pow(R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()), 2);
884 drho -= ((R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3()))
885 * (rho * R * (1 + delta * AS.dalphar_dDelta() - tau * delta * AS.d2alphar_dDelta_dTau()))
886 + (T * R / rho * (tau * delta * AS.d2alphar_dDelta_dTau() + delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()))
887 * (R
888 * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()
889 - 2 * delta * tau * AS.d2alphar_dDelta_dTau() - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau())))
890 / (R * T * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2()));
891 if (index == iCpmass) {
892 drho /= AS.molar_mass();
893 dT /= AS.molar_mass();
894 }
895 break;
896 }
897 case ispeed_sound: {
898 //dwdT
899 double aa = 1.0 + delta * AS.dalphar_dDelta() - delta * tau * AS.d2alphar_dDelta_dTau();
900 double bb = pow(tau, 2) * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
901 double daa_dTau = -delta * tau * AS.d3alphar_dDelta_dTau2();
902 double dbb_dTau = pow(tau, 2) * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2.0 * tau * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2());
903 double w = AS.speed_sound();
904 dT = 1.0 / 2.0 / w / T
905 * (pow(w, 2)
906 - R * Tr / AS.molar_mass()
907 * (2.0 * delta * AS.d2alphar_dDelta_dTau() + pow(delta, 2) * AS.d3alphar_dDelta2_dTau()
908 - (2 * aa / bb * daa_dTau - pow(aa / bb, 2) * dbb_dTau)));
909 //dwdrho
910 double daa_dDelta =
911 AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2() - tau * (AS.d2alphar_dDelta_dTau() + delta * AS.d3alphar_dDelta2_dTau());
912 double dbb_dDelta = pow(tau, 2) * (AS.d3alpha0_dDelta_dTau2() + AS.d3alphar_dDelta_dTau2());
913 drho = R * T / 2.0 / AS.molar_mass() / w / rhor
914 * (2.0 * (AS.dalphar_dDelta() + delta * AS.d2alphar_dDelta2())
915 + (2.0 * delta * AS.d2alphar_dDelta2() + pow(delta, 2) * AS.d3alphar_dDelta3())
916 - (2 * aa / bb * daa_dDelta - pow(aa / bb, 2) * dbb_dDelta));
917 break;
918 }
919 default:
920 throw ValueError(format("input to get_dT_drho[%s] is invalid", get_parameter_information(index, "short").c_str()));
921 }
922}
923void get_dT_drho_second_derivatives(AbstractState& AS, int index, CoolPropDbl& dT2, CoolPropDbl& drho_dT, CoolPropDbl& drho2) {
924 CoolPropDbl T = AS.T(), rho = AS.rhomolar(), rhor = AS.rhomolar_reducing(), Tr = AS.T_reducing(), R = AS.gas_constant(), delta = rho / rhor,
925 tau = Tr / T;
926
927 // Here we use T and rho as independent variables since derivations are already done by Thorade, 2013,
928 // Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z
929
930 switch (index) {
931 case iT:
932 case iDmass:
933 case iDmolar:
934 dT2 = 0; // d2rhomolar_dtau2
935 drho2 = 0;
936 drho_dT = 0;
937 break;
938 case iTau:
939 dT2 = 2 * Tr / pow(T, 3);
940 drho_dT = 0;
941 drho2 = 0;
942 break;
943 case iDelta:
944 dT2 = 0;
945 drho_dT = 0;
946 drho2 = 0;
947 break;
948 case iP: {
949 drho2 =
950 T * R / rho * (2 * delta * AS.dalphar_dDelta() + 4 * pow(delta, 2) * AS.d2alphar_dDelta2() + pow(delta, 3) * AS.d3alphar_dDelta3());
951 dT2 = rho * R / T * (pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
952 drho_dT = R
953 * (1 + 2 * delta * AS.dalphar_dDelta() + pow(delta, 2) * AS.d2alphar_dDelta2() - 2 * delta * tau * AS.d2alphar_dDelta_dTau()
954 - tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
955 break;
956 }
957 case iHmass:
958 case iHmolar: {
959 // d2h/drho2|T
960 drho2 = R * T * pow(delta / rho, 2) * (tau * AS.d3alphar_dDelta2_dTau() + 2 * AS.d2alphar_dDelta2() + delta * AS.d3alphar_dDelta3());
961 // d2h/dT2|rho
962 dT2 = R / T * pow(tau, 2)
963 * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2())
964 + delta * AS.d3alphar_dDelta_dTau2());
965 // d2h/drho/dT
966 drho_dT = R / rho * delta
967 * (delta * AS.d2alphar_dDelta2() - pow(tau, 2) * AS.d3alphar_dDelta_dTau2() + AS.dalphar_dDelta()
968 - tau * delta * AS.d3alphar_dDelta2_dTau() - tau * AS.d2alphar_dDelta_dTau());
969 if (index == iHmass) {
970 drho2 /= AS.molar_mass();
971 drho_dT /= AS.molar_mass();
972 dT2 /= AS.molar_mass();
973 }
974 break;
975 }
976 case iSmass:
977 case iSmolar: {
978 // d2s/rho2|T
979 drho2 = R / pow(rho, 2) * (1 - pow(delta, 2) * AS.d2alphar_dDelta2() + tau * pow(delta, 2) * AS.d3alphar_dDelta2_dTau());
980 // d2s/dT2|rho
981 dT2 = R * pow(tau / T, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 3 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
982 // d2s/drho/dT
983 drho_dT = R / (T * rho) * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
984 if (index == iSmass) {
985 drho2 /= AS.molar_mass();
986 drho_dT /= AS.molar_mass();
987 dT2 /= AS.molar_mass();
988 }
989 break;
990 }
991 case iUmass:
992 case iUmolar: {
993 // d2u/rho2|T
994 drho2 = R * T * tau * pow(delta / rho, 2) * AS.d3alphar_dDelta2_dTau();
995 // d2u/dT2|rho
996 dT2 = R / T * pow(tau, 2) * (tau * (AS.d3alpha0_dTau3() + AS.d3alphar_dTau3()) + 2 * (AS.d2alpha0_dTau2() + AS.d2alphar_dTau2()));
997 // d2u/drho/dT
998 drho_dT = R / rho * (-pow(tau, 2) * delta * AS.d3alphar_dDelta_dTau2());
999 if (index == iUmass) {
1000 drho2 /= AS.molar_mass();
1001 drho_dT /= AS.molar_mass();
1002 dT2 /= AS.molar_mass();
1003 }
1004 break;
1005 }
1006 default:
1007 throw ValueError(format("input to get_dT_drho_second_derivatives[%s] is invalid", get_parameter_information(index, "short").c_str()));
1008 }
1009}
1010CoolPropDbl AbstractState::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant) {
1011 CoolPropDbl dOf_dT, dOf_drho, dWrt_dT, dWrt_drho, dConstant_dT, dConstant_drho;
1012
1013 get_dT_drho(*this, Of, dOf_dT, dOf_drho);
1014 get_dT_drho(*this, Wrt, dWrt_dT, dWrt_drho);
1015 get_dT_drho(*this, Constant, dConstant_dT, dConstant_drho);
1016
1017 return (dOf_dT * dConstant_drho - dOf_drho * dConstant_dT) / (dWrt_dT * dConstant_drho - dWrt_drho * dConstant_dT);
1018}
1019CoolPropDbl AbstractState::calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2) {
1020 CoolPropDbl dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2, d2Wrt1_dT2,
1021 d2Wrt1_drhodT, d2Wrt1_drho2, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2, dWrt2_dT, dWrt2_drho, dConstant2_dT, dConstant2_drho, N, D,
1022 dNdrho__T, dDdrho__T, dNdT__rho, dDdT__rho, dderiv1_drho, dderiv1_dT, second;
1023
1024 // First and second partials needed for terms involved in first derivative
1025 get_dT_drho(*this, Of1, dOf1_dT, dOf1_drho);
1026 get_dT_drho(*this, Wrt1, dWrt1_dT, dWrt1_drho);
1027 get_dT_drho(*this, Constant1, dConstant1_dT, dConstant1_drho);
1028 get_dT_drho_second_derivatives(*this, Of1, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2);
1029 get_dT_drho_second_derivatives(*this, Wrt1, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2);
1030 get_dT_drho_second_derivatives(*this, Constant1, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2);
1031
1032 // First derivatives of terms involved in the second derivative
1033 get_dT_drho(*this, Wrt2, dWrt2_dT, dWrt2_drho);
1034 get_dT_drho(*this, Constant2, dConstant2_dT, dConstant2_drho);
1035
1036 // Numerator and denominator of first partial derivative term
1037 N = dOf1_dT * dConstant1_drho - dOf1_drho * dConstant1_dT;
1038 D = dWrt1_dT * dConstant1_drho - dWrt1_drho * dConstant1_dT;
1039
1040 // Derivatives of the numerator and denominator of the first partial derivative term with respect to rho, T held constant
1041 // They are of similar form, with Of1 and Wrt1 swapped
1042 dNdrho__T = dOf1_dT * d2Constant1_drho2 + d2Of1_drhodT * dConstant1_drho - dOf1_drho * d2Constant1_drhodT - d2Of1_drho2 * dConstant1_dT;
1043 dDdrho__T = dWrt1_dT * d2Constant1_drho2 + d2Wrt1_drhodT * dConstant1_drho - dWrt1_drho * d2Constant1_drhodT - d2Wrt1_drho2 * dConstant1_dT;
1044
1045 // Derivatives of the numerator and denominator of the first partial derivative term with respect to T, rho held constant
1046 // They are of similar form, with Of1 and Wrt1 swapped
1047 dNdT__rho = dOf1_dT * d2Constant1_drhodT + d2Of1_dT2 * dConstant1_drho - dOf1_drho * d2Constant1_dT2 - d2Of1_drhodT * dConstant1_dT;
1048 dDdT__rho = dWrt1_dT * d2Constant1_drhodT + d2Wrt1_dT2 * dConstant1_drho - dWrt1_drho * d2Constant1_dT2 - d2Wrt1_drhodT * dConstant1_dT;
1049
1050 // First partial of first derivative term with respect to T
1051 dderiv1_drho = (D * dNdrho__T - N * dDdrho__T) / pow(D, 2);
1052
1053 // First partial of first derivative term with respect to rho
1054 dderiv1_dT = (D * dNdT__rho - N * dDdT__rho) / pow(D, 2);
1055
1056 // Complete second derivative
1057 second = (dderiv1_dT * dConstant2_drho - dderiv1_drho * dConstant2_dT) / (dWrt2_dT * dConstant2_drho - dWrt2_drho * dConstant2_dT);
1058
1059 return second;
1060}
1061// // ----------------------------------------
1062// // Smoothing functions for density
1063// // ----------------------------------------
1064// /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
1065// virtual double AbstractState::drhodh_constp_smoothed(double xend);
1066// /// A smoothed version of the derivative using a spline curve in the region of x=0 to x=xend
1067// virtual double AbstractState::drhodp_consth_smoothed(double xend);
1068// /// Density corresponding to the smoothed derivatives in the region of x=0 to x=xend
1069// virtual void AbstractState::rho_smoothed(double xend, double *rho_spline, double *dsplinedh, double *dsplinedp);
1070
1071} /* namespace CoolProp */
1072
1073#ifdef ENABLE_CATCH
1074
1075#include <catch2/catch_all.hpp>
1076
1077TEST_CASE("Check AbstractState", "[AbstractState]") {
1078 SECTION("bad backend") {
1079 CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("DEFINITELY_A_BAD_BACKEND", "Water")));
1080 }
1081 SECTION("good backend - bad fluid") {
1082 CHECK_THROWS(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "DEFINITELY_A_BAD_FLUID")));
1083 }
1084 SECTION("good backend - helmholtz") {
1085 CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("HEOS", "Water")));
1086 }
1087 SECTION("good backend - incomp") {
1088 CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("INCOMP", "DEB")));
1089 }
1090 SECTION("good backend - REFPROP") {
1091 CHECK_NOTHROW(shared_ptr<CoolProp::AbstractState>(CoolProp::AbstractState::factory("REFPROP", "Water")));
1092 }
1093}
1094
1095TEST_CASE("Check derivatives in first_partial_deriv", "[derivs_in_first_partial_deriv]") {
1096 shared_ptr<CoolProp::AbstractState> Water(CoolProp::AbstractState::factory("HEOS", "Water"));
1097 shared_ptr<CoolProp::AbstractState> WaterplusT(CoolProp::AbstractState::factory("HEOS", "Water"));
1098 shared_ptr<CoolProp::AbstractState> WaterminusT(CoolProp::AbstractState::factory("HEOS", "Water"));
1099 shared_ptr<CoolProp::AbstractState> Waterplusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
1100 shared_ptr<CoolProp::AbstractState> Waterminusrho(CoolProp::AbstractState::factory("HEOS", "Water"));
1101
1102 double dT = 1e-3, drho = 1;
1103 Water->update(CoolProp::PT_INPUTS, 101325, 300);
1104 WaterplusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 + dT);
1105 WaterminusT->update(CoolProp::DmolarT_INPUTS, Water->rhomolar(), 300 - dT);
1106 Waterplusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() + drho, 300);
1107 Waterminusrho->update(CoolProp::DmolarT_INPUTS, Water->rhomolar() - drho, 300);
1108
1109 // Numerical derivatives
1110 CoolPropDbl dP_dT_num = (WaterplusT->p() - WaterminusT->p()) / (2 * dT);
1111 CoolPropDbl dP_drho_num = (Waterplusrho->p() - Waterminusrho->p()) / (2 * drho);
1112
1113 CoolPropDbl dHmolar_dT_num = (WaterplusT->hmolar() - WaterminusT->hmolar()) / (2 * dT);
1114 CoolPropDbl dHmolar_drho_num = (Waterplusrho->hmolar() - Waterminusrho->hmolar()) / (2 * drho);
1115 CoolPropDbl dHmass_dT_num = (WaterplusT->hmass() - WaterminusT->hmass()) / (2 * dT);
1116 CoolPropDbl dHmass_drho_num = (Waterplusrho->hmass() - Waterminusrho->hmass()) / (2 * drho);
1117
1118 CoolPropDbl dSmolar_dT_num = (WaterplusT->smolar() - WaterminusT->smolar()) / (2 * dT);
1119 CoolPropDbl dSmolar_drho_num = (Waterplusrho->smolar() - Waterminusrho->smolar()) / (2 * drho);
1120 CoolPropDbl dSmass_dT_num = (WaterplusT->smass() - WaterminusT->smass()) / (2 * dT);
1121 CoolPropDbl dSmass_drho_num = (Waterplusrho->smass() - Waterminusrho->smass()) / (2 * drho);
1122
1123 CoolPropDbl dUmolar_dT_num = (WaterplusT->umolar() - WaterminusT->umolar()) / (2 * dT);
1124 CoolPropDbl dUmolar_drho_num = (Waterplusrho->umolar() - Waterminusrho->umolar()) / (2 * drho);
1125 CoolPropDbl dUmass_dT_num = (WaterplusT->umass() - WaterminusT->umass()) / (2 * dT);
1126 CoolPropDbl dUmass_drho_num = (Waterplusrho->umass() - Waterminusrho->umass()) / (2 * drho);
1127
1128 CoolPropDbl dGmolar_dT_num = (WaterplusT->gibbsmolar() - WaterminusT->gibbsmolar()) / (2 * dT);
1129 CoolPropDbl dGmolar_drho_num = (Waterplusrho->gibbsmolar() - Waterminusrho->gibbsmolar()) / (2 * drho);
1130 CoolPropDbl dGmass_dT_num = (WaterplusT->gibbsmass() - WaterminusT->gibbsmass()) / (2 * dT);
1131 CoolPropDbl dGmass_drho_num = (Waterplusrho->gibbsmass() - Waterminusrho->gibbsmass()) / (2 * drho);
1132
1133 CoolPropDbl dCvmolar_dT_num = (WaterplusT->cvmolar() - WaterminusT->cvmolar()) / (2 * dT);
1134 CoolPropDbl dCvmolar_drho_num = (Waterplusrho->cvmolar() - Waterminusrho->cvmolar()) / (2 * drho);
1135 CoolPropDbl dCvmass_dT_num = (WaterplusT->cvmass() - WaterminusT->cvmass()) / (2 * dT);
1136 CoolPropDbl dCvmass_drho_num = (Waterplusrho->cvmass() - Waterminusrho->cvmass()) / (2 * drho);
1137
1138 CoolPropDbl dCpmolar_dT_num = (WaterplusT->cpmolar() - WaterminusT->cpmolar()) / (2 * dT);
1139 CoolPropDbl dCpmolar_drho_num = (Waterplusrho->cpmolar() - Waterminusrho->cpmolar()) / (2 * drho);
1140 CoolPropDbl dCpmass_dT_num = (WaterplusT->cpmass() - WaterminusT->cpmass()) / (2 * dT);
1141 CoolPropDbl dCpmass_drho_num = (Waterplusrho->cpmass() - Waterminusrho->cpmass()) / (2 * drho);
1142
1143 CoolPropDbl dspeed_sound_dT_num = (WaterplusT->speed_sound() - WaterminusT->speed_sound()) / (2 * dT);
1144 CoolPropDbl dspeed_sound_drho_num = (Waterplusrho->speed_sound() - Waterminusrho->speed_sound()) / (2 * drho);
1145
1146 // Pressure
1147 CoolPropDbl dP_dT_analyt, dP_drho_analyt;
1148 CoolProp::get_dT_drho(*Water, CoolProp::iP, dP_dT_analyt, dP_drho_analyt);
1149 // Enthalpy
1150 CoolPropDbl dHmolar_dT_analyt, dHmolar_drho_analyt;
1151 CoolProp::get_dT_drho(*Water, CoolProp::iHmolar, dHmolar_dT_analyt, dHmolar_drho_analyt);
1152 CoolPropDbl dHmass_dT_analyt, dHmass_drho_analyt;
1153 CoolProp::get_dT_drho(*Water, CoolProp::iHmass, dHmass_dT_analyt, dHmass_drho_analyt);
1154 // Entropy
1155 CoolPropDbl dSmolar_dT_analyt, dSmolar_drho_analyt;
1156 CoolProp::get_dT_drho(*Water, CoolProp::iSmolar, dSmolar_dT_analyt, dSmolar_drho_analyt);
1157 CoolPropDbl dSmass_dT_analyt, dSmass_drho_analyt;
1158 CoolProp::get_dT_drho(*Water, CoolProp::iSmass, dSmass_dT_analyt, dSmass_drho_analyt);
1159 // Internal energy
1160 CoolPropDbl dUmolar_dT_analyt, dUmolar_drho_analyt;
1161 CoolProp::get_dT_drho(*Water, CoolProp::iUmolar, dUmolar_dT_analyt, dUmolar_drho_analyt);
1162 CoolPropDbl dUmass_dT_analyt, dUmass_drho_analyt;
1163 CoolProp::get_dT_drho(*Water, CoolProp::iUmass, dUmass_dT_analyt, dUmass_drho_analyt);
1164 // Gibbs
1165 CoolPropDbl dGmolar_dT_analyt, dGmolar_drho_analyt;
1166 CoolProp::get_dT_drho(*Water, CoolProp::iGmolar, dGmolar_dT_analyt, dGmolar_drho_analyt);
1167 CoolPropDbl dGmass_dT_analyt, dGmass_drho_analyt;
1168 CoolProp::get_dT_drho(*Water, CoolProp::iGmass, dGmass_dT_analyt, dGmass_drho_analyt);
1169 // Isochoric heat capacity
1170 CoolPropDbl dCvmolar_dT_analyt, dCvmolar_drho_analyt;
1171 CoolProp::get_dT_drho(*Water, CoolProp::iCvmolar, dCvmolar_dT_analyt, dCvmolar_drho_analyt);
1172 CoolPropDbl dCvmass_dT_analyt, dCvmass_drho_analyt;
1173 CoolProp::get_dT_drho(*Water, CoolProp::iCvmass, dCvmass_dT_analyt, dCvmass_drho_analyt);
1174 // Isobaric heat capacity
1175 CoolPropDbl dCpmolar_dT_analyt, dCpmolar_drho_analyt;
1176 CoolProp::get_dT_drho(*Water, CoolProp::iCpmolar, dCpmolar_dT_analyt, dCpmolar_drho_analyt);
1177 CoolPropDbl dCpmass_dT_analyt, dCpmass_drho_analyt;
1178 CoolProp::get_dT_drho(*Water, CoolProp::iCpmass, dCpmass_dT_analyt, dCpmass_drho_analyt);
1179 // Speed of sound
1180 CoolPropDbl dspeed_sound_dT_analyt, dspeed_sound_drho_analyt;
1181 CoolProp::get_dT_drho(*Water, CoolProp::ispeed_sound, dspeed_sound_dT_analyt, dspeed_sound_drho_analyt);
1182
1183 double eps = 1e-3;
1184
1185 CHECK(std::abs(dP_dT_analyt / dP_dT_num - 1) < eps);
1186 CHECK(std::abs(dP_drho_analyt / dP_drho_num - 1) < eps);
1187
1188 CHECK(std::abs(dHmolar_dT_analyt / dHmolar_dT_num - 1) < eps);
1189 CHECK(std::abs(dHmolar_drho_analyt / dHmolar_drho_num - 1) < eps);
1190 CHECK(std::abs(dHmass_dT_analyt / dHmass_dT_num - 1) < eps);
1191 CHECK(std::abs(dHmass_drho_analyt / dHmass_drho_num - 1) < eps);
1192
1193 CHECK(std::abs(dSmolar_dT_analyt / dSmolar_dT_num - 1) < eps);
1194 CHECK(std::abs(dSmolar_drho_analyt / dSmolar_drho_num - 1) < eps);
1195 CHECK(std::abs(dSmass_dT_analyt / dSmass_dT_num - 1) < eps);
1196 CHECK(std::abs(dSmass_drho_analyt / dSmass_drho_num - 1) < eps);
1197
1198 CHECK(std::abs(dUmolar_dT_analyt / dUmolar_dT_num - 1) < eps);
1199 CHECK(std::abs(dUmolar_drho_analyt / dUmolar_drho_num - 1) < eps);
1200 CHECK(std::abs(dUmass_dT_analyt / dUmass_dT_num - 1) < eps);
1201 CHECK(std::abs(dUmass_drho_analyt / dUmass_drho_num - 1) < eps);
1202
1203 CHECK(std::abs(dGmolar_dT_analyt / dGmolar_dT_num - 1) < eps);
1204 CHECK(std::abs(dGmolar_drho_analyt / dGmolar_drho_num - 1) < eps);
1205 CHECK(std::abs(dGmass_dT_analyt / dGmass_dT_num - 1) < eps);
1206 CHECK(std::abs(dGmass_drho_analyt / dGmass_drho_num - 1) < eps);
1207
1208 CHECK(std::abs(dCvmolar_dT_analyt / dCvmolar_dT_num - 1) < eps);
1209 CHECK(std::abs(dCvmolar_drho_analyt / dCvmolar_drho_num - 1) < eps);
1210 CHECK(std::abs(dCvmass_dT_analyt / dCvmass_dT_num - 1) < eps);
1211 CHECK(std::abs(dCvmass_drho_analyt / dCvmass_drho_num - 1) < eps);
1212
1213 CHECK(std::abs(dCpmolar_dT_analyt / dCpmolar_dT_num - 1) < eps);
1214 CHECK(std::abs(dCpmolar_drho_analyt / dCpmolar_drho_num - 1) < eps);
1215 CHECK(std::abs(dCpmass_dT_analyt / dCpmass_dT_num - 1) < eps);
1216 CHECK(std::abs(dCpmass_drho_analyt / dCpmass_drho_num - 1) < eps);
1217
1218 CHECK(std::abs(dspeed_sound_dT_analyt / dspeed_sound_dT_num - 1) < eps);
1219 CHECK(std::abs(dspeed_sound_drho_analyt / dspeed_sound_drho_num - 1) < eps);
1220}
1221
1222#endif
Generated on Mon Apr 21 2025 21:53:40 for CoolProp by doxygen 1.9.4