diff --git a/src/HumidAirProp.cpp b/src/HumidAirProp.cpp index 157650b7..2e23de71 100644 --- a/src/HumidAirProp.cpp +++ b/src/HumidAirProp.cpp @@ -13,12 +13,14 @@ #include "crossplatform_shared_ptr.h" #include "Exceptions.h" +#include // std::next_permutation #include #include "math.h" #include "time.h" #include "stdio.h" #include #include +#include /// This is a stub overload to help with all the strcmp calls below and avoid needing to rewrite all of them std::size_t strcmp(const std::string &s, const std::string e){ @@ -166,7 +168,8 @@ static double Brent_HAProps_T(const std::string &OutputName, const std::string & ~BrentSolverResids(){}; double call(double T){ - return HAPropsSI(OutputName,"T",T,Input1Name,Input1,Input2Name,Input2)-TargetVal; + double val = HAPropsSI(OutputName, "T", T, Input1Name, Input1, Input2Name, Input2); + return val - TargetVal; } }; @@ -198,8 +201,19 @@ static double Brent_HAProps_T(const std::string &OutputName, const std::string & T_min_valid = ValidNumber(r_min); } } - T = CoolProp::Brent(BSR, T_min, T_max, 1e-7, 1e-4, 50, errstr); - + // We will do a secant call if the values at T_min and T_max have the same sign + if (r_min*r_max > 0){ + if (std::abs(r_min) < std::abs(r_max)){ + T = CoolProp::Secant(BSR, T_min, 0.01*T_min, 1e-7, 50, errstr); + } + else{ + T = CoolProp::Secant(BSR, T_max, -0.01*T_max, 1e-7, 50, errstr); + } + } + else{ + T = CoolProp::Brent(BSR, T_min, T_max, 1e-7, 1e-4, 50, errstr); + } + return T; } static double Secant_Tdb_at_saturated_W(double psi_w, double p, double T_guess) @@ -215,17 +229,14 @@ static double Secant_Tdb_at_saturated_W(double psi_w, double p, double T_guess) double call(double T){ double p_ws; - if (T>=273.16) - { + if (T>=273.16){ // Saturation pressure [Pa] Water->update(CoolProp::QT_INPUTS, 0, T); p_ws= Water->keyed_output(CoolProp::iP); } - else - { + else{ // Sublimation pressure [Pa] p_ws=psub_Ice(T); - } double f = f_factor(T, p); double pp_water_calc = f*p_ws; @@ -238,7 +249,7 @@ static double Secant_Tdb_at_saturated_W(double psi_w, double p, double T_guess) BrentSolverResids Resids(psi_w, p); std::string errstr; - T = CoolProp::Secant(Resids, 150, 100, 1e-7, 100, errstr); + T = CoolProp::Brent(Resids, 150, 350, 1e-16, 1e-7, 100, errstr); return T; } @@ -1154,7 +1165,6 @@ double MoleFractionWater(double T, double p, int HumInput, double InVal) { // Sublimation pressure [Pa] p_ws=psub_Ice(T); - } // Enhancement Factor [-] f=f_factor(T,p); @@ -1276,8 +1286,8 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d { // Add a check to make sure that Air and Water fluid states have been properly instantiated check_fluid_instantiation(); - - int In1Type, In2Type, In3Type,iT,iW,iTdp,iRH,ip,Type1,Type2; + int iT, iW, iTdp, iRH, ip; + givens In1Type, In2Type, In3Type, Type1, Type2, OutputType; double vals[3],p,T,RH,W,Tdp,psi_w,M_ha,v_bar,h_bar,s_bar,MainInputValue,SecondaryInputValue,T_guess; double Value1,Value2,W_guess; std::string MainInputName, SecondaryInputName, Name1, Name2; @@ -1285,11 +1295,17 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d vals[0]=Input1; vals[1]=Input2; vals[2]=Input3; + + OutputType=Name2Type(OutputName.c_str()); // First figure out what kind of inputs you have, convert names to Macro expansions In1Type=Name2Type(Input1Name.c_str()); In2Type=Name2Type(Input2Name.c_str()); In3Type=Name2Type(Input3Name.c_str()); + + if (OutputType == In1Type){return Input1;} + if (OutputType == In2Type){return Input2;} + if (OutputType == In3Type){return Input3;} // Pressure must be included ip=TypeMatch(GIVEN_P,Input1Name,Input2Name,Input3Name); @@ -1331,17 +1347,17 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d // Temperature and pressure are known, figure out which variable holds the other value if (In1Type!=GIVEN_T && In1Type!=GIVEN_P) { - strcpy(SecondaryInputName,Input1Name); + strcpy(SecondaryInputName, Input1Name); SecondaryInputValue=Input1; } else if (In2Type!=GIVEN_T && In2Type!=GIVEN_P) { - strcpy(SecondaryInputName,Input2Name); + strcpy(SecondaryInputName, Input2Name); SecondaryInputValue=Input2; } else if (In3Type!=GIVEN_T && In3Type!=GIVEN_P) { - strcpy(SecondaryInputName,Input3Name); + strcpy(SecondaryInputName, Input3Name); SecondaryInputValue=Input3; } else{ @@ -1397,11 +1413,9 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d if (Type1 == GIVEN_HUMRAT) { // MainInput is the one that you are using in the call to HAProps - MainInputValue=Value1; - strcpy(MainInputName,Name1); + MainInputValue=Value1; strcpy(MainInputName,Name1); // SecondaryInput is the one that you are trying to match - SecondaryInputValue=Value2; - strcpy(SecondaryInputName,Name2); + SecondaryInputValue=Value2; strcpy(SecondaryInputName,Name2); } else if (Type2 == GIVEN_HUMRAT) { @@ -1428,11 +1442,11 @@ double HAPropsSI(const std::string &OutputName, const std::string &Input1Name, d } else { - printf("Sorry, but currently at least one of the variables as an input to HAPropsSI() must be temperature, relative humidity, humidity ratio, or dewpoint\n Eventually will add a 2-D NR solver to find T and psi_w simultaneously, but not included now\n"); + CoolProp::ValueError("Sorry, but currently at least one of the variables as an input to HAPropsSI() must be temperature, relative humidity, humidity ratio, or dewpoint\n Eventually will add a 2-D NR solver to find T and psi_w simultaneously, but not included now\n"); return _HUGE; } - double T_min = 210; + double T_min = 200; double T_max = 450; if (Name2Type(MainInputName) == GIVEN_RH){ @@ -2007,13 +2021,88 @@ TEST_CASE_METHOD(HAPropsConsistencyFixture, "ASHRAE RP1485 Tables", "[RP1485]") } } } - TEST_CASE("Assorted tests","[HAPropsSI]") { CHECK(ValidNumber(HumidAir::HAPropsSI("T", "H", 267769, "P", 104300, "W", 0.0))); CHECK(ValidNumber(HumidAir::HAPropsSI("T", "B", 252.84, "W", 5.097e-4, "P", 101325))); CHECK(ValidNumber(HumidAir::HAPropsSI("T", "B",290, "R", 1, "P", 101325))); } +// a predicate implemented as a function: +bool is_not_a_pair (const std::set &item) { return item.size() != 2; } + +const int number_of_inputs = 6; +std::string inputs[number_of_inputs] = {"W","D","B","R","T","V"};//,"H","S"}; + +class ConsistencyTestData +{ +public: + bool is_built; + std::vector data; + std::list > inputs_list; + ConsistencyTestData(){ + is_built = false; + }; + void build(){ + if (is_built){return;} + std::vector indices(number_of_inputs); + for (std::size_t i = 0; i < number_of_inputs; ++i){ indices[i] = i;} + // Generate a powerset of all the permutations of all lengths of inputs + std::set indices_set(indices.begin(), indices.end()); + std::set > inputs_powerset = powerset(indices_set); + inputs_list = std::list >(inputs_powerset.begin(), inputs_powerset.end()); + inputs_list.remove_if(is_not_a_pair); + + const int NT = 20, NW = 20; + double p = 101325; + for (double T = 210; T < 350; T += (350-210)/(NT-1)) + { + double Wsat = HumidAir::HAPropsSI("W", "T", T, "P", p, "R", 1.0); + for (double W = 1e-10; W < Wsat; W += (Wsat-1e-10)/(NW-1)){ + Dictionary vals; + // Calculate all the values using T, W + for (int i = 0; i < number_of_inputs; ++i){ + double v = HumidAir::HAPropsSI(inputs[i], "T", T, "P", p, "W", W); + vals.add_number(inputs[i], v); + } + data.push_back(vals); + std::cout << format("T %g W %g\n",T,W); + } + } + is_built = true; + }; +} consistency_data; + +TEST_CASE("HAPropsSI", "[HAPropsSI]") +{ + consistency_data.build(); + double p = 101325; + for (std::size_t i = 0; i < consistency_data.data.size(); ++i) + { + for (std::list >::iterator iter = consistency_data.inputs_list.begin(); iter != consistency_data.inputs_list.end(); ++iter) + { + std::vector pair(iter->begin(), iter->end()); + std::string i0 = inputs[pair[0]], i1 = inputs[pair[1]]; + double v0 = consistency_data.data[i].get_double(i0), v1 = consistency_data.data[i].get_double(i1); + std::ostringstream ss2; + ss2 << "Inputs: \"" << i0 << "\"," << v0 << ",\"" << i1 << "\"," << v1; + SECTION(ss2.str(), ""){ + + double T = consistency_data.data[i].get_double("T"); + double W = consistency_data.data[i].get_double("W"); + double Wcalc = HumidAir::HAPropsSI("W", i0, v0, i1, v1, "P", p); + double Tcalc = HumidAir::HAPropsSI("T", i0, v0, i1, v1, "P", p); + std::string err = CoolProp::get_global_param_string("errstring"); + CAPTURE(T); + CAPTURE(W); + CAPTURE(Tcalc); + CAPTURE(Wcalc); + CAPTURE(err); + CHECK(std::abs(Tcalc - T) < 1e-1); + CHECK(std::abs((Wcalc - W)/W) < 1e-3); + } + } + } +} #endif /* CATCH_ENABLED */