From 5ff28fc302372d6efbd96e174fa3b899285d177f Mon Sep 17 00:00:00 2001 From: Jorrit Wronski Date: Wed, 7 Jan 2015 12:16:00 +0100 Subject: [PATCH] Clears the cached transport properties and fixes #394 --- src/AbstractState.cpp | 29 ++++---- src/CoolProp.cpp | 159 ++++++++++++++++++++++-------------------- 2 files changed, 102 insertions(+), 86 deletions(-) diff --git a/src/AbstractState.cpp b/src/AbstractState.cpp index 09ec08fc..71cacf6d 100644 --- a/src/AbstractState.cpp +++ b/src/AbstractState.cpp @@ -121,7 +121,7 @@ bool AbstractState::clear() { this->_gibbsmolar.clear(); this->_logp.clear(); this->_logrhomolar.clear(); - + ///// Smoothing values //this->rhospline = -_HUGE; //this->dsplinedp = -_HUGE; @@ -154,6 +154,11 @@ bool AbstractState::clear() { this->_d2alphar_dDelta_dTau_lim.clear(); this->_d3alphar_dDelta2_dTau_lim.clear(); + /// Transport properties + this->_viscosity.clear(); + this->_conductivity.clear(); + this->_surface_tension.clear(); + return true; } double AbstractState::trivial_keyed_output(int key) @@ -419,7 +424,7 @@ void get_dT_drho(AbstractState &AS, parameters index, long double &dT, long doub R = AS.gas_constant(), delta = rho/rhor, tau = Tr/T; - + switch (index) { case iT: @@ -497,9 +502,9 @@ void get_dT_drho_second_derivatives(AbstractState &AS, int index, long double &d delta = rho/rhor, tau = Tr/T; - // Here we use T and rho as independent variables since derivations are already done by Thorade, 2013, + // Here we use T and rho as independent variables since derivations are already done by Thorade, 2013, // Partial derivatives of thermodynamic state propertiesfor dynamic simulation, DOI 10.1007/s12665-013-2394-z - + switch (index) { case iT: @@ -584,7 +589,7 @@ long double AbstractState::calc_first_partial_deriv(parameters Of, parameters Wr } long double AbstractState::calc_second_partial_deriv(parameters Of1, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2) { - long double dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, + long double dOf1_dT, dOf1_drho, dWrt1_dT, dWrt1_drho, dConstant1_dT, dConstant1_drho, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2, dWrt2_dT, dWrt2_drho, dConstant2_dT, dConstant2_drho, N, D, dNdrho__T, dDdrho__T, dNdT__rho, dDdT__rho, dderiv1_drho, dderiv1_dT, second; @@ -596,34 +601,34 @@ long double AbstractState::calc_second_partial_deriv(parameters Of1, parameters get_dT_drho_second_derivatives(*this, Of1, d2Of1_dT2, d2Of1_drhodT, d2Of1_drho2); get_dT_drho_second_derivatives(*this, Wrt1, d2Wrt1_dT2, d2Wrt1_drhodT, d2Wrt1_drho2); get_dT_drho_second_derivatives(*this, Constant1, d2Constant1_dT2, d2Constant1_drhodT, d2Constant1_drho2); - + // First derivatives of terms involved in the second derivative get_dT_drho(*this, Wrt2, dWrt2_dT, dWrt2_drho); get_dT_drho(*this, Constant2, dConstant2_dT, dConstant2_drho); - + // Numerator and denominator of first partial derivative term N = dOf1_dT*dConstant1_drho - dOf1_drho*dConstant1_dT; D = dWrt1_dT*dConstant1_drho - dWrt1_drho*dConstant1_dT; - + // Derivatives of the numerator and denominator of the first partial derivative term with respect to rho, T held constant // They are of similar form, with Of1 and Wrt1 swapped dNdrho__T = dOf1_dT*d2Constant1_drho2 + d2Of1_drhodT*dConstant1_drho - dOf1_drho*d2Constant1_drhodT - d2Of1_drho2*dConstant1_dT; dDdrho__T = dWrt1_dT*d2Constant1_drho2 + d2Wrt1_drhodT*dConstant1_drho - dWrt1_drho*d2Constant1_drhodT - d2Wrt1_drho2*dConstant1_dT; - + // Derivatives of the numerator and denominator of the first partial derivative term with respect to T, rho held constant // They are of similar form, with Of1 and Wrt1 swapped dNdT__rho = dOf1_dT*d2Constant1_drhodT + d2Of1_dT2*dConstant1_drho - dOf1_drho*d2Constant1_dT2 - d2Of1_drhodT*dConstant1_dT; dDdT__rho = dWrt1_dT*d2Constant1_drhodT + d2Wrt1_dT2*dConstant1_drho - dWrt1_drho*d2Constant1_dT2 - d2Wrt1_drhodT*dConstant1_dT; - + // First partial of first derivative term with respect to T dderiv1_drho = (D*dNdrho__T - N*dDdrho__T)/pow(D, 2); - + // First partial of first derivative term with respect to rho dderiv1_dT = (D*dNdT__rho - N*dDdT__rho)/pow(D, 2); // Complete second derivative second = (dderiv1_dT*dConstant2_drho - dderiv1_drho*dConstant2_dT)/(dWrt2_dT*dConstant2_drho - dWrt2_drho*dConstant2_dT); - + return second; } // // ---------------------------------------- diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp index 570648da..0168f159 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -77,7 +77,7 @@ void extract_backend(const std::string &fluid_string, std::string &backend, std: { std::size_t i; std::string _fluid_string = fluid_string; - // For backwards compatibility reasons, if "REFPROP-" or "REFPROP-MIX:" start + // For backwards compatibility reasons, if "REFPROP-" or "REFPROP-MIX:" start // the fluid_string, replace them with "REFPROP::" if (_fluid_string.find("REFPROP-MIX:") == 0) { @@ -177,21 +177,21 @@ std::string extract_fractions(const std::string &fluid_string, std::vector &fluid_names, - const std::vector &z, +void _PropsSI_initialize(const std::string &backend, + const std::vector &fluid_names, + const std::vector &z, shared_ptr &State){ - + if (fluid_names.empty()){throw ValueError("fluid_names cannot be empty");} - + std::vector fractions(1, 1.0); // Default to one component, unity fraction const std::vector *fractions_ptr = NULL; // Pointer to the array to be used; - + if (fluid_names.size() > 1){ // Set the pointer - we are going to use the supplied fractions; they must be provided fractions_ptr = &z; @@ -220,7 +220,7 @@ void _PropsSI_initialize(const std::string &backend, State.reset(AbstractState::factory(backend, fluid_names)); } } - + // Set the fraction for the state if (State->using_mole_fractions()){ // If a predefined mixture or a pure fluid, the fractions will already be set @@ -268,25 +268,36 @@ struct output_parameter{ }; }; -void _PropsSI_outputs(shared_ptr &State, - std::vector output_parameters, - CoolProp::input_pairs input_pair, - const std::vector &in1, +void _PropsSI_outputs(shared_ptr &State, + std::vector output_parameters, + CoolProp::input_pairs input_pair, + const std::vector &in1, const std::vector &in2, std::vector > &IO){ - + // Check the inputs if (in1.size() != in2.size()){ throw ValueError(format("lengths of in1 [%d] and in2 [%d] are not the same", in1.size(), in2.size()));} bool one_input_one_output = (in1.size() == 1 && in2.size() == 1 && output_parameters.size() == 1); - + + if (get_debug_level() > 100) + { + std::cout << format("%s (%d): input pair = %d ",__FILE__,__LINE__, input_pair) << std::endl; + std::cout << format("%s (%d): in1 = %s ",__FILE__,__LINE__, vec_to_string(in1).c_str()) << std::endl; + std::cout << format("%s (%d): in2 = %s ",__FILE__,__LINE__, vec_to_string(in2).c_str()) << std::endl; + } + // Resize the output matrix std::size_t N1 = std::max(static_cast(1), in1.size()); std::size_t N2 = std::max(static_cast(1), output_parameters.size()); IO.resize(N1, std::vector(N2, _HUGE)); - + // Throw an error if at the end, there were no successes bool success = false; - + + if (get_debug_level() > 100) + { + std::cout << format("%s (%d): Iterating over %d input value pairs.",__FILE__,__LINE__,IO.size()) << std::endl; + } // Iterate over the state variable inputs for (std::size_t i = 0; i < IO.size(); ++i){ try{ @@ -301,7 +312,7 @@ void _PropsSI_outputs(shared_ptr &State, for (std::size_t j = 0; j < IO[i].size(); ++j){ IO[i][j] = _HUGE; } continue; } - + for (std::size_t j = 0; j < IO[i].size(); ++j){ try{ output_parameter &output = output_parameters[j]; @@ -328,13 +339,13 @@ void _PropsSI_outputs(shared_ptr &State, if (success == false) { IO.clear(); throw ValueError(format("No outputs were able to be calculated"));} } -void _PropsSImulti(const std::vector &Outputs, - const std::string &Name1, +void _PropsSImulti(const std::vector &Outputs, + const std::string &Name1, const std::vector &Prop1, - const std::string &Name2, - const std::vector &Prop2, - const std::string &backend, - const std::vector &fluids, + const std::string &Name2, + const std::vector &Prop2, + const std::string &backend, + const std::vector &fluids, const std::vector &fractions, std::vector > &IO) { @@ -342,7 +353,7 @@ void _PropsSImulti(const std::vector &Outputs, CoolProp::input_pairs input_pair; std::vector output_parameters; std::vector v1, v2; - + try{ // Initialize the State class _PropsSI_initialize(backend, fluids, fractions, State); @@ -362,9 +373,9 @@ void _PropsSImulti(const std::vector &Outputs, } catch (std::exception &e){ // Input parameter parsing failed. Stop - throw ValueError(format("Input pair parsing failed for Name1: \"%s\", Name2: \"%s\"; err: %s", Name1.c_str(), Name2.c_str(), e.what())); + throw ValueError(format("Input pair parsing failed for Name1: \"%s\", Name2: \"%s\"; err: %s", Name1.c_str(), Name2.c_str(), e.what())); } - + try{ output_parameters = output_parameter::get_output_parameters(Outputs); } @@ -376,34 +387,34 @@ void _PropsSImulti(const std::vector &Outputs, // Calculate the output(s). In the case of a failure, all values will be filled with _HUGE _PropsSI_outputs(State, output_parameters, input_pair, v1, v2, IO); } - -std::vector > PropsSImulti(const std::vector &Outputs, - const std::string &Name1, + +std::vector > PropsSImulti(const std::vector &Outputs, + const std::string &Name1, const std::vector &Prop1, - const std::string &Name2, - const std::vector &Prop2, - const std::string &backend, - const std::vector &fluids, + const std::string &Name2, + const std::vector &Prop2, + const std::string &backend, + const std::vector &fluids, const std::vector &fractions) { std::vector > IO; - + #if !defined(NO_ERROR_CATCHING) try{ #endif - + // Call the subfunction that can bubble errors _PropsSImulti(Outputs, Name1, Prop1, Name2, Prop2, backend, fluids, fractions, IO); - + // Return the value(s) return IO; - + #if !defined(NO_ERROR_CATCHING) } catch(const std::exception& e){ - set_error_string(e.what()); + set_error_string(e.what()); #if defined (PROPSSI_ERROR_STDOUT) - std::cout << e.what() << std::endl; + std::cout << e.what() << std::endl; #endif if (get_debug_level() > 1){std::cout << e.what() << std::endl;} } @@ -418,12 +429,12 @@ double PropsSI(const std::string &Output, const std::string &Name1, double Prop1 #if !defined(NO_ERROR_CATCHING) try{ #endif - + // BEGIN OF TRY // Here is the real code that is inside the try block - - - + + + extract_backend(Ref, backend, fluid); if (has_fractions_in_string(fluid)){ extract_fractions(fluid, fractions); @@ -432,49 +443,49 @@ double PropsSI(const std::string &Output, const std::string &Name1, double Prop1 _PropsSImulti(strsplit(Output,'&'), Name1, std::vector(1, Prop1), Name2, std::vector(1, Prop2), backend, std::vector(1, fluid), fractions, IO); if (IO.empty()){ throw ValueError(get_global_param_string("errstring").c_str()); } if (IO.size()!= 1 || IO[0].size() != 1){ throw ValueError(format("output should be 1x1; error was %s", get_global_param_string("errstring").c_str())); } - + double val = IO[0][0]; - + if (get_debug_level() > 1){ std::cout << format("_PropsSI will return %g",val) << std::endl; } return val; // END OF TRY #if !defined(NO_ERROR_CATCHING) } catch(const std::exception& e){ - set_error_string(e.what() + format(" : PropsSI(\"%s\",\"%s\",%0.10g,\"%s\",%0.10g,\"%s\")",Output.c_str(),Name1.c_str(), Prop1, Name2.c_str(), Prop2, Ref.c_str())); + set_error_string(e.what() + format(" : PropsSI(\"%s\",\"%s\",%0.10g,\"%s\",%0.10g,\"%s\")",Output.c_str(),Name1.c_str(), Prop1, Name2.c_str(), Prop2, Ref.c_str())); #if defined (PROPSSI_ERROR_STDOUT) - std::cout << e.what() << std::endl; + std::cout << e.what() << std::endl; #endif if (get_debug_level() > 1){std::cout << e.what() << std::endl;} - return _HUGE; + return _HUGE; } catch(...){ - return _HUGE; + return _HUGE; } #endif } #if defined(ENABLE_CATCH) TEST_CASE("Check inputs to PropsSI","[PropsSI]") { - SECTION("Single state, single output"){ + SECTION("Single state, single output"){ CHECK(ValidNumber(CoolProp::PropsSI("T","P",101325,"Q",0,"Water"))); }; - SECTION("Single state, single output, pure incompressible"){ + SECTION("Single state, single output, pure incompressible"){ CHECK(ValidNumber(CoolProp::PropsSI("D","P",101325,"T",300,"INCOMP::DowQ"))); }; - SECTION("Bad input pair"){ + SECTION("Bad input pair"){ CHECK(!ValidNumber(CoolProp::PropsSI("D","Q",0,"Q",0,"Water"))); }; - SECTION("Single state, single output, 40% incompressible"){ + SECTION("Single state, single output, 40% incompressible"){ CHECK(ValidNumber(CoolProp::PropsSI("D","P",101325,"T",300,"INCOMP::MEG[0.40]"))); }; - SECTION("Single state, single output, predefined CoolProp mixture"){ + SECTION("Single state, single output, predefined CoolProp mixture"){ CHECK(ValidNumber(CoolProp::PropsSI("T","Q",1,"P",3e6,"HEOS::R125[0.7]&R32[0.3]"))); }; - SECTION("Single state, single output"){ + SECTION("Single state, single output"){ CHECK(ValidNumber(CoolProp::PropsSI("T","P",101325,"Q",0,"HEOS::Water"))); }; - SECTION("Single state, single output, predefined mixture"){ + SECTION("Single state, single output, predefined mixture"){ CHECK(ValidNumber(CoolProp::PropsSI("T","P",101325,"Q",0,"R410A.mix"))); }; SECTION("Predefined mixture"){ @@ -587,7 +598,7 @@ double Props1SI(const std::string &FluidName, const std::string &Output) // They are backwards, swap std::swap(_Output, _FluidName); } - + // First input is the fluid, second input is the input parameter double val1 = PropsSI(_Output, "", 0, "", 0, _FluidName); if (!ValidNumber(val1)){ @@ -601,19 +612,19 @@ double Props1SI(const std::string &FluidName, const std::string &Output) #if defined(ENABLE_CATCH) TEST_CASE("Check inputs to Props1SI","[Props1SI],[PropsSI]") { - SECTION("Good fluid, good parameter"){ + SECTION("Good fluid, good parameter"){ CHECK(ValidNumber(CoolProp::Props1SI("Tcrit","Water"))); }; - SECTION("Good fluid, good parameter"){ + SECTION("Good fluid, good parameter"){ CHECK(ValidNumber(CoolProp::PropsSI("Tcrit","",0,"",0,"Water"))); }; - SECTION("Good fluid, good parameter, inverted"){ + SECTION("Good fluid, good parameter, inverted"){ CHECK(ValidNumber(CoolProp::Props1SI("Water","Tcrit"))); }; - SECTION("Good fluid, bad parameter"){ + SECTION("Good fluid, bad parameter"){ CHECK(!ValidNumber(CoolProp::Props1SI("Water","????????????"))); }; - SECTION("Bad fluid, good parameter"){ + SECTION("Bad fluid, good parameter"){ CHECK(!ValidNumber(CoolProp::Props1SI("?????","Tcrit"))); }; }; @@ -637,14 +648,14 @@ bool is_valid_fluid_string(std::string &input_fluid_string) } } double saturation_ancillary(const std::string &fluid_name, const std::string &output, int Q, const std::string &input, double value){ - + // Generate the state instance std::vector names(1, fluid_name); shared_ptr HEOS(new CoolProp::HelmholtzEOSMixtureBackend(names)); - + parameters iInput = get_parameter_index(input); parameters iOutput = get_parameter_index(output); - + return HEOS->saturation_ancillary(iOutput, Q, iInput, value); } void set_reference_stateS(std::string Ref, std::string reference_state) @@ -709,7 +720,7 @@ void set_reference_stateD(std::string Ref, double T, double rhomolar, double h0, shared_ptr HEOS; std::vector _comps(1, Ref); HEOS.reset(new CoolProp::HelmholtzEOSMixtureBackend(_comps)); - + HEOS->update(DmolarT_INPUTS, rhomolar, T); // Get current values for the enthalpy and entropy @@ -779,8 +790,8 @@ TEST_CASE("Check inputs to get_global_param_string","[get_global_param_string]") std::ostringstream ss3c; for (int i = 0; i INCOMP(new CoolProp::IncompressibleBackend(fluid)); - + if (!ParamName.compare("long_name")){ return INCOMP->calc_name(); } @@ -805,12 +816,12 @@ std::string get_fluid_param_string(std::string FluidName, std::string ParamName) catch(std::exception &e){ throw ValueError(format("CoolProp error: %s", e.what())); } catch(...){ throw ValueError("CoolProp error: Indeterminate error"); } } - + try{ std::vector comps(1, FluidName); shared_ptr HEOS(new CoolProp::HelmholtzEOSMixtureBackend(comps)); CoolProp::CoolPropFluid *fluid = HEOS->get_components()[0]; - + if (!ParamName.compare("aliases")){ return strjoin(fluid->aliases, ", "); } @@ -844,7 +855,7 @@ TEST_CASE("Check inputs to get_fluid_param_string", "[get_fluid_param_string]") std::ostringstream ss3c; for (int i = 0; i < num_good_inputs; ++i){ ss3c << "Test for" << good_inputs[i]; - SECTION(ss3c.str(), ""){ + SECTION(ss3c.str(), ""){ CHECK_NOTHROW(CoolProp::get_fluid_param_string("Water", good_inputs[i])); }; } @@ -886,7 +897,7 @@ std::string PhaseSI(const std::string &Name1, double Prop1, const std::string &N std::size_t Phase_int = static_cast(Phase_double); return phase_lookup_string(static_cast(Phase_int)); } - + /* std::string PhaseSI(const std::string &Name1, double Prop1, const std::string &Name2, double Prop2, const std::string &FluidName, const std::vector &z) {