5 #if defined(ENABLE_CATCH)
8 # include <catch2/catch_all.hpp>
16 if (!type.compare(
"rational_polynomial")) {
17 this->type = TYPE_RATIONAL_POLYNOMIAL;
29 if (!type.compare(
"rhoLnoexp"))
30 this->type = TYPE_NOT_EXPONENTIAL;
32 this->type = TYPE_EXPONENTIAL;
46 if (type == TYPE_NOT_SET) {
48 }
else if (type == TYPE_RATIONAL_POLYNOMIAL) {
52 double THETA = 1 - T /
T_r;
54 for (std::size_t i = 0; i <
N; ++i) {
55 s[i] = n[i] * pow(THETA, t[i]);
57 double summer = std::accumulate(s.begin(), s.end(), 0.0);
59 if (type == TYPE_NOT_EXPONENTIAL) {
64 tau_r_value =
T_r / T;
82 double call(
double T) {
84 return current_value - value;
87 solver_resid resid(
this, value);
88 std::string errstring;
90 min_bound = Tmin - 0.01;
109 for (std::size_t i = 0; i <
simon.
parts.size(); ++i) {
147 throw ValueError(
"Melting line curve not set");
151 }
else if (OF ==
iP_min) {
153 }
else if (OF ==
iT_max) {
155 }
else if (OF ==
iT_min) {
157 }
else if (OF ==
iP && GIVEN ==
iT) {
161 for (std::size_t i = 0; i <
simon.
parts.size(); ++i) {
164 return part.
p_0 + part.
a * (pow(T / part.
T_0, part.
c) - 1);
167 throw ValueError(
"unable to calculate melting line (p,T) for Simon curve");
176 throw ValueError(
"unable to calculate melting line (p,T) for polynomial_in_Tr curve");
185 throw ValueError(
"unable to calculate melting line (p,T) for polynomial_in_Theta curve");
192 for (std::size_t i = 0; i <
simon.
parts.size(); ++i) {
196 if (T >= part.
T_0 && T <= part.
T_max) {
200 throw ValueError(
format(
"unable to calculate melting line T(p) for Simon curve for p=%Lg; bounds are %Lg,%Lg Pa", value,
pmin,
pmax));
208 double call(
double T) {
213 return given_p - calc_p;
221 solver_resid resid(&part, value);
227 format(
"unable to calculate melting line T(p) for polynomial_in_Theta curve for p=%Lg; bounds are %Lg,%Lg Pa", value,
pmin,
pmax));
236 double call(
double T) {
241 return given_p - calc_p;
249 solver_resid resid(&part, value);
256 format(
"unable to calculate melting line T(p) for polynomial_in_Theta curve for p=%Lg; bounds are %Lg,%Lg Pa", value,
pmin,
pmax));
265 #if defined(ENABLE_CATCH)
266 TEST_CASE(
"Water melting line",
"[melting]") {
269 SECTION(
"Ice Ih-liquid") {
270 double actual = AS->melting_line(
iT,
iP, 138.268e6);
271 double expected = 260.0;
274 CHECK(std::abs(actual - expected) < 0.01);
276 SECTION(
"Ice III-liquid") {
277 double actual = AS->melting_line(
iT,
iP, 268.685e6);
278 double expected = 254;
281 CHECK(std::abs(actual - expected) < 0.01);
283 SECTION(
"Ice V-liquid") {
284 double actual = AS->melting_line(
iT,
iP, 479.640e6);
285 double expected = 265;
288 CHECK(std::abs(actual - expected) < 0.01);
290 SECTION(
"Ice VI-liquid") {
291 double actual = AS->melting_line(
iT,
iP, 1356.76e6);
292 double expected = 320;
295 CHECK(std::abs(actual - expected) < 1);
299 TEST_CASE(
"Tests for values from melting lines",
"[melting]") {
302 for (std::size_t i = 0; i < fluids.size(); ++i) {
306 if (!AS->has_melting_line() || !fluids[i].compare(
"Water")) {
316 std::ostringstream ss0;
317 ss0 <<
"Check melting line limits for fluid " << fluids[i];
318 SECTION(ss0.str(),
"") {
327 for (
double p = 0.1 * (pmax - pmin) + pmin; p < pmax; p += 0.2 * (pmax - pmin)) {
329 std::ostringstream ss1;
330 ss1 <<
"Melting line for " << fluids[i] <<
" at p=" << p;
331 SECTION(ss1.str(),
"") {
332 double actual_T = AS->melting_line(
iT,
iP, p);
336 CHECK(actual_T > Tmin);
337 CHECK(actual_T < Tmax);
341 std::ostringstream ss2;
342 ss2 <<
"Ensure melting line valid for " << fluids[i] <<
" @ EOS pmax";
343 SECTION(ss2.str(),
"") {
345 double EOS_pmax = AS->pmax();
347 CHECK_NOTHROW(actual_T = AS->melting_line(
iT,
iP, EOS_pmax));
353 TEST_CASE(
"Test that hs_anchor enthalpy/entropy agrees with EOS",
"[ancillaries]") {
355 for (std::size_t i = 0; i < fluids.size(); ++i) {
361 std::ostringstream ss1;
362 ss1 <<
"Check hs_anchor for " << fluids[i];
363 SECTION(ss1.str(),
"") {
365 "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS";
367 double EOS_hmolar = AS->
hmolar();
368 double EOS_smolar = AS->smolar();
369 CAPTURE(hs_anchor.
hmolar);
370 CAPTURE(hs_anchor.
smolar);
373 CHECK(std::abs(EOS_hmolar - hs_anchor.
hmolar) < 1e-3);
374 CHECK(std::abs(EOS_smolar - hs_anchor.
smolar) < 1e-3);
379 TEST_CASE(
"Surface tension",
"[surface_tension]") {
380 SECTION(
"from PropsSI") {
383 SECTION(
"from saturation_ancillary") {
386 SECTION(
"from AbstractState") {
389 CHECK_NOTHROW(AS->surface_tension());