CoolProp/dev/Tickets/60.cpp

#include "../../CoolProp/CoolPropTools.h"
#include "../../CoolProp/IncompBase.h"
#include "../../CoolProp/IncompLiquid.h"
#include "../../CoolProp/IncompSolution.h"

#if defined(__ISWINDOWS__)
#include <Windows.h>
#include "stdafx.h"
#else
#include <sys/time.h>
#include <ctime>
#include <stdint.h>
#endif

#include <vector>
#include <math.h>
#include <string.h>
#include <iostream>
#include <map>
#include <utility>


/** A class to test the different implementations,
 *  exposes more functions than the other ones.
 */
class IncompressibleTest : public IncompressibleClass {
public:
	IncompressibleTest(){};
    ~IncompressibleTest(){};
    // Some functions need top be overwritten!

public:
	double testSi(std::vector<double> const& coefficients, double x){
		return simplePolynomial(coefficients,x);
	}

	double testHo(std::vector<double> const& coefficients, double x){
		return baseHorner(coefficients,x);
	}

	double testSiInt(std::vector<double> const& coefficients, double T){
		return simplePolynomialInt(coefficients,T);
	}

	double testHoInt(std::vector<double> const& coefficients, double T){
		return baseHornerInt(coefficients,T);
	}

	double test2sInt(std::vector<double> const& coefficients, double T){
		return integrateIn2Steps(coefficients,T);
	}

	double testSiInt(std::vector<double> const& coefficients, double T1, double T0){
		return simplePolynomialInt(coefficients,T1,T0);
	}

	double testHoInt(std::vector<double> const& coefficients, double T1, double T0){
		return baseHornerInt(coefficients,T1,T0);
	}

	double test2sInt(std::vector<double> const& coefficients, double T1, double T0){
		return integrateIn2Steps(coefficients,T1,T0);
	}

	double testSiFra(std::vector<double> const& coefficients, double T){
		return simpleFracInt(coefficients,T);
	}

	double testHoFra(std::vector<double> const& coefficients, double T){
		return baseHornerFra(coefficients,T);
	}

	double test2sFra(std::vector<double> const& coefficients, double T){
		return fracIntIn2Steps(coefficients,T);
	}

	double testSiFra(std::vector<double> const& coefficients, double T1, double T0){
		return simpleFracInt(coefficients,T1,T0);
	}

	double testHoFra(std::vector<double> const& coefficients, double T1, double T0){
		return baseHornerFra(coefficients,T1,T0);
	}

	double test2sFra(std::vector<double> const& coefficients, double T1, double T0){
		return fracIntIn2Steps(coefficients,T1,T0);
	}

	// And the same in 2D
	double testSi(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return simplePolynomial(coefficients,x,T);
	}

	double testHo(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return baseHorner(coefficients,x,T);
	}

	double testSiInt(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return simplePolynomialInt(coefficients,x,T);
	}

	double testHoInt(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return baseHornerInt(coefficients,x,T);
	}

	double test2sInt(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return integrateIn2Steps(coefficients,x,T,false);
	}

	double testSiInt(std::vector<std::vector<double> > const& coefficients, double x, double T1, double T0){
		return simplePolynomialInt(coefficients,x,T1,T0);
	}

	double testHoInt(std::vector<std::vector<double> > const& coefficients, double x, double T1, double T0){
		return baseHornerInt(coefficients,x,T1,T0);
	}

	double test2sInt(std::vector<std::vector<double> > const& coefficients, double x, double T1, double T0){
		return integrateIn2Steps(coefficients,x,T1,T0);
	}

	double testSiFra(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return simpleFracInt(coefficients,x,T);
	}

	double testHoFra(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return baseHornerFra(coefficients,x,T);
	}

	double test2sFra(std::vector<std::vector<double> > const& coefficients, double x, double T){
		return fracIntIn2Steps(coefficients,x,T);
	}

	double testSiFra(std::vector<std::vector<double> > const& coefficients, double x, double T1, double T0){
		return simpleFracInt(coefficients,x,T1,T0);
	}

	double testHoFra(std::vector<std::vector<double> > const& coefficients, double x, double T1, double T0){
		return baseHornerFra(coefficients,x,T1,T0);
	}

	double test2sFra(std::vector<std::vector<double> > const& coefficients, double x, double T1, double T0){
		return fracIntIn2Steps(coefficients,x,T1,T0);
	}
};

/* Returns the amount of milliseconds elapsed since the UNIX epoch. Works on both
 * windows and linux.
 * Taken from http://stackoverflow.com/questions/1861294/how-to-calculate-execution-time-of-a-code-snippet-in-c
 * */
uint64_t getTimeValue() {
#if defined(__ISWINDOWS__)
	/* Windows */
	FILETIME ft;
	LARGE_INTEGER li;

	/* Get the amount of 100 nano seconds intervals elapsed since January 1, 1601 (UTC) and copy it
	 * to a LARGE_INTEGER structure. */
	GetSystemTimeAsFileTime(&ft);
	li.LowPart = ft.dwLowDateTime;
	li.HighPart = ft.dwHighDateTime;

	uint64_t ret = li.QuadPart;
	ret -= 116444736000000000LL; /* Convert from file time to UNIX epoch time. */
	ret /= 10; /* From 100 nano seconds (10^-7) to 1 microsecond (10^-6) intervals */

	return ret;
#else
	/* Linux */
	struct timeval tv;

	gettimeofday(&tv, NULL);

	uint64_t ret = tv.tv_usec;
	/* In micro seconds (10^-6), add the seconds (10^0)
	 * after converting them to microseconds (10^-6) */
	ret += (tv.tv_sec * 1000000);

	return ret;
#endif
}

std::vector< std::vector<double> > makeMatrix(std::vector<double> const& coefficients){
	//IncompressibleClass::checkCoefficients(coefficients,18);
	std::vector< std::vector<double> > matrix;
	std::vector<double> tmpVector;

	tmpVector.clear();
	tmpVector.push_back(coefficients[0]);
	tmpVector.push_back(coefficients[6]);
	tmpVector.push_back(coefficients[11]);
	tmpVector.push_back(coefficients[15]);
	matrix.push_back(tmpVector);

	tmpVector.clear();
	tmpVector.push_back(coefficients[1]);
	tmpVector.push_back(coefficients[7]);
	tmpVector.push_back(coefficients[12]);
	tmpVector.push_back(coefficients[16]);
	matrix.push_back(tmpVector);

	tmpVector.clear();
	tmpVector.push_back(coefficients[2]);
	tmpVector.push_back(coefficients[8]);
	tmpVector.push_back(coefficients[13]);
	tmpVector.push_back(coefficients[17]);
	matrix.push_back(tmpVector);

	tmpVector.clear();
	tmpVector.push_back(coefficients[3]);
	tmpVector.push_back(coefficients[9]);
	tmpVector.push_back(coefficients[14]);
	tmpVector.push_back(0.0);
	matrix.push_back(tmpVector);

	tmpVector.clear();
	tmpVector.push_back(coefficients[4]);
	tmpVector.push_back(coefficients[10]);
	tmpVector.push_back(0.0);
	tmpVector.push_back(0.0);
	matrix.push_back(tmpVector);

	tmpVector.clear();
	tmpVector.push_back(coefficients[5]);
	tmpVector.push_back(0.0);
	tmpVector.push_back(0.0);
	tmpVector.push_back(0.0);
	matrix.push_back(tmpVector);

	tmpVector.clear();
	return matrix;
}

std::map<std::string, double> testObject(IncompressibleTest* fluid, std::vector<std::vector<double> > coeffs, int exponent, bool print){

	std::map<std::string, double> results;

	uint64_t runs = 10;
	for (int i=1; i<exponent; i++){
		runs *= 10;
	}

	if (print) std::cout << "Runs: " << runs << std::endl;

	std::vector<double> simCoeffs;
	if (coeffs.size()>0) {
		simCoeffs = coeffs[1];
	} else {
		results["error"] = -_HUGE;
		return results;
	}

	if (print) {
		std::cout << "simCoeffs: " << simCoeffs[0];
		for (int i=1; i < simCoeffs.size(); i++) {
			std::cout << ", " << simCoeffs[i];
		}
		std::cout << std::endl;
	}


	uint64_t start;
	uint64_t end;
	double time;

	double Tin = 280;
	double T0  = Tin-10;
	double xin = 0.25;

	// Testing the 1D functions first
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
	  	fluid->testSi(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D std sim simple"] = time;
	fluid->setDebug(print);
	fluid->testSi(simCoeffs, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHo(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D std sim Horner"] = time;
	fluid->setDebug(print);
	fluid->testHo(simCoeffs, Tin);

	// Testing the 1D integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiInt(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D int ind simple"] = time;
	fluid->setDebug(print);
	fluid->testSiInt(simCoeffs, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoInt(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D int ind Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoInt(simCoeffs, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sInt(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D int ind 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sInt(simCoeffs, Tin);

	// Testing the 1D definite integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiInt(simCoeffs, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D int def simple"] = time;
	fluid->setDebug(print);
	fluid->testSiInt(simCoeffs, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoInt(simCoeffs, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D int def Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoInt(simCoeffs, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sInt(simCoeffs, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D int def 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sInt(simCoeffs, Tin, T0);

	// Testing the 1D fraction integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiFra(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D fra ind simple"] = time;
	fluid->setDebug(print);
	fluid->testSiFra(simCoeffs, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoFra(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D fra ind Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoFra(simCoeffs, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sFra(simCoeffs, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D fra ind 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sFra(simCoeffs, Tin);

	// Testing the 1D definite fraction integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiFra(simCoeffs, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D fra def simple"] = time;
	fluid->setDebug(print);
	fluid->testSiFra(simCoeffs, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoFra(simCoeffs, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D fra def Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoFra(simCoeffs, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sFra(simCoeffs, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["1D fra def 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sFra(simCoeffs, Tin, T0);

	// Testing the 2D functions
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSi(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D std sim simple"] = time;
	fluid->setDebug(print);
	fluid->testSi(coeffs, xin, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHo(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D std sim Horner"] = time;
	fluid->setDebug(print);
	fluid->testHo(coeffs, xin, Tin);

	// Testing the 2D integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiInt(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D int ind simple"] = time;
	fluid->setDebug(print);
	fluid->testSiInt(coeffs, xin, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoInt(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D int ind Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoInt(coeffs, xin, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sInt(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D int ind 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sInt(coeffs, xin, Tin);

	// Testing the 2D definite integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiInt(coeffs, xin, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D int def simple"] = time;
	fluid->setDebug(print);
	fluid->testSiInt(coeffs, xin, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoInt(coeffs, xin, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D int def Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoInt(coeffs, xin, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sInt(coeffs, xin, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D int def 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sInt(coeffs, xin, Tin, T0);

	// Testing the 2D fraction integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiFra(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D fra ind simple"] = time;
	fluid->setDebug(print);
	fluid->testSiFra(coeffs, xin, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoFra(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D fra ind Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoFra(coeffs, xin, Tin);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sFra(coeffs, xin, Tin+i/runs);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D fra ind 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sFra(coeffs, xin, Tin);

	// Testing the 2D definite fraction integrators
	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testSiFra(coeffs, xin, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D fra def simple"] = time;
	fluid->setDebug(print);
	fluid->testSiFra(coeffs, xin, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->testHoFra(coeffs, xin, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D fra def Horner"] = time;
	fluid->setDebug(print);
	fluid->testHoFra(coeffs, xin, Tin, T0);

	fluid->setDebug(false);
	start = getTimeValue();
	for (int i=0; i<runs; i++){
		fluid->test2sFra(coeffs, xin, Tin+i/runs, T0);
	}
	end = getTimeValue();
	time = (end-start)*1e3/runs;
	results["2D fra def 2Steps"] = time;
	fluid->setDebug(print);
	fluid->test2sFra(coeffs, xin, Tin, T0);

	return results;
}

int main(int argc, const char* argv[]) {



	std::vector<double> tmpVector;
	std::vector<std::vector<double> > coeffs;
	int exponent = 6;

	IncompressibleTest* fluid = new IncompressibleTest();

	tmpVector.clear();
	tmpVector.push_back( 1081.6353100);
	tmpVector.push_back(-2.4559523700);
	tmpVector.push_back( 0.0058152057);
	tmpVector.push_back(-7.500013E-05);
	tmpVector.push_back(-7.575759E-07);
	tmpVector.push_back( 1.666671E-07);
	tmpVector.push_back(-5.6609963900);
	tmpVector.push_back( 0.1002726190);
	tmpVector.push_back(-0.0004797330);
	tmpVector.push_back( 1.333333E-06);
	tmpVector.push_back( 3.636364E-08);
	tmpVector.push_back(-0.0852857143);
	tmpVector.push_back( 0.0007904762);
	tmpVector.push_back( 1.428571E-06);
	tmpVector.push_back( 6.666668E-07);
	tmpVector.push_back(-0.0037650794);
	tmpVector.push_back( 3.333333E-05);
	tmpVector.push_back( 6.984127E-07);
	coeffs.clear();
	coeffs = makeMatrix(tmpVector);


	std::map<std::string, double> time = testObject(fluid, coeffs, exponent, true);

	char buff[100];
	std::map<std::string, double>::iterator iter;
	for (iter = time.begin(); iter != time.end(); ++iter) {
	    sprintf(buff, "%8.3f", iter->second);
	    std::cout << "Time consumption in " << iter->first << ": " << buff << " ns per call for 1e" << exponent << " calls." << std::endl;
	}

	//std::cout << "Time consumption liquid: "   << time_liq << " µs per call from 1e" << exponent << " calls." << std::endl;

	//std::cout << "Time consumption liquid: "   << time_liq << " µs per call from 1e" << exponent << " calls." << std::endl;
	//std::cout << "Time consumption solution: " << time_sol << " µs per call from 1e" << exponent << " calls." << std::endl;


//	SecCoolSolution* obj = new MethanolSolution();
//	double x = 0.25;
//	double T = 5.0 + 273.15;
//	double p = 3e5;
//
//	obj->testInputs(T + 00, p, x);
//	obj->testInputs(T + 05, p, x);
//	obj->testInputs(T + 10, p, x);
//	obj->testInputs(T + 15, p, x);

}








//double result = coefficients[0] * log(T);
//if (coefficients.size() > 1) {
//	for (unsigned int i=1; i<coefficients.size(); i++){
//		result += 1/(i) * coefficients[i] * pow(T,(int)(i));
//	}
//}
//
//
//double result = coefficients[0] * log(T);
//if (coefficients.size() > 1) {
//	std::vector<double> newCoeffs(coefficients.begin() + 1, coefficients.end());
//	result += polyint(newCoeffs,T);
//}