From 6eaea3e642dc41668ba748b4c5ffeed9b833cd70 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Thu, 26 Feb 2015 20:54:30 -0700 Subject: [PATCH] Bicubic coefficients are calculated in all cells Signed-off-by: Ian Bell --- include/Exceptions.h | 9 +++ src/AbstractState.cpp | 8 +++ src/Backends/Tabular/BicubicBackend.cpp | 92 ++++++++++++++++++++++--- src/Backends/Tabular/BicubicBackend.h | 54 ++++++++++++++- src/Backends/Tabular/TTSEBackend.h | 1 + 5 files changed, 153 insertions(+), 11 deletions(-) diff --git a/include/Exceptions.h b/include/Exceptions.h index 81c0c4e7..ce7c38fc 100644 --- a/include/Exceptions.h +++ b/include/Exceptions.h @@ -46,6 +46,15 @@ public: virtual const char* what() const throw(){ return err.c_str(); } }; +class KeyError: public CoolPropBaseError +{ +public: + KeyError() throw() {}; + KeyError(std::string errstring) throw() {err=errstring;}; + ~KeyError() throw() {}; + virtual const char* what() const throw(){ return err.c_str(); } +}; + class AttributeError: public CoolPropBaseError { public: diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 9c0aab1a..163a1c95 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -15,6 +15,7 @@ #include "Backends/Incompressible/IncompressibleBackend.h" #include "Backends/Helmholtz/Fluids/FluidLibrary.h" #include "Backends/Tabular/TTSEBackend.h" +#include "Backends/Tabular/BicubicBackend.h" namespace CoolProp { @@ -51,6 +52,13 @@ AbstractState * AbstractState::factory(const std::string &backend, const std::ve shared_ptr AS(factory(backend.substr(5), fluid_names[0])); return new TTSEBackend(AS); } + else if (backend.find("BICUBIC&") == 0) + { + if (fluid_names.size() != 1){throw ValueError(format("For backend [%s], name vector must be one element long", backend.c_str()));} + // Will throw if there is a problem with this backend + shared_ptr AS(factory(backend.substr(8), fluid_names[0])); + return new BicubicBackend(AS); + } else if (!backend.compare("TREND")) { throw ValueError("TREND backend not yet implemented"); diff --git a/src/Backends/Tabular/BicubicBackend.cpp b/src/Backends/Tabular/BicubicBackend.cpp index ccd8a050..fadf4581 100644 --- a/src/Backends/Tabular/BicubicBackend.cpp +++ b/src/Backends/Tabular/BicubicBackend.cpp @@ -1,5 +1,5 @@ #include "BicubicBackend.h" -#include "Eigen/Core" +#include "MatrixMath.h" /// The inverse of the A matrix for the bicubic interpolation (http://en.wikipedia.org/wiki/Bicubic_interpolation) static const double Ainv_data[16*16] = { @@ -19,17 +19,89 @@ 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); -namespace CoolProp +void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double val1, double val2) { -double do_one(){ - Eigen::Matrix f; - f.setRandom(); - Eigen::Matrix alpha = Ainv*f; - Eigen::Matrix a(alpha.data()); - return a(2,0); -} } +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 + std::vector > *f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL; + + SinglePhaseGriddedTableData &table = single_phase_logph; + std::vector > &coeffs = coeffs_ph; + + clock_t t1 = clock(); + // Resize the coefficient structures + coeffs_ph.resize(table.Nx - 1, std::vector(table.Ny - 1)); + int valid_cell_count = 0; + for (std::size_t k = 0; k < param_count; ++k){ + parameters param = param_list[k]; + if (param == table.xkey || param == table.ykey){continue;} // Skip tables that match either of the input variables + + switch(param){ + case iT: + f = &(table.T); fx = &(table.dTdx); fy = &(table.dTdy); fxy = &(table.d2Tdxdy); + break; + case iP: + 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); + break; + case iSmolar: + 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); + break; + default: + throw ValueError(); + } + for (std::size_t i = 0; i < table.Nx-1; ++i) // -1 since we have one fewer cells than nodes + { + for (std::size_t j = 0; j < table.Ny-1; ++j) // -1 since we have one fewer cells than nodes + { + if (ValidNumber((*f)[i][j]) && ValidNumber((*f)[i+1][j]) && ValidNumber((*f)[i][j+1]) && ValidNumber((*f)[i+1][j+1])){ + + // This will hold the scaled f values for the cell + Eigen::Matrix F; + // The output values (do not require scaling + F(0) = (*f)[i][j]; F(1) = (*f)[i+1][j]; F(2) = (*f)[i][j+1]; F(3) = (*f)[i+1][j+1]; + // Scaling parameter + // d(f)/dxhat = df/dx * dx/dxhat, where xhat = (x-x_i)/(x_{i+1}-x_i) + coeffs[i][j].dx_dxhat = table.xvec[i+1]-table.xvec[i]; + double dx_dxhat = coeffs[i][j].dx_dxhat; + F(4) = (*fx)[i][j]*dx_dxhat; F(5) = (*fx)[i+1][j]*dx_dxhat; + F(6) = (*fx)[i][j+1]*dx_dxhat; F(7) = (*fx)[i+1][j+1]*dx_dxhat; + // Scaling parameter + // d(f)/dyhat = df/dy * dy/dyhat, where yhat = (y-y_j)/(y_{j+1}-y_j) + coeffs[i][j].dy_dyhat = table.yvec[j+1]-table.yvec[j]; + double dy_dyhat = coeffs[i][j].dy_dyhat; + F(8) = (*fy)[i][j]*dy_dyhat; F(9) = (*fy)[i+1][j]*dy_dyhat; + F(10) = (*fy)[i][j+1]*dy_dyhat; F(11) = (*fy)[i+1][j+1]*dy_dyhat; + // Cross derivatives are doubly scaled following the examples above + 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 + std::vector valpha = eigen_to_vec1D(alpha); + coeffs[i][j].set(param, valpha); + coeffs[i][j].set_valid(); + valid_cell_count++; + } + else{ + coeffs[i][j].set_invalid(); + } + } + } + } + double elapsed = (clock() - t1)/((double)CLOCKS_PER_SEC); + if (debug){ + std::cout << format("Calculated bicubic coefficients for %d good cells in %g sec.", valid_cell_count, elapsed); + } +} diff --git a/src/Backends/Tabular/BicubicBackend.h b/src/Backends/Tabular/BicubicBackend.h index 65383362..7ed9591c 100644 --- a/src/Backends/Tabular/BicubicBackend.h +++ b/src/Backends/Tabular/BicubicBackend.h @@ -2,9 +2,52 @@ #define BICUBICBACKEND_H #include "TabularBackends.h" +#include "Exceptions.h" +#include "Eigen/Core" + namespace CoolProp { + +/// This structure holds the coefficients for one cell, the coefficients are stored in matrices +/// and can be obtained by the get() function. +class CellCoeffs{ + private: + bool _valid; + public: + + double dx_dxhat, dy_dyhat; + CellCoeffs(){_valid = false;} + std::vector T, rhomolar, hmolar, p, smolar, umolar; + /// Return a const reference to the desired matrix + const std::vector & get(parameters params){ + switch(params){ + case iT: return T; + case iP: return p; + case iDmolar: return rhomolar; + case iHmolar: return hmolar; + case iSmolar: return smolar; + default: throw KeyError(format("Invalid key to get() function of CellCoeffs")); + } + }; + /// Set one of the matrices in this class + void set(parameters params, const std::vector &mat){ + switch(params){ + case iT: T = mat; break; + case iP: p = mat; break; + case iDmolar: rhomolar = mat; break; + case iHmolar: hmolar = mat; break; + case iSmolar: smolar = mat; break; + default: throw KeyError(format("Invalid key to set() function of CellCoeffs")); + } + }; + /// Returns true if the cell coefficients seem to have been calculated properly + bool valid(){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 + void set_invalid(){_valid = false;}; +}; /** \brief This class implements bicubic interpolation, as very clearly laid out by * the page on wikipedia: http://en.wikipedia.org/wiki/Bicubic_interpolation @@ -55,11 +98,20 @@ A^{-1} = \left[ \begin{array}{*{16}c} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \f] \ */ +typedef std::vector > mat; class BicubicBackend : public TabularBackend { + protected: + std::vector > coeffs_ph, coeffs_pT; public: + /// Instantiator; base class loads or makes tables + BicubicBackend(shared_ptr AS) : TabularBackend (AS){build_coeffs_ph();}; std::string backend_name(void){return "BicubicBackend";} - BicubicBackend(shared_ptr AS) : TabularBackend (AS) {} + /// 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 do_one(); diff --git a/src/Backends/Tabular/TTSEBackend.h b/src/Backends/Tabular/TTSEBackend.h index b4dcb89d..1f3e17fb 100644 --- a/src/Backends/Tabular/TTSEBackend.h +++ b/src/Backends/Tabular/TTSEBackend.h @@ -10,6 +10,7 @@ class TTSEBackend : public TabularBackend { public: std::string backend_name(void){return "TTSEBackend";} + /// Instantiator; base class loads or makes tables TTSEBackend(shared_ptr AS) : TabularBackend (AS) {} 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);