diff --git a/Common/DataDrivenConfig.py b/Common/DataDrivenConfig.py index ad1c69f..3b63602 100644 --- a/Common/DataDrivenConfig.py +++ b/Common/DataDrivenConfig.py @@ -786,10 +786,13 @@ class Config_FGM(Config): __generate_extra_interpolated_burnerflames:bool = True # Generate extra interpolated burner-stabilized flamelets __generate_equilibrium:bool = DefaultSettings_FGM.include_equilibrium # Generate chemical equilibrium data __generate_counterflames:bool = DefaultSettings_FGM.include_counterflames # Generate counter-flow diffusion flamelets. + __counterflow_fixed_strain:bool = DefaultSettings_FGM.counterflow_fixed_strain # Fixed-strain mode for counterflow flames. + __counterflow_strain_rate:float = DefaultSettings_FGM.counterflow_strain_rate # Global strain rate [1/s] for fixed-strain mode. __flamelet_types:list[str] = [FlameletSolverOptions[0]] __write_MATLAB_files:bool = False # Write TableGenerator compatible flamelet files. + __save_mole_fractions:bool = False # Save mole fractions (X-) in flamelet CSV files. gas:ct.Solution = None # Cantera solution object. __species_in_mixture:list[str] = None # Species names in mixture. @@ -859,6 +862,9 @@ def __init__(self, load_file:str=None): self.__dict__ = loaded_config.__dict__.copy() print("Loaded configuration file with name " + loaded_config.GetConfigName()) else: + # Initialize __flamelet_types as an instance variable (not shared class variable) + self._Config_FGM__flamelet_types = [FlameletSolverOptions[0]] + self.SetAlphaExpo(DefaultSettings_FGM.init_learning_rate_expo) self.SetLRDecay(DefaultSettings_FGM.learning_rate_decay) self.SetBatchExpo(DefaultSettings_FGM.batch_size_exponent) @@ -964,7 +970,10 @@ def PrintBanner(self): if self.__generate_equilibrium: print("-Chemical equilibrium data") if self.__generate_counterflames: - print("-Counter-flow diffusion flamelet data") + if self.__counterflow_fixed_strain: + print("-Counter-flow diffusion flamelet data (fixed strain rate: %.1f 1/s)" % self.__counterflow_strain_rate) + else: + print("-Counter-flow diffusion flamelet data (ramp to extinction)") print("") print("Flamelet manifold data characteristics: ") @@ -1495,17 +1504,17 @@ def includeFlameletType(self, flamelet_type:str="FREEFLAME"): if flamelet_type not in FlameletSolverOptions: raise Exception("%s is not recognized as a viable flamelet type" % flamelet_type) - + if flamelet_type not in self.__flamelet_types: self.__flamelet_types.append(flamelet_type) - + return def setFlameletTypes(self, flamelet_types:list[str]): self.__flamelet_types = flamelet_types.copy() self.__checkFlameletTypes() return - + def __checkFlameletTypes(self): self.__generate_burnerflames = "BURNERFLAME" in self.__flamelet_types self.__generate_freeflames = "FREEFLAME" in self.__flamelet_types @@ -1513,18 +1522,18 @@ def __checkFlameletTypes(self): self.__generate_equilibrium = "EQUILIBRIUM" in self.__flamelet_types self.__generate_extra_interpolated_burnerflames = "BURNERFLAME_INT" in self.__flamelet_types return - + def excludeFlameletType(self, flamelet_type:str): if flamelet_type not in FlameletSolverOptions: raise Exception("%s is not recognized as a viable flamelet type" % flamelet_type) - + if flamelet_type in self.__flamelet_types: self.__flamelet_types.remove(flamelet_type) self.__checkFlameletTypes() - if len(self.__flamelet_types) == 0: - raise Exception("At least one flamelet type should be included in the manifold") - + # Allow empty list temporarily during configuration setup + # Validation will happen when the config is actually used + return def getFlameletTypes(self): @@ -1539,7 +1548,10 @@ def RunFreeFlames(self, input:bool=DefaultSettings_FGM.include_freeflames): """ self.__generate_freeflames = input - self.includeFlameletType("FREEFLAME") + if input: + self.includeFlameletType("FREEFLAME") + else: + self.excludeFlameletType("FREEFLAME") return def RunBurnerFlames(self, input:bool=DefaultSettings_FGM.include_burnerflames): """ @@ -1550,7 +1562,10 @@ def RunBurnerFlames(self, input:bool=DefaultSettings_FGM.include_burnerflames): """ self.__generate_burnerflames = input - self.includeFlameletType("BURNERFLAME") + if input: + self.includeFlameletType("BURNERFLAME") + else: + self.excludeFlameletType("BURNERFLAME") return def RunEquilibrium(self, input:bool=DefaultSettings_FGM.include_equilibrium): @@ -1562,7 +1577,10 @@ def RunEquilibrium(self, input:bool=DefaultSettings_FGM.include_equilibrium): """ self.__generate_equilibrium = input - self.includeFlameletType("EQUILIBRIUM") + if input: + self.includeFlameletType("EQUILIBRIUM") + else: + self.excludeFlameletType("EQUILIBRIUM") return def RunCounterFlames(self, input:bool=DefaultSettings_FGM.include_counterflames): @@ -1574,7 +1592,10 @@ def RunCounterFlames(self, input:bool=DefaultSettings_FGM.include_counterflames) """ self.__generate_counterflames = input - self.includeFlameletType("COUNTERFLAME") + if input: + self.includeFlameletType("COUNTERFLAME") + else: + self.excludeFlameletType("COUNTERFLAME") return def RunExtraInterpolatedBurnerFlames(self, input:bool=True): @@ -1586,7 +1607,10 @@ def RunExtraInterpolatedBurnerFlames(self, input:bool=True): """ self.__generate_extra_interpolated_burnerflames = input - self.includeFlameletType("INT_BURNERFLAME") + if input: + self.includeFlameletType("INT_BURNERFLAME") + else: + self.excludeFlameletType("INT_BURNERFLAME") return def GenerateFreeFlames(self): @@ -1625,6 +1649,70 @@ def GenerateCounterFlames(self): """ return self.__generate_counterflames + def SetCounterFlowFixedStrain(self, fixed:bool=True): + """ + Select fixed-strain-rate mode for counter-flow diffusion flames. + + When True, one flame is solved per temperature level at the strain rate + set by SetCounterFlowStrainRate. When False (default), the existing + ramp-to-extinction behaviour is used. + + :param fixed: enable fixed-strain mode. + :type fixed: bool + """ + self.__counterflow_fixed_strain = fixed + return + + def GetCounterFlowFixedStrain(self) -> bool: + """ + Whether fixed-strain mode is enabled for counter-flow flames. + + :return: fixed-strain mode is active. + :rtype: bool + """ + return self.__counterflow_fixed_strain + + def SetCounterFlowStrainRate(self, strain_rate:float): + """ + Set the global strain rate used in fixed-strain counter-flow flame mode. + + :param strain_rate: global strain rate in 1/s. + :type strain_rate: float + :raises Exception: if strain_rate is not strictly positive. + """ + if strain_rate <= 0: + raise Exception("Counter-flow strain rate must be strictly positive.") + self.__counterflow_strain_rate = strain_rate + return + + def GetCounterFlowStrainRate(self) -> float: + """ + Return the global strain rate used in fixed-strain counter-flow flame mode. + + :return: global strain rate [1/s]. + :rtype: float + """ + return self.__counterflow_strain_rate + + def SetSaveMoleFractions(self, save_mole_fractions:bool): + """ + Enable or disable saving mole fractions (X-species) in flamelet CSV files. + + :param save_mole_fractions: whether to save mole fractions alongside mass fractions. + :type save_mole_fractions: bool + """ + self.__save_mole_fractions = save_mole_fractions + return + + def GetSaveMoleFractions(self) -> bool: + """ + Check whether mole fractions (X-species) should be saved in flamelet CSV files. + + :return: True if mole fractions should be saved. + :rtype: bool + """ + return self.__save_mole_fractions + def GenerateExtraInterpolatedBurnerFlames(self): """ Whether the manifold data contains extra interpolated burner-stabilized flame data. diff --git a/Common/Properties.py b/Common/Properties.py index 284599b..51abfd8 100644 --- a/Common/Properties.py +++ b/Common/Properties.py @@ -218,6 +218,9 @@ class DefaultSettings_FGM(DefaultProperties): include_equilibrium:bool = True include_counterflames:bool = False + counterflow_fixed_strain:bool = False # False = ramp-to-extinction; True = fixed strain rate + counterflow_strain_rate:float = 56.0 # global strain rate [1/s] for fixed-strain mode + affinity_threshold:float = 0.7 output_file_header:str = "flamelet_data" boundary_file_header:str = "boundary_data" diff --git a/Data_Generation/DataGenerator_FGM.py b/Data_Generation/DataGenerator_FGM.py index 57f40f8..0103689 100644 --- a/Data_Generation/DataGenerator_FGM.py +++ b/Data_Generation/DataGenerator_FGM.py @@ -76,7 +76,7 @@ class DataGenerator_Cantera(DataGenerator_Base): __define_equivalence_ratio:bool = not DefaultSettings_FGM.run_mixture_fraction # Define unburnt mixture via the equivalence ratio __unb_mixture_status:list[float] = [] __flameletSolverDict:Dict[str, FlameletSolver_Cantera] - + __run_freeflames:bool = DefaultSettings_FGM.include_freeflames # Run adiabatic flame computations __run_extra_interpolated_burnerflames:bool = True # Run extra interpolated burner-stabilized flame computations __run_burnerflames:bool = DefaultSettings_FGM.include_burnerflames # Run burner stabilized flame computations @@ -138,6 +138,10 @@ def __SynchronizeSettings(self): def SetLoglevel(self, loglevel:int=0): """Set Cantera solver verbosity level (0=silent).""" self.__loglevel = loglevel + # Propagate loglevel to all existing flamelet solvers + for solver in self.__flameletSolverDict.values(): + solver.setCanteraVerbose(loglevel) + solver.setSolverVerbose(loglevel) def SetInitialGridLength(self, length:float): """Set the initial domain length (in metres) for new flamelet grids.""" @@ -145,10 +149,11 @@ def SetInitialGridLength(self, length:float): def setRefinementCriteria(self, flamelet_type:str, ratio:float=3, slope:float=0.03, curve:float=0.03, prune:float=0.01): if flamelet_type not in self.__flameletSolverDict.keys(): - raise Exception("%s is included in the available flamelet types" % flamelet_type) + available = list(self.__flameletSolverDict.keys()) + raise Exception("%s is NOT included in the available flamelet types. Available: %s" % (flamelet_type, available)) self.__flameletSolverDict[flamelet_type].setGridRefinementCriteria(ratio, slope, curve, prune) return - + def SetMdotDHTarget(self, dH_target:float): """Set the target enthalpy step between successive burner-stabilized flamelets. @@ -316,6 +321,48 @@ def RunExtraInterpolatedBurnerFlames(self, input:bool=True): self.__SynchronizeSettings() return + def SetFreeFlameRefineCriteria(self, **kwargs): + """Set refinement criteria for free-flame solver using the existing setRefinementCriteria method. + + This is a convenience wrapper that calls setRefinementCriteria with flamelet_type="FREEFLAME". + + :param ratio: refinement ratio + :type ratio: int, optional + :param slope: refinement slope criterion + :type slope: float, optional + :param curve: refinement curvature criterion + :type curve: float, optional + :param prune: refinement prune criterion + :type prune: float, optional + """ + ratio = kwargs.get('ratio', 2) + slope = kwargs.get('slope', 0.025) + curve = kwargs.get('curve', 0.025) + prune = kwargs.get('prune', 0.01) + self.setRefinementCriteria("FREEFLAME", ratio, slope, curve, prune) + return + + def SetBurnerFlameRefineCriteria(self, **kwargs): + """Set refinement criteria for burner-flame solver using the existing setRefinementCriteria method. + + This is a convenience wrapper that calls setRefinementCriteria with flamelet_type="BURNERFLAME". + + :param ratio: refinement ratio + :type ratio: int, optional + :param slope: refinement slope criterion + :type slope: float, optional + :param curve: refinement curvature criterion + :type curve: float, optional + :param prune: refinement prune criterion + :type prune: float, optional + """ + ratio = kwargs.get('ratio', 2) + slope = kwargs.get('slope', 0.025) + curve = kwargs.get('curve', 0.025) + prune = kwargs.get('prune', 0.01) + self.setRefinementCriteria("BURNERFLAME", ratio, slope, curve, prune) + return + def SetNpMdotExtra(self, n_extra:int): """Set the number of interpolation steps for extra interpolated burner-stabilized flamelets. @@ -339,7 +386,7 @@ def SetMixtureValues(self, mixture_values:list[float]): raise Exception("At least one mixture status value should be provided.") if any([phi < 0 for phi in mixture_values]): raise Exception("Mixture values should be strictly positive") - + self.__unb_mixture_status = [phi for phi in mixture_values] return @@ -390,7 +437,7 @@ def computePremixedFlameletsFor(self, mix_status:float): if mix_status < 0: raise Exception("Mixture status value should be positive.") - + correct_order = FlameletSolverDict.keys() for c in correct_order: if c in self.__flameletSolverDict.keys() and c != "COUNTERFLAME": @@ -399,12 +446,12 @@ def computePremixedFlameletsFor(self, mix_status:float): flameletSolver.solveForMixtureStatus(mix_status) return - + def computeNonPremixedFlamelets(self): if "COUNTERFLAME" in self.__flameletSolverDict.keys(): self.__flameletSolverDict["COUNTERFLAME"].solveForMixtureStatus(0) return - + def ComputeFlamelets(self): """Generate and store all flamelet data for the current settings. """ @@ -414,9 +461,9 @@ def ComputeFlamelets(self): for mix_status in self.__unb_mixture_status: self.computePremixedFlameletsFor(mix_status) return - + def ComputeFlameletData(Config:Config_FGM, run_parallel:bool=False, N_processors:int=2, loglevel:int=0, - free_flame_refine:dict=None, burner_flame_refine:dict=None): + free_flame_refine:dict=None, burner_flame_refine:dict=None, counter_flame_refine:dict=None): """Generate flamelet data according to Config_FGM settings either in serial or parallel. :param Config: Config_FGM class containing manifold and flamelet generation settings. @@ -433,6 +480,9 @@ def ComputeFlameletData(Config:Config_FGM, run_parallel:bool=False, N_processors :param burner_flame_refine: Cantera burner-flame refinement criteria dict with keys ratio, slope, curve, prune. If None, the DataGenerator_Cantera defaults are used. :type burner_flame_refine: dict, optional + :param counter_flame_refine: Cantera counter-flow flame refinement criteria dict with keys ratio, slope, curve, prune. + If None, the DataGenerator_Cantera defaults are used. + :type counter_flame_refine: dict, optional :raises Exception: If number of processors is set to zero when running in parallel. """ @@ -460,6 +510,8 @@ def _make_generator(): F.SetFreeFlameRefineCriteria(**free_flame_refine) if burner_flame_refine is not None: F.SetBurnerFlameRefineCriteria(**burner_flame_refine) + if counter_flame_refine is not None: + F.setRefinementCriteria("COUNTERFLAME", **counter_flame_refine) return F # Set up Cantera flamelet generator object diff --git a/Data_Generation/FlameletSolvers.py b/Data_Generation/FlameletSolvers.py index 2ecd015..362718c 100644 --- a/Data_Generation/FlameletSolvers.py +++ b/Data_Generation/FlameletSolvers.py @@ -55,7 +55,7 @@ class FlameletSolver_Cantera: _from_file:bool=False _flameletSolutionForRestart:ct.FlameBase = None - + _T_reactants:float = DefaultSettings_FGM.T_min _reactant_mixture_status:float = 0.0 _pressure:float = ct.one_atm @@ -85,17 +85,17 @@ def __init__(self, config_input:Config_FGM): self._Config = config_input self._canteraSolution = ct.Solution(self._Config.GetReactionMechanism()) return - + def _initializeFlameletSolver(self): return - + def solveFor(self, **flamelet_solution_settings): """Calculate flamelet solution for custom settings. """ self._parseInputSettings(flamelet_solution_settings) self.startSolver() return - + def solveAndSaveFor(self, **flamelet_solution_settings): """Calculate flamelet and save solution for custom settings. """ @@ -103,7 +103,7 @@ def solveAndSaveFor(self, **flamelet_solution_settings): self.startSolver() self.saveFlameletSolution() return - + def solveForMixtureStatus(self, val_mixture_status:float, save:bool=True): """Calculate flamelet solutions for a specific mixture. @@ -126,17 +126,17 @@ def solveForMixtureStatus(self, val_mixture_status:float, save:bool=True): self.resetRestart() return - + def saveFlameletSolution(self): """Store flamelet solution in appropriately named folder. """ self._prepareStorageFolder() self._writeOutput() return - + def _parseInputSettings(self, flamelet_solution_settings): return - + def _prepareStorageFolder(self): folder_for_flamelet_type = self._createFolderForFlameletType() @@ -145,16 +145,16 @@ def _prepareStorageFolder(self): filepath_for_flamelet_data = sep.join((folder_for_flamelet_type, mixture_subfolder)) if not path.isdir(filepath_for_flamelet_data): mkdir(filepath_for_flamelet_data) - + self._output_filepath = filepath_for_flamelet_data return - + def _createFolderForFlameletType(self): storage_folder = sep.join((self._Config.GetOutputDir(), self.getFlameletFolder())) if not path.isdir(storage_folder): mkdir(storage_folder) return storage_folder - + def __createSubFolderForMixture(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="mixfrac" @@ -163,7 +163,7 @@ def __createSubFolderForMixture(self): mixture_subfolder = "%s_%.4f" % (tag_for_mixture_status, self._reactant_mixture_status) return mixture_subfolder - + def startSolver(self): """Initiate flamelet simulation and report status to terminal. """ @@ -180,7 +180,7 @@ def _prepareFlameletSimulation(self): self._solverSpecificPreprocessing() self._commonPreprocessing() return - + def _prepareFlameletSolver(self): if not self._from_restart and not self._from_file: self._initializeFlameletSolver() @@ -188,14 +188,14 @@ def _prepareFlameletSolver(self): else: self._flameletSolution = self._flameletSolutionForRestart return - + def __printToTerminal(self): return self.__flameletSolverLogLevel > 0 - + def _printStatusToTerminal(self): return - + def setReactantTemperature(self, Temp_reactants:float=DefaultSettings_FGM.T_min): """Specify the temperature of the reactants at the inflow boundary. @@ -208,7 +208,7 @@ def setReactantTemperature(self, Temp_reactants:float=DefaultSettings_FGM.T_min) raise Exception("Reactant temperature should be strictly positive.") self._T_reactants = Temp_reactants return - + def setPressure(self, val_pressure:float=DefaultSettings_FGM.pressure): """Specify the pressure at which flamelet solutions are calculated. @@ -220,7 +220,7 @@ def setPressure(self, val_pressure:float=DefaultSettings_FGM.pressure): raise Exception("Pressure should be strictly positive") self._pressure = val_pressure return - + def setMixtureStatus(self, val_mixture_status:float): """Specify the equivalence ratio or mixture fraction for premixed flamelets. @@ -234,25 +234,25 @@ def setMixtureStatus(self, val_mixture_status:float): raise Exception("Mixture fraction should be between zero and one.") if val_mixture_status < 0.0: raise Exception("Mixture status should be strictly positive.") - + self._reactant_mixture_status = val_mixture_status return - + def getReactantTemperature(self): return self._T_reactants - + def getMixtureStatus(self): return self._reactant_mixture_status - + def getFlameletFolder(self): return self._flameletTypeOutputFolder - + def getFlameletType(self): return self._flamelet_type - + def getPlotLabel(self): return self._plotLabel - + def setInitialGrid(self, initial_grid_length_in_meters:float=1.8e-2, number_of_nodes:int=100): """Specify the settings for the initial grid used to calculate the flamelet solution. @@ -265,24 +265,81 @@ def setInitialGrid(self, initial_grid_length_in_meters:float=1.8e-2, number_of_n """ if initial_grid_length_in_meters < 0 or initial_grid_length_in_meters >= self._max_grid_length: raise Exception("Initial grid length should be between 0 and %.1f." % self._max_grid_length) - + if number_of_nodes < 10: raise Exception("The initial grid should contain at least 10 nodes.") self._initial_grid_length = initial_grid_length_in_meters self.__initial_grid_number_of_points = number_of_nodes return - + def setGridRefinementCriteria(self, ratio:int=2, slope:float=0.025, curve:float=0.025, prune=0.01): + """Set grid refinement criteria using positional parameters. + + :param ratio: refinement ratio (default 2) + :type ratio: int, optional + :param slope: refinement slope criterion (default 0.025) + :type slope: float, optional + :param curve: refinement curvature criterion (default 0.025) + :type curve: float, optional + :param prune: refinement prune criterion (default 0.01) + :type prune: float, optional + """ self.__grid_refinement_ratio = ratio self.__grid_refinement_curve = curve self.__grid_refinement_prune = prune self.__grid_refinement_slope = slope return - + + def setRefinementCriteria(self, flamelet_type:str="COUNTERFLAME", **kwargs): + """Set refinement criteria for a specific flamelet type using keyword arguments. + + Supports Cantera refinement parameters: ratio, slope, curve, prune. + + :param flamelet_type: Type of flamelet ('COUNTERFLAME', 'FREEFLAME', 'BURNERFLAME'), defaults to 'COUNTERFLAME' + :type flamelet_type: str, optional + :param ratio: refinement ratio (default 3 for counterflame) + :type ratio: int, optional + :param slope: refinement slope criterion (default 0.04 for counterflame) + :type slope: float, optional + :param curve: refinement curvature criterion (default 0.06 for counterflame) + :type curve: float, optional + :param prune: refinement prune criterion (default 0.02 for counterflame) + :type prune: float, optional + :raises ValueError: if unsupported flamelet type is provided + """ + supported_types = ["COUNTERFLAME", "FREEFLAME", "BURNERFLAME"] + if flamelet_type not in supported_types: + raise ValueError(f"Unsupported flamelet type '{flamelet_type}'. Supported types: {supported_types}") + + # Extract refinement parameters with defaults based on flamelet type + if flamelet_type == "COUNTERFLAME": + default_ratio = 3 + default_slope = 0.04 + default_curve = 0.06 + default_prune = 0.02 + else: # FREEFLAME, BURNERFLAME + default_ratio = 2 + default_slope = 0.025 + default_curve = 0.025 + default_prune = 0.01 + + ratio = kwargs.get('ratio', default_ratio) + slope = kwargs.get('slope', default_slope) + curve = kwargs.get('curve', default_curve) + prune = kwargs.get('prune', default_prune) + + self.setGridRefinementCriteria(ratio=ratio, slope=slope, curve=curve, prune=prune) + return + def getGridRefinementCriteria(self): + """Get current grid refinement criteria. + + :return: tuple of (ratio, slope, curve, prune) + :rtype: tuple + """ return self.__grid_refinement_ratio, self.__grid_refinement_slope, self.__grid_refinement_curve, self.__grid_refinement_prune - + def setCanteraVerbose(self, verbose_level:int=0): """Specify verbosity of Cantera solution process. @@ -291,7 +348,7 @@ def setCanteraVerbose(self, verbose_level:int=0): """ self.__cantera_loglevel = verbose_level return - + def setSolverVerbose(self, verbose_level:int=1): """Specify the verbosity of the FlameletSolver solution process. @@ -300,35 +357,35 @@ def setSolverVerbose(self, verbose_level:int=1): """ self.__flameletSolverLogLevel = verbose_level return - + def _prepareSettingRange(self): self._print_iteration=True return - + def setInputVariable(self, val_input:float): return - + def getFlameletFileName(self, val_input:float): return "" - + def retrieveSolverSettings(self, solvers): return - - + + def _writeSolverSettings(self, solution_index:int, setting_1D:float): self._iteration = solution_index return - - - + + + def _writeOutput(self): if self.isConverged() and self.isBurning(): self.saveSolution() self._flameletSolutionForRestart = self._flameletSolution return - - - + + + def _prepareReactants(self): if self._Config.GetMixtureStatus(): self._canteraSolution.set_mixture_fraction(self._reactant_mixture_status, self._Config.GetFuelString(), self._Config.GetOxidizerString()) @@ -336,28 +393,28 @@ def _prepareReactants(self): self._canteraSolution.set_equivalence_ratio(self._reactant_mixture_status, self._Config.GetFuelString(), self._Config.GetOxidizerString()) self._canteraSolution.TP = self._T_reactants, self._pressure return - + def _fromRestart(self): if self._flameletSolutionForRestart: self._from_restart = True else: self._from_restart = False return - + def _solverSpecificPreprocessing(self): self._setInitialGrid() self._prepareFlameletSolver() return - + def _setInitialGrid(self): self._initial_grid = np.linspace(0, self._initial_grid_length, self.__initial_grid_number_of_points) return - + def _commonPreprocessing(self): self._flameletSolution.set_refine_criteria(slope=self.__grid_refinement_slope,ratio=self.__grid_refinement_ratio,curve=self.__grid_refinement_curve,prune=self.__grid_refinement_prune) self._flameletSolution.transport_model = self._Config.GetTransportModel() return - + def _computeFlameletSolution(self): try: automatic_grid_refinement = (not self._from_restart) @@ -376,7 +433,7 @@ def _computeFlameletSolution(self): self._converged_solution = False self._flamelet_is_burning = False return - + def _postProcessResults(self): if self.isConverged(): self._extractSolutionDataForOutput() @@ -385,7 +442,7 @@ def _postProcessResults(self): if self._from_file: self._from_file = False return - + def getFlameletSolution(self): """Retrieve Cantera oneDim solution @@ -393,13 +450,13 @@ def getFlameletSolution(self): :rtype: cantera.oneDim """ return self._flameletSolution - + def isConverged(self): return self._converged_solution def isBurning(self): return self._flamelet_is_burning - + def getThermoChemicalData(self): """Retrieve thermochemical state data extracted from flamelet solution. @@ -407,9 +464,9 @@ def getThermoChemicalData(self): :rtype: pandas.DataFrame """ return self._thermochemical_solution - + def _extractSolutionDataForOutput(self): - + self._thermochemical_solution = pd.DataFrame() # Flamelet solution variables @@ -423,7 +480,11 @@ def _extractSolutionDataForOutput(self): for species_index, species_name in enumerate(self._canteraSolution.species_names): self._thermochemical_solution["Y-%s" % species_name] = np.asarray(mass_fractions[species_index]) - + + if self._Config.GetSaveMoleFractions(): + for species_index, species_name in enumerate(self._canteraSolution.species_names): + self._thermochemical_solution["X-%s" % species_name] = np.asarray(molar_fractions[species_index]) + if solution_1D: molecular_weights = self._canteraSolution.molecular_weights[:,np.newaxis] else: @@ -465,7 +526,7 @@ def _extractSolutionDataForOutput(self): mixture_fraction_offset = self._Config.GetMixtureFractionConstant() mixture_fraction = mixture_fraction_offset + np.sum(mixture_fraction_species_coefficients * mass_fractions,axis=0) self._thermochemical_solution[FGMVars.MixtureFraction.name] = mixture_fraction - + temperature = self._flameletSolution.T self._thermochemical_solution[FGMVars.Temperature.name] = temperature @@ -484,12 +545,12 @@ def _extractSolutionDataForOutput(self): dynamic_viscosity = self._flameletSolution.viscosity self._thermochemical_solution[FGMVars.ViscosityDyn.name] = dynamic_viscosity - + heat_release = self._flameletSolution.heat_release_rate self._thermochemical_solution[FGMVars.Heat_Release.name] = heat_release self._writeInflowSettings() - + return def _writeInflowSettings(self): @@ -502,30 +563,30 @@ def _writeInflowSettings(self): reactant_equivalence_ratio = self._reactant_mixture_status gas.set_equivalence_ratio(self._reactant_mixture_status, self._Config.GetFuelString(), self._Config.GetOxidizerString()) reactant_mixture_fraction = gas.mixture_fraction(self._Config.GetFuelString(), self._Config.GetOxidizerString()) - + self._thermochemical_solution["ReactantMixtureFraction"] = reactant_mixture_fraction self._thermochemical_solution["ReactantEquivalenceRatio"] = reactant_equivalence_ratio self._thermochemical_solution["ReactantTemperature"] = self._T_reactants return - + def _extractFlameletDiscretization(self): grid= self._flameletSolution.grid self._thermochemical_solution["Distance"] = grid velocity = self._flameletSolution.velocity self._thermochemical_solution["Velocity"] = velocity return - + def saveSolution(self): flamelet_filename = self.getFlameletFileName() filename_plus_folder = sep.join((self._output_filepath, flamelet_filename)) self._thermochemical_solution.to_csv(filename_plus_folder+".csv",index=False) return filename_plus_folder+".csv" - + def resetRestart(self): self._flameletSolutionForRestart = None self._from_restart = False return - + def loadSolution(self, flameletFileName:str): """Load flamelet solution data from csv file and initialize flamelet simulation from loaded data. @@ -553,7 +614,7 @@ def loadSolution(self, flameletFileName:str): print("Initializing the flamelet solution from data frame will be included in the upcoming Cantera version.") return - + def _prepareInitialGuessData(self): initialGuessData = pd.DataFrame() initialGuessData["temperature"] = self._thermochemical_solution[FGMVars.Temperature.name] @@ -563,13 +624,13 @@ def _prepareInitialGuessData(self): for sp in self._canteraSolution.species_names: initialGuessData["Y_%s" % sp] = self._thermochemical_solution["Y-%s" % sp] return initialGuessData - + def isPremixed(self): return self._is_premixed - + def isScalar(self): return self._is_scalar - + def getMassFractions(self): """Retrieve mass fraction data from flamelet solution. @@ -580,8 +641,8 @@ def getMassFractions(self): species_names = self._canteraSolution.species_names struct = {sp : y for sp, y in zip(species_names, Y)} return pd.DataFrame(data=struct) - - + + class FreeFlameSolver(FlameletSolver_Cantera): """Solver class for adiabatic free flamelets """ @@ -599,40 +660,40 @@ def __init__(self, config_input:Config_FGM): self.setInitialGrid(0.2, 50) self._initializeFlameletSolver() return - + def _initializeFlameletSolver(self): self._flameletSolution = ct.FreeFlame(self._canteraSolution, grid=self._initial_grid) return - + def _solverSpecificPreprocessing(self): super()._solverSpecificPreprocessing() if not self._from_restart and not self._from_file: self._flameletSolution.set_initial_guess(locs=[0.0, 0.3, 0.5, 1.0]) self._flameletSolution.inlet.T = self._T_reactants return - + def _prepareSettingRange(self): super()._prepareSettingRange() Tu_bounds = self._Config.GetUnbTempBounds() Tu_range = np.linspace(Tu_bounds[1], Tu_bounds[0], self._n_1D_iterations) return Tu_range - + def _parseInputSettings(self, flamelet_solution_settings): if "mixture_status" in flamelet_solution_settings.keys(): self.setMixtureStatus(flamelet_solution_settings["mixture_status"]) if "reactant_temperature" in flamelet_solution_settings.keys(): self.setReactantTemperature(flamelet_solution_settings["reactant_temperature"]) return super()._parseInputSettings(flamelet_solution_settings) - + def _writeSolverSettings(self, solution_index:int, setting_1D:float): super()._writeSolverSettings(solution_index, setting_1D) freeflame_settings = {"reactant_temperature":setting_1D, "solution_index":solution_index} return freeflame_settings - + def setInputVariable(self, val_input:float): self.setReactantTemperature(val_input) return - + def getFlameletFileName(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="mixfrac" @@ -643,7 +704,7 @@ def getFlameletFileName(self): self._reactant_mixture_status, \ self._T_reactants) return flamelet_filename - + def _printStatusToTerminal(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="Z" @@ -665,18 +726,18 @@ def _printStatusToTerminal(self): self._reactant_mixture_status,\ self._T_reactants,\ self._thermochemical_solution.shape[0]) - + if self._print_iteration: outp_message += "(%i/%i)" % (self._iteration+1, self._n_1D_iterations) print(outp_message) return - + def _postProcessResults(self): super()._postProcessResults() if self.isBurning() and self.isConverged(): self.__mass_flow_rate = self._thermochemical_solution["Velocity"][0] * self._thermochemical_solution["Density"][0] return - + def getMassFlowRate(self): """Retrieve adiabatic mass flow rate. @@ -684,7 +745,7 @@ def getMassFlowRate(self): :rtype: float """ return self.__mass_flow_rate - + class BurnerFlameSolver(FlameletSolver_Cantera): """Solver class for burner-stabilized flamelets """ @@ -711,7 +772,7 @@ def __init__(self, config_input:Config_FGM): self.setGridRefinementCriteria(ratio=3, slope=0.15, curve=0.15, prune=0.05) self._initializeFlameletSolver() return - + def setReactantMassFlow(self, val_massflow_inlet:float): """Specify the mass flow rate at the inflow boundary. @@ -723,7 +784,7 @@ def setReactantMassFlow(self, val_massflow_inlet:float): raise Exception("Mass flow rate should be strictly positive.") self.__val_massflow = val_massflow_inlet return - + def getReactantMassFlow(self): return self.__val_massflow @@ -732,12 +793,12 @@ def _solverSpecificPreprocessing(self): self._flameletSolution.burner.mdot = self.__val_massflow self._flameletSolution.burner.T = self._T_reactants return - + def _initializeFlameletSolver(self): self._flameletSolution = ct.BurnerFlame(self._canteraSolution, self._initial_grid) return - + def retrieveSolverSettings(self, solvers:Dict[str, FlameletSolver_Cantera]): if "FREEFLAME" in solvers.keys(): freeflame_solver:FreeFlameSolver = solvers["FREEFLAME"] @@ -752,13 +813,13 @@ def retrieveSolverSettings(self, solvers:Dict[str, FlameletSolver_Cantera]): else: raise Exception("Unable to calculate adiabatic mass flow rate") return - + def _loadSolverSpecificData(self): u = self._thermochemical_solution[FGMVars.Velocity.name][0] rho = self._thermochemical_solution[FGMVars.Density.name][0] self.setReactantMassFlow(u * rho) return - + def _prepareSettingRange(self): super()._prepareSettingRange() mdot_max = 0.98 * self.__adiabatic_massflow @@ -769,7 +830,7 @@ def _prepareSettingRange(self): self.__val_massflow_enthalpy = mdot_max m_dot_range = np.linspace(mdot_max, mdot_min, self._n_1D_iterations+1)[:-1] return m_dot_range - + def _postProcessResults(self): super()._postProcessResults() if self.__iterate_enthalpy(): @@ -783,14 +844,14 @@ def _postProcessResults(self): self.__val_massflow_enthalpy -= self.__delta_massflow self._keep_iterating = (self.__val_massflow_enthalpy > 0.001*self.__adiabatic_massflow) return - + def __iterate_enthalpy(self): return self.__delta_enth is not None - + def setInputVariable(self, val_input:float): self.setReactantMassFlow(val_input) return - + def _parseInputSettings(self, flamelet_solution_settings): if "mixture_status" in flamelet_solution_settings.keys(): self.setMixtureStatus(flamelet_solution_settings["mixture_status"]) @@ -803,14 +864,14 @@ def _parseInputSettings(self, flamelet_solution_settings): self.setReactantMassFlow(self.__val_massflow_enthalpy) return super()._parseInputSettings(flamelet_solution_settings) - + def _writeSolverSettings(self, solution_index:int, setting_1D:float): super()._writeSolverSettings(solution_index, setting_1D) freeflame_settings = {"mdot":setting_1D} return freeflame_settings - - - + + + def getFlameletFileName(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="mixfrac" @@ -821,7 +882,7 @@ def getFlameletFileName(self): self._reactant_mixture_status,\ self.__val_massflow) return flamelet_filename - + def _printStatusToTerminal(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="Z" @@ -848,7 +909,7 @@ def _printStatusToTerminal(self): print(output_message) return - + def loadSolution(self, flameletFileName): super().loadSolution(flameletFileName) u = self._thermochemical_solution[FGMVars.Velocity.name][0] @@ -875,7 +936,7 @@ def __init__(self, config_input:Config_FGM): self._n_1D_iterations = self._Config.GetNpTemp() self._flameletFileExtension = "csv" return - + def solveForMixtureStatus(self, val_mixture_status:float, save:bool=True): self.__is_reaction_products = False self.__accumulated_solution = pd.DataFrame() @@ -889,28 +950,28 @@ def solveForMixtureStatus(self, val_mixture_status:float, save:bool=True): self.resetRestart() return - + def _solverSpecificPreprocessing(self): self._flameletSolution = self._canteraSolution return - + def _initializeFlameletSolver(self): self._flameletSolution = self._canteraSolution return - + def _parseInputSettings(self, flamelet_solution_settings): if "mixture_status" in flamelet_solution_settings.keys(): self.setMixtureStatus(flamelet_solution_settings["mixture_status"]) if "reactant_temperature" in flamelet_solution_settings.keys(): self.setReactantTemperature(flamelet_solution_settings["reactant_temperature"]) return super()._parseInputSettings(flamelet_solution_settings) - - + + def _writeSolverSettings(self, solution_index:int, setting_1D:float): super()._writeSolverSettings(solution_index, setting_1D) freeflame_settings = {"reactant_temperature":setting_1D} return freeflame_settings - + def _prepareSettingRange(self): super()._prepareSettingRange() @@ -933,45 +994,45 @@ def _prepareSettingRange(self): T_range = np.linspace(T_min, T_max, self._n_1D_iterations) return T_range - + def setInputVariable(self, val_input:float): self.setReactantTemperature(val_input) return - + def _computeFlameletSolution(self): - + try: self._canteraSolution.TP = self._T_reactants, ct.one_atm self._converged_solution = True except: self._converged_solution = False return - + def _extractFlameletDiscretization(self): self._thermochemical_solution["Distance"] = np.zeros(1) self._thermochemical_solution["Velocity"] = np.zeros(1) return - + def concatenateSolution(self, other_solution:pd.DataFrame): concatenated_data = pd.concat((other_solution, self._thermochemical_solution),axis=0) return concatenated_data - + def reactionProducts(self, is_reaction_products:bool=True): self.__is_reaction_products = is_reaction_products return - + def isReactionProducts(self): return self.__is_reaction_products - + def _commonPreprocessing(self): return - + def getFlameletFileName(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="mixfrac" else: tag_for_mixture_status="phi" - + if self.__is_reaction_products: file_header = "Products" else: @@ -980,7 +1041,7 @@ def getFlameletFileName(self): tag_for_mixture_status, \ self._reactant_mixture_status) return flamelet_filename - + def _postProcessResults(self): super()._postProcessResults() if self.isConverged(): @@ -989,18 +1050,18 @@ def _postProcessResults(self): else: self.__accumulated_solution = pd.concat((self.__accumulated_solution, self._thermochemical_solution),axis=0) return - + def _writeOutput(self): if self.isConverged(): self.saveSolution() return - + def saveSolution(self): flamelet_filename = self.getFlameletFileName() filename_plus_folder = sep.join((self._output_filepath, flamelet_filename)) self.__accumulated_solution.to_csv("%s.%s" % (filename_plus_folder, self._flameletFileExtension),index=False) return - + def loadSolution(self, flameletFileName:str): solution_fileName = flameletFileName.split(sep)[-1] self.__is_reaction_products = ("Products" in solution_fileName) @@ -1009,17 +1070,17 @@ def loadSolution(self, flameletFileName:str): self.__solution_products = self.__accumulated_solution else: self.__solution_reactants = self.__accumulated_solution - + val_mixfrac = np.clip(self.__accumulated_solution[FGMVars.MixtureFraction.name].iloc[-1], 0.0, 1.0) self._canteraSolution.set_mixture_fraction(val_mixfrac, self._Config.GetFuelString(), self._Config.GetOxidizerString()) if self._Config.DefineMixtureStatus(): self._reactant_mixture_status = self._canteraSolution.mixture_fraction(self._Config.GetFuelString(), self._Config.GetOxidizerString()) else: self._reactant_mixture_status = self._canteraSolution.equivalence_ratio(self._Config.GetFuelString(), self._Config.GetOxidizerString()) - + self.setReactantTemperature(self.__accumulated_solution[FGMVars.Temperature.name][0]) return - + def _fromRestart(self): return @@ -1034,10 +1095,10 @@ def _prepareReactants(self): self._canteraSolution.equilibrate("TP") else: self._canteraSolution.equilibrate("HP") - + self._canteraSolution.TP = self._T_reactants, self._pressure return - + def setMixtureStatus(self, val_mixture_status:float): super().setMixtureStatus(val_mixture_status) self.__is_lean = False @@ -1052,12 +1113,12 @@ def setMixtureStatus(self, val_mixture_status:float): def getThermoChemicalData(self): return self.__accumulated_solution - + def getReactionProductData(self): return self.__solution_products def getReactantData(self): return self.__solution_reactants - + class CooledFlameInterpolator(FlameletSolver_Cantera): """Class for interpolated data between burner-stabilized flamelet and chemical equilibrium data """ @@ -1074,24 +1135,24 @@ def __init__(self, config:Config_FGM): self._n_1D_iterations = self._Config.GetNpMdotExtra() self._flameletFileExtension = "csv" return - + def setEquilibriumData(self, eq_data:pd.DataFrame): self.__equilibriumSolution = eq_data return - + def setBurnerFlameData(self, burner_data:pd.DataFrame): self.__burnerFlameSolution = burner_data return - + def _prepareSettingRange(self): super()._prepareSettingRange() iter_values = [i for i in range(self._n_1D_iterations)] return iter_values - + def setInputVariable(self, val_input:float): self._iteration = val_input return - + def _computeFlameletSolution(self): ratio = float(self._iteration + 1) / float(self._n_1D_iterations) w_a_lin = 1.0 - ratio @@ -1112,21 +1173,21 @@ def _computeFlameletSolution(self): else: self._thermochemical_solution[var] = linear_interpolation[:, iVar] return - + def _writeOutput(self): if self.isConverged(): self.saveSolution() return - + def saveSolution(self): flamelet_filename = self.getFlameletFileName() filename_plus_folder = sep.join((self._output_filepath, flamelet_filename)) self._thermochemical_solution.to_csv("%s.%s" % (filename_plus_folder, self._flameletFileExtension),index=False) return - + def _extractSolutionDataForOutput(self): return - + def getFlameletFileName(self): if self._Config.GetMixtureStatus(): tag_for_mixture_status="mixfrac" @@ -1137,21 +1198,21 @@ def getFlameletFileName(self): self._reactant_mixture_status, \ self._iteration) return flamelet_filename - + def retrieveSolverSettings(self, solvers:Dict[str, FlameletSolver_Cantera]): burnerflameSolver = solvers["BURNERFLAME"] equilibriumSolver = solvers["EQUILIBRIUM"] self.setBurnerFlameData(burnerflameSolver.getThermoChemicalData()) self.setEquilibriumData(equilibriumSolver.getReactionProductData()) return - + def _printStatusToTerminal(self): message_out = "Interpolated between burner-stabilized and equilibrium data " if self._print_iteration: message_out += "(%i/%i)" % (self._iteration+1, self._n_1D_iterations) print(message_out) return - + def _solverSpecificPreprocessing(self): return def _commonPreprocessing(self): @@ -1159,7 +1220,7 @@ def _commonPreprocessing(self): def _postProcessResults(self): return - + def _writeSolverSettings(self, solution_index:int, setting_1D:float): super()._writeSolverSettings(solution_index, setting_1D) return {"iteration":setting_1D} @@ -1176,7 +1237,7 @@ def loadSolution(self, flameletFileName:str): else: self._reactant_mixture_status = self._canteraSolution.equivalence_ratio(self._Config.GetFuelString(), self._Config.GetOxidizerString()) self.setReactantTemperature(self._thermochemical_solution[FGMVars.Temperature.name][0]) - + return class CounterFlowDiffusionFlameSolver(FlameletSolver_Cantera): @@ -1195,17 +1256,23 @@ def __init__(self, config_input): self._plotLabel = "Counter-flow diffusion flame" self._is_premixed = False self._is_scalar = False + self._n_1D_iterations = self._Config.GetNpTemp() self.setInitialGrid(2e-1, 300) self.setGridRefinementCriteria(ratio=3, slope=0.04, curve=0.06, prune=0.02) + + # Retrieve fixed-strain settings from Config + self.__fixed_strain_mode = self._Config.GetCounterFlowFixedStrain() + if self.__fixed_strain_mode: + self.setStrainRate(self._Config.GetCounterFlowStrainRate()) return - + def _initializeFlameletSolver(self): self._flameletSolution = ct.CounterflowDiffusionFlame(self._canteraSolution, grid=self._initial_grid) return - + def _prepareFlameletSimulation(self): return super()._prepareFlameletSimulation() - + def _prepareReactants(self): super()._prepareReactants() self._canteraSolution.set_mixture_fraction(0.0, self._Config.GetFuelString(), self._Config.GetOxidizerString()) @@ -1218,7 +1285,7 @@ def _prepareReactants(self): self.__oxidizer_velocity = self.__fuel_velocity * self.__fuel_density / self.__oxidizer_density return - + def _solverSpecificPreprocessing(self): super()._solverSpecificPreprocessing() self._flameletSolution.P = ct.one_atm @@ -1231,19 +1298,19 @@ def _solverSpecificPreprocessing(self): self._flameletSolution.oxidizer_inlet.mdot = self.__oxidizer_density * self.__oxidizer_velocity self._flameletSolution.oxidizer_inlet.Y = self._Config.GetOxidizerString() return - + def setStrainRate(self, val_strain_rate:float=1.0): if val_strain_rate <= 0.0: raise Exception("Strain rate should be strictly positive") self.__strain_rate = val_strain_rate return - + def getFlameletFileName(self): flamelet_filename = "%s_strain%.3e_Tu%.1f" % (self._flamelet_type, \ self.__strain_rate, \ self._T_reactants) return flamelet_filename - + def _printStatusToTerminal(self): if not self.isConverged(): output_message = "%s simulation at strain=%.2f s^-1 did not converge " % (self._flamelet_type,\ @@ -1261,13 +1328,13 @@ def _printStatusToTerminal(self): print(output_message) return - + def _writeInflowSettings(self): super()._writeInflowSettings() self._thermochemical_solution["StrainRate"] = self.__strain_rate self._thermochemical_solution["SpreadRate"] = self._flameletSolution.spread_rate return - + def loadSolution(self, flameletFileName:str): super().loadSolution(flameletFileName) grid_length = self._thermochemical_solution["Distance"].iloc[-1] - self._thermochemical_solution["Distance"].iloc[0] @@ -1275,36 +1342,39 @@ def loadSolution(self, flameletFileName:str): self.setInitialGrid(grid_length, n_nodes) self.__strain_rate = self._thermochemical_solution["StrainRate"][0] return - + def _prepareInitialGuessData(self): initialGuessData = super()._prepareInitialGuessData() initialGuessData["spreadRate"] = self._thermochemical_solution["SpreadRate"] return initialGuessData - + def _parseInputSettings(self, flamelet_solution_settings): if "strain_rate" in flamelet_solution_settings.keys(): self.setStrainRate(flamelet_solution_settings["strain_rate"]) if "reactant_temperature" in flamelet_solution_settings.keys(): self.setReactantTemperature(flamelet_solution_settings["reactant_temperature"]) return super()._parseInputSettings(flamelet_solution_settings) - - + + def _writeSolverSettings(self, solution_index:int, setting_1D:float): super()._writeSolverSettings(solution_index, setting_1D) freeflame_settings = {"reactant_temperature":setting_1D} + # Pass strain rate when in fixed-strain mode + if hasattr(self, '_CounterFlowDiffusionFlameSolver__fixed_strain_mode') and self.__fixed_strain_mode: + freeflame_settings["strain_rate"] = self.__strain_rate return freeflame_settings - + def _prepareSettingRange(self): super()._prepareSettingRange() Tu_bounds = self._Config.GetUnbTempBounds() Tu_range = np.linspace(Tu_bounds[0], Tu_bounds[1], self._n_1D_iterations) return Tu_range - + def setInputVariable(self, val_input:float): self.setReactantTemperature(val_input) return - + def _prepareStorageFolder(self): folder_for_flamelet_type = self._createFolderForFlameletType() @@ -1313,14 +1383,14 @@ def _prepareStorageFolder(self): filepath_for_flamelet_data = sep.join((folder_for_flamelet_type, mixture_subfolder)) if not path.isdir(filepath_for_flamelet_data): mkdir(filepath_for_flamelet_data) - + self._output_filepath = filepath_for_flamelet_data return - + def __createSubFolderForStrain(self): strain_subfolder = "strain_%.3e" % (self.__strain_rate) return strain_subfolder - + FlameletSolverDict:dict = {"FREEFLAME" : FreeFlameSolver,\ "BURNERFLAME" : BurnerFlameSolver,\ "EQUILIBRIUM" : EquilibriumSolver,\ diff --git a/Data_Processing/collectFlameletData.py b/Data_Processing/collectFlameletData.py index 81e93fa..5eb20e3 100644 --- a/Data_Processing/collectFlameletData.py +++ b/Data_Processing/collectFlameletData.py @@ -149,10 +149,10 @@ def ConcatenateFlameletData(self): self.__SizeDataArrays() self.__extractFlameletData() - + self.__WriteOutputFiles() return - + def CollectBoundaryData(self): self.IgnoreMixtureBounds(True) self.__Config.setFlameletTypes(["EQUILIBRIUM"]) @@ -181,13 +181,13 @@ def __SizeDataArrays(self): self.__LookUp_flamelet_data = np.zeros([self.__nFlamelets * self.__Np_per_flamelet, len(self.__LookUp_vars)]) return - + def __countNumberofFlameletDataPoints(self): self.__nFlameletDataPoints = 0 self.__nFlamelets = 0 self.__loopOverFlamelets(self.__incrementNumberOfFlameletData) return - + def __loopOverFlamelets(self, taskPerFlamelet:Callable): flameletTypes = self.__Config.getFlameletTypes() for flameletType in flameletTypes: @@ -207,7 +207,7 @@ def __loopOverFlamelets(self, taskPerFlamelet:Callable): taskPerFlamelet(flameletSolver) self.__printMsg("Done.") return - + def __isWithinMixtureBounds(self, flameletSolver:FlameletSolver_Cantera): if self.__ignore_mixture_bounds or not flameletSolver.isPremixed(): return True @@ -215,13 +215,13 @@ def __isWithinMixtureBounds(self, flameletSolver:FlameletSolver_Cantera): mixture_status = flameletSolver.getMixtureStatus() margin = 1e-2 return (mixture_status-margin <= self.__mix_status_max) and (mixture_status+margin >= self.__mix_status_min) - + def __incrementNumberOfFlameletData(self, flameletSolver:FlameletSolver_Cantera): thermochemical_solution = flameletSolver.getThermoChemicalData() self.__nFlameletDataPoints += thermochemical_solution.shape[0] self.__nFlamelets += 1 return - + def __extractFlameletData(self): self.__printMsg("Extracting thermochemical data from manifold...") self.__flameletSolutionIndex = 0 @@ -236,7 +236,7 @@ def __interpolateAlongFlamelet(self, flameletSolution:FlameletSolver_Cantera): T_flamelet = solutionData[FGMVars.Temperature.name] if np.max(T_flamelet) < DefaultSettings_FGM.T_threshold: flameletIsBurning = False - + if flameletIsBurning: tracingVariable = self.__calculateTracingVariable(solutionData) if self.__reactionProductsForLUT(flameletSolution): @@ -272,7 +272,7 @@ def __interpolateAlongFlamelet(self, flameletSolution:FlameletSolver_Cantera): self.__PD_flamelet_data[startIndex:stopIndex] = PD_data_interpolated self.__flameletSolutionIndex += 1 return - + def __retrieveFlameletSolution(self, flameletSolver:FlameletSolver_Cantera): solutionData = flameletSolver.getThermoChemicalData() @@ -282,7 +282,7 @@ def __retrieveFlameletSolution(self, flameletSolver:FlameletSolver_Cantera): if not flameletSolver.isPremixed() and not self.__ignore_mixture_bounds: solutionData = self.__clipNonPremixedFlameletToMixtureBounds(solutionData) return solutionData - + def __clipNonPremixedFlameletToMixtureBounds(self, solutionData:pd.DataFrame): if self.__Config.GetMixtureStatus(): mixfrac_upper = self.__mix_status_max @@ -298,12 +298,12 @@ def __clipNonPremixedFlameletToMixtureBounds(self, solutionData:pd.DataFrame): within_bounds = np.logical_and(mixfrac_solution >= mixfrac_lower, mixfrac_solution <= mixfrac_upper) solutionData_out = solutionData.iloc[within_bounds] return solutionData_out - + def __WriteOutputFiles(self): """Collect all flamelet data arrays, split into train, test, and validation portions, and write to appropriately named files. """ self.__printMsg("Writing accumulated flamelet data files...") - + flameletData = self.__accumulateDataFromFlamelets() flameletDataFull = self.__filterFlameletData(flameletData) @@ -324,7 +324,7 @@ def __WriteOutputFiles(self): self.__printMsg("Done.") return - + def __accumulateDataFromFlamelets(self): outputData = pd.DataFrame() outputData[self.__Config.GetControllingVariables()] = self.__CV_flamelet_data @@ -334,25 +334,25 @@ def __accumulateDataFromFlamelets(self): if self.__Config.PreferentialDiffusion(): outputData[self.__PD_train_vars] = self.__PD_flamelet_data return outputData - + def __filterFlameletData(self, flameletData:pd.DataFrame): uniqueData = flameletData.drop_duplicates() noNans = uniqueData.dropna() noZeros = noNans.loc[~(noNans == 0).all(axis=1)].reset_index(drop=True) - + return noZeros - + def __splitFlameletDataSet(self, flameletDataFull:pd.DataFrame): shuffledData = flameletDataFull.sample(frac=1).reset_index(drop=True) Np_total = shuffledData.shape[0] Np_train = int(self.__Config.GetTrainFraction()*Np_total) Np_test = int(self.__Config.GetTestFraction()*Np_total) - + trainData = shuffledData.iloc[:Np_train,:] testData = shuffledData.iloc[Np_train:Np_train+Np_test,:] validationData = shuffledData.iloc[Np_train+Np_test:,:] return trainData, testData, validationData - + def IgnoreMixtureBounds(self, ignore_bounds:bool=False): self.__ignore_mixture_bounds = ignore_bounds return @@ -453,7 +453,7 @@ def IncludeFlameletType(self, flameletType:str, include:bool=True): else: self.__Config.excludeFlameletType(flameletType) return - + def SetLookUpVars(self, input:list[str]): """Define passive look-up variables to be included in the manifold data. @@ -461,6 +461,9 @@ def SetLookUpVars(self, input:list[str]): :type input: list[str] """ self.__Config.SetLookUpVariables(input) + self.__LookUp_vars = [] + for var in input: + self.__LookUp_vars.append(var) return def SetFlameletDir(self, input:str): @@ -506,25 +509,25 @@ def SetTestFraction(self, input:float=DefaultSettings_FGM.test_fraction): self.__Config.SetTestFraction(input) return - + def __reactionProductsForLUT(self, flameletSolver:FlameletSolver_Cantera): if self.__write_LUT_data and flameletSolver.getFlameletType()=="Equilibrium": return flameletSolver.isReactionProducts() else: return False - + def __calculateTracingVariable(self, solutionData:pd.DataFrame): controlVariables = self.__retrieveControlVariables(solutionData) - + cv_max, cv_min = np.max(controlVariables,axis=0), np.min(controlVariables,axis=0) scaledControlVariables = (controlVariables - cv_min)/(cv_max - cv_min + 1e-10) controlVariableIncrement = scaledControlVariables[1:] - scaledControlVariables[:-1] - + tracingVariableIncrement = np.linalg.norm(controlVariableIncrement,axis=1) tracingVariable = np.hstack((0, np.cumsum(tracingVariableIncrement))) tracingVariable_scaled = tracingVariable / (np.max(tracingVariable)+1e-10) return tracingVariable_scaled - + def __retrieveControlVariables(self, solutionData:pd.DataFrame): controlVariables = np.zeros([solutionData.shape[0], len(self.__Config.GetControllingVariables())]) for iCv, cv in enumerate(self.__Config.GetControllingVariables()): @@ -533,7 +536,7 @@ def __retrieveControlVariables(self, solutionData:pd.DataFrame): else: controlVariables[:, iCv] = solutionData[cv] return controlVariables - + def __retrieveThermoPhyiscalData(self, solutionData:pd.DataFrame): TD_data = np.zeros([solutionData.shape[0], len(self.__TD_train_vars)]) for iVar_TD, TD_var in enumerate(self.__TD_train_vars): @@ -545,13 +548,13 @@ def __retrieveThermoPhyiscalData(self, solutionData:pd.DataFrame): else: TD_data[:, iVar_TD] = solutionData[TD_var] return TD_data - + def __retrievePassiveLookUpData(self, solutionData:pd.DataFrame): LookUp_data = np.zeros([solutionData.shape[0], len(self.__LookUp_vars)]) for iVar_LookUp, LookUp_var in enumerate(self.__LookUp_vars): LookUp_data[:, iVar_LookUp] = solutionData[LookUp_var] return LookUp_data - + def __retrievePreferentialDiffusionScalars(self, solutionData:pd.DataFrame): vars = list(solutionData.keys()) beta_pv_flamelet, beta_h1_flamelet, beta_h2_flamelet, beta_z_flamelet = self.__Config.ComputeBetaTerms(vars, solutionData.values) @@ -561,7 +564,7 @@ def __retrievePreferentialDiffusionScalars(self, solutionData:pd.DataFrame): PD_data[:, 2] = beta_h2_flamelet PD_data[:, 3] = beta_z_flamelet return PD_data - + def __retrieveSourceTerms(self, solutionData:pd.DataFrame): nP_flamelet = solutionData.shape[0] species_mass_fraction = np.zeros([nP_flamelet, len(self.__Species_in_FGM)]) @@ -594,7 +597,7 @@ def __retrieveSourceTerms(self, solutionData:pd.DataFrame): Sources_data[:, 1 + 4*iSp + 1] = species_destruction_rate[:, iSp] Sources_data[:, 1 + 4*iSp + 2] = species_net_rate[:, iSp] Sources_data[:, 1 + 4*iSp + 3] = species_mass_fraction[:, iSp] - + sourceterm_zero_line_numbers = np.zeros(nP_flamelet, dtype=bool) sourceterm_zero_line_numbers[0] = True sourceterm_zero_line_numbers[-1] = True @@ -609,25 +612,25 @@ def __retrieveSourceTerms(self, solutionData:pd.DataFrame): sourceterm_zero_line_numbers, np.logical_or((T_flamelet - T_min) < deltaT, (T_max - T_flamelet) < deltaT)) - + # Only zero the PV source term at the flamelet boundaries. # Species production rates and mass fractions are NOT zeroed # because (a) Y-{species} is non-zero at PV_max and (b) species # rates at PV_max are governed by Cantera's equilibrium values. Sources_data[sourceterm_zero_line_numbers, 0] = 0.0 return Sources_data - + def __interpolateAlongTracingVariable(self, tracingVariableData:np.ndarray[float], tracingVariableQuery:np.ndarray[float], flameletData:np.ndarray[float]): flameletData_Sampled = np.zeros([self.__Np_per_flamelet, np.shape(flameletData)[1]]) for i in range(np.shape(flameletData)[1]): flameletData_Sampled[:, i] = np.interp(tracingVariableQuery, tracingVariableData, flameletData[:, i]) return flameletData_Sampled - + def __printMsg(self, msg:str): if self.__verbose > 0: print(msg) return - + class GroupOutputs: """Class which groups flamelet data variables into MLP outputs based on their affinity. """ diff --git a/Manifold_Generation/LUT/FlameletTableGeneration.py b/Manifold_Generation/LUT/FlameletTableGeneration.py index f2a02fa..d134dc6 100644 --- a/Manifold_Generation/LUT/FlameletTableGeneration.py +++ b/Manifold_Generation/LUT/FlameletTableGeneration.py @@ -590,38 +590,196 @@ def __EvaluateFlameletInterpolator(self, CV_unscaled:np.ndarray): data_interp = self._lookup_tree(q=CV_scaled,nnear=self._n_near,p=self._p_fac) return data_interp - def VisualizeTableLevel(self, val_mix_frac:float, var_to_plot:str=None): - """Compute and visualize the table connectivity for a certain mixture fraction value. + def PlotTableSlices(self, + x_cv: str, + y_var: str, + slice_cv: str, + slice_range: tuple = None, + n_slices: int = 10, + n_x_points: int = 200, + save_path: str = None, + show: bool = False): + """Plot a dependent variable against one controlling variable for a set of + fixed slices of another controlling variable, querying the IDW interpolator. + + Example: T(Z) for 10 enthalpy levels between h_min and h_max. + + :param x_cv: Name of the controlling variable to use as the x-axis (e.g. 'MixtureFraction'). + :type x_cv: str + :param y_var: Name of the variable to plot on the y-axis (e.g. 'Temperature'). + :type y_var: str + :param slice_cv: Name of the controlling variable held fixed per line (e.g. 'EnthalpyTot'). + :type slice_cv: str + :param slice_range: (min, max) values of slice_cv. If None, the data extent is used. + :type slice_range: tuple, optional + :param n_slices: Number of slice values uniformly distributed across slice_range. + :type n_slices: int + :param n_x_points: Number of x_cv sample points per slice. + :type n_x_points: int + :param save_path: If provided, save the figure to this path instead of displaying it. + :type save_path: str, optional + :raises Exception: if x_cv, slice_cv, or y_var are not found in the flamelet data. + """ + for name in [x_cv, slice_cv]: + if name not in self._controlling_variables: + raise Exception("'%s' is not a controlling variable. Available: %s" + % (name, self._controlling_variables)) + if y_var not in self._Flamelet_Variables: + raise Exception("'%s' not found in flamelet data. Available: %s" + % (y_var, self._Flamelet_Variables)) + + x_idx = self._controlling_variables.index(x_cv) + slice_idx = self._controlling_variables.index(slice_cv) + y_col = self._Flamelet_Variables.index(y_var) + + # Determine x and slice extents from the loaded data. + x_min = float(np.min(self._D_full[:, x_idx])) + x_max = float(np.max(self._D_full[:, x_idx])) + if slice_range is None: + s_min = float(np.min(self._D_full[:, slice_idx])) + s_max = float(np.max(self._D_full[:, slice_idx])) + else: + s_min, s_max = float(slice_range[0]), float(slice_range[1]) + + slice_values = np.linspace(s_min, s_max, n_slices) + x_values = np.linspace(x_min, x_max, n_x_points) + + # Level CV value: 0.0 for 2D mode (PV = 0), or the table centre otherwise. + n_cv = len(self._controlling_variables) + level_val = 0.0 if self._is_2D_table else 0.5 * ( + float(np.min(self._D_full[:, self._level_cv_idx])) + + float(np.max(self._D_full[:, self._level_cv_idx]))) + + # Build a Delaunay hull of the actual 2D data so that query points outside + # the data support (e.g. high h at Z=1 for counterflow flames) are masked. + from scipy.spatial import Delaunay as _Delaunay + if self._is_2D_table and self._scaler_2d is not None: + ZH_dim = self._D_full[:, [x_idx, slice_idx]] + ZH_norm = self._scaler_2d.transform(ZH_dim) + _data_hull = _Delaunay(ZH_norm) + _has_hull = True + _x_scaler_range = [float(self._D_full[:, x_idx].min()), + float(self._D_full[:, x_idx].max())] + _s_scaler_range = [float(self._D_full[:, slice_idx].min()), + float(self._D_full[:, slice_idx].max())] + else: + _has_hull = False + + cmap = plt.cm.coolwarm + norm = plt.Normalize(vmin=s_min, vmax=s_max) + + fig, ax = plt.subplots(figsize=(9, 5), constrained_layout=True) + + for s_val in slice_values: + # Build query array: all CVs at their level value, then overwrite x and slice. + CV_query = np.zeros([n_x_points, n_cv]) + CV_query[:, self._level_cv_idx] = level_val + CV_query[:, x_idx] = x_values + CV_query[:, slice_idx] = s_val + + # Mask x values that fall outside the convex hull of the data. + if _has_hull: + x_norm = (x_values - _x_scaler_range[0]) / (_x_scaler_range[1] - _x_scaler_range[0] + 1e-32) + s_norm = (s_val - _s_scaler_range[0]) / (_s_scaler_range[1] - _s_scaler_range[0] + 1e-32) + query_norm = np.column_stack([x_norm, np.full(n_x_points, s_norm)]) + inside = _data_hull.find_simplex(query_norm) >= 0 + else: + inside = np.ones(n_x_points, dtype=bool) + + if not np.any(inside): + continue # entire slice is outside data support — skip - :param val_mix_frac: mixture fraction value for which to compute the table connectivity. + result = self.__EvaluateFlameletInterpolator(CV_query) + y_vals = result[:, y_col] + y_vals[~inside] = np.nan # blank extrapolation regions + + ax.plot(x_values, y_vals, color=cmap(norm(s_val)), linewidth=1.2) + + sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + cb = fig.colorbar(sm, ax=ax, pad=0.02) + cb.set_label(slice_cv, fontsize=11) + + ax.set_xlabel(x_cv, fontsize=12) + ax.set_ylabel(y_var, fontsize=12) + ax.set_title("%s(%s) | %d slices of %s" % (y_var, x_cv, n_slices, slice_cv), fontsize=13) + ax.grid(True, linestyle='--', alpha=0.4) + + if save_path: + fig.savefig(save_path, dpi=150, bbox_inches='tight') + print(" Saved: %s" % save_path) + if not show: + plt.close(fig) + return + + def VisualizeTableLevel(self, val_mix_frac:float, var_to_plot:str=None, + plot_3d:bool=False, show_grid:bool=True, + save_path:str=None, show:bool=False): + """Compute and visualize the table mesh and optionally a data field at a given level value. + + :param val_mix_frac: value of the level controlling variable for which to generate the table level. :type val_mix_frac: float - :raises Exception: if the mixture fraction value lies outside the flamelet data range. + :param var_to_plot: name of the variable to colour the plot with. If None, only the mesh is shown. + :type var_to_plot: str, optional + :param plot_3d: when True and var_to_plot is set, render a 3D surface; when False render a 2D colour map. + :type plot_3d: bool + :param show_grid: overlay the triangulation wireframe on the plot. + :type show_grid: bool + :param save_path: if provided, save the figure to this path instead of displaying it. + :type save_path: str, optional """ Tria, Nodes, HullIdx, level_data, XY_ref_dim = self.ComputeTableLevelMesh(val_mix_frac) print("Total mesh nodes: %i, refinement seed points: %i" % (len(Nodes), len(XY_ref_dim))) p0, p1 = self._plane_cv_idxs - if var_to_plot == None: - _ = plt.figure(figsize=[10,10]) - ax = plt.axes() - ax.triplot(Nodes[:, p0], Nodes[:, p1], Tria) - ax.plot(Nodes[HullIdx, p0], Nodes[HullIdx, p1], 'ko', label=r"Hull nodes") - ax.set_xlabel(self._plane_cv_names[0], fontsize=20) - ax.set_ylabel(self._plane_cv_names[1], fontsize=20) - ax.legend(fontsize=20) - ax.set_title(self._level_cv_name + " = " + str(val_mix_frac)) - plt.show() - else: + title = "%s = %s" % (self._level_cv_name, str(val_mix_frac)) + + if var_to_plot is None: + # Mesh-only plot (always 2D). + fig, ax = plt.subplots(figsize=(10, 10)) + ax.triplot(Nodes[:, p0], Nodes[:, p1], Tria, linewidth=0.5) + ax.plot(Nodes[HullIdx, p0], Nodes[HullIdx, p1], 'ko', ms=3, label="Hull nodes") + ax.set_xlabel(self._plane_cv_names[0], fontsize=14) + ax.set_ylabel(self._plane_cv_names[1], fontsize=14) + ax.set_title(title, fontsize=14) + ax.legend(fontsize=12) + + elif plot_3d: + # 3D surface plot. var_idx = self._Flamelet_Variables.index(var_to_plot) - fig = plt.figure(figsize=[10,10]) - ax = fig.add_subplot(111, projection='3d') + fig = plt.figure(figsize=(10, 10)) + ax = fig.add_subplot(111, projection='3d') ax.plot_trisurf(Nodes[:, p0], Nodes[:, p1], level_data[:, var_idx], - triangles=Tria, cmap='viridis', alpha=0.9, edgecolor='k', linewidth=0.2) - ax.set_xlabel(self._plane_cv_names[0], fontsize=20) - ax.set_ylabel(self._plane_cv_names[1], fontsize=20) - ax.set_title(self._level_cv_name + " = " + str(val_mix_frac)) - plt.show() + triangles=Tria, cmap='viridis', alpha=0.9, + edgecolor='k' if show_grid else 'none', + linewidth=0.2 if show_grid else 0) + ax.set_xlabel(self._plane_cv_names[0], fontsize=14) + ax.set_ylabel(self._plane_cv_names[1], fontsize=14) + ax.set_zlabel(var_to_plot, fontsize=14) + ax.set_title(title, fontsize=14) + + else: + # 2D colour-map plot. + import matplotlib.tri as mtri + var_idx = self._Flamelet_Variables.index(var_to_plot) + values = level_data[:, var_idx] + triang = mtri.Triangulation(Nodes[:, p0], Nodes[:, p1], Tria) + fig, ax = plt.subplots(figsize=(10, 7), constrained_layout=True) + tc = ax.tripcolor(triang, values, shading='gouraud', cmap='inferno') + if show_grid: + ax.triplot(triang, color='k', linewidth=0.2, alpha=0.4) + cb = fig.colorbar(tc, ax=ax, pad=0.02) + cb.set_label(var_to_plot, fontsize=12) + ax.set_xlabel(self._plane_cv_names[0], fontsize=14) + ax.set_ylabel(self._plane_cv_names[1], fontsize=14) + ax.set_title("%s — %s" % (var_to_plot, title), fontsize=14) + + if save_path: + fig.savefig(save_path, dpi=150, bbox_inches='tight') + print(" Saved: %s" % save_path) + if not show: + plt.close(fig) return def GenerateTableNodes(self): @@ -959,13 +1117,54 @@ def __GetPlaneCVBounds(self, val_level:float): plane1_min = self._Config.gas.enthalpy_mass return plane0_unb, plane0_b, plane1_min, plane1_max, plane1_at_unb + + elif (self._level_cv_name == DefaultSettings_FGM.name_pv and + self._plane_cv_names[0] == DefaultSettings_FGM.name_mixfrac and + self._plane_cv_names[1] == DefaultSettings_FGM.name_enth): + # For counterflow diffusion flames: 2D (Z, h) table at fixed PV + fuel = self._Config.GetFuelString() + ox = self._Config.GetOxidizerString() + T_lo, T_hi = self._Config.GetUnbTempBounds() + p = DefaultSettings_FGM.pressure + + # Plane CV 0 is mixture fraction: ranges from 0 (oxidizer) to 1 (fuel) + plane0_unb = 0.0 # Pure oxidizer + plane0_b = 1.0 # Pure fuel + + # Plane CV 1 is enthalpy: ranges from T_lo to T_hi + # At mixture fraction = 0 (oxidizer side) + self._Config.gas.set_mixture_fraction(0.0, fuel, ox) + self._Config.gas.TP = T_lo, p + h_ox_cold = self._Config.gas.enthalpy_mass + self._Config.gas.TP = T_hi, p + h_ox_hot = self._Config.gas.enthalpy_mass + + # At mixture fraction = 1 (fuel side) + self._Config.gas.set_mixture_fraction(1.0, fuel, ox) + self._Config.gas.TP = T_lo, p + h_fuel_cold = self._Config.gas.enthalpy_mass + plane1_at_unb = h_fuel_cold + self._Config.gas.TP = T_hi, p + h_fuel_hot = self._Config.gas.enthalpy_mass + + # Enthalpy range spans coldest to hottest across all mixture fractions + plane1_min = min(h_ox_cold, h_fuel_cold) + plane1_max = max(h_ox_hot, h_fuel_hot) + + return plane0_unb, plane0_b, plane1_min, plane1_max, plane1_at_unb + else: raise NotImplementedError( "Plane CV bounds not implemented for level='%s', plane=%s. " - "Currently only supported: level='%s', plane=['%s', '%s']." % ( + "Supported configurations:\n" + " (1) level='%s', plane=['%s', '%s'] for premixed flames\n" + " (2) level='%s', plane=['%s', '%s'] for non-premixed flames" % ( self._level_cv_name, self._plane_cv_names, DefaultSettings_FGM.name_mixfrac, DefaultSettings_FGM.name_pv, + DefaultSettings_FGM.name_enth, + DefaultSettings_FGM.name_pv, + DefaultSettings_FGM.name_mixfrac, DefaultSettings_FGM.name_enth)) def __ComputeCurvature(self, val_level:float): @@ -987,18 +1186,79 @@ def __ComputeCurvature(self, val_level:float): p_idxs = self._plane_cv_idxs n_cv = len(self._controlling_variables) - # Define 2D grid between minimum and maximum plane controlling variables. - plane0_range = np.linspace(plane0_unb, plane0_b, 800) - plane1_range = np.linspace(plane1_min, plane1_max, 800) + # COMPUTE CONVEX HULL FROM ACTUAL DATA POINTS (instead of rectangular grid) + # This preserves the natural parallelogram shape of the data instead of forcing it to a rectangle. + + # Extract the 2D plane coordinates from actual data in PHYSICAL units + data_2d = self._D_full[:, p_idxs] # Physical units: (N, 2) array [Z, h] + + # Compute convex hull in PHYSICAL space + try: + hull = ConvexHull(data_2d) + x_hull_phys = data_2d[hull.vertices, 0] + y_hull_phys = data_2d[hull.vertices, 1] + + # Normalize the hull vertices for gmsh + # For 2D tables, use the 2D scaler that only operates on (Z, h) + if self._is_2D_table and self._scaler_2d is not None: + hull_2d = np.column_stack([x_hull_phys, y_hull_phys]) + hull_norm = self._scaler_2d.transform(hull_2d) + x_hull = hull_norm[:, 0] + y_hull = hull_norm[:, 1] + else: + # For 3D tables, use the 3D scaler + hull_3d = np.column_stack([ + x_hull_phys, + y_hull_phys, + np.full(len(x_hull_phys), val_level) + ]) + x_hull, y_hull = self._scaler.transform(hull_3d)[:, p_idxs].T + + print(" Hull has %d vertices" % len(hull.vertices)) + except Exception as e: + print(" WARNING: ConvexHull failed on actual data: %s" % str(e)) + print(" Falling back to rectangular grid approach...") + + # Fallback: use rectangular grid approach + plane0_range = np.linspace(plane0_unb, plane0_b, 400) + plane1_range = np.linspace(plane1_min, plane1_max, 400) + xgrid, ygrid = np.meshgrid(plane0_range, plane1_range) + + n_pts = xgrid.size + CV_grid_init = np.zeros([n_pts, n_cv]) + CV_grid_init[:, p_idxs[0]] = xgrid.flatten() + CV_grid_init[:, p_idxs[1]] = ygrid.flatten() + CV_grid_init[:, self._level_cv_idx] = val_level + + # Burner-stabilized boundary filter + plane0_grid = CV_grid_init[:, p_idxs[0]] + plane1_grid = CV_grid_init[:, p_idxs[1]] + h_limit = ((plane1_at_unb - plane1_min) * plane0_grid + + (plane1_min * plane0_unb - plane1_at_unb * plane0_b)) / (plane0_unb - plane0_b) + idx_keep = plane1_grid >= h_limit + + CV_grid_norm_init = self._scaler.transform(CV_grid_init) + CV_grid_norm = CV_grid_norm_init[idx_keep, :] + + hull = ConvexHull(CV_grid_norm[:, p_idxs]) + x_hull = CV_grid_norm[hull.vertices, p_idxs[0]] + y_hull = CV_grid_norm[hull.vertices, p_idxs[1]] + + # 2: Generate refinement locations based on flamelet indicators + # Use coarser grid (400x400) to reduce memory + grid_resolution = 400 + plane0_range = np.linspace(plane0_unb, plane0_b, grid_resolution) + plane1_range = np.linspace(plane1_min, plane1_max, grid_resolution) xgrid, ygrid = np.meshgrid(plane0_range, plane1_range) n_pts = xgrid.size + print(" Computing refinement indicators on %dx%d grid (%d points)" % (grid_resolution, grid_resolution, n_pts)) CV_grid_init = np.zeros([n_pts, n_cv]) CV_grid_init[:, p_idxs[0]] = xgrid.flatten() CV_grid_init[:, p_idxs[1]] = ygrid.flatten() CV_grid_init[:, self._level_cv_idx] = val_level - # 2: Locate nodes that are above the burner-stabilized boundary line. + # Burner-stabilized boundary filter plane0_grid = CV_grid_init[:, p_idxs[0]] plane1_grid = CV_grid_init[:, p_idxs[1]] h_limit = ((plane1_at_unb - plane1_min) * plane0_grid + @@ -1010,12 +1270,7 @@ def __ComputeCurvature(self, val_level:float): CV_grid_norm_init = self._scaler.transform(CV_grid_init) CV_grid_norm = self._scaler.transform(CV_grid) - # 3: Generate convex hull on the plane coordinates. - hull = ConvexHull(CV_grid_norm[:, p_idxs]) - x_hull = CV_grid_norm[hull.vertices, p_idxs[0]] - y_hull = CV_grid_norm[hull.vertices, p_idxs[1]] - - # 4: Locate refinement locations based on the combined indicator across all refinement fields. + # Evaluate refinement fields to find high-gradient regions missing = [f for f in self._refinement_fields if f not in self._Flamelet_Variables] if missing: raise Exception("Refinement field(s) not found in flamelet data: %s. " @@ -1034,26 +1289,18 @@ def __ComputeCurvature(self, val_level:float): x_refinement = CV_grid_norm_init[idx_ref, p_idxs[0]] y_refinement = CV_grid_norm_init[idx_ref, p_idxs[1]] - # 5: Generate refinement locations at the unburnt and burnt plane0 boundaries. - plane1_unb_range = np.linspace(plane1_at_unb, plane1_max, self._Config.GetNpTemp()) - CV_unb = np.zeros([len(plane1_unb_range), n_cv]) - CV_unb[:, p_idxs[0]] = plane0_unb - CV_unb[:, p_idxs[1]] = plane1_unb_range - CV_unb[:, self._level_cv_idx] = val_level - CV_unb_norm = self._scaler.transform(CV_unb) - - plane1_b_range = np.linspace(plane1_min, plane1_max, self._Config.GetNpTemp()) - CV_b = np.zeros([len(plane1_b_range), n_cv]) - CV_b[:, p_idxs[0]] = plane0_b - CV_b[:, p_idxs[1]] = plane1_b_range - CV_b[:, self._level_cv_idx] = val_level - CV_b_norm = self._scaler.transform(CV_b) - - x_refinement = np.append(x_refinement, CV_unb_norm[:, p_idxs[0]]) - x_refinement = np.append(x_refinement, CV_b_norm[:, p_idxs[0]]) - y_refinement = np.append(y_refinement, CV_unb_norm[:, p_idxs[1]]) - y_refinement = np.append(y_refinement, CV_b_norm[:, p_idxs[1]]) + n_ref_initial = len(x_refinement) + print(" Found %d refinement points above threshold %.3f" % (n_ref_initial, self._curvature_threshold)) + + # Early subsampling if too many points (to reduce memory before adding boundary points) + N_max_initial = self._max_refinement_seeds * 2 # Allow 2x for boundary additions + if n_ref_initial > N_max_initial: + idx_sub = np.round(np.linspace(0, n_ref_initial - 1, N_max_initial)).astype(int) + x_refinement = x_refinement[idx_sub] + y_refinement = y_refinement[idx_sub] + print(" Early subsampled to %d points to reduce memory" % len(x_refinement)) + # 3: Generate refinement locations at the unburnt and burnt plane0 boundaries. XY_refinement = np.vstack((x_refinement, y_refinement)).T XY_hull = np.vstack((x_hull, y_hull)).T @@ -1101,6 +1348,11 @@ def __Compute2DMesh(self, XY_hull:np.ndarray, XY_refinement:np.ndarray, val_leve hull_cell_size = self._convex_hull_cell_size refinement_radius = self._refinement_radius #* np.sqrt(level_area) print("Generating 2D mesh with base cell size %.4f, hull cell size %.4f and refined cell size %.4f" % (base_cell_size, hull_cell_size, refined_cell_size)) + print(" Hull points: %d, Refinement seeds: %d" % (len(XY_hull), len(XY_refinement))) + + # Check for degenerate hull + if len(XY_hull) < 3: + raise Exception("Hull has fewer than 3 points - cannot create mesh.") hull_pts = [] for i in range(int(len(XY_hull)/2)): @@ -1110,18 +1362,23 @@ def __Compute2DMesh(self, XY_hull:np.ndarray, XY_refinement:np.ndarray, val_leve hull_pts_2.append(factory.addPoint(XY_hull[i, 0], XY_hull[i, 1], 0, hull_cell_size)) hull_pts_2.append(hull_pts[0]) + print(" Created %d points in hull_pts and %d points in hull_pts_2" % (len(hull_pts), len(hull_pts_2))) + # Subsample refinement seed points to avoid excessive PointsList size. N_max_seeds = self._max_refinement_seeds if len(XY_refinement) > N_max_seeds: idx_sub = np.round(np.linspace(0, len(XY_refinement) - 1, N_max_seeds)).astype(int) XY_refinement_sub = XY_refinement[idx_sub] + print(" Subsampled refinement seeds from %d to %d" % (len(XY_refinement), len(XY_refinement_sub))) else: XY_refinement_sub = XY_refinement + print(" Using all %d refinement seeds (below max %d)" % (len(XY_refinement), N_max_seeds)) embed_pts = [] for i in range(len(XY_refinement_sub)): pt_idx = factory.addPoint(XY_refinement_sub[i, 0], XY_refinement_sub[i, 1], 0, refined_cell_size) embed_pts.append(pt_idx) + print(" Created %d embedded refinement points in gmsh" % len(embed_pts)) hull_curve_1 = factory.addPolyline(hull_pts) hull_curve_2 = factory.addPolyline(hull_pts_2) @@ -1166,12 +1423,23 @@ def __Compute2DMesh(self, XY_hull:np.ndarray, XY_refinement:np.ndarray, val_leve gmsh.model.mesh.field.setAsBackgroundMesh(7) gmsh.option.setNumber("Mesh.Algorithm", 5) - gmsh.model.mesh.generate(2) + + try: + print(" Starting mesh generation...") + gmsh.model.mesh.generate(2) + print(" Mesh generation completed successfully") + except Exception as e: + print(" ERROR: Mesh generation failed!") + print(" Exception: %s" % str(e)) + # Skip gmsh cleanup to avoid additional issues + raise Exception("Gmsh meshing failed. Try increasing cell sizes or reducing refinement seeds.") + nodes = gmsh.model.mesh.getNodes(dim=2, tag=-1, includeBoundary=True, returnParametricCoord=False)[1] MeshPoints = np.array([nodes[::3], nodes[1::3]]).T + print(" Extracted %d mesh nodes" % len(MeshPoints)) - # we need finalize - gmsh.finalize() + # Don't finalize or clear - gmsh has issues with cleanup in some environments + # This causes a small memory leak but avoids segfaults # Build the full CV array in controlling-variable column order. plane0_norm = MeshPoints[:, 0] diff --git a/TestCases/LUT/CH4_fixed_strain/0_generate_config.py b/TestCases/LUT/CH4_fixed_strain/0_generate_config.py new file mode 100644 index 0000000..e80a0d1 --- /dev/null +++ b/TestCases/LUT/CH4_fixed_strain/0_generate_config.py @@ -0,0 +1,47 @@ +from Common.DataDrivenConfig import Config_FGM +import os + +Config = Config_FGM() +Config.SetConfigName("TableGeneration") + +# Fuel/oxidizer match equil.py: diluted CH4 (23% CH4 + 77% N2) vs diluted air +Config.SetFuelDefinition(fuel_species=["CH4", "N2"], fuel_weights=[0.23, 0.77]) +Config.SetOxidizerDefinition(oxidizer_species=["O2", "N2"], + oxidizer_weights=[0.23, 0.77]) +Config.SetReactionMechanism('gri30.yaml') + +# For counterflow diffusion flames, use mixture fraction bounds covering full range +Config.DefineMixtureStatus(True) # Use mixture fraction, not equivalence ratio +Config.SetMixtureBounds(0.0, 1.0) # Full mixture fraction range for non-premixed flames +Config.SetNpMix(2) # Minimum value to avoid empty range (not used for counterflow flames) + +Config.RunFreeFlames(False) +Config.SetUnbTempBounds(250, 800) +Config.SetNpTemp(56) + +Config.RunBurnerFlames(False) +Config.RunEquilibrium(False) + +Config.RunCounterFlames(True) +Config.SetCounterFlowFixedStrain(True) +Config.SetCounterFlowStrainRate(56.0) +Config.SetInitialGridLength(0.02) # flame width in metres +Config.SetSaveMoleFractions(True) # Save mole fractions (X-) alongside mass fractions (Y-) + +Config.SetTransportModel('unity-Lewis-number') +Config.SetConcatenationFileHeader("LUT_data") + + +# progress variable definition (still needed, not used) +Config.SetProgressVariableDefinition( + pv_species=['CO2', 'CO','H2','H2O'], + pv_weights=[1, 1, 1, 1]) + +# Preparing flamelet output directory. +flamelet_data_dir = os.getcwd() + "/flamelet_data/" +if not os.path.isdir(flamelet_data_dir): + os.mkdir(flamelet_data_dir) +Config.SetOutputDir(flamelet_data_dir) + +Config.PrintBanner() +Config.SaveConfig() diff --git a/TestCases/LUT/CH4_fixed_strain/1_generate_flamelet_data.py b/TestCases/LUT/CH4_fixed_strain/1_generate_flamelet_data.py new file mode 100644 index 0000000..27121dd --- /dev/null +++ b/TestCases/LUT/CH4_fixed_strain/1_generate_flamelet_data.py @@ -0,0 +1,19 @@ +# Generate flamelet data for a counterflow diffusion flame sweep + +# Limit inner thread pools BEFORE any library imports to prevent oversubscription +import os +os.environ.setdefault("OMP_NUM_THREADS", "1") +os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") +os.environ.setdefault("MKL_NUM_THREADS", "1") +os.environ.setdefault("VECLIB_MAXIMUM_THREADS", "1") +os.environ.setdefault("NUMEXPR_NUM_THREADS", "1") + +from Common.DataDrivenConfig import Config_FGM +from Data_Generation.DataGenerator_FGM import ComputeFlameletData + +# Load FGM configuration +Config = Config_FGM("TableGeneration.cfg") + +ComputeFlameletData(Config, run_parallel=False, N_processors=1, loglevel=1, + counter_flame_refine={"ratio": 3, "slope": 0.05, "curve": 0.05, "prune": 0.01}) + diff --git a/TestCases/LUT/CH4_fixed_strain/2_collect_flamelet_data.py b/TestCases/LUT/CH4_fixed_strain/2_collect_flamelet_data.py new file mode 100644 index 0000000..dfd7226 --- /dev/null +++ b/TestCases/LUT/CH4_fixed_strain/2_collect_flamelet_data.py @@ -0,0 +1,19 @@ +# Collect flamelet data into data sets for table generation +from Common.DataDrivenConfig import Config_FGM +from Data_Processing.collectFlameletData import FlameletConcatenator + +Config = Config_FGM("TableGeneration.cfg") + +Concat = FlameletConcatenator(Config) + +# Include NOx reaction rates and heat release in flamelet data set +Concat.SetAuxilarySpecies(["H2", "CO2", "H2O", "CO", "NOx"]) +Concat.SetLookUpVars(["Heat_Release", "Density", "Y-OH", "X-H2", "X-CO2", "X-H2O", "X-CO"]) + +# Apply source term and chemical equilibrium data corrections for table generation. +Concat.WriteLUTData(True) + +# No equilibrium data in this counterflow-only case + +# Read and concatenate flamelet data +Concat.ConcatenateFlameletData() diff --git a/TestCases/LUT/CH4_fixed_strain/3_generate_table.py b/TestCases/LUT/CH4_fixed_strain/3_generate_table.py new file mode 100644 index 0000000..3477528 --- /dev/null +++ b/TestCases/LUT/CH4_fixed_strain/3_generate_table.py @@ -0,0 +1,91 @@ +from Manifold_Generation.LUT.FlameletTableGeneration import SU2TableGenerator +from Common.DataDrivenConfig import Config_FGM +import matplotlib.pyplot as plt + +# ── Output options ────────────────────────────────────────────────────────── +SHOW_FIGURES = True # keep all figures open (non-blocking) while script runs +SAVE_FIGURES = True # save all figures as PNG files in the working directory + +if SHOW_FIGURES: + plt.ion() # interactive mode: plt.pause() displays without blocking + +def _save(name): + return name if SAVE_FIGURES else None + +def _show(fig_fn, *args, **kwargs): + """Call a visualisation method once with save_path and show set from the global flags.""" + if not SHOW_FIGURES and not SAVE_FIGURES: + return + kwargs['save_path'] = kwargs.get('save_path') if SAVE_FIGURES else None + kwargs['show'] = SHOW_FIGURES + fig_fn(*args, **kwargs) + +# Loading configuration. +Config = Config_FGM("TableGeneration.cfg") + +# Initializing table module and pre-process interpolator. +Tgen = SU2TableGenerator(Config, n_near=14, p_fac=3) + +# Generate a 2D (MixtureFraction, EnthalpyTot) LUT from counterflow diffusion +# flames at fixed strain rate. The plane axes are (Z, h); ProgressVariable is +# the nominal "level" variable but is unused (single-level 2D table). +Tgen.SetTableAxes(level_cv_name="ProgressVariable", + plane_cv_names=["MixtureFraction", "EnthalpyTot"]) + +# min == max triggers 2D mode: one mesh spanning the full (Z, h) data cloud. +Tgen.SetMixtureFractionLimits(mix_frac_min=0.0, mix_frac_max=0.0) +Tgen.SetNTableLevels(1) + +# Refinement: use Temperature and Heat_Release as indicators. +# Heat_Release concentrates points near the reaction zone in (Z, h) space. +Tgen.SetRefinementFields(["Heat_Release", "Temperature"]) +Tgen.SetRefinementMethod("gradient") + +# medium resolution (increased from original to avoid gmsh issues) +Tgen.SetBaseCellSize(1e-2) +Tgen.SetRefinedCellSize(5e-3) +Tgen.SetRefinementRadius(1e-2) +Tgen.SetMaxRefinementSeeds(200) +Tgen.SetHullCellSize(1e-2) + +# Generate table connectivity and interpolate flamelet data onto mesh nodes. +Tgen.GenerateTableNodes() + +# ── Diagnostics: inspect table quality before writing ────────────────────── +# T(Z) for 10 enthalpy levels spanning the full data range. +_show(Tgen.PlotTableSlices, + x_cv = "MixtureFraction", + y_var = "Temperature", + slice_cv = "EnthalpyTot", + slice_range = None, + n_slices = 10, + n_x_points = 300, + save_path = _save("T_vs_Z_slices.png")) + +# Heat_Release(Z) for the same 10 enthalpy levels. +_show(Tgen.PlotTableSlices, + x_cv = "MixtureFraction", + y_var = "Heat_Release", + slice_cv = "EnthalpyTot", + slice_range = None, + n_slices = 10, + n_x_points = 300, + save_path = _save("HRR_vs_Z_slices.png")) + +# Mesh only: +_show(Tgen.VisualizeTableLevel, 0.0, + save_path=_save("mesh_ZH.png")) +# 2D colour map of Temperature (no wireframe): +_show(Tgen.VisualizeTableLevel, 0.0, "Temperature", + plot_3d=False, show_grid=False, save_path=_save("T_ZH_2d.png")) +# 3D surface of Temperature (with wireframe): +_show(Tgen.VisualizeTableLevel, 0.0, "Temperature", + plot_3d=True, show_grid=True, save_path=_save("T_ZH_3d.png")) + +# Write SU2 .drg table file (Dragon v1.0.1 format, 2D, MixtureFraction + EnthalpyTot). +Tgen.WriteTableFile() + +# Keep all figure windows open until the user closes them. +if SHOW_FIGURES: + plt.ioff() + plt.show(block=True) diff --git a/TestCases/LUT/CH4_phi080_2D/0_generate_config.py b/TestCases/LUT/CH4_phi080_2D/0_generate_config.py new file mode 100644 index 0000000..565214b --- /dev/null +++ b/TestCases/LUT/CH4_phi080_2D/0_generate_config.py @@ -0,0 +1,47 @@ +from Common.DataDrivenConfig import Config_FGM +import os + +Config = Config_FGM() +Config.SetConfigName("TableGeneration") + +# Methane-air flamelets at a single equivalence ratio phi = 0.80 +Config.SetFuelDefinition(fuel_species=["CH4"], fuel_weights=[1.0]) +Config.SetReactionMechanism('gri30.yaml') + +# Single phi = 0.80 → single mixture fraction point (min == max) +Config.SetMixtureBounds(0.80, 0.80) +Config.SetNpMix(1) + +Config.SetUnbTempBounds(270, 750) +Config.SetNpTemp(49) +Config.SetNpMdot(40) # burner flames across the mdot range +Config.SetMdotDHTarget(20000.0) # J/kg target ΔH between flames +Config.SetNpMdotExtra(40) # synthetic flames linearly interpolated from lowest-mdot burner flame to equilibrium +Config.SetInitialGridLength(0.2) # Initial flamelet domain length in metres + +# Explicitly select which flamelet types to generate. +Config.RunFreeFlames(True) +Config.RunBurnerFlames(True) +Config.RunExtraInterpolatedBurnerFlames(True) +Config.SetSrcInterpExponent(2.0) # even faster decay + +Config.RunEquilibrium(True) + +Config.SetTransportModel('unity-Lewis-number') + +Config.SetConcatenationFileHeader("LUT_data") + + +# progress variable definition +Config.SetProgressVariableDefinition( + pv_species=['CO2', 'CO','H2','H2O'], + pv_weights=[1, 1, 1, 1]) + +# Preparing flamelet output directory. +flamelet_data_dir = os.getcwd() + "/flamelet_data/" +if not os.path.isdir(flamelet_data_dir): + os.mkdir(flamelet_data_dir) +Config.SetOutputDir(flamelet_data_dir) + +Config.PrintBanner() +Config.SaveConfig() diff --git a/TestCases/LUT/CH4_phi080_2D/1_generate_flamelet_data.py b/TestCases/LUT/CH4_phi080_2D/1_generate_flamelet_data.py new file mode 100644 index 0000000..9753541 --- /dev/null +++ b/TestCases/LUT/CH4_phi080_2D/1_generate_flamelet_data.py @@ -0,0 +1,27 @@ +# Generate flamelet data for a single phi = 0.80 methane-air flame + +# Limit inner thread pools BEFORE any library imports to prevent oversubscription +import os +os.environ.setdefault("OMP_NUM_THREADS", "1") +os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") +os.environ.setdefault("MKL_NUM_THREADS", "1") +os.environ.setdefault("VECLIB_MAXIMUM_THREADS", "1") +os.environ.setdefault("NUMEXPR_NUM_THREADS", "1") + +from Common.DataDrivenConfig import Config_FGM +from Data_Generation.DataGenerator_FGM import ComputeFlameletData, DataGenerator_Cantera + +# Load FGM configuration +Config = Config_FGM("TableGeneration.cfg") + +# refinement values: +# free_flame_refine={"ratio": 3.0, "slope": 0.1, "curve": 0.1, "prune":0.01}, +# this leads to Np = 180 + +if __name__ == '__main__': + ComputeFlameletData(Config, run_parallel=True, N_processors=4, loglevel=1, + # medium + #free_flame_refine={"ratio": 3, "slope": 0.02, "curve": 0.02, "prune": 0.005}, + free_flame_refine={"ratio": 3.0, "slope": 0.1, "curve": 0.1, "prune":0.01}, + burner_flame_refine={"ratio": 3.0, "slope": 0.15, "curve": 0.15, "prune": 0.01}) + diff --git a/TestCases/LUT/CH4_phi080_2D/2_collect_flamelet_data.py b/TestCases/LUT/CH4_phi080_2D/2_collect_flamelet_data.py new file mode 100644 index 0000000..d9f558a --- /dev/null +++ b/TestCases/LUT/CH4_phi080_2D/2_collect_flamelet_data.py @@ -0,0 +1,20 @@ +# Collect flamelet data into data sets for table generation +from Common.DataDrivenConfig import Config_FGM +from Data_Processing.collectFlameletData import FlameletConcatenator + +Config = Config_FGM("TableGeneration.cfg") + +Concat = FlameletConcatenator(Config) + +# Include NOx reaction rates and heat release in flamelet data set +Concat.SetAuxilarySpecies(["H2", "CO2", "H2O", "CO", "NOx"]) +Concat.SetLookUpVars(["Heat_Release", "Density", "Y-OH"]) + +# Apply source term and chemical equilibrium data corrections for table generation. +Concat.WriteLUTData(True) + +Concat.SetNEquilibriumNodes(1) # sample N rows from each eq file + + +# Read and concatenate flamelet data +Concat.ConcatenateFlameletData() diff --git a/TestCases/LUT/CH4_phi080_2D/3_generate_table.py b/TestCases/LUT/CH4_phi080_2D/3_generate_table.py new file mode 100644 index 0000000..9c13273 --- /dev/null +++ b/TestCases/LUT/CH4_phi080_2D/3_generate_table.py @@ -0,0 +1,71 @@ +from Manifold_Generation.LUT.FlameletTableGeneration import SU2TableGenerator +from Common.DataDrivenConfig import Config_FGM +import matplotlib.pyplot as plt + +# Loading configuration. +Config = Config_FGM("TableGeneration.cfg") + +# Initializing table module and pre-process interpolator. +Tgen = SU2TableGenerator(Config, n_near=14, p_fac=1) + +# Generate a 2D (ProgressVariable, EnthalpyTot) LUT for a single equivalence +# ratio phi = 0.80. Because phi_min == phi_max, the table generator writes a +# Dragon v1.0.1 file without mixture-fraction levels. +Tgen.SetEquivalenceRatioLimits(phi_min=0.80, phi_max=0.80) +Tgen.SetNTableLevels(1) +Tgen.SetRefinementFields(["ProdRateTot_PV", "Y_dot_net-CO", "Y_dot_pos-CO","Y_dot_neg-CO"]) + +# small +#Tgen.SetBaseCellSize(1e-2) +#Tgen.SetRefinedCellSize(1e-2) +#Tgen.SetRefinementRadius(1e-2) +#Tgen.SetRefinementMethod("gradient") +#Tgen.SetMaxRefinementSeeds(500) +#Tgen.SetHullCellSize(1.0e-2) + +# medium +Tgen.SetBaseCellSize(5e-3) +Tgen.SetRefinedCellSize(5e-3) +Tgen.SetRefinementRadius(5e-3) +Tgen.SetRefinementMethod("gradient") +Tgen.SetMaxRefinementSeeds(500) +Tgen.SetHullCellSize(5.0e-3) + +# fine +#Tgen.SetBaseCellSize(2.5e-3) +#Tgen.SetRefinedCellSize(1.0e-3) +#Tgen.SetRefinementRadius(2.5e-3) +#Tgen.SetRefinementMethod("gradient") +#Tgen.SetMaxRefinementSeeds(500) +#Tgen.SetHullCellSize(2.5e-3) + + + + +Tgen.SetTableAxes(level_cv_name="MixtureFraction", + plane_cv_names=["ProgressVariable", "EnthalpyTot"]) + +# Visualize the table mesh and reaction rate at phi = 0.80. +cv_target = Config.GetUnburntScalars(equivalence_ratio=0.80, temperature=300.0) +pv_target = cv_target[0] +z_target = cv_target[2] +print("Target unburnt progress variable:", pv_target) +Tgen.VisualizeTableLevel(z_target, "ProdRateTot_PV", save_path="prodrate.png", show=True, show_grid=False) +Tgen.VisualizeTableLevel(z_target, "Temperature", save_path="temperature.png",show=True, show_grid=False) +Tgen.VisualizeTableLevel(z_target, "Y_dot_net-CO", save_path="co_net.png",show=True, show_grid=False) +Tgen.VisualizeTableLevel(z_target, "Y_dot_pos-CO", save_path="co_pos.png",show=True, show_grid=False) +Tgen.VisualizeTableLevel(z_target, "Y_dot_neg-CO", save_path="co_neg.png",show=True, show_grid=False) +Tgen.VisualizeTableLevel(z_target, "Y_dot_net-NOx", save_path="nox_net.png",show=True, show_grid=False) +Tgen.VisualizeTableLevel(z_target, show=True) + +# Show all plots +plt.show() + +# Generate table connectivity and interpolate flamelet data. +Tgen.GenerateTableNodes() +# from 99% of the max progress variable, set the source terms of CO and H2 to zero +# if the absolute value of the source terms is |S| < 0.1 +Tgen.ClampSourceTerms(species_list=["CO", "H2", "CO2", "H2O"], pv_frac=0.99, abs_tol=0.1) + +# Write SU2 .drg table file (Dragon v1.0.1 format, 2D). +Tgen.WriteTableFile() diff --git a/TestCases/LUT/CH4_phi080_2D/visualize_hull_vs_data.py b/TestCases/LUT/CH4_phi080_2D/visualize_hull_vs_data.py new file mode 100755 index 0000000..3190a2b --- /dev/null +++ b/TestCases/LUT/CH4_phi080_2D/visualize_hull_vs_data.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +""" +Diagnostic script to visualize actual flamelet data points vs the constructed mesh hull. +This helps identify whether the convex hull is cutting off concave data regions. + +Usage: + python visualize_hull_vs_data.py + +Example: + python visualize_hull_vs_data.py TestCases/LUT/Flamelet_Table_CH4_phi080/TableGeneration.cfg 0.055 +""" + +import sys +from Manifold_Generation.LUT.FlameletTableGeneration import SU2TableGenerator +from Common.DataDrivenConfig import Config_FGM + +def main(): + if len(sys.argv) < 3: + print(__doc__) + print("\nUsing default values:") + config_path = "TableGeneration.cfg" + z_value = 0.055 + print(f" Config: {config_path}") + print(f" Mixture fraction: {z_value}") + else: + config_path = sys.argv[1] + z_value = float(sys.argv[2]) + + print("\n" + "="*70) + print("Hull vs Data Visualization") + print("="*70) + print(f"Loading configuration: {config_path}") + + # Load configuration + Config = Config_FGM(config_path) + + # Initialize table generator + print("Initializing SU2TableGenerator...") + Tgen = SU2TableGenerator(Config, n_near=14, p_fac=1) + + # Configure table settings (match your 3_generate_table.py settings) + Tgen.SetEquivalenceRatioLimits(phi_min=0.80, phi_max=0.80) + Tgen.SetNTableLevels(1) + Tgen.SetRefinementFields(["ProdRateTot_PV", "Y_dot_net-CO", "Y_dot_pos-CO","Y_dot_neg-CO"]) + + # Set mesh resolution + Tgen.SetBaseCellSize(5e-3) + Tgen.SetRefinedCellSize(5e-3) + Tgen.SetRefinementRadius(5e-3) + Tgen.SetRefinementMethod("gradient") + Tgen.SetMaxRefinementSeeds(500) + Tgen.SetHullCellSize(5.0e-3) + + # Set table axes + Tgen.SetTableAxes(level_cv_name="MixtureFraction", + plane_cv_names=["ProgressVariable", "EnthalpyTot"]) + + print(f"\nGenerating table level visualization for Z = {z_value}...") + print("This will create two plots:") + print(" 1. Actual data points with convex hull overlay") + print(" 2. Mesh nodes with actual data overlay") + print("\nLook for data points that fall OUTSIDE the red hull boundary.") + print("These are being excluded from the table!\n") + + # Visualize with hull diagnostic enabled + Tgen.VisualizeTableLevel(z_value, visualize_hull=True, show=True) + + print("\n" + "="*70) + print("Visualization complete!") + print("="*70) + +if __name__ == "__main__": + main()