github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-21 03:48:08 -05:00
Code Issues Packages Projects Releases Wiki Activity

Get tabular backends working with consistency plots; closes #675

Browse Source 
Fixed a bunch of bugs with the Tabular methods
This commit is contained in:
Ian Bell
2015-05-16 20:48:41 -06:00
parent 4271f00f4f
commit 44bf43a65f
6 changed files with 142 additions and 69 deletions
Download Patch File Download Diff File Expand all files Collapse all files

72
src/Backends/Tabular/BicubicBackend.cpp
View File

@@ -63,7 +63,7 @@ void CoolProp::BicubicBackend::build_coeffs(SinglePhaseGriddedTableData &table,
f = &(table.umolar); fx = &(table.dumolardx); fy = &(table.dumolardy); fxy = &(table.d2umolardxdy);
break;
default:
throw ValueError();
throw ValueError("Invalid variable type to build_coeffs");
}
for (std::size_t i = 0; i < table.Nx-1; ++i) // -1 since we have one fewer cells than nodes
{
@@ -248,6 +248,16 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
// Find and cache the indices i, j
selected_table = SELECTED_PH_TABLE;
single_phase_logph.find_nearest_neighbor(iP, _p, otherkey, otherval, cached_single_phase_i, cached_single_phase_j);
CellCoeffs &cell = coeffs_ph[cached_single_phase_i][cached_single_phase_j];
if (!cell.valid()){
if (cell.has_valid_neighbor()){
// Get new good neighbor
cell.get_alternate(cached_single_phase_i, cached_single_phase_j);
}
else{
if (!cell.valid()){throw ValueError(format("Cell is invalid and has no good neighbors for p = %g Pa, T= %g K",val1,val2));}
}
}
// Now find hmolar given P, X for X in Hmolar, Smolar, Umolar
invert_single_phase_x(single_phase_logph, coeffs_ph, otherkey, otherval, _p, cached_single_phase_i, cached_single_phase_j);
}
@@ -291,6 +301,17 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
// Find and cache the indices i, j
selected_table = SELECTED_PT_TABLE;
single_phase_logpT.find_native_nearest_good_cell(_T, _p, cached_single_phase_i, cached_single_phase_j);
CellCoeffs &cell = coeffs_pT[cached_single_phase_i][cached_single_phase_j];
if (!cell.valid()){
if (cell.has_valid_neighbor()){
// Get new good neighbor
cell.get_alternate(cached_single_phase_i, cached_single_phase_j);
}
else{
if (!cell.valid()){throw ValueError(format("Cell is invalid and has no good neighbors for p = %g Pa, T= %g K",val1,val2));}
}
}
// If p < pc, you might be getting a liquid solution when you want a vapor solution or vice versa
// if you are very close to the saturation curve, so we figure out what the saturation temperature
// is for the given pressure
@@ -341,12 +362,23 @@ void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double v
else{
cached_saturation_iL = iL; cached_saturation_iV = iV;
}
_p = pure_saturation.evaluate(iP, _T, _Q, iL, iV);
}
else{
// Find and cache the indices i, j
selected_table = SELECTED_PT_TABLE;
single_phase_logpT.find_nearest_neighbor(iT, _T, otherkey, otherval, cached_single_phase_i, cached_single_phase_j);
// Now find the y variable (Dmolar in this case)
CellCoeffs &cell = coeffs_pT[cached_single_phase_i][cached_single_phase_j];
if (!cell.valid()){
if (cell.has_valid_neighbor()){
// Get new good neighbor
cell.get_alternate(cached_single_phase_i, cached_single_phase_j);
}
else{
if (!cell.valid()){throw ValueError(format("Cell is invalid and has no good neighbors for p = %g Pa, T= %g K",val1,val2));}
}
}
// Now find the y variable (Dmolar or Smolar in this case)
invert_single_phase_y(single_phase_logpT, coeffs_pT, otherkey, otherval, _T, cached_single_phase_i, cached_single_phase_j);
}
break;
@@ -416,7 +448,7 @@ double CoolProp::BicubicBackend::evaluate_single_phase_transport(SinglePhaseGrid
switch(output){
case iconductivity: _conductivity = val; break;
case iviscosity: _viscosity = val; break;
default: throw ValueError();
default: throw ValueError("Invalid output variable in evaluate_single_phase_transport");
}
return val;
}
@@ -449,8 +481,8 @@ double CoolProp::BicubicBackend::evaluate_single_phase(const SinglePhaseGriddedT
case iDmolar: _rhomolar = val; break;
case iSmolar: _smolar = val; break;
case iHmolar: _hmolar = val; break;
//case iUmolar:
default: throw ValueError();
case iUmolar: _umolar = val; break;
default: throw ValueError("Invalid output variable in evaluate_single_phase");
}
return val;
}
@@ -530,10 +562,19 @@ void CoolProp::BicubicBackend::invert_single_phase_x(const SinglePhaseGriddedTab
xhat = xhat0;
}
else if (N == 2){
xhat = std::min(xhat0, xhat1);
xhat = std::abs(xhat0) < std::abs(xhat1) ? xhat0 : xhat1;
}
else if (N == 3){
xhat = min3(xhat0, xhat1, xhat2);
if (std::abs(xhat0) < std::abs(xhat1) && std::abs(xhat0) < std::abs(xhat2)){
xhat = xhat0;
}
// Already know that xhat1 < xhat0 (xhat0 is not the minimum)
else if (std::abs(xhat1) < std::abs(xhat2)){
xhat = xhat1;
}
else{
xhat = xhat2;
}
}
else if (N == 0){
throw ValueError("Could not find a solution in invert_single_phase_x");
@@ -547,7 +588,7 @@ void CoolProp::BicubicBackend::invert_single_phase_x(const SinglePhaseGriddedTab
switch(table.xkey){
case iHmolar: _hmolar = val; break;
case iT: _T = val; break;
default: throw ValueError();
default: throw ValueError("Invalid output variable in invert_single_phase_x");
}
}
@@ -578,10 +619,19 @@ void CoolProp::BicubicBackend::invert_single_phase_y(const SinglePhaseGriddedTab
yhat = yhat0;
}
else if (N == 2){
yhat = std::min(yhat0, yhat1);
yhat = std::abs(yhat0) < std::abs(yhat1) ? yhat0 : yhat1;
}
else if (N == 3){
yhat = min3(yhat0, yhat1, yhat2);
if (std::abs(yhat0) < std::abs(yhat1) && std::abs(yhat0) < std::abs(yhat2)){
yhat = yhat0;
}
// Already know that yhat1 < yhat0 (yhat0 is not the minimum)
else if (std::abs(yhat1) < std::abs(yhat2)){
yhat = yhat1;
}
else{
yhat = yhat2;
}
}
else if (N == 0){
throw ValueError("Could not find a solution in invert_single_phase_x");
@@ -594,7 +644,7 @@ void CoolProp::BicubicBackend::invert_single_phase_y(const SinglePhaseGriddedTab
// Cache the output value calculated
switch(table.ykey){
case iP: _p = val; break;
default: throw ValueError();
default: throw ValueError("Invalid output variable in invert_single_phase_x");
}
}

6
src/Backends/Tabular/BicubicBackend.h
View File

@@ -49,7 +49,7 @@ class CellCoeffs{
}
};
/// Returns true if the cell coefficients seem to have been calculated properly
bool valid(){return _valid;};
bool valid() const {return _valid;};
/// Call this function to set the valid flag to true
void set_valid(){_valid = true;};
/// Call this function to set the valid flag to false
@@ -57,7 +57,7 @@ class CellCoeffs{
/// Set the neighboring (alternate) cell to be used if the cell is invalid
void set_alternate(std::size_t i, std::size_t j){alt_i = i; alt_j = j; _has_valid_neighbor = true;}
/// Get neighboring(alternate) cell to be used if this cell is invalid
void get_alternate(std::size_t &i, std::size_t &j){
void get_alternate(std::size_t &i, std::size_t &j) const {
if (_has_valid_neighbor){
i = alt_i; j = alt_j;
}
@@ -66,7 +66,7 @@ class CellCoeffs{
}
}
/// Returns true if cell is invalid and it has valid neighbor
bool has_valid_neighbor(){
bool has_valid_neighbor() const{
return _has_valid_neighbor;
}
};

10
src/Backends/Tabular/TabularBackends.cpp
View File

@@ -28,15 +28,15 @@ template <typename T> void load_table(T &table, const std::string &path_to_table
throw UnableToLoadError(err);
}
std::vector<char> newBuffer(raw.size()*5);
uLong newBufferSize = newBuffer.size();
uLong newBufferSize = static_cast<uLong>(newBuffer.size());
int code;
do{
code = uncompress((unsigned char *)(&(newBuffer[0])), &newBufferSize,
(unsigned char *)(&(raw[0])), raw.size());
(unsigned char *)(&(raw[0])), static_cast<mz_ulong>(raw.size()));
if (code == Z_BUF_ERROR){
// Output buffer is too small, make it bigger and try again
newBuffer.resize(newBuffer.size()*5);
newBufferSize = newBuffer.size();
newBufferSize = static_cast<uLong>(newBuffer.size());
}
else if (code != 0){ // Something else, a big problem
std::string err = format("Unable to uncompress file %s with miniz code %d", path_to_table.c_str(), code);
@@ -68,9 +68,9 @@ template <typename T> void write_table(const T &table, const std::string &path_t
std::string tabPath = std::string(path_to_tables + "/" + name + ".bin");
std::string zPath = tabPath + ".z";
std::vector<char> buffer(sbuf.size());
uLong outSize = buffer.size();
uLong outSize = static_cast<uLong>(buffer.size());
compress((unsigned char *)(&(buffer[0])), &outSize,
(unsigned char*)(sbuf.data()), sbuf.size());
(unsigned char*)(sbuf.data()), static_cast<mz_ulong>(sbuf.size()));
std::ofstream ofs2(zPath.c_str(), std::ofstream::binary);
ofs2.write(&buffer[0], outSize);

25
src/Backends/Tabular/TabularBackends.h
View File

@@ -395,7 +395,8 @@ class SinglePhaseGriddedTableData{
};
/// Check that the native inputs (the inputs the table is based on) are in range
bool native_inputs_are_in_range(double x, double y){
return x >= xmin && x <= xmax && y >= ymin && y <= ymax;
double e = 10*DBL_EPSILON;
return x >= xmin-e && x <= xmax+e && y >= ymin-e && y <= ymax+e;
}
/// @brief Find the nearest neighbor for native inputs (the inputs the table is based on)
/// Does not check whether this corresponds to a valid node or not
@@ -433,9 +434,6 @@ class SinglePhaseGriddedTableData{
// This one is fine because we now end up with a vector<double> in the other variable
const std::vector<std::vector<double> > & v = get(otherkey);
bisect_vector(v[i], otherval, j);
if (j < v[i].size()-1 && std::abs(v[i][j+1] - otherval) < std::abs(v[i][j] - otherval)){
j++;
}
}
}
/// Find the nearest good neighbor node for inputs that are the same as the grid inputs
@@ -627,7 +625,11 @@ class TabularBackend : public AbstractState
LogPTTable single_phase_logpT;
PureFluidSaturationTableData pure_saturation; // This will ultimately be split into pure and mixture backends which derive from this backend
PhaseEnvelopeData phase_envelope;
CoolPropDbl calc_T_critical(void){return this->AS->T_critical();};
CoolPropDbl calc_Ttriple(void){return this->AS->Ttriple();};
CoolPropDbl calc_p_triple(void){return this->AS->p_triple();};
CoolPropDbl calc_pmax(void){return this->AS->pmax();};
CoolPropDbl calc_Tmax(void){return this->AS->Tmax();};
bool using_mole_fractions(void){return true;}
bool using_mass_fractions(void){return false;}
bool using_volu_fractions(void){return false;}
@@ -725,6 +727,19 @@ class TabularBackend : public AbstractState
return pure_saturation.evaluate(iSmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
}
}
CoolPropDbl calc_umolar(void){
if (using_single_phase_table){
switch(selected_table){
case SELECTED_PH_TABLE: return evaluate_single_phase_phmolar(iUmolar, cached_single_phase_i, cached_single_phase_j);
case SELECTED_PT_TABLE: return evaluate_single_phase_pT(iUmolar, cached_single_phase_i, cached_single_phase_j);
case SELECTED_NO_TABLE: throw ValueError("table not selected");
}
return _HUGE; // not needed, will never be hit, just to make compiler happy
}
else{
return pure_saturation.evaluate(iUmolar, _p, _Q, cached_saturation_iL, cached_saturation_iV);
}
}
CoolPropDbl calc_cpmolar(void){
if (using_single_phase_table){
return calc_first_partial_deriv(iHmolar, iT, iP);

12
src/Solvers.cpp
View File

@@ -221,16 +221,8 @@ double Secant(FuncWrapper1D* f, double x0, double dx, double tol, int maxiter, s
if (iter>1)
{
double deltax = x2-x1;
if (std::abs(deltax)<1e-14)
{
if (std::abs(fval) < tol*10)
{
return x;
}
else
{
throw ValueError("Step is small but not solved to tolerance");
}
if (std::abs(deltax)<1e-14){
return x;
}
y2=fval;
x3=x2-y2/(y2-y1)*(x2-x1);

86
wrappers/Python/CoolProp/Plots/ConsistencyPlots.py
View File

@@ -5,7 +5,6 @@ import numpy as np
import time
import CoolProp as CP
from CoolProp.CoolProp import PropsSI
from CoolProp.Plots import PropsPlot
CP.CoolProp.set_debug_level(00)
@@ -64,18 +63,20 @@ def split_pair_xy(pair):
raise ValueError(pair)
class ConsistencyFigure(object):
def __init__(self, fluid, figsize = (15, 23)):
def __init__(self, fluid, figsize = (15, 23), backend = 'HEOS', additional_skips = []):
self.fluid = fluid
self.backend = backend
self.fig, self.axes = plt.subplots(nrows = 5, ncols = 3, figsize = figsize)
self.pairs = all_solvers
pairs_generator = iter(self.pairs)
states = [CP.AbstractState(backend, fluid) for _ in range(3)]
self.axes_list = []
for row in self.axes:
for ax in row:
pair = pairs_generator.next()
self.axes_list.append(ConsistencyAxis(ax, self, pair, self.fluid))
self.axes_list.append(ConsistencyAxis(ax, self, pair, self.fluid, self.backend, *states))
ax.set_title(pair)
self.calc_saturation_curves()
@@ -83,14 +84,14 @@ class ConsistencyFigure(object):
self.calc_Tmax_curve()
self.plot_Tmax_curve()
self.calc_melting_curve()
self.plot_melting_curve()
self.tight_layout()
for i, (ax, pair) in enumerate(zip(self.axes_list, self.pairs)):
if pair not in not_implemented_solvers:
if pair not in not_implemented_solvers and pair not in additional_skips:
ax.consistency_check_singlephase()
if pair not in no_two_phase_solvers:
ax.consistency_check_twophase()
@@ -104,7 +105,7 @@ class ConsistencyFigure(object):
"""
Calculate all the saturation curves in one shot using the state class to save computational time
"""
HEOS = CP.AbstractState('HEOS', self.fluid)
HEOS = CP.AbstractState(self.backend, self.fluid)
self.dictL, self.dictV = {}, {}
for Q, dic in zip([0, 1], [self.dictL, self.dictV]):
rhomolar,smolar,hmolar,T,p,umolar = [],[],[],[],[],[]
@@ -112,6 +113,9 @@ class ConsistencyFigure(object):
try:
HEOS.update(CP.QT_INPUTS, Q, _T)
if (HEOS.p() < 0): raise ValueError('P is negative:'+str(HEOS.p()))
HEOS.T(), HEOS.p(), HEOS.rhomolar(), HEOS.hmolar(), HEOS.smolar()
HEOS.umolar()
T.append(HEOS.T())
p.append(HEOS.p())
rhomolar.append(HEOS.rhomolar())
@@ -134,12 +138,17 @@ class ConsistencyFigure(object):
ax.plot_saturation_curves()
def calc_Tmax_curve(self):
HEOS = CP.AbstractState('HEOS', self.fluid)
HEOS = CP.AbstractState(self.backend, self.fluid)
rhomolar,smolar,hmolar,T,p,umolar = [],[],[],[],[],[]
for _p in np.logspace(np.log10(HEOS.keyed_output(CP.iP_min)*1.01), np.log10(HEOS.keyed_output(CP.iP_max)), 300):
try:
HEOS.update(CP.PT_INPUTS, _p, HEOS.keyed_output(CP.iT_max))
except ValueError as VE:
print('Tmax',_p, VE)
continue
try:
T.append(HEOS.T())
p.append(HEOS.p())
rhomolar.append(HEOS.rhomolar())
@@ -147,7 +156,7 @@ class ConsistencyFigure(object):
smolar.append(HEOS.smolar())
umolar.append(HEOS.umolar())
except ValueError as VE:
print('Tmax',VE)
print('Tmax access', VE)
self.Tmax = dict(T = np.array(T),
P = np.array(p),
@@ -204,11 +213,15 @@ class ConsistencyFigure(object):
self.fig.savefig(fname, **kwargs)
class ConsistencyAxis(object):
def __init__(self, axis, fig, pair, fluid):
def __init__(self, axis, fig, pair, fluid, backend, state1, state2, state3):
self.ax = axis
self.fig = fig
self.pair = pair
self.fluid = fluid
self.backend = backend
self.state = state1
self.state_PT = state2
self.state_QT = state3
#self.saturation_curves()
def label_axes(self):
@@ -255,7 +268,6 @@ class ConsistencyAxis(object):
def consistency_check_singlephase(self):
tic = time.time()
state = CP.AbstractState('HEOS', self.fluid)
# Update the state given the desired set of inputs
param1, param2 = split_pair(self.pair)
@@ -272,51 +284,53 @@ class ConsistencyAxis(object):
xbad, ybad = [], []
xexcep, yexcep = [], []
for p in np.logspace(np.log10(state.keyed_output(CP.iP_min)*1.01), np.log10(state.keyed_output(CP.iP_max)), 40):
for p in np.logspace(np.log10(self.state.keyed_output(CP.iP_min)*1.01), np.log10(self.state.keyed_output(CP.iP_max)), 40):
Tmin = state.keyed_output(CP.iT_triple)
if state.has_melting_line():
Tmin = self.state.keyed_output(CP.iT_triple)
if self.state.has_melting_line():
try:
pmelt_min = state.melting_line(CP.iP_min, -1, -1)
pmelt_min = self.state.melting_line(CP.iP_min, -1, -1)
if p < pmelt_min:
T0 = Tmin
else:
T0 = state.melting_line(CP.iT, CP.iP, p)
T0 = self.state.melting_line(CP.iT, CP.iP, p)
except Exception as E:
T0 = Tmin + 1.1
print('MeltingLine:', E)
else:
T0 = Tmin+1.1
for T in np.linspace(T0, state.keyed_output(CP.iT_max), 40):
state_PT = CP.AbstractState('HEOS', self.fluid)
for T in np.linspace(T0, self.state.keyed_output(CP.iT_max), 40):
try:
# Update the state using PT inputs in order to calculate all the remaining inputs
state_PT.update(CP.PT_INPUTS, p, T)
self.state_PT.update(CP.PT_INPUTS, p, T)
except ValueError as VE:
print('consistency',VE)
continue
_exception = False
try:
state.update(pairkey, state_PT.keyed_output(key1), state_PT.keyed_output(key2))
print(self.pair, self.state_PT.keyed_output(key1), self.state_PT.keyed_output(key2))
self.state.update(pairkey, self.state_PT.keyed_output(key1), self.state_PT.keyed_output(key2))
except ValueError as VE:
print('update', VE)
print('update', self.state_PT.keyed_output(key1), self.state_PT.keyed_output(key2), VE)
_exception = True
x = self.to_axis_units(xparam, state_PT.keyed_output(xkey))
y = self.to_axis_units(yparam, state_PT.keyed_output(ykey))
x = self.to_axis_units(xparam, self.state_PT.keyed_output(xkey))
y = self.to_axis_units(yparam, self.state_PT.keyed_output(ykey))
if _exception:
xexcep.append(x)
yexcep.append(y)
else:
# Check the error on the density
if abs(state_PT.rhomolar()/state.rhomolar()-1) < 1e-3 and abs(state_PT.p()/state.p()-1) < 1e-3 and abs(state_PT.T() - state.T()) < 1e-3:
if abs(self.state_PT.rhomolar()/self.state.rhomolar()-1) < 1e-3 and abs(self.state_PT.p()/self.state.p()-1) < 1e-3 and abs(self.state_PT.T() - self.state.T()) < 1e-3:
xgood.append(x)
ygood.append(y)
else:
print('bad', x, y, abs(self.state_PT.rhomolar()/self.state.rhomolar()-1), abs(self.state_PT.p()/self.state.p()-1), abs(self.state_PT.T() - self.state.T()))
xbad.append(x)
ybad.append(y)
@@ -330,7 +344,7 @@ class ConsistencyAxis(object):
def consistency_check_twophase(self):
tic = time.time()
state = CP.AbstractState('HEOS', self.fluid)
state = self.state
# Update the state given the desired set of inputs
param1, param2 = split_pair(self.pair)
@@ -351,32 +365,32 @@ class ConsistencyAxis(object):
Tmin = state.keyed_output(CP.iT_triple)+1
for T in np.linspace(Tmin, state.keyed_output(CP.iT_critical)-0.5, 20):
state_QT = CP.AbstractState('HEOS', self.fluid)
for T in np.linspace(Tmin, state.keyed_output(CP.iT_critical)-1, 20):
try:
# Update the state using QT inputs in order to calculate all the remaining inputs
state_QT.update(CP.QT_INPUTS, q, T)
self.state_QT.update(CP.QT_INPUTS, q, T)
except ValueError as VE:
print('consistency',VE)
continue
_exception = False
try:
state.update(pairkey, state_QT.keyed_output(key1), state_QT.keyed_output(key2))
state.update(pairkey, self.state_QT.keyed_output(key1), self.state_QT.keyed_output(key2))
except ValueError as VE:
print('update', state_QT.keyed_output(key1), state_QT.keyed_output(key2), VE)
print('update_QT', T, q)
print('update', self.state_QT.keyed_output(key1), self.state_QT.keyed_output(key2), VE)
_exception = True
x = self.to_axis_units(xparam, state_QT.keyed_output(xkey))
y = self.to_axis_units(yparam, state_QT.keyed_output(ykey))
x = self.to_axis_units(xparam, self.state_QT.keyed_output(xkey))
y = self.to_axis_units(yparam, self.state_QT.keyed_output(ykey))
if _exception:
xexcep.append(x)
yexcep.append(y)
else:
# Check the error on the density
if abs(state_QT.rhomolar()/state.rhomolar()-1) < 1e-3 and abs(state_QT.p()/state.p()-1) < 1e-3 and abs(state_QT.T() - state.T()) < 1e-3:
if abs(self.state_QT.rhomolar()/self.state.rhomolar()-1) < 1e-3 and abs(self.state_QT.p()/self.state.p()-1) < 1e-3 and abs(self.state_QT.T() - self.state.T()) < 1e-3:
xgood.append(x)
ygood.append(y)
else:
@@ -405,13 +419,15 @@ class ConsistencyAxis(object):
if __name__=='__main__':
PVT = PdfPages('Consistency.pdf')
for fluid in ['Ethanol','Water']:#CP.__fluids__:
for fluid in ['Water']:#CP.__fluids__:
print('************************************************')
print(fluid)
print('************************************************')
ff = ConsistencyFigure(fluid)
skips = ['DmolarHmolar','DmolarSmolar','DmolarUmolar','HmolarSmolar']
ff = ConsistencyFigure(fluid, backend = 'BICUBIC&HEOS', additional_skips = skips)
ff.add_to_pdf(PVT)
ff.savefig(fluid + '.png')
ff.savefig(fluid + '.pdf')
plt.close()
del ff
PVT.close()