Implement calling a flash routine with guess values provided at the low-level (PQ only for now)

src/Backends/Helmholtz/FlashRoutines.cpp

@@ -443,6 +443,27 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
     }
 }
 
+void FlashRoutines::PQ_flash_with_guesses(HelmholtzEOSMixtureBackend &HEOS, const GuessesStructure &guess)
+{
+	SaturationSolvers::newton_raphson_saturation NR;
+    SaturationSolvers::newton_raphson_saturation_options IO;
+	IO.rhomolar_liq = guess.rhomolar_liq;
+	IO.rhomolar_vap = guess.rhomolar_vap;
+	IO.x = std::vector<CoolPropDbl>(guess.x.begin(), guess.x.end());
+	IO.y = std::vector<CoolPropDbl>(guess.y.begin(), guess.y.end());
+	IO.T = guess.T;
+	IO.p = guess.p;
+	IO.bubble_point = false;
+	IO.imposed_variable = SaturationSolvers::newton_raphson_saturation_options::P_IMPOSED;
+	
+	NR.call(HEOS, IO.y, IO.x, IO);
+
+	// Load the other outputs
+    HEOS._phase = iphase_twophase;
+	HEOS._rhomolar = 1/(HEOS._Q/IO.rhomolar_vap + (1 - HEOS._Q)/IO.rhomolar_liq);
+    HEOS._T = IO.T;
+}
+
 void FlashRoutines::PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS, parameters other, CoolPropDbl value)
 {
     

src/Backends/Helmholtz/FlashRoutines.h

@@ -31,7 +31,11 @@ public:
     /// Flash for given pressure and (molar) quality
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     static void PQ_flash(HelmholtzEOSMixtureBackend &HEOS);
-    
+
+	/// Flash for given pressure and (molar) quality with guess values provided
+    /// @param HEOS The HelmholtzEOSMixtureBackend to be used
+	static void PQ_flash_with_guesses(HelmholtzEOSMixtureBackend &HEOS, const GuessesStructure &guess);
+
     /// Flash for given temperature and (molar) quality
     /// @param HEOS The HelmholtzEOSMixtureBackend to be used
     static void QT_flash(HelmholtzEOSMixtureBackend &HEOS);

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp

@@ -942,6 +942,27 @@ void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double
     
 }
 
+void HelmholtzEOSMixtureBackend::update_with_guesses(CoolProp::input_pairs input_pair, double value1, double value2, const GuessesStructure &guesses)
+{
+	if (get_debug_level() > 10){std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)",__FILE__,__LINE__, input_pair, get_input_pair_short_desc(input_pair).c_str(), value1, value2) << std::endl;}
+    
+    CoolPropDbl ld_value1 = value1, ld_value2 = value2;
+    pre_update(input_pair, ld_value1, ld_value2);
+    value1 = ld_value1; value2 = ld_value2;
+
+    switch(input_pair)
+    {
+        case PQ_INPUTS:
+            _p = value1; _Q = value2; FlashRoutines::PQ_flash_with_guesses(*this, guesses); break;
+        default:
+            throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
+    }
+
+
+    
+    post_update();
+}
+
 void HelmholtzEOSMixtureBackend::post_update()
 {
     // Check the values that must always be set

src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h

@@ -101,6 +101,11 @@ public:
      * @param value2 The second input value
      */
     void update(CoolProp::input_pairs input_pair, double value1, double value2);
+
+	/** \brief Update the state using guess values
+	 * 
+	 */
+	void update_with_guesses(CoolProp::input_pairs input_pair, double Value1, double Value2, const GuessesStructure &guesses);
     
     /** \brief Update all the internal variables for a state by copying from another state
      */

src/Backends/Helmholtz/MixtureParameters.cpp

@@ -194,9 +194,15 @@ void set_mixture_binary_pair_data(const std::string &CAS1, const std::string &CA
         std::vector<Dictionary> &v = mixturebinarypairlibrary.binary_pair_map[CAS];
         try{
             v[0].add_number(key, value);
+            double got = v[0].get_double(key);
+            if (std::abs(got-value) > 1e-10){
+                throw ValueError("Did not set value properly");
+            }
+        }
+        catch(std::exception &e){ 
+            throw ValueError(format("Could not set the parameter [%s] for the binary pair [%s,%s] - for now this is an error; error: %s", 
+                                    key.c_str(), CAS1.c_str(), CAS2.c_str(), e.what()));
         }
-        catch(...){ }
-        throw ValueError(format("Could not match the parameter [%s] for the binary pair [%s,%s] - for now this is an error.", key.c_str(), CAS1.c_str(), CAS2.c_str()));
     }
     else{
         throw ValueError(format("Could not match the binary pair [%s,%s] - for now this is an error.",CAS1.c_str(), CAS2.c_str()));
@@ -361,6 +367,12 @@ void MixtureParameters::set_mixture_parameters(HelmholtzEOSMixtureBackend &HEOS)
                                                                                                         dict_red.get_string("Name1").c_str(),
                                                                                                         dict_red.get_string("Name2").c_str()   ));
             }
+            /*
+            if (i == 0){
+                std::cout << format("betaT %10.9Lg gammaT %10.9Lg betaV %10.9Lg gammaV %10.9Lg %s",  
+                                    beta_T[i][j], gamma_T[i][j], beta_v[i][j], gamma_v[i][j], get_mixture_binary_pair_data(CAS[0],CAS[1],"gammaT").c_str()) << std::endl;
+            }
+            */
 
             // ***************************************************
             //     Departure functions used in excess term

src/Backends/Helmholtz/VLERoutines.cpp

@@ -1156,10 +1156,7 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS
         
         double omega = 1.0;
         if (options.sstype == imposed_p){
-            if (std::abs(change) > 0.05*T){
-                omega = 0.1;
-            }
-            T += omega*change;
+			T += change;
         }
         else if (options.sstype == imposed_T){
             if (std::abs(change) > 0.05*p){
@@ -1308,8 +1305,9 @@ void SaturationSolvers::newton_raphson_saturation::check_Jacobian()
 void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<CoolPropDbl> &z, std::vector<CoolPropDbl> &z_incipient, newton_raphson_saturation_options &IO)
 {
     int iter = 0;
+	bool debug = get_debug_level() > 9 || false;
     
-    if (get_debug_level() > 9){std::cout << " NRsat::call:  p" << IO.p << " T" << IO.T << " dl" << IO.rhomolar_liq << " dv" << IO.rhomolar_vap << std::endl;}
+    if (debug){std::cout << " NRsat::call:  p " << IO.p << " T " << IO.T << " dl " << IO.rhomolar_liq << " dv " << IO.rhomolar_vap << std::endl;}
     
     // Reset all the variables and resize
     pre_call();
@@ -1375,7 +1373,9 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
         else{
             throw ValueError("invalid imposed_variable");
         }
-        //std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg")  << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
+        if(debug){
+			std::cout << format("\t%Lg ", this->error_rms) << T << " " << rhomolar_liq << " " << rhomolar_vap << " v " << vec_to_string(v, "%0.10Lg")  << " x " << vec_to_string(x, "%0.10Lg") << " r " << vec_to_string(r, "%0.10Lg") << std::endl;
+		}
         
         min_rel_change = min_abs_value(err_rel);
         iter++;