CoolProp 8.0.0
An open-source fluid property and humid air property database
VLERoutines.cpp
Go to the documentation of this file.
1
2#include "HelmholtzEOSMixtureBackend.h"
3#include "VLERoutines.h"
4#include <numeric>
5
6#include <cmath>
7#include "CoolProp/numerics/MatrixMath.h"
8#include "MixtureDerivatives.h"
9#include "CoolProp/Configuration.h"
10#include "FlashRoutines.h"
11#include <Eigen/Dense>
12
13namespace CoolProp {
14
15void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend& HEOS, parameters ykey, CoolPropDbl y) {
16
17 class inner_resid : public FuncWrapper1D
18 {
19 public:
20 HelmholtzEOSMixtureBackend* HEOS;
21 CoolPropDbl T, desired_p;
22
23 inner_resid(HelmholtzEOSMixtureBackend* HEOS, CoolPropDbl T, CoolPropDbl desired_p) : HEOS(HEOS), T(T), desired_p(desired_p) {};
24 double call(double rhomolar_liq) override {
25 HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
26 CoolPropDbl calc_p = HEOS->SatL->p();
27 std::cout << format("inner p: %0.16Lg; res: %0.16Lg", calc_p, calc_p - desired_p) << '\n';
28 return calc_p - desired_p;
29 }
30 };
31
32 // Define the outer residual to be driven to zero - this is the equality of
33 // Gibbs function for both co-existing phases
34 class outer_resid : public FuncWrapper1D
35 {
36 public:
37 HelmholtzEOSMixtureBackend* HEOS;
38 parameters ykey;
39 CoolPropDbl y;
40 CoolPropDbl rhomolar_crit;
41
42 outer_resid(HelmholtzEOSMixtureBackend& HEOS, CoolProp::parameters ykey, CoolPropDbl y)
43 : HEOS(&HEOS), ykey(ykey), y(y), rhomolar_crit(HEOS.rhomolar_critical()) {};
44 double call(double rhomolar_vap) override {
45 // Calculate the other variable (T->p or p->T) for given vapor density
46 CoolPropDbl T = NAN, p = NAN, rhomolar_liq = NAN;
47 switch (ykey) {
48 case iT: {
49 T = y;
50 HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, y);
51 p = HEOS->SatV->p();
52 std::cout << format("outer p: %0.16Lg", p) << '\n';
53 inner_resid inner(HEOS, T, p);
54 rhomolar_liq = Brent(inner, rhomolar_crit * 1.5, rhomolar_crit * (1 + 1e-8), LDBL_EPSILON, 1e-10, 100);
55 break;
56 }
57 default:
58 throw ValueError("Wrong input for outer_resid");
59 }
60 HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
61 HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, T);
62
63 // Calculate the Gibbs functions for liquid and vapor
64 //CoolPropDbl gL = HEOS->SatL->gibbsmolar();
65 //CoolPropDbl gV = HEOS->SatV->gibbsmolar();
66
67 // Residual is difference in Gibbs function
68 // r = gL - gV;
69
70 return p;
71 };
72 };
73 outer_resid resid(HEOS, iT, y);
74
75 double rhomolar_crit = HEOS.rhomolar_critical();
76
77 Brent(&resid, rhomolar_crit * (1 - 1e-8), rhomolar_crit * 0.5, DBL_EPSILON, 1e-9, 20);
78}
79
80void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options) {
81
82 // Define the residual to be driven to zero
83 class solver_resid : public FuncWrapper1D
84 {
85 public:
86 HelmholtzEOSMixtureBackend* HEOS;
87 CoolPropDbl T, rhomolar_liq, rhomolar_vap;
88
89 solver_resid(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl rhomolar_liq_guess, CoolPropDbl rhomolar_vap_guess)
90 : HEOS(&HEOS), T(T), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess) {};
91 double call(double p) override {
92 // Recalculate the densities using the current guess values
93 HEOS->SatL->update_TP_guessrho(T, p, rhomolar_liq);
94 HEOS->SatV->update_TP_guessrho(T, p, rhomolar_vap);
95
96 // Calculate the Gibbs functions for liquid and vapor
97 CoolPropDbl gL = HEOS->SatL->gibbsmolar();
98 CoolPropDbl gV = HEOS->SatV->gibbsmolar();
99
100 // Residual is difference in Gibbs function
101 return gL - gV;
102 };
103 };
104 solver_resid resid(HEOS, T, options.rhoL, options.rhoV);
105
106 if (!ValidNumber(options.p)) {
107 throw ValueError(format("options.p is not valid in saturation_T_pure_1D_P for T = %Lg", T));
108 };
109 if (!ValidNumber(options.rhoL)) {
110 throw ValueError(format("options.rhoL is not valid in saturation_T_pure_1D_P for T = %Lg", T));
111 };
112 if (!ValidNumber(options.rhoV)) {
113 throw ValueError(format("options.rhoV is not valid in saturation_T_pure_1D_P for T = %Lg", T));
114 };
115
116 try {
117 Secant(resid, options.p, options.p * 1.1, 1e-10, 100);
118 } catch (...) {
119 CoolPropDbl pmax = std::min(options.p * 1.03, static_cast<CoolPropDbl>(HEOS.p_critical() + 1e-6));
120 CoolPropDbl pmin = std::max(options.p * 0.97, static_cast<CoolPropDbl>(HEOS.p_triple() - 1e-6));
121 Brent(resid, pmin, pmax, LDBL_EPSILON, 1e-8, 100);
122 }
123}
124
125void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, saturation_PHSU_pure_options& options) {
126
127 // Define the residual to be driven to zero
128 class solver_resid : public FuncWrapper1D
129 {
130 public:
131 HelmholtzEOSMixtureBackend* HEOS;
132 CoolPropDbl p, rhomolar_liq, rhomolar_vap;
133
134 solver_resid(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, CoolPropDbl rhomolar_liq_guess, CoolPropDbl rhomolar_vap_guess)
135 : HEOS(&HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess) {};
136 double call(double T) override {
137 // Recalculate the densities using the current guess values
138 HEOS->SatL->update_TP_guessrho(T, p, rhomolar_liq);
139 HEOS->SatV->update_TP_guessrho(T, p, rhomolar_vap);
140
141 // Calculate the Gibbs functions for liquid and vapor
142 CoolPropDbl gL = HEOS->SatL->gibbsmolar();
143 CoolPropDbl gV = HEOS->SatV->gibbsmolar();
144
145 // Residual is difference in Gibbs function
146 return gL - gV;
147 };
148 };
149 solver_resid resid(HEOS, p, options.rhoL, options.rhoV);
150
151 if (!ValidNumber(options.T)) {
152 throw ValueError("options.T is not valid in saturation_P_pure_1D_T");
153 };
154 if (!ValidNumber(options.rhoL)) {
155 throw ValueError("options.rhoL is not valid in saturation_P_pure_1D_T");
156 };
157 if (!ValidNumber(options.rhoV)) {
158 throw ValueError("options.rhoV is not valid in saturation_P_pure_1D_T");
159 };
160
161 CoolPropDbl Tmax = std::min(options.T + 2, static_cast<CoolPropDbl>(HEOS.T_critical() - 1e-6));
162 CoolPropDbl Tmin = std::max(options.T - 2, static_cast<CoolPropDbl>(HEOS.Ttriple() + 1e-6));
163 Brent(resid, Tmin, Tmax, LDBL_EPSILON, 1e-11, 100);
164}
165
166void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl specified_value, saturation_PHSU_pure_options& options) {
167 /*
168 This function is inspired by the method of Akasaka:
169
170 R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
171 Helmholtz Energy Equations of State",
172 Journal of Thermal Science and Technology v3 n3,2008
173
174 Ancillary equations are used to get a sensible starting point
175 */
176 std::vector<CoolPropDbl> negativer(3, _HUGE), v;
177 std::vector<std::vector<CoolPropDbl>> J(3, std::vector<CoolPropDbl>(3, _HUGE));
178
179 HEOS.calc_reducing_state();
180 const SimpleState& reduce = HEOS.get_reducing_state();
181 CoolProp::SimpleState crit = HEOS.get_state("reducing");
182 shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
183
184 CoolPropDbl T = NAN, rhoL = NAN, rhoV = NAN, pL = NAN, pV = NAN, hL = NAN, sL = NAN, hV = NAN, sV = NAN;
185 CoolPropDbl deltaL = 0, deltaV = 0, tau = 0, error = NAN;
186 int iter = 0, specified_parameter = 0;
187
188 // Use the density ancillary function as the starting point for the solver
189 try {
190 if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_PL
191 || options.specified_variable == saturation_PHSU_pure_options::IMPOSED_PV) {
192 // Invert liquid density ancillary to get temperature
193 // TODO: fit inverse ancillaries too
194 try {
195 T = HEOS.get_components()[0].ancillaries.pL.invert(specified_value);
196 } catch (...) {
197 throw ValueError("Unable to invert ancillary equation");
198 }
199 } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HL) {
200 CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
201 // Ancillary is deltah = h - hs_anchor.h
202 try {
203 T = HEOS.get_components()[0].ancillaries.hL.invert(specified_value - hs_anchor.hmolar);
204 } catch (...) {
205 throw ValueError("Unable to invert ancillary equation for hL");
206 }
207 } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HV) {
208 class Residual : public FuncWrapper1D
209 {
210 public:
211 CoolPropFluid* component;
212 double h;
213 Residual(CoolPropFluid& component, double h) : component(&component), h(h) {}
214 double call(double T) override {
215 CoolPropDbl h_liq = component->ancillaries.hL.evaluate(T) + component->EOS().hs_anchor.hmolar;
216 return h_liq + component->ancillaries.hLV.evaluate(T) - h;
217 };
218 };
219 Residual resid(HEOS.get_components()[0], HEOS.hmolar());
220
221 // Ancillary is deltah = h - hs_anchor.h
222 CoolPropDbl Tmin_satL = NAN, Tmin_satV = NAN;
223 HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
224 double Tmin = Tmin_satL;
225 double Tmax = HEOS.calc_Tmax_sat();
226 try {
227 T = Brent(resid, Tmin - 3, Tmax + 1, DBL_EPSILON, 1e-10, 50);
228 } catch (...) {
229 shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy = std::make_shared<HelmholtzEOSMixtureBackend>(HEOS.get_components());
230 HEOS_copy->update(QT_INPUTS, 1, Tmin);
231 double hTmin = HEOS_copy->hmolar();
232 HEOS_copy->update(QT_INPUTS, 1, Tmax);
233 double hTmax = HEOS_copy->hmolar();
234 T = (Tmax - Tmin) / (hTmax - hTmin) * (HEOS.hmolar() - hTmin) + Tmin;
235 }
236 } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SL) {
237 CoolPropFluid& component = HEOS.get_components()[0];
238 CoolProp::SaturationAncillaryFunction& anc = component.ancillaries.sL;
239 CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
240 // If near the critical point, use a near critical guess value for T
241 if (std::abs(HEOS.smolar() - crit.smolar) < std::abs(component.ancillaries.sL.get_max_abs_error())) {
242 T = std::max(0.99 * crit.T, crit.T - 0.1);
243 } else {
244 CoolPropDbl Tmin = NAN, Tmax = NAN, Tmin_satV = NAN;
245 HEOS.calc_Tmin_sat(Tmin, Tmin_satV);
246 Tmax = HEOS.calc_Tmax_sat();
247 // Ancillary is deltas = s - hs_anchor.s
248 // First try a conventional call
249 try {
250 T = anc.invert(specified_value - hs_anchor.smolar, Tmin, Tmax);
251 } catch (...) {
252 try {
253 T = anc.invert(specified_value - hs_anchor.smolar, Tmin - 3, Tmax + 3);
254 } catch (...) {
255 double vmin = anc.evaluate(Tmin);
256 double vmax = anc.evaluate(Tmax);
257 if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)) {
258 T = Tmax - 0.1;
259 } else {
260 throw ValueError(format("Unable to invert ancillary equation for sL for value %Lg with Tminval %g and Tmaxval %g ",
261 specified_value - hs_anchor.smolar, vmin, vmax));
262 }
263 }
264 }
265 }
266 } else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SV) {
267 CoolPropFluid& component = HEOS.get_components()[0];
268 CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
269 class Residual : public FuncWrapper1D
270 {
271 public:
272 CoolPropFluid* component;
273 double s;
274 Residual(CoolPropFluid& component, double s) : component(&component), s(s) {}
275 double call(double T) override {
276 CoolPropDbl s_liq = component->ancillaries.sL.evaluate(T) + component->EOS().hs_anchor.smolar;
277 CoolPropDbl resid = s_liq + component->ancillaries.sLV.evaluate(T) - s;
278
279 return resid;
280 };
281 };
282 Residual resid(component, HEOS.smolar());
283
284 // Ancillary is deltas = s - hs_anchor.s
285 CoolPropDbl Tmin_satL = NAN, Tmin_satV = NAN;
286 HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
287 double Tmin = Tmin_satL;
288 double Tmax = HEOS.calc_Tmax_sat();
289 try {
290 T = Brent(resid, Tmin - 3, Tmax, DBL_EPSILON, 1e-10, 50);
291 } catch (...) {
292 CoolPropDbl vmax = resid.call(Tmax);
293 // If near the critical point, use a near critical guess value for T
294 if (std::abs(specified_value - hs_anchor.smolar) < std::abs(vmax)) {
295 T = std::max(0.99 * crit.T, crit.T - 0.1);
296 } else {
297 shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy = std::make_shared<HelmholtzEOSMixtureBackend>(HEOS.get_components());
298 HEOS_copy->update(QT_INPUTS, 1, Tmin);
299 double sTmin = HEOS_copy->smolar();
300 HEOS_copy->update(QT_INPUTS, 1, Tmax);
301 double sTmax = HEOS_copy->smolar();
302 T = (Tmax - Tmin) / (sTmax - sTmin) * (HEOS.smolar() - sTmin) + Tmin;
303 }
304 }
305 } else {
306 throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
307 }
308 // If T from the ancillaries is above the critical temp, this will cause failure
309 // in ancillaries for rhoV and rhoL, decrease if needed
310 T = std::min(T, static_cast<CoolPropDbl>(HEOS.T_critical() - 0.1));
311
312 // Evaluate densities from the ancillary equations
313 rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
314 rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
315
316 // Apply a single step of Newton's method to improve guess value for liquid
317 // based on the error between the gas pressure (which is usually very close already)
318 // and the liquid pressure, which can sometimes (especially at low pressure),
319 // be way off, and often times negative
320 SatL->update(DmolarT_INPUTS, rhoL, T);
321 SatV->update(DmolarT_INPUTS, rhoV, T);
322 double rhoL_updated = rhoL - (SatL->p() - SatV->p()) / SatL->first_partial_deriv(iP, iDmolar, iT);
323
324 // Accept the update if the liquid density is greater than the vapor density
325 if (rhoL_updated > rhoV) {
326 rhoL = rhoL_updated;
327 }
328
329 // Update the state again with the better guess for the liquid density
330 SatL->update(DmolarT_INPUTS, rhoL, T);
331 SatV->update(DmolarT_INPUTS, rhoV, T);
332
333 deltaL = rhoL / reduce.rhomolar;
334 deltaV = rhoV / reduce.rhomolar;
335 tau = reduce.T / T;
336 } catch (NotImplementedError&) {
337 throw; // ??? What is this try...catch for?
338 }
339
340 do {
341 /*if (get_debug_level()>8){
342 std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
343 }*/
344
345 pL = SatL->p();
346 hL = SatL->hmolar();
347 sL = SatL->smolar();
348 pV = SatV->p();
349 hV = SatV->hmolar();
350 sV = SatV->smolar();
351
352 // These derivatives are needed for both cases
353 CoolPropDbl alpharL = SatL->alphar();
354 CoolPropDbl alpharV = SatV->alphar();
355 CoolPropDbl dalphar_dtauL = SatL->dalphar_dTau();
356 CoolPropDbl dalphar_dtauV = SatV->dalphar_dTau();
357 CoolPropDbl d2alphar_ddelta_dtauL = SatL->d2alphar_dDelta_dTau();
358 CoolPropDbl d2alphar_ddelta_dtauV = SatV->d2alphar_dDelta_dTau();
359 CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
360 CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
361 CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
362 CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
363
364 // -r_1 (equate the pressures)
365 negativer[0] = -(deltaV * (1 + deltaV * dalphar_ddeltaV) - deltaL * (1 + deltaL * dalphar_ddeltaL));
366 // -r_2 (equate the gibbs energy)
367 negativer[1] = -(deltaV * dalphar_ddeltaV + alpharV + log(deltaV) - deltaL * dalphar_ddeltaL - alpharL - log(deltaL));
368 switch (options.specified_variable) {
369 case saturation_PHSU_pure_options::IMPOSED_PL:
370 // -r_3 (equate calculated pressure and specified liquid pressure)
371 negativer[2] = -(pL / specified_value - 1);
372 break;
373 case saturation_PHSU_pure_options::IMPOSED_PV:
374 // -r_3 (equate calculated pressure and specified vapor pressure)
375 negativer[2] = -(pV / specified_value - 1);
376 break;
377 case saturation_PHSU_pure_options::IMPOSED_HL:
378 // -r_3 (equate calculated liquid enthalpy and specified liquid enthalpy)
379 negativer[2] = -(hL - specified_value);
380 break;
381 case saturation_PHSU_pure_options::IMPOSED_HV:
382 // -r_3 (equate calculated vapor enthalpy and specified vapor enthalpy)
383 negativer[2] = -(hV - specified_value);
384 break;
385 case saturation_PHSU_pure_options::IMPOSED_SL:
386 // -r_3 (equate calculated liquid entropy and specified liquid entropy)
387 negativer[2] = -(sL - specified_value);
388 break;
389 case saturation_PHSU_pure_options::IMPOSED_SV:
390 // -r_3 (equate calculated vapor entropy and specified vapor entropy)
391 negativer[2] = -(sV - specified_value);
392 break;
393 default:
394 throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
395 }
396
397 // dr1_dtau
398 J[0][0] = pow(deltaV, 2) * d2alphar_ddelta_dtauV - pow(deltaL, 2) * d2alphar_ddelta_dtauL;
399 // dr2_dtau
400 J[1][0] = deltaV * d2alphar_ddelta_dtauV + dalphar_dtauV - deltaL * d2alphar_ddelta_dtauL - dalphar_dtauL;
401
402 if (options.use_logdelta) {
403 // dr_1/d_log(delta'')
404 J[0][1] = -deltaL - 2 * pow(deltaL, 2) * dalphar_ddeltaL - pow(deltaL, 3) * d2alphar_ddelta2L;
405 // dr_2/d_log(delta'')
406 J[1][1] = -pow(deltaL, 2) * d2alphar_ddelta2L - 2 * deltaL * dalphar_ddeltaL - 1;
407 } else {
408 // dr_1/ddelta''
409 J[0][1] = -1 - 2 * deltaL * dalphar_ddeltaL - pow(deltaL, 2) * d2alphar_ddelta2L;
410 // dr_2/ddelta''
411 J[1][1] = -1 / deltaL - 2 * dalphar_ddeltaL - deltaL * d2alphar_ddelta2L;
412 }
413
414 if (options.use_logdelta) {
415 // dr_1/d_log(delta'')
416 J[0][2] = deltaV + 2 * pow(deltaV, 2) * dalphar_ddeltaV + pow(deltaV, 3) * d2alphar_ddelta2V;
417 // dr_2/d_log(delta'')
418 J[1][2] = 1 + 2 * deltaV * dalphar_ddeltaV + 1 + pow(deltaV, 2) * d2alphar_ddelta2V;
419 } else {
420 // dr_1/ddelta''
421 J[0][2] = 1 + 2 * deltaV * dalphar_ddeltaV + pow(deltaV, 2) * d2alphar_ddelta2V;
422 // dr_2/ddelta''
423 J[1][2] = deltaV * d2alphar_ddelta2V + 2 * dalphar_ddeltaV + 1 / deltaV;
424 }
425
426 // Derivatives of the specification equation
427 switch (options.specified_variable) {
428 case saturation_PHSU_pure_options::IMPOSED_PL:
429 // dr_3/dtau
430 J[2][0] = SatL->first_partial_deriv(iP, iTau, iDelta) / specified_value;
431 if (options.use_logdelta) {
432 // dr_3/d(log(delta'))
433 J[2][1] = deltaL * SatL->first_partial_deriv(iP, iDelta, iTau) / specified_value;
434 } else {
435 // dr_3/ddelta'
436 J[2][1] = SatL->first_partial_deriv(iP, iDelta, iTau) / specified_value;
437 }
438 // dr_3/ddelta'' (liquid pressure not a function of vapor density)
439 J[2][2] = 0;
440 specified_parameter = CoolProp::iP;
441 break;
442 case saturation_PHSU_pure_options::IMPOSED_PV:
443 // dr_3/dtau
444 J[2][0] = SatV->first_partial_deriv(iP, iTau, iDelta) / specified_value;
445 // dr_3/ddelta' (vapor pressure not a function of liquid density)
446 J[2][1] = 0;
447 if (options.use_logdelta) {
448 // dr_3/d(log(delta'')
449 J[2][2] = deltaV * SatV->first_partial_deriv(iP, iDelta, iTau) / specified_value;
450 } else {
451 // dr_3/ddelta''
452 J[2][2] = SatV->first_partial_deriv(iP, iDelta, iTau) / specified_value;
453 }
454 specified_parameter = CoolProp::iP;
455 break;
456 case saturation_PHSU_pure_options::IMPOSED_HL:
457 // dr_3/dtau
458 J[2][0] = SatL->first_partial_deriv(iHmolar, iTau, iDelta);
459 // dr_3/ddelta'
460 J[2][1] = SatL->first_partial_deriv(iHmolar, iDelta, iTau);
461 if (options.use_logdelta) {
462 J[2][1] *= deltaL;
463 }
464 // dr_3/ddelta''
465 J[2][2] = 0; //(liquid enthalpy not a function of vapor density)
466 specified_parameter = CoolProp::iHmolar;
467 break;
468 case saturation_PHSU_pure_options::IMPOSED_HV:
469 // dr_3/dtau
470 J[2][0] = SatV->first_partial_deriv(iHmolar, iTau, iDelta);
471 // dr_3/ddelta'
472 J[2][1] = 0; //(vapor enthalpy not a function of liquid density)
473 // dr_3/ddelta''
474 J[2][2] = SatV->first_partial_deriv(iHmolar, iDelta, iTau);
475 if (options.use_logdelta) {
476 J[2][2] *= deltaV;
477 }
478 specified_parameter = CoolProp::iHmolar;
479 break;
480 case saturation_PHSU_pure_options::IMPOSED_SL:
481 // dr_3/dtau
482 J[2][0] = SatL->first_partial_deriv(iSmolar, iTau, iDelta);
483 // dr_3/ddelta'
484 J[2][1] = SatL->first_partial_deriv(iSmolar, iDelta, iTau);
485 if (options.use_logdelta) {
486 J[2][1] *= deltaL;
487 }
488 // dr_3/ddelta''
489 J[2][2] = 0; //(liquid entropy not a function of vapor density)
490 specified_parameter = CoolProp::iSmolar;
491 break;
492 case saturation_PHSU_pure_options::IMPOSED_SV:
493 // dr_3/dtau
494 J[2][0] = SatV->first_partial_deriv(iSmolar, iTau, iDelta);
495 // dr_3/ddelta'
496 J[2][1] = 0; //(vapor enthalpy not a function of liquid density)
497 // dr_3/ddelta''
498 J[2][2] = SatV->first_partial_deriv(iSmolar, iDelta, iTau);
499 if (options.use_logdelta) {
500 J[2][2] *= deltaV;
501 }
502 specified_parameter = CoolProp::iSmolar;
503 break;
504 default:
505 throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
506 }
507
508 v = linsolve(J, negativer);
509
510 // Conditions for an acceptable step are:
511 // a) tau > 1
512 // b) rhoL > rhoV or deltaL > deltaV
513 double tau0 = tau, deltaL0 = deltaL, deltaV0 = deltaV;
514 // Geometric damping search (~25 iters) — no FP accumulation.
515 for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) { // NOLINT(cert-flp30-c)
516 tau = tau0 + omega_local * options.omega * v[0];
517 if (options.use_logdelta) {
518 deltaL = exp(log(deltaL0) + omega_local * options.omega * v[1]);
519 deltaV = exp(log(deltaV0) + omega_local * options.omega * v[2]);
520 } else {
521 deltaL = deltaL0 + omega_local * options.omega * v[1];
522 deltaV = deltaV0 + omega_local * options.omega * v[2];
523 }
524 if (tau > 1 && deltaL > deltaV) {
525 break;
526 }
527 }
528
529 rhoL = deltaL * reduce.rhomolar;
530 rhoV = deltaV * reduce.rhomolar;
531 T = reduce.T / tau;
532
533 SatL->update(DmolarT_INPUTS, rhoL, T);
534 SatV->update(DmolarT_INPUTS, rhoV, T);
535
536 error = sqrt(pow(negativer[0], 2) + pow(negativer[1], 2) + pow(negativer[2], 2));
537 iter++;
538 if (T < 0) {
539 throw SolutionError(format("saturation_PHSU_pure solver T < 0"));
540 }
541 // If the change is very small, stop
542 if (max_abs_value(v) < 1e-10) {
543 break;
544 }
545 if (iter > 50) {
546 // Set values back into the options structure for use in next solver
547 options.rhoL = rhoL;
548 options.rhoV = rhoV;
549 options.T = T;
550 // Error out
551 std::string info = get_parameter_information(specified_parameter, "short");
552 throw SolutionError(format("saturation_PHSU_pure solver did not converge after 50 iterations for %s=%Lg current error is %Lg",
553 info.c_str(), specified_value, error));
554 }
555 } while (error > 1e-9);
556 // Recalculate error
557 // The result has changed since the last error calculation.
558 // In rare scenarios, the final step can become unstable due to solving a singular
559 // J matrix. This final error check verifies that the solution is still good.
560 // Furthermore, the forced phase of SatL and SatV may have caused errors. We will recalculate them without this assumption.
561 SatL->unspecify_phase();
562 SatV->unspecify_phase();
563 SatL->update(DmolarT_INPUTS, rhoL, T);
564 SatV->update(DmolarT_INPUTS, rhoV, T);
565 negativer[0] = -(deltaV * (1 + deltaV * SatV->dalphar_dDelta()) - deltaL * (1 + deltaL * SatL->dalphar_dDelta()));
566 negativer[1] = -(deltaV * SatV->dalphar_dDelta() + SatV->alphar() + log(deltaV) - deltaL * SatL->dalphar_dDelta() - SatL->alphar() - log(deltaL));
567 switch (options.specified_variable) {
568 case saturation_PHSU_pure_options::IMPOSED_PL:
569 // -r_3 (equate calculated pressure and specified liquid pressure)
570 negativer[2] = -(SatL->p() / specified_value - 1);
571 break;
572 case saturation_PHSU_pure_options::IMPOSED_PV:
573 // -r_3 (equate calculated pressure and specified vapor pressure)
574 negativer[2] = -(SatV->p() / specified_value - 1);
575 break;
576 case saturation_PHSU_pure_options::IMPOSED_HL:
577 // -r_3 (equate calculated liquid enthalpy and specified liquid enthalpy)
578 negativer[2] = -(SatL->hmolar() - specified_value);
579 break;
580 case saturation_PHSU_pure_options::IMPOSED_HV:
581 // -r_3 (equate calculated vapor enthalpy and specified vapor enthalpy)
582 negativer[2] = -(SatV->hmolar() - specified_value);
583 break;
584 case saturation_PHSU_pure_options::IMPOSED_SL:
585 // -r_3 (equate calculated liquid entropy and specified liquid entropy)
586 negativer[2] = -(SatL->smolar() - specified_value);
587 break;
588 case saturation_PHSU_pure_options::IMPOSED_SV:
589 // -r_3 (equate calculated vapor entropy and specified vapor entropy)
590 negativer[2] = -(SatV->smolar() - specified_value);
591 break;
592 default:
593 throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
594 }
595 error = sqrt(pow(negativer[0], 2) + pow(negativer[1], 2) + pow(negativer[2], 2));
596 // reset the phase for the next update.
597 SatL->specify_phase(iphase_liquid);
598 SatV->specify_phase(iphase_gas);
599 if (error > 1e-8 && max_abs_value(v) > 1e-9) {
600 throw SolutionError(format("saturation_PHSU_pure solver was good, but went bad. Current error is %Lg", error));
601 }
602}
603void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options) {
604 /*
605 This function is inspired by the method of Akasaka:
606
607 R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
608 Helmholtz Energy Equations of State",
609 Journal of Thermal Science and Technology v3 n3,2008
610
611 Ancillary equations are used to get a sensible starting point
612 */
613 std::vector<CoolPropDbl> r(2, _HUGE), v;
614 std::vector<std::vector<CoolPropDbl>> J(2, std::vector<CoolPropDbl>(2, _HUGE));
615
616 HEOS.calc_reducing_state();
617 const SimpleState& reduce = HEOS.get_reducing_state();
618 shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
619
620 CoolPropDbl T = NAN, rhoL = NAN, rhoV = NAN;
621 CoolPropDbl deltaL = 0, deltaV = 0, tau = 0, error = NAN, p_error = NAN;
622 int iter = 0;
623
624 // Use the density ancillary function as the starting point for the solver
625 try {
626 if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
627 // Invert liquid density ancillary to get temperature
628 // TODO: fit inverse ancillaries too
629 T = HEOS.get_components()[0].ancillaries.rhoL.invert(rhomolar);
630 rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
631 rhoL = rhomolar;
632 } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) {
633 // Invert vapor density ancillary to get temperature
634 // TODO: fit inverse ancillaries too
635 T = HEOS.get_components()[0].ancillaries.rhoV.invert(rhomolar);
636 rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
637 rhoV = rhomolar;
638 } else {
639 throw ValueError(format("imposed rho to saturation_D_pure [%d%] is invalid", options.imposed_rho));
640 }
641
642 deltaL = rhoL / reduce.rhomolar;
643 deltaV = rhoV / reduce.rhomolar;
644 tau = reduce.T / T;
645 } catch (NotImplementedError&) {
646 throw; // ??? What is this try...catch for?
647 }
648
649 do {
650 /*if (get_debug_level()>8){
651 std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
652 }*/
653
654 // Calculate once to save on calls to EOS
655 SatL->update(DmolarT_INPUTS, rhoL, T);
656 SatV->update(DmolarT_INPUTS, rhoV, T);
657
658 CoolPropDbl pL = SatL->p();
659 CoolPropDbl pV = SatV->p();
660
661 // These derivatives are needed for both cases
662 CoolPropDbl dalphar_dtauL = SatL->dalphar_dTau();
663 CoolPropDbl dalphar_dtauV = SatV->dalphar_dTau();
664 CoolPropDbl d2alphar_ddelta_dtauL = SatL->d2alphar_dDelta_dTau();
665 CoolPropDbl d2alphar_ddelta_dtauV = SatV->d2alphar_dDelta_dTau();
666 CoolPropDbl alpharL = SatL->alphar();
667 CoolPropDbl alpharV = SatV->alphar();
668 CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
669 CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
670
671 // -r_1
672 r[0] = -(deltaV * (1 + deltaV * dalphar_ddeltaV) - deltaL * (1 + deltaL * dalphar_ddeltaL));
673 // -r_2
674 r[1] = -(deltaV * dalphar_ddeltaV + alpharV + log(deltaV) - deltaL * dalphar_ddeltaL - alpharL - log(deltaL));
675
676 // dr1_dtau
677 J[0][0] = pow(deltaV, 2) * d2alphar_ddelta_dtauV - pow(deltaL, 2) * d2alphar_ddelta_dtauL;
678 // dr2_dtau
679 J[1][0] = deltaV * d2alphar_ddelta_dtauV + dalphar_dtauV - deltaL * d2alphar_ddelta_dtauL - dalphar_dtauL;
680
681 if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
682 CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
683 if (options.use_logdelta) {
684 J[0][1] = deltaV + 2 * pow(deltaV, 2) * dalphar_ddeltaV + pow(deltaV, 3) * d2alphar_ddelta2V;
685 J[1][1] = pow(deltaV, 2) * d2alphar_ddelta2V + 2 * deltaV * dalphar_ddeltaV + 1;
686 } else {
687 J[0][1] = 1 + 2 * deltaV * dalphar_ddeltaV + pow(deltaV, 2) * d2alphar_ddelta2V;
688 J[1][1] = deltaV * d2alphar_ddelta2V + 2 * dalphar_ddeltaV + 1 / deltaV;
689 }
690 } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) {
691 CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
692 if (options.use_logdelta) {
693 J[0][1] = -deltaL - 2 * pow(deltaL, 2) * dalphar_ddeltaL - pow(deltaL, 3) * d2alphar_ddelta2L;
694 J[1][1] = -pow(deltaL, 2) * d2alphar_ddelta2L - 2 * deltaL * dalphar_ddeltaL - 1;
695 } else {
696 J[0][1] = -1 - 2 * deltaL * dalphar_ddeltaL - pow(deltaL, 2) * d2alphar_ddelta2L;
697 J[1][1] = -deltaL * d2alphar_ddelta2L - 2 * dalphar_ddeltaL - 1 / deltaL;
698 }
699 }
700
701 //double DET = J[0][0]*J[1][1]-J[0][1]*J[1][0];
702
703 v = linsolve(J, r);
704
705 double omega = options.omega;
706
707 if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOL) {
708 if (options.use_logdelta)
709 deltaV = exp(log(deltaV) + omega * v[1]);
710 else {
711 if (deltaV + omega * v[1] <= 0) {
712 omega = 0.5 * deltaV / v[1];
713 } // gone off track, take a smaller step
714 if (tau + omega * v[0] <= 0) {
715 omega = 0.5 * tau / v[0];
716 }
717 deltaV += omega * v[1];
718 }
719 } else {
720 if (options.use_logdelta)
721 deltaL = exp(log(deltaL) + omega * v[1]);
722 else {
723 if (deltaL + omega * v[1] <= 0) {
724 omega = 0.5 * deltaL / v[1];
725 } // gone off track, take a smaller step
726 if (tau + omega * v[0] <= 0) {
727 omega = 0.5 * tau / v[0];
728 }
729 deltaL += omega * v[1];
730 }
731 }
732
733 tau += omega * v[0];
734
735 rhoL = deltaL * reduce.rhomolar;
736 rhoV = deltaV * reduce.rhomolar;
737 T = reduce.T / tau;
738
739 p_error = (pL - pV) / pL;
740
741 error = sqrt(pow(r[0], 2) + pow(r[1], 2));
742 iter++;
743 if (T < 0) {
744 throw SolutionError(format("saturation_D_pure solver T < 0"));
745 }
746 if (iter > options.max_iterations) {
747 throw SolutionError(
748 format("saturation_D_pure solver did not converge after %d iterations with rho: %g mol/m^3", options.max_iterations, rhomolar));
749 }
750 } while (error > 1e-9);
751 CoolPropDbl p_error_limit = 1e-3;
752 if (std::abs(p_error) > p_error_limit) {
753 throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit));
754 }
755}
756void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options) {
757 // Set some input options
758 SaturationSolvers::saturation_T_pure_Akasaka_options _options(false);
759 _options.omega = 1.0;
760 try {
761 // Actually call the solver
762 SaturationSolvers::saturation_T_pure_Maxwell(HEOS, T, _options);
763 } catch (...) {
764 try {
765 // Actually call the solver
766 SaturationSolvers::saturation_T_pure_Akasaka(HEOS, T, _options);
767 } catch (...) {
768 // If there was an error, store values for use in later solvers
769 options.pL = _options.pL;
770 options.pV = _options.pV;
771 options.rhoL = _options.rhoL;
772 options.rhoV = _options.rhoV;
773 options.p = _options.pL;
774 SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options);
775 }
776 }
777}
778void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options) {
779 // Start with the method of Akasaka
780
781 /*
782 This function implements the method of Akasaka
783
784 R. Akasaka,"A Reliable and Useful Method to Determine the Saturation State from
785 Helmholtz Energy Equations of State",
786 Journal of Thermal Science and Technology v3 n3,2008
787
788 Ancillary equations are used to get a sensible starting point
789 */
790
791 HEOS.calc_reducing_state();
792 const SimpleState& reduce = HEOS.get_reducing_state();
793 CoolPropDbl R_u = HEOS.gas_constant();
794 shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
795
796 CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, JL = NAN, JV = NAN, KL = NAN, KV = NAN, dJL = NAN, dJV = NAN, dKL = NAN, dKV = NAN;
797 CoolPropDbl DELTA = NAN, deltaL = 0, deltaV = 0, error = NAN, PL = NAN, PV = NAN, stepL = NAN, stepV = NAN;
798 int iter = 0;
799
800 try {
801 if (options.use_guesses) {
802 // Use the guesses provided in the options structure
803 rhoL = options.rhoL;
804 rhoV = options.rhoV;
805 } else {
806 // Use the density ancillary function as the starting point for the solver
807
808 // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
809 if (T > 0.99 * HEOS.get_reducing_state().T) {
810 rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T - 0.1);
811 rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T - 0.1);
812 } else {
813 rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
814 rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
815
816 // Apply a single step of Newton's method to improve guess value for liquid
817 // based on the error between the gas pressure (which is usually very close already)
818 // and the liquid pressure, which can sometimes (especially at low pressure),
819 // be way off, and often times negative
820 SatL->update(DmolarT_INPUTS, rhoL, T);
821 SatV->update(DmolarT_INPUTS, rhoV, T);
822
823 // Update the guess for liquid density using density solver with vapor pressure
824 // and liquid density guess from ancillaries
825 HEOS.specify_phase(iphase_liquid);
826 rhoL = HEOS.solver_rho_Tp(T, SatV->p(), rhoL);
827 HEOS.unspecify_phase();
828 }
829 }
830
831 deltaL = rhoL / reduce.rhomolar;
832 deltaV = rhoV / reduce.rhomolar;
833 } catch (NotImplementedError&) { // NOLINT(bugprone-empty-catch)
834 // Backend doesn't implement the saturation-density ancillaries
835 // (e.g. PCSAFT, incompressible) — keep the deltaL/deltaV initial
836 // guess from the caller and let the Newton iteration below
837 // converge from there. The commented-out Soave fallback was an
838 // earlier attempt at a guess-from-Tc/pc/omega path; left in
839 // place as a hint if anyone revisits this.
840 /*double Tc = crit.T;
841 double pc = crit.p.Pa;
842 double w = 6.67228479e-09*Tc*Tc*Tc-7.20464352e-06*Tc*Tc+3.16947758e-03*Tc-2.88760012e-01;
843 double q = -6.08930221451*w -5.42477887222;
844 double pt = exp(q*(Tc/T-1))*pc;*/
845
846 //double rhoL = density_Tp_Soave(T, pt, 0), rhoV = density_Tp_Soave(T, pt, 1);
847
848 //deltaL = rhoL/reduce.rhomolar;
849 //deltaV = rhoV/reduce.rhomolar;
850 //tau = reduce.T/T;
851 }
852 //if (get_debug_level()>5){
853 // std::cout << format("%s:%d: Akasaka guess values deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
854 // }
855
856 do {
857 /*if (get_debug_level()>8){
858 std::cout << format("%s:%d: right before the derivs with deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
859 }*/
860
861 // Calculate once to save on calls to EOS
862 SatL->update(DmolarT_INPUTS, rhoL, T);
863 SatV->update(DmolarT_INPUTS, rhoV, T);
864 CoolPropDbl alpharL = SatL->alphar();
865 CoolPropDbl alpharV = SatV->alphar();
866 CoolPropDbl dalphar_ddeltaL = SatL->dalphar_dDelta();
867 CoolPropDbl dalphar_ddeltaV = SatV->dalphar_dDelta();
868 CoolPropDbl d2alphar_ddelta2L = SatL->d2alphar_dDelta2();
869 CoolPropDbl d2alphar_ddelta2V = SatV->d2alphar_dDelta2();
870
871 JL = deltaL * (1 + deltaL * dalphar_ddeltaL);
872 JV = deltaV * (1 + deltaV * dalphar_ddeltaV);
873 KL = deltaL * dalphar_ddeltaL + alpharL + log(deltaL);
874 KV = deltaV * dalphar_ddeltaV + alpharV + log(deltaV);
875
876 PL = R_u * reduce.rhomolar * T * JL;
877 PV = R_u * reduce.rhomolar * T * JV;
878
879 // At low pressure, the magnitude of d2alphar_ddelta2L and d2alphar_ddelta2V are enormous, truncation problems arise for all the partials
880 dJL = 1 + 2 * deltaL * dalphar_ddeltaL + deltaL * deltaL * d2alphar_ddelta2L;
881 dJV = 1 + 2 * deltaV * dalphar_ddeltaV + deltaV * deltaV * d2alphar_ddelta2V;
882 dKL = 2 * dalphar_ddeltaL + deltaL * d2alphar_ddelta2L + 1 / deltaL;
883 dKV = 2 * dalphar_ddeltaV + deltaV * d2alphar_ddelta2V + 1 / deltaV;
884
885 DELTA = dJV * dKL - dJL * dKV;
886
887 error = sqrt((KL - KV) * (KL - KV) + (JL - JV) * (JL - JV));
888
889 // Get the predicted step
890 stepL = options.omega / DELTA * ((KV - KL) * dJV - (JV - JL) * dKV);
891 stepV = options.omega / DELTA * ((KV - KL) * dJL - (JV - JL) * dKL);
892
893 CoolPropDbl deltaL0 = deltaL, deltaV0 = deltaV;
894 // Conditions for an acceptable step are:
895 // a) rhoL > rhoV or deltaL > deltaV
896 // Geometric damping search (~25 iters) — no FP accumulation.
897 for (double omega_local = 1.0; omega_local > 0.1; omega_local /= 1.1) { // NOLINT(cert-flp30-c)
898 deltaL = deltaL0 + omega_local * stepL;
899 deltaV = deltaV0 + omega_local * stepV;
900
901 if (deltaL > 1 && deltaV < 1 && deltaV > 0) {
902 break;
903 }
904 }
905
906 rhoL = deltaL * reduce.rhomolar;
907 rhoV = deltaV * reduce.rhomolar;
908 iter++;
909 if (iter > 100) {
910 throw SolutionError(format("Akasaka solver did not converge after 100 iterations"));
911 }
912 } while (error > 1e-10 && std::abs(stepL) > 10 * DBL_EPSILON * std::abs(stepL) && std::abs(stepV) > 10 * DBL_EPSILON * std::abs(stepV));
913
914 CoolPropDbl p_error_limit = 1e-3;
915 CoolPropDbl p_error = (PL - PV) / PL;
916 if (std::abs(p_error) > p_error_limit) {
917 options.pL = PL;
918 options.pV = PV;
919 options.rhoL = rhoL;
920 options.rhoV = rhoV;
921 throw SolutionError(format("saturation_T_pure_Akasaka solver abs error on p [%g] > limit [%g]", std::abs(p_error), p_error_limit));
922 }
923}
924
925CoolPropDbl sign(CoolPropDbl x) {
926 if (x > 0) {
927 return 1;
928 } else {
929 return -1;
930 }
931}
932
933void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options) {
934
935 /*
936 This function implements the method of
937
938 Ancillary equations are used to get a sensible starting point
939 */
940
941 HEOS.calc_reducing_state();
942 shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL, SatV = HEOS.SatV;
943 CoolProp::SimpleState& crit = HEOS.get_components()[0].crit;
944 CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL = NAN, DeltavV = NAN, pL = NAN, pV = NAN, p = NAN, last_error = NAN;
945 int iter = 0, small_step_count = 0,
946 backwards_step_count = 0; // Counter for the number of times you have taken a step that increases error
947
948 try {
949 if (options.use_guesses) {
950 // Use the guesses provided in the options structure
951 rhoL = options.rhoL;
952 rhoV = options.rhoV;
953 } else {
954 // Use the density ancillary function as the starting point for the solver
955
956 // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
957 if (T > 0.9999 * HEOS.get_reducing_state().T) {
958 rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T - 0.1);
959 rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T - 0.1);
960 } else {
961 rhoL = HEOS.get_components()[0].ancillaries.rhoL.evaluate(T);
962 rhoV = HEOS.get_components()[0].ancillaries.rhoV.evaluate(T);
963 p = HEOS.get_components()[0].ancillaries.pV.evaluate(T);
964
965 CoolProp::SimpleState& tripleL = HEOS.get_components()[0].triple_liquid;
966 CoolProp::SimpleState& tripleV = HEOS.get_components()[0].triple_vapor;
967
968 // If the guesses are terrible, apply a simple correction
969 // but only if the limits are being checked
970 if ((rhoL < crit.rhomolar * 0.8 || rhoL > tripleL.rhomolar * 1.2 || rhoV > crit.rhomolar * 1.2 || rhoV < tripleV.rhomolar * 0.8)
971 && !get_config_bool(DONT_CHECK_PROPERTY_LIMITS)) {
972
973 if (get_debug_level() > 5) {
974 std::cout << format(
975 "[Maxwell] ancillaries correction T: %0.16Lg rhoL: %0.16Lg rhoV: %0.16Lg rhoc: %g rhoLtrip: %g rhoVtrip: %g\n", T, rhoL,
976 rhoV, crit.rhomolar, tripleL.rhomolar, tripleV.rhomolar);
977 }
978
979 // Lets assume that liquid density is more or less linear with T
980 rhoL = (crit.rhomolar - tripleL.rhomolar) / (crit.T - tripleL.T) * (T - tripleL.T) + tripleL.rhomolar;
981 // Then we calculate pressure from this density
982 SatL->update_DmolarT_direct(rhoL, T);
983 // Then we assume vapor to be ideal gas
984 if (SatL->p() > 0) {
985 rhoV = SatL->p() / (SatL->gas_constant() * T);
986 } else {
987 rhoV = p / (SatL->gas_constant() * T);
988 }
989 // Update the vapor state
990 SatV->update_DmolarT_direct(rhoV, T);
991 } else {
992 SatL->update_DmolarT_direct(rhoL, T);
993 SatV->update_DmolarT_direct(rhoV, T);
994 }
995 if (get_debug_level() > 5) {
996 std::cout << format("[Maxwell] ancillaries T: %0.16Lg rhoL: %0.16Lg rhoV: %0.16Lg pL: %g pV: %g\n", T, rhoL, rhoV, SatL->p(),
997 SatV->p());
998 }
999
1000 // Update the guess for liquid density using density solver with vapor pressure
1001 // and liquid density guess from ancillaries, but only if the pressures are not
1002 // close to each other
1003 if (std::abs((SatL->p() - p) / p) > 0.1) {
1004 for (int iii = 0; iii < 6; ++iii) {
1005 // Use Halley's method to update the liquid density (http://en.wikipedia.org/wiki/Halley%27s_method)
1006 CoolPropDbl f = SatL->p() - SatV->p();
1007 CoolPropDbl dpdrho = SatL->first_partial_deriv(iP, iDmolar, iT);
1008 CoolPropDbl d2pdrho2 = SatL->second_partial_deriv(iP, iDmolar, iT, iDmolar, iT);
1009 CoolPropDbl deltarhoLHalley = -(2 * f * dpdrho) / (2 * POW2(dpdrho) - f * d2pdrho2);
1010 rhoL += deltarhoLHalley;
1011 if (std::abs(deltarhoLHalley / rhoL) < DBL_EPSILON) {
1012 break;
1013 }
1014 SatL->update_DmolarT_direct(rhoL, T);
1015 }
1016 }
1017
1018 SatL->update_DmolarT_direct(rhoL, T);
1019 SatV->update_DmolarT_direct(rhoV, T);
1020
1021 // Update the guess for vapor density using density solver with vapor pressure
1022 // and density guess from ancillaries, but only if the pressures are not
1023 // close to each other
1024 if (std::abs((SatV->p() - p) / p) > 0.1) {
1025 HEOS.specify_phase(iphase_gas);
1026 rhoV = SatV->solver_rho_Tp(T, p, rhoV);
1027 HEOS.unspecify_phase();
1028 }
1029 }
1030 }
1031 } catch (NotImplementedError&) { // NOLINT(bugprone-empty-catch)
1032 // SatL/SatV ancillaries not implemented for this backend — keep
1033 // the input rhoL/rhoV initial guesses and let the caller's
1034 // crit.rhomolar clamps below handle the dome edge.
1035 }
1036
1037 if (rhoL < crit.rhomolar) {
1038 rhoL = 1.01 * crit.rhomolar;
1039 }
1040 if (rhoV > crit.rhomolar) {
1041 rhoV = 0.99 * crit.rhomolar;
1042 }
1043 last_error = _HUGE;
1044 SatL->update_DmolarT_direct(rhoL, T);
1045 SatV->update_DmolarT_direct(rhoV, T);
1046 if (get_debug_level() > 5) {
1047 std::cout << format("[Maxwell] starting T: %0.16Lg rhoL: %Lg rhoV: %Lg pL: %Lg pV: %g\n", T, rhoL, rhoV, SatL->p(), SatV->p());
1048 }
1049 do {
1050 pL = SatL->p();
1051 pV = SatV->p();
1052 CoolPropDbl vL = 1 / SatL->rhomolar(), vV = 1 / SatV->rhomolar();
1053 // Get alpha, the pressure derivative with volume at constant T
1054 // Given by (dp/drho|T)*drhodv
1055 CoolPropDbl alphaL = SatL->first_partial_deriv(iP, iDmolar, iT) * (-POW2(SatL->rhomolar()));
1056 CoolPropDbl alphaV = SatV->first_partial_deriv(iP, iDmolar, iT) * (-POW2(SatV->rhomolar()));
1057
1058 // Total helmholtz energy for liquid and vapor
1059 CoolPropDbl RT = SatL->gas_constant() * T;
1060 CoolPropDbl helmholtzL = (SatL->calc_alpha0() + SatL->calc_alphar()) * RT;
1061 CoolPropDbl helmholtzV = (SatV->calc_alpha0() + SatV->calc_alphar()) * RT;
1062
1063 // Calculate the mean pressure
1064 CoolPropDbl pM = (helmholtzL - helmholtzV) / (vV - vL);
1065
1066 // Coefficients for the quadratic in the step
1067 CoolPropDbl A = 0.5 * alphaL * (alphaL - alphaV);
1068 CoolPropDbl B = alphaL * (pL - pV - alphaV * (vL - vV));
1069 CoolPropDbl C = alphaV * (vL - vV) * (pM - pL) + 0.5 * POW2(pL - pV);
1070
1071 // Argument to square root
1072 CoolPropDbl sqrt_arg = std::abs(B * B / (4 * A * A) - C / A);
1073
1074 // If the argument to sqrt is very small, we multiply it by a large factor to make it
1075 // larger, and then also divide the sqrt by the sqrt of the factor
1076 if (std::abs(sqrt_arg) > 1e-10) {
1077 DeltavL = -0.5 * B / A + sign((alphaL - alphaV) / alphaV) * sqrt(sqrt_arg);
1078 } else {
1079 // Scale the argument to sqrt() function to make it about 1.0, and divide by sqrt(factor) to yield a factor of 1
1080 CoolPropDbl powerfactor = -log10(sqrt_arg);
1081 DeltavL = -0.5 * B / A + sign((alphaL - alphaV) / alphaV) * sqrt(sqrt_arg * powerfactor) / sqrt(powerfactor);
1082 }
1083 DeltavV = (pL - pV + alphaL * DeltavL) / alphaV;
1084
1085 // Update the densities of liquid and vapor
1086 rhoL = 1 / (vL + DeltavL);
1087 rhoV = 1 / (vV + DeltavV);
1088 if (B * B / (4 * A * A) - C / A < -10 * DBL_EPSILON) {
1089 rhoL *= 1.01;
1090 rhoV /= 1.01;
1091 }
1092
1093 // Update the states again
1094 SatL->update_DmolarT_direct(rhoL, T);
1095 SatV->update_DmolarT_direct(rhoV, T);
1096
1097 // Calculate the error (here the relative error in pressure)
1098 error = std::abs((SatL->p() - SatV->p()) / SatL->p());
1099
1100 if (get_debug_level() > 5) {
1101 std::cout << format("[Maxwell] rhoL: %0.16Lg rhoV: %0.16Lg error: %Lg dvL/vL: %Lg dvV/vV: %Lg pL: %Lg pV: %Lg\n", rhoL, rhoV, error,
1102 DeltavL / vL, DeltavV / vV, pL, pV);
1103 }
1104
1105 // If the step size is small, start a counter to allow the other density
1106 // to be corrected a few times
1107 if (std::abs(DeltavL * rhoL) < 1e-13 || std::abs(DeltavV * rhoV) < 1e-13) {
1108 small_step_count++;
1109 }
1110 // If you are not continuing to march towards the solution, after a couple of times, stop
1111 // This is especially a problem for water
1112 if (std::abs(error) > std::abs(last_error)) {
1113 backwards_step_count++;
1114 }
1115
1116 iter++;
1117 last_error = error;
1118 if (iter > 30) {
1119 throw SolutionError(format("Maxwell solver did not converge after 30 iterations; rhoL: %0.16Lg rhoV: %0.16Lg error: %Lg dvL/vL: %Lg "
1120 "dvV/vV: %Lg pL: %Lg pV: %Lg\n",
1121 rhoL, rhoV, error, DeltavL / vL, DeltavV / vV, pL, pV));
1122 }
1123 } while ((SatL->p() < 0) || (error > 1e-10 && small_step_count < 4 && backwards_step_count < 6));
1124 if (get_debug_level() > 5) {
1125 std::cout << format("[Maxwell] pL: %g pV: %g\n", SatL->p(), SatV->p());
1126 }
1127}
1128
1129void SaturationSolvers::x_and_y_from_K(CoolPropDbl beta, const std::vector<CoolPropDbl>& K, const std::vector<CoolPropDbl>& z,
1130 std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y) {
1131 for (unsigned int i = 0; i < K.size(); i++) {
1132 double denominator = (1 - beta + beta * K[i]); // Common denominator
1133 x[i] = z[i] / denominator;
1134 y[i] = K[i] * z[i] / denominator;
1135 }
1136}
1137
1138void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend& HEOS, const CoolPropDbl beta, CoolPropDbl T, CoolPropDbl p,
1139 const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K, mixture_VLE_IO& options) {
1140 int iter = 1;
1141 CoolPropDbl change = NAN, f = NAN, df = NAN, deriv_liq = NAN, deriv_vap = NAN;
1142 std::size_t N = z.size();
1143 std::vector<CoolPropDbl> ln_phi_liq, ln_phi_vap;
1144 ln_phi_liq.resize(N);
1145 ln_phi_vap.resize(N);
1146
1147 std::vector<CoolPropDbl>&x = HEOS.SatL->get_mole_fractions_ref(), &y = HEOS.SatV->get_mole_fractions_ref();
1148 x_and_y_from_K(beta, K, z, x, y);
1149
1150 HEOS.SatL->specify_phase(iphase_liquid);
1151 HEOS.SatV->specify_phase(iphase_gas);
1152
1153 normalize_vector(x);
1154 normalize_vector(y);
1155
1156 HEOS.SatL->set_mole_fractions(x);
1157 HEOS.SatV->set_mole_fractions(y);
1158 HEOS.SatL->calc_reducing_state();
1159 HEOS.SatV->calc_reducing_state();
1160 CoolPropDbl rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(T, p, iphase_liquid); // [mol/m^3]
1161 CoolPropDbl rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(T, p, iphase_gas); // [mol/m^3]
1162
1163 // Use Peneloux volume translation to shift liquid volume
1164 // As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
1165 double summer_c = 0, v_SRK = 1 / rhomolar_liq;
1166 const std::vector<CoolPropFluid>& components = HEOS.get_components();
1167 for (std::size_t i = 0; i < components.size(); ++i) {
1168 // Get the parameters for the cubic EOS
1169 CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical);
1170 CoolPropDbl pc = HEOS.get_fluid_constant(i, iP_critical);
1171 CoolPropDbl rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
1172 CoolPropDbl R = 8.3144598;
1173
1174 summer_c += z[i] * (0.40768 * R * Tc / pc * (0.29441 - pc / (rhomolarc * R * Tc)));
1175 }
1176 rhomolar_liq = 1 / (v_SRK - summer_c);
1177 HEOS.SatL->update_TP_guessrho(T, p, rhomolar_liq);
1178 HEOS.SatV->update_TP_guessrho(T, p, rhomolar_vap);
1179
1180 do {
1181 HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
1182 HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
1183
1184 f = 0;
1185 df = 0;
1186
1187 x_N_dependency_flag xN_flag = XN_INDEPENDENT;
1188 for (std::size_t i = 0; i < N; ++i) {
1189 ln_phi_liq[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i, xN_flag);
1190 ln_phi_vap[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i, xN_flag);
1191
1192 if (options.sstype == imposed_p) {
1193 deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatL.get()), i, xN_flag);
1194 deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatV.get()), i, xN_flag);
1195 } else if (options.sstype == imposed_T) {
1196 deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatL.get()), i, xN_flag);
1197 deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatV.get()), i, xN_flag);
1198 } else {
1199 throw ValueError();
1200 }
1201
1202 K[i] = exp(ln_phi_liq[i] - ln_phi_vap[i]);
1203
1204 f += z[i] * (K[i] - 1) / (1 - beta + beta * K[i]);
1205
1206 double dfdK = K[i] * z[i] / pow(1 - beta + beta * K[i], (int)2);
1207 df += dfdK * (deriv_liq - deriv_vap);
1208 }
1209
1210 if (std::abs(df) <= 1e-14) { // To avoid dividing by 0
1211 if (std::abs(f) <= 1e-12) // 1e-12 is the loop convergence criterion
1212 {
1213 change = -f; // Should be converged. f <= e-12, so change will have nearly no impact.
1214 } else {
1215 throw ValueError(format("df very small (df = %g) in successive_substitution but f is not converged (f = %g > 1e-12).", df, f));
1216 }
1217 } else {
1218 change = -f / df;
1219 }
1220
1221 double omega = 1.0;
1222 if (options.sstype == imposed_p) {
1223 T += change;
1224 } else if (options.sstype == imposed_T) {
1225 if (std::abs(change) > 0.05 * p) {
1226 omega = 0.1;
1227 }
1228 p += omega * change;
1229 }
1230
1231 x_and_y_from_K(beta, K, z, x, y);
1232 normalize_vector(x);
1233 normalize_vector(y);
1234 HEOS.SatL->set_mole_fractions(x);
1235 HEOS.SatV->set_mole_fractions(y);
1236
1237 iter += 1;
1238 if (iter > 50) {
1239 throw ValueError(format("saturation_p was unable to reach a solution within 50 iterations"));
1240 }
1241 } while (std::abs(f) > 1e-12 && iter < options.Nstep_max);
1242
1243 HEOS.SatL->update_TP_guessrho(T, p, HEOS.SatL->rhomolar());
1244 HEOS.SatV->update_TP_guessrho(T, p, HEOS.SatV->rhomolar());
1245
1246 options.p = HEOS.SatL->p();
1247 options.T = HEOS.SatL->T();
1248 options.rhomolar_liq = HEOS.SatL->rhomolar();
1249 options.rhomolar_vap = HEOS.SatV->rhomolar();
1250 options.x = x;
1251 options.y = y;
1252}
1253void SaturationSolvers::newton_raphson_saturation::resize(std::size_t N) {
1254 this->N = N;
1255 x.resize(N);
1256 y.resize(N);
1257
1258 if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1259 r.resize(N + 1);
1260 err_rel.resize(N + 1);
1261 J.resize(N + 1, N + 1);
1262 } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1263 r.resize(N);
1264 err_rel.resize(N);
1265 J.resize(N, N);
1266 } else {
1267 throw ValueError();
1268 }
1269}
1270void SaturationSolvers::newton_raphson_saturation::check_Jacobian() {
1271 // References to the classes for concision
1272 HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1273
1274 // Build the Jacobian and residual vectors
1275 build_arrays();
1276
1277 // Make copies of the base
1278 CoolPropDbl T0 = T;
1279 std::vector<CoolPropDbl> x0 = x;
1280 Eigen::VectorXd r0 = r;
1281 Eigen::MatrixXd J0 = J;
1282 CoolPropDbl rhomolar_liq0 = rSatL.rhomolar();
1283 CoolPropDbl rhomolar_vap0 = rSatV.rhomolar();
1284
1285 {
1286 // Derivatives with respect to T
1287 double dT = 1e-3, T1 = T + dT, T2 = T - dT;
1288 this->T = T1;
1289 this->rhomolar_liq = rhomolar_liq0;
1290 this->rhomolar_vap = rhomolar_vap0;
1291 build_arrays();
1292 Eigen::VectorXd r1 = r;
1293 this->T = T2;
1294 this->rhomolar_liq = rhomolar_liq0;
1295 this->rhomolar_vap = rhomolar_vap0;
1296 build_arrays();
1297 Eigen::VectorXd r2 = r;
1298
1299 Eigen::VectorXd diffn = (r1 - r2) / (2 * dT);
1300 std::cout << format("For T\n");
1301 //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1302 //std::cout << "analytic: " << vec_to_string(J0.col(N-1), "%0.11Lg") << std::endl;
1303 }
1304 {
1305 // Derivatives with respect to rho'
1306 double drho = 1;
1307 this->T = T0;
1308 this->rhomolar_liq = rhomolar_liq0 + drho;
1309 this->rhomolar_vap = rhomolar_vap0;
1310 build_arrays();
1311 Eigen::VectorXd rr1 = r;
1312 this->T = T0;
1313 this->rhomolar_liq = rhomolar_liq0 - drho;
1314 this->rhomolar_vap = rhomolar_vap0;
1315 build_arrays();
1316 Eigen::VectorXd rr2 = r;
1317
1318 Eigen::VectorXd diffn = (rr1 - rr2) / (2 * drho);
1319 std::cout << format("For rho\n");
1320 //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1321 //std::cout << "analytic: " << vec_to_string(J0.col(N-1), "%0.11Lg") << std::endl;
1322 }
1323 for (std::size_t i = 0; i < x.size() - 1; ++i) {
1324 // Derivatives with respect to x[i]
1325 double dx = 1e-5;
1326 this->x = x0;
1327 this->x[i] += dx;
1328 this->x[x.size() - 1] -= dx;
1329 this->T = T0;
1330 this->rhomolar_liq = rhomolar_liq0;
1331 this->rhomolar_vap = rhomolar_vap0;
1332 build_arrays();
1333 Eigen::VectorXd r1 = this->r;
1334
1335 this->x = x0;
1336 this->x[i] -= dx;
1337 this->x[x.size() - 1] += dx;
1338 this->T = T0;
1339 this->rhomolar_liq = rhomolar_liq0;
1340 this->rhomolar_vap = rhomolar_vap0;
1341 build_arrays();
1342 Eigen::VectorXd r2 = this->r;
1343
1344 Eigen::VectorXd diffn = (r1 - r2) / (2 * dx);
1345 std::cout << format("For x%d N %d\n", i, N);
1346 //std::cout << "numerical: " << vec_to_string(diffn, "%0.11Lg") << std::endl;
1347 //std::cout << "analytic: " << vec_to_string(J0.col(i), "%0.11Lg") << std::endl;
1348 }
1349}
1350void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& z,
1351 std::vector<CoolPropDbl>& z_incipient, newton_raphson_saturation_options& IO) {
1352 int iter = 0;
1353 bool debug = get_debug_level() > 9 || false;
1354
1355 if (debug) {
1356 std::cout << " NRsat::call: p " << IO.p << " T " << IO.T << " dl " << IO.rhomolar_liq << " dv " << IO.rhomolar_vap << '\n';
1357 }
1358
1359 // Reset all the variables and resize
1360 pre_call();
1361
1362 this->bubble_point = IO.bubble_point;
1363 rhomolar_liq = IO.rhomolar_liq;
1364 rhomolar_vap = IO.rhomolar_vap;
1365 T = IO.T;
1366 p = IO.p;
1367 imposed_variable = IO.imposed_variable;
1368
1369 resize(z.size());
1370
1371 if (bubble_point) {
1372 // Bubblepoint, vapor (y) is the incipient phase
1373 x = z;
1374 y = z_incipient;
1375 } else {
1376 // Dewpoint, liquid (x) is the incipient phase
1377 y = z;
1378 x = z_incipient;
1379 }
1380
1381 // Hold a pointer to the backend
1382 this->HEOS = &HEOS;
1383
1384 //check_Jacobian();
1385
1386 do {
1387 // Build the Jacobian and residual vectors
1388 build_arrays();
1389
1390 // Solve for the step; v is the step with the contents
1391 // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
1392 Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1393
1394 if (bubble_point) {
1395 for (unsigned int i = 0; i < N - 1; ++i) {
1396 err_rel[i] = v[i] / y[i];
1397 y[i] += v[i];
1398 }
1399 y[N - 1] = 1 - std::accumulate(y.begin(), y.end() - 1, 0.0);
1400 } else {
1401 for (unsigned int i = 0; i < N - 1; ++i) {
1402 err_rel[i] = v[i] / x[i];
1403 x[i] += v[i];
1404 }
1405 x[N - 1] = 1 - std::accumulate(x.begin(), x.end() - 1, 0.0);
1406 }
1407 if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED) {
1408 T += v[N - 1];
1409 err_rel[N - 1] = v[N - 1] / T;
1410 } else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1411 p += v[N - 1];
1412 err_rel[N - 1] = v[N - 1] / p;
1413 } else if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1414 T += v[N - 1];
1415 err_rel[N - 1] = v[N - 1] / T;
1416 rhomolar_liq += v[N];
1417 err_rel[N] = v[N] / rhomolar_liq;
1418 } else {
1419 throw ValueError("invalid imposed_variable");
1420 }
1421 if (debug) {
1422 //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
1423 }
1424
1425 min_rel_change = err_rel.cwiseAbs().minCoeff();
1426 iter++;
1427
1428 if (iter == IO.Nstep_max) {
1429 throw ValueError(format("newton_raphson_saturation::call reached max number of iterations [%d]", IO.Nstep_max));
1430 }
1431 } while (this->error_rms > 1e-7 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1432
1433 IO.Nsteps = iter;
1434 IO.p = p;
1435 IO.x = x; // Mole fractions in liquid
1436 IO.y = y; // Mole fractions in vapor
1437 IO.T = T;
1438 IO.rhomolar_liq = rhomolar_liq;
1439 IO.rhomolar_vap = rhomolar_vap;
1440 const std::vector<CoolPropFluid>& fluidsL = HEOS.SatL->get_components();
1441 const std::vector<CoolPropFluid>& fluidsV = HEOS.SatV->get_components();
1442 if (!fluidsL.empty() && !fluidsV.empty()) {
1443 IO.hmolar_liq = HEOS.SatL->hmolar();
1444 IO.hmolar_vap = HEOS.SatV->hmolar();
1445 IO.smolar_liq = HEOS.SatL->smolar();
1446 IO.smolar_vap = HEOS.SatV->smolar();
1447 }
1448}
1449
1450void SaturationSolvers::newton_raphson_saturation::build_arrays() {
1451 // References to the classes for concision
1452 HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1453
1454 // Step 0:
1455 // -------
1456 // Set mole fractions for the incipient phase
1457 if (bubble_point) {
1458 // Vapor is incipient phase, set its composition
1459 rSatV.set_mole_fractions(y);
1460 rSatL.set_mole_fractions(x);
1461 } else {
1462 // Liquid is incipient phase, set its composition
1463 rSatL.set_mole_fractions(x);
1464 rSatV.set_mole_fractions(y);
1465 }
1466
1467 if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1468 rSatL.update(DmolarT_INPUTS, rhomolar_liq, T);
1469 rSatV.update(DmolarT_INPUTS, rhomolar_vap, T);
1470 } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED || imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1471 rSatL.update_TP_guessrho(T, p, rhomolar_liq);
1472 rhomolar_liq = rSatL.rhomolar();
1473 rSatV.update_TP_guessrho(T, p, rhomolar_vap);
1474 rhomolar_vap = rSatV.rhomolar();
1475 } else {
1476 throw ValueError("imposed variable not set for NR VLE");
1477 }
1478
1479 // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
1480 CoolPropDbl p_liq = rSatL.p();
1481 CoolPropDbl p_vap = rSatV.p();
1482 p = 0.5 * (p_liq + p_vap);
1483
1484 // Step 2:
1485 // -------
1486 // Build the residual vector and the Jacobian matrix
1487
1488 x_N_dependency_flag xN_flag = XN_DEPENDENT;
1489
1490 if (imposed_variable == newton_raphson_saturation_options::RHOV_IMPOSED) {
1491 // For the residuals F_i (equality of fugacities)
1492 for (std::size_t i = 0; i < N; ++i) {
1493 // Equate the liquid and vapor fugacities
1494 CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1495 CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1496 r(i) = ln_f_liq - ln_f_vap;
1497
1498 for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1499 if (bubble_point) {
1500 J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatV, i, j, xN_flag);
1501 } else {
1502 J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(rSatL, i, j, xN_flag);
1503 }
1504 }
1505 J(i, N - 1) = MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatL, i, xN_flag)
1506 - MixtureDerivatives::dln_fugacity_i_dT__constrho_n(rSatV, i, xN_flag);
1507 J(i, N) = MixtureDerivatives::dln_fugacity_i_drho__constT_n(rSatL, i, xN_flag);
1508 }
1509 // ---------------------------------------------------------------
1510 // Derivatives of pL(T,rho',x)-p(T,rho'',y) with respect to inputs
1511 // ---------------------------------------------------------------
1512 r(N) = p_liq - p_vap;
1513 for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1514 J(N, j) = MixtureDerivatives::dpdxj__constT_V_xi(rSatL, j, xN_flag); // p'' not a function of x0
1515 }
1516 // Fixed composition derivatives
1517 J(N, N - 1) = rSatL.first_partial_deriv(iP, iT, iDmolar) - rSatV.first_partial_deriv(iP, iT, iDmolar);
1518 J(N, N) = rSatL.first_partial_deriv(iP, iDmolar, iT);
1519 } else if (imposed_variable == newton_raphson_saturation_options::P_IMPOSED) {
1520 // Independent variables are N-1 mole fractions of incipient phase and T
1521
1522 // For the residuals F_i (equality of fugacities)
1523 for (std::size_t i = 0; i < N; ++i) {
1524 // Equate the liquid and vapor fugacities
1525 CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1526 CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1527 r(i) = ln_f_liq - ln_f_vap;
1528
1529 for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1530 if (bubble_point) {
1531 J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1532 } else {
1533 J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1534 }
1535 }
1536 J(i, N - 1) =
1537 MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
1538 }
1539 } else if (imposed_variable == newton_raphson_saturation_options::T_IMPOSED) {
1540 // Independent variables are N-1 mole fractions of incipient phase and p
1541
1542 // For the residuals F_i (equality of fugacities)
1543 for (std::size_t i = 0; i < N; ++i) {
1544 // Equate the liquid and vapor fugacities
1545 CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1546 CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1547 r(i) = ln_f_liq - ln_f_vap;
1548
1549 for (std::size_t j = 0; j < N - 1; ++j) { // j from 0 to N-2
1550 if (bubble_point) {
1551 J(i, j) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1552 } else {
1553 J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1554 }
1555 }
1556 J(i, N - 1) =
1557 MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
1558 }
1559 } else {
1560 throw ValueError();
1561 }
1562
1563 error_rms = r.norm();
1564
1565 // Calculate derivatives along phase boundary;
1566 // Gernert thesis 3.96 and 3.97
1567 CoolPropDbl dQ_dPsat = 0, dQ_dTsat = 0;
1568 for (std::size_t i = 0; i < N; ++i) {
1569 dQ_dPsat += x[i]
1570 * (MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatL, i, xN_flag)
1571 - MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(rSatV, i, xN_flag));
1572 dQ_dTsat += x[i]
1573 * (MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatL, i, xN_flag)
1574 - MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(rSatV, i, xN_flag));
1575 }
1576 dTsat_dPsat = -dQ_dPsat / dQ_dTsat;
1577 dPsat_dTsat = -dQ_dTsat / dQ_dPsat;
1578}
1579
1580void SaturationSolvers::newton_raphson_twophase::call(HelmholtzEOSMixtureBackend& HEOS, newton_raphson_twophase_options& IO) {
1581 int iter = 0;
1582
1583 if (get_debug_level() > 9) {
1584 std::cout << " NRsat::call: p" << IO.p << " T" << IO.T << " dl" << IO.rhomolar_liq << " dv" << IO.rhomolar_vap << '\n';
1585 }
1586
1587 // Reset all the variables and resize
1588 pre_call();
1589
1590 rhomolar_liq = IO.rhomolar_liq;
1591 rhomolar_vap = IO.rhomolar_vap;
1592 T = IO.T;
1593 p = IO.p;
1594 imposed_variable = IO.imposed_variable;
1595 x = IO.x;
1596 y = IO.y;
1597 z = IO.z;
1598 beta = IO.beta;
1599
1600 this->N = z.size();
1601 x.resize(N);
1602 y.resize(N);
1603 r.resize(2 * N - 1);
1604 J.resize(2 * N - 1, 2 * N - 1);
1605 err_rel.resize(2 * N - 1);
1606
1607 // Hold a pointer to the backend
1608 this->HEOS = &HEOS;
1609
1610 do {
1611 // Build the Jacobian and residual vectors
1612 build_arrays();
1613
1614 // A previous step may have left the composition/density infeasible, making the
1615 // fugacity (hence the residual) non-finite. Fail loudly so the caller falls back
1616 // to the blind solver, rather than letting "NaN > tol == false" silently exit the
1617 // loop and return a stale, unconverged state (GH #3192 follow-up).
1618 if (!ValidNumber(error_rms)) {
1619 throw ValueError("newton_raphson_twophase::call produced a non-finite residual");
1620 }
1621
1622 // Solve for the step; v is the step with the contents
1623 // [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]
1624 Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
1625
1626 // Damp the Newton step so neither phase composition leaves the open interval (0,1).
1627 // An undamped step routinely overshoots (x_i < 0 or y_i > 1); the next fugacity
1628 // evaluation then returns NaN and the solve diverges (GH #3192 follow-up). tau is
1629 // the largest fraction of the full step that keeps every stepped mole fraction
1630 // inside (0,1), backed off by step_safety so it stays strictly interior.
1631 double tau = 1.0;
1632 const double step_safety = 0.8;
1633 // The dependent mole fractions x[N-1] = 1 - sum(x[0..N-2]) and y[N-1] likewise change by
1634 // minus the sum of the independent steps; accumulate those so the dependent component is
1635 // bounded to (0,1) too (otherwise it can still overshoot for N>=3 and reintroduce the
1636 // non-finite fugacity path this damping is meant to prevent).
1637 double dx_last = 0.0, dy_last = 0.0;
1638 for (std::size_t i = 0; i < N - 1; ++i) {
1639 const double dx = v[i], dy = v[i + (N - 1)];
1640 dx_last -= dx;
1641 dy_last -= dy;
1642 if (x[i] + dx <= 0.0) {
1643 tau = std::min(tau, step_safety * (-x[i] / dx));
1644 } else if (x[i] + dx >= 1.0) {
1645 tau = std::min(tau, step_safety * ((1.0 - x[i]) / dx));
1646 }
1647 if (y[i] + dy <= 0.0) {
1648 tau = std::min(tau, step_safety * (-y[i] / dy));
1649 } else if (y[i] + dy >= 1.0) {
1650 tau = std::min(tau, step_safety * ((1.0 - y[i]) / dy));
1651 }
1652 }
1653 // Bound the dependent (Nth) mole fraction with the same rule.
1654 const double x_last = x[N - 1], y_last = y[N - 1];
1655 if (x_last + dx_last <= 0.0) {
1656 tau = std::min(tau, step_safety * (-x_last / dx_last));
1657 } else if (x_last + dx_last >= 1.0) {
1658 tau = std::min(tau, step_safety * ((1.0 - x_last) / dx_last));
1659 }
1660 if (y_last + dy_last <= 0.0) {
1661 tau = std::min(tau, step_safety * (-y_last / dy_last));
1662 } else if (y_last + dy_last >= 1.0) {
1663 tau = std::min(tau, step_safety * ((1.0 - y_last) / dy_last));
1664 }
1665
1666 for (unsigned int i = 0; i < N - 1; ++i) {
1667 err_rel[i] = tau * v[i] / x[i];
1668 x[i] += tau * v[i];
1669 err_rel[i + (N - 1)] = tau * v[i + (N - 1)] / y[i];
1670 y[i] += tau * v[i + (N - 1)];
1671 }
1672 x[N - 1] = 1 - std::accumulate(x.begin(), x.end() - 1, 0.0);
1673 y[N - 1] = 1 - std::accumulate(y.begin(), y.end() - 1, 0.0);
1674
1675 if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED) {
1676 T += tau * v[2 * N - 2];
1677 err_rel[2 * N - 2] = tau * v[2 * N - 2] / T;
1678 } else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED) {
1679 p += tau * v[2 * N - 2];
1680 err_rel[2 * N - 2] = tau * v[2 * N - 2] / p;
1681 } else {
1682 throw ValueError("invalid imposed_variable");
1683 }
1684
1685 // Stop only when the LARGEST relative change is small (maxCoeff, not minCoeff: a
1686 // single near-stationary variable must not terminate the iteration prematurely).
1687 min_rel_change = err_rel.cwiseAbs().maxCoeff();
1688 iter++;
1689
1690 if (iter == IO.Nstep_max) {
1691 throw ValueError(format("newton_raphson_twophase::call reached max number of iterations [%d]", IO.Nstep_max));
1692 }
1693 } while (this->error_rms > 1e-9 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
1694
1695 // Refresh the residual at the final iterate (this also leaves SatL/SatV holding the
1696 // converged state) and require genuine convergence; otherwise signal failure so the
1697 // caller can fall back to the blind solver instead of returning the last guess.
1698 build_arrays();
1699 if (!ValidNumber(error_rms) || error_rms > 1e-7) {
1700 throw ValueError(format("newton_raphson_twophase::call did not converge (error_rms = %g)", static_cast<double>(error_rms)));
1701 }
1702
1703 IO.Nsteps = iter;
1704 IO.p = p;
1705 IO.x = x; // Mole fractions in liquid
1706 IO.y = y; // Mole fractions in vapor
1707 IO.T = T;
1708 IO.rhomolar_liq = rhomolar_liq;
1709 IO.rhomolar_vap = rhomolar_vap;
1710 IO.hmolar_liq = HEOS.SatL.get()->hmolar();
1711 IO.hmolar_vap = HEOS.SatV.get()->hmolar();
1712 IO.smolar_liq = HEOS.SatL.get()->smolar();
1713 IO.smolar_vap = HEOS.SatV.get()->smolar();
1714}
1715
1716void SaturationSolvers::newton_raphson_twophase::build_arrays() {
1717 // References to the classes for concision
1718 HelmholtzEOSMixtureBackend &rSatL = *(HEOS->SatL.get()), &rSatV = *(HEOS->SatV.get());
1719
1720 // Zero the Jacobian first: the beta-constraint rows (k = N .. 2N-2) below only write the
1721 // two mole-fraction columns the constraint depends on; the imposed-variable (T/p) column,
1722 // and for N>=3 the off-diagonal mole-fraction columns, are left untouched and their true
1723 // value is 0. Eigen::resize() does not zero storage, so without this those entries are
1724 // uninitialized memory feeding the Newton solve (GH #3192 follow-up).
1725 J.setZero();
1726
1727 // Step 0:
1728 // -------
1729 // Set mole fractions
1730 rSatL.set_mole_fractions(x);
1731 rSatV.set_mole_fractions(y);
1732
1733 //std::vector<CoolPropDbl> &x = rSatL.get_mole_fractions();
1734 //std::vector<CoolPropDbl> &y = rSatV.get_mole_fractions();
1735
1736 rSatL.update_TP_guessrho(T, p, rhomolar_liq);
1737 rhomolar_liq = rSatL.rhomolar();
1738 rSatV.update_TP_guessrho(T, p, rhomolar_vap);
1739 rhomolar_vap = rSatV.rhomolar();
1740
1741 // For diagnostic purposes calculate the pressures (no derivatives are evaluated)
1742 CoolPropDbl p_liq = rSatL.p();
1743 CoolPropDbl p_vap = rSatV.p();
1744 p = 0.5 * (p_liq + p_vap);
1745
1746 // Step 2:
1747 // -------
1748 // Build the residual vector and the Jacobian matrix
1749
1750 x_N_dependency_flag xN_flag = XN_DEPENDENT;
1751
1752 // Form of residuals do not depend on which variable is imposed
1753 for (std::size_t i = 0; i < N; ++i) {
1754 // Equate the liquid and vapor fugacities
1755 CoolPropDbl ln_f_liq = log(MixtureDerivatives::fugacity_i(rSatL, i, xN_flag));
1756 CoolPropDbl ln_f_vap = log(MixtureDerivatives::fugacity_i(rSatV, i, xN_flag));
1757 r[i] = ln_f_liq - ln_f_vap; // N of these
1758
1759 if (i != N - 1) {
1760 // Equate the specified vapor mole fraction and that given defined by the ith component
1761 r[i + N] = (z[i] - x[i]) / (y[i] - x[i]) - beta; // N-1 of these
1762 }
1763 }
1764
1765 // First part of derivatives with respect to ln f_i
1766 for (std::size_t i = 0; i < N; ++i) {
1767 for (std::size_t j = 0; j < N - 1; ++j) {
1768 J(i, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatL, i, j, xN_flag);
1769 J(i, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(rSatV, i, j, xN_flag);
1770 }
1771
1772 // Last derivative with respect to either T or p depending on what is imposed
1773 if (imposed_variable == newton_raphson_twophase_options::P_IMPOSED) {
1774 J(i, 2 * N - 2) =
1775 MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dT__constp_n(rSatV, i, xN_flag);
1776 } else if (imposed_variable == newton_raphson_twophase_options::T_IMPOSED) {
1777 J(i, 2 * N - 2) =
1778 MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatL, i, xN_flag) - MixtureDerivatives::dln_fugacity_i_dp__constT_n(rSatV, i, xN_flag);
1779 } else {
1780 throw ValueError();
1781 }
1782 }
1783 // Derivatives with respect to the vapor mole fractions residual
1784 for (std::size_t i = 0; i < N - 1; ++i) {
1785 std::size_t k = i + N; // N ln f_i residuals
1786 J(k, i) = (z[i] - y[i]) / pow(y[i] - x[i], 2);
1787 J(k, i + (N - 1)) = -(z[i] - x[i]) / pow(y[i] - x[i], 2);
1788 }
1789
1790 error_rms = r.norm(); // Square-root (The R in RMS)
1791}
1792
1793class RachfordRiceResidual : public FuncWrapper1DWithDeriv
1794{
1795 private:
1796 const std::vector<double>&z, &lnK;
1797
1798 public:
1799 RachfordRiceResidual(const std::vector<double>& z, const std::vector<double>& lnK) : z(z), lnK(lnK) {};
1800 double call(double beta) override {
1801 return FlashRoutines::g_RachfordRice(z, lnK, beta);
1802 }
1803 double deriv(double beta) override {
1804 return FlashRoutines::dgdbeta_RachfordRice(z, lnK, beta);
1805 }
1806};
1807
1808// Solve the Rachford-Rice residual sum_i z_i (K_i - 1) / (1 - beta + beta K_i) = 0
1809// for beta on [0, 1] by bisection. The caller ensures the residual does not change
1810// sign in the wrong direction (residual >= 0 at beta = 0 and <= 0 at beta = 1); the
1811// residual is strictly decreasing with no interior pole (1 + beta (K_i - 1) > 0 for
1812// K_i > 0 and beta in [0, 1]), so bisection converges to the (possibly boundary) root.
1813// Used in place of the double-typed RachfordRiceResidual + Brent so the
1814// CoolPropDbl (incl. long-double) build is type-consistent and no solver throws here.
1815static CoolPropDbl rachford_rice_beta_bisect(const std::vector<CoolPropDbl>& z, const std::vector<CoolPropDbl>& K) {
1816 CoolPropDbl a = 0.0, b = 1.0;
1817 for (int it = 0; it < 100; ++it) {
1818 CoolPropDbl beta = 0.5 * (a + b);
1819 CoolPropDbl g = 0;
1820 for (std::size_t i = 0; i < z.size(); ++i)
1821 g += z[i] * (K[i] - 1.0) / (1.0 - beta + beta * K[i]);
1822 if (g > 0)
1823 a = beta; // residual is decreasing: g > 0 means beta is too small
1824 else
1825 b = beta;
1826 }
1827 return 0.5 * (a + b);
1828}
1829
1830void SaturationSolvers::successive_substitution_guessrho(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y,
1831 CoolPropDbl& rhomolar_liq, CoolPropDbl& rhomolar_vap, const std::vector<CoolPropDbl>& z,
1832 int num_steps, double tol) {
1833 // Successive-substitution refinement of a two-phase guess at fixed (T, p).
1834 // Adapted from jakobreichert's PR #2720: re-solve each phase density near its
1835 // current guess, recompute K from the fugacity-coefficient ratio, and re-split
1836 // via Rachford-Rice.
1837 const std::size_t N = z.size();
1838 std::vector<CoolPropDbl> lnK(N, 0.0), K(N);
1839 for (int ss = 0; ss < num_steps; ++ss) {
1840 // Preserve the current (last-good) densities so a bad density root this step can be
1841 // rolled back: update_TP_guessrho overwrites rhomolar_liq/vap before the finite check
1842 // below, and on !finite the caller's Newton solver must still receive a valid seed.
1843 const CoolPropDbl rho_liq_prev = rhomolar_liq, rho_vap_prev = rhomolar_vap;
1844 HEOS.SatL->set_mole_fractions(x);
1845 HEOS.SatL->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_liq);
1846 HEOS.SatV->set_mole_fractions(y);
1847 HEOS.SatV->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_vap);
1848 rhomolar_liq = HEOS.SatL->rhomolar();
1849 rhomolar_vap = HEOS.SatV->rhomolar();
1850
1851 CoolPropDbl g0 = 0, g1 = 0, max_lnK_change = 0;
1852 bool finite = true;
1853 for (std::size_t i = 0; i < N; ++i) {
1854 CoolPropDbl lnK_new = log(HEOS.SatL->fugacity_coefficient(i) / HEOS.SatV->fugacity_coefficient(i));
1855 if (!ValidNumber(lnK_new)) {
1856 finite = false;
1857 break;
1858 }
1859 max_lnK_change = std::max(max_lnK_change, std::abs(lnK_new - lnK[i]));
1860 lnK[i] = lnK_new;
1861 K[i] = exp(lnK[i]);
1862 if (!ValidNumber(K[i])) { // exp overflow at an extreme lnK -> step is unusable
1863 finite = false;
1864 break;
1865 }
1866 g0 += z[i] * (K[i] - 1.0);
1867 g1 += z[i] * (1.0 - 1.0 / K[i]);
1868 }
1869 // A non-finite fugacity coefficient (e.g. a bad density root) makes this SS step
1870 // meaningless; restore the last-good densities (x/y are not yet updated this step) and
1871 // stop, so the Newton solver receives the last valid x/y/rho.
1872 if (!finite) {
1873 rhomolar_liq = rho_liq_prev;
1874 rhomolar_vap = rho_vap_prev;
1875 break;
1876 }
1877
1878 CoolPropDbl beta;
1879 if (g0 < 0)
1880 beta = 0;
1881 else if (g1 > 0)
1882 beta = 1;
1883 else
1884 beta = rachford_rice_beta_bisect(z, K);
1885 x_and_y_from_K(beta, K, z, x, y);
1886 normalize_vector(x);
1887 normalize_vector(y);
1888
1889 if (ss > 0 && max_lnK_change < tol) break;
1890 }
1891}
1892
1893bool SaturationSolvers::guess_split_from_wilson(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y,
1894 CoolPropDbl& rhomolar_liq, CoolPropDbl& rhomolar_vap, const std::vector<CoolPropDbl>& z,
1895 CoolPropDbl T, CoolPropDbl p, int num_steps, bool require_bracket) {
1896 // Seed a two-phase guess from the ideal (Wilson) K-factor estimate and refine it
1897 // by successive substitution. Returns false when the Wilson estimate does not
1898 // bracket a two-phase Rachford-Rice root (the feed is outside its bubble/dew
1899 // points) or when a phase density cannot be obtained. May also throw from the
1900 // underlying density solver; the blind-flash caller treats any throw as "not
1901 // two-phase" and falls back to the single-phase path. Used to recover a genuinely
1902 // two-phase state that the TPD stability test reported as single-phase, e.g. for
1903 // cubic mixtures at high vapor fraction near the dew point (CoolProp-zgpy).
1904 const std::size_t N = z.size();
1905 std::vector<CoolPropDbl> K(N), lnK(N);
1906 CoolPropDbl g0 = 0, g1 = 0;
1907 for (std::size_t i = 0; i < N; ++i) {
1908 lnK[i] = Wilson_lnK_factor(HEOS, T, p, i);
1909 if (!ValidNumber(lnK[i])) return false;
1910 K[i] = exp(lnK[i]);
1911 if (!ValidNumber(K[i])) return false; // exp overflow at an extreme lnK -> no usable estimate
1912 g0 += z[i] * (K[i] - 1.0); // Rachford-Rice residual at beta = 0
1913 g1 += z[i] * (1.0 - 1.0 / K[i]); // Rachford-Rice residual at beta = 1
1914 }
1915 // Two-phase iff the residual changes sign on (0, 1): g0 > 0 (z above its bubble
1916 // point) AND g1 < 0 (z below its dew point). When require_bracket is false the
1917 // caller (a non-conclusive stability verdict) forces an attempt even though the
1918 // IDEAL Wilson estimate places the feed outside (0, 1); the EOS-based SS refinement
1919 // below then decides, and the flash's verify/fallback rejects it if it does not
1920 // converge to a genuine equilibrium split.
1921 bool bracketed = (g0 > 0 && g1 < 0);
1922 if (require_bracket && !bracketed) return false;
1923
1924 CoolPropDbl beta = bracketed ? rachford_rice_beta_bisect(z, K) : 0.5;
1925 x.resize(N);
1926 y.resize(N);
1927 x_and_y_from_K(beta, K, z, x, y);
1928 normalize_vector(x);
1929 normalize_vector(y);
1930
1931 // Density guesses at the (heavy-rich) liquid and (light-rich) vapor compositions --
1932 // NOT the feed, whose single (vapor) root at high vapor fraction would collapse the
1933 // split to the trivial solution.
1934 HEOS.SatL->set_mole_fractions(x);
1935 rhomolar_liq = HEOS.SatL->solver_rho_Tp_global(T, p, HEOS.SatL->calc_rhomolar_max_bound());
1936 HEOS.SatV->set_mole_fractions(y);
1937 rhomolar_vap = HEOS.SatV->solver_rho_Tp_global(T, p, HEOS.SatV->calc_rhomolar_max_bound());
1938 if (!ValidNumber(rhomolar_liq) || rhomolar_liq <= 0 || !ValidNumber(rhomolar_vap) || rhomolar_vap <= 0) return false;
1939
1940 successive_substitution_guessrho(HEOS, x, y, rhomolar_liq, rhomolar_vap, z, num_steps);
1941 return true;
1942}
1943
1944void StabilityRoutines::StabilityEvaluationClass::trial_compositions() {
1945
1946 x.resize(z.size());
1947 y.resize(z.size());
1948 lnK.resize(z.size());
1949 K.resize(z.size());
1950 double g0 = 0, g1 = 0, beta = -1;
1951
1952 for (int i = 0; i < static_cast<int>(z.size()); ++i) {
1953 // Calculate the K-factor
1954 if (m_T < 0 && m_p < 0) {
1955 // Using T&P from the class
1956 lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, HEOS.T(), HEOS.p(), i);
1957 } else {
1958 // Using specified T&P
1959 lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, m_T, m_p, i);
1960 }
1961 K[i] = exp(lnK[i]);
1962 g0 += z[i] * (K[i] - 1); // The summation for beta = 0
1963 g1 += z[i] * (1 - 1 / K[i]); // The summation for beta = 1
1964 }
1965 // Copy K-factors for later use
1966 K0 = K;
1967 // Now see what to do about the g(0) and g(1) values
1968 // -----
1969 //
1970 if (g0 < 0) {
1971 beta = 0; // Assumed to be at bubble-point temperature
1972 } else if (g1 > 0) {
1973 beta = 1; // Assumed to be at the dew-point temperature
1974 } else {
1975 // Need to iterate to find beta that makes g of Rachford-Rice zero
1976 RachfordRiceResidual resid(z, lnK);
1977 beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
1978 }
1979 // Get the compositions from given value for beta, K, z
1980 SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
1981 normalize_vector(x);
1982 normalize_vector(y);
1983 if (debug) {
1984 std::cout << format("1) T: %g p: %g beta: %g\n", HEOS.T(), HEOS.p(), beta);
1985 }
1986}
1987void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int num_steps) {
1988 // ----
1989 // Do a few steps of successive substitution
1990 // ----
1991
1992 HEOS.SatL->set_mole_fractions(x);
1993 HEOS.SatL->calc_reducing_state();
1994 HEOS.SatV->set_mole_fractions(y);
1995 HEOS.SatV->calc_reducing_state();
1996
1997 if (debug) {
1998 std::cout << format("2) SS1: i beta K x y rho' rho''\n");
1999 }
2000 for (int step_count = 0; step_count < num_steps; ++step_count) {
2001 // Set the composition
2002 HEOS.SatL->set_mole_fractions(x);
2003 HEOS.SatV->set_mole_fractions(y);
2004 HEOS.SatL->calc_reducing_state();
2005 HEOS.SatV->calc_reducing_state();
2006
2007 this->rho_TP_global();
2008
2009 // Calculate the new K-factors from the fugacity coefficients
2010 double g0 = 0, g1 = 0;
2011 for (std::size_t i = 0; i < z.size(); ++i) {
2012 lnK[i] = log(HEOS.SatL->fugacity_coefficient(i) / HEOS.SatV->fugacity_coefficient(i));
2013 K[i] = exp(lnK[i]);
2014 g0 += z[i] * (K[i] - 1); // The summation for beta = 0
2015 g1 += z[i] * (1 - 1 / K[i]); // The summation for beta = 1
2016 }
2017 RachfordRiceResidual resid(z, lnK);
2018 if (g0 < 0) {
2019 beta = 0;
2020 } else if (g1 > 0) {
2021 beta = 1;
2022 } else {
2023 // Need to iterate to find beta that makes g of Rachford-Rice zero
2024 beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
2025 }
2026
2027 // Get the compositions from given values for beta, K, z
2028 SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
2029 normalize_vector(x);
2030 normalize_vector(y);
2031 if (debug) {
2032 std::cout << format("2) %d %g %s %s %s %g %g\n", step_count, beta, vec_to_string(K, "%0.6f").c_str(), vec_to_string(x, "%0.6f").c_str(),
2033 vec_to_string(y, "%0.6f").c_str(), rhomolar_liq, rhomolar_vap);
2034 }
2035 }
2036}
2047void StabilityRoutines::StabilityEvaluationClass::check_stability() {
2048 if (use_michelsen) {
2049 check_stability_michelsen();
2050 } else {
2051 check_stability_legacy();
2052 }
2053}
2054
2055// Solve a trial phase's density at fixed (T, p) for its current composition, WARM-STARTED
2056// from the previous root (rho_warm). The stability TPD trajectory moves the composition
2057// gradually, so the stable density root tracks continuously; a local Newton from the prior
2058// root (update_TP_guessrho) reaches it in a few EOS evaluations, vs the full spinodal scan
2059// + double Brent solve of solver_rho_Tp_global. Falls back to the global solver on the
2060// first call (rho_warm <= 0), or if the warm solve fails or returns a non-physical root,
2061// so the stable-root contract is preserved and the throw-on-total-failure contract (which
2062// callers rely on) is unchanged. Updates the backend state and rho_warm to the new root.
2063static CoolPropDbl solve_trial_rho_warm(HelmholtzEOSMixtureBackend& phase, CoolPropDbl T, CoolPropDbl p, CoolPropDbl& rho_warm) {
2064 if (rho_warm > 0) {
2065 CoolPropDbl r = -1;
2066 bool warm_ok = false;
2067 try {
2068 phase.update_TP_guessrho(T, p, rho_warm);
2069 r = phase.rhomolar();
2070 // Accept the warm (local) root only if it is finite, positive, and has NOT
2071 // jumped branches. The composition moves gradually along these trajectories,
2072 // so the stable root's density changes only slightly per step; a large jump
2073 // (the liquid and vapor branches differ by ~10x or more) signals the local
2074 // Newton landing on the OTHER, now-metastable branch after a spinodal crossing.
2075 // (Near the critical point the two branches merge, so a sub-2x change there is
2076 // genuinely the same root.)
2077 warm_ok = ValidNumber(r) && r > 0 && r < 2.0 * rho_warm && r > 0.5 * rho_warm;
2078 } catch (...) {
2079 warm_ok = false; // warm solve threw -> fall back to the global solver below
2080 }
2081 if (warm_ok) {
2082 rho_warm = r;
2083 return r;
2084 }
2085 }
2086 // Cold start, branch jump, or warm-solve failure: re-confirm with the global,
2087 // lowest-Gibbs solver so a metastable root is never silently accepted.
2088 CoolPropDbl rg = phase.solver_rho_Tp_global(T, p, phase.calc_rhomolar_max_bound());
2089 phase.update_DmolarT_direct(rg, T);
2090 rho_warm = rg;
2091 return rg;
2092}
2093
2094void StabilityRoutines::StabilityEvaluationClass::check_stability_michelsen() {
2095 const std::size_t N = z.size();
2096 CoolPropDbl the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
2097 CoolPropDbl the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
2098 _stable = true;
2099 _uncertain = false;
2100 bool any_uncertain = false; // a trial's minimize_tpd was non-conclusive (step/density fail, max-iter)
2101
2102 // Evaluate feed fugacities: d_i = ln(z_i) + ln(phi_i(z))
2103 HEOS.SatL->set_mole_fractions(z);
2104 CoolPropDbl rho_b;
2105 try {
2106 rho_b = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, HEOS.SatL->calc_rhomolar_max_bound());
2107 } catch (...) {
2108 // solver_rho_Tp_global can fail for multiparameter mixtures when the pressure
2109 // lies between the spinodal pressures. Fall back to SRK-seeded solver.
2110 HEOS.SatL->specify_phase(iphase_gas);
2111 rho_b = HEOS.SatL->solver_rho_Tp(the_T, the_p);
2112 HEOS.SatL->unspecify_phase();
2113 }
2114 HEOS.SatL->update_DmolarT_direct(rho_b, the_T);
2115
2116 std::vector<CoolPropDbl> ln_f_z(N);
2117 for (std::size_t i = 0; i < N; ++i) {
2118 if (z[i] > 0)
2119 ln_f_z[i] = std::log(z[i]) + std::log(HEOS.SatL->fugacity_coefficient(i));
2120 else
2121 ln_f_z[i] = -1e30; // Effectively -inf for absent components
2122 }
2123
2124 // Build two trial compositions from Wilson K-factors ([Michelsen1982a] Eq. 28):
2125 // K_i = (Pc_i/P) * exp(5.373*(1+omega_i)*(1-Tc_i/T))
2126 std::vector<CoolPropDbl> yV(N), xL(N);
2127 for (std::size_t i = 0; i < N; ++i) {
2128 double Ki = std::exp(SaturationSolvers::Wilson_lnK_factor(HEOS, the_T, the_p, i));
2129 yV[i] = z[i] * Ki;
2130 xL[i] = z[i] / Ki;
2131 }
2132
2133 // Test both trial directions (vapor-like and liquid-like)
2134 std::vector<std::vector<CoolPropDbl>> trials = {yV, xL};
2135 for (std::size_t t = 0; t < trials.size(); ++t) {
2136 auto& Y = trials[t];
2137 // Warm-start density root for this trial's composition trajectory. Reset per trial:
2138 // the vapor-like and liquid-like trials live on different density branches.
2139 CoolPropDbl rho_warm = -1;
2140
2141 // --- Phase 1: Successive substitution with GDEM acceleration ---
2142 // Fixed-point map: ln(Y_i^new) = d_i - ln(phi_i(y_norm))
2143 // See [Michelsen1982a] Eq. 19 and [M&M2007] Ch. 12, Sec. 12.6
2144 const int max_ss_loops = 4; // Each loop does 2 SS steps + 1 GDEM step
2145 const double cntol = 1e-7;
2146 bool ss_decided = false;
2147
2148 for (int loop = 0; loop < max_ss_loops && !ss_decided; ++loop) {
2149 std::array<double, 2> esq_pair = {0, 0};
2150 std::vector<CoolPropDbl> err(N);
2151
2152 for (int kk = 0; kk < 2 && !ss_decided; ++kk) {
2153 // Normalize Y to get trial composition
2154 CoolPropDbl sumY = 0;
2155 for (std::size_t i = 0; i < N; ++i)
2156 sumY += Y[i];
2157 std::vector<CoolPropDbl> y_norm(N);
2158 for (std::size_t i = 0; i < N; ++i)
2159 y_norm[i] = Y[i] / sumY;
2160
2161 // Evaluate fugacity coefficients at trial composition
2162 HEOS.SatV->set_mole_fractions(y_norm);
2163 try {
2164 solve_trial_rho_warm(*HEOS.SatV, the_T, the_p, rho_warm);
2165 } catch (...) {
2166 ss_decided = true;
2167 break;
2168 }
2169
2170 // Compute new Y, objective function, and convergence metrics
2171 CoolPropDbl tm = 1.0; // Modified TPD: tm = 1 + sum Y_i*(ln Y_i + ln phi_i - d_i - 1)
2172 CoolPropDbl gmax = 0;
2173 double esq = 0;
2174 for (std::size_t i = 0; i < N; ++i) {
2175 double ln_phi_y = std::log(HEOS.SatV->fugacity_coefficient(i));
2176 double ln_Y_new = ln_f_z[i] - ln_phi_y;
2177 double ln_Y_old = std::log(std::max(Y[i], 1e-300));
2178 double diff = ln_Y_new - ln_Y_old;
2179 err[i] = diff;
2180 esq += Y[i] * diff * diff;
2181 gmax = std::max(gmax, std::abs(diff));
2182
2183 double s_i = ln_Y_old + ln_phi_y - ln_f_z[i];
2184 tm += Y[i] * (s_i - 1.0);
2185
2186 Y[i] = std::exp(ln_Y_new);
2187 }
2188
2189 // Early exit: tm < 0 means unstable
2190 if (tm < -cntol) {
2191 _stable = false;
2192 CoolPropDbl sY = 0;
2193 for (std::size_t i = 0; i < N; ++i)
2194 sY += Y[i];
2195 for (std::size_t i = 0; i < N; ++i)
2196 y_norm[i] = Y[i] / sY;
2197 if (t == 0) {
2198 this->y = y_norm;
2199 this->x = z;
2200 } else {
2201 this->x = y_norm;
2202 this->y = z;
2203 }
2204 return;
2205 }
2206
2207 // Converged to a stationary point
2208 if (gmax < cntol) {
2209 ss_decided = true;
2210 break;
2211 }
2212
2213 // Proximity test: if trial is close to feed and curvature is positive, stable
2214 CoolPropDbl distance_sq = 0, curvature = 0;
2215 for (std::size_t i = 0; i < N; ++i) {
2216 double zysq = std::sqrt(Y[i] * z[i]);
2217 distance_sq += Y[i] + z[i] - 2.0 * zysq;
2218 curvature -= (Y[i] - zysq) * err[i];
2219 }
2220 if (distance_sq < 0) distance_sq = -distance_sq;
2221 if (std::sqrt(distance_sq) < 0.1 && curvature > 0 && tm / curvature > 0.8) {
2222 ss_decided = true;
2223 break;
2224 }
2225
2226 esq_pair[kk] = esq;
2227 }
2228
2229 if (ss_decided) break;
2230
2231 // GDEM extrapolation step ([M&M2007] Ch. 12, Sec. 12.6)
2232 if (esq_pair[0] > 0) {
2233 double ratio = std::sqrt(esq_pair[1] / esq_pair[0]);
2234 // A NaN ratio (non-finite esq) passes both < 0 and >= 0.95, so guard it
2235 // explicitly -- otherwise it poisons the GDEM lnK/Y update (CoolProp-1tbe.8 finding 4c).
2236 if (!ValidNumber(ratio) || ratio < 0 || ratio >= 0.95) ratio = 0.95;
2237 double factor = ratio / (1.0 - ratio);
2238 for (std::size_t i = 0; i < N; ++i) {
2239 double ln_Y = std::log(std::max(Y[i], 1e-300));
2240 ln_Y += factor * err[i];
2241 Y[i] = std::exp(ln_Y);
2242 }
2243 }
2244 } // end SS+GDEM loops
2245
2246 // --- Phase 2: Second-order TPD minimization ---
2247 // If SS did not conclusively decide stability, use quasi-Newton
2248 // minimization in alpha variables. See [Michelsen1982a] Eq. 25-27.
2249 bool trial_unstable = false;
2250 bool trial_ok = minimize_tpd(Y, ln_f_z, the_T, the_p, trial_unstable);
2251 if (!trial_ok) any_uncertain = true; // could not conclude this trial -> not a clean "stable"
2252 if (trial_ok) {
2253 if (trial_unstable) {
2254 _stable = false;
2255 CoolPropDbl sY = 0;
2256 for (std::size_t i = 0; i < N; ++i)
2257 sY += Y[i];
2258 std::vector<CoolPropDbl> y_norm(N);
2259 for (std::size_t i = 0; i < N; ++i)
2260 y_norm[i] = Y[i] / sY;
2261 if (t == 0) {
2262 this->y = y_norm;
2263 this->x = z;
2264 } else {
2265 this->x = y_norm;
2266 this->y = z;
2267 }
2268 return;
2269 }
2270 }
2271 // If minimize_tpd fails or finds tm >= 0, try the other trial direction
2272 }
2273 // Both trials indicate stability -- but flag a non-conclusive verdict so the caller
2274 // attempts a verified split instead of trusting a fail-open (CoolProp-zgpy).
2275 _uncertain = any_uncertain;
2276 _stable = true;
2277}
2278
2288bool StabilityRoutines::StabilityEvaluationClass::minimize_tpd(std::vector<CoolPropDbl>& Y, const std::vector<CoolPropDbl>& ln_f_z, CoolPropDbl the_T,
2289 CoolPropDbl the_p, bool& is_unstable) {
2290
2291 const std::size_t N = Y.size();
2292 const double cntol = 1e-7;
2293 const int max_iter = 20;
2294 is_unstable = false;
2295
2296 // Alpha variables: alpha_i = 2*sqrt(Y_i)
2297 std::vector<CoolPropDbl> alpha(N), alpha_old(N);
2298 for (std::size_t i = 0; i < N; ++i)
2299 alpha[i] = 2.0 * std::sqrt(std::max(Y[i], 1e-300));
2300
2301 double trust_radius = 0.25; // Initial trust-region size
2302 double diagonal_shift = 0.0;
2303 CoolPropDbl rho_warm = -1; // warm-start density root, tracked across iterations/inner steps
2304
2305 for (int iter = 0; iter < max_iter; ++iter) {
2306 // Update Y from alpha
2307 for (std::size_t i = 0; i < N; ++i) {
2308 Y[i] = (0.5 * alpha[i]) * (0.5 * alpha[i]);
2309 alpha_old[i] = alpha[i];
2310 }
2311
2312 // Evaluate fugacity coefficients at normalized trial composition
2313 CoolPropDbl sumY = 0;
2314 for (std::size_t i = 0; i < N; ++i)
2315 sumY += Y[i];
2316 std::vector<CoolPropDbl> y_norm(N);
2317 for (std::size_t i = 0; i < N; ++i)
2318 y_norm[i] = Y[i] / sumY;
2319
2320 HEOS.SatV->set_mole_fractions(y_norm);
2321 try {
2322 solve_trial_rho_warm(*HEOS.SatV, the_T, the_p, rho_warm);
2323 } catch (...) {
2324 return false; // Density solve failed
2325 }
2326
2327 // Compute gradient and Hessian in alpha variables
2328 // s_i = ln(Y_i) + ln(phi_i(y)) - d_i (scaled gradient component)
2329 // grad_i = -alpha_i/2 * s_i
2330 // H_ij = delta_ij * (1 + s_i/2) + (alpha_i*alpha_j)/(4*N_tot) * d(ln phi_i)/dn_j
2331 std::vector<CoolPropDbl> scaled_grad(N), grad(N), half_alpha(N);
2332 double max_gradient = 0, obj_old = 1.0;
2333
2334 for (std::size_t i = 0; i < N; ++i) {
2335 half_alpha[i] = alpha[i] * 0.5;
2336 if (z[i] > 0) {
2337 double ln_Y = std::log(std::max(Y[i], 1e-300));
2338 double ln_phi = std::log(HEOS.SatV->fugacity_coefficient(i));
2339 scaled_grad[i] = ln_Y + ln_phi - ln_f_z[i];
2340 } else {
2341 scaled_grad[i] = 0;
2342 }
2343 grad[i] = -scaled_grad[i] * half_alpha[i];
2344 max_gradient = std::max(max_gradient, std::abs(grad[i]));
2345 obj_old += Y[i] * (scaled_grad[i] - 1.0);
2346 }
2347
2348 // Converged?
2349 if (max_gradient < cntol) {
2350 is_unstable = (obj_old < -cntol);
2351 return true;
2352 }
2353
2354 // Build Hessian
2355 Eigen::MatrixXd H(N, N);
2356 for (std::size_t i = 0; i < N; ++i) {
2357 double ahi = half_alpha[i] / sumY;
2358 for (std::size_t j = 0; j <= i; ++j) {
2359 double dln_phi_dnj = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), i, j, XN_INDEPENDENT);
2360 double term = ahi * half_alpha[j] * dln_phi_dnj;
2361 H(i, j) = term;
2362 H(j, i) = term;
2363 }
2364 H(i, i) += 1.0 + 0.5 * scaled_grad[i];
2365 }
2366
2367 // Solve (H + lambda*I) * delta_alpha = -grad using trust-region method
2368 // Simplified Hebden: try full Newton step, shrink trust region if objective increases
2369 Eigen::VectorXd g_vec(N), delta_alpha(N);
2370 for (std::size_t i = 0; i < N; ++i)
2371 g_vec(i) = grad[i];
2372
2373 // Inner loop: adjust diagonal_shift until step is accepted
2374 const int max_inner = 20;
2375 bool step_accepted = false;
2376 for (int inner = 0; inner < max_inner; ++inner) {
2377 Eigen::MatrixXd H_shifted = H;
2378 H_shifted.diagonal().array() += diagonal_shift;
2379
2380 delta_alpha = H_shifted.colPivHouseholderQr().solve(g_vec);
2381
2382 // Clamp to prevent negative alpha
2383 double step_norm_sq = 0;
2384 for (std::size_t i = 0; i < N; ++i) {
2385 double da = delta_alpha(i);
2386 if (da + alpha_old[i] <= 0) da = -0.9 * alpha_old[i];
2387 delta_alpha(i) = da;
2388 alpha[i] = alpha_old[i] + da;
2389 Y[i] = (0.5 * alpha[i]) * (0.5 * alpha[i]);
2390 step_norm_sq += da * da;
2391 }
2392 double step_size = std::sqrt(step_norm_sq);
2393
2394 // Check if step is within trust region
2395 if (step_size > trust_radius && diagonal_shift == 0) {
2396 // Step too large, add regularization
2397 diagonal_shift = step_size / trust_radius - 1.0;
2398 // Restore alpha
2399 for (std::size_t i = 0; i < N; ++i) {
2400 alpha[i] = alpha_old[i];
2401 Y[i] = (0.5 * alpha[i]) * (0.5 * alpha[i]);
2402 }
2403 continue;
2404 }
2405
2406 // Evaluate new objective
2407 sumY = 0;
2408 for (std::size_t i = 0; i < N; ++i)
2409 sumY += Y[i];
2410 for (std::size_t i = 0; i < N; ++i)
2411 y_norm[i] = Y[i] / sumY;
2412
2413 HEOS.SatV->set_mole_fractions(y_norm);
2414 try {
2415 solve_trial_rho_warm(*HEOS.SatV, the_T, the_p, rho_warm);
2416 } catch (...) {
2417 // Density solve failed, shrink trust region
2418 trust_radius = step_size / 3.0;
2419 diagonal_shift = 0;
2420 for (std::size_t i = 0; i < N; ++i) {
2421 alpha[i] = alpha_old[i];
2422 Y[i] = (0.5 * alpha[i]) * (0.5 * alpha[i]);
2423 }
2424 continue;
2425 }
2426
2427 double obj_new = 1.0;
2428 for (std::size_t i = 0; i < N; ++i) {
2429 double ln_Y = std::log(std::max(Y[i], 1e-300));
2430 double ln_phi = std::log(HEOS.SatV->fugacity_coefficient(i));
2431 obj_new += Y[i] * (ln_Y + ln_phi - ln_f_z[i] - 1.0);
2432 }
2433
2434 // Quick exit if already found unstable
2435 if (obj_new < -cntol) {
2436 is_unstable = true;
2437 return true;
2438 }
2439
2440 // Check actual vs predicted reduction
2441 if (obj_new > obj_old + 1e-12) {
2442 // Objective increased, shrink trust region
2443 trust_radius = step_size / 3.0;
2444 diagonal_shift = 0;
2445 for (std::size_t i = 0; i < N; ++i) {
2446 alpha[i] = alpha_old[i];
2447 Y[i] = (0.5 * alpha[i]) * (0.5 * alpha[i]);
2448 }
2449 continue;
2450 }
2451
2452 // Compute predicted reduction for trust-region update
2453 Eigen::VectorXd Hd = H * delta_alpha;
2454 double predicted = 0.5 * delta_alpha.dot(Hd) - delta_alpha.dot(g_vec);
2455 double actual = obj_new - obj_old;
2456 double ratio = (predicted != 0) ? actual / predicted : 1.0;
2457
2458 if (ratio < 0.25) {
2459 trust_radius = step_size / 2.0;
2460 } else if (ratio > 0.75 && diagonal_shift > 0) {
2461 trust_radius = step_size * 2.0;
2462 } else {
2463 trust_radius = step_size;
2464 }
2465 diagonal_shift = 0;
2466 step_accepted = true;
2467 break;
2468 }
2469
2470 if (!step_accepted) {
2471 // Could not find an acceptable step
2472 return false;
2473 }
2474 }
2475 return false; // Reached max iterations without convergence -> non-conclusive (caller decides)
2476}
2477
2478void StabilityRoutines::StabilityEvaluationClass::check_stability_legacy() {
2479 std::vector<double> tpdL, tpdH;
2480
2481 // Calculate the temperature and pressure to be used
2482 double the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
2483 double the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
2484
2485 // If beta value is between epsilon and 1-epsilon, check the TPD
2486 if (beta > DBL_EPSILON && beta < 1 - DBL_EPSILON) {
2487
2488 // Set the composition back to the bulk composition for both liquid and vapor phases
2489 HEOS.SatL->set_mole_fractions(z);
2490 HEOS.SatV->set_mole_fractions(z);
2491 HEOS.SatL->calc_reducing_state();
2492 HEOS.SatV->calc_reducing_state();
2493
2494 // Update the densities in each class
2495 double rhoL = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, HEOS.SatL->calc_rhomolar_max_bound());
2496 double rhoV = HEOS.SatV->solver_rho_Tp_global(the_T, the_p, HEOS.SatV->calc_rhomolar_max_bound());
2497 HEOS.SatL->update_DmolarT_direct(rhoL, the_T);
2498 HEOS.SatV->update_DmolarT_direct(rhoV, the_T);
2499
2500 // Calculate the tpd and the Gibbs energy difference (Gernert, 2014, Eqs. 20-22)
2501 // The trial compositions are the phase compositions from before
2502 this->tpd_liq = HEOS.SatL->tangent_plane_distance(the_T, the_p, x, rhomolar_liq);
2503 this->tpd_vap = HEOS.SatV->tangent_plane_distance(the_T, the_p, y, rhomolar_vap);
2504
2505 this->DELTAG_nRT = (1 - beta) * tpd_liq + beta * (tpd_vap);
2506 if (debug) {
2507 std::cout << format("3) tpd': %g tpd'': %g DELTAG/nRT: %g\n", tpd_liq, tpd_vap, DELTAG_nRT);
2508 }
2509
2510 // If any of these cases are met, feed is conclusively unstable, stop!
2511 if (this->tpd_liq < -DBL_EPSILON || this->tpd_vap < -DBL_EPSILON || this->DELTAG_nRT < -DBL_EPSILON) {
2512 if (debug) {
2513 std::cout << format("3) PHASE SPLIT beta in (eps,1-eps) \n");
2514 }
2515 _stable = false;
2516 return;
2517 }
2518 }
2519
2520 // Ok, we aren't sure about stability, need to keep going with the full tpd analysis
2521
2522 // Use the global density solver to obtain the density root (or the lowest Gibbs energy root if more than one)
2523 CoolPropDbl rho_bulk = HEOS.solver_rho_Tp_global(the_T, the_p, HEOS.calc_rhomolar_max_bound());
2524 HEOS.update_DmolarT_direct(rho_bulk, the_T);
2525
2526 // Calculate the fugacity coefficient at initial composition of the bulk phase.
2527 std::vector<double> fugacity_coefficient0(z.size()), fugacity0(z.size());
2528 for (std::size_t i = 0; i < z.size(); ++i) {
2529 fugacity_coefficient0[i] = exp(MixtureDerivatives::ln_fugacity_coefficient(HEOS, i, XN_DEPENDENT));
2530 fugacity0[i] = MixtureDerivatives::fugacity_i(HEOS, i, XN_DEPENDENT);
2531 }
2532
2533 // Generate light and heavy test compositions (Gernert, 2014, Eq. 23)
2534 xL.resize(z.size());
2535 xH.resize(z.size());
2536 for (std::size_t i = 0; i < z.size(); ++i) {
2537 xL[i] = z[i] * K0[i]; // Light-phase composition
2538 xH[i] = z[i] / K0[i]; // Heavy-phase composition
2539 }
2540 normalize_vector(xL);
2541 normalize_vector(xH);
2542
2543 // For each composition, use successive substitution to try to evaluate stability
2544 if (debug) {
2545 std::cout << format("3) SS2: i x' x'' rho' rho'' tpd' tpd''\n");
2546 }
2547
2548 // We got this far, we assume stable phases
2549 _stable = true;
2550
2551 double diffbulkL = 0, diffbulkH = 0;
2552 for (int step_count = 0; step_count < 100; ++step_count) {
2553
2554 // Set the composition
2555 HEOS.SatL->set_mole_fractions(xH);
2556 HEOS.SatV->set_mole_fractions(xL);
2557 HEOS.SatL->calc_reducing_state();
2558 HEOS.SatV->calc_reducing_state();
2559
2560 // Do the global density solver for both phases
2561 rho_TP_global();
2562
2563 double tpd_L = 0, tpd_H = 0;
2564 for (std::size_t i = 0; i < xL.size(); ++i) {
2565 tpd_L += xL[i] * (log(MixtureDerivatives::fugacity_i(*HEOS.SatV, i, XN_DEPENDENT)) - log(fugacity0[i]));
2566 tpd_H += xH[i] * (log(MixtureDerivatives::fugacity_i(*HEOS.SatL, i, XN_DEPENDENT)) - log(fugacity0[i]));
2567 }
2568 tpdL.push_back(tpd_L);
2569 tpdH.push_back(tpd_H);
2570
2571 // Calculate the new composition from the fugacity coefficients
2572 diffbulkL = 0, diffbulkH = 0;
2573 for (std::size_t i = 0; i < z.size(); ++i) {
2574 xL[i] = z[i] * fugacity_coefficient0[i] / HEOS.SatV->fugacity_coefficient(i);
2575 diffbulkL += std::abs(xL[i] - z[i]);
2576 xH[i] = z[i] * fugacity_coefficient0[i] / HEOS.SatL->fugacity_coefficient(i);
2577 diffbulkH += std::abs(xH[i] - z[i]);
2578 }
2579 normalize_vector(xL);
2580 normalize_vector(xH);
2581 if (debug) {
2582 std::cout << format("3) %d %s %s %g %g %g %g\n", step_count, vec_to_string(xL, "%0.6f").c_str(), vec_to_string(xH, "%0.6f").c_str(),
2583 rhomolar_liq, rhomolar_vap, tpd_L, tpd_H);
2584 }
2585
2586 // Check if either of the phases have the bulk composition. If so, no phase split
2587 if (diffbulkL < 1e-2 || diffbulkH < 1e-2) {
2588 _stable = true;
2589 return;
2590 }
2591
2592 // Check if either tpd is negative, if so, phases definitively split, quit
2593 if (tpd_L < -1e-12 || tpd_H < -1e-12) {
2594 _stable = false;
2595 return;
2596 }
2597 }
2598 if (diffbulkH > 0.25 || diffbulkL > 0.25) {
2599 // At least one test phase is definitely not the bulk composition, so phase split predicted
2600 _stable = false;
2601 }
2602}
2603
2604void StabilityRoutines::StabilityEvaluationClass::rho_TP_global() {
2605
2606 // Calculate the temperature and pressure to be used
2607 double the_T = (m_T > 0 && m_p > 0) ? m_T : HEOS.T();
2608 double the_p = (m_T > 0 && m_p > 0) ? m_p : HEOS.p();
2609
2610 // Calculate covolume of SRK, use it as the maximum density
2611 double rhoL = HEOS.SatL->solver_rho_Tp_global(the_T, the_p, HEOS.SatL->calc_rhomolar_max_bound());
2612 double rhoV = HEOS.SatV->solver_rho_Tp_global(the_T, the_p, HEOS.SatV->calc_rhomolar_max_bound());
2613 HEOS.SatL->update_DmolarT_direct(rhoL, the_T);
2614 HEOS.SatV->update_DmolarT_direct(rhoV, the_T);
2615
2616 rhomolar_liq = HEOS.SatL->rhomolar();
2617 rhomolar_vap = HEOS.SatV->rhomolar();
2618}
2619
2620void StabilityRoutines::StabilityEvaluationClass::rho_TP_w_guesses() {
2621
2622 // Re-calculate the density
2623 if (m_T > 0 && m_p > 0) {
2624 HEOS.SatL->update_TP_guessrho(m_T, m_p, rhomolar_liq);
2625 HEOS.SatV->update_TP_guessrho(m_T, m_p, rhomolar_vap);
2626 } else {
2627 HEOS.SatL->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_liq);
2628 HEOS.SatV->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_vap);
2629 }
2630 rhomolar_liq = HEOS.SatL->rhomolar();
2631 rhomolar_vap = HEOS.SatV->rhomolar();
2632}
2633
2634void StabilityRoutines::StabilityEvaluationClass::rho_TP_SRK_translated() {
2635
2636 // First use cubic as a guess for the density of liquid and vapor phases
2637 if (m_T > 0 && m_p > 0) {
2638 rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(m_T, m_p, iphase_liquid); // [mol/m^3]
2639 rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(m_T, m_p, iphase_gas); // [mol/m^3]
2640 } else {
2641 rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_liquid); // [mol/m^3]
2642 rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_gas); // [mol/m^3]
2643 }
2644
2645 // Apply volume translation to liquid density only
2646 if (HEOS.backend_name().find("Helmholtz") == 0) {
2647 // Use Peneloux volume translation to shift liquid volume
2648 // As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
2649 double summer_c = 0, v_SRK = 1 / rhomolar_liq;
2650 for (std::size_t i = 0; i < z.size(); ++i) {
2651 // Get the parameters for the cubic EOS
2652 CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical), pc = HEOS.get_fluid_constant(i, iP_critical),
2653 rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
2654 CoolPropDbl R = 8.3144598;
2655 summer_c += z[i] * (0.40768 * R * Tc / pc * (0.29441 - pc / (rhomolarc * R * Tc)));
2656 }
2657 rhomolar_liq = 1 / (v_SRK - summer_c);
2658 }
2659}
2660
2670void SaturationSolvers::PTflash_twophase::solve() {
2671 if (get_config_int(MIXTURE_STABILITY_ALGORITHM) != 0) {
2672 solve_michelsen();
2673 } else {
2674 solve_legacy();
2675 }
2676}
2677void SaturationSolvers::PTflash_twophase::solve_michelsen() {
2678 const std::size_t N = IO.x.size();
2679 if (!ValidNumber(IO.p)) IO.p = HEOS.p();
2680 if (!ValidNumber(IO.T)) IO.T = HEOS.T();
2681 // Reset the convergence-failure flag so it reflects only THIS solve attempt. IO is held
2682 // by reference, so a reused options object could otherwise carry a stale nonconvergence=true
2683 // from a prior attempt into the single-phase gate in PT_flash_mixtures.
2684 IO.nonconvergence = false;
2685 // Warm-start density roots for the liquid/vapor phases (tracked across SS + Newton
2686 // iterations; first solve per phase falls back to the global solver inside the helper).
2687 CoolPropDbl rho_warm_L = -1, rho_warm_V = -1;
2688
2689 // Store K-factors in log space to prevent overflow for wide-boiling mixtures.
2690 // lnK[i] = ln(phi_i^L / phi_i^V) = ln(K_i). See [Michelsen1982b] Eq. 5.
2691 std::vector<CoolPropDbl> lnK(N);
2692 for (std::size_t i = 0; i < N; ++i) {
2693 double ratio = IO.y[i] / std::max(IO.x[i], 1e-300);
2694 lnK[i] = std::log(std::max(ratio, 1e-300));
2695 }
2696 CoolPropDbl beta = IO.beta;
2697
2698 // Reject a non-finite seed up front. A stability false-positive at a single-phase point
2699 // (e.g. below the bubble) can hand in a NaN trial composition; without this guard the NaN
2700 // propagates through Rachford-Rice into x/y and surfaces as a misleading "lost a phase
2701 // density solve" error instead of a clean single-phase fallback (#3192). Reporting
2702 // nonconvergence routes PT_flash_mixtures to its single-phase branch.
2703 for (std::size_t i = 0; i < N; ++i) {
2704 if (!ValidNumber(lnK[i]) || !ValidNumber(IO.x[i]) || !ValidNumber(IO.y[i])) {
2705 IO.nonconvergence = true;
2706 throw SolutionError(format("PTflash_twophase::solve_michelsen got a non-finite seed at T = %g K, p = %g Pa", static_cast<double>(IO.T),
2707 static_cast<double>(IO.p)));
2708 }
2709 }
2710
2711 // Helper: solve Rachford-Rice in log-K space with Newton + bisection safeguards
2712 auto solve_rachford_rice = [&]() {
2713 CoolPropDbl beta_min = 0, beta_max = 1.0;
2714 // Guard a non-finite / out-of-range Newton start (e.g. a stale beta carried in from a
2715 // prior failed iterate): an invalid beta makes denom=1+beta*term hit a pole, poisoning
2716 // r/dr -> NaN below (#3192).
2717 if (!ValidNumber(beta) || beta < 0.0 || beta > 1.0) beta = 0.5;
2718 for (int rr_iter = 0; rr_iter < 50; ++rr_iter) {
2719 CoolPropDbl r = 0, dr = 0;
2720 for (std::size_t i = 0; i < N; ++i) {
2721 // Clamp lnK before exp() so a wide-boiling step (lnK grows unbounded via the
2722 // GDEM acceleration that feeds this solve) cannot overflow Ki -- nor term*term /
2723 // denom*denom below -- to +inf and poison r/dr with NaN. At the clamp Ki
2724 // dominates, so term/denom -> 1/beta and term^2/denom^2 -> 1/beta^2, i.e. the
2725 // exact Ki->inf asymptotic limits; 350 keeps term*term (~exp(2*lnK)) in range.
2726 // CoolPropDbl (not double) matches lnK's storage type.
2727 CoolPropDbl Ki = std::exp(std::min(lnK[i], static_cast<CoolPropDbl>(350.0)));
2728 CoolPropDbl term = Ki - 1.0;
2729 CoolPropDbl denom = 1.0 + beta * term;
2730 r += IO.z[i] * term / denom;
2731 dr -= IO.z[i] * term * term / (denom * denom);
2732 }
2733 if (r > 0)
2734 beta_min = beta;
2735 else
2736 beta_max = beta;
2737 if (std::abs(r) < 1e-11) break;
2738 CoolPropDbl beta_new = beta - r / dr;
2739 // A non-finite Newton step (dr -> 0, or a denom pole when beta strays to a bracket
2740 // edge) must fall back to bisection, not slip through: NaN passes BOTH the <= and >=
2741 // comparisons, so without the ValidNumber guard a NaN beta would poison x/y and the
2742 // downstream density solve (#3192; same failure class as the #3167 GDEM-ratio guard).
2743 if (!ValidNumber(beta_new) || beta_new <= beta_min || beta_new >= beta_max) beta_new = 0.5 * (beta_min + beta_max);
2744 if (std::abs(beta_new - beta) < 1e-11) break;
2745 beta = beta_new;
2746 }
2747 for (std::size_t i = 0; i < N; ++i) {
2748 // Same clamp as the residual loop: keeps a heavy component's Ki finite so x_i -> 0
2749 // and y_i -> z_i/beta (the correct Ki->inf limits) instead of inf/NaN.
2750 CoolPropDbl Ki = std::exp(std::min(lnK[i], static_cast<CoolPropDbl>(350.0)));
2751 IO.x[i] = IO.z[i] / (1.0 + beta * (Ki - 1.0));
2752 IO.y[i] = Ki * IO.x[i];
2753 }
2754 normalize_vector(IO.x);
2755 normalize_vector(IO.y);
2756 };
2757
2758 // Helper: evaluate phase densities and fugacities
2759 auto evaluate_phases = [&]() -> bool {
2760 HEOS.SatL->set_mole_fractions(IO.x);
2761 try {
2762 IO.rhomolar_liq = solve_trial_rho_warm(*HEOS.SatL, IO.T, IO.p, rho_warm_L);
2763 } catch (...) {
2764 return false;
2765 }
2766 HEOS.SatV->set_mole_fractions(IO.y);
2767 try {
2768 IO.rhomolar_vap = solve_trial_rho_warm(*HEOS.SatV, IO.T, IO.p, rho_warm_V);
2769 } catch (...) {
2770 return false;
2771 }
2772 return true;
2773 };
2774
2775 // --- Phase 1: Successive substitution with GDEM acceleration ---
2776 // K_i^new = phi_i^L(x) / phi_i^V(y). Exit when K-factors converge or
2777 // after max_ss_loops. See [Michelsen1982b] Eq. 5-9 and [M&M2007] Sec. 12.6.
2778 const int max_ss_loops = 4; // Each loop: 2 SS steps + 1 GDEM extrapolation
2779 const double ss_tol = 1e-7;
2780 bool ss_converged = false;
2781
2782 for (int loop = 0; loop < max_ss_loops && !ss_converged; ++loop) {
2783 std::array<double, 2> esq_pair = {0, 0};
2784 std::vector<CoolPropDbl> err(N);
2785
2786 for (int kk = 0; kk < 2 && !ss_converged; ++kk) {
2787 solve_rachford_rice();
2788 if (!evaluate_phases()) {
2789 throw SolutionError("PT flash lost a phase density solve during successive substitution");
2790 }
2791 double gmax = 0;
2792 double esq = 0;
2793 for (std::size_t i = 0; i < N; ++i) {
2794 double lnK_new = std::log(HEOS.SatL->fugacity_coefficient(i)) - std::log(HEOS.SatV->fugacity_coefficient(i));
2795 double diff = lnK_new - lnK[i];
2796 err[i] = diff;
2797 esq += IO.z[i] * diff * diff;
2798 gmax = std::max(gmax, std::abs(diff));
2799 lnK[i] = lnK_new;
2800 }
2801 esq_pair[kk] = esq;
2802 if (gmax < ss_tol) {
2803 ss_converged = true;
2804 }
2805 }
2806
2807 if (ss_converged) break;
2808
2809 // GDEM extrapolation ([M&M2007] Ch. 12, Sec. 12.6)
2810 if (esq_pair[0] > 0) {
2811 double ratio = std::sqrt(esq_pair[1] / esq_pair[0]);
2812 // A NaN ratio (non-finite esq) passes both < 0 and >= 0.95, so guard it
2813 // explicitly -- otherwise it poisons the GDEM lnK/Y update (CoolProp-1tbe.8 finding 4c).
2814 if (!ValidNumber(ratio) || ratio < 0 || ratio >= 0.95) ratio = 0.95;
2815 double factor = ratio / (1.0 - ratio);
2816 for (std::size_t i = 0; i < N; ++i) {
2817 lnK[i] += factor * err[i];
2818 }
2819 }
2820 }
2821
2822 // Ensure phases are up-to-date after SS
2823 solve_rachford_rice();
2824 if (!evaluate_phases()) {
2825 throw SolutionError("PT flash lost a phase density solve after successive substitution");
2826 }
2827
2828 // --- Phase 2: Second-order Gibbs energy minimization ---
2829 // See [Michelsen1982b] Appendix B. The reduced gradient is
2830 // g_i = beta_L*beta_V*(ln f_i^V - ln f_i^L) and the reduced Hessian
2831 // uses d(ln phi)/dn derivatives from both phases.
2832
2833 // Compute Gibbs energy of current state for objective-decrease checking
2834 auto compute_gibbs = [&]() -> double {
2835 double G = 0;
2836 for (std::size_t i = 0; i < N; ++i) {
2837 if (IO.x[i] > 0) G += (1.0 - beta) * IO.x[i] * (std::log(IO.x[i]) + std::log(HEOS.SatL->fugacity_coefficient(i)));
2838 if (IO.y[i] > 0) G += beta * IO.y[i] * (std::log(IO.y[i]) + std::log(HEOS.SatV->fugacity_coefficient(i)));
2839 }
2840 return G;
2841 };
2842
2843 double G_old = compute_gibbs();
2844
2845 // Convergence target on the equal-fugacity residual max_i |ln f_i^V - ln f_i^L|.
2846 // The post-loop gate below (GitHub #3168) throws when we exit worse than the SS
2847 // tolerance, restoring solve_legacy's "no silent wrong answer" contract.
2848 const double gibbs_tol = 1e-9;
2849 const int max_gibbs_iter = 50;
2850 const int max_restart = 2;
2851 const int max_inner = 40;
2852 bool converged = false;
2853 double last_max_g = 1e300;
2854
2855 for (int restart = 0; restart < max_restart && !converged; ++restart) {
2856 // Trust-region radius in scaled-variable (u = D*v) units; restart with a
2857 // tighter radius if a previous pass stalled.
2858 double trust_radius = (restart == 0) ? 1.0 : 0.2;
2859
2860 for (int gibbs_iter = 0; gibbs_iter < max_gibbs_iter && !converged; ++gibbs_iter) {
2861 const CoolPropDbl L_frac = 1.0 - beta;
2862 const CoolPropDbl V_frac = beta;
2863
2864 // Re-sync the SatL/SatV backends to the current (IO.x, IO.y) before reading
2865 // their fugacity coefficients and composition derivatives below. A rejected
2866 // or density-failed trial in the previous iteration's inner loop may have
2867 // left them on a trial composition; re-syncing here keeps the gradient,
2868 // Hessian, and the convergence residual (max_g) consistent with the state
2869 // that will actually be published.
2870 HEOS.SatL->set_mole_fractions(IO.x);
2871 HEOS.SatL->update_DmolarT_direct(IO.rhomolar_liq, IO.T);
2872 HEOS.SatV->set_mole_fractions(IO.y);
2873 HEOS.SatV->update_DmolarT_direct(IO.rhomolar_vap, IO.T);
2874
2875 // Reduced gradient g_i = beta*(1-beta)*(ln f_i^V - ln f_i^L), reduced
2876 // Hessian H ([Michelsen1982b] Eq. B6), and the diagonal conditioning
2877 // scale dia_i = sqrt(z_i/(x_i*y_i)) ([Michelsen1982b] Appendix B). The
2878 // raw Newton system is badly conditioned near the mixture critical point
2879 // and for wide-boiling mixtures, so we scale and then take a Hebden
2880 // restricted (trust-region) step.
2881 Eigen::VectorXd g(N), dia(N);
2882 Eigen::MatrixXd H(N, N);
2883 CoolPropDbl max_g = 0;
2884 Eigen::MatrixXd DL(N, N), DV(N, N);
2885 for (std::size_t i = 0; i < N; ++i) {
2886 for (std::size_t j = 0; j < N; ++j) {
2887 DL(i, j) = CoolProp::MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), i, j, CoolProp::XN_INDEPENDENT);
2888 DV(i, j) = CoolProp::MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), i, j, CoolProp::XN_INDEPENDENT);
2889 }
2890 }
2891 for (std::size_t i = 0; i < N; ++i) {
2892 CoolPropDbl sum_x_DL = 0, sum_y_DV = 0;
2893 for (std::size_t k = 0; k < N; ++k) {
2894 sum_x_DL += IO.x[k] * DL(i, k);
2895 sum_y_DV += IO.y[k] * DV(i, k);
2896 }
2897 for (std::size_t j = 0; j < N; ++j) {
2898 CoolPropDbl dln_phi_L_dnj = DL(i, j) - sum_x_DL;
2899 CoolPropDbl dln_phi_V_dnj = DV(i, j) - sum_y_DV;
2900 H(i, j) = V_frac * dln_phi_L_dnj + L_frac * dln_phi_V_dnj - 1.0;
2901 }
2902 H(i, i) += (V_frac / IO.x[i]) + (L_frac / IO.y[i]);
2903 CoolPropDbl l_act = std::log(IO.x[i]) + std::log(HEOS.SatL->fugacity_coefficient(i));
2904 CoolPropDbl v_act = std::log(IO.y[i]) + std::log(HEOS.SatV->fugacity_coefficient(i));
2905 g(i) = V_frac * L_frac * (v_act - l_act);
2906 max_g = std::max(max_g, std::abs(v_act - l_act));
2907 // Conditioning scale. The fallback to 1.0 only keeps the scaled system
2908 // D^{-1} H D^{-1} finite if dia underflows to ~0; it does not make a
2909 // zero-feed component (z_i = 0) usable here (the H and ln(x_i) terms
2910 // above already require z_i > 0, as in the original solver).
2911 CoolPropDbl dia_i = std::sqrt(IO.z[i] / std::max(IO.x[i] * IO.y[i], 1e-300));
2912 dia(i) = (dia_i > 1e-300) ? dia_i : 1.0;
2913 }
2914 last_max_g = max_g;
2915 if (max_g < gibbs_tol) {
2916 converged = true;
2917 break;
2918 }
2919
2920 // Scale the Newton system: Hs = D^{-1} H D^{-1}, gs = D^{-1} g (D = diag(dia)).
2921 Eigen::VectorXd gs(N);
2922 Eigen::MatrixXd Hs(N, N);
2923 for (std::size_t i = 0; i < N; ++i) {
2924 gs(i) = g(i) / dia(i);
2925 for (std::size_t j = 0; j < N; ++j) {
2926 Hs(i, j) = H(i, j) / (dia(i) * dia(j));
2927 }
2928 }
2929
2930 // Hebden restricted step: pick a diagonal shift so that (Hs + shift*I) is
2931 // positive definite (guaranteeing a descent direction; see [Michelsen1982b]
2932 // "Positive definiteness of the Hessian matrix") and the scaled step stays
2933 // within the trust radius. Accept on a Gibbs decrease, adapting the radius.
2934 double diagonal_shift = 0.0;
2935 bool step_ok = false;
2936 for (int inner = 0; inner < max_inner; ++inner) {
2937 Eigen::MatrixXd Hl = Hs;
2938 Hl.diagonal().array() += diagonal_shift;
2939
2940 // Positive-definiteness guard via the smallest eigenvalue (N is small).
2941 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(Hl, Eigen::EigenvaluesOnly);
2942 double min_eig = es.eigenvalues().minCoeff();
2943 if (!ValidNumber(min_eig)) break;
2944 if (min_eig < 1e-8) {
2945 diagonal_shift += (1e-8 - min_eig);
2946 continue;
2947 }
2948
2949 Eigen::VectorXd ds = Hl.ldlt().solve(-gs);
2950 double step_norm = ds.norm();
2951
2952 // Enforce the trust region by adding more shift if the step is too long.
2953 if (step_norm > trust_radius) {
2954 diagonal_shift = (diagonal_shift > 0.0) ? diagonal_shift * 3.0 : (step_norm / trust_radius - 1.0) * std::max(min_eig, 1e-3);
2955 continue;
2956 }
2957
2958 // Unscale to the mole-number step delta_v_i = ds_i / dia_i and clamp to
2959 // keep vapor/liquid mole numbers positive.
2960 Eigen::VectorXd delta_v(N);
2961 for (std::size_t i = 0; i < N; ++i) {
2962 delta_v(i) = ds(i) / dia(i);
2963 }
2964 CoolPropDbl pos_scale = 1.0;
2965 for (std::size_t i = 0; i < N; ++i) {
2966 CoolPropDbl v_old = beta * IO.y[i];
2967 if (delta_v(i) > 0 && v_old + delta_v(i) > IO.z[i]) pos_scale = std::min(pos_scale, 0.99 * (IO.z[i] - v_old) / delta_v(i));
2968 if (delta_v(i) < 0 && v_old + delta_v(i) < 0) pos_scale = std::min(pos_scale, 0.99 * (-v_old) / delta_v(i));
2969 }
2970
2971 CoolPropDbl V_new = 0, L_new = 0;
2972 std::vector<CoolPropDbl> v_new(N), l_new(N);
2973 for (std::size_t i = 0; i < N; ++i) {
2974 v_new[i] = beta * IO.y[i] + pos_scale * delta_v(i);
2975 l_new[i] = IO.z[i] - v_new[i];
2976 V_new += v_new[i];
2977 L_new += l_new[i];
2978 }
2979 std::vector<CoolPropDbl> x_trial(N), y_trial(N);
2980 for (std::size_t i = 0; i < N; ++i) {
2981 y_trial[i] = v_new[i] / V_new;
2982 x_trial[i] = l_new[i] / L_new;
2983 }
2984
2985 bool eval_ok = false;
2986 try {
2987 HEOS.SatL->set_mole_fractions(x_trial);
2988 CoolPropDbl rL = HEOS.SatL->solver_rho_Tp_global(IO.T, IO.p, HEOS.SatL->calc_rhomolar_max_bound());
2989 HEOS.SatV->set_mole_fractions(y_trial);
2990 CoolPropDbl rV = HEOS.SatV->solver_rho_Tp_global(IO.T, IO.p, HEOS.SatV->calc_rhomolar_max_bound());
2991 if (ValidNumber(rL) && ValidNumber(rV) && rL > 0 && rV > 0) {
2992 HEOS.SatL->update_DmolarT_direct(rL, IO.T);
2993 HEOS.SatV->update_DmolarT_direct(rV, IO.T);
2994
2995 // Tentatively adopt the trial state and check the Gibbs energy.
2996 double beta_save = beta;
2997 auto x_save = IO.x;
2998 auto y_save = IO.y;
2999 CoolPropDbl rL_save = IO.rhomolar_liq, rV_save = IO.rhomolar_vap;
3000 beta = V_new;
3001 IO.x = x_trial;
3002 IO.y = y_trial;
3003 IO.rhomolar_liq = rL;
3004 IO.rhomolar_vap = rV;
3005 double G_new = compute_gibbs();
3006
3007 if (G_new < G_old + 1e-12) {
3008 G_old = G_new;
3009 step_ok = true;
3010 eval_ok = true;
3011 // Grow the trust region if the step rode its boundary.
3012 if (step_norm > 0.8 * trust_radius) trust_radius = std::min(2.0 * trust_radius, 1e3);
3013 } else {
3014 // Gibbs increased: restore composition and density. The
3015 // SatL/SatV backends are left on the rejected trial, but the
3016 // next inner attempt overwrites them and the next outer
3017 // iteration re-syncs them to IO.x/IO.y at its top.
3018 beta = beta_save;
3019 IO.x = x_save;
3020 IO.y = y_save;
3021 IO.rhomolar_liq = rL_save;
3022 IO.rhomolar_vap = rV_save;
3023 }
3024 }
3025 } catch (...) {
3026 // density solve threw — treated as a rejected step below
3027 }
3028 if (eval_ok) break;
3029 // Density failure or no decrease: shrink the trust region, reset shift.
3030 trust_radius = 0.5 * step_norm;
3031 if (trust_radius < 1e-10) break;
3032 diagonal_shift = 0.0;
3033 }
3034 if (!step_ok) break;
3035 }
3036 }
3037 IO.beta = beta;
3038
3039 // Recompute the equal-fugacity residual on the FINAL published (IO.x, IO.y) state
3040 // so the gate below tests exactly what is returned, not a pre-step value captured
3041 // at the top of the last iteration (which could over-throw a split that converged
3042 // on its final accepted step).
3043 {
3044 HEOS.SatL->set_mole_fractions(IO.x);
3045 HEOS.SatL->update_DmolarT_direct(IO.rhomolar_liq, IO.T);
3046 HEOS.SatV->set_mole_fractions(IO.y);
3047 HEOS.SatV->update_DmolarT_direct(IO.rhomolar_vap, IO.T);
3048 double final_max_g = 0;
3049 for (std::size_t i = 0; i < N; ++i) {
3050 double l_act = std::log(IO.x[i]) + std::log(HEOS.SatL->fugacity_coefficient(i));
3051 double v_act = std::log(IO.y[i]) + std::log(HEOS.SatV->fugacity_coefficient(i));
3052 final_max_g = std::max(final_max_g, std::abs(v_act - l_act));
3053 }
3054 last_max_g = final_max_g;
3055 converged = (final_max_g < gibbs_tol);
3056 }
3057
3058 // Convergence gate (GitHub #3168, refined for #3192): never publish a grossly-unconverged
3059 // or trivial split -- the stability-false-positive signature #3168 targets -- but DO accept
3060 // a genuine, near-converged equilibrium. The original hard 1e-7 throw also discarded real
3061 // wide-boiling splits that converge only to ~1e-6 (SS converges linearly and the second-order
3062 // stage stalls for stiff mixtures), so a true two-phase state was silently misclassified as
3063 // single-phase. Distinguish the two: a genuine split has a non-trivial composition spread
3064 // AND an interior phase fraction AND an equal-fugacity residual at engineering tolerance; a
3065 // false positive is trivial (x==y), collapsed (beta -> 0/1), or grossly unconverged.
3066 if (!converged) {
3067 CoolPropDbl spread = 0;
3068 for (std::size_t i = 0; i < N; ++i)
3069 spread = std::max(spread, std::abs(IO.x[i] - IO.y[i]));
3070 const bool genuine = ValidNumber(last_max_g) && last_max_g <= 1e-5 // near-converged equilibrium
3071 && spread >= 1e-4 // not a trivial (x==y) split
3072 && beta > 1e-8 && beta < 1.0 - 1e-8; // not a collapsed phase
3073 if (!genuine) {
3074 IO.nonconvergence = true;
3075 throw SolutionError(format("PTflash_twophase::solve_michelsen failed to converge: max|ln f_V - ln f_L| = %g at T = %g K, p = %g Pa",
3076 last_max_g, static_cast<double>(IO.T), static_cast<double>(IO.p)));
3077 }
3078 }
3079}
3080
3081void SaturationSolvers::PTflash_twophase::solve_legacy() {
3082 const std::size_t N = IO.x.size();
3083 int iter = 0;
3084 double min_rel_change = NAN;
3085 do {
3086 // Build the Jacobian and residual vectors
3087 build_arrays();
3088
3089 // Solve for the step; v is the step with the contents
3090 // [delta(x'_0), delta(x'_1), ..., delta(x'_{N-1}), delta(x''_0), delta(x''_1), ..., delta(x''_{N-1})]
3091
3092 // Uncomment to see Jacobian and residual at every step
3093 //std::cout << vec_to_string(J, "%0.12Lg") << std::endl;
3094 //std::cout << vec_to_string(-r, "%0.12Lg") << std::endl;
3095
3096 Eigen::VectorXd v = J.colPivHouseholderQr().solve(-r);
3097
3098 for (unsigned int i = 0; i < N - 1; ++i) {
3099 err_rel[i] = v[i] / IO.x[i];
3100 IO.x[i] += v[i];
3101 err_rel[i + (N - 1)] = v[i + (N - 1)] / IO.y[i];
3102 IO.y[i] += v[i + (N - 1)];
3103 }
3104 IO.x[N - 1] = 1 - std::accumulate(IO.x.begin(), IO.x.end() - 1, 0.0);
3105 IO.y[N - 1] = 1 - std::accumulate(IO.y.begin(), IO.y.end() - 1, 0.0);
3106
3107 //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg") << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
3108
3109 min_rel_change = err_rel.cwiseAbs().minCoeff();
3110 iter++;
3111
3112 if (iter == IO.Nstep_max) {
3113 throw ValueError(format("PTflash_twophase::call reached max number of iterations [%d]", IO.Nstep_max));
3114 }
3115 } while (this->error_rms > 1e-9 && min_rel_change > 1000 * DBL_EPSILON && iter < IO.Nstep_max);
3116
3117 // Recover the vapor molar fraction (beta) from the converged overall mass
3118 // balance z_i = (1 - beta) x_i + beta y_i. The amount-of-substance
3119 // residuals in build_arrays() drive (z_i - x_i)/(y_i - x_i) to a single
3120 // common value across all components, so any component yields the same
3121 // beta; pick the one with the largest phase separation |y_i - x_i| for
3122 // numerical robustness. Without this the caller (PT_flash_mixtures) reads
3123 // the default beta = 0.5 and reports a wrong vapor quality and bulk density
3124 // for every two-phase mixture (CoolProp-1tbe.1); solve_michelsen() sets
3125 // IO.beta with the same convention.
3126 std::size_t i_best = 0;
3127 CoolPropDbl sep_best = -1;
3128 for (std::size_t i = 0; i < N; ++i) {
3129 CoolPropDbl sep = std::abs(IO.y[i] - IO.x[i]);
3130 if (sep > sep_best) {
3131 sep_best = sep;
3132 i_best = i;
3133 }
3134 }
3135 if (sep_best > 10 * DBL_EPSILON) {
3136 IO.beta = (IO.z[i_best] - IO.x[i_best]) / (IO.y[i_best] - IO.x[i_best]);
3137 } else {
3138 // Degenerate (trivial) split: the phases are numerically identical, so
3139 // beta is undefined. Signal a non-split so the caller falls back to
3140 // the single-phase branch rather than dividing by ~zero.
3141 IO.beta = 0;
3142 }
3143}
3144void SaturationSolvers::PTflash_twophase::build_arrays() {
3145 const std::size_t N = IO.x.size();
3146
3147 r.resize(2 * N - 2);
3148 J.resize(2 * N - 2, 2 * N - 2);
3149 err_rel.resize(2 * N - 2);
3150
3151 HEOS.SatL->set_mole_fractions(IO.x);
3152 HEOS.SatL->update_TP_guessrho(IO.T, IO.p, IO.rhomolar_liq);
3153
3154 HEOS.SatV->set_mole_fractions(IO.y);
3155 HEOS.SatV->update_TP_guessrho(IO.T, IO.p, IO.rhomolar_vap);
3156
3157 // Independent variables are
3158 // [delta(x'_0), delta(x'_1), ..., delta(x'_{N-1}), delta(x''_0), delta(x''_1), ..., delta(x''_{N-1})]
3159
3160 // First N residuals are the iso-fugacity condition
3161 for (std::size_t k = 0; k < N; ++k) {
3162 r(k) = log(HEOS.SatL->fugacity(k) / HEOS.SatV->fugacity(k));
3163 for (std::size_t j = 0; j < N - 1; ++j) {
3164 if (k == N - 1) {
3165 J(k, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), k, j, XN_DEPENDENT);
3166 J(k, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), k, j, XN_DEPENDENT);
3167 } else {
3168 J(k, j) = MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatL.get()), k, j, XN_DEPENDENT);
3169 J(k, j + N - 1) = -MixtureDerivatives::dln_fugacity_dxj__constT_p_xi(*(HEOS.SatV.get()), k, j, XN_DEPENDENT);
3170 }
3171 }
3172 }
3173 // Next N-2 residuals are amount of substance balances
3174 for (std::size_t i = 0; i < N - 2; ++i) {
3175 std::size_t k = i + N;
3176 r(k) = (IO.z[i] - IO.x[i]) / (IO.y[i] - IO.x[i]) - (IO.z[N - 1] - IO.x[N - 1]) / (IO.y[N - 1] - IO.x[N - 1]);
3177 for (std::size_t j = 0; j < N - 2; ++j) {
3178 J(k, j) = (IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
3179 J(k, j + N - 1) = -(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
3180 }
3181 std::size_t j = N - 2;
3182 J(k, j) = -(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
3183 J(k, j + N - 1) = +(IO.z[j] - IO.x[j]) / POW2(IO.y[j] - IO.x[j]);
3184 }
3185 this->error_rms = r.norm();
3186}
3187} /* namespace CoolProp*/
3188
3189#if defined(ENABLE_CATCH)
3190# include <catch2/catch_all.hpp>
3191
3192TEST_CASE("Check the PT flash calculation for two-phase inputs", "[PTflash_twophase]") {
3193 shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", "Propane&Ethane"));
3194 AS->set_mole_fractions(std::vector<double>(2, 0.5));
3195 // Dewpoint calculation
3196 AS->update(CoolProp::PQ_INPUTS, 101325, 1);
3197
3198 // Do a dummy calculation at the dewpoint state - make sure all the residuals are zero as they should be
3199 CoolProp::SaturationSolvers::PTflash_twophase_options o;
3200 o.x = AS->mole_fractions_liquid();
3201 o.y = AS->mole_fractions_vapor();
3202 o.z = AS->get_mole_fractions();
3203 o.rhomolar_liq = AS->saturated_liquid_keyed_output(CoolProp::iDmolar);
3204 o.rhomolar_vap = AS->saturated_vapor_keyed_output(CoolProp::iDmolar);
3205 o.T = AS->T();
3206 o.p = AS->p();
3207 o.omega = 1.0;
3208 CoolProp::SaturationSolvers::PTflash_twophase solver(*static_cast<CoolProp::HelmholtzEOSMixtureBackend*>(AS.get()), o);
3209 solver.build_arrays();
3210 double err = solver.r.norm();
3211 REQUIRE(std::abs(err) < 1e-10);
3212
3213 // Now, perturb the solution a little bit and actually do the solve
3214 std::vector<double> x0 = o.x;
3215 o.x[0] *= 1.1;
3216 o.x[1] = 1 - o.x[0];
3217 solver.solve();
3218 // Make sure we end up with the same liquid composition
3219 double diffx0 = o.x[0] - x0[0];
3220 REQUIRE(std::abs(diffx0) < 1e-10);
3221
3222 // Now do the blind flash call with PT as inputs
3223 AS->update(CoolProp::PT_INPUTS, AS->p(), AS->T() - 2);
3224 REQUIRE(AS->phase() == CoolProp::iphase_twophase);
3225}
3226
3227TEST_CASE("Legacy PT flash recovers the vapor quality (not pinned at 0.5)", "[PTflash_twophase]") {
3228 // Regression for CoolProp-1tbe.1: solve_legacy() never wrote IO.beta, so
3229 // a blind PT flash with MIXTURE_STABILITY_ALGORITHM=0 read the default
3230 // beta = 0.5 and reported a wrong vapor quality (and bulk density) for
3231 // every two-phase mixture. Build a reference two-phase point at a known,
3232 // non-trivial quality, then confirm both stability algorithms recover it.
3233 const std::string fluids = "Methane&Ethane";
3234 const std::vector<double> z = {0.6, 0.4};
3235 const double Qtarget = 0.25;
3236
3237 // RAII so the global config is restored even if a CHECK below throws.
3238 struct StabilityAlgoGuard
3239 {
3240 int prev;
3241 StabilityAlgoGuard() : prev(CoolProp::get_config_int(MIXTURE_STABILITY_ALGORITHM)) {}
3242 StabilityAlgoGuard(const StabilityAlgoGuard&) = delete;
3243 StabilityAlgoGuard& operator=(const StabilityAlgoGuard&) = delete;
3244 StabilityAlgoGuard(StabilityAlgoGuard&&) = delete;
3245 StabilityAlgoGuard& operator=(StabilityAlgoGuard&&) = delete;
3246 ~StabilityAlgoGuard() {
3247 CoolProp::set_config_int(MIXTURE_STABILITY_ALGORITHM, prev);
3248 }
3249 } stability_guard;
3250
3251 // Reference point at a known quality, built with the default algorithm.
3252 double p_ref, T_ref;
3253 {
3254 shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", fluids));
3255 AS->set_mole_fractions(z);
3256 AS->update(CoolProp::PQ_INPUTS, 2e6, Qtarget);
3257 p_ref = AS->p();
3258 T_ref = AS->T();
3259 }
3260
3261 auto flash_Q = [&](int algo) {
3262 CoolProp::set_config_int(MIXTURE_STABILITY_ALGORITHM, algo);
3263 shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS", fluids));
3264 AS->set_mole_fractions(z);
3265 AS->update(CoolProp::PT_INPUTS, p_ref, T_ref);
3266 REQUIRE(AS->phase() == CoolProp::iphase_twophase);
3267 return AS->Q();
3268 };
3269
3270 const double Q_michelsen = flash_Q(1);
3271 const double Q_legacy = flash_Q(0);
3272
3273 // Both algorithms recover the target quality...
3274 CHECK(Q_michelsen == Catch::Approx(Qtarget).margin(0.02));
3275 CHECK(Q_legacy == Catch::Approx(Qtarget).margin(0.02));
3276 // ...and the legacy path agrees with Michelsen rather than being pinned
3277 // at the default 0.5 (the bug).
3278 CHECK(Q_legacy == Catch::Approx(Q_michelsen).margin(0.01));
3279 CHECK(std::abs(Q_legacy - 0.5) > 0.1);
3280}
3281
3282#endif
Generated by doxygen 1.9.4