Skip to content

Commit

Permalink
Fix some bugs with cubics and critical point evaluation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ibell
ibell committed Feb 16, 2016
1 parent d51d47e commit 4ad9e02
Show file tree
Hide file tree
Showing 5 changed files with 46 additions and 25 deletions.
17 changes: 9 additions & 8 deletions src/Backends/Cubics/CubicBackend.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ void CoolProp::AbstractCubicBackend::rho_Tp_cubic(CoolPropDbl T, CoolPropDbl p,
A + B*B*(Delta_1*Delta_2-Delta_1-Delta_2) - B*(Delta_1+Delta_2),
-A*B-Delta_1*Delta_2*(POW2(B)+POW3(B)),
Nsolns, Z0, Z1, Z2);
if (Nsolns == 0){ rho0 = p/(Z0*R*T); }
if (Nsolns == 1){ rho0 = p/(Z0*R*T); }
else if (Nsolns == 3){
rho0 = p/(Z0*R*T);
rho1 = p/(Z1*R*T);
Expand Down Expand Up @@ -274,11 +274,7 @@ CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp(CoolPropDbl T, CoolPro
rho = rho0;
}
else if (Nsoln == 3){
if (rho_guess > 0 && imposed_phase_index == iphase_not_imposed){
// Use guessed density to select root

}
else if (rho_guess < 0 && imposed_phase_index != iphase_not_imposed){
if (imposed_phase_index != iphase_not_imposed){
// Use imposed phase to select root
if (imposed_phase_index == iphase_gas || imposed_phase_index == iphase_supercritical_gas){
rho = rho0;
Expand All @@ -297,8 +293,13 @@ CoolPropDbl CoolProp::AbstractCubicBackend::solver_rho_Tp(CoolPropDbl T, CoolPro
else{
throw ValueError("Obtained neither 1 nor three roots");
}
// Set some variables at the end
this->recalculate_singlephase_phase();
if (is_pure_or_pseudopure){
// Set some variables at the end
this->recalculate_singlephase_phase();
}
else{
_phase = iphase_gas; // TODO: fix this
}
_Q = -1;
return rho;
}
4 changes: 2 additions & 2 deletions src/Backends/Cubics/GeneralizedCubic.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -344,7 +344,7 @@ double AbstractCubic::d3_PI_12_dxidxjdxk(double delta, const std::vector<double>
}
double AbstractCubic::psi_plus(double delta, const std::vector<double> &x, std::size_t idelta)
{
double b = bm_term(x);
//double b = bm_term(x);
switch(idelta){
case 0:
return A_term(delta, x)*c_term(x)/(Delta_1-Delta_2);
Expand Down Expand Up @@ -468,7 +468,7 @@ double AbstractCubic::d2_psi_plus_dxidxj(double delta, const std::vector<double>
}
double AbstractCubic::d3_psi_plus_dxidxjdxk(double delta, const std::vector<double> &x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k, bool xN_independent)
{
double bracket = 0;
//double bracket = 0;
double PI12 = PI_12(delta, x, 0);
switch (idelta){
case 0:
Expand Down
6 changes: 3 additions & 3 deletions src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3279,8 +3279,8 @@ std::vector<CoolProp::CriticalState> HelmholtzEOSMixtureBackend::calc_all_critic
tau0 *= 1.1;
}
double tau_L0 = Halley(resid_L0, tau0, 1e-10, 100, errstr);
double T0 = T_reducing()/tau_L0;
double rho0 = delta0*rhomolar_reducing();
//double T0 = T_reducing()/tau_L0;
//double rho0 = delta0*rhomolar_reducing();

L0CurveTracer tracer(*this, tau_L0, delta0);

Expand All @@ -3303,11 +3303,11 @@ double HelmholtzEOSMixtureBackend::calc_tangent_plane_distance(const double T, c
}
update_TPD_state();
TPD_state->set_mole_fractions(std::vector<CoolPropDbl>(w.begin(), w.end()));
TPD_state->specify_phase(iphase_liquid); // Something homogeneous
if (rhomolar_guess < 0){
TPD_state->update(PT_INPUTS, p, T);
}
else{
TPD_state->specify_phase(iphase_gas); // Something homogeneous
TPD_state->update_TP_guessrho(T, p, rhomolar_guess);
}
double summer = 0;
Expand Down
4 changes: 4 additions & 0 deletions src/Backends/Helmholtz/ReducingFunctions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ CoolPropDbl ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector<CoolPropD
}
else if (xN_flag == XN_DEPENDENT){
if (j == N-1){ return 0;}
if (N == 0){ return 0; }
CoolPropDbl s = 0;
for (std::size_t k = 0; k < N-1; k++)
{
Expand All @@ -37,6 +38,7 @@ CoolPropDbl ReducingFunction::d2_ndTrdni_dxj_dxk__constxi(const std::vector<Cool
return d3Trdxidxjdxk(x, i, j, k, xN_flag) - 2*d2Trdxidxj(x, j, k, xN_flag)-s;
}
else if (xN_flag == XN_DEPENDENT){
if (N == 0){ return 0; }
if (j == N-1){ return 0; }
CoolPropDbl s = 0;
for (std::size_t m = 0; m < N-1; m++)
Expand Down Expand Up @@ -80,6 +82,7 @@ CoolPropDbl ReducingFunction::ndrhorbardni__constnj(const std::vector<CoolPropDb
}
else if (xN_flag == XN_DEPENDENT){
CoolPropDbl summer_term1 = 0;
if (N == 0){ return 0; }
for (std::size_t k = 0; k < N-1; ++k)
{
summer_term1 += x[k]*drhormolardxi__constxj(x, k, xN_flag);
Expand All @@ -103,6 +106,7 @@ CoolPropDbl ReducingFunction::ndTrdni__constnj(const std::vector<CoolPropDbl> &x
}
else if (xN_flag == XN_DEPENDENT){
CoolPropDbl summer_term1 = 0;
if (N==0){ return 0; }
for (std::size_t k = 0; k < N-1; ++k)
{
summer_term1 += x[k]*dTrdxi__constxj(x, k, xN_flag);
Expand Down
40 changes: 28 additions & 12 deletions src/Backends/Helmholtz/ReducingFunctions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,9 +64,9 @@ class ReducingFunction
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right)
* \f]
*/
CoolPropDbl d_ndTrdni_dxj__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d_ndTrdni_dxj__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);

CoolPropDbl d2_ndTrdni_dxj_dxk__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d2_ndTrdni_dxj_dxk__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);


/** \brief
Expand All @@ -82,18 +82,18 @@ class ReducingFunction
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-2}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right)
* \f]
*/
CoolPropDbl d_ndrhorbardni_dxj__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d_ndrhorbardni_dxj__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);

CoolPropDbl d2_ndrhorbardni_dxj_dxk__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d2_ndrhorbardni_dxj_dxk__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);

CoolPropDbl ndrhorbardni__constnj(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
CoolPropDbl ndTrdni__constnj(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
CoolPropDbl PSI_rho(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
CoolPropDbl d_PSI_rho_dxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
CoolPropDbl d2_PSI_rho_dxj_dxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
CoolPropDbl PSI_T(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
CoolPropDbl d_PSI_T_dxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
CoolPropDbl d2_PSI_T_dxj_dxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
virtual CoolPropDbl ndrhorbardni__constnj(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
virtual CoolPropDbl ndTrdni__constnj(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
virtual CoolPropDbl PSI_rho(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d_PSI_rho_dxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d2_PSI_rho_dxj_dxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
virtual CoolPropDbl PSI_T(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d_PSI_T_dxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
virtual CoolPropDbl d2_PSI_T_dxj_dxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
};

/** \brief The reducing function model of GERG-2008
Expand Down Expand Up @@ -389,6 +389,22 @@ class ConstantReducingFunction : public ReducingFunction
CoolPropDbl d2rhormolardxidxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return 0; };
/// \brief Derivative of the molar reducing density with respect to component i, j, and k mole fractions
CoolPropDbl d3rhormolardxidxjdxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ return 0; };

//virtual CoolPropDbl d_ndTrdni_dxj__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return 0; };
//virtual CoolPropDbl d2_ndTrdni_dxj_dxk__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ return 0; };
//virtual CoolPropDbl d_ndrhorbardni_dxj__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return 0; };
//virtual CoolPropDbl d2_ndrhorbardni_dxj_dxk__constxi(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ return 0; };
//virtual CoolPropDbl ndrhorbardni__constnj(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag){ return 0; };
//virtual CoolPropDbl ndTrdni__constnj(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag){ return 0; };

/// Note: this one is one, not zero
virtual CoolPropDbl PSI_rho(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag){ return 1; };
virtual CoolPropDbl d_PSI_rho_dxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return 0; };
virtual CoolPropDbl d2_PSI_rho_dxj_dxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ return 0; };
virtual CoolPropDbl PSI_T(const std::vector<CoolPropDbl> &x, std::size_t i, x_N_dependency_flag xN_flag){ return 0; };
virtual CoolPropDbl d_PSI_T_dxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ return 0; };
virtual CoolPropDbl d2_PSI_T_dxj_dxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ return 0; };

};

/** \brief Reducing function converter for dry air and HFC blends
Expand Down

0 comments on commit 4ad9e02

Please sign in to comment.