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

Derivatives work for both TTSE and BICUBIC, as do cv and cp

Browse Source
This commit is contained in:
Ian Bell
2015-03-21 20:44:43 -06:00
parent 55aade4776
commit 0fc326b886
6 changed files with 130 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -3,6 +3,12 @@
#include "BicubicBackend.h"
#include "MatrixMath.h"
CoolPropDbl pow(CoolPropDbl x, std::size_t n){
return pow(x, static_cast<int>(n));
}
double pow(double x, std::size_t n){
return pow(x, static_cast<int>(n));
}
/// 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] = {
@@ -27,8 +33,8 @@ static Eigen::Matrix<double, 16, 16> Ainv(Ainv_data);
void CoolProp::BicubicBackend::build_coeffs(SinglePhaseGriddedTableData &table, std::vector<std::vector<CellCoeffs> > &coeffs)
{
const bool debug = get_debug_level() > 5 || false;
const int param_count = 5;
parameters param_list[param_count] = {iDmolar, iT, iSmolar, iHmolar, iP}; //iUmolar
const int param_count = 6;
parameters param_list[param_count] = {iDmolar, iT, iSmolar, iHmolar, iP, iUmolar};
std::vector<std::vector<double> > *f = NULL, *fx = NULL, *fy = NULL, *fxy = NULL;
clock_t t1 = clock();
@@ -56,6 +62,9 @@ void CoolProp::BicubicBackend::build_coeffs(SinglePhaseGriddedTableData &table,
break;
case iHmolar:
f = &(table.hmolar); fx = &(table.dhmolardx); fy = &(table.dhmolardy); fxy = &(table.d2hmolardxdy);
break;
case iUmolar:
f = &(table.umolar); fx = &(table.dumolardx); fy = &(table.dumolardy); fxy = &(table.d2umolardxdy);
break;
default:
throw ValueError();
@@ -303,6 +312,7 @@ double CoolProp::BicubicBackend::evaluate_single_phase(SinglePhaseGriddedTableDa
/// Use the single_phase table to evaluate an output
double CoolProp::BicubicBackend::evaluate_single_phase_derivative(SinglePhaseGriddedTableData &table, std::vector<std::vector<CellCoeffs> > &coeffs, parameters output, double x, double y, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny)
{
// Get the cell
CellCoeffs &cell = coeffs[i][j];
@@ -318,6 +328,8 @@ double CoolProp::BicubicBackend::evaluate_single_phase_derivative(SinglePhaseGri
// Calculate the output value desired
double val = 0;
if (Nx == 1 && Ny == 0){
if (output == table.xkey) { return 1.0; }
if (output == table.ykey) { return 0.0; }
for (std::size_t l = 1; l < 4; ++l)
{
for(std::size_t m = 0; m < 4; ++m)
@@ -329,6 +341,8 @@ double CoolProp::BicubicBackend::evaluate_single_phase_derivative(SinglePhaseGri
return val*dxhatdx;
}
else if (Ny == 1 && Nx == 0){
if (output == table.ykey) { return 1.0; }
if (output == table.xkey) { return 0.0; }
for (std::size_t l = 0; l < 4; ++l)
{
for(std::size_t m = 1; m < 4; ++m)

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

@@ -27,6 +27,7 @@ class CellCoeffs{
case iDmolar: return rhomolar;
case iHmolar: return hmolar;
case iSmolar: return smolar;
case iUmolar: return umolar;
default: throw KeyError(format("Invalid key to get() function of CellCoeffs"));
}
};
@@ -38,6 +39,7 @@ class CellCoeffs{
case iDmolar: rhomolar = mat; break;
case iHmolar: hmolar = mat; break;
case iSmolar: smolar = mat; break;
case iUmolar: umolar = mat; break;
default: throw KeyError(format("Invalid key to set() function of CellCoeffs"));
}
};

71
src/Backends/Tabular/TTSEBackend.cpp
View File

@@ -149,7 +149,10 @@ double CoolProp::TTSEBackend::evaluate_single_phase(SinglePhaseGriddedTableData
z = &table.hmolar; dzdx = &table.dhmolardx; dzdy = &table.dhmolardy;
d2zdxdy = &table.d2hmolardxdy; d2zdx2 = &table.d2hmolardx2; d2zdy2 = &table.d2hmolardy2;
break;
//case iUmolar:
case iUmolar:
z = &table.umolar; dzdx = &table.dumolardx; dzdy = &table.dumolardy;
d2zdxdy = &table.d2umolardxdy; d2zdx2 = &table.d2umolardx2; d2zdy2 = &table.d2umolardy2;
break;
case iviscosity:
z = &table.visc; break;
case iconductivity:
@@ -162,10 +165,6 @@ double CoolProp::TTSEBackend::evaluate_single_phase(SinglePhaseGriddedTableData
double deltax = x - table.xvec[i];
double deltay = y - table.yvec[j];
if (output == iconductivity || output == iviscosity){
// Linear interpolation
}
// Calculate the output value desired
double val = (*z)[i][j]+deltax*(*dzdx)[i][j]+deltay*(*dzdy)[i][j]+0.5*deltax*deltax*(*d2zdx2)[i][j]+0.5*deltay*deltay*(*d2zdy2)[i][j]+deltay*deltax*(*d2zdxdy)[i][j];
@@ -181,4 +180,66 @@ double CoolProp::TTSEBackend::evaluate_single_phase(SinglePhaseGriddedTableData
return val;
}
/// Use the single-phase table to evaluate an output
double CoolProp::TTSEBackend::evaluate_single_phase_derivative(SinglePhaseGriddedTableData &table, parameters output, double x, double y, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny)
{
if (Nx == 1 && Ny == 0){
if (output == table.xkey) { return 1.0; }
if (output == table.ykey) { return 0.0; }
}
else if (Ny == 1 && Nx == 0){
if (output == table.ykey) { return 1.0; }
if (output == table.xkey) { return 0.0; }
}
// Define pointers for the matrices to be used;
std::vector<std::vector<double> > *z = NULL, *dzdx = NULL, *dzdy = NULL, *d2zdxdy = NULL, *d2zdy2 = NULL, *d2zdx2 = NULL;
// Connect the pointers based on the output variable desired
switch(output){
case iT:
z = &table.T; dzdx = &table.dTdx; dzdy = &table.dTdy;
d2zdxdy = &table.d2Tdxdy; d2zdx2 = &table.d2Tdx2; d2zdy2 = &table.d2Tdy2;
break;
case iDmolar:
z = &table.rhomolar; dzdx = &table.drhomolardx; dzdy = &table.drhomolardy;
d2zdxdy = &table.d2rhomolardxdy; d2zdx2 = &table.d2rhomolardx2; d2zdy2 = &table.d2rhomolardy2;
break;
case iSmolar:
z = &table.smolar; dzdx = &table.dsmolardx; dzdy = &table.dsmolardy;
d2zdxdy = &table.d2smolardxdy; d2zdx2 = &table.d2smolardx2; d2zdy2 = &table.d2smolardy2;
break;
case iHmolar:
z = &table.hmolar; dzdx = &table.dhmolardx; dzdy = &table.dhmolardy;
d2zdxdy = &table.d2hmolardxdy; d2zdx2 = &table.d2hmolardx2; d2zdy2 = &table.d2hmolardy2;
break;
case iUmolar:
z = &table.umolar; dzdx = &table.dumolardx; dzdy = &table.dumolardy;
d2zdxdy = &table.d2umolardxdy; d2zdx2 = &table.d2umolardx2; d2zdy2 = &table.d2umolardy2;
break;
default:
throw ValueError();
}
// Distances from the node
double deltax = x - table.xvec[i];
double deltay = y - table.yvec[j];
double val;
// Calculate the output value desired
if (Nx == 1 && Ny == 0){
if (output == table.xkey) { return 1.0; }
if (output == table.ykey) { return 0.0; }
val = (*dzdx)[i][j] + deltax*(*d2zdx2)[i][j] + deltay*(*d2zdxdy)[i][j];
}
else if (Ny == 1 && Nx == 0){
if (output == table.ykey) { return 1.0; }
if (output == table.xkey) { return 0.0; }
val = (*dzdy)[i][j] + deltay*(*d2zdy2)[i][j] + deltax*(*d2zdxdy)[i][j];
}
else{
throw NotImplementedError("only first derivatives currently supported");
}
return val;
}
#endif // !defined(NO_TABULAR_BACKENDS)

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

@@ -41,9 +41,7 @@ class TTSEBackend : public TabularBackend
* @param Ny The number of derivatives with respect to y with x held constant
* @return
*/
double evaluate_single_phase_derivative(SinglePhaseGriddedTableData &table, parameters output, double x, double y, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny){
throw NotImplementedError("Not yet implemented");
}
double evaluate_single_phase_derivative(SinglePhaseGriddedTableData &table, parameters output, double x, double y, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny);
double evaluate_single_phase_phmolar_derivative(parameters output, std::size_t i, std::size_t j, std::size_t Nx, std::size_t Ny){
return evaluate_single_phase_derivative(single_phase_logph, output, _hmolar, _p, i, j, Nx, Ny);
};

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

@@ -212,6 +212,7 @@ void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractS
rhomolar[i][j] = AS->rhomolar();
hmolar[i][j] = AS->hmolar();
smolar[i][j] = AS->smolar();
umolar[i][j] = AS->umolar();
// -------------------------
// Transport properties
@@ -237,6 +238,8 @@ void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractS
dhmolardy[i][j] = AS->first_partial_deriv(iHmolar, ykey, xkey);
dsmolardx[i][j] = AS->first_partial_deriv(iSmolar, xkey, ykey);
dsmolardy[i][j] = AS->first_partial_deriv(iSmolar, ykey, xkey);
dumolardx[i][j] = AS->first_partial_deriv(iUmolar, xkey, ykey);
dumolardy[i][j] = AS->first_partial_deriv(iUmolar, ykey, xkey);
// ----------------------------------------
// Second derivatives of state variables
@@ -256,6 +259,9 @@ void CoolProp::SinglePhaseGriddedTableData::build(shared_ptr<CoolProp::AbstractS
d2smolardx2[i][j] = AS->second_partial_deriv(iSmolar, xkey, ykey, xkey, ykey);
d2smolardxdy[i][j] = AS->second_partial_deriv(iSmolar, xkey, ykey, ykey, xkey);
d2smolardy2[i][j] = AS->second_partial_deriv(iSmolar, ykey, xkey, ykey, xkey);
d2umolardx2[i][j] = AS->second_partial_deriv(iUmolar, xkey, ykey, xkey, ykey);
d2umolardxdy[i][j] = AS->second_partial_deriv(iUmolar, xkey, ykey, ykey, xkey);
d2umolardy2[i][j] = AS->second_partial_deriv(iUmolar, ykey, xkey, ykey, xkey);
}
}
}

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

@@ -15,7 +15,7 @@
* See http://stackoverflow.com/a/148610
* See http://stackoverflow.com/questions/147267/easy-way-to-use-variables-of-enum-types-as-string-in-c#202511
*/
#define LIST_OF_MATRICES X(T) X(p) X(rhomolar) X(hmolar) X(smolar) X(dTdx) X(dTdy) X(dpdx) X(dpdy) X(drhomolardx) X(drhomolardy) X(dhmolardx) X(dhmolardy) X(dsmolardx) X(dsmolardy) X(d2Tdx2) X(d2Tdxdy) X(d2Tdy2) X(d2pdx2) X(d2pdxdy) X(d2pdy2) X(d2rhomolardx2) X(d2rhomolardxdy) X(d2rhomolardy2) X(d2hmolardx2) X(d2hmolardxdy) X(d2hmolardy2) X(d2smolardx2) X(d2smolardxdy) X(d2smolardy2) X(visc) X(cond)
#define LIST_OF_MATRICES X(T) X(p) X(rhomolar) X(hmolar) X(smolar) X(umolar) X(dTdx) X(dTdy) X(dpdx) X(dpdy) X(drhomolardx) X(drhomolardy) X(dhmolardx) X(dhmolardy) X(dsmolardx) X(dsmolardy) X(dumolardx) X(dumolardy) X(d2Tdx2) X(d2Tdxdy) X(d2Tdy2) X(d2pdx2) X(d2pdxdy) X(d2pdy2) X(d2rhomolardx2) X(d2rhomolardxdy) X(d2rhomolardy2) X(d2hmolardx2) X(d2hmolardxdy) X(d2hmolardy2) X(d2smolardx2) X(d2smolardxdy) X(d2smolardy2) X(d2umolardx2) X(d2umolardxdy) X(d2umolardy2) X(visc) X(cond)
/** ***MAGIC WARNING***!! X Macros in use
* See http://stackoverflow.com/a/148610
@@ -156,6 +156,12 @@ class PureFluidSaturationTableData{
double hV = CubicInterp(logpV, hmolarV, iV-2, iV-1, iV, iV+1, logp);
double hL = CubicInterp(logpL, hmolarL, iL-2, iL-1, iL, iL+1, logp);
return Q*hV + (1-Q)*hL;
}
case iUmolar:
{
double uV = CubicInterp(logpV, umolarV, iV-2, iV-1, iV, iV+1, logp);
double uL = CubicInterp(logpL, umolarL, iL-2, iL-1, iL, iL+1, logp);
return Q*uV + (1-Q)*uL;
}
case iDmolar:
{
@@ -502,7 +508,15 @@ class TabularBackend : public AbstractState
return calc_first_partial_deriv(iHmolar, iT, iP);
}
else{
throw ValueError("table not selected");
throw ValueError("Two-phase not possible for cpmolar currently");
}
}
CoolPropDbl calc_cvmolar(void){
if (using_single_phase_table){
return calc_first_partial_deriv(iUmolar, iT, iDmolar);
}
else{
throw ValueError("Two-phase not possible for cvmolar currently");
}
}
@@ -534,20 +548,34 @@ class TabularBackend : public AbstractState
}
CoolPropDbl calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant){
if (using_single_phase_table){
CoolPropDbl dOf_dx, dOf_dy, dWrt_dx, dWrt_dy, dConstant_dx, dConstant_dy;
switch(selected_table){
case SELECTED_PH_TABLE: {
if (Of == iHmolar && Wrt == iT && Constant == iP){
double val = 1/(evaluate_single_phase_phmolar_derivative(iT,cached_single_phase_i, cached_single_phase_j,1,0));
return val;
}
else{
throw ValueError("This derivative output not yet supported");
}
dOf_dx = evaluate_single_phase_phmolar_derivative(Of, cached_single_phase_i, cached_single_phase_j,1,0);
dOf_dy = evaluate_single_phase_phmolar_derivative(Of, cached_single_phase_i, cached_single_phase_j,0,1);
dWrt_dx = evaluate_single_phase_phmolar_derivative(Wrt, cached_single_phase_i, cached_single_phase_j,1,0);
dWrt_dy = evaluate_single_phase_phmolar_derivative(Wrt, cached_single_phase_i, cached_single_phase_j,0,1);
dConstant_dx = evaluate_single_phase_phmolar_derivative(Constant, cached_single_phase_i, cached_single_phase_j,1,0);
dConstant_dy = evaluate_single_phase_phmolar_derivative(Constant, cached_single_phase_i, cached_single_phase_j,0,1);
break;
}
case SELECTED_PT_TABLE:{
dOf_dx = evaluate_single_phase_pT_derivative(Of, cached_single_phase_i, cached_single_phase_j,1,0);
dOf_dy = evaluate_single_phase_pT_derivative(Of, cached_single_phase_i, cached_single_phase_j,0,1);
dWrt_dx = evaluate_single_phase_pT_derivative(Wrt, cached_single_phase_i, cached_single_phase_j,1,0);
dWrt_dy = evaluate_single_phase_pT_derivative(Wrt, cached_single_phase_i, cached_single_phase_j,0,1);
dConstant_dx = evaluate_single_phase_pT_derivative(Constant, cached_single_phase_i, cached_single_phase_j,1,0);
dConstant_dy = evaluate_single_phase_pT_derivative(Constant, cached_single_phase_i, cached_single_phase_j,0,1);
break;
}
case SELECTED_PT_TABLE: throw ValueError("No P-T derivatives yet");
case SELECTED_NO_TABLE: throw ValueError("table not selected");
}
return _HUGE; // not needed, will never be hit, just to make compiler happy
return (dOf_dx*dConstant_dy-dOf_dy*dConstant_dx)/(dWrt_dx*dConstant_dy-dWrt_dy*dConstant_dx);
//return 1/(evaluate_single_phase_phmolar_derivative(iT,cached_single_phase_i, cached_single_phase_j,1,0));
}
else{
return pure_saturation.evaluate(iconductivity, _p, _Q, cached_saturation_iL, cached_saturation_iV);