diff --git a/applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp b/applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp index 943cc72ef9c0..60e36ef86e54 100644 --- a/applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp +++ b/applications/RomApplication/custom_python/add_custom_utilities_to_python.cpp @@ -54,6 +54,8 @@ void AddCustomUtilitiesToPython(pybind11::module& m) .def_static("GetElementIdsNotInHRomModelPart", &RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart) .def_static("GetHRomMinimumConditionsIds", &RomAuxiliaryUtilities::GetHRomMinimumConditionsIds) .def_static("ProjectRomSolutionIncrementToNodes", &RomAuxiliaryUtilities::ProjectRomSolutionIncrementToNodes) + .def_static("GetElementIdsInModelPart", &RomAuxiliaryUtilities::GetElementIdsInModelPart) + .def_static("GetConditionIdsInModelPart", &RomAuxiliaryUtilities::GetConditionIdsInModelPart) ; } diff --git a/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp b/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp index 4a1064f62209..917bd6a43597 100644 --- a/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp +++ b/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.cpp @@ -342,7 +342,7 @@ std::vector RomAuxiliaryUtilities::GetNodalNeighbouringElementIdsNotI { std::vector new_element_ids; const auto& r_elem_weights = rHRomWeights.at("Elements"); - + FindGlobalNodalEntityNeighboursProcess find_nodal_elements_neighbours_process(rModelPart); find_nodal_elements_neighbours_process.Execute(); @@ -373,7 +373,7 @@ std::vector RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart( for (const auto& r_elem : rModelPartWithElementsToInclude.Elements()) { IndexType element_id = r_elem.Id(); - + // Check if the element is already added if (r_elem_weights.find(element_id - 1) == r_elem_weights.end()) { new_element_ids.push_back(element_id - 1); @@ -383,6 +383,7 @@ std::vector RomAuxiliaryUtilities::GetElementIdsNotInHRomModelPart( return new_element_ids; } + std::vector RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart( const ModelPart& rModelPart, const ModelPart& rModelPartWithConditionsToInclude, @@ -393,7 +394,7 @@ std::vector RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart( for (const auto& r_cond : rModelPartWithConditionsToInclude.Conditions()) { IndexType condition_id = r_cond.Id(); - + // Check if the condition is already added if (r_cond_weights.find(condition_id - 1) == r_cond_weights.end()) { new_condition_ids.push_back(condition_id - 1); @@ -403,6 +404,28 @@ std::vector RomAuxiliaryUtilities::GetConditionIdsNotInHRomModelPart( return new_condition_ids; } +std::vector RomAuxiliaryUtilities::GetElementIdsInModelPart( + const ModelPart& rModelPart) +{ + std::vector element_ids; + + for (const auto& r_elem : rModelPart.Elements()) { + element_ids.push_back(r_elem.Id()); + } + return element_ids; +} + +std::vector RomAuxiliaryUtilities::GetConditionIdsInModelPart( + const ModelPart& rModelPart) +{ + std::vector condition_ids; + + for (const auto& r_cond : rModelPart.Conditions()) { + condition_ids.push_back(r_cond.Id()); + } + return condition_ids; +} + std::vector RomAuxiliaryUtilities::GetHRomMinimumConditionsIds( const ModelPart& rModelPart, const std::map& rHRomConditionWeights) @@ -567,6 +590,6 @@ void RomAuxiliaryUtilities::GetPsiElemental( noalias(row(rPsiElemental, i)) = row(r_nodal_rom_basis, row_id); } } - } + } } // namespace Kratos diff --git a/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h b/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h index 5bd8364cdfd0..0e52a38caa49 100644 --- a/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h +++ b/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h @@ -107,14 +107,14 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities static std::vector GetHRomConditionParentsIds( const ModelPart& rModelPart, const std::map>& rHRomWeights); - + /** * @brief Retrieve the IDs of elements neighboring nodes in a given sub-model part but not present in HRom weights. * - * This function iterates over all the nodes in the provided sub-model part (`rGivenModelPart`) and collects - * the IDs of the elements neighboring each node. The neighboring elements are determined using the - * 'NEIGHBOUR_ELEMENTS' value attached to each node. The function then checks if these elements are already - * present in `rHRomWeights`. If not, their IDs are added to the return vector. It's important to note that + * This function iterates over all the nodes in the provided sub-model part (`rGivenModelPart`) and collects + * the IDs of the elements neighboring each node. The neighboring elements are determined using the + * 'NEIGHBOUR_ELEMENTS' value attached to each node. The function then checks if these elements are already + * present in `rHRomWeights`. If not, their IDs are added to the return vector. It's important to note that * this function assumes that the 'NEIGHBOUR_ELEMENTS' values are already computed for the nodes in the model part. * * @param rModelPart The complete model part which houses all the elements. @@ -168,6 +168,15 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities const ModelPart& rModelPart, const std::map& rHRomConditionWeights); + + static std::vector GetElementIdsInModelPart( + const ModelPart& rModelPart + ); + + static std::vector GetConditionIdsInModelPart( + const ModelPart& rModelPart + ); + /** * @brief Project the ROM solution increment onto the nodal basis * For a given model part this function takes the ROM_SOLUTION_INCREMENT, which is assumed to be diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 62639521aa48..b1573e424ea2 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -1,5 +1,5 @@ import numpy as np - +from KratosMultiphysics import Logger try: from matplotlib import pyplot as plt @@ -14,31 +14,42 @@ class EmpiricalCubatureMethod(): Reference: Hernandez 2020. "A multiscale method for periodic structures using domain decomposition and ECM-hyperreduction" """ - - """ - Constructor setting up the parameters for the Element Selection Strategy - ECM_tolerance: approximation tolerance for the element selection algorithm - Filter_tolerance: parameter limiting the number of candidate points (elements) to those above this tolerance - Plotting: whether to plot the error evolution of the element selection algorithm - """ - def __init__(self, ECM_tolerance = 0, Filter_tolerance = 1e-16, Plotting = False): + def __init__( + self, + ECM_tolerance = 0, + Filter_tolerance = 1e-16, + Plotting = False, + MaximumNumberUnsuccesfulIterations = 100 + ): + """ + Constructor setting up the parameters for the Element Selection Strategy + ECM_tolerance: approximation tolerance for the element selection algorithm + Filter_tolerance: parameter limiting the number of candidate points (elements) to those above this tolerance + Plotting: whether to plot the error evolution of the element selection algorithm + """ self.ECM_tolerance = ECM_tolerance self.Filter_tolerance = Filter_tolerance self.Name = "EmpiricalCubature" self.Plotting = Plotting - - - - """ - Method for setting up the element selection - input: - ResidualsBasis: numpy array containing a basis to the residuals projected - - constrain_sum_of_weights: enable the user to constrain weights to be the sum of the number of entities. - - constrain_conditions: enable the user to enforce weights to consider conditions (for specific boundary conditions). - """ - def SetUp(self, ResidualsBasis, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = 0): - + self.MaximumNumberUnsuccesfulIterations = MaximumNumberUnsuccesfulIterations + + def SetUp( + self, + ResidualsBasis, + InitialCandidatesSet = None, + constrain_sum_of_weights=True, + constrain_conditions = False, + number_of_conditions = 0 + ): + """ + Method for setting up the element selection + input: - ResidualsBasis: numpy array containing a basis to the residuals projected + - constrain_sum_of_weights: enable the user to constrain weights to be the sum of the number of entities. + - constrain_conditions: enable the user to enforce weights to consider conditions (for specific boundary conditions). + """ self.W = np.ones(np.shape(ResidualsBasis)[0]) self.G = ResidualsBasis.T + self.y = InitialCandidatesSet self.add_constrain_count = None total_number_of_entities = np.shape(self.G)[1] elements_constraint = np.ones(total_number_of_entities) @@ -66,47 +77,89 @@ def SetUp(self, ResidualsBasis, constrain_sum_of_weights=True, constrain_conditi self.G = np.vstack([ self.G , projection_of_constant_vector_conditions ] ) self.add_constrain_count = -2 self.b = self.G @ self.W + self.UnsuccesfulIterations = 0 - - """ - Method performing calculations required before launching the Calculate method - """ def Initialize(self): - self.Gnorm = np.sqrt(sum(np.multiply(self.G, self.G), 0)) + """ + Method performing calculations required before launching the Calculate method + """ + self.GnormNOONE = np.linalg.norm(self.G[:self.add_constrain_count,:], axis = 0) M = np.shape(self.G)[1] normB = np.linalg.norm(self.b) - self.y = np.arange(0,M,1) # Set of candidate points (those whose associated column has low norm are removed) - GnormNOONE = np.sqrt(sum(np.multiply(self.G[:self.add_constrain_count,:], self.G[:self.add_constrain_count,:]), 0)) - if self.Filter_tolerance > 0: - TOL_REMOVE = self.Filter_tolerance * normB - rmvpin = np.where(GnormNOONE[self.y] < TOL_REMOVE) - self.y = np.delete(self.y,rmvpin) + + if self.y is None: + self.y = np.arange(0,M,1) # Set of candidate points (those whose associated column has low norm are removed) + + if self.Filter_tolerance > 0: + TOL_REMOVE = self.Filter_tolerance * normB + rmvpin = np.where(self.GnormNOONE[self.y] < TOL_REMOVE) + #self.y_complement = self.y[rmvpin] + self.y = np.delete(self.y,rmvpin) + else: + self.y_complement = np.arange(0, M, 1) # Initialize complement with all points + self.y_complement = np.delete(self.y_complement, self.y) # Remove candidates from complement + + if self.Filter_tolerance > 0: + TOL_REMOVE = self.Filter_tolerance * normB # Compute removal tolerance + + # Filter out low-norm columns from complement + rmvpin_complement = np.where(self.GnormNOONE[self.y_complement] < TOL_REMOVE) + self.y_complement = np.delete(self.y_complement, rmvpin_complement) + + # Filter out low-norm columns from candidates + rmvpin = np.where(self.GnormNOONE[self.y] < TOL_REMOVE) + removed_count = np.size(rmvpin) + self.y = np.delete(self.y, rmvpin) + + # Warning if some candidates were removed + if removed_count > 0: + Logger.PrintWarning("EmpiricalCubatureMethod", f"Some of the candidates were removed ({removed_count} removed). To include all candidates (with 0 weights in the HROM model part) for visualization and projection, consider using 'include_elements_model_parts_list' and 'include_conditions_model_parts_list' in the 'hrom_settings'.") + + # Warning if all candidates were removed + if np.size(self.y) == 0: + Logger.PrintWarning("EmpiricalCubatureMethod", "All candidates were removed because they have no contribution to the residual. To include them all (with 0 weights in the HROM model part) for visualization and projection, use 'include_elements_model_parts_list' and 'include_conditions_model_parts_list' in the 'hrom_settings'.") + self.y = self.y_complement # Set candidates to complement + self.z = {} # Set of intergration points self.mPOS = 0 # Number of nonzero weights - self.r = self.b # residual vector + self.r = self.b.copy() # residual vector self.m = len(self.b) # Default number of points self.nerror = np.linalg.norm(self.r)/normB self.nerrorACTUAL = self.nerror - def Run(self): self.Initialize() self.Calculate() + def expand_candidates_with_complement(self): + self.y = np.r_[self.y,self.y_complement] + print('expanding set to include the complement...') + ExpandedSetFlag = True + return ExpandedSetFlag - """ - Method launching the element selection algorithm to find a set of elements: self.z, and wiegths: self.w - """ def Calculate(self): - + """ + Method launching the element selection algorithm to find a set of elements: self.z, and wiegths: self.w + """ + MaximumLengthZ = 0 + ExpandedSetFlag = False k = 1 # number of iterations - while self.nerrorACTUAL > self.ECM_tolerance and self.mPOS < self.m and len(self.y) != 0: + self.success = True + while self.nerrorACTUAL > self.ECM_tolerance and self.mPOS < self.m and np.size(self.y) != 0: + + if self.UnsuccesfulIterations > self.MaximumNumberUnsuccesfulIterations and not ExpandedSetFlag: + ExpandedSetFlag = self.expand_candidates_with_complement() #Step 1. Compute new point - ObjFun = self.G[:,self.y].T @ self.r.T - ObjFun = ObjFun.T / self.Gnorm[self.y] - indSORT = np.argmax(ObjFun) - i = self.y[indSORT] + if np.size(self.y)==1: + #candidate set consists of a single element + indSORT = 0 + i = int(self.y) + else: + ObjFun = self.G[:,self.y].T @ self.r.T + ObjFun = ObjFun.T / self.GnormNOONE[self.y] + indSORT = np.argmax(ObjFun) + i = self.y[indSORT] if k==1: alpha = np.linalg.lstsq(self.G[:, [i]], self.b)[0] H = 1/(self.G[:,i] @ self.G[:,i].T) @@ -118,7 +171,13 @@ def Calculate(self): self.z = i else: self.z = np.r_[self.z,i] - self.y = np.delete(self.y,indSORT) + + #self.y = np.delete(self.y,indSORT) + if np.size(self.y)==1: + self.expand_candidates_with_complement() + self.y = np.delete(self.y,indSORT) + else: + self.y = np.delete(self.y,indSORT) # Step 4. Find possible negative weights if any(alpha < 0): @@ -130,9 +189,15 @@ def Calculate(self): alpha = H @ (self.G[:, self.z].T @ self.b) alpha = alpha.reshape(len(alpha),1) + if np.size(self.z) > MaximumLengthZ : + self.UnsuccesfulIterations = 0 + else: + self.UnsuccesfulIterations += 1 + #Step 6 Update the residual - if len(alpha)==1: - self.r = self.b - (self.G[:,self.z] * alpha) + if np.size(alpha)==1: + self.r = self.b.reshape(-1,1) - (self.G[:,self.z] * alpha).reshape(-1,1) + self.r = np.squeeze(self.r) else: Aux = self.G[:,self.z] @ alpha self.r = np.squeeze(self.b - Aux.T) @@ -150,8 +215,18 @@ def Calculate(self): ERROR_GLO = np.c_[ ERROR_GLO , self.nerrorACTUAL] NPOINTS = np.c_[ NPOINTS , np.size(self.z)] + MaximumLengthZ = max(MaximumLengthZ, np.size(self.z)) k = k+1 + if k>1000: + """ + this means using the initial candidate set, it was impossible to obtain a set of positive weights. + Try again without constraints!!! + TODO: incorporate this into greater workflow + """ + self.success = False + break + self.w = alpha.T * np.sqrt(self.W[self.z]) #TODO FIXME cope with weights vectors different from 1 print(f'Total number of iterations = {k}') @@ -163,17 +238,16 @@ def Calculate(self): plt.ylabel('Error %') plt.show() - - - """ - Method for the quick update of weights (self.w), whenever a negative weight is found - """ def _UpdateWeightsInverse(self, A,Aast,a,xold): + """ + Method for the quick update of weights (self.w), whenever a negative weight is found + """ c = np.dot(A.T, a) d = np.dot(Aast, c).reshape(-1, 1) s = np.dot(a.T, a) - np.dot(c.T, d) aux1 = np.hstack([Aast + np.outer(d, d) / s, -d / s]) if np.shape(-d.T / s)[1]==1: + s = s.reshape(1,-1) aux2 = np.squeeze(np.hstack([-d.T / s, 1 / s])) else: aux2 = np.hstack([np.squeeze(-d.T / s), 1 / s]) @@ -182,24 +256,20 @@ def _UpdateWeightsInverse(self, A,Aast,a,xold): x = np.vstack([(xold - d * v), v]) return Bast, x - - - """ - Method for the quick update of weights (self.w), whenever a negative weight is found - """ def _MultiUpdateInverseHermitian(self, invH, neg_indexes): + """ + Method for the quick update of weights (self.w), whenever a negative weight is found + """ neg_indexes = np.sort(neg_indexes) for i in range(np.size(neg_indexes)): neg_index = neg_indexes[i] - i invH = self._UpdateInverseHermitian(invH, neg_index) return invH - - - """ - Method for the quick update of weights (self.w), whenever a negative weight is found - """ def _UpdateInverseHermitian(self, invH, neg_index): + """ + Method for the quick update of weights (self.w), whenever a negative weight is found + """ if neg_index == np.shape(invH)[1]: aux = (invH[0:-1, -1] * invH[-1, 0:-1]) / invH(-1, -1) invH_new = invH[:-1, :-1] - aux diff --git a/applications/RomApplication/python_scripts/hrom_training_utility.py b/applications/RomApplication/python_scripts/hrom_training_utility.py index 6afab61db463..9b0a00a5ec1e 100644 --- a/applications/RomApplication/python_scripts/hrom_training_utility.py +++ b/applications/RomApplication/python_scripts/hrom_training_utility.py @@ -61,6 +61,32 @@ def __init__(self, solver, custom_settings): if not self.solver.model.HasModelPart(model_part_name): raise Exception('The model part named "' + model_part_name + '" does not exist in the model') + + # getting initial candidate ids for empirical cubature + initial_candidate_elements_model_part_list = settings["initial_candidate_elements_model_part_list"].GetStringArray() + initial_candidate_conditions_model_part_list = settings["initial_candidate_conditions_model_part_list"].GetStringArray() + + candidate_ids = np.empty(0) + for model_part_name in initial_candidate_elements_model_part_list: + if not self.solver.model.HasModelPart(model_part_name): + raise Exception('The model part named "' + model_part_name + '" does not exist in the model') + this_modelpart_element_ids = KratosROM.RomAuxiliaryUtilities.GetElementIdsInModelPart(self.solver.model.GetModelPart(model_part_name)) + if len(this_modelpart_element_ids)>0: + candidate_ids = np.r_[candidate_ids, np.array(this_modelpart_element_ids)] + + number_of_elements = self.solver.GetComputingModelPart().GetRootModelPart().NumberOfElements() + for model_part_name in initial_candidate_conditions_model_part_list: + if not self.solver.model.HasModelPart(model_part_name): + raise Exception('The model part named "' + model_part_name + '" does not exist in the model') + this_modelpart_condition_ids = KratosROM.RomAuxiliaryUtilities.GetConditionIdsInModelPart(self.solver.model.GetModelPart(model_part_name)) + if len(this_modelpart_condition_ids)>0: + candidate_ids = np.r_[candidate_ids, np.array(this_modelpart_condition_ids)+number_of_elements] + + if np.size(candidate_ids)>0: + self.candidate_ids = np.unique(candidate_ids).astype(int) - 1 # this -1 takes into account the id difference in numpy and Kratos + else: + self.candidate_ids = None + # Rom settings files self.rom_basis_output_name = Path(custom_settings["rom_basis_output_name"].GetString()) self.rom_basis_output_folder = Path(custom_settings["rom_basis_output_folder"].GetString()) @@ -97,9 +123,13 @@ def CalculateAndSaveHRomWeights(self): # Calculate the residuals basis and compute the HROM weights from it residual_basis = self.__CalculateResidualBasis() n_conditions = self.solver.GetComputingModelPart().NumberOfConditions() # Conditions must be included as an extra restriction to enforce ECM to capture all BC's regions. - self.hyper_reduction_element_selector.SetUp(residual_basis, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = n_conditions) + self.hyper_reduction_element_selector.SetUp(residual_basis, InitialCandidatesSet = self.candidate_ids, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = n_conditions) self.hyper_reduction_element_selector.Run() - + if not self.hyper_reduction_element_selector.success: + KratosMultiphysics.Logger.PrintWarning("HRomTrainingUtility", "The Empirical Cubature Method did not converge using the initial set of candidates. Launching again without initial candidates.") + #Imposing an initial candidate set can lead to no convergence. Restart without imposing the initial candidate set + self.hyper_reduction_element_selector.SetUp(residual_basis, InitialCandidatesSet = None, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = n_conditions) + self.hyper_reduction_element_selector.Run() # Save the HROM weights in the RomParameters.json # Note that in here we are assuming this naming convention for the ROM json file self.AppendHRomWeightsToRomParameters() @@ -172,9 +202,11 @@ def __GetHRomTrainingDefaultSettings(cls): "projection_strategy": "galerkin", "include_conditions_model_parts_list": [], "include_elements_model_parts_list": [], + "initial_candidate_elements_model_part_list" : [], + "initial_candidate_conditions_model_part_list" : [], "include_nodal_neighbouring_elements_model_parts_list":[], "include_minimum_condition": false, - "include_condition_parents": true + "include_condition_parents": false }""") return default_settings @@ -250,7 +282,7 @@ def AppendHRomWeightsToRomParameters(self): # If needed, update your weights and indexes using __AddSelectedElementsWithZeroWeights function with the new_elements weights, indexes = self.__AddSelectedElementsWithZeroWeights(weights, indexes, new_elements) - + # Add nodal neighbouring elements for model_part_name in self.include_nodal_neighbouring_elements_model_parts_list: # Check if the sub model part exists @@ -286,6 +318,8 @@ def AppendHRomWeightsToRomParameters(self): # If required, add the HROM conditions parent elements # Note that we add these with zero weight so their future assembly will have no effect if self.include_condition_parents: + KratosMultiphysics.Logger.PrintWarning("HRomTrainingUtility", 'Make sure you set "assign_neighbour_elements_to_conditions": true in the solver_settings to have a parent element for each condition.') + # Get the HROM condition parents from the current HROM weights missing_condition_parents = KratosROM.RomAuxiliaryUtilities.GetHRomConditionParentsIds( self.solver.GetComputingModelPart().GetRootModelPart(), #TODO: I think this one should be the root diff --git a/applications/RomApplication/python_scripts/rom_manager.py b/applications/RomApplication/python_scripts/rom_manager.py index b254e81976b9..6698289b181c 100644 --- a/applications/RomApplication/python_scripts/rom_manager.py +++ b/applications/RomApplication/python_scripts/rom_manager.py @@ -309,8 +309,13 @@ def __LaunchTrainHROM(self, mu_train, store_residuals_projected=False): self._StoreSnapshotsMatrix('residuals_projected',RedidualsSnapshotsMatrix) u,_,_,_ = RandomizedSingularValueDecomposition(COMPUTE_V=False).Calculate(RedidualsSnapshotsMatrix, self.hrom_training_parameters["element_selection_svd_truncation_tolerance"].GetDouble()) - simulation.GetHROM_utility().hyper_reduction_element_selector.SetUp(u) + simulation.GetHROM_utility().hyper_reduction_element_selector.SetUp(u, InitialCandidatesSet = simulation.GetHROM_utility().candidate_ids) simulation.GetHROM_utility().hyper_reduction_element_selector.Run() + if not simulation.GetHROM_utility().hyper_reduction_element_selector.success: + KratosMultiphysics.Logger.PrintWarning("HRomTrainingUtility", "The Empirical Cubature Method did not converge using the initial set of candidates. Launching again without initial candidates.") + #Imposing an initial candidate set can lead to no convergence. Restart without imposing the initial candidate set + self.hyper_reduction_element_selector.SetUp(u, InitialCandidatesSet = None) + self.hyper_reduction_element_selector.Run() simulation.GetHROM_utility().AppendHRomWeightsToRomParameters() simulation.GetHROM_utility().CreateHRomModelParts() @@ -482,6 +487,16 @@ def __LaunchRunHROM(self, mu_run, use_full_model_part): simulation.Run() + def _AddHromParametersToRomParameters(self,f): + f["rom_settings"]['rom_bns_settings']['monotonicity_preserving'] = self.general_rom_manager_parameters["ROM"]["galerkin_rom_bns_settings"]["monotonicity_preserving"].GetBool() if self.general_rom_manager_parameters["ROM"]["lspg_rom_bns_settings"].Has("monotonicity_preserving") else False + f["hrom_settings"]["element_selection_type"] = self.general_rom_manager_parameters["HROM"]["element_selection_type"].GetString() + f["hrom_settings"]["element_selection_svd_truncation_tolerance"] = self.general_rom_manager_parameters["HROM"]["element_selection_svd_truncation_tolerance"].GetDouble() + f["hrom_settings"]["create_hrom_visualization_model_part"] = self.general_rom_manager_parameters["HROM"]["create_hrom_visualization_model_part"].GetBool() + f["hrom_settings"]["echo_level"] = self.general_rom_manager_parameters["HROM"]["echo_level"].GetInt() + f["hrom_settings"]["include_condition_parents"] = self.general_rom_manager_parameters["HROM"]["include_condition_parents"].GetBool() if self.general_rom_manager_parameters["HROM"].Has("include_condition_parents") else False + f["hrom_settings"]["initial_candidate_elements_model_part_list"] = self.general_rom_manager_parameters["HROM"]["initial_candidate_elements_model_part_list"].GetStringArray() if self.general_rom_manager_parameters["HROM"].Has("initial_candidate_elements_model_part_list") else [] + f["hrom_settings"]["initial_candidate_conditions_model_part_list"] = self.general_rom_manager_parameters["HROM"]["initial_candidate_conditions_model_part_list"].GetStringArray() if self.general_rom_manager_parameters["HROM"].Has("initial_candidate_conditions_model_part_list") else [] + def _ChangeRomFlags(self, simulation_to_run = 'ROM'): """ This method updates the Flags present in the RomParameters.json file @@ -500,6 +515,7 @@ def _ChangeRomFlags(self, simulation_to_run = 'ROM'): with parameters_file_path.open('r+') as parameter_file: f=json.load(parameter_file) f['assembling_strategy'] = self.general_rom_manager_parameters['assembling_strategy'].GetString() if self.general_rom_manager_parameters.Has('assembling_strategy') else 'global' + self._AddHromParametersToRomParameters(f) if simulation_to_run=='GalerkinROM': f['projection_strategy']="galerkin" f['train_hrom']=False @@ -663,10 +679,19 @@ def _GetDefaulRomBasisOutputParameters(self): def _GetDefaulHromTrainingParameters(self): hrom_training_parameters = KratosMultiphysics.Parameters("""{ + "hrom_format": "numpy", "element_selection_type": "empirical_cubature", - "element_selection_svd_truncation_tolerance": 0, + "element_selection_svd_truncation_tolerance": 1.0e-6, "echo_level" : 0, - "create_hrom_visualization_model_part" : true + "create_hrom_visualization_model_part" : true, + "projection_strategy": "galerkin", + "include_conditions_model_parts_list": [], + "include_elements_model_parts_list": [], + "initial_candidate_elements_model_part_list" : [], + "initial_candidate_conditions_model_part_list" : [], + "include_nodal_neighbouring_elements_model_parts_list":[], + "include_minimum_condition": false, + "include_condition_parents": false }""") return hrom_training_parameters diff --git a/applications/RomApplication/tests/test_empirical_cubature_method.py b/applications/RomApplication/tests/test_empirical_cubature_method.py index 69ce1e0538fd..1cb7cd0bf5f7 100644 --- a/applications/RomApplication/tests/test_empirical_cubature_method.py +++ b/applications/RomApplication/tests/test_empirical_cubature_method.py @@ -27,11 +27,12 @@ def test_empirical_cubature_method(self): TestMatrix = synthetic_matrix(degree) #Get a synthetic matrix (During the training of a ROM model, this is a matrix of residuals projected onto a basis) TestMatrixBasis = calculate_basis(TestMatrix) ElementSelector = EmpiricalCubatureMethod() - ElementSelector.SetUp(TestMatrixBasis, constrain_sum_of_weights=False) - ElementSelector.Initialize() - ElementSelector.Calculate() - - self.assertEqual( len(ElementSelector.z) , degree + 1) #for a polynomial of degree n, n+1 points are to be selected + ElementSelector.SetUp(TestMatrixBasis, InitialCandidatesSet = np.array([0]), constrain_sum_of_weights=False) + ElementSelector.Run() + if ElementSelector.success is not True: + ElementSelector.SetUp(TestMatrixBasis, InitialCandidatesSet = None, constrain_sum_of_weights=False) + ElementSelector.Run() + self.assertEqual( np.size(ElementSelector.z) , degree + 1) #for a polynomial of degree n, n+1 points are to be selected # Cleaning kratos_utilities.DeleteDirectoryIfExisting("__pycache__")