diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml
index 882264a..774ac58 100644
--- a/.github/workflows/regression.yml
+++ b/.github/workflows/regression.yml
@@ -18,18 +18,18 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: Pre Cleanup
- uses: docker://ghcr.io/su2code/su2_dataminer_test
+ uses: docker://ghcr.io/su2code/su2_dataminer_test:v2
with:
entrypoint: /bin/rm
args: -rf src/
- name: Run regression tests
- uses: docker://ghcr.io/su2code/su2_dataminer_test
+ uses: docker://ghcr.io/su2code/su2_dataminer_test:v2
with:
args: -b ${{ github.ref }} -s run_regression.py
- name: Post Cleanup
- uses: docker://ghcr.io/su2code/su2_dataminer_test
+ uses: docker://ghcr.io/su2code/su2_dataminer_test:v2
with:
entrypoint: /bin/rm
args: -rf src/
@@ -39,18 +39,18 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: Pre Cleanup
- uses: docker://ghcr.io/su2code/su2_dataminer_test
+ uses: docker://ghcr.io/su2code/su2_dataminer_test:v2
with:
entrypoint: /bin/rm
args: -rf src/
- name: MLPCpp integration tests
- uses: docker://ghcr.io/su2code/su2_dataminer_test
+ uses: docker://ghcr.io/su2code/su2_dataminer_test:v2
with:
args: -b ${{ github.ref }} -s MLPCpp_tests.py -m main
- name: Post Cleanup
- uses: docker://ghcr.io/su2code/su2_dataminer_test
+ uses: docker://ghcr.io/su2code/su2_dataminer_test:v2
with:
entrypoint: /bin/rm
args: -rf src/
\ No newline at end of file
diff --git a/Common/DataDrivenConfig.py b/Common/DataDrivenConfig.py
index fb1a1c7..52f3c51 100644
--- a/Common/DataDrivenConfig.py
+++ b/Common/DataDrivenConfig.py
@@ -57,6 +57,17 @@ class Config_NICFD(Config):
# Fluid definition settings
__fluid_names:list[str] = ["MM"] # List of fluid names used for data generation.
__fluid_string:str="MM" # Fluid string for defining the abstract state in CoolProp
+
+ __calc_transport_properties:bool = False
+ __viscosity_model:str = DefaultSettings_NICFD.viscosity_model
+ __conductivity_model:str = DefaultSettings_NICFD.conductivity_model
+
+ # Phases to include in fluid data.
+ __gasphase:bool = True
+ __twophase:bool = False
+ __liquidphase:bool = False
+ __supercritical:bool = True
+
__EOS_type:str=DefaultSettings_NICFD.EOS_type # Equation of state used by CoolProp
__fluid_mole_fractions:list[float] = [1.0] # Mole fractions for components in fluid mixture.
__use_PT:bool = DefaultSettings_NICFD.use_PT_grid # Use a pressure-temperature based grid for fluid training data.
@@ -78,7 +89,7 @@ class Config_NICFD(Config):
_state_vars:list[str] = ["s", "T","p","c2"] # State variable names for which the physics-informed MLP is trained.
# Table Generation Settings
-
+ __Table_discretization:str = DefaultSettings_NICFD.tabulation_method
__Table_base_cell_size:float = None # Table base cell size per table level.
__Table_ref_cell_size:float = None # Refined cell size per table level.
__Table_ref_radius:float = None # Refinement radius within which refined cell size is applied.
@@ -151,6 +162,15 @@ def PrintBanner(self):
print("Density range: %.2f kg/m3 -> %.2f kg/m3 (%i steps)" % (self.__Rho_lower, self.__Rho_upper, self.__Np_P))
print("")
print("State variables considered during physics-informed learning: "+", ".join((v for v in self._state_vars)))
+
+ if self.__twophase:
+ print("Including two-phase fluid data.")
+ if self.__calc_transport_properties:
+ print("Including transport properties in fluid data.")
+ if self.__twophase:
+ print("Two-phase conductivity model: %s" % self.__conductivity_model)
+ print("Two-phase viscosity model: %s" % self.__viscosity_model)
+
return
@@ -208,6 +228,137 @@ def SetEquationOfState(self, EOS_type_in:str=DefaultSettings_NICFD.EOS_type):
self.__EOS_type=EOS_type_in.upper()
return
+ def IncludeTransportProperties(self, calc_transport_properties:bool=False):
+ """Include transport properties in fluid data calculation
+
+ :param calc_transport_properties: evaluate transport properties, defaults to False
+ :type calc_transport_properties: bool, optional
+ """
+ self.__calc_transport_properties = calc_transport_properties
+ return
+
+ def CalcTransportProperties(self):
+ """Whether transport properties are included in the fluid data set.
+
+ :return: calculation of transport properties.
+ :rtype: bool
+ """
+ return self.__calc_transport_properties
+
+ def SetConductivityModel(self, conductivity_model:str=DefaultSettings_NICFD.conductivity_model):
+ """Specify the two-phase conductivity model.
+
+ :param conductivity_model: two-phase conductivity model option, defaults to "volume"
+ :type conductivity_model: str, optional
+ :raises Exception: if the specified option is not supported.
+ """
+ if conductivity_model not in DefaultSettings_NICFD.conductivity_models:
+ raise Exception("Two-phase conductivity model should be one of the following: " + ",".join(c for c in DefaultSettings_NICFD.conductivity_models))
+ self.__conductivity_model = conductivity_model
+ if not self.TwoPhase():
+ self.EnableTwophase(True)
+ print("Two-phase conductivity model specified, including two-phase fluid data.")
+ return
+
+ def GetConductivityModel(self):
+ """Get two-phase conductivity model.
+
+ :return: two-phase conductivity model option
+ :rtype: str
+ """
+ return self.__conductivity_model
+
+ def SetViscosityModel(self, viscosity_model:str=DefaultSettings_NICFD.viscosity_model):
+ """Specify the two-phase viscosity model.
+
+ :param viscosity_model: two-phase viscosity model option, defaults to "mcadams"
+ :type viscosity_model: str, optional
+ :raises Exception: if the specified option is not supported.
+ """
+ if viscosity_model not in DefaultSettings_NICFD.viscosity_models:
+ raise Exception("Two-phase viscosity model should be one of the following: " + ",".join(c for c in DefaultSettings_NICFD.viscosity_models))
+ self.__viscosity_model = viscosity_model
+ if not self.TwoPhase():
+ self.EnableTwophase(True)
+ print("Two-phase viscosity model specified, including two-phase fluid data.")
+ return
+
+ def GetViscosityModel(self):
+ """Get two-phase viscosity model.
+
+ :return: two-phase viscosity model option
+ :rtype: str
+ """
+ return self.__viscosity_model
+
+ def EnableGasPhase(self, gas_phase:bool=True):
+ """Include thermophysical state data from the fluid in gas phase
+
+ :param gas_phase: include gas phase data, defaults to True
+ :type gas_phase: bool, optional
+ """
+ self.__gasphase = gas_phase
+ return
+
+ def GasPhase(self):
+ """Whether thermophysical state data from the fluid in the gaseous phase are included.
+
+ :return: inclusion of gas phase data.
+ :rtype: bool
+ """
+ return self.__gasphase
+
+ def EnableSuperCritical(self, supercritical:bool=True):
+ """Include thermophysical state data of the fluid in supercritial, and of the supercritical gas and liquid phase if specified.
+
+ :param supercritical: include supercritical phase data, defaults to True
+ :type supercritical: bool, optional
+ """
+ self.__supercritical = supercritical
+ return
+
+ def EnableTwophase(self, two_phase:bool=False):
+ """Include two-phase thermophysical data of the fluid.
+
+ :param two_phase: include two-phase data in fluid data, defaults to False
+ :type two_phase: bool, optional
+ """
+ self.__twophase = two_phase
+ return
+
+ def EnableLiquidPhase(self, liquid_phase:bool=False):
+ """Include thermodynamic state data from fluid in the liquid phase.
+
+ :param liquid_phase: include liquid-phase data, defaults to False
+ :type liquid_phase: bool, optional
+ """
+ self.__liquidphase = liquid_phase
+ return
+
+ def TwoPhase(self):
+ """Whether thermophysical state data from the fluid in the two-phase are included.
+
+ :return: inclusion of two-phase data.
+ :rtype: bool
+ """
+ return self.__twophase
+
+ def LiquidPhase(self):
+ """Whether thermophysical state data from the fluid in the liquid phase are included.
+
+ :return: inclusion of liquid phase data.
+ :rtype: bool
+ """
+ return self.__liquidphase
+
+ def SuperCritical(self):
+ """Whether thermophysical state data from the fluid in the supercritical phases are included.
+
+ :return: inclusion of supercritical phase data.
+ :rtype: bool
+ """
+ return self.__supercritical
+
def GetEquationOfState(self):
"""Retrieve the equation of state backend used by CoolProp for fluid data calculations.
@@ -262,9 +413,9 @@ def UseAutoRange(self, use_auto_range:bool=True):
return
def GetAutoRange(self):
- """Whether fluid data ranges are automatically determined by CoolProp
+ """The span of the thermodynamic state space is determined automatically.
- :return: boolean if range is automatically determined
+ :return: whether automatic ranging is used.
:rtype: bool
"""
return self.__use_auto_range
@@ -358,9 +509,9 @@ def GetNpEnergy(self):
def SetDensityBounds(self, Rho_lower:float=DefaultSettings_NICFD.Rho_min, Rho_upper:float=DefaultSettings_NICFD.Rho_max):
"""Define the density bounds of the density-energy based fluid data grid.
- :param Rho_lower: lower limit density value in kg/m3, defaults to DefaultSettings_NICFD.Rho_min
+ :param Rho_lower: lower limit density value, defaults to DefaultSettings_NICFD.Rho_min
:type Rho_lower: float, optional
- :param Rho_upper: upper limit for density in kg/m3, defaults to DefaultSettings_NICFD.Rho_max
+ :param Rho_upper: upper limit for density, defaults to DefaultSettings_NICFD.Rho_max
:type Rho_upper: float, optional
:raises Exception: if lower value for density exceeds upper value.
"""
@@ -489,6 +640,26 @@ def GetNpPressure(self):
"""
return self.__Np_P
+ def SetTableDiscretization(self, table_method:str=DefaultSettings_NICFD.tabulation_method):
+ """Specify how thermodynamic state space is discretized for look-up table.
+
+ :param table_method: discretization method, defaults to "cartesian"
+ :type table_method: str, optional
+ :raises Exception: if method is not supported
+ """
+ if table_method not in DefaultSettings_NICFD.tabulation_options:
+ raise Exception("Table discretization method should be one of the following: %s" % ",".join(t for t in DefaultSettings_NICFD.tabulation_options))
+ self.__Table_discretization = table_method
+ return
+
+ def GetTableDiscretization(self):
+ """Get thermodynamic state space discretization method.
+
+ :return: table discretization method.
+ :rtype: str
+ """
+ return self.__Table_discretization
+
def SetTableCellSize(self, base_cell_size:float, refined_cell_size:float=None):
"""Define the base and optional refined 2D table cell sizes.
@@ -1019,10 +1190,11 @@ def SetReactionMechanism(self, mechanism_input:str=DefaultSettings_FGM.reaction_
return
def GetReactionMechanism(self):
- """Get the reaction mechanism used for flamelet generation.
-
+ """
+ Get the reaction mechanism used for flamelet generation.
:return: reaction mechanism name.
:rtype: str
+
"""
return self.__reaction_mechanism
@@ -1079,10 +1251,11 @@ def SetMixtureBounds(self, mix_lower:float=DefaultSettings_FGM.eq_ratio_min, mix
return
def GetMixtureBounds(self):
- """Get the mixture status bounds.
-
+ """
+ Get the mixture status bounds.
:return: List containing lower and upper mixture status values.
:rtype: list[float]
+
"""
return [self.__mix_status_lower, self.__mix_status_upper]
@@ -1102,7 +1275,8 @@ def SetNpMix(self, input:int=DefaultSettings_FGM.Np_eq):
return
def GetNpMix(self):
- """Get the number of divisions between the lean and rich mixture status for flamelet generation.
+ """
+ Get the number of divisions between the lean and rich mixture status for flamelet generation.
:return: number of divisions between rich and lean.
:rtype: int
@@ -1381,8 +1555,8 @@ def SetPassiveSpecies(self, passive_species:list[str]=[]):
"""
Set the passive transported species for which source terms should be saved in the manifold.
- :param passive_species: list of species names.
- :type passive_species: list[str]
+ :param __passive_species: list of species names.
+ :type __passive_species: list[str]
"""
self.__passive_species = passive_species
return
diff --git a/Common/Properties.py b/Common/Properties.py
index fc1c85f..a39f235 100644
--- a/Common/Properties.py
+++ b/Common/Properties.py
@@ -65,6 +65,11 @@ class EntropicVars(Enum):
dsdp_rho=auto()
dsdrho_p=auto()
cp=auto()
+ cv=auto()
+ Enthalpy=auto()
+ Conductivity=auto()
+ ViscosityDyn=auto()
+ VaporQuality=auto()
N_STATE_VARS=auto()
class FGMVars(Enum):
@@ -137,6 +142,11 @@ class DefaultSettings_NICFD(DefaultProperties):
fluid_name:str = "Air"
EOS_type:str = "HEOS"
+ conductivity_models:list[str] = ["volume", "mass"]
+ viscosity_models:list[str] = ["mcadams", "cicchitti", "dukler"]
+ conductivity_model:str = "volume"
+ viscosity_model:str = "mcadams"
+
use_PT_grid:bool = False
controlling_variables:list[str] = [EntropicVars.Density.name, \
@@ -152,6 +162,9 @@ class DefaultSettings_NICFD(DefaultProperties):
config_type:str = "EntropicAI"
supported_state_vars:list[str] = ["s","T","p","c2","dTdrho_e","dTde_rho","dpdrho_e","dpde_rho"]
supported_backends:list[str] = ["HEOS","PR", "SRK", "IF97","REFPROP"]
+
+ tabulation_options:list[str] = ["cartesian","adaptive"]
+ tabulation_method:str="cartesian"
class DefaultSettings_FGM(DefaultProperties):
config_name:str = "config_FGM"
@@ -213,4 +226,5 @@ class DefaultSettings_FGM(DefaultProperties):
"exponential" : tf.keras.activations.exponential,\
"gelu" : tf.keras.activations.gelu,\
"sigmoid" : tf.keras.activations.sigmoid,\
- "swish" : tf.keras.activations.swish}
\ No newline at end of file
+ "swish" : tf.keras.activations.swish}
+
diff --git a/Data_Generation/DataGenerator_NICFD.py b/Data_Generation/DataGenerator_NICFD.py
index be4ec4e..d00e63e 100644
--- a/Data_Generation/DataGenerator_NICFD.py
+++ b/Data_Generation/DataGenerator_NICFD.py
@@ -47,7 +47,7 @@ class DataGenerator_CoolProp(DataGenerator_Base):
"""
_Config:Config_NICFD
fluid = None
- __accepted_phases:list[int] = [CoolP.iphase_gas, CoolP.iphase_supercritical_gas, CoolP.iphase_supercritical]
+ __accepted_phases:list[int] = []
# Pressure and temperature limits
__use_PT:bool = DefaultSettings_NICFD.use_PT_grid
__T_min:float = DefaultSettings_NICFD.T_min
@@ -69,17 +69,13 @@ class DataGenerator_CoolProp(DataGenerator_Base):
# Entropy derivatives.
__StateVars_fluid:np.ndarray[float] = None
- __StateVars_additional:np.ndarray[float] = None
-
__success_locations:np.ndarray[bool] = None
__mixture:bool = False
- def __init__(self, Config_in:Config_NICFD=None):
- """Constructor, load settings from SU2 DataMiner configuration.
+ __fd_step_size_rho:float = 7e-3
+ __fd_step_size_e:float = 7e-3
- :param Config_in: SU2 DataMiner configuration class, defaults to None
- :type Config_in: Config_NICFD, optional
- """
+ def __init__(self, Config_in:Config_NICFD=None):
DataGenerator_Base.__init__(self, Config_in=Config_in)
if Config_in is None:
@@ -91,6 +87,12 @@ def __init__(self, Config_in:Config_NICFD=None):
if len(self._Config.GetFluidNames()) > 1:
self.__mixture = True
+ # Define phases for which to generate fluid data.
+ self.EnableTwophase(self._Config.TwoPhase())
+ self.EnableLiquidPhase(self._Config.LiquidPhase())
+ self.EnableGasPhase(self._Config.GasPhase())
+ self.EnableSuperCritical(self._Config.SuperCritical())
+
self.fluid = CP.AbstractState(self._Config.GetEquationOfState(), self._Config.GetFluid())
self.__auto_range = self._Config.GetAutoRange()
@@ -168,16 +170,12 @@ def PreprocessData(self):
X_max = self.__rho_max
Y_min = self.__e_min
Y_max = self.__e_max
-
X_range = (X_min - X_max) * np.cos(np.linspace(0, 0.5*np.pi, self.__Np_X)) + X_max
Y_range = np.linspace(Y_min, Y_max, self.__Np_Y)
-
self.__X_grid, self.__Y_grid = np.meshgrid(X_range, Y_range)
return
def UpdateConfig(self):
- """Update the bounds of the thermodynamic state data stored in the configuration class.
- """
if self.__use_PT:
self._Config.SetPressureBounds(self.__P_min, self.__P_max)
self._Config.SetTemperatureBounds(self.__T_min, self.__T_max)
@@ -323,6 +321,166 @@ def GetPressureBounds(self):
"""
return [self.__P_min, self.__P_max]
+ def GetDensityBounds(self):
+ """Get fluid density bounds values.
+
+ :return: minimum and maximum fluid density in data set.
+ :rtype: tuple
+ """
+ return self.__rho_min, self.__rho_max
+
+ def GetEnergyBounds(self):
+ """Get fluid static energy bounds values.
+
+ :return: minimum and maximum fluid static energy in data set.
+ :rtype: tuple
+ """
+ return self.__e_min, self.__e_max
+
+ def IncludeTransportProperties(self, calc_transport_properties:bool=False):
+ """Include transport properties in fluid data calculation
+
+ :param calc_transport_properties: evaluate transport properties, defaults to False
+ :type calc_transport_properties: bool, optional
+ """
+ self._Config.IncludeTransportProperties(calc_transport_properties)
+ return
+
+ def CalcTransportProperties(self):
+ return self._Config.CalcTransportProperties()
+
+ def EnableGasPhase(self, gas_phase:bool=True):
+ """Include thermophysical data of the fluid in gas phase.
+
+ :param gas_phase: include gas phase data, defaults to True
+ :type gas_phase: bool, optional
+ """
+ self._Config.EnableGasPhase(gas_phase)
+ if self._Config.GasPhase() and CoolP.iphase_gas not in self.__accepted_phases:
+ self.__accepted_phases.append(CoolP.iphase_gas)
+ if self._Config.SuperCritical():
+ self.__accepted_phases.append(CoolP.iphase_supercritical_gas)
+ if not self._Config.GasPhase() and CoolP.iphase_gas in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_gas)
+ if CoolP.iphase_supercritical_gas in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_supercritical_gas)
+ return
+
+ def GasPhase(self):
+ """Whether gas phase data are included in the fluid data set.
+
+ :return: gas phase data included in fluid data set.
+ :rtype: bool
+ """
+ return self._Config.GasPhase()
+
+ def EnableTwophase(self, two_phase:bool=False):
+ """Include two-phase region in fluid data.
+
+ :param two_phase: include two-phase data in fluid data, defaults to False
+ :type two_phase: bool, optional
+ """
+ self._Config.EnableTwophase(two_phase)
+ if self._Config.TwoPhase() and CoolP.iphase_twophase not in self.__accepted_phases:
+ self.__accepted_phases.append(CoolP.iphase_twophase)
+ if not self._Config.TwoPhase() and CoolP.iphase_twophase in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_twophase)
+ return
+
+ def TwoPhase(self):
+ return self._Config.TwoPhase()
+
+ def EnableLiquidPhase(self, liquid_phase:bool=False):
+ """Include thermophysical data of the fluid in liquid phase.
+
+ :param liquid_phase: include liquid phase fluid data, defaults to False
+ :type liquid_phase: bool, optional
+ """
+ self._Config.EnableLiquidPhase(liquid_phase)
+ if self._Config.LiquidPhase() and CoolP.iphase_liquid not in self.__accepted_phases:
+ self.__accepted_phases.append(CoolP.iphase_liquid)
+ if self._Config.SuperCritical():
+ self.__accepted_phases.append(CoolP.iphase_supercritical_liquid)
+ if not self._Config.LiquidPhase() and CoolP.iphase_liquid in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_liquid)
+ if CoolP.iphase_supercritical_liquid in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_supercritical_liquid)
+ return
+
+ def LiquidPhase(self):
+ """Whether to include fluid data in liquid phase.
+
+ :return: liquid phase data are included in the data set.
+ :rtype: bool
+ """
+ return self._Config.LiquidPhase()
+
+ def EnableSuperCritical(self, supercritical:bool=True):
+ """Whether to include fluid data in the supercritical phase in the data set.
+
+ :param supercritical: include data in the supercritical phase, defaults to True
+ :type supercritical: bool, optional
+ """
+ self._Config.EnableSuperCritical(supercritical)
+ if self._Config.SuperCritical():
+ if CoolP.iphase_supercritical not in self.__accepted_phases:
+ self.__accepted_phases.append(CoolP.iphase_supercritical)
+ if CoolP.iphase_gas in self.__accepted_phases and CoolP.iphase_supercritical_gas not in self.__accepted_phases:
+ self.__accepted_phases.append(CoolP.iphase_supercritical_gas)
+ if (CoolP.iphase_liquid in self.__accepted_phases) and (CoolP.iphase_supercritical_liquid not in self.__accepted_phases):
+ self.__accepted_phases.append(CoolP.iphase_supercritical_liquid)
+
+ else:
+ if CoolP.iphase_supercritical in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_supercritical)
+ if CoolP.iphase_supercritical_gas in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_supercritical_gas)
+ if CoolP.iphase_supercritical_liquid in self.__accepted_phases:
+ self.__accepted_phases.remove(CoolP.iphase_supercritical_liquid)
+ return
+
+ def SetConductivityModel(self, conductivity_model:str=DefaultSettings_NICFD.conductivity_model):
+ """Specify the two-phase conductivity model.
+
+ :param conductivity_model: two-phase conductivity model option, defaults to "volume"
+ :type conductivity_model: str, optional
+ :raises Exception: if the specified option is not supported.
+ """
+ self._Config.SetConductivityModel(conductivity_model)
+ if not self.TwoPhase():
+ self.EnableTwophase(True)
+ print("Two-phase conductivity model specified, including two-phase fluid data.")
+ return
+
+ def GetConductivityModel(self):
+ """Get two-phase conductivity model.
+
+ :return: two-phase conductivity model option
+ :rtype: str
+ """
+ return self._Config.GetConductivityModel()
+
+ def SetViscosityModel(self, viscosity_model:str=DefaultSettings_NICFD.viscosity_model):
+ """Specify the two-phase viscosity model.
+
+ :param viscosity_model: two-phase viscosity model option, defaults to "mcadams"
+ :type viscosity_model: str, optional
+ :raises Exception: if the specified option is not supported.
+ """
+ self._Config.SetViscosityModel(viscosity_model)
+ if not self.TwoPhase():
+ self.EnableTwophase(True)
+ print("Two-phase viscosity model specified, including two-phase fluid data.")
+ return
+
+ def GetViscosityModel(self):
+ """Get two-phase viscosity model.
+
+ :return: two-phase viscosity model option
+ :rtype: str
+ """
+ return self._Config.GetViscosityModel()
+
def ComputeData(self):
super().ComputeData()
@@ -333,7 +491,6 @@ def ComputeData(self):
self.__StateVars_fluid = np.zeros([self.__Np_X, self.__Np_Y, EntropicVars.N_STATE_VARS.value])
self.__success_locations = np.ones([self.__Np_X, self.__Np_Y],dtype=bool)
-
# Loop over density-based or pressure-based grid.
for i in tqdm(range(self.__Np_X)):
for j in range(self.__Np_Y):
@@ -346,45 +503,599 @@ def ComputeData(self):
T = self.fluid.T()
self.fluid.update(CP.PT_INPUTS, p, T)
# Check if fluid phase is not vapor or liquid
- self.__StateVars_fluid[i,j,:], self.__success_locations[i,j] = self.__GetStateVector()
+ self.__StateVars_fluid[i,j,:], self.__success_locations[i,j] = self.GetStateVector()
except:
self.__success_locations[i,j] = False
self.__StateVars_fluid[i, j, :] = None
-
return
- def __GetStateVector(self):
+ def __TransportProperties(self, q:float):
+ """Evaluate fluid transport properties
+
+ :param q: vapor quality
+ :type q: float
+ :return: dynamic viscosity, thermal conductivity, vapor quality
+ :rtype: float, float, float
+ """
+ p = self.fluid.p()
+ rho = self.fluid.rhomass()
+ e = self.fluid.umass()
+ phase = self.fluid.phase()
+ is_twophase = (phase == CoolP.iphase_twophase)
+ if is_twophase:
+ vapor_q = q
+ sat_props = self.__get_saturated_transport_properties(p)
+
+ if sat_props is not None:
+ # Compute void fraction for volume-weighted models
+ alpha = self.__compute_void_fraction(
+ q, sat_props['rho_l'], sat_props['rho_g'])
+ # Mixture viscosity
+ viscosity = self.__compute_twophase_viscosity(
+ q, sat_props['mu_l'], sat_props['mu_g'],
+ sat_props['rho_l'], sat_props['rho_g'], alpha)
+
+ # Mixture thermal conductivity
+ conductivity = self.__compute_twophase_conductivity(
+ q, sat_props['k_l'], sat_props['k_g'],
+ sat_props['rho_l'], sat_props['rho_g'], alpha)
+ else:
+ # Fallback: try to get properties at current state (might fail)
+ try:
+ self.fluid.update(CP.DmassUmass_INPUTS, rho, e)
+ viscosity = self.fluid.viscosity()
+ conductivity = self.fluid.conductivity()
+ except:
+ # Last resort: interpolate based on quality
+ viscosity = 1e-5 # Reasonable default
+ conductivity = 0.1 # Reasonable default
+ else:
+ vapor_q = -1.0
+ viscosity = self.fluid.viscosity()
+ conductivity = self.fluid.conductivity()
+ return viscosity, conductivity, vapor_q
+
+ def __EntropicEoS(self, rho, e, s, derivs, state_vector_struct:dict):
+ state_vector_struct[EntropicVars.Density.name] = rho
+ state_vector_struct[EntropicVars.Energy.name] = e
+ state_vector_struct[EntropicVars.s.name] = s
+ dsdrho_e = derivs["dsdrho_e"]
+ state_vector_struct[EntropicVars.dsdrho_e.name] = dsdrho_e
+ dsde_rho = derivs["dsde_rho"]
+ state_vector_struct[EntropicVars.dsde_rho.name] = dsde_rho
+ d2sdrho2 = derivs["d2sdrho2"]
+ state_vector_struct[EntropicVars.d2sdrho2.name] = d2sdrho2
+ d2sde2 = derivs["d2sde2"]
+ state_vector_struct[EntropicVars.d2sde2.name] = d2sde2
+ d2sdedrho = derivs["d2sdedrho"]
+ state_vector_struct[EntropicVars.d2sdedrho.name] = d2sdedrho
+ Temperature = pow(dsde_rho, -1)
+ state_vector_struct[EntropicVars.T.name] = Temperature
+ Pressure = -rho * rho * Temperature * dsdrho_e
+ state_vector_struct[EntropicVars.p.name] = Pressure
+ dPde_rho = -rho*rho * Temperature * (-Temperature * (d2sde2 * dsdrho_e) + d2sdedrho)
+ state_vector_struct[EntropicVars.dpde_rho.name] = dPde_rho
+ dPdrho_e = - rho * Temperature * (dsdrho_e * (2 - rho * Temperature * d2sdedrho) + rho * d2sdrho2)
+ state_vector_struct[EntropicVars.dpdrho_e.name] = dPdrho_e
+ SoundSpeed2 = dPdrho_e - (dsdrho_e / dsde_rho) * dPde_rho
+ state_vector_struct[EntropicVars.c2.name] = SoundSpeed2
+ dTde_rho = -Temperature * Temperature * d2sde2
+ state_vector_struct[EntropicVars.dTde_rho.name] = dTde_rho
+ dTdrho_e = -Temperature * Temperature * d2sdedrho
+ state_vector_struct[EntropicVars.dTdrho_e.name] = dTdrho_e
+
+ dhde_rho = 1 + dPde_rho / rho
+ dhdrho_e = -Pressure * np.power(rho, -2) + dPdrho_e / rho
+ state_vector_struct[EntropicVars.dhde_rho.name] = dhdrho_e
+ state_vector_struct[EntropicVars.dhdrho_e.name] = dhde_rho
+
+ # drhode_p = -dPde_rho / dPdrho_e
+ # dTde_p = dTde_rho + dTdrho_e * drhode_p
+ # dhde_p = dhde_rho + drhode_p*dhdrho_e
+ # Cp = dhde_p / (dTde_p+np.sign(dTde_p)*1e-8)
+ # Cv = 1 / (dTde_rho+np.sign(dTde_rho)*1e-8)
+
+ X = self.fluid.Q()
+ self.fluid.update(CoolP.PQ_INPUTS, Pressure, 1)
+ cp_vap, cv_vap,rho_vap = self.fluid.cpmass(), self.fluid.cvmass(), self.fluid.rhomass()
+ self.fluid.update(CoolP.PQ_INPUTS, Pressure, 0)
+ cp_liq, cv_liq = self.fluid.cpmass(), self.fluid.cvmass()
+ alpha=X*rho/rho_vap
+ Cp =alpha*cp_vap+(1-alpha)*cp_liq
+ Cv =alpha*cv_vap+(1-alpha)*cv_liq
+ self.fluid.update(CP.DmassUmass_INPUTS, rho, e)
+
+ state_vector_struct[EntropicVars.cp.name] = Cp
+ state_vector_struct[EntropicVars.cv.name] = Cv
+ dhdrho_P = dhdrho_e - dhde_rho * (1 / dPde_rho) * dPdrho_e
+ dhdP_rho = dhde_rho * (1 / dPde_rho)
+ dsdrho_P = dsdrho_e - dPdrho_e * (1 / dPde_rho) * dsde_rho
+ dsdP_rho = dsde_rho / dPde_rho
+ state_vector_struct[EntropicVars.dhdrho_p.name] = dhdrho_P
+ state_vector_struct[EntropicVars.dhdp_rho.name] = dhdP_rho
+ state_vector_struct[EntropicVars.dsdrho_p.name] = dsdrho_P
+ state_vector_struct[EntropicVars.dsdp_rho.name] = dsdP_rho
+ return
+
+
+
+ def __EquationofState(self, s:float, rho:float, e:float, derivs:dict, state_vector_struct:dict):
+ """Calculate thermodynamic state variables with entropic equation of state
+
+ :param s: fluid entropy
+ :type s: float
+ :param rho: fluid density
+ :type rho: float
+ :param e: fluid static energy
+ :type e: float
+ :param derivs: entropy Jacobian and Hessian
+ :type derivs: dict
+ :param state_vector_struct: state struct
+ :type state_vector_struct: dict
+ """
+
+ # Retrieve entropy Jacobian and Hessian components
+ dsdrho_e = derivs["dsdrho_e"]
+ dsde_rho = derivs["dsde_rho"]
+ d2sdrho2 = derivs["d2sdrho2"]
+ d2sde2 = derivs["d2sde2"]
+ d2sdedrho = derivs["d2sdedrho"]
+
+ # Store entropy, Jacobian, and Hessian components
+ state_vector_struct[EntropicVars.s.name] = s
+ state_vector_struct[EntropicVars.Density.name] = rho
+ state_vector_struct[EntropicVars.Energy.name] = e
+ state_vector_struct[EntropicVars.dsdrho_e.name] = dsdrho_e
+ state_vector_struct[EntropicVars.dsde_rho.name] = dsde_rho
+ state_vector_struct[EntropicVars.d2sdrho2.name] = d2sdrho2
+ state_vector_struct[EntropicVars.d2sde2.name] = d2sde2
+ state_vector_struct[EntropicVars.d2sdedrho.name] = d2sdedrho
+
+ # Compute and store thermodynamic state variables
+ Temperature = self.fluid.T()#pow(dsde_rho, -1)
+ state_vector_struct[EntropicVars.T.name] = Temperature
+ Pressure = self.fluid.p()#-rho * rho * Temperature * dsdrho_e
+ state_vector_struct[EntropicVars.p.name] = Pressure
+ Enthalpy = self.fluid.hmass()
+ state_vector_struct[EntropicVars.Enthalpy.name] = Enthalpy
+
+ X = self.fluid.Q()
+ if X<=0 or X>=1:
+ state_vector_struct[EntropicVars.dpde_rho.name] = self.fluid.first_partial_deriv(CP.iP, CP.iUmass, CP.iDmass)
+ state_vector_struct[EntropicVars.dpdrho_e.name] = self.fluid.first_partial_deriv(CP.iP, CP.iDmass, CP.iUmass)
+ state_vector_struct[EntropicVars.c2.name] = self.fluid.speed_sound()**2
+ state_vector_struct[EntropicVars.dTde_rho.name] = self.fluid.first_partial_deriv(CP.iT, CP.iUmass, CP.iDmass)
+ state_vector_struct[EntropicVars.dTdrho_e.name] = self.fluid.first_partial_deriv(CP.iT, CP.iDmass, CP.iUmass)
+ state_vector_struct[EntropicVars.cp.name] = self.fluid.cpmass()
+ state_vector_struct[EntropicVars.cv.name] = self.fluid.cvmass()
+ state_vector_struct[EntropicVars.dhdrho_e.name] = self.fluid.first_partial_deriv(CP.iHmass, CP.iDmass, CP.iUmass)
+ state_vector_struct[EntropicVars.dhde_rho.name] = self.fluid.first_partial_deriv(CP.iHmass, CP.iUmass, CP.iDmass)
+ state_vector_struct[EntropicVars.dhdrho_p.name] = self.fluid.first_partial_deriv(CP.iHmass, CP.iDmass, CP.iP)
+ state_vector_struct[EntropicVars.dhdp_rho.name] = self.fluid.first_partial_deriv(CP.iHmass, CP.iP, CP.iDmass)
+ state_vector_struct[EntropicVars.dsdrho_p.name] = self.fluid.first_partial_deriv(CP.iSmass, CP.iDmass, CP.iP)
+ state_vector_struct[EntropicVars.dsdp_rho.name] = self.fluid.first_partial_deriv(CP.iSmass, CP.iP, CP.iDmass)
+ else:
+ self.__EntropicEoS(rho, e, s, derivs, state_vector_struct)
+ return
+
+ def __ThermodynamicState(self):
+ """Evaluate the thermodynamic state variables
+
+ :return: thermodynamic state variables.
+ :rtype: struct
+ """
+ state_vector_vals = {v.name : None for v in EntropicVars}
+
+ phase = self.fluid.phase()
+ rho = self.fluid.rhomass()
+ e = self.fluid.umass()
+ s = self.fluid.smass()
+ is_twophase = (phase == CoolP.iphase_twophase)
+ # Evaluate entropy Jacobian and Hessian components depending on phase
+ derivs = {"dsdrho_e":0,"dsde_rho":0,"d2sdrho2":0,"d2sde2":0,"d2sdedrho":0}
+ if is_twophase:
+ derivs_twophase = self.__compute_derivatives_fd(rho, e, s)
+ derivs["dsdrho_e"] = derivs_twophase["dsdrho_e"]
+ derivs["dsde_rho"] = derivs_twophase["dsde_rho"]
+ derivs["d2sdrho2"] = derivs_twophase["d2sdrho2"]
+ derivs["d2sde2"] = derivs_twophase["d2sde2"]
+ derivs["d2sdedrho"] = derivs_twophase["d2sdedrho"]
+ self.UpdateFluid(rho, e)
+ else:
+ derivs["dsdrho_e"] = self.fluid.first_partial_deriv(CP.iSmass, CP.iDmass, CP.iUmass)
+ derivs["dsde_rho"] = self.fluid.first_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass)
+ derivs["d2sdrho2"] = self.fluid.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iUmass, CP.iDmass, CP.iUmass)
+ derivs["d2sde2"] = self.fluid.second_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass, CP.iUmass, CP.iDmass)
+ derivs["d2sdedrho"] = self.fluid.second_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass, CP.iDmass, CP.iUmass)
+
+ # Calculate thermodynamic state variables with entropy-based equation of state
+ self.__EquationofState(s, rho, e, derivs, state_vector_vals)
+
+ return state_vector_vals
+
+
+ def GetStateVector(self):
+ """Retrieve thermodynamic or thermophysical state information
+
+ :return: array with fluid properties, whether phase is supported
+ :rtype: np.ndarray[float], bool
+ """
state_vector_vals = np.ones(EntropicVars.N_STATE_VARS.value)
- correct_phase = True
- if self.fluid.phase() in self.__accepted_phases:
- state_vector_vals[EntropicVars.s.value] = self.fluid.smass()
- if not self.__mixture:
- state_vector_vals[EntropicVars.dsde_rho.value] = self.fluid.first_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass)
- state_vector_vals[EntropicVars.dsdrho_e.value] = self.fluid.first_partial_deriv(CP.iSmass, CP.iDmass, CP.iUmass)
- state_vector_vals[EntropicVars.d2sde2.value] = self.fluid.second_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass, CP.iUmass, CP.iDmass)
- state_vector_vals[EntropicVars.d2sdedrho.value] = self.fluid.second_partial_deriv(CP.iSmass, CP.iUmass, CP.iDmass, CP.iDmass, CP.iUmass)
- state_vector_vals[EntropicVars.d2sdrho2.value] = self.fluid.second_partial_deriv(CP.iSmass, CP.iDmass, CP.iUmass, CP.iDmass, CP.iUmass)
- state_vector_vals[EntropicVars.Density.value] = self.fluid.rhomass()
- state_vector_vals[EntropicVars.Energy.value] = self.fluid.umass()
- state_vector_vals[EntropicVars.T.value] = self.fluid.T()
- state_vector_vals[EntropicVars.p.value] = self.fluid.p()
- state_vector_vals[EntropicVars.c2.value] = self.fluid.speed_sound()**2
- state_vector_vals[EntropicVars.dTde_rho.value] = self.fluid.first_partial_deriv(CP.iT, CP.iUmass, CP.iDmass)
- state_vector_vals[EntropicVars.dTdrho_e.value] = self.fluid.first_partial_deriv(CP.iT, CP.iDmass, CP.iUmass)
- state_vector_vals[EntropicVars.dpde_rho.value] = self.fluid.first_partial_deriv(CP.iP, CP.iUmass, CP.iDmass)
- state_vector_vals[EntropicVars.dpdrho_e.value] = self.fluid.first_partial_deriv(CP.iP, CP.iDmass, CP.iUmass)
- state_vector_vals[EntropicVars.dhde_rho.value] = self.fluid.first_partial_deriv(CP.iHmass, CP.iUmass, CP.iDmass)
- state_vector_vals[EntropicVars.dhdrho_e.value] = self.fluid.first_partial_deriv(CP.iHmass, CP.iDmass, CP.iUmass)
- state_vector_vals[EntropicVars.dhdp_rho.value] = self.fluid.first_partial_deriv(CP.iHmass, CP.iP, CP.iDmass)
- state_vector_vals[EntropicVars.dhdrho_p.value] = self.fluid.first_partial_deriv(CP.iHmass, CP.iDmass, CP.iP)
- state_vector_vals[EntropicVars.dsdp_rho.value] = self.fluid.first_partial_deriv(CP.iSmass, CP.iP, CP.iDmass)
- state_vector_vals[EntropicVars.dsdrho_p.value] = self.fluid.first_partial_deriv(CP.iSmass, CP.iDmass, CP.iP)
- state_vector_vals[EntropicVars.cp.value] = self.fluid.cpmass()
+
+ phase = self.fluid.phase()
+ q = self.fluid.Q()
+ correct_phase = True
+ if phase in self.__accepted_phases:
+ state_vals = self.__ThermodynamicState()
+ if self._Config.CalcTransportProperties():
+ viscosity, conductivity,vapor_quality = self.__TransportProperties(q)
+ state_vals[EntropicVars.Conductivity.name] = conductivity
+ state_vals[EntropicVars.ViscosityDyn.name] = viscosity
+ state_vals[EntropicVars.VaporQuality.name] = vapor_quality
+ for s in EntropicVars:
+ if s.value < EntropicVars.N_STATE_VARS.value:
+ state_vector_vals[s.value] = state_vals[s.name]
+
else:
correct_phase = False
state_vector_vals[:] = None
return state_vector_vals, correct_phase
+
+ def __get_entropy_safe(self, rho:float, e:float):
+ """Retrieve fluid entropy and check if CoolProp converges
+
+ :param rho: fluid density
+ :type rho: float
+ :param e: fluid static energy
+ :type e: float
+ :return: fluid entropy
+ :rtype: float
+ """
+ try:
+ self.fluid.update(CP.DmassUmass_INPUTS, rho, e)
+ phase = self.fluid.phase()
+ if phase not in self.__accepted_phases:
+ return None
+ return self.fluid.smass()
+ except:
+ return None
+
+ def SetFDStepSizes(self, fd_step_rho:float=1e-5, fd_step_e:float=1e-5):
+ """Set relative step sizes for density and static energy used for finite-difference analysis.
+
+ :param fd_step_rho: relative density FD step, defaults to 1e-5
+ :type fd_step_rho: float, optional
+ :param fd_step_e: relative static energy FD step, defaults to 1e-5
+ :type fd_step_e: float, optional
+ :raises Exception: if negative step sizes are provided
+ """
+ if fd_step_rho <= 0 or fd_step_e <= 0:
+ raise Exception("Relative step sizes for finite-differences should be positive")
+ self.__fd_step_size_rho = fd_step_rho
+ self.__fd_step_size_e = fd_step_e
+ return
+
+ def ComputeSaturationCurve(self):
+ """Compute the density and static energy along the fluid saturation curve
+
+ :return: density and static energy array of saturation curve
+ :rtype: np.ndarray[float]
+ """
+ eos_fluid = "%s::%s" % (self._Config.GetEquationOfState(), self._Config.GetFluid())
+ ptriple = CP.PropsSI("PTRIPLE", eos_fluid)
+ pcrit = CP.PropsSI("PCRIT", eos_fluid)
+ Psat = np.linspace(ptriple, pcrit, 2000)
+
+ rhoLiq=CP.PropsSI("D","P",Psat,"Q",0,eos_fluid)
+ rhoVap=CP.PropsSI("D","P",Psat,"Q",1,eos_fluid)
+
+ eLiq=CP.PropsSI("U","P",Psat,"Q",0,eos_fluid)
+ eVap=CP.PropsSI("U","P",Psat,"Q",1,eos_fluid)
+
+ rho_sat=np.concatenate((rhoLiq[:-1],np.flip(rhoVap)))
+ e_sat=np.concatenate((eLiq[:-1],np.flip(eVap)))
+ sat_curve=np.column_stack((rho_sat,e_sat))
+ return sat_curve
+
+ def GetFDStepSizes(self):
+ return self.__fd_step_size_rho, self.__fd_step_size_e
+
+ def __get_saturated_transport_properties(self, p:float):
+ """Get transport properties at saturation conditions for both liquid and vapor phases.
+
+ :param p: pressure [Pa]
+ :type p: float
+ :return: dict with keys: 'mu_l', 'mu_g', 'k_l', 'k_g', 'rho_l', 'rho_g'
+ :rtype: dict
+ """
+ try:
+ # Create temporary fluid states for saturated liquid and vapor
+ # Use pressure-quality inputs to get saturation properties
+
+ # Saturated liquid (Q=0)
+ self.fluid.update(CP.PQ_INPUTS, p, 0.0)
+ mu_l = self.fluid.viscosity()
+ k_l = self.fluid.conductivity()
+ rho_l = self.fluid.rhomass()
+
+ # Saturated vapor (Q=1)
+ self.fluid.update(CP.PQ_INPUTS, p, 1.0)
+ mu_g = self.fluid.viscosity()
+ k_g = self.fluid.conductivity()
+ rho_g = self.fluid.rhomass()
+
+ return {
+ 'mu_l': mu_l, 'mu_g': mu_g,
+ 'k_l': k_l, 'k_g': k_g,
+ 'rho_l': rho_l, 'rho_g': rho_g
+ }
+ except Exception as ex:
+ return None
+
+ def __compute_void_fraction(self, quality:float, rho_l:float, rho_g:float):
+ """Compute void fractio from quality using the homogeneous model.
+
+ :param quality: vapor quality (mass fraction of vapor)
+ :type quality: float
+ :param rho_l: Saturated liquid density [kg/m3]
+ :type rho_l: float
+ :param rho_g: Saturated vapor density [kg/m3]
+ :type rho_g: float
+ :return: Void fraction (volume fraction of vapor)
+ :rtype: float
+ """
+
+ if quality <= 0:
+ return 0.0
+ elif quality >= 1:
+ return 1.0
+ else:
+ # Homogeneous model void fraction
+ return 1.0 / (1.0 + (1.0 - quality) / quality * rho_g / rho_l)
+
+
+ def __compute_twophase_viscosity(self, quality:float, mu_l:float, mu_g:float, rho_l:float=None, rho_g:float=None, alpha:float=None):
+ """Compute two-phase mixture viscosity using selected mixing model.
+
+ :param quality: Vapor quality (mass fraction)
+ :type quality: float
+ :param mu_l: Saturated liquid viscosity [Pas]
+ :type mu_l: float
+ :param mu_g: Saturated vapor viscosity [Pas]
+ :type mu_g: float
+ :param rho_l: Saturated liquid density [kg/m3] (needed for some models), defaults to None
+ :type rho_l: float, optional
+ :param rho_g: Saturated vapor density [kg/m3] (needed for some models), defaults to None
+ :type rho_g: float, optional
+ :param alpha: Void fraction, defaults to None
+ :type alpha: float, optional
+ :return: Two-phase mixture viscosity [Pas]
+ :rtype: float
+ """
+
+ x = quality
+ if self._Config.GetViscosityModel() == "mcadams":
+ # McAdams et al. (1942) - Reciprocal average (recommended for most cases)
+ if x <= 0:
+ return mu_l
+ elif x >= 1:
+ return mu_g
+ else:
+ return 1.0 / (x / mu_g + (1.0 - x) / mu_l)
+
+ elif self._Config.GetViscosityModel() == "cicchitti":
+ # Cicchitti et al. (1960) - Mass-weighted average
+ return x * mu_g + (1.0 - x) * mu_l
+
+ elif self._Config.GetViscosityModel() == "dukler":
+ # Dukler et al. (1964) - Volume-weighted (needs void fraction)
+ if alpha is None:
+ alpha = self.__compute_void_fraction(x, rho_l, rho_g)
+ return alpha * mu_g + (1.0 - alpha) * mu_l
+
+ else:
+ # Default to McAdams
+ if x <= 0:
+ return mu_l
+ elif x >= 1:
+ return mu_g
+ else:
+ return 1.0 / (x / mu_g + (1.0 - x) / mu_l)
+
+ def __compute_twophase_conductivity(self, quality:float, k_l:float, k_g:float, rho_l:float=None, rho_g:float=None, alpha:float=None):
+ """Compute two-phase mixture thermal conductivity using selected mixing model.
+
+ :param quality: Vapor quality (mass fraction)
+ :type quality: float
+ :param k_l: Saturated liquid thermal conductivity [W/(mK)]
+ :type k_l: float
+ :param k_g: Saturated vapor thermal conductivity [W/(mK)]
+ :type k_g: float
+ :param rho_l: Saturated liquid density [kg/m3], defaults to None
+ :type rho_l: float, optional
+ :param rho_g: Saturated vapor density [kg/m3], defaults to None
+ :type rho_g: float, optional
+ :param alpha: Void fraction, defaults to None
+ :type alpha: float, optional
+ :return: Two-phase mixture thermal conductivity [W/(mK)]
+ :rtype: float
+ """
+
+ x = quality
+
+ if alpha is None and rho_l is not None and rho_g is not None:
+ alpha = self.__compute_void_fraction(x, rho_l, rho_g)
+
+ if self._Config.GetConductivityModel() == "volume":
+ # Volume-weighted average (recommended)
+ if alpha is None:
+ # Fallback to mass-weighted if void fraction unavailable
+ return x * k_g + (1.0 - x) * k_l
+ return alpha * k_g + (1.0 - alpha) * k_l
+
+ elif self._Config.GetConductivityModel() == "mass":
+ # Mass-weighted average
+ return x * k_g + (1.0 - x) * k_l
+
+ else:
+ # Default to volume-weighted
+ if alpha is None:
+ return x * k_g + (1.0 - x) * k_l
+ return alpha * k_g + (1.0 - alpha) * k_l
+
+
+ def DiscretizationError(self, rho, e):
+ self.UpdateFluid(rho, e)
+ state_vector,_ = self.GetStateVector()
+ self.UpdateFluid(rho, e)
+ dPde_rho_ref = self.fluid.first_partial_deriv(CP.iP, CP.iUmass, CP.iDmass)
+ dPdrho_e_ref = self.fluid.first_partial_deriv(CP.iP, CP.iDmass, CP.iUmass)
+ dTde_rho_ref = self.fluid.first_partial_deriv(CP.iT, CP.iUmass, CP.iDmass)
+ dTdrho_e_ref = self.fluid.first_partial_deriv(CP.iT, CP.iDmass, CP.iUmass)
+ P_ref = self.fluid.p()
+ T_ref = self.fluid.T()
+ ref_data = np.array([P_ref, T_ref])
+
+ dPde_rho_test = state_vector[EntropicVars.dpde_rho.value]
+ dPdrho_e_test = state_vector[EntropicVars.dpdrho_e.value]
+ dTde_rho_test = state_vector[EntropicVars.dTde_rho.value]
+ dTdrho_e_test = state_vector[EntropicVars.dTdrho_e.value]
+ P_test = state_vector[EntropicVars.p.value]
+ T_test = state_vector[EntropicVars.T.value]
+ test_data = np.array([P_test, T_test])
+
+ e = np.sum(np.power((test_data - ref_data)/test_data,2))
+ return e
+
+ def __compute_derivatives_fd(self, rho:float, e:float, s_center:float):
+ """Compute all entropy derivatives using central finite differences. Works for both single-phase and two-phase regions.
+
+ :param rho: fluid density [kg/m3]
+ :type rho: float
+ :param e: fluid static energy [J/kg]
+ :type e: float
+ :param s_center: fluid entropy [J/kg]
+ :type s_center: float
+ :return: dict with keys: 'dsdrho_e', 'dsde_rho', 'd2sdrho2', 'd2sde2', 'd2sdedrho'
+ :rtype: dict
+ """
+
+ # Compute step sizes
+ drho = rho * self.__fd_step_size_rho#max(rho * self.__fd_step_size_rho, 1e-6) # Minimum absolute step
+ de = abs(e) * self.__fd_step_size_e#max(abs(e) * self.__fd_step_size_e, 1.0) # Minimum 1 J/kg step
+
+ # Get entropy at stencil points for first derivatives (central difference)
+ s_rho_plus = self.__get_entropy_safe(rho + drho, e)
+ s_rho_minus = self.__get_entropy_safe(rho - drho, e)
+ s_e_plus = self.__get_entropy_safe(rho, e + de)
+ s_e_minus = self.__get_entropy_safe(rho, e - de)
+
+ # Check if first derivative stencil is valid
+ if any(s is None for s in [s_rho_plus, s_rho_minus, s_e_plus, s_e_minus]):
+ # Try one-sided differences as fallback
+ return self.__compute_derivatives_fd_onesided(rho, e, s_center, drho, de)
+
+ # First derivatives (central difference: O(h^2) accuracy)
+ dsdrho_e = (s_rho_plus - s_rho_minus) / (2 * drho)
+ dsde_rho = (s_e_plus - s_e_minus) / (2 * de)
+
+ # Second derivatives
+ d2sdrho2 = (s_rho_plus - 2*s_center + s_rho_minus) / (drho**2)
+ d2sde2 = (s_e_plus - 2*s_center + s_e_minus) / (de**2)
+
+ # Mixed derivative d2s/drhode
+ # Use central difference: [s(rho+,e+) - s(rho+,e-) - s(rho-,e+) + s(rho-,e-)] / (4 * drho * de)
+ s_pp = self.__get_entropy_safe(rho + drho, e + de)
+ s_pm = self.__get_entropy_safe(rho + drho, e - de)
+ s_mp = self.__get_entropy_safe(rho - drho, e + de)
+ s_mm = self.__get_entropy_safe(rho - drho, e - de)
+
+ if any(s is None for s in [s_pp, s_pm, s_mp, s_mm]):
+ # Mixed derivative failed - use simpler approximation
+ d2sdedrho = 0.0 # This is less critical than first derivatives
+ else:
+ d2sdedrho = (s_pp - s_pm - s_mp + s_mm) / (4 * drho * de)
+
+ return {
+ 'dsdrho_e': dsdrho_e,
+ 'dsde_rho': dsde_rho,
+ 'd2sdrho2': d2sdrho2,
+ 'd2sde2': d2sde2,
+ 'd2sdedrho': d2sdedrho
+ }
+
+ def __compute_derivatives_fd_onesided(self, rho:float, e:float, s_center:float, drho:float, de:float):
+ """Fallback: compute derivatives using one-sided finite differences. Less accurate but more robust near boundaries.
+
+ :param rho: fluid density [kg/m3]
+ :type rho: float
+ :param e: fluid static energy [J/kg]
+ :type e: float
+ :param s_center: fluid entropy [J/kg]
+ :type s_center: float
+ :param drho: density increment [kg/m3]
+ :type drho: float
+ :param de: static energy increment [J/kg]
+ :type de: float
+ :return: dict with keys: 'dsdrho_e', 'dsde_rho', 'd2sdrho2', 'd2sde2', 'd2sdedrho'
+ :rtype: dict
+ """
+
+ # Fallback: compute derivatives using one-sided finite differences.
+ # Less accurate but more robust near boundaries.
+
+ derivs = {}
+
+ # Try forward difference for ds/drho
+ s_rho_plus = self.__get_entropy_safe(rho + drho, e)
+ s_rho_minus = self.__get_entropy_safe(rho - drho, e)
+
+ if s_rho_plus is not None and s_rho_minus is not None:
+ derivs['dsdrho_e'] = (s_rho_plus - s_rho_minus) / (2 * drho)
+ derivs['d2sdrho2'] = (s_rho_plus - 2*s_center + s_rho_minus) / (drho**2)
+ elif s_rho_plus is not None:
+ derivs['dsdrho_e'] = (s_rho_plus - s_center) / drho
+ derivs['d2sdrho2'] = 0.0
+ elif s_rho_minus is not None:
+ derivs['dsdrho_e'] = (s_center - s_rho_minus) / drho
+ derivs['d2sdrho2'] = 0.0
+ else:
+ return None # Cannot compute
+
+ # Try for ds/de
+ s_e_plus = self.__get_entropy_safe(rho, e + de)
+ s_e_minus = self.__get_entropy_safe(rho, e - de)
+
+ if s_e_plus is not None and s_e_minus is not None:
+ derivs['dsde_rho'] = (s_e_plus - s_e_minus) / (2 * de)
+ derivs['d2sde2'] = (s_e_plus - 2*s_center + s_e_minus) / (de**2)
+ elif s_e_plus is not None:
+ derivs['dsde_rho'] = (s_e_plus - s_center) / de
+ derivs['d2sde2'] = 0.0
+ elif s_e_minus is not None:
+ derivs['dsde_rho'] = (s_center - s_e_minus) / de
+ derivs['d2sde2'] = 0.0
+ else:
+ return None
+
+ # Mixed derivative - just set to zero if we're already using fallback
+ derivs['d2sdedrho'] = 0.0
+
+ return derivs
+
+ def UpdateFluid(self, val_rho:float, val_e:float):
+ """Update the CoolProp fluid object with a given density and static energy.
+
+ :param val_rho: fluid density [kg/m3]
+ :type val_rho: float
+ :param val_e: fluid static energy [J/kg]
+ :type val_e: float
+ """
+
+ self.fluid.update(CP.DmassUmass_INPUTS, val_rho, val_e)
+ return
+
def VisualizeFluidData(self):
"""Visualize computed fluid data.
"""
diff --git a/Documentation/source/Config.rst b/Documentation/source/Config.rst
deleted file mode 100644
index 234d547..0000000
--- a/Documentation/source/Config.rst
+++ /dev/null
@@ -1,239 +0,0 @@
-.. _configurations:
-
-***************************
-SU2 DataMiner Configuration
-***************************
-*SU2 DataMiner* uses a configuration class in order to store important information regarding the data generation, data mining, and
-manifold generation processes. This page lists some of the important functions of the *Config* class which acts as the :ref:`base`
-class for configurations specific to the application such as :ref:`NICFD` and FGM :ref:`FGM`, for which additional settings can be specified.
-
-.. contents:: :depth: 2
-
-.. _base:
-
-Base Config Class
-=================
-
-*SU2 DataMiner* uses a configuration class in order to store important information regarding the data generation, data mining, and
-manifold generation processes. This page lists some of the important functions of the *Config* class which acts as the base
-class for configurations specific to the application such as :ref:`NICFD` and FGM :ref:`FGM`, for which additional settings can be specified.
-
-
-Storage Location and Configuration Information
-----------------------------------------------
-
-During the various processes in *SU2 DataMiner*, data are generated, processed, and analyzed. All information regarding these
-processes is stored on the current hardware in a user-defined location. *SU2 DataMiner* configurations can be saved locally under
-different names in order to keep track of various data sets and manifolds at once.
-The following functions can be used to manipulate and access the storage location for fluid data and manifolds of the *SU2 DataMiner* configuration
-and save and load configurations.
-
-.. autofunction:: Common.Config_base.Config.SetConfigName
-
-.. autofunction:: Common.Config_base.Config.GetConfigName
-
-.. autofunction:: Common.Config_base.Config.SaveConfig
-
-
-.. code-block::
-
- from su2dataminer.config import Config
-
- c = Config()
- c.SetConfigName("test")
- c.SaveConfig()
-
-
-.. autofunction:: Common.Config_base.Config.SetOutputDir
-
-.. autofunction:: Common.Config_base.Config.GetOutputDir
-
-.. autofunction:: Common.Config_base.Config.PrintBanner
-
-.. autofunction:: Common.Config_base.Config.SetConcatenationFileHeader
-
-.. autofunction:: Common.Config_base.Config.GetConcatenationFileHeader
-
-
-Settings for Machine Learning Applications
-------------------------------------------
-
-The data-driven fluid modeling applications of *SU2 DataMiner* involve the use of multi-layer perceptrons to calculate the thermodynamic state of the fluid during flow calculations in SU2.
-The values of the weights and biases, the hidden layer architecture(s) and other parameters needed to train the network can be stored in and retrieved from the *SU2 DataMiner* configuration.
-
-.. autofunction:: Common.Config_base.Config.SetControllingVariables
-
-.. autofunction:: Common.Config_base.Config.GetControllingVariables
-
-.. autofunction:: Common.Config_base.Config.SetTrainFraction
-
-.. autofunction:: Common.Config_base.Config.GetTrainFraction
-
-.. autofunction:: Common.Config_base.Config.SetTestFraction
-
-.. autofunction:: Common.Config_base.Config.GetTestFraction
-
-The following functions are used to specify the parameters used for training the networks on the fluid data. Currently, *SU2 DataMiner* uses supervised, gradient-based methods for training where the learning
-rate follows the exponential decay method. The value of the initial learning rate and decay parameter can be accessed through
-
-.. autofunction:: Common.Config_base.Config.SetAlphaExpo
-
-.. autofunction:: Common.Config_base.Config.GetAlphaExpo
-
-.. autofunction:: Common.Config_base.Config.SetLRDecay
-
-.. autofunction:: Common.Config_base.Config.GetLRDecay
-
-The networks are trained using batches of training data. The size of the training batches and the maximum number of epochs for which the networks are trained can be accessed through
-
-.. autofunction:: Common.Config_base.Config.SetBatchExpo
-
-.. autofunction:: Common.Config_base.Config.GetBatchExpo
-
-.. autofunction:: Common.Config_base.Config.SetNEpochs
-
-.. autofunction:: Common.Config_base.Config.GetNEpochs
-
-Additionally, the number of nodes in the hidden layers of the networks and activation function can be accessed. *SU2 DataMiner* currently supports the use of a single activation function to all hidden layers which can be selected from a list of options and the input and output layers use a linear function.
-
-.. autofunction:: Common.Config_base.Config.SetHiddenLayerArchitecture
-
-.. autofunction:: Common.Config_base.Config.GetHiddenLayerArchitecture
-
-.. autofunction:: Common.Config_base.Config.SetActivationFunction
-
-.. autofunction:: Common.Config_base.Config.GetActivationFunction
-
-Finally, the weights and biases data can be accessed in the configuration using
-
-.. autofunction:: Common.Config_base.Config.SetWeights
-
-.. autofunction:: Common.Config_base.Config.SetBiases
-
-.. autofunction:: Common.Config_base.Config.GetWeightsBiases
-
-All the settings regarding the training parameters, architecture, weights, and biases of the trained network can be stored automatically after training a network from the :ref:`Trainer` class
-
-.. autofunction:: Common.Config_base.Config.UpdateMLPHyperParams
-
-
-To use the network stored in the configuration in SU2 simulations, the network information needs to be written to an ASCII file such that it can be loaded in SU2 through the `MLPCpp`_ module.
-All relevant information about the network is automatically written to a properly formatted ASCII file using
-
-.. autofunction:: Common.Config_base.Config.WriteSU2MLP
-
-.. _NICFD:
-
-Configuration for real-gas applications
-=======================================
-
-The following section describes the most important functions of the Config_NICFD class.
-
-.. autofunction:: Common.DataDrivenConfig.Config_NICFD.__init__
-
-.. autofunction:: Common.DataDrivenConfig.Config_NICFD.SetFluid
-
-.. autofunction:: Common.DataDrivenConfig.Config_NICFD.SetEquationOfState
-
-Example
--------
-The code snippet below demonstrates the
-::
-
- from su2dataminer.config import Config_NICFD
-
- config = Config_NICFD()
- config.SetFluid("MM")
- config.SetEquationOfState("HEOS")
- config.SetName("siloxane_MM_heos")
- config.UseAutoRange(True)
- config.PrintBanner()
-
-
-.. _FGM:
-
-Configuration for combustion applications
-=========================================
-
-The following section describes the most important functions of the SU2 DataMiner configuration class for FGM applications.
-
-SU2 DataMiner supports the generation of premixed flamelets using `Cantera `_ for FGM applications.
-
-
-Reaction Mechanism and Species Transport Model
-----------------------------------------------
-
-
-The reaction mechanism used for flamelet calculations is specified through the following function:
-
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetReactionMechanism
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.GetReactionMechanism
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetTransportModel
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.GetTransportModel
-
-
-Reactants and Flamelet Types
-----------------------------
-
-
-The fuel and oxidizer used for flamelet calculations can be defined through the following functions:
-
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetFuelDefinition
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetOxidizerDefinition
-
-
-The types of flamelet data that are currently supported by SU2 DataMiner are adiabatic flamelets, burner-stabilized flamelets, and chemical equilibrium data.
-Each of these flamelet types can be included or excluded from the manifold through the following functions:
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.RunFreeFlames
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.RunBurnerFlames
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.RunEquilibrium
-
-
-FGM Controlling Variables
--------------------------
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetControllingVariables
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.GetControllingVariables
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetProgressVariableDefinition
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.ResetProgressVariableDefinition
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.GetUnburntScalars
-
-
-Preferential Diffusion
-----------------------
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.EnablePreferentialDiffusion
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.SetAverageLewisNumbers
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.ComputeBetaTerms
-
-
-Thermochemical State Variables
-------------------------------
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.AddOutputGroup
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.RemoveOutputGroup
-
-
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.EditOutputGroup
-
-
-
-.. _MLPCpp : https://github.com/EvertBunschoten/MLPCpp.git
-
-.. _CanteraLink: https://cantera.org/
\ No newline at end of file
diff --git a/Documentation/source/Setup.rst b/Documentation/source/Setup.rst
index bceaabb..7a12b0a 100644
--- a/Documentation/source/Setup.rst
+++ b/Documentation/source/Setup.rst
@@ -3,7 +3,7 @@
Set-up and Requirements
=======================
-*SU2 DataMiner* is an open-source, python-based software suite that is available on `GitHub `_.
+*SU2 DataMiner* is an open-source, python-based software suite that is available on `GitHub `_.
.. important::
@@ -32,7 +32,7 @@ The *SU2 DataMiner* source code can be downloaded from GitHub by cloning the rep
.. code-block::
- >>> git clone https://github.com/Propulsion-Power-TU-Delft/SU2_DataMiner.git
+ >>> git clone https://github.com/su2code/SU2_DataMiner.git
where `` refers to the target location where the *SU2 DataMiner* source code will be stored.
diff --git a/Documentation/source/Tutorials/FGM/index.rst b/Documentation/source/Tutorials/FGM/index.rst
new file mode 100644
index 0000000..c54899a
--- /dev/null
+++ b/Documentation/source/Tutorials/FGM/index.rst
@@ -0,0 +1,7 @@
+Tutorials for Flamelet-Generated Manifolds
+==========================================
+
+.. toctree::
+ methane_flamelets.rst
+ :maxdepth: 1
+
diff --git a/Documentation/source/Tutorials/FGM/methane_flamelets.rst b/Documentation/source/Tutorials/FGM/methane_flamelets.rst
new file mode 100644
index 0000000..c9a24c8
--- /dev/null
+++ b/Documentation/source/Tutorials/FGM/methane_flamelets.rst
@@ -0,0 +1,18 @@
+.. sectionauthor:: Evert Bunschoten
+
+.. _tutorial_methane_flamelets:
+
+
+
+Generating a methane flamelet manifold
+======================================
+
+This tutorial shows how to generate a manifold of methane-air flamelets using *SU2 DataMiner*.
+The manifold contains adiabatic flamelets, burner-stabilized flamelets, and chemical equilibrium data generated for equivalence ratios between 0.8 and 1.2.
+
+
+.. contents:: :depth: 2
+
+
+Set-up
+------
diff --git a/Documentation/source/Tutorials/NICFD/index.rst b/Documentation/source/Tutorials/NICFD/index.rst
new file mode 100644
index 0000000..3ced445
--- /dev/null
+++ b/Documentation/source/Tutorials/NICFD/index.rst
@@ -0,0 +1,7 @@
+Tutorials for Non-Ideal Compressible Fluid Dynamics
+===================================================
+
+.. toctree::
+ tablegeneration
+ :maxdepth: 1
+
diff --git a/Documentation/source/Tutorials/NICFD/tablegeneration.rst b/Documentation/source/Tutorials/NICFD/tablegeneration.rst
new file mode 100644
index 0000000..f6abad3
--- /dev/null
+++ b/Documentation/source/Tutorials/NICFD/tablegeneration.rst
@@ -0,0 +1,64 @@
+.. _NICFD_LUT:
+
+.. sectionauthor:: Evert Bunschoten
+
+Table Generation for NICFD Applications
+=======================================
+
+SU2 DataMiner supports the creation of look-up table methods for thermophyscial state evaluations in NICFD simulations in SU2.
+This tutorial demonstrates how to generate thermophyisical tables for NICFD applications in SU2.
+
+To get started, you will need to have installed SU2 DataMiner according to the :ref:`installation instructions `.
+
+
+.. contents:: :depth: 2
+
+
+1. Config Generation
+--------------------
+
+As for any process within the SU2 DataMiner workflow, all settings regarding the setup of the fluid data generation and tabulation are stored in an SU2 DataMiner :ref:`configuration object `.
+The tutorial for setting up a generic SU2 DataMiner configuration can be found :ref:`here `.
+
+In this example, a look-up table will be created for the application of modeling fluid properties of carbondioxide in supercritical conditions.
+The following Python code snippet shows the initial set-up.
+
+
+.. code-block::
+
+ #!/usr/bin/env python3
+ from su2dataminer.config import Config_NICFD
+
+ config = Config_NICFD()
+ config.SetFluid("CarbonDioxide")
+ config.SetEquationOfState("HEOS")
+
+
+
+
+
+
+Step 2:
+-------
+
+
+
+.. image:: /Tutorials/SU2DataMiner_logo.png
+ :height: 200 px
+ :width: 200 px
+ :scale: 50 %
+ :alt: this is a detailed caption of the image
+ :align: left
+ :loading: embed
+
+
+and some text here
+
+Step 3:
+-------
+
+
+.. _literature_link_1:
+
+
+.. _literature_link_2:
diff --git a/Documentation/source/Tutorials/SU2DataMiner_logo.png b/Documentation/source/Tutorials/SU2DataMiner_logo.png
new file mode 100644
index 0000000..dff639e
Binary files /dev/null and b/Documentation/source/Tutorials/SU2DataMiner_logo.png differ
diff --git a/Documentation/source/Tutorials/create_configs.rst b/Documentation/source/Tutorials/create_configs.rst
index 4ebe7f6..b0b166e 100644
--- a/Documentation/source/Tutorials/create_configs.rst
+++ b/Documentation/source/Tutorials/create_configs.rst
@@ -1,3 +1,7 @@
+.. _tutorialconfigs:
+
+.. sectionauthor:: Evert Bunschoten
+
Creating SU2 DataMiner Configurations
=====================================
diff --git a/Documentation/source/Tutorials/tutorial_template.rst b/Documentation/source/Tutorials/tutorial_template.rst
new file mode 100644
index 0000000..22d61de
--- /dev/null
+++ b/Documentation/source/Tutorials/tutorial_template.rst
@@ -0,0 +1,58 @@
+.. _tutorial_label:
+
+.. sectionauthor:: your name
+
+Tutorial Template
+=================
+
+
+Use this file as a template for writing your own tutorial when creating new functionalities in SU2 DataMiner. Copy this file to the folder of the respective category and include the file name in the toctree of the ``index.rst`` file.
+
+Introduction text in which you describe the goal of this tutorial. Create hyperlinks to other tutorials and documentation with ``:ref:`this command ``.
+For example, :ref:`this page` leads to the installation instructions page. Consult the `Sphinx documentation page `_ for formatting documentation.
+You can also add images and reference them with :ref:`hyperlinks `.
+
+
+.. contents:: :depth: 2
+
+
+Step 1:
+-------
+
+Some information on how to get started.
+
+.. code-block::
+
+ >>> Explain your setup with code blocks or command line instructions.
+
+
+.. code-block::
+
+ #!/usr/bin/env python3
+ print(Or code snippets.)
+
+
+Step 2:
+-------
+
+
+
+.. image:: /Tutorials/SU2DataMiner_logo.png
+ :height: 200 px
+ :width: 200 px
+ :scale: 50 %
+ :alt: this is a detailed caption of the image
+ :align: left
+ :loading: embed
+
+
+and some text here
+
+Step 3:
+-------
+
+
+.. _literature_link_1:
+
+
+.. _literature_link_2:
diff --git a/Documentation/source/Tutorials/tutorials.rst b/Documentation/source/Tutorials/tutorials.rst
new file mode 100644
index 0000000..18096ed
--- /dev/null
+++ b/Documentation/source/Tutorials/tutorials.rst
@@ -0,0 +1,10 @@
+Tutorials
+=========
+
+.. toctree::
+ create_configs
+ FGM/index
+ NICFD/index
+ tutorial_template
+ :maxdepth: 1
+
diff --git a/Documentation/source/conf.py b/Documentation/source/conf.py
index 34b6780..e3d9eb2 100644
--- a/Documentation/source/conf.py
+++ b/Documentation/source/conf.py
@@ -23,7 +23,7 @@
exclude_patterns = []
add_module_names = False
-
+show_authors=True
numpydoc_show_class_members = False
# -- Options for HTML output -------------------------------------------------
# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output
diff --git a/Documentation/source/documentation/configs/Config_FGM.rst b/Documentation/source/documentation/configs/Config_FGM.rst
index f17a0cf..39639d4 100644
--- a/Documentation/source/documentation/configs/Config_FGM.rst
+++ b/Documentation/source/documentation/configs/Config_FGM.rst
@@ -299,7 +299,7 @@ The composition of the output groups can be manually adapted and accessed with t
.. autofunction:: Common.DataDrivenConfig.Config_FGM.AddOutputGroup
-.. autofunction:: Common.DataDrivenConfig.Config_FGM.DefineOutputGroup
+.. autofunction:: Common.DataDrivenConfig.Config_FGM.EditOutputGroup
.. autofunction:: Common.DataDrivenConfig.Config_FGM.RemoveOutputGroup
diff --git a/Documentation/source/documentation/configs/Config_NICFD.rst b/Documentation/source/documentation/configs/Config_NICFD.rst
index 60c38db..48fa550 100644
--- a/Documentation/source/documentation/configs/Config_NICFD.rst
+++ b/Documentation/source/documentation/configs/Config_NICFD.rst
@@ -1,3 +1,6 @@
+.. sectionauthor:: Evert Bunschoten
+
+
.. _NICFD:
SU2 DataMiner Configuration for NICFD
@@ -35,6 +38,7 @@ retrieved by running the following code snippet:
.. autofunction:: Common.DataDrivenConfig.Config_NICFD.SetFluid
+.. _nicfd_data_limits:
Range of the Thermodynamic State Data
-------------------------------------
@@ -70,6 +74,34 @@ The limits of the fluid data grid are set with the following functions
.. autofunction:: Common.DataDrivenConfig.Config_NICFD.SetDensityBounds
+Phases and transport properties
+-------------------------------
+
+Thermophysical state data can be generated at various phases. The following functions can be used to specify in which phases thermophysical state data are included in the data generation process.
+*SU2 DataMiner* currently only supports the inclusion of fluid data in the gaseous, liquid, two-phase, supercritical, supercritical gas, and supercritical liquid phases.
+Information about the phase data can be accessed with the following functions.
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.EnableGasPhase
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.EnableLiquidPhase
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.EnableSuperCritical
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.GasPhase
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.LiquidPhase
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.SuperCritical
+
+Transport properties can also be included in the fluid data set. Currently, these are limited to the fluid *thermal conductivity*, *dynamic viscosity*, and *vapor quality*.
+Information on whether transport properties are included in the fluid data set can be accessed with the following functions:
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.IncludeTransportProperties
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.CalcTransportProperties
+
+
+.. _resolutionsettings:
Resolution
----------
@@ -116,7 +148,22 @@ By default, these are entropy, pressure, temperature, and the speed of sound, bu
-
+Tabulation
+----------
+
+In addition to the EEoS-PINN method, thermophysical state information can also be stored in look-up tables.
+The thermodynamic state space can be discretized with a *Cartesian* method or with an *adaptive* method.
+With the Cartesian method, the thermodynamic state space is divided into equally sized, rectangular cells in the density-static energy space.
+The number of cells along the density and energy axes are specified with the :ref:`setters ` for thermodynamic state space resolution.
+With the adaptive method, the thermodynamic state space is divided into triangular elements and can be locally refined in areas of interest such as around the saturation curve.
+For both methods, the limits of the table can be specified with the :ref:`setters ` for the density and static energy ranges.
+
+Information on the table discretization method can be accessed with the following functions:
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.SetTableDiscretization
+
+.. autofunction:: Common.DataDrivenConfig.Config_NICFD.GetTableDiscretization
+
References
----------
diff --git a/Documentation/source/documentation/index.rst b/Documentation/source/documentation/index.rst
index ed81d26..04ccd0a 100644
--- a/Documentation/source/documentation/index.rst
+++ b/Documentation/source/documentation/index.rst
@@ -6,5 +6,6 @@ Documentation
datageneration/index.rst
datamining/index.rst
machinelearning/index.rst
+ tabulation/index.rst
:maxdepth: 1
diff --git a/Documentation/source/documentation/tabulation/FGM.rst b/Documentation/source/documentation/tabulation/FGM.rst
new file mode 100644
index 0000000..e69de29
diff --git a/Documentation/source/documentation/tabulation/NICFD.rst b/Documentation/source/documentation/tabulation/NICFD.rst
new file mode 100644
index 0000000..dbe511d
--- /dev/null
+++ b/Documentation/source/documentation/tabulation/NICFD.rst
@@ -0,0 +1,38 @@
+.. _doc_nicfd_tabulation:
+
+.. sectionauthor:: Evert Bunschoten
+
+Tabulation methods for NICFD applications
+=========================================
+
+
+This page documents the tabulation methods for NICFD applications in *SU2 DataMiner*
+
+.. contents:: :depth: 2
+
+
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.__init__
+
+
+Refinement settings
+-------------------
+
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.SetTableDiscretization
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.SetCellSize_Coarse
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.SetCellSize_Refined
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.SetRefinement_Radius
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.AddRefinementCriterion
+
+Table Generation
+----------------
+
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.SetTableVars
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.GenerateTable
+
+Output of tabulation files
+--------------------------
+
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.WriteTableFile
+
+.. autofunction:: Manifold_Generation.LUT.LUTGenerators.SU2TableGenerator_NICFD.WriteOutParaview
+
diff --git a/Documentation/source/documentation/tabulation/index.rst b/Documentation/source/documentation/tabulation/index.rst
new file mode 100644
index 0000000..63e92f8
--- /dev/null
+++ b/Documentation/source/documentation/tabulation/index.rst
@@ -0,0 +1,7 @@
+Tabulation Methods
+==================
+
+.. toctree::
+ NICFD
+ :caption: Tabulation Methods
+ :maxdepth: 1
\ No newline at end of file
diff --git a/Documentation/source/index.rst b/Documentation/source/index.rst
index 1862bbd..85dcb77 100644
--- a/Documentation/source/index.rst
+++ b/Documentation/source/index.rst
@@ -24,7 +24,7 @@ using physics-informed machine learning methods.
.. toctree::
- Tutorials/create_configs.rst
+ Tutorials/tutorials
:maxdepth: 2
:caption: Tutorials:
diff --git a/Manifold_Generation/LUT/LUTGenerators.py b/Manifold_Generation/LUT/LUTGenerators.py
index c7ff8e8..fb2edc6 100644
--- a/Manifold_Generation/LUT/LUTGenerators.py
+++ b/Manifold_Generation/LUT/LUTGenerators.py
@@ -23,302 +23,724 @@
#=============================================================================================#
import numpy as np
-from scipy.spatial import ConvexHull, Delaunay
-from sklearn.preprocessing import MinMaxScaler,RobustScaler,StandardScaler, QuantileTransformer
-import matplotlib.pyplot as plt
+from Common.Properties import EntropicVars,DefaultSettings_NICFD
+from su2dataminer.generate_data import DataGenerator_CoolProp
+from scipy.spatial import Delaunay
+from sklearn.preprocessing import MinMaxScaler
from tqdm import tqdm
-import sys,os
from Common.DataDrivenConfig import Config_NICFD
-import cantera as ct
import gmsh
-import pickle
-from multiprocessing import Pool
-from sklearn.metrics import mean_squared_error
-from Common.Interpolators import Invdisttree
-from random import sample
-
-class SU2TableGenerator:
+from concave_hull import concave_hull, concave_hull_indexes
+import meshio
+import matplotlib.pyplot as plt
+
+
+def FiniteDifferenceDerivative(y:np.ndarray[float], x:np.ndarray[float]):
+ """Calculate second-order accurate, one-dimensional finite-difference derivatives of y with respect to x.
+
+ :param y: data to calculate the finite-differences for.
+ :type y: np.ndarray[float]
+ :param x: axial coordinates.
+ :type x: np.ndarray[float]
+ :return: finite-difference derivatives of y with respect to x.
+ :rtype: np.ndarray[float]
+ """
+ Np = len(x)
+ dydx = np.zeros(Np)
+ for i in range(1, Np-1):
+ y_m = y[i-1]
+ y_p = y[i+1]
+ y_0 = y[i]
+ x_m = x[i-1]
+ x_p = x[i+1]
+ x_0 = x[i]
+ dx_1 = x_p - x_0
+ dx_2 = x_0 - x_m
+ dx2_1 = dx_1*dx_1
+ dx2_2 = dx_2*dx_2
+ if (dx_1==0) or (dx_2==0):
+ dydx[i] = 0.0
+ else:
+ dydx[i] = (dx2_2 * y_p + (dx2_1 - dx2_2)*y_0 - dx2_1*y_m)/(dx_1*dx_2*(dx_1+dx_2))
+ dx_1 = x[1] - x[0]
+ dx_2 = x[2] - x[0]
+ dx2_1 = dx_1*dx_1
+ dx2_2 = dx_2*dx_2
+ y_0 = y[0]
+ y_p = y[1]
+ y_pp = y[2]
+ if (dx_1==0) or (dx_2==0):
+ dydx[0] = 0.0
+ else:
+ dydx[0] = (dx2_1 * y_pp + (dx2_2 - dx2_1)*y_0 - dx2_2*y_p)/(dx_1*dx_2*(dx_1 - dx_2))
+
+ dx_1 = x[-2] - x[-1]
+ dx_2 = x[-3] - x[-1]
+ dx2_1 = dx_1*dx_1
+ dx2_2 = dx_2*dx_2
+ y_0 = y[-1]
+ y_p = y[-2]
+ y_pp = y[-3]
+ if (dx_1==0) or (dx_2==0):
+ dydx[-1] = 0.0
+ else:
+ dydx[-1] = (dx2_1 * y_pp + (dx2_2 - dx2_1)*y_0 - dx2_2*y_p)/(dx_1*dx_2*(dx_1 - dx_2))
+ return dydx
+
+class SU2TableGenerator_NICFD:
_Config:Config_NICFD = None # Config_FGM class from which to read settings.
+ _DataGenerator:DataGenerator_CoolProp = None
_savedir:str
- _mixfrac_min:float = None # Minimum mixture fraction value of the flamelet data.
- _mixfrac_max:float = None # Maximum mixture fraction value of the flamelet data.
-
- _pv_full_norm:np.ndarray[float] = None # Normalized progress variable values of the flamelet data.
- _enth_full_norm:np.ndarray[float] = None # Normalized total enthalpy values of the flamelet data.
- _mixfrac_full_norm:np.ndarray[float] = None # Normalized mixture fraction values of the flamelet data.
-
- _Fluid_Variables:list[str] = None # Variable names in the concatenated flamelet data file.
- _Flamelet_Data:np.ndarray[float] = None # Concatenated flamelet data.
-
- _custom_table_limits_set:bool = False
- _mixfrac_min_table:float = None # Lower mixture fraction limit of the table.
- _mixfrac_max_table:float = None # Upper mixture fraction limit of the table.
-
- __run_parallel:bool = False
- __Np_cores:int = 1
-
- _N_table_levels:int = 100 # Number of table levels.
- _mixfrac_range_table:np.ndarray[float] = None # Mixture fraction values of the table levels.
- _base_cell_size:float = 1e-2#3.7e-3 # Table level base cell size.
-
- _refined_cell_size:float = 3e-4#2.5e-3#1.5e-3 # Table level refined cell size.
- _refinement_radius:float = 5e-3#5e-2 # Table level radius within which refinement is applied.
- _curvature_threshold:float = 0.15 # Curvature threshold above which refinement is applied.
+ _base_cell_size:float = 2e-2 # Table level base cell size.
+ _refined_cell_size:float = 5e-3#2.5e-3#1.5e-3 # Table level refined cell size.
+ _refinement_radius:float = 1e-2#5e-2 # Table level radius within which refinement is applied.
+ __refinement_vars:list[str] = []
+ __refinement_norm_min:list[float] = []
+ __refinement_norm_max:list[float] = []
+ _table_vars:list[str] = [s.name for s in EntropicVars][:-1]
_table_nodes = [] # Progress variable, total enthalpy, and mixture fraction node values for each table level.
_table_nodes_norm = [] # Normalized table nodes for each level.
_table_connectivity = [] # Table node connectivity per table level.
_table_hullnodes = [] # Hull node indices per table level.
- __table_insert_levels:list[float] = []
_controlling_variables:list[str]=["Density",\
"Energy"] # FGM controlling variables
- _lookup_tree:Invdisttree = None # KD tree with inverse distance weighted interpolation for flamelet data interpolation.
- _flamelet_data_scaler:MinMaxScaler = None # Scaler for flamelet data controlling variables.
- _n_near:int = 9 # Number of nearest neighbors from which to evaluate flamelet data.
- _p_fac:int = 3 # Power by which to weigh distances from query point.
+ _fluid_data_scaler:MinMaxScaler= MinMaxScaler() # Scaler for flamelet data controlling variables.
- def __init__(self, Config:Config_NICFD, load_file:str=None):
+ def __init__(self, Config:Config_NICFD):
"""
- Initiate table generator class.
+ Initiate table generator class. Settings regarding the fluid data generation and table resolution are automatically retrieved from the configuration object.
:param Config: Config_FGM object.
:type Config: Config_FGM
"""
+ self._Config = Config
+ self._controlling_variables= [c for c in self._Config.GetControllingVariables()]
+
+ self._DataGenerator = DataGenerator_CoolProp(self._Config)
+ entropic_vars = [a.name for a in EntropicVars][:-1]
+ self._table_vars = entropic_vars.copy()
+ if not self._Config.TwoPhase():
+ self._table_vars.remove(EntropicVars.VaporQuality.name)
+ if not self._Config.CalcTransportProperties():
+ self._table_vars.remove(EntropicVars.ViscosityDyn.name)
+ self._table_vars.remove(EntropicVars.Conductivity.name)
+ return
+
+ def SetFDStepSize(self, val_step_size:float=3e-7):
+ """Set the relative step size for density and static energy for evaluating fluid properties in the two-phase region.
- if load_file:
- # Load an existing TableGenerator object.
- with open(load_file, "rb") as fid:
- loaded_table_generator = pickle.load(fid)
- self.__dict__ = loaded_table_generator.__dict__.copy()
- else:
- # Create new TableGenerator object.
- self._Config = Config
-
- self.__DefineFluidDataInterpolator()
-
- self._savedir = self._Config.GetOutputDir()
-
- self._table_nodes, self.table_data, self._table_connectivity, self._table_hullnodes = self.__ComputeCurvature()
-
- fig = plt.figure()
- ax = plt.axes(projection='3d')
- ax.plot3D(self._table_nodes[:,0],self._table_nodes[:,1], self.table_data[:, self._Fluid_Variables.index("d2sdrho2")],'k.')
- plt.show()
- self.table_vars = ['s','dsdrho_e','dsde_rho','d2sdrho2','d2sdedrho','d2sde2']
+ :param val_step_size: relative finite-difference step size, defaults to 3e-7
+ :type val_step_size: float, optional
+ :raises Exception: if the provided value is negative or zero.
+ """
+ if val_step_size <= 0:
+ raise Exception("Relative step size for finite-differences should be positive.")
+ self._DataGenerator.SetFDStepSizes(val_step_size,val_step_size)
+ return
- self.WriteTableFile(self._Config.GetOutputDir()+"/LUT_"+self._Config.GetConfigName()+".drg")
+ def SetCellSize_Coarse(self, cell_size_coarse:float=1e-2):
+ """Specify the coarse level cell size of the table
- def __DefineFluidDataInterpolator(self):
- print("Configuring KD-tree for most accurate lookups")
+ :param cell_size_coarse: coarse cell size, defaults to 1e-2
+ :type cell_size_coarse: float, optional
+ :raises Exception: if specified cell size is negative or zero
+ """
+ if cell_size_coarse <= 0:
+ raise Exception("Cell size value should be positive")
+ self._base_cell_size = cell_size_coarse
+ return
+
+ def SetCellSize_Refined(self, cell_size_ref:float=5e-3):
+ """Specify the refined level cell size of the table
- print("Loading fluid data...")
- # Define scaler for FGM controlling variables.
- full_data_file = self._Config.GetOutputDir()+"/"+self._Config.GetConcatenationFileHeader()+"_full.csv"
- with open(full_data_file,'r') as fid:
- self._Fluid_Variables = fid.readline().strip().split(',')
- D_full = np.loadtxt(full_data_file,delimiter=',',skiprows=1)
- self._scaler = MinMaxScaler()
- self.CV_full = np.vstack(tuple(D_full[:, self._Fluid_Variables.index(c)] for c in self._controlling_variables)).T
- self.__min_CV, self.__max_CV = np.min(self.CV_full,axis=0), np.max(self.CV_full,axis=0)
+ :param cell_size_ref: refined cell size, defaults to 1e-2
+ :type cell_size_ref: float, optional
+ :raises Exception: if specified cell size is negative or zero
+ """
+ if cell_size_ref <= 0:
+ raise Exception("Cell size value should be positive")
+ self._refined_cell_size = cell_size_ref
+ return
+
+ def SetRefinement_Radius(self, refinement_radius:float=1e-2):
+ """Specify the radius around each refinement point within which the refined cell size is applied
- CV_full_scaled = self._scaler.fit_transform(self.CV_full)
+ :param refinement_radius: refinement radius, defaults to 1e-2
+ :type refinement_radius: float, optional
+ :raises Exception: if specified value is negative or zero
+ """
+ if refinement_radius <= 0:
+ raise Exception("Refinement radius should be positive")
+ self._refinement_radius = refinement_radius
+ return
+
+ def SetTableDiscretization(self, method:str=DefaultSettings_NICFD.tabulation_method):
+ """Overwrite the thermodynamic state space discretization method from the configuration.
- # Exctract train and test data
- train_data_file = self._Config.GetOutputDir()+"/"+self._Config.GetConcatenationFileHeader()+"_full.csv"
- test_data_file = self._Config.GetOutputDir()+"/"+self._Config.GetConcatenationFileHeader()+"_test.csv"
+ :param method: discratization method, defaults to 'cartesian'
+ :type method: str, optional
+ """
+ self._Config.SetTableDiscretization(method)
+ return
+
+ def SetTableVars(self, table_vars_in:list[str]):
+ """Specify the thermophysical variables to be included in the table file. All quantities are included by default. The list shoud at least contain "Density" and "Energy".
+
+ :param table_vars_in: list with thermophysical variables to be included in the table.
+ :type table_vars_in: list[str]
+ :raises Exception: if any of the specified variables are not supported by SU2 DataMiner.
+ """
+ self._table_vars = []
+ if EntropicVars.Density.name not in table_vars_in:
+ print("Density should always be included in table variables")
+ self._table_vars.append(EntropicVars.Density.name)
+
+ if EntropicVars.Energy.name not in table_vars_in:
+ print("Energy should always be included in table variables")
+ self._table_vars.append(EntropicVars.Energy.name)
+
+ if self._Config.EnableTwophase() and EntropicVars.VaporQuality.name in table_vars_in:
+ print("Table generator not configured for two-phase, ignoring vapor quality from table data.")
+ table_vars_in.remove(EntropicVars.VaporQuality.name)
- var_to_test_for = "d2sdrho2"
+ if not self._Config.CalcTransportProperties():
+ if EntropicVars.Conductivity.name in table_vars_in:
+ print("Table generator not configured for transport properties, ignoring conductivity data")
+ if EntropicVars.ViscosityDyn.name in table_vars_in:
+ print("Table generator not configured for transport properties, ignoring viscosity data")
+
+ valid_vars = True
+ for v in table_vars_in:
+ found_var = False
+ for q in EntropicVars:
+ if v.lower() == q.name.lower():
+ found_var = True
+ self._table_vars.append(q.name)
+ if not found_var:
+ print("Error, \"%s\" is not supported by SU2 DataMiner" % v)
+ valid_vars = False
+ if not valid_vars:
+ raise Exception("Some specified thermophysical variables are not supported.")
+ return
+
+ def __Compute2DMesh(self, points:np.ndarray[float], ref_pts:np.ndarray[float]=[],show:bool=False,sat_curve_pts:np.ndarray[float]=[]):
+ """Populate two-dimensional thermodynamic state space with table nodes according to refinement settings.
+
+ :param points: initial point cloud.
+ :type points: np.ndarray[float]
+ :param ref_pts: locations to apply refinement to, defaults to []
+ :type ref_pts: np.ndarray[float], optional
+ :param show: show the discretization generated by Gmesh, defaults to False
+ :type show: bool, optional
+ :param sat_curve_pts: saturation curve points, defaults to []
+ :type sat_curve_pts: np.ndarray[float], optional
+ :return: table nodes, table connectivity.
+ :rtype: tuple
+ """
+ # Create concave hull of normalized table coordinates.
+ add_sat_curve = (len(sat_curve_pts) > 0)
+ XY_hull = concave_hull(np.unique(points,axis=0), length_threshold=self._base_cell_size)
+ nans = np.isnan(XY_hull).all(1)
+ XY_hull = XY_hull[np.invert(nans), :]
+
+ i = 0
+ hull_indices = [i]
+ while i < (len(XY_hull)-1):
+ i_next = i+1
+ found_next_pt = False
+ while not found_next_pt:
+ dist = np.sqrt(np.sum(np.power(XY_hull[i_next, :] - XY_hull[i, :], 2)))
+ if (dist >= self._base_cell_size) or (i_next == len(XY_hull)-1):
+ found_next_pt = True
+ else:
+ i_next += 1
+ i = i_next
+ hull_indices.append(i_next)
+ XY_hull = XY_hull[hull_indices, :]
+
+ # Initiate gmsh
+ gmsh.initialize()
+ gmsh.model.add("table_level")
+ gmsh.option.setNumber("General.Verbosity", 0)
+ factory = gmsh.model.geo
+
+ # Create hull points
+ hull_pts = []
+ for i in range(int(len(XY_hull))):
+ hull_pts.append(factory.addPoint(XY_hull[i, 0], XY_hull[i, 1], 0, self._base_cell_size))
- D_train = np.loadtxt(train_data_file,delimiter=',',skiprows=1)
- D_test = np.loadtxt(test_data_file,delimiter=',',skiprows=1)
+ # Connect hull points to a closed multi-component curve
+ hull_lines = []
+ for i in range(len(hull_pts)-1):
+ hull_lines.append(factory.addLine(hull_pts[i], hull_pts[i+1]))
+ hull_lines.append(factory.addLine(hull_pts[-1], hull_pts[0]))
+
+ # Create a 2D plane of the enclosed space
+ curvloop = factory.addCurveLoop(hull_lines)
+
+ # Apply refinement points
+ ref_pt_ids = []
+ if len(ref_pts)>0:
+ for i in range(len(ref_pts)):
+ ref_pt_ids.append(factory.addPoint(ref_pts[i,0], ref_pts[i, 1], 0.0))
+
+ factory.synchronize()
- CV_train = np.vstack(tuple(D_train[:, self._Fluid_Variables.index(c)] for c in self._controlling_variables)).T
- CV_test = np.vstack(tuple(D_test[:, self._Fluid_Variables.index(c)] for c in self._controlling_variables)).T
+ if add_sat_curve:
+
+ # Create normal vector to saturation curve.
+ dedrho_sat_norm = FiniteDifferenceDerivative(sat_curve_pts[:,0], sat_curve_pts[:,1])
+ norm_vector = np.column_stack((-1.0 / dedrho_sat_norm, np.ones(len(dedrho_sat_norm))))
+ norm_vector = norm_vector / np.sqrt(np.sum(np.power(norm_vector, 2), axis=1))[:,np.newaxis]
+
+ # Create offset curves to ensure that no nodes are generated on the saturation curve itself.
+ sat_curve_upper_pts = []
+ sat_curve_lower_pts = []
+ dx = 1e-3*self._refined_cell_size
+
+ sat_curve_rhoe_upper = sat_curve_pts + dx * norm_vector
+ sat_curve_rhoe_lower = sat_curve_pts - dx * norm_vector
+ inside_upper = np.logical_and(sat_curve_rhoe_upper > 0, sat_curve_rhoe_upper < 1).all(1)
+ inside_lower = np.logical_and(sat_curve_rhoe_lower > 0, sat_curve_rhoe_lower< 1).all(1)
+ valid_sat_curve_pts = np.logical_and(inside_upper, inside_lower)
+ nans_lower = np.isnan(sat_curve_rhoe_lower).all(1)
+ nans_upper = np.isnan(sat_curve_rhoe_upper).all(1)
+ valid_nans = np.logical_and(np.invert(nans_lower), np.invert(nans_upper))
+ valid_pts = np.logical_and(valid_sat_curve_pts, valid_nans)
+
+ hull_DT = Delaunay(XY_hull)
+ within_hull = np.zeros(len(sat_curve_pts),dtype=np.bool)
+ for i in range(len(sat_curve_pts)):
+ rhoe_sat_upper = sat_curve_rhoe_upper[i,:]
+ within_hull_upper = hull_DT.find_simplex(rhoe_sat_upper)>=0
+ rhoe_sat_lower = sat_curve_rhoe_lower[i,:]
+ within_hull_lower = hull_DT.find_simplex(rhoe_sat_lower)>=0
+ within_hull[i] = (within_hull_upper and within_hull_lower)
+ valid_pts = np.logical_and(valid_pts, within_hull)
+
+ sat_curve_pts = sat_curve_pts[valid_pts, :]
+ norm_vector = norm_vector[valid_pts, :]
+
+ i = 0
+ j = 1
+ sat_curve_upper_pts.append(factory.addPoint(sat_curve_pts[i,0] + dx*norm_vector[i, 0],\
+ sat_curve_pts[i,1] + dx*norm_vector[i, 1],0, 0.6*self._refined_cell_size))
+ sat_curve_lower_pts.append(factory.addPoint(sat_curve_pts[i,0] - dx*norm_vector[i, 0],\
+ sat_curve_pts[i,1] - dx*norm_vector[i, 1],0, 0.6*self._refined_cell_size))
+ dists = []
+ while j < len(sat_curve_pts):
+ dist = np.sqrt(np.sum(np.power(sat_curve_pts[j,:] - sat_curve_pts[i,:],2)))
+ if dist < 0.6*self._refined_cell_size:
+ j += 1
+ else:
+ i = j
+ j += 1
+ dists.append(dist)
+ sat_curve_upper_pts.append(factory.addPoint(sat_curve_pts[i,0] + dx*norm_vector[i, 0],\
+ sat_curve_pts[i,1] + dx*norm_vector[i, 1],0, 0.6*self._refined_cell_size))
+ sat_curve_lower_pts.append(factory.addPoint(sat_curve_pts[i,0] - dx*norm_vector[i, 0],\
+ sat_curve_pts[i,1] - dx*norm_vector[i, 1],0, 0.6*self._refined_cell_size))
+ sat_curve_connecting_lines = []
+ for i in range(len(sat_curve_upper_pts)):
+ sat_curve_connecting_lines.append(factory.addLine(sat_curve_lower_pts[i], sat_curve_upper_pts[i]))
+ sat_curve_upper_lines = []
+ sat_curve_lower_lines = []
+
+ sat_curve_crvloops = []
+ sat_curve_cornertags = []
+ for i in range(len(sat_curve_lower_pts)-1):
+ if dists[i] < self._base_cell_size:
+ c_upper=factory.addLine(sat_curve_upper_pts[i],sat_curve_upper_pts[i+1])
+ c_lower=factory.addLine(sat_curve_lower_pts[i],sat_curve_lower_pts[i+1])
+ sat_curve_upper_lines.append(c_upper)
+ sat_curve_lower_lines.append(c_lower)
+ sat_curve_crvloops.append(factory.addCurveLoop([c_upper, -sat_curve_connecting_lines[i+1], -c_lower, sat_curve_connecting_lines[i]]))
+ sat_curve_cornertags.append([sat_curve_lower_pts[i], sat_curve_lower_pts[i+1], sat_curve_upper_pts[i+1], sat_curve_upper_pts[i]])
+ factory.synchronize()
+ fluid_plane_crvloops = [curvloop] + [-c for c in sat_curve_crvloops]
+ fluid_surf = factory.addPlaneSurface(fluid_plane_crvloops)
+ sat_surfs = [factory.addPlaneSurface([c]) for c in sat_curve_crvloops]
+ factory.addPhysicalGroup(2, [fluid_surf] + sat_surfs)
+ else:
+ factory.synchronize()
+ fluid_surf = factory.addPlaneSurface([curvloop])
+ factory.addPhysicalGroup(2, [fluid_surf])
+ # Apply conditional refinement, where the refined cell size is applied in proximity to the refinement points
+ factory.synchronize()
+ threshold_fields = []
+ dist_field_ref_pt = gmsh.model.mesh.field.add("Distance")
+ gmsh.model.mesh.field.setNumbers(dist_field_ref_pt, "PointsList", ref_pt_ids)
+ gmsh.model.mesh.field.setNumber(dist_field_ref_pt, "Sampling", 100)
+ t_field_ref_pt = gmsh.model.mesh.field.add("Threshold")
+ gmsh.model.mesh.field.setNumber(t_field_ref_pt, "InField", dist_field_ref_pt)
+ gmsh.model.mesh.field.setNumber(t_field_ref_pt, "SizeMin", self._refined_cell_size)
+ gmsh.model.mesh.field.setNumber(t_field_ref_pt, "SizeMax", self._base_cell_size)
+ gmsh.model.mesh.field.setNumber(t_field_ref_pt, "DistMin", 0.5*self._refinement_radius)
+ gmsh.model.mesh.field.setNumber(t_field_ref_pt, "DistMax", 1.5*self._refinement_radius)
+ threshold_fields.append(t_field_ref_pt)
+
+ if add_sat_curve:
+ dist_field_sat_crv = gmsh.model.mesh.field.add("Distance")
+ gmsh.model.mesh.field.setNumbers(dist_field_sat_crv, \
+ "PointsList", sat_curve_upper_pts + sat_curve_lower_pts)
+ gmsh.model.mesh.field.setNumber(dist_field_sat_crv, "Sampling", 10)
+ t_field_sat_crv = gmsh.model.mesh.field.add("Threshold")
+ gmsh.model.mesh.field.setNumber(t_field_sat_crv, "InField", dist_field_sat_crv)
+ gmsh.model.mesh.field.setNumber(t_field_sat_crv, "SizeMin", self._refined_cell_size)
+ gmsh.model.mesh.field.setNumber(t_field_sat_crv, "SizeMax", self._base_cell_size)
+ gmsh.model.mesh.field.setNumber(t_field_sat_crv, "DistMin", 0.5*self._refinement_radius)
+ gmsh.model.mesh.field.setNumber(t_field_sat_crv, "DistMax", 1.5*self._refinement_radius)
+ threshold_fields.append(t_field_sat_crv)
+
+ back_field = gmsh.model.mesh.field.add("Min")
+ gmsh.model.mesh.field.setNumbers(back_field, "FieldsList", threshold_fields)
+ gmsh.model.mesh.field.setAsBackgroundMesh(back_field)
+
+ factory.synchronize()
+ if add_sat_curve:
+ for i in range(len(sat_curve_connecting_lines)):
+ gmsh.model.mesh.setTransfiniteCurve(sat_curve_connecting_lines[i], 2)
+ for i in range(len(sat_curve_lower_lines)):
+ gmsh.model.mesh.setTransfiniteCurve(sat_curve_lower_lines[i], 2)
+ gmsh.model.mesh.setTransfiniteCurve(sat_curve_upper_lines[i], 2)
+
+ for s, c in zip(sat_surfs, sat_curve_cornertags):
+ gmsh.model.mesh.setTransfiniteSurface(s, cornerTags=c)
+ gmsh.model.mesh.setRecombine(2, s)
+ # Generate 2D mesh and extract table nodes
+ gmsh.model.mesh.generate(2)
- CV_train_scaled = self._scaler.transform(CV_train)
- CV_test_scaled = self._scaler.transform(CV_test)
+ if show:
+ gmsh.fltk.run()
+ # Global nodes
+ nodeTags, coords, _ = gmsh.model.mesh.getNodes()
+ nodeTags = np.asarray(nodeTags, dtype=np.int64)
+ MeshPoints = np.asarray(coords, dtype=float).reshape(-1, 3)[:, :2]
- PPV_test = D_test[:, self._Fluid_Variables.index(var_to_test_for)]
- print("Done!")
+ order = np.argsort(nodeTags)
+ nodeTags_sorted = nodeTags[order]
- self._lookup_tree = Invdisttree(X=CV_train_scaled,z=D_train)
+ # 2) 2D elements
+ if fluid_surf is None:
+ elemTypes, _, elemNodeTags = gmsh.model.mesh.getElements(2)
+ else:
+ elemTypes, _, elemNodeTags = gmsh.model.mesh.getElements(2, fluid_surf)
+
+ tris = []
+ quads = []
+
+ for et, nodes_flat in zip(elemTypes, elemNodeTags):
+ if et == 2: # triangles with 3 nodes
+ tri_tags = np.asarray(nodes_flat, dtype=np.int64).reshape(-1, 3)
+ tris.append(self.__map_tags(tri_tags, nodeTags_sorted,order).reshape(-1, 3))
+ elif et == 3: # quad with 4 nodes
+ quad_tags = np.asarray(nodes_flat, dtype=np.int64).reshape(-1, 4)
+ quads.append(self.__map_tags(quad_tags, nodeTags_sorted,order).reshape(-1, 4))
+
+ tris = np.vstack(tris) if tris else np.zeros((0, 3), dtype=np.int64)
+
+ if quads:
+ quads = np.vstack(quads)
+ # split quad -> 2 tri: (0,1,2) + (0,2,3)
+ tris = np.vstack([
+ tris,
+ quads[:, [0, 1, 2]],
+ quads[:, [0, 2, 3]],
+ ])
+ gmsh.finalize()
- print("Search for best tree parameters...")
- # Do brute-force search to get the optimum number of nearest neighbors and distance power.
- n_near_range = range(1, 20)
- p_range = range(2, 6)
- RMS_ppv = np.zeros([len(n_near_range), len(p_range)])
- for i in tqdm(range(len(n_near_range))):
- for j in range(len(p_range)):
- PPV_predicted = self._lookup_tree(q=CV_test_scaled, nnear=n_near_range[i], p=p_range[j])[:, self._Fluid_Variables.index(var_to_test_for)]
- rms_local = mean_squared_error(y_true=PPV_test, y_pred=PPV_predicted)
- RMS_ppv[i,j] = rms_local
- [imin,jmin] = divmod(RMS_ppv.argmin(), RMS_ppv.shape[1])
- self._n_near = n_near_range[imin]
- self._p_fac = p_range[jmin]
- print("Done!")
- print("Best found number of nearest neighbors: "+str(self._n_near))
- print("Best found distance power: "+str(self._p_fac))
+ return MeshPoints, tris
+
+ def __map_tags(self,tags,nodeTags_sorted,order):
+ tags = np.asarray(tags, dtype=np.int64).ravel()
+ pos = np.searchsorted(nodeTags_sorted, tags)
+ ok = (pos < len(nodeTags_sorted)) & (nodeTags_sorted[pos] == tags)
+ if not np.all(ok):
+ missing = np.unique(tags[~ok])
+ raise RuntimeError(f"Node tags non trovati in getNodes(): {missing[:20]} (tot missing={len(missing)})")
+ return order[pos]
+
+ def __CalcMeshData(self, fluid_data_mesh:np.ndarray[float]):
+ """Calculate the fluid thermodynamic state variables for the table nodes
- def __EvaluateFluidInterpolator(self, CV_unscaled:np.ndarray):
- CV_scaled = self._scaler.transform(CV_unscaled)
- data_interp = self._lookup_tree(q=CV_scaled,nnear=self._n_near,p=self._p_fac)
- return data_interp
+ :param fluid_data_mesh: table mesh nodes of density and static energy
+ :type fluid_data_mesh: np.ndarray[float]
+ :return: filtered thermodynamic state data at the table nodes
+ :rtype: np.ndarray[float]
+ """
+ fluid_data_out = fluid_data_mesh.copy()
+ self.valid_mask = np.zeros(len(fluid_data_mesh),dtype=np.bool)
+ for i in tqdm(range(len(fluid_data_mesh)),desc="Evaluating fluid properties on table nodes..."):
+ try:
+ self._DataGenerator.UpdateFluid(fluid_data_mesh[i, EntropicVars.Density.value], fluid_data_mesh[i, EntropicVars.Energy.value])
+ state_vector, correct_phase = self._DataGenerator.GetStateVector()
+ if correct_phase:
+ fluid_data_out[i, :] = state_vector
+ self.valid_mask[i] = True
+ else:
+ fluid_data_out[i, :] = None
+ except:
+ fluid_data_out[i, :] = None
+ fluid_data_out = fluid_data_out[self.valid_mask,:]
+ return fluid_data_out
- def ComputeTableMesh(self):
- return
+ def __CartesianTableData(self):
+ print("Generating table on Cartesian grid")
+ Np_rho = self._Config.GetNpDensity()
+ Np_e = self._Config.GetNpEnergy()
+ self._DataGenerator.PreprocessData()
+ if self._Config.GetAutoRange():
+ rho_min, rho_max = self._DataGenerator.GetDensityBounds()
+ e_min, e_max = self._DataGenerator.GetEnergyBounds()
+ else:
+ rho_minmax = self._Config.GetDensityBounds()
+ rho_min = rho_minmax[0]
+ rho_max = rho_minmax[1]
+ e_minmax = self._Config.GetEnergyBounds()
+ e_min = e_minmax[0]
+ e_max = e_minmax[1]
+ rho_range = np.linspace(rho_min, rho_max, Np_rho)
+ e_range = np.linspace(e_min, e_max, Np_e)
+ self.rho_grid, self.e_grid = np.meshgrid(rho_range, e_range)
+
+ print(f"Grid Configuration:")
+ print(f" Density: [{rho_min:.2f}, {rho_max:.2f}] kg/m3 ({Np_rho} points)")
+ print(f" Energy: [{e_min:.0f}, {e_max:.0f}] J/kg ({Np_e} points)")
+ print(f" Total grid points: {Np_rho * Np_e:,}")
+ print()
+
+ shape = self.rho_grid.shape
+ n_points = shape[0] * shape[1]
+
+ # Initialize storage arrays
+ self.state_data = np.zeros([shape[0], shape[1], EntropicVars.N_STATE_VARS.value])
+
+ # Validity mask
+ self.valid_mask = np.zeros(shape, dtype=bool)
+
+ # Flatten for iteration
+ rho_flat = self.rho_grid.flatten()
+ e_flat = self.e_grid.flatten()
+
+ success_count = 0
+ for i in tqdm(range(n_points), desc="Evaluating"):
+ rho = rho_flat[i]
+ e = e_flat[i]
+ idx_2d = np.unravel_index(i, shape)
+ try:
+ self._DataGenerator.UpdateFluid(rho, e)
+ state_data, correct_phase = self._DataGenerator.GetStateVector()
+ if correct_phase:
+ self.state_data[idx_2d[0], idx_2d[1], :] = state_data
+ success_count += 1
+ self.valid_mask[idx_2d] = True
+ else:
+ self.state_data[idx_2d[0], idx_2d[1], :] = None
+ except:
+ self.state_data[idx_2d[0], idx_2d[1], :] = None
+ return
- def __ComputeCurvature(self):
- rho_min = self.__min_CV[self._controlling_variables.index("Density")]
- rho_max = self.__max_CV[self._controlling_variables.index("Density")]
- e_min = self.__min_CV[self._controlling_variables.index("Energy")]
- e_max = self.__max_CV[self._controlling_variables.index("Energy")]
-
- rho_range = (rho_min - rho_max)* (np.cos(np.linspace(0, 0.5*np.pi, 800))) + rho_max
- e_range = np.linspace(e_min, e_max, 100)
- xgrid, ygrid = np.meshgrid(rho_range, e_range)
-
- CV_probe = self.CV_full
- probe_data = self.__EvaluateFluidInterpolator(CV_probe)
-
- DT = Delaunay(self.CV_full)
- Tria = DT.simplices
- HullNodes = DT.convex_hull[:,0]
+ def __CartesianTriangulation(self):
+ """
+ Create Delaunay triangulation of valid grid points.
+ """
+ print("Creating Delaunay triangulation...")
+
+ # Extract valid points
+ rho_table = self.state_data[:,:,EntropicVars.Density.value]
+ e_table = self.state_data[:,:,EntropicVars.Energy.value]
+ rho_valid = rho_table[self.valid_mask].flatten()
+ e_valid = e_table[self.valid_mask].flatten()
+
+ # Stack as (N, 2) array
+ cv_table = np.column_stack([rho_valid, e_valid])
+
+ #self._table_nodes = np.column_stack(tuple(self.state_data[:,:,EntropicVars[v].value][self.valid_mask].flatten() for v in self._table_vars))
+ self._table_nodes = np.column_stack(tuple(self.state_data[:,:,i][self.valid_mask].flatten() for i in range(EntropicVars.N_STATE_VARS.value)))
+
+ # Create Delaunay triangulation
+ tri = Delaunay(cv_table)
+ self._table_connectivity = tri.simplices
+
+ # Identify hull nodes
+ edges = np.vstack([tri.simplices[:, [0, 1]],
+ tri.simplices[:, [1, 2]],
+ tri.simplices[:, [2, 0]]])
+ edges = np.sort(edges, axis=1)
+ unique_edges, counts = np.unique(edges, axis=0, return_counts=True)
+ boundary_edges = unique_edges[counts == 1]
+ self._table_hullnodes= np.unique(boundary_edges.flatten())
+
+ print(f" Triangulation nodes: {len(self._table_nodes):,}")
+ print(f" Triangles: {len(self._table_connectivity):,}")
+ print(f" Hull nodes: {len(self._table_hullnodes):,}")
+ print()
+ return
+
+ def __CreateSaturationCurve(self):
+ rhoe_sat_curve = self._DataGenerator.ComputeSaturationCurve()
+ rho_min, rho_max = self._DataGenerator.GetDensityBounds()
+ e_min, e_max = self._DataGenerator.GetEnergyBounds()
+ within_bounds_density = np.logical_and(rhoe_sat_curve[:,0] > rho_min, rhoe_sat_curve[:,0] < rho_max)
+ within_bounds_energy = np.logical_and(rhoe_sat_curve[:,1] > e_min, rhoe_sat_curve[:,1] < e_max)
+ within_bounds = np.logical_and(within_bounds_density, within_bounds_energy)
+
+ state_sat_curve = np.zeros([len(rhoe_sat_curve), EntropicVars.N_STATE_VARS.value])
+ state_sat_curve[:, EntropicVars.Density.value] = rhoe_sat_curve[:,0]
+ state_sat_curve[:, EntropicVars.Energy.value] = rhoe_sat_curve[:,1]
+
+ state_sat_curve_norm = self._fluid_data_scaler.transform(state_sat_curve[within_bounds, :])
+
+ sat_curve_pts_norm = state_sat_curve_norm[:, [EntropicVars.Density.value,EntropicVars.Energy.value]]
+ return sat_curve_pts_norm
+
+ def __GenerateMeshAndData(self, rhoe_norm:np.ndarray[float], sat_curve_pts:np.ndarray[float],ix_ref=[]):
- return CV_probe, probe_data, Tria, HullNodes
-
- # fig = plt.figure()
- # ax = plt.axes(projection='3d')
- # ax.plot3D(CV_probe[:,0],CV_probe[:,1],probe_data[:,self._Fluid_Variables.index("d2sdrho2")],'k.')
- # ax.plot3D(CV_probe[HullNodes, 0], CV_probe[HullNodes, 1], probe_data[HullNodes,self._Fluid_Variables.index("d2sdrho2")],'r.')
- # plt.show()
- # plt.triplot(CV_probe[:,0],CV_probe[:,1],DT.simplices)
- # plt.show()
+ rhoe_norm_mesh_nodes,tria = self.__Compute2DMesh(rhoe_norm, ref_pts=rhoe_norm[ix_ref,:], show=False, sat_curve_pts=sat_curve_pts)
+ # Calculate thermodynamic state variables of initial table nodes
+ fluid_data_norm = np.zeros([len(rhoe_norm_mesh_nodes), EntropicVars.N_STATE_VARS.value])
+ fluid_data_norm[:, EntropicVars.Density.value] = rhoe_norm_mesh_nodes[:,0]
+ fluid_data_norm[:, EntropicVars.Energy.value] = rhoe_norm_mesh_nodes[:,1]
+ fluid_data_mesh = self._fluid_data_scaler.inverse_transform(fluid_data_norm)
+ fluid_data_mesh = self.__CalcMeshData(fluid_data_mesh)
+ return fluid_data_mesh, tria
+
+ def GenerateTable(self):
+ """Initiate table generation process
+ """
- # rho_probe = probe_data[:, self._Fluid_Variables.index("Density")]
- # e_probe = probe_data[:, self._Fluid_Variables.index("Energy")]
+ self.__CartesianTableData()
+ self._table_nodes = np.column_stack(tuple(self.state_data[:,:,i][self.valid_mask].flatten() for i in range(EntropicVars.N_STATE_VARS.value)))
+ self._fluid_data_scaler.fit(self._table_nodes)
- # test_points = np.hstack((rho_probe[:,np.newaxis], e_probe[:,np.newaxis]))
- # CV_unique, idx_unique = np.unique(test_points,axis=0,return_index=True)
+ # Load initial fluid data and scale it
+ if self._Config.GetTableDiscretization()=="cartesian":
+ self.__CartesianTriangulation()
+ else:
+ print("Generating table with adaptive refinement")
+
+ self._table_nodes = np.column_stack(tuple(self.state_data[:,:,i][self.valid_mask].flatten() for i in range(EntropicVars.N_STATE_VARS.value)))
+
+ fluid_data_norm = self._fluid_data_scaler.fit_transform(self._table_nodes)
- # CV_vals_norm = self._scaler.transform(CV_probe)
+ sat_curve_pts_norm = self.__CreateSaturationCurve()
- # hull = ConvexHull(CV_vals_norm)
- # s_grid = probe_data[:, self._Fluid_Variables.index("s")]
- # s_grid = np.reshape(s_grid, np.shape(xgrid))
+ # Generate initial coarse table of fluid data
+ rhoe_norm = fluid_data_norm[:, [EntropicVars.Density.value, EntropicVars.Energy.value]]
+ print("Generating coarse thermodynamic mesh...")
+ fluid_data_coarse, _ = self.__GenerateMeshAndData(rhoe_norm, sat_curve_pts_norm)
+ print("Done!")
+ # Identify refinement locations
+ fluid_data_norm = self._fluid_data_scaler.transform(fluid_data_coarse)
+ ix_ref = self.__ApplyRefinement(fluid_data_norm)
- # idx_ref = self.__ComputeEntropyCurvature(s_grid)
- # x_refinement = CV_vals_norm[idx_ref, 0]
- # y_refinement = CV_vals_norm[idx_ref, 1]
- # x_hull = CV_vals_norm[hull.vertices, 0]
- # y_hull = CV_vals_norm[hull.vertices, 1]
+ # # Regenerate table including refinement locations
+ rhoe_norm_mesh = fluid_data_norm[:, [EntropicVars.Density.value, EntropicVars.Energy.value]]
- # XY_refinement = np.vstack((x_refinement, y_refinement)).T
- # XY_hull = np.vstack((x_hull, y_hull)).T
+ # Create triangulation of filtered thermodynamic state data
+ print("Generating refined thermodynamic mesh...")
+ fluid_data_ref, _ = self.__GenerateMeshAndData(rhoe_norm_mesh, sat_curve_pts_norm, ix_ref=ix_ref)
+ print("Done!")
+ fluid_data_norm_ref = self._fluid_data_scaler.transform(fluid_data_ref)
- # return XY_refinement, XY_hull, hull.area
-
- def __ComputeEntropyCurvature(self, s_interp:np.ndarray[float]):
- Q_norm = (s_interp - np.min(s_interp))/(np.max(s_interp) - np.min(s_interp))
- dQdy, dQdx = np.gradient(Q_norm)
- dQ_mag = np.sqrt(np.power(dQdy, 2) + np.power(dQdx, 2))
- dQ_norm = dQ_mag / np.max(dQ_mag)
- d2Qdy2, d2Qdx2 = np.gradient(dQ_norm)
- d2Q_mag = np.sqrt(np.power(d2Qdy2, 2) + np.power(d2Qdx2, 2))
- d2Q_norm = d2Q_mag / np.max(d2Q_mag)
- d2Q_norm = d2Q_norm.flatten()
- idx_ref = np.where(d2Q_norm > self._curvature_threshold/10)
- return idx_ref
+ DT = Delaunay(fluid_data_norm_ref[:, [EntropicVars.Density.value,EntropicVars.Energy.value]])
+
+ # Extract triangulation, hull nodes, and table data
+ Tria = DT.simplices
+ HullNodes = concave_hull_indexes(fluid_data_norm_ref[:, [EntropicVars.Density.value,EntropicVars.Energy.value]])
+
+ self._table_nodes = fluid_data_ref
+ self._table_connectivity = Tria
+ self._table_hullnodes = HullNodes
+
+ return
- def __Compute2DMesh(self, XY_hull:np.ndarray, XY_refinement:np.ndarray, level_area:float):
- """
- Generate a 2D mesh for the current table level.
-
- :param XY_hull: Array containing normalized pv and enth coordinates of the outline of the table level.
- :type XY_hull: NDArray
- :param XY_refinement: Array containing normalized pv and enth coordinates where refinement should be applied.
- :type XY_refinement: NDArray
- :return: mesh nodes of the 2D table mesh.
- :rtype: NDArray
- """
- gmsh.initialize()
+ def __remove_invalid_nodes_from_mesh(self, connectivity, valid_mask, rhoe_mesh_norm):
- gmsh.option.setNumber("General.Terminal", 0)
- gmsh.option.setNumber("General.Verbosity", 1)
- gmsh.model.add("table_level")
- factory = gmsh.model.geo
+ conn = np.asarray(connectivity, dtype=np.int64)
- base_cell_size = self._base_cell_size * level_area
- refined_cell_size = self._refined_cell_size * level_area
- refinement_radius = self._refinement_radius * np.sqrt(level_area)
+ tri_keep = np.all(valid_mask[conn], axis=1)
+ conn2 = conn[tri_keep]
- hull_pts = []
- for i in range(int(len(XY_hull)/2)):
- hull_pts.append(factory.addPoint(XY_hull[i, 0], XY_hull[i, 1], 0, base_cell_size))
- hull_pts_2 = [hull_pts[-1]]
- for i in range(int(len(XY_hull)/2), len(XY_hull)):
- hull_pts_2.append(factory.addPoint(XY_hull[i, 0], XY_hull[i, 1], 0, base_cell_size))
- hull_pts_2.append(hull_pts[0])
- embed_pts = []
- for i in range(len(XY_refinement)):
- pt_idx = factory.addPoint(XY_refinement[i, 0], XY_refinement[i, 1], 0, refined_cell_size)
- embed_pts.append(pt_idx)
-
-
- hull_curve_1 = factory.addPolyline(hull_pts)
- hull_curve_2 = factory.addPolyline(hull_pts_2)
-
- CL = factory.addCurveLoop([hull_curve_1, hull_curve_2])
-
- surf = factory.addPlaneSurface([CL])
- gmsh.model.addPhysicalGroup(1, [hull_curve_1], name="hull_curve_1")
- gmsh.model.addPhysicalGroup(1, [hull_curve_2], name="hull_curve_2")
- gmsh.model.addPhysicalGroup(2, [surf], name="table_level")
- gmsh.model.geo.synchronize()
-
- gmsh.model.mesh.field.add("Distance", 1)
- gmsh.model.mesh.field.setNumbers(1, "PointsList", embed_pts)
- gmsh.model.mesh.field.setNumber(1, "Sampling", 100)
- gmsh.model.mesh.field.add("Threshold", 2)
- gmsh.model.mesh.field.setNumber(2, "InField", 1)
- gmsh.model.mesh.field.setNumber(2, "SizeMin", refined_cell_size)
- gmsh.model.mesh.field.setNumber(2, "SizeMax", base_cell_size)
- gmsh.model.mesh.field.setNumber(2, "DistMin", refinement_radius)
- gmsh.model.mesh.field.setNumber(2, "DistMax", 1.5*refinement_radius)
-
- gmsh.model.mesh.field.add("Min", 7)
- gmsh.model.mesh.field.setNumbers(7, "FieldsList", [2])
- gmsh.model.mesh.field.setAsBackgroundMesh(7)
+ keep_nodes = np.flatnonzero(valid_mask)
+ old_to_new = -np.ones(len(rhoe_mesh_norm), dtype=np.int64)
+ old_to_new[keep_nodes] = np.arange(len(keep_nodes), dtype=np.int64)
+
+ conn2 = old_to_new[conn2]
+
+ return conn2
+ def WriteOutParaview(self,file_name_out:str="vtk_table"):
+ """
+ write a file containing all the LuT data that can be opened with Paraview
- lc = base_cell_size
- def meshSizeCallback(dim,tag,x,y,z,lc):
- return lc
+ :param file_name_out: string indicating the name and extension of the saved file
+ """
+
+ #x, y = self._table_nodes[:, EntropicVars.Density.value], self._table_nodes[:, EntropicVars.Energy.value]
+ table_data_norm = self._fluid_data_scaler.transform(self._table_nodes)
+ x, y = table_data_norm[:, EntropicVars.Density.value], table_data_norm[:, EntropicVars.Energy.value]
+ # scale_x= self._fluid_data_scaler.data_max_[EntropicVars.Density.value] - self._fluid_data_scaler.data_min_[EntropicVars.Density.value]
+ # scale_y= self._fluid_data_scaler.data_max_[EntropicVars.Energy.value] - self._fluid_data_scaler.data_min_[EntropicVars.Energy.value]
- gmsh.model.mesh.setSizeCallback(meshSizeCallback)
- gmsh.option.setNumber("Mesh.Algorithm", 5)
- gmsh.model.mesh.generate(2)
- nodes = gmsh.model.mesh.getNodes(dim=2, tag=-1, includeBoundary=True, returnParametricCoord=False)[1]
- MeshPoints = np.array([nodes[::3], nodes[1::3]]).T
- gmsh.finalize()
+ pts = np.column_stack([x, y, np.zeros_like(x)]) # z=0
- # Remove mesh nodes that are out of bounds.
- pv_norm, enth_norm = MeshPoints[:, 0], MeshPoints[:, 1]
+ conn = np.asarray(self._table_connectivity, dtype=np.int64)
+ if conn.min() == 1:
+ conn = conn - 1
- CV_level_norm = np.vstack((pv_norm, enth_norm)).T
- CV_level_dim = self._scaler.inverse_transform(CV_level_norm)
+ point_data = {}
+ for var in self._table_vars:
+ ivar = EntropicVars[var].value
+ point_data[var] = np.asarray(self._table_nodes[:, ivar])
- MeshPoints = np.zeros([np.shape(MeshPoints)[0], 2])
- MeshPoints[:, 0] = pv_norm
- MeshPoints[:, 1] = enth_norm
+ mesh = meshio.Mesh(
+ points=pts,
+ cells=[("triangle", conn)],
+ point_data=point_data
+ )
+ mesh.write("%s.vtk" % file_name_out)
- table_level_data = self.__EvaluateFluidInterpolator(CV_level_dim)
+ return
- return MeshPoints, table_level_data
+ def AddRefinementCriterion(self, TD_variable:str, norm_val_min:float=np.inf, norm_val_max:float=-np.inf):
+ """Apply refinement in the table where the normalized value of the thermodynamic variable lies between the specified bounds.
+
+ :param TD_variable: name of the thermodynamic variable for which to apply refinement
+ :type TD_variable: str
+ :param norm_val_min: lower bound of the normalized thermodynamic variable, defaults to np.inf
+ :type norm_val_min: float, optional
+ :param norm_val_max: upper bound of the normalized thermodynamic variable, defaults to -np.inf
+ :type norm_val_max: float, optional
+ :raises Exception: if thermodynamic state variable is unknown to SU2 DataMiner
+ """
+ if TD_variable not in self._table_vars:
+ raise Exception("%s is not present in fluid data" % TD_variable)
+
+ self.__refinement_vars.append(TD_variable)
+ self.__refinement_norm_min.append(norm_val_min)
+ self.__refinement_norm_max.append(norm_val_max)
+ return
+ def __ApplyRefinement(self, fluid_data_norm_ref:np.ndarray[float]):
+ ix_ref = np.array([],dtype=np.int64)
+ fluid_vars = [a.name for a in EntropicVars][:-1]
+ fluid_data_inv = self._fluid_data_scaler.inverse_transform(fluid_data_norm_ref)
+ for TD_var, val_min, val_max in zip(self.__refinement_vars, self.__refinement_norm_min, self.__refinement_norm_max):
+ norm_data_var = fluid_data_inv[:, fluid_vars.index(TD_var)]
+
+ ix = np.argwhere(np.logical_and(norm_data_var>=val_min, norm_data_var<=val_max))[:,0]
+ ix_ref = np.append(ix_ref, ix)
+ if len(ix_ref) > 0:
+ return np.unique(ix_ref)
+ else:
+ return []
+
+
def WriteTableFile(self, output_filepath:str=None):
"""
Save the table data and connectivity as a Dragon library file. If no file name is provided, the table file will be named according to the Config_FGM class name.
@@ -350,9 +772,9 @@ def WriteTableFile(self, output_filepath:str=None):
fid.write("%i\n" % np.shape(self._table_hullnodes)[0])
fid.write("\n")
- fid.write("[Number of variables]\n%i\n\n" % (len(self.table_vars)+2))
+ fid.write("[Number of variables]\n%i\n\n" % (len(self._table_vars)))
fid.write("[Variable names]\n")
- for iVar, Var in enumerate(self._controlling_variables + self.table_vars):
+ for iVar, Var in enumerate(self._table_vars):
fid.write(str(iVar + 1)+":"+Var+"\n")
fid.write("\n")
@@ -360,10 +782,10 @@ def WriteTableFile(self, output_filepath:str=None):
print("Writing table data...")
fid.write("\n")
- for iNode in range(np.shape(self._table_nodes)[0]):
- fid.write("\t".join("%+.14e" % cv for cv in self._table_nodes[iNode, :]))
- for var in self.table_vars:
- fid.write("\t%+.14e" % self.table_data[:, self._Fluid_Variables.index(var)][iNode])
+ for iNode in range(len(self._table_nodes)):
+ for var in self._table_vars:
+ ivar = EntropicVars[var].value
+ fid.write("\t%+.14e" % self._table_nodes[iNode, ivar])
fid.write("\n")
fid.write("\n\n")
print("Done!")
@@ -383,3 +805,6 @@ def WriteTableFile(self, output_filepath:str=None):
print("Done!")
fid.close()
+
+ return
+
\ No newline at end of file
diff --git a/README.md b/README.md
index 7bdaa41..77a7590 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@
# SU2 DataMiner
-This repository describes the workflow for manifold generation for data-driven fluid modeling in SU2. The workflow allows the user to generate fluid data and convert these into tables and train multi-layer perceptrons in order to retrieve thermo-chemical quantities during simulations in SU2. The applications are currently limited to non-ideal computational fluid dynamics and flamelet-generated manifold simulations for arbitrary fluids and reactants respectively. Documentation of the source code, the theory guide, and tutorial collection can be found on the [SU2 DataMiner web page](https://propulsion-power-tu-delft.github.io/SU2_DataMiner/).
+This repository describes the workflow for manifold generation for data-driven fluid modeling in SU2. The workflow allows the user to generate fluid data and convert these into tables and train multi-layer perceptrons in order to retrieve thermo-chemical quantities during simulations in SU2. The applications are currently limited to non-ideal computational fluid dynamics and flamelet-generated manifold simulations for arbitrary fluids and reactants respectively. Documentation of the source code, the theory guide, and tutorial collection can be found on the [SU2 DataMiner web page](https://su2code.github.io/SU2_DataMiner/).
## Capabilities
The SU2 DataMiner workflow allows the user to generate fluid data and convert these into look-up tables (LUT) or multi-layer perceptrons (MLP) for usage in SU2 simulations. The types of simulations for which this workflow is suitable are flamelet-generated manifold (FGM) and non-ideal computational fluid dynamics (NICFD) simulations. This tool allows the user to start from scratch and end up with a table input file or a set of MLP input files which can immediately be used within SU2.
@@ -33,17 +33,4 @@ export PATH=$PATH:$PINNTRAINING_HOME/bin
where `````` is the path to where you cloned the repository.
-Tutorials can be found under ```TestCases```, proper documentation will follow soon.
-
-## Getting Started
-
-Generating a fluid data manifold for FGM or NICFD applications consists of the following steps:
-
-1. *Generate a configuration*: The manifold settings such as the type of fluid, storage directory and range are stored in a configuration class. Configurations can be defined through python (example scripts are found in the TestCases folder, named ```generate_config.py```) or interactively through terminal inputs. In order to generate a configuration interactively, run the command ```GenerateConfig.py``` in the terminal.
-
-2. *Generate fluid data*: Raw fluid data can be generated once the configuration is defined. Similarly to the configuration set-up, fluid data can be generated through a python interface enabling more flexibility, or through the terminal. Run the command ```GenerateFlameletData.py -h``` to see the available options. Optionally, flamelet data can be visualized through the ```PlotFlamelets.py``` command.
-
-3. *Process fluid data*: Raw fluid data needs to be processed in order to be converted into a manifold usable in SU2. Especially a flamelet-based manifold requires additional steps to convert raw flamelet data into a usable manifold. These steps include optimizing the progress variable, homogenizing the flamelet data, and grouping the various flamelet data into groups of high correlation. TODO: write executable for this step.
-
-4. *Generate manifold*: The processed fluid data is ready for conversion into a manifold for SU2 simulations. The available formats are the look-up table (LUT) and multi-layer perceptron (MLP).
diff --git a/TestCases/NICFD/LUT_Twophase/generate_table.py b/TestCases/NICFD/LUT_Twophase/generate_table.py
new file mode 100644
index 0000000..7c878a4
--- /dev/null
+++ b/TestCases/NICFD/LUT_Twophase/generate_table.py
@@ -0,0 +1,46 @@
+from su2dataminer.config import Config_NICFD
+from su2dataminer.manifold import SU2TableGenerator_NICFD
+
+# Generate properties of MM with REFPROP library.
+config = Config_NICFD()
+config.SetFluid("MM")
+config.SetEquationOfState("REFPROP")
+
+# Include gas, liquid, two-phase, and supercritical fluid data.
+config.EnableTwophase(True)
+config.EnableLiquidPhase(True)
+config.EnableGasPhase(True)
+config.EnableSuperCritical(True)
+
+# Include visosity, conductivity, and vapor quality.
+config.IncludeTransportProperties(True)
+config.UseAutoRange(False)
+config.SetDensityBounds(0.1, 460)
+config.SetEnergyBounds(200e3, 360e3)
+config.SetNpDensity(50)
+config.SetNpEnergy(50)
+
+# Initiate table generator with adaptive triangulation.
+tablegen = SU2TableGenerator_NICFD(config)
+tablegen.SetTableDiscretization("adaptive")
+
+# Specify table resolution for coarse and refined sections.
+tablegen.SetCellSize_Coarse(2e-2)
+tablegen.SetCellSize_Refined(2e-3)
+
+# Relative step size for finite-differences.
+tablegen.SetFDStepSize(7e-3)
+
+# Optionally, specify thermophysical variables to be included in the table. By default, all variables are included.
+# tablegen.SetTableVars(["Density","Energy","s","p","T", "dsdrho_e","dsde_rho", "d2sdrho2","d2sde2","d2sdedrho","VaporQuality","ViscosityDyn","Conductivity"])
+
+# Specify custom refinement regions (low density, around Trova isentrope)
+tablegen.AddRefinementCriterion("Density", 0.0, 10.0)
+tablegen.AddRefinementCriterion("s", 729.13-2, 729.13+40)
+
+# Generate table.
+tablegen.GenerateTable()
+
+# Write SU2 DRG and vtk table.
+tablegen.WriteTableFile("LUT_adaptive.drg")
+tablegen.WriteOutParaview("vtktable_adaptive")
\ No newline at end of file
diff --git a/docker/build/Dockerfile b/docker/build/Dockerfile
index d66e747..673e2f4 100644
--- a/docker/build/Dockerfile
+++ b/docker/build/Dockerfile
@@ -10,6 +10,7 @@ RUN apt-get update && apt-get install -y \
build-essential \
swig \
ccache \
+ gmsh \
&& rm -rf /var/lib/apt/lists/* \
&& update-alternatives --install /usr/bin/python python /usr/bin/python3 10 \
&& /usr/sbin/update-ccache-symlinks \
diff --git a/environment.yml b/environment.yml
index 3af31e9..48ba895 100644
--- a/environment.yml
+++ b/environment.yml
@@ -25,3 +25,6 @@ dependencies:
- pygad
- paretoset
- pymoo
+- pydata-sphinx-theme
+- meshio
+- concave_hull
\ No newline at end of file
diff --git a/required_packages.txt b/required_packages.txt
index d2a13bf..10f3d7f 100644
--- a/required_packages.txt
+++ b/required_packages.txt
@@ -17,4 +17,6 @@ tqdm==4.67.1
pygad==3.5.0
paretoset==1.2.5
pymoo==0.6.1.6
-pydata-sphinx-theme
\ No newline at end of file
+pydata-sphinx-theme
+meshio==5.3.5
+concave_hull==0.1.2
\ No newline at end of file
diff --git a/su2dataminer/manifold.py b/su2dataminer/manifold.py
index d97cd20..14605d8 100644
--- a/su2dataminer/manifold.py
+++ b/su2dataminer/manifold.py
@@ -1,3 +1,4 @@
from Manifold_Generation.MLP.Trainers_FGM.Trainers import *
from Manifold_Generation.MLP.Trainers_NICFD.Trainers import *
-from Manifold_Generation.MLP.optimizeHP import *
\ No newline at end of file
+from Manifold_Generation.MLP.optimizeHP import *
+from Manifold_Generation.LUT.LUTGenerators import *
\ No newline at end of file