Improved resolution along saturation curve in two-phase tables

Common/DataDrivenConfig.py

@@ -1902,8 +1902,7 @@ def EnablePreferentialDiffusion(self, use_PD:bool=DefaultSettings_FGM.preferenti
         """
 
         self.__preferential_diffusion = use_PD
-        if not self.__custom_Le_av_set:
-            self.SetAverageLewisNumbers()
+        self.SetAverageLewisNumbers()
 
         return
 

Data_Generation/DataGenerator_NICFD.py

@@ -32,6 +32,7 @@
 from tqdm import tqdm
 import csv
 import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
 np.random.seed(2)
 
 #---------------------------------------------------------------------------------------------#
@@ -74,6 +75,7 @@ class DataGenerator_CoolProp(DataGenerator_Base):
 
     __fd_step_size_rho:float = 7e-3
     __fd_step_size_e:float = 7e-3
+    __saturation_curve_interpolator:interp1d = False
 
     def __init__(self, Config_in:Config_NICFD=None):
         DataGenerator_Base.__init__(self, Config_in=Config_in)
@@ -795,6 +797,14 @@ def ComputeSaturationCurve(self, N_samples:int=4000):
         sat_curve=np.column_stack((rho_sat,e_sat))
         return sat_curve
 
+    def GenerateSaturationCurveInterpolator(self):
+        rhoe_sat_curve = self.ComputeSaturationCurve(4000)
+        self.__saturation_curve_interpolator = interp1d(rhoe_sat_curve[:,0], rhoe_sat_curve[:,1],kind="linear",bounds_error=False,fill_value="extrapolate")
+        return
+
+    def GetSaturationCurveStaticEnergy(self, val_rho:float):
+        return self.__saturation_curve_interpolator(val_rho)
+
     def GetFDStepSizes(self):
         return self.__fd_step_size_rho, self.__fd_step_size_e
 

Manifold_Generation/LUT/LUTGenerators.py

@@ -91,10 +91,21 @@ def _checkIfVariableIsValid(self, var_to_check:str):
 
     def _passRefinementOptions(self, mesher:MeshThermodynamicPlane):
         if self._Config.TwoPhase():
-            rhoe_saturation_curve = self.__datagenerator.ComputeSaturationCurve(N_samples=2000)
-            saturation_curve_pts_scaled = self._scaler_controlling_variables.transform(rhoe_saturation_curve)
+            saturation_curve_pts_scaled = self.__calculateSaturationCurvePoints()
             mesher.setSaturationCurvePoints(saturation_curve_pts_scaled)
+
         return super()._passRefinementOptions(mesher)
+
+    def __calculateSaturationCurvePoints(self):
+        self.__datagenerator.GenerateSaturationCurveInterpolator()
+        n_samples = 5000
+        rho_saturation_curve = self.__datagenerator.ComputeSaturationCurve(N_samples=n_samples)[:,0]
+        rho_min, rho_max = np.min(rho_saturation_curve), np.max(rho_saturation_curve)
+        rho_saturation_curve = np.linspace(rho_min, rho_max, n_samples)
+        e_saturation_curve = self.__datagenerator.GetSaturationCurveStaticEnergy(rho_saturation_curve)
+        rhoe_saturation_curve = np.column_stack((rho_saturation_curve, e_saturation_curve))
+        saturation_curve_pts_scaled = self._scaler_controlling_variables.transform(rhoe_saturation_curve)
+        return saturation_curve_pts_scaled
 
     def _initiateMesher(self):
         return MeshThermodynamicPlane()

Manifold_Generation/LUT/MeshTools.py

@@ -3,6 +3,7 @@
 from concave_hull import concave_hull, concave_hull_indexes
 from collections.abc import Callable
 from Common.CommonMethods import shoelace, FiniteDifferenceDerivative
+from shapely.geometry import Point, Polygon
 
 def default_refinement_function(this, x:float, y:float, z:float):
     return 1.0
@@ -424,19 +425,13 @@ def __createSaturationCurveOffsets(self, normal_vector:np.ndarray[float]):
     def __clipSaturationCurveToPlane(self, pts_offset_lower:np.ndarray[float], pts_offset_upper:np.ndarray[float]):
 
         hull_pts_orig = concave_hull(self._pointCloud_hullNodes, length_threshold=self._base_cell_size)
-        ref_area = shoelace(hull_pts_orig)
+        polygon_thermodynamic_plane = Polygon(hull_pts_orig)
         within_hull = np.zeros(len(self.__saturation_curve_points),dtype=bool)
 
         for i in range(len(self.__saturation_curve_points)):
-            XY_with_pt = np.vstack((self._pointCloud_hullNodes[:,:2], pts_offset_upper[i,:]))
-            hull_n = concave_hull(XY_with_pt, length_threshold=self._base_cell_size)
-            area_n = shoelace(hull_n)
-            within_hull_upper = (area_n <= ref_area )
-            XY_with_pt = np.vstack((self._pointCloud_hullNodes[:,:2], pts_offset_lower[i,:]))
-            hull_n = concave_hull(XY_with_pt, length_threshold=self._base_cell_size)
-            area_n = shoelace(hull_n)
-            within_hull_lower = (area_n <= ref_area )
-            within_hull[i] = (within_hull_upper and within_hull_lower)
+            pt_upper = Point(pts_offset_upper[i])
+            pt_lower = Point(pts_offset_lower[i])
+            within_hull[i] = polygon_thermodynamic_plane.contains(pt_upper) and polygon_thermodynamic_plane.contains(pt_lower)
 
         return within_hull
 

Manifold_Generation/MLP/Trainers_NICFD/Trainers.py

@@ -41,12 +41,11 @@
 import matplotlib.pyplot as plt
 from matplotlib import ticker
 
-import CoolProp as CoolP
-import CoolProp.CoolProp as CP
 
 from Common.DataDrivenConfig import Config_NICFD
 from Common.CommonMethods import readStateDataFromFile
 from Common.Properties import DefaultSettings_NICFD, EntropicVars
+from Data_Generation.DataGenerator_NICFD import DataGenerator_CoolProp
 from Manifold_Generation.MLP.Trainer_Base import TensorFlowFit,PhysicsInformedTrainer,TrainMLP
 
 LabelPairing = {EntropicVars.s.name:r"Entropy $(s)[J/kg]$",\
@@ -58,45 +57,6 @@
                 EntropicVars.dpdrho_e.name:r"Pressure-density derivative $\left(\left.\frac{\partial p}{\partial \rho}\right|_e\right)$",\
                 EntropicVars.dpde_rho.name:r"Pressure-energy derivative $\left(\left.\frac{\partial p}{\partial e}\right|_\rho\right)$"}
 
-def GetStateVector(fluid:CP.AbstractState):
-    """Extract the fluid thermodynamic properties from the CoolProp abstract state object.
-
-    :param fluid: CoolProp abstract state of the equation of state.
-    :type fluid: CoolProp.CoolProp.AbstractState
-    :return: thermodynamic state properties, phase
-    :rtype: np.ndarray, bool
-    """
-    state_vector_vals = np.ones(EntropicVars.N_STATE_VARS.value)
-    correct_phase = True
-    accepted_phases:list[int] = [CoolP.iphase_gas, CoolP.iphase_supercritical_gas, CoolP.iphase_supercritical]
-    if fluid.phase() in accepted_phases:
-        state_vector_vals[EntropicVars.s.value] = fluid.smass()
-        state_vector_vals[EntropicVars.dsde_rho.value] = fluid.first_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass)
-        state_vector_vals[EntropicVars.dsdrho_e.value] = fluid.first_partial_deriv(CP.iSmass, CP.iDmass, CP.iUmass)
-        state_vector_vals[EntropicVars.d2sde2.value] = fluid.second_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass, CP.iUmass, CP.iDmass)
-        state_vector_vals[EntropicVars.d2sdedrho.value] = fluid.second_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass, CP.iDmass, CP.iUmass)
-        state_vector_vals[EntropicVars.d2sdrho2.value] = fluid.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iUmass, CP.iDmass, CP.iUmass)
-        state_vector_vals[EntropicVars.Density.value] = fluid.rhomass()
-        state_vector_vals[EntropicVars.Energy.value] = fluid.umass()
-        state_vector_vals[EntropicVars.T.value] = fluid.T()
-        state_vector_vals[EntropicVars.p.value] = fluid.p()
-        state_vector_vals[EntropicVars.c2.value] = fluid.speed_sound()**2
-        state_vector_vals[EntropicVars.dTde_rho.value] = fluid.first_partial_deriv(CP.iT, CP.iUmass, CP.iDmass)
-        state_vector_vals[EntropicVars.dTdrho_e.value] = fluid.first_partial_deriv(CP.iT, CP.iDmass, CP.iUmass)
-        state_vector_vals[EntropicVars.dpde_rho.value] = fluid.first_partial_deriv(CP.iP, CP.iUmass, CP.iDmass)
-        state_vector_vals[EntropicVars.dpdrho_e.value] = fluid.first_partial_deriv(CP.iP, CP.iDmass, CP.iUmass)
-        state_vector_vals[EntropicVars.dhde_rho.value] = fluid.first_partial_deriv(CP.iHmass, CP.iUmass, CP.iDmass)
-        state_vector_vals[EntropicVars.dhdrho_e.value] = fluid.first_partial_deriv(CP.iHmass, CP.iDmass, CP.iUmass)
-        state_vector_vals[EntropicVars.dhdp_rho.value] = fluid.first_partial_deriv(CP.iHmass, CP.iP, CP.iDmass)
-        state_vector_vals[EntropicVars.dhdrho_p.value] = fluid.first_partial_deriv(CP.iHmass, CP.iDmass, CP.iP)
-        state_vector_vals[EntropicVars.dsdp_rho.value] = fluid.first_partial_deriv(CP.iSmass, CP.iP, CP.iDmass)
-        state_vector_vals[EntropicVars.dsdrho_p.value] = fluid.first_partial_deriv(CP.iSmass, CP.iDmass, CP.iP)
-        state_vector_vals[EntropicVars.cp.value] = fluid.cpmass()
-    else:
-        correct_phase = False
-        state_vector_vals[:] = None
-    return state_vector_vals, correct_phase
-
 def ComputeRhoEGridData(config:Config_NICFD):
     """Compute the fluid thermodynamic properties for a coarse rho-e grid for visualization purposes.
 
@@ -117,16 +77,18 @@ def ComputeRhoEGridData(config:Config_NICFD):
 
     state_cp = np.empty([Np_x, Np_y, EntropicVars.N_STATE_VARS.value])
     state_cp[:] = np.nan
-    fluid = CP.AbstractState(config.GetEquationOfState(), config.GetFluid())
+    dg = DataGenerator_CoolProp(config)
+
     for i in range(Np_x):
         for j in range(Np_y):
             rho = xx[i,j]
             e = yy[i,j]
             state_cp[i,j,EntropicVars.Density.value] = rho
             state_cp[i,j,EntropicVars.Energy.value] = e
+
             try:
-                fluid.update(CP.DmassUmass_INPUTS, rho, e)
-                state_vector_vals, correct_phase = GetStateVector(fluid)
+                dg.UpdateFluid(rho, e)
+                state_vector_vals, correct_phase = dg.GetStateVector()
                 if correct_phase:
                     state_cp[i,j,:] = state_vector_vals
             except: