github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-08 21:05:14 -05:00
Code Issues Packages Projects Releases Wiki Activity

Added build directory to ignore, more on polymath and incompressibles

Browse Source
This commit is contained in:
jowr
2014-05-18 17:23:15 +02:00
parent 8987cf3954
commit 70fdcf2613
16 changed files with 540 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
src/Backends/Helmholtz/Fluids/Incompressible.cpp
View File

@@ -7,15 +7,124 @@
#include "Incompressible.h"
namespace CoolProp {
Incompressible::Incompressible() {
// TODO Auto-generated constructor stub
std::vector<std::vector<double> > Incompressible::changeAxis(const std::vector<double> &input){
std::vector<std::vector<double> > output;
std::vector<double> tmp;
size_t sizeX = input.size();
for (size_t i = 0; i < sizeX; ++i){
tmp.clear();
tmp.push_back(input[i]);
output.push_back(tmp);
}
return output;
}
Incompressible::~Incompressible() {
// TODO Auto-generated destructor stub
/* All functions need T and p as input. Might not
* be necessary, but gives a clearer structure.
*/
/// Density as a function of temperature, pressure and composition.
double Incompressible::rho(double T_K, double p) {
return poly.polyval(cRho, getxInput(x), getTInput(T_K));
}
/// Heat capacities as a function of temperature, pressure and composition.
double Incompressible::c(double T_K, double p) {
return poly.polyval(cHeat, getxInput(x), getTInput(T_K));
}
/// Enthalpy as a function of temperature, pressure and composition.
double Incompressible::h(double T_K, double p) {
return h_u(T_K, p);
}
/// Entropy as a function of temperature, pressure and composition.
double Incompressible::s(double T_K, double p) {
return poly.polyfracintcentral(cHeat, getxInput(x), T_K, Tbase)
- poly.polyfracintcentral(cHeat, getxInput(x), Tref, Tbase);
}
/// Viscosity as a function of temperature, pressure and composition.
double Incompressible::visc(double T_K, double p) {
return expo.expval(cVisc, getxInput(x), getTInput(T_K), 2) / 1e3;
}
/// Thermal conductivity as a function of temperature, pressure and composition.
double Incompressible::cond(double T_K, double p) {
return poly.polyval(cCond, getxInput(x), getTInput(T_K));
}
/// Internal energy as a function of temperature, pressure and composition.
double Incompressible::u(double T_K, double p) {
return poly.polyint(cHeat, getxInput(x), getTInput(T_K))
- poly.polyint(cHeat, getxInput(x), getTInput(Tref));
}
/// Saturation pressure as a function of temperature and composition.
double Incompressible::psat(double T_K ){throw NotImplementedError("Psat is not available");};
/// Freezing temperature as a function of pressure and composition.
double Incompressible::Tfreeze( double p){throw NotImplementedError("Tfreeze is not available");};
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not
* cause errors. */
bool Incompressible::checkT(double T_K, double p){
if( Tmin < 0. ) {
throw ValueError("Please specify the minimum temperature.");
} else if( Tmax < 0.) {
throw ValueError("Please specify the maximum temperature.");
} else if ( (Tmin>T_K) || (T_K>Tmax) ) {
throw ValueError(format("Your temperature %f is not between %f and %f.",T_K,Tmin,Tmax));
} else if (T_K < Tfreeze(p)) {
throw ValueError(format("Your temperature %f is below the freezing point of %f.",T_K,Tfreeze(p)));
} else {
return true;
}
return false;
}
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool Incompressible::checkP(double T_K, double p) {
double ps = psat(T_K);
if (p<ps) {
throw ValueError(format("Equations are valid for solution phase only: %f < %f. ",p,ps));
} else {
return true;
}
}
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool Incompressible::checkX(double x){
if( xmin < 0. ) {
throw ValueError("Please specify the minimum concentration.");
} else if( xmax < 0.) {
throw ValueError("Please specify the maximum concentration.");
} else if ( (xmin>x) || (x>xmax) ) {
throw ValueError(format("Your composition %f is not between %f and %f.",x,xmin,xmax));
} else {
return true;
}
return false;
}
/// Check validity of temperature, pressure and composition input.
bool Incompressible::checkTPX(double T, double p, double x) {
return (checkT(T,p) && checkP(T,p) && checkX(x));
}
} /* namespace CoolProp */

157
src/Backends/Helmholtz/Fluids/Incompressible.h
View File

@@ -8,14 +8,169 @@
#ifndef INCOMPRESSIBLE_H_
#define INCOMPRESSIBLE_H_
//#include "AbstractFluid.h"
#include "PolyMath.h"
namespace CoolProp {
class Incompressible{
protected:
std::string name;
std::string description;
std::string reference;
double Tmin, Tmax;
double xmin, xmax;
double TminPsat;
double xref, Tref;
double xbase, Tbase;
double x;
std::vector< std::vector<double> > cRho;
std::vector< std::vector<double> > cHeat;
std::vector< std::vector<double> > cVisc;
std::vector< std::vector<double> > cCond;
std::vector< std::vector<double> > cPsat;
std::vector<double> cTfreeze;
std::vector<std::vector<double> > changeAxis(const std::vector<double> &input);
BasePolynomial poly;
PolynomialSolver solver;
BaseExponential expo;
public:
Incompressible();
virtual ~Incompressible();
std::string getName() const {return name;}
std::string get_name() const {return getName();}// For backwards-compatibility.
std::string getDescription() const {return description;}
std::string getReference() const {return reference;}
double getTmax() const {return Tmax;}
double getTmin() const {return Tmin;}
double getxmax() const {return xmax;}
double getxmin() const {return xmin;}
double getTminPsat() const {return TminPsat;}
double getTref() const {return Tref;}
double getxref() const {return xref;}
double getTbase() const {return Tbase;}
double getxbase() const {return xbase;}
double getx() const {return x;}
void setName(std::string name) {this->name = name;}
void setDescription(std::string description) {this->description = description;}
void setReference(std::string reference) {this->reference = reference;}
void setTmax(double Tmax) {this->Tmax = Tmax;}
void setTmin(double Tmin) {this->Tmin = Tmin;}
void setxmax(double xmax) {this->xmax = xmax;}
void setxmin(double xmin) {this->xmin = xmin;}
void setTminPsat(double TminPsat) {this->TminPsat = TminPsat;}
void setTref(double Tref) {this->Tref = Tref;}
void setxref(double xref) {this->xref = xref;}
void setTbase(double Tbase) {this->Tbase = Tbase;}
void setxbase(double xbase) {this->xbase = xbase;}
void setx(double x) {this->x = x;}
// Setters for the coefficients
void setcRho(std::vector<std::vector<double> > cRho){this->cRho = cRho;}
void setcHeat(std::vector<std::vector<double> > cHeat){this->cHeat = cHeat;}
void setcVisc(std::vector<std::vector<double> > cVisc){this->cVisc = cVisc;}
void setcCond(std::vector<std::vector<double> > cCond){this->cCond = cCond;}
void setcPsat(std::vector<std::vector<double> > cPsat){this->cPsat = cPsat;}
void setcTfreeze(std::vector<double> cTfreeze){this->cTfreeze = cTfreeze;}
void setcRho(std::vector<double> cRho){this->cRho = changeAxis(cRho);}
void setcHeat(std::vector<double> cHeat){this->cHeat = changeAxis(cHeat);}
void setcVisc(std::vector<double> cVisc){this->cVisc = changeAxis(cVisc);}
void setcCond(std::vector<double> cCond){this->cCond = changeAxis(cCond);}
void setcPsat(std::vector<double> cPsat){this->cPsat = changeAxis(cPsat);}
double getTInput(double curTValue){return curTValue-Tbase;}
double getxInput(double curxValue){return (curxValue-xbase)*100.0;}
/* All functions need T and p as input. Might not
* be necessary, but gives a clearer structure.
*/
/// Density as a function of temperature, pressure and composition.
virtual double rho (double T_K, double p);
/// Heat capacities as a function of temperature, pressure and composition.
virtual double c (double T_K, double p);
virtual double cp (double T_K, double p){return c(T_K,p);};
virtual double cv (double T_K, double p){return c(T_K,p);};
/// Entropy as a function of temperature, pressure and composition.
virtual double s (double T_K, double p);
/// Internal energy as a function of temperature, pressure and composition.
virtual double u (double T_K, double p);
/// Enthalpy as a function of temperature, pressure and composition.
virtual double h (double T_K, double p);
/// Viscosity as a function of temperature, pressure and composition.
virtual double visc(double T_K, double p);
/// Thermal conductivity as a function of temperature, pressure and composition.
virtual double cond(double T_K, double p);
/// Saturation pressure as a function of temperature and composition.
virtual double psat(double T_K );
/// Freezing temperature as a function of pressure and composition.
virtual double Tfreeze( double p);
protected:
/* Define internal energy and enthalpy as functions of the
* other properties to provide data in case there are no
* coefficients.
*/
/// Enthalpy from u, p and rho.
/** Calculate enthalpy as a function of temperature and
* pressure employing functions for internal energy and
* density. Provides consistent formulations. */
double h_u(double T_K, double p) {
return u(T_K,p)+p/rho(T_K,p);
};
/// Internal energy from h, p and rho.
/** Calculate internal energy as a function of temperature
* and pressure employing functions for enthalpy and
* density. Provides consistent formulations. */
double u_h(double T_K, double p) {
return h(T_K,p)-p/rho(T_K,p);
};
/*
* Some more functions to provide a single implementation
* of important routines.
* We start with the check functions that can validate input
* in terms of pressure p, temperature T and composition x.
*/
/// Check validity of temperature input.
/** Compares the given temperature T to the result of a
* freezing point calculation. This is not necessarily
* defined for all fluids, default values do not cause errors. */
bool checkT(double T_K, double p);
/// Check validity of pressure input.
/** Compares the given pressure p to the saturation pressure at
* temperature T and throws and exception if p is lower than
* the saturation conditions.
* The default value for psat is -1 yielding true if psat
* is not redefined in the subclass.
* */
bool checkP(double T_K, double p);
/// Check validity of composition input.
/** Compares the given composition x to a stored minimum and
* maximum value. Enforces the redefinition of xmin and
* xmax since the default values cause an error. */
bool checkX(double x);
/// Check validity of temperature, pressure and composition input.
bool checkTPX(double T, double p, double x);
};
} /* namespace CoolProp */

98
src/Backends/Incompressible/IncompressibleBackend.cpp Normal file
View File

@@ -0,0 +1,98 @@
#if defined(_MSC_VER)
#define _CRTDBG_MAP_ALLOC
#define _CRT_SECURE_NO_WARNINGS
#include <crtdbg.h>
#include <sys/stat.h>
#else
#include <sys/stat.h>
#endif
#include <string>
//#include "CoolProp.h"
#include "IncompressibleBackend.h"
#include "../Fluids/FluidLibrary.h"
#include "Solvers.h"
#include "MatrixMath.h"
namespace CoolProp {
IncompressibleBackend::IncompressibleBackend(const std::string &fluid_name) {
throw CoolProp::NotImplementedError("Mixture-style constructor is not implemented yet for incompressible fluids");
}
IncompressibleBackend::IncompressibleBackend(const std::vector<std::string> &component_names) {
throw CoolProp::NotImplementedError("Mixture-style constructor is not implemented yet for incompressible fluids");
}
bool IncompressibleBackend::using_mole_fractions(){return false;}
void IncompressibleBackend::update(long input_pair, double value1, double value2) {
switch (input_pair) {
case PT_INPUTS: {
break;
}
// case DmassP_INPUTS: {
//
// }
// break;
// }
// case HmassP_INPUTS: {
// // Call again, but this time with molar units
// // H: [J/kg] * [kg/mol] -> [J/mol]
// update(HmolarP_INPUTS, value1 * (double) _molar_mass, value2);
// return;
// }
// case PUmass_INPUTS: {
// // Call again, but this time with molar units
// // U: [J/kg] * [kg/mol] -> [J/mol]
// update(PUmolar_INPUTS, value1, value2 * (double) _molar_mass);
// return;
// }
// case PSmass_INPUTS: {
// // Call again, but this time with molar units
// // U: [J/kg] * [kg/mol] -> [J/mol]
// update(PUmolar_INPUTS, value1, value2 * (double) _molar_mass);
// return;
// }
default: {
throw ValueError(
format("This pair of inputs [%s] is not yet supported",
get_input_pair_short_desc(input_pair).c_str()));
}
}
}
/// Set the mole fractions
/**
@param mole_fractions The vector of mole fractions of the components
*/
void IncompressibleBackend::set_mole_fractions(const std::vector<long double> &mole_fractions) {
throw CoolProp::NotImplementedError("Cannot set mole fractions for incompressible fluid");
}
/// Set the mass fractions
/**
@param mass_fractions The vector of mass fractions of the components
*/
void IncompressibleBackend::set_mass_fractions(const std::vector<long double> &mass_fractions) {
this->mass_fractions = mass_fractions;
}
/// Check if the mole fractions have been set, etc.
void IncompressibleBackend::check_status() {
throw CoolProp::NotImplementedError("Cannot check status for incompressible fluid");
}
/// Get the viscosity [Pa-s]
long double IncompressibleBackend::calc_viscosity(void){
throw NotImplementedError();
}
/// Get the thermal conductivity [W/m/K] (based on the temperature and density in the state class)
long double IncompressibleBackend::calc_conductivity(void){
throw NotImplementedError();
}
}

23
src/Backends/Incompressible/IncompressibleBackend.h
View File

@@ -20,44 +20,45 @@ public:
virtual ~IncompressibleBackend(){};
/// The instantiator
/// @param fluid_names The vector of strings of the fluid components, without file ending
IncompressibleBackend(const std::string & fluid_names){ };
/// @param fluid_name the string with the fluid name
IncompressibleBackend(const std::string &fluid_name);
/// The instantiator
/// @param fluid_names The vector of strings of the fluid components, without file ending
IncompressibleBackend(const std::vector<std::string> &component_names);
// Incompressible backend uses mole fractions
bool using_mole_fractions(){return false;}
bool using_mole_fractions();
/// Updating function for incompressible fluid
/**
In this function we take a pair of thermodynamic states, those defined in the input_pairs
enumeration and update all the internal variables that we can. REFPROP calculates
a lot of other state variables each time you use a flash routine so we cache all the
outputs that we can, which saves on computational time.
enumeration and update all the internal variables that we can.
@param input_pair Integer key from CoolProp::input_pairs to the two inputs that will be passed to the function
@param value1 First input value
@param value2 Second input value
*/
void update(long input_pair, double value1, double value2){throw CoolProp::NotImplementedError("Update not implemented yet for incompressible fluids");};
void update(long input_pair, double value1, double value2);
/// Set the mole fractions
/**
@param mole_fractions The vector of mole fractions of the components
*/
void set_mole_fractions(const std::vector<long double> &mole_fractions){throw CoolProp::NotImplementedError("Cannot set mole fractions for incompressible fluid");};
void set_mole_fractions(const std::vector<long double> &mole_fractions);
/// Set the mass fractions
/**
@param mass_fractions The vector of mass fractions of the components
*/
void set_mass_fractions(const std::vector<long double> &mass_fractions){this->mass_fractions = mass_fractions;};
void set_mass_fractions(const std::vector<long double> &mass_fractions);
/// Check if the mole fractions have been set, etc.
void check_status();
/// Get the viscosity [Pa-s]
long double calc_viscosity(void){throw NotImplementedError();};
long double calc_viscosity(void);
/// Get the thermal conductivity [W/m/K] (based on the temperature and density in the state class)
long double calc_conductivity(void){throw NotImplementedError();};
long double calc_conductivity(void);
};
} /* namespace CoolProp */

154
src/PolyMath.cpp
View File

@@ -656,6 +656,55 @@ double PolyIntResidual::deriv(double x){
return polyRes;
}
double PolyFracIntResidual::call(double x){
double polyRes = -1;
switch (this->dim) {
case i1D:
polyRes = this->poly.polyfracint(this->coefficients[0], x);
break;
case i2D:
polyRes = this->poly.polyfracint(this->coefficients, this->firstDim, x);
break;
default:
throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
}
return polyRes - this->output;
}
double PolyFracIntResidual::deriv(double x){
double polyRes = -1;
switch (this->dim) {
case i1D:
polyRes = this->poly.polyfracval(this->coefficients[0], x);
break;
case i2D:
polyRes = this->poly.polyfracval(this->coefficients, this->firstDim, x);
break;
default:
throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
}
return polyRes;
}
double PolyFracIntCentralResidual::call(double x){
double polyRes = -1;
switch (this->dim) {
case i1D:
polyRes = this->poly.polyfracintcentral(this->coefficients[0], x, this->baseVal);
break;
case i2D:
polyRes = this->poly.polyfracintcentral(this->coefficients, this->firstDim, x, this->baseVal);
break;
default:
throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
}
return polyRes - this->output;
}
double PolyFracIntCentralResidual::deriv(double x){
throw CoolProp::NotImplementedError("Derivative of a polynomial frac int is not defined.");
}
double PolyDerResidual::call(double x){
double polyRes = -1;
switch (this->dim) {
@@ -678,7 +727,7 @@ double PolyDerResidual::deriv(double x){
/** Implements the same public functions as the
/** Implements the same public functions as the BasePolynomial
* but solves the polynomial for the given value
* instead of evaluating it.
* TODO: This class does not check for bijective
@@ -689,7 +738,7 @@ PolynomialSolver::PolynomialSolver(){
this->DEBUG = false;
this->macheps = DBL_EPSILON;
this->tol = DBL_EPSILON*1e3;
this->maxiter = 100;
this->maxiter = 50;
}
/** Everything related to the normal polynomials goes in this
@@ -781,7 +830,8 @@ double PolynomialSolver::polyfracval(const std::vector< std::vector<double> > &c
/// @param coefficients vector containing the ordered coefficients
/// @param y double value that represents the current output
double PolynomialSolver::polyfracint(const std::vector<double> &coefficients, double y){
throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function
PolyFracIntResidual residual = PolyFracIntResidual(coefficients, y);
return this->solve(residual);
}
/// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable
@@ -789,7 +839,8 @@ double PolynomialSolver::polyfracint(const std::vector<double> &coefficients, do
/// @param x double value that represents the current input in the 1st dimension
/// @param z double value that represents the current output
double PolynomialSolver::polyfracint(const std::vector< std::vector<double> > &coefficients, double x, double z){
throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function
PolyFracIntResidual residual = PolyFracIntResidual(coefficients, x, z);
return this->solve(residual);
}
/// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable
@@ -797,7 +848,8 @@ double PolynomialSolver::polyfracint(const std::vector< std::vector<double> > &c
/// @param y double value that represents the current output
/// @param xbase central x-value for fitted function
double PolynomialSolver::polyfracintcentral(const std::vector<double> &coefficients, double y, double xbase){
throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function
PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, y, xbase);
return this->solve(residual);
}
/// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable
@@ -806,7 +858,8 @@ double PolynomialSolver::polyfracintcentral(const std::vector<double> &coefficie
/// @param z double value that represents the current output
/// @param ybase central y-value for fitted function
double PolynomialSolver::polyfracintcentral(const std::vector< std::vector<double> > &coefficients, double x, double z, double ybase){
throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function
PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, x, z, ybase);
return this->solve(residual);
}
@@ -1024,6 +1077,30 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]")
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyint(cHeat, T);
solver.setGuess(T+100);
val2 = solver.polyint(cHeat, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
// val1 = poly.polyder(cHeat, T);
// solver.setGuess(T+100);
// val2 = solver.polyder(cHeat, val1);
// CAPTURE(T);
// CAPTURE(val1);
// CAPTURE(val2);
// CHECK(fabs(T-val2) < 1e-1);
//
// val1 = poly.polyfracint(cHeat, T);
// solver.setGuess(T+100);
// val2 = solver.polyfracint(cHeat, val1);
// CAPTURE(T);
// CAPTURE(val1);
// CAPTURE(val2);
// CHECK(fabs(T-val2) < 1e-1);
}
}
SECTION("Solve1DBrent") {
@@ -1035,6 +1112,39 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]")
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyint(cHeat, T);
solver.setLimits(T-300,T+300);
val2 = solver.polyint(cHeat, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyder(cHeat, T);
solver.setLimits(T-300,T+300);
val2 = solver.polyder(cHeat, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyfracint(cHeat, T);
solver.setLimits(T-100,T+100);
val2 = solver.polyfracint(cHeat, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyfracintcentral(cHeat, T, 250.0);
solver.setLimits(T-100,T+100);
val2 = solver.polyfracintcentral(cHeat, val1, 250.0);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
}
}
@@ -1094,6 +1204,38 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]")
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyint(cHeat2D, xDim1, T);
solver.setLimits(T-300,T+300);
val2 = solver.polyint(cHeat2D, xDim1, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyder(cHeat2D, xDim1, T);
solver.setLimits(T-300,T+300);
val2 = solver.polyder(cHeat2D, xDim1, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyfracint(cHeat2D, xDim1, T);
solver.setLimits(T-100,T+100);
val2 = solver.polyfracint(cHeat2D, xDim1, val1);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
val1 = poly.polyfracintcentral(cHeat2D, xDim1, T, 250);
solver.setLimits(T-100,T+100);
val2 = solver.polyfracintcentral(cHeat2D, xDim1, val1, 250);
CAPTURE(T);
CAPTURE(val1);
CAPTURE(val2);
CHECK(fabs(T-val2) < 1e-1);
}
}