Merge pull request #29521 from GiudGiud/PR_fp_tab_vh

Add more options to tabulated FPs with (v,h)

modules/fluid_properties/include/fluidproperties/SinglePhaseFluidProperties.h

@@ -438,12 +438,15 @@ class SinglePhaseFluidProperties : public FluidProperties
 
   /**
    * Newton's method may be used to convert between variable sets
-   * _tolerance, _T_initial_guess, and _p_initial_guess are the parameters for these
-   * iterative solves
    */
+  /// Relative tolerance of the solves
   const Real _tolerance;
+  /// Initial guess for temperature (or temperature used to compute the initial guess)
   const Real _T_initial_guess;
+  /// Initial guess for pressure (or pressure used to compute the initial guess)
   const Real _p_initial_guess;
+  /// Maximum number of iterations for the variable conversion newton solves
+  const unsigned int _max_newton_its;
 
 private:
   void unimplementedDerivativeMethod(const std::string & property_function_name) const

modules/fluid_properties/include/fluidproperties/TabulatedFluidProperties.h

@@ -226,7 +226,15 @@ class TabulatedFluidProperties : public SinglePhaseFluidProperties
   /// - if user-specified, use _v_min/max and _e_min/max
   /// - if reading a (v,e) interpolation, the bounds of that range
   /// - if a _fp exist find the min/max v/e from T_min/max and p_min/max
+  void createVGridVector();
   void createVEGridVectors();
+  /// Create (or reset) the grid vectors for the specific volume and enthalpy interpolation
+  /// The order of priority for determining the range boundaries in v and h:
+  /// - if user-specified, use _v_min/max and _e_min/max
+  /// - if a _fp exist find the min/max v/e from T_min/max and p_min/max
+  /// - if reading a (p,T) tabulation, the bounds of the enthalpy grid
+  /// - if reading a (v,e) tabulation, the bounds of the enthalpy grid
+  void createVHGridVectors();
 
   /// Standardized error message for missing interpolation
   void missingVEInterpolationError(const std::string & function_name) const;
@@ -336,6 +344,8 @@ class TabulatedFluidProperties : public SinglePhaseFluidProperties
   bool _log_space_v;
   /// log-space the internal energy interpolation grid axis instead of linear
   bool _log_space_e;
+  /// log-space the enthalpy interpolation grid axis instead of linear
+  bool _log_space_h;
 
   /// User-selected out-of-bounds behavior
   MooseEnum _OOBBehavior;
@@ -345,6 +355,7 @@ class TabulatedFluidProperties : public SinglePhaseFluidProperties
     Ignore,
     Throw,
     DeclareInvalid,
+    WarnInvalid,
     SetToClosestBound
   };
 

modules/fluid_properties/src/fluidproperties/NaKFluidProperties.C

@@ -69,8 +69,13 @@ NaKFluidProperties::T_from_p_rho(Real pressure, Real density) const
   // NOTE we could also invert analytically the third degree polynomial, see Cardan's method
   auto lambda = [&](Real p, Real current_T, Real & new_rho, Real & drho_dp, Real & drho_dT)
   { rho_from_p_T(p, current_T, new_rho, drho_dp, drho_dT); };
-  Real T = FluidPropertiesUtils::NewtonSolve(
-               pressure, density, _T_initial_guess, _tolerance, lambda, name() + "::T_from_p_rho")
+  Real T = FluidPropertiesUtils::NewtonSolve(pressure,
+                                             density,
+                                             _T_initial_guess,
+                                             _tolerance,
+                                             lambda,
+                                             name() + "::T_from_p_rho",
+                                             _max_newton_its)
                .first;
   // check for nans
   if (std::isnan(T))

modules/fluid_properties/src/fluidproperties/SimpleFluidProperties.C

@@ -67,7 +67,8 @@ SimpleFluidProperties::molarMass() const
   return _molar_mass;
 }
 
-Real SimpleFluidProperties::beta_from_p_T(Real /*pressure*/, Real /*temperature*/) const
+Real
+SimpleFluidProperties::beta_from_p_T(Real /*pressure*/, Real /*temperature*/) const
 {
   return _thermal_expansion;
 }
@@ -81,7 +82,8 @@ SimpleFluidProperties::beta_from_p_T(
   dbeta_dT = 0.0;
 }
 
-Real SimpleFluidProperties::cp_from_p_T(Real /*pressure*/, Real /*temperature*/) const
+Real
+SimpleFluidProperties::cp_from_p_T(Real /*pressure*/, Real /*temperature*/) const
 {
   return _cp;
 }
@@ -95,7 +97,11 @@ SimpleFluidProperties::cp_from_p_T(
   dcp_dT = 0.0;
 }
 
-Real SimpleFluidProperties::cp_from_v_e(Real /*v*/, Real /*e*/) const { return _cp; }
+Real
+SimpleFluidProperties::cp_from_v_e(Real /*v*/, Real /*e*/) const
+{
+  return _cp;
+}
 
 void
 SimpleFluidProperties::cp_from_v_e(Real v, Real e, Real & cp, Real & dcp_dv, Real & dcp_de) const
@@ -105,7 +111,8 @@ SimpleFluidProperties::cp_from_v_e(Real v, Real e, Real & cp, Real & dcp_dv, Rea
   dcp_de = 0.0;
 }
 
-Real SimpleFluidProperties::cv_from_p_T(Real /*pressure*/, Real /*temperature*/) const
+Real
+SimpleFluidProperties::cv_from_p_T(Real /*pressure*/, Real /*temperature*/) const
 {
   return _cv;
 }
@@ -119,7 +126,11 @@ SimpleFluidProperties::cv_from_p_T(
   dcv_dT = 0.0;
 }
 
-Real SimpleFluidProperties::cv_from_v_e(Real /*v*/, Real /*e*/) const { return _cv; }
+Real
+SimpleFluidProperties::cv_from_v_e(Real /*v*/, Real /*e*/) const
+{
+  return _cv;
+}
 
 void
 SimpleFluidProperties::cv_from_v_e(Real v, Real e, Real & cv, Real & dcv_dv, Real & dcv_de) const
@@ -171,7 +182,8 @@ SimpleFluidProperties::c_from_v_e(Real v, Real e, Real & c, Real & dc_dv, Real &
   dc_de = 0.0;
 }
 
-Real SimpleFluidProperties::k_from_p_T(Real /*pressure*/, Real /*temperature*/) const
+Real
+SimpleFluidProperties::k_from_p_T(Real /*pressure*/, Real /*temperature*/) const
 {
   return _thermal_conductivity;
 }
@@ -185,7 +197,8 @@ SimpleFluidProperties::k_from_p_T(
   dk_dT = 0;
 }
 
-Real SimpleFluidProperties::k_from_v_e(Real /*v*/, Real /*e*/) const
+Real
+SimpleFluidProperties::k_from_v_e(Real /*v*/, Real /*e*/) const
 {
   return _thermal_conductivity;
 }
@@ -199,7 +212,8 @@ SimpleFluidProperties::k_from_v_e(
   dk_de = 0;
 }
 
-Real SimpleFluidProperties::s_from_p_T(Real /*pressure*/, Real /*temperature*/) const
+Real
+SimpleFluidProperties::s_from_p_T(Real /*pressure*/, Real /*temperature*/) const
 {
   return _specific_entropy;
 }
@@ -213,12 +227,17 @@ SimpleFluidProperties::s_from_p_T(
   ds_dT = 0;
 }
 
-Real SimpleFluidProperties::s_from_h_p(Real /*enthalpy*/, Real /*pressure*/) const
+Real
+SimpleFluidProperties::s_from_h_p(Real /*enthalpy*/, Real /*pressure*/) const
 {
   return _specific_entropy;
 }
 
-Real SimpleFluidProperties::s_from_v_e(Real /*v*/, Real /*e*/) const { return _specific_entropy; }
+Real
+SimpleFluidProperties::s_from_v_e(Real /*v*/, Real /*e*/) const
+{
+  return _specific_entropy;
+}
 
 void
 SimpleFluidProperties::s_from_v_e(
@@ -326,7 +345,7 @@ SimpleFluidProperties::T_from_p_h(Real p, Real h) const
   auto lambda = [&](Real p, Real current_T, Real & new_rho, Real & dh_dp, Real & dh_dT)
   { h_from_p_T(p, current_T, new_rho, dh_dp, dh_dT); };
   return FluidPropertiesUtils::NewtonSolve(
-             p, h, T_initial, _tolerance, lambda, name() + "::T_from_p_h")
+             p, h, T_initial, _tolerance, lambda, name() + "::T_from_p_h", _max_newton_its)
       .first;
 }
 
@@ -429,7 +448,8 @@ SimpleFluidProperties::e_from_v_h(Real v, Real h, Real & e, Real & de_dv, Real &
   de_dh = de_dp * dp_dh + de_dT * dT_dh;
 }
 
-Real SimpleFluidProperties::mu_from_p_T(Real /*pressure*/, Real /*temperature*/) const
+Real
+SimpleFluidProperties::mu_from_p_T(Real /*pressure*/, Real /*temperature*/) const
 {
   return _viscosity;
 }
@@ -443,7 +463,11 @@ SimpleFluidProperties::mu_from_p_T(
   dmu_dT = 0.0;
 }
 
-Real SimpleFluidProperties::mu_from_v_e(Real /*v*/, Real /*e*/) const { return _viscosity; }
+Real
+SimpleFluidProperties::mu_from_v_e(Real /*v*/, Real /*e*/) const
+{
+  return _viscosity;
+}
 
 void
 SimpleFluidProperties::mu_from_v_e(Real v, Real e, Real & mu, Real & dmu_dv, Real & dmu_de) const

modules/fluid_properties/src/fluidproperties/SinglePhaseFluidProperties.C

@@ -28,6 +28,8 @@ SinglePhaseFluidProperties::validParams()
       2e5,
       "p_initial_guess > 0",
       "Pressure initial guess for Newton Method variable set conversion");
+  params.addParam<unsigned int>(
+      "max_newton_its", 100, "Maximum number of Newton iterations for variable set conversions");
   params.addParamNamesToGroup("tolerance T_initial_guess p_initial_guess",
                               "Variable set conversions Newton solve");
 
@@ -39,7 +41,8 @@ SinglePhaseFluidProperties::SinglePhaseFluidProperties(const InputParameters & p
     // downstream apps are creating fluid properties without their parameters, hence the workaround
     _tolerance(isParamValid("tolerance") ? getParam<Real>("tolerance") : 1e-8),
     _T_initial_guess(isParamValid("T_initial_guess") ? getParam<Real>("T_initial_guess") : 400),
-    _p_initial_guess(isParamValid("p_initial_guess") ? getParam<Real>("p_initial_guess") : 2e5)
+    _p_initial_guess(isParamValid("p_initial_guess") ? getParam<Real>("p_initial_guess") : 2e5),
+    _max_newton_its(getParam<unsigned int>("max_newton_its"))
 {
 }
 

modules/fluid_properties/src/fluidproperties/TabulatedBicubicFluidProperties.C

@@ -69,8 +69,8 @@ TabulatedBicubicFluidProperties::constructInterpolation()
   bool conversion_succeeded = true;
   unsigned int fail_counter_ve = 0;
   unsigned int fail_counter_vh = 0;
-  // Create interpolations of (p,T) from (v,e) and (v,h)
-  if (_construct_pT_from_ve || _construct_pT_from_vh)
+  // Create interpolations of (p,T) from (v,e)
+  if (_construct_pT_from_ve)
   {
     // Grids in specific volume and internal energy can be either linear or logarithmic
     // NOTE: this could have been called already when generating tabulated data
@@ -159,31 +159,12 @@ TabulatedBicubicFluidProperties::constructInterpolation()
         std::make_unique<BicubicInterpolation>(_specific_volume, _internal_energy, T_from_v_e);
   }
 
+  // Create interpolations of (p,T) from (v,h)
   if (_construct_pT_from_vh)
   {
-    if (_fp)
-    {
-      // extreme values of enthalpy for the grid bounds
-      Real h1 = h_from_p_T(_pressure_min, _temperature_min);
-      Real h2 = h_from_p_T(_pressure_max, _temperature_min);
-      Real h3 = h_from_p_T(_pressure_min, _temperature_max);
-      Real h4 = h_from_p_T(_pressure_max, _temperature_max);
-      _h_min = std::min({h1, h2, h3, h4});
-      _h_max = std::max({h1, h2, h3, h4});
-    }
-    // if csv exists, get max and min values from csv file
-    else
-    {
-      _h_max = *max_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
-      _h_min = *min_element(_properties[_enthalpy_idx].begin(), _properties[_enthalpy_idx].end());
-    }
-    Real dh = (_h_max - _h_min) / ((Real)_num_e - 1);
-
-    // Create h grid for interpolation
-    // enthalpy & internal energy use same # grid points
-    _enthalpy.resize(_num_e);
-    for (unsigned int j = 0; j < _num_e; ++j)
-      _enthalpy[j] = _h_min + j * dh;
+    // Grids in specific volume and enthalpy can be either linear or logarithmic
+    // NOTE: the specific volume grid could have been created when generating tabulated data
+    createVHGridVectors();
 
     // initialize vectors for interpolation
     std::vector<std::vector<Real>> p_from_v_h(_num_v);