From cf7ea4c49d36af9169c4db32413e8d47235bb7ff Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 26 Feb 2015 22:00:04 -0700 Subject: [PATCH] Bicubic tables work and are VERY fast Edge cases need handling still Signed-off-by: Ian Bell --- src/Backends/Tabular/BicubicBackend.cpp | 100 +++++++++++++++++++++--- src/Backends/Tabular/BicubicBackend.h | 9 ++- src/Backends/Tabular/TabularBackends.h | 2 + 3 files changed, 99 insertions(+), 12 deletions(-) diff --git a/src/Backends/Tabular/BicubicBackend.cpp b/src/Backends/Tabular/BicubicBackend.cpp index fadf4581..e11b748f 100644 --- a/src/Backends/Tabular/BicubicBackend.cpp +++ b/src/Backends/Tabular/BicubicBackend.cpp @@ -2,6 +2,7 @@ #include "MatrixMath.h" /// The inverse of the A matrix for the bicubic interpolation (http://en.wikipedia.org/wiki/Bicubic_interpolation) +/// NOTE: The matrix is transposed below static const double Ainv_data[16*16] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, @@ -19,18 +20,16 @@ static const double Ainv_data[16*16] = { 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 1, 0, 1, 0, -6, 6, 6, -6, -4, -2, 4, 2, -3, 3, -3, 3, -2, -1, -2, -1, 4, -4, -4, 4, 2, 2, -2, -2, 2, -2, 2, -2, 1, 1, 1, 1}; -static Eigen::Matrix Ainv(Ainv_data); - -void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double val1, double val2) -{ -} +static Eigen::Matrix Ainv(Ainv_data); + void CoolProp::BicubicBackend::build_coeffs_ph() { const bool debug = get_debug_level() > 5 || true; const int param_count = 5; - parameters param_list[param_count] = {iT, iDmolar, iSmolar, iHmolar, iP}; //iUmolar + parameters param_list[param_count] = {iDmolar, iT, iSmolar, iHmolar, iP}; //iUmolar std::vector > *f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL; + // To be passed as argument SinglePhaseGriddedTableData &table = single_phase_logph; std::vector > &coeffs = coeffs_ph; @@ -51,13 +50,13 @@ void CoolProp::BicubicBackend::build_coeffs_ph() f = &(table.p); fx = &(table.dpdx); fy = &(table.dpdy); fxy = &(table.d2pdxdy); break; case iDmolar: - f = &(table.rhomolar); fx = &(table.drhomolardx); *fy = (table.drhomolardy); fxy = &(table.d2rhomolardxdy); + f = &(table.rhomolar); fx = &(table.drhomolardx); fy = &(table.drhomolardy); fxy = &(table.d2rhomolardxdy); break; case iSmolar: - f = &(table.smolar); fx = &(table.dsmolardx); *fy = (table.dsmolardy); fxy = &(table.d2smolardxdy); + f = &(table.smolar); fx = &(table.dsmolardx); fy = &(table.dsmolardy); fxy = &(table.d2smolardxdy); break; case iHmolar: - f = &(table.hmolar); fx = &(table.dhmolardx); *fy = (table.dhmolardy); fxy = &(table.d2hmolardxdy); + f = &(table.hmolar); fx = &(table.dhmolardx); fy = &(table.dhmolardy); fxy = &(table.d2hmolardxdy); break; default: throw ValueError(); @@ -88,7 +87,7 @@ void CoolProp::BicubicBackend::build_coeffs_ph() F(12) = (*fxy)[i][j]*dy_dyhat*dx_dxhat; F(13) = (*fxy)[i+1][j]*dy_dyhat*dx_dxhat; F(14) = (*fxy)[i][j+1]*dy_dyhat*dx_dxhat; F(15) = (*fxy)[i+1][j+1]*dy_dyhat*dx_dxhat; // Calculate the alpha coefficients - Eigen::MatrixXd alpha = Ainv*F; // 16x1 + Eigen::MatrixXd alpha = Ainv.transpose()*F; // 16x1; Watch out for the transpose! std::vector valpha = eigen_to_vec1D(alpha); coeffs[i][j].set(param, valpha); coeffs[i][j].set_valid(); @@ -105,3 +104,84 @@ void CoolProp::BicubicBackend::build_coeffs_ph() std::cout << format("Calculated bicubic coefficients for %d good cells in %g sec.", valid_cell_count, elapsed); } } + +void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double val1, double val2) +{ + // Flush the cached indices (set to large number) + cached_single_phase_i = std::numeric_limits::max(); + cached_single_phase_j = std::numeric_limits::max(); + cached_saturation_iL = std::numeric_limits::max(); + cached_saturation_iV = std::numeric_limits::max(); + + switch(input_pair){ + case HmolarP_INPUTS:{ + _hmolar = val1; _p = val2; + if (!single_phase_logph.native_inputs_are_in_range(_hmolar, _p)){ + // Use the AbstractState instance + using_single_phase_table = false; + if (get_debug_level() > 5){ std::cout << "inputs are not in range"; } + throw ValueError(format("inputs are not in range, hmolar=%Lg, p=%Lg", static_cast(_hmolar), _p)); + } + else{ + using_single_phase_table = true; // Use the table ! + std::size_t iL, iV; + CoolPropDbl hL = 0, hV = 0; + if (pure_saturation.is_inside(_p, iHmolar, _hmolar, iL, iV, hL, hV)){ + using_single_phase_table = false; + _Q = (static_cast(_hmolar)-hL)/(hV-hL); + if(!is_in_closed_range(0.0,1.0,static_cast(_Q))){ + throw ValueError("vapor quality is not in (0,1)"); + } + else{ + cached_saturation_iL = iL; cached_saturation_iV = iV; + } + } + else{ + // Find and cache the indices i, j + selected_table = SELECTED_PH_TABLE; + single_phase_logph.find_native_nearest_good_cell(_hmolar, _p, cached_single_phase_i, cached_single_phase_j); + } + } + break; + } + case HmassP_INPUTS:{ + update(HmolarP_INPUTS, val1 * AS->molar_mass(), val2); // H: [J/kg] * [kg/mol] -> [J/mol] + return; + } + } +} + +/// Use the single_phase table to evaluate an output +double CoolProp::BicubicBackend::evaluate_single_phase(SinglePhaseGriddedTableData &table, std::vector > &coeffs, parameters output, double x, double y, std::size_t i, std::size_t j) +{ + CellCoeffs &cell = coeffs[i][j]; + + if (!cell.valid()){throw ValueError("Cell is invalid");} + + // Get the alpha coefficients + const std::vector &alpha = cell.get(output); + + // Normalized value in the range (0, 1) + double xhat = (x - table.xvec[i])/(table.xvec[i+1] - table.xvec[i]); + double yhat = (y - table.yvec[j])/(table.yvec[j+1] - table.yvec[j]); + + // Calculate the output value desired + double val = 0; + for (std::size_t l = 0; l < 4; ++l) + { + for(std::size_t m = 0; m < 4; ++m) + { + val += alpha[m*4+l]*pow(xhat, static_cast(l))*pow(yhat, static_cast(m)); + } + } + + // Cache the output value calculated + switch(output){ + case iT: _T = val; break; + case iDmolar: _rhomolar = val; break; + case iSmolar: _smolar = val; break; + //case iUmolar: + default: throw ValueError(); + } + return val; +} \ No newline at end of file diff --git a/src/Backends/Tabular/BicubicBackend.h b/src/Backends/Tabular/BicubicBackend.h index 7ed9591c..bb402604 100644 --- a/src/Backends/Tabular/BicubicBackend.h +++ b/src/Backends/Tabular/BicubicBackend.h @@ -110,8 +110,13 @@ class BicubicBackend : public TabularBackend /// Build the \f$a_{i,j}\f$ coefficients for bicubic interpolation void build_coeffs_ph(); void update(CoolProp::input_pairs input_pair, double val1, double val2); - double evaluate_single_phase_phmolar(parameters output, std::size_t i, std::size_t j){return _HUGE;}; - double evaluate_single_phase_pT(parameters output, std::size_t i, std::size_t j){return _HUGE;}; + double evaluate_single_phase(SinglePhaseGriddedTableData &table, std::vector > &coeffs, parameters output, double x, double y, std::size_t i, std::size_t j); + double evaluate_single_phase_phmolar(parameters output, std::size_t i, std::size_t j){ + return evaluate_single_phase(single_phase_logph, coeffs_ph, output, _hmolar, _p, i, j); + }; + double evaluate_single_phase_pT(parameters output, std::size_t i, std::size_t j){ + return evaluate_single_phase(single_phase_logpT, coeffs_pT, output, _T, _p, i, j); + }; }; double do_one(); diff --git a/src/Backends/Tabular/TabularBackends.h b/src/Backends/Tabular/TabularBackends.h index d9406b13..95d2abb1 100644 --- a/src/Backends/Tabular/TabularBackends.h +++ b/src/Backends/Tabular/TabularBackends.h @@ -331,6 +331,8 @@ class SinglePhaseGriddedTableData{ /// Find the nearest cell with lower left coordinate (i,j) where (i,j) is a good node, and so are (i+1,j), (i,j+1), (i+1,j+1) /// This is needed for bicubic interpolation void find_native_nearest_good_cell(double x, double y, std::size_t &i, std::size_t &j){ + bisect_vector(xvec, x, i); + bisect_vector(yvec, y, j); } };