github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-23 03:00:17 -04:00
Code Issues Packages Projects Releases Wiki Activity

Bicubic coefficients are calculated in all cells

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2015-02-26 20:54:30 -07:00
parent 52a123437c
commit 6eaea3e642
5 changed files with 153 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
include/Exceptions.h
View File

@@ -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:

8
src/AbstractState.cpp
View File

@@ -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<AbstractState> 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<AbstractState> AS(factory(backend.substr(8), fluid_names[0]));
return new BicubicBackend(AS);
}
else if (!backend.compare("TREND"))
{
throw ValueError("TREND backend not yet implemented");

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

@@ -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<double, 16, 16> Ainv(Ainv_data);
namespace CoolProp
void CoolProp::BicubicBackend::update(CoolProp::input_pairs input_pair, double val1, double val2)
{
double do_one(){
Eigen::Matrix<double, 16, 1> f;
f.setRandom();
Eigen::Matrix<double, 16, 1> alpha = Ainv*f;
Eigen::Matrix<double, 4, 4> 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<std::vector<double> > *f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL;
SinglePhaseGriddedTableData &table = single_phase_logph;
std::vector<std::vector<CellCoeffs> > &coeffs = coeffs_ph;
clock_t t1 = clock();
// Resize the coefficient structures
coeffs_ph.resize(table.Nx - 1, std::vector<CellCoeffs>(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<double, 16, 1> 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<double> 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);
}
}

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

@@ -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<double> T, rhomolar, hmolar, p, smolar, umolar;
/// Return a const reference to the desired matrix
const std::vector<double> & 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<double> &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<std::vector<double> > mat;
class BicubicBackend : public TabularBackend
{
protected:
std::vector<std::vector<CellCoeffs> > coeffs_ph, coeffs_pT;
public:
/// Instantiator; base class loads or makes tables
BicubicBackend(shared_ptr<CoolProp::AbstractState> AS) : TabularBackend (AS){build_coeffs_ph();};
std::string backend_name(void){return "BicubicBackend";}
BicubicBackend(shared_ptr<CoolProp::AbstractState> 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();

1
src/Backends/Tabular/TTSEBackend.h
View File

@@ -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<CoolProp::AbstractState> 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);