diff --git a/src/Backends/Tabular/BicubicBackend.cpp b/src/Backends/Tabular/BicubicBackend.cpp index dcf9ea62..0d2c1699 100644 --- a/src/Backends/Tabular/BicubicBackend.cpp +++ b/src/Backends/Tabular/BicubicBackend.cpp @@ -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"); } } diff --git a/src/Backends/Tabular/BicubicBackend.h b/src/Backends/Tabular/BicubicBackend.h index 7dee908b..777b5477 100644 --- a/src/Backends/Tabular/BicubicBackend.h +++ b/src/Backends/Tabular/BicubicBackend.h @@ -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; } }; diff --git a/src/Backends/Tabular/TabularBackends.cpp b/src/Backends/Tabular/TabularBackends.cpp index 95229451..3146453d 100644 --- a/src/Backends/Tabular/TabularBackends.cpp +++ b/src/Backends/Tabular/TabularBackends.cpp @@ -28,15 +28,15 @@ template void load_table(T &table, const std::string &path_to_table throw UnableToLoadError(err); } std::vector newBuffer(raw.size()*5); - uLong newBufferSize = newBuffer.size(); + uLong newBufferSize = static_cast(newBuffer.size()); int code; do{ code = uncompress((unsigned char *)(&(newBuffer[0])), &newBufferSize, - (unsigned char *)(&(raw[0])), raw.size()); + (unsigned char *)(&(raw[0])), static_cast(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(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 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 buffer(sbuf.size()); - uLong outSize = buffer.size(); + uLong outSize = static_cast(buffer.size()); compress((unsigned char *)(&(buffer[0])), &outSize, - (unsigned char*)(sbuf.data()), sbuf.size()); + (unsigned char*)(sbuf.data()), static_cast(sbuf.size())); std::ofstream ofs2(zPath.c_str(), std::ofstream::binary); ofs2.write(&buffer[0], outSize); diff --git a/src/Backends/Tabular/TabularBackends.h b/src/Backends/Tabular/TabularBackends.h index 02df1678..fb8921b6 100644 --- a/src/Backends/Tabular/TabularBackends.h +++ b/src/Backends/Tabular/TabularBackends.h @@ -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 in the other variable const std::vector > & 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); diff --git a/src/Solvers.cpp b/src/Solvers.cpp index 81aaf433..cce22a8a 100644 --- a/src/Solvers.cpp +++ b/src/Solvers.cpp @@ -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); diff --git a/wrappers/Python/CoolProp/Plots/ConsistencyPlots.py b/wrappers/Python/CoolProp/Plots/ConsistencyPlots.py index 8cf16919..66dff731 100644 --- a/wrappers/Python/CoolProp/Plots/ConsistencyPlots.py +++ b/wrappers/Python/CoolProp/Plots/ConsistencyPlots.py @@ -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() \ No newline at end of file