From 49c7fcfa00afd9395ca7bc3d1c96068de7b51a33 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Mon, 18 Sep 2023 10:43:11 +0200 Subject: [PATCH 01/21] Placing docstrings in the correct place --- .../empirical_cubature_method.py | 81 +++++++++---------- 1 file changed, 38 insertions(+), 43 deletions(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 62639521aa48..da3f8143860d 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -14,29 +14,34 @@ 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): + """ + 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): - + def SetUp( + self, + ResidualsBasis, + 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.add_constrain_count = None @@ -67,11 +72,10 @@ def SetUp(self, ResidualsBasis, constrain_sum_of_weights=True, constrain_conditi self.add_constrain_count = -2 self.b = self.G @ self.W - - """ - Method performing calculations required before launching the Calculate method - """ def Initialize(self): + """ + Method performing calculations required before launching the Calculate method + """ self.Gnorm = np.sqrt(sum(np.multiply(self.G, self.G), 0)) M = np.shape(self.G)[1] normB = np.linalg.norm(self.b) @@ -88,17 +92,14 @@ def Initialize(self): self.nerror = np.linalg.norm(self.r)/normB self.nerrorACTUAL = self.nerror - def Run(self): self.Initialize() self.Calculate() - - """ - 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 + """ k = 1 # number of iterations while self.nerrorACTUAL > self.ECM_tolerance and self.mPOS < self.m and len(self.y) != 0: @@ -163,12 +164,10 @@ 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) @@ -182,24 +181,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 From 8027324e50d84f566ea6ea762e25fecf6026c37e Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Mon, 18 Sep 2023 10:46:44 +0200 Subject: [PATCH 02/21] adding new flag --- .../python_scripts/empirical_cubature_method.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index da3f8143860d..5b6bd3fab581 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -18,7 +18,9 @@ def __init__( self, ECM_tolerance = 0, Filter_tolerance = 1e-16, - Plotting = False): + Plotting = False, + MaximumNumberUnsuccesfulIterations = 100 + ): """ Constructor setting up the parameters for the Element Selection Strategy ECM_tolerance: approximation tolerance for the element selection algorithm @@ -35,7 +37,8 @@ def SetUp( ResidualsBasis, constrain_sum_of_weights=True, constrain_conditions = False, - number_of_conditions = 0): + number_of_conditions = 0 + ): """ Method for setting up the element selection input: - ResidualsBasis: numpy array containing a basis to the residuals projected From 31984bd742758cba64a4fbc775a40ba904569e6d Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Mon, 18 Sep 2023 10:52:04 +0200 Subject: [PATCH 03/21] adding initial candidate set option in ECM --- .../RomApplication/python_scripts/empirical_cubature_method.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 5b6bd3fab581..189fc32de6f6 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -31,10 +31,12 @@ def __init__( self.Filter_tolerance = Filter_tolerance self.Name = "EmpiricalCubature" self.Plotting = Plotting + self.MaximumNumberUnsuccesfulIterations = MaximumNumberUnsuccesfulIterations def SetUp( self, ResidualsBasis, + InitialCandidatesSet = None, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = 0 @@ -47,6 +49,7 @@ def SetUp( """ 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) From 65a9747025935374859ee2ef7c020d0c49ca902c Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Mon, 18 Sep 2023 18:41:38 +0200 Subject: [PATCH 04/21] Modifications on ecm to accept initial candadate set --- .../empirical_cubature_method.py | 67 +++++++++++++++---- 1 file changed, 53 insertions(+), 14 deletions(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 189fc32de6f6..f2ca761804c9 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -82,18 +82,29 @@ def Initialize(self): """ Method performing calculations required before launching the Calculate method """ - self.Gnorm = np.sqrt(sum(np.multiply(self.G, self.G), 0)) + 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) + self.y_complement = np.delete(self.y_complement, self.y)# 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_complement] < TOL_REMOVE) + self.y_complement = np.delete(self.y_complement,rmvpin) + 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 @@ -102,18 +113,31 @@ def Run(self): self.Initialize() self.Calculate() + def expand_candidates_with_complement(self): + self.y = np.union1d(self.y, self.y_complement) + print('expanding set to include the complement...') + ExpandedSetFlag = True + return ExpandedSetFlag + 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: #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 isinstance(self.y, np.int64) or isinstance(self.y, np.int32): + #candidate set consists of a single element + indSORT = 0 + i = self.y + else: + 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 k==1: alpha = np.linalg.lstsq(self.G[:, [i]], self.b)[0] H = 1/(self.G[:,i] @ self.G[:,i].T) @@ -137,9 +161,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) @@ -159,6 +189,15 @@ def Calculate(self): 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}') From 975574c123a209dc815dc216d4d9eda05093c3e8 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Mon, 18 Sep 2023 18:44:41 +0200 Subject: [PATCH 05/21] bug --- .../RomApplication/python_scripts/empirical_cubature_method.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index f2ca761804c9..7906a7a2ce7f 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -135,7 +135,7 @@ def Calculate(self): i = self.y else: ObjFun = self.G[:,self.y].T @ self.r.T - ObjFun = ObjFun.T / self.Gnorm[self.y] + ObjFun = ObjFun.T / self.GnormNOONE[self.y] indSORT = np.argmax(ObjFun) i = self.y[indSORT] if k==1: From c2d8ce9168efa6cf004484a36af49f39d0efbd8a Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 11:40:05 +0200 Subject: [PATCH 06/21] Adding to cpp functions to retrive elements and conditions ids --- .../rom_auxiliary_utilities.cpp | 31 ++++++++++++++++--- .../rom_auxiliary_utilities.h | 19 +++++++++--- 2 files changed, 41 insertions(+), 9 deletions(-) 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..c4a1f89ba1d1 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); + + std::vector GetElementIdsInModelPart( + const ModelPart& rModelPart + ); + + 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 From 1900ef74cc4362bbdaac5b2c054c7caf9c088a72 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 11:49:39 +0200 Subject: [PATCH 07/21] Exporting functions to Python --- .../custom_python/add_custom_utilities_to_python.cpp | 2 ++ .../RomApplication/custom_utilities/rom_auxiliary_utilities.h | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) 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.h b/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h index c4a1f89ba1d1..0e52a38caa49 100644 --- a/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h +++ b/applications/RomApplication/custom_utilities/rom_auxiliary_utilities.h @@ -169,11 +169,11 @@ class KRATOS_API(ROM_APPLICATION) RomAuxiliaryUtilities const std::map& rHRomConditionWeights); - std::vector GetElementIdsInModelPart( + static std::vector GetElementIdsInModelPart( const ModelPart& rModelPart ); - std::vector GetConditionIdsInModelPart( + static std::vector GetConditionIdsInModelPart( const ModelPart& rModelPart ); From 25da04c0733743e27444565078318a1fb12b5876 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 14:15:13 +0200 Subject: [PATCH 08/21] ECM with candidates --- .../python_scripts/empirical_cubature_method.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 7906a7a2ce7f..da2bea4edeaf 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -77,6 +77,7 @@ def SetUp( 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 def Initialize(self): """ @@ -128,6 +129,9 @@ def Calculate(self): k = 1 # number of iterations while self.nerrorACTUAL > self.ECM_tolerance and self.mPOS < self.m and len(self.y) != 0: + if self.UnsuccesfulIterations > self.MaximumNumberUnsuccesfulIterations and not ExpandedSetFlag: + ExpandedSetFlag = self.expand_candidates_with_complement() + #Step 1. Compute new point if isinstance(self.y, np.int64) or isinstance(self.y, np.int32): #candidate set consists of a single element From 667e83bb2e122b144542d15382a66d67d4db4aac Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 14:15:49 +0200 Subject: [PATCH 09/21] Adding initial candidate sets to ECM in hrom_training_utility --- .../python_scripts/hrom_training_utility.py | 34 +++++++++++++++++-- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/applications/RomApplication/python_scripts/hrom_training_utility.py b/applications/RomApplication/python_scripts/hrom_training_utility.py index 6afab61db463..61a7c8dce94a 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) - 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,7 +123,7 @@ 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() # Save the HROM weights in the RomParameters.json @@ -172,9 +198,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 +278,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 From 43c61b5560429d4df23b3c0bc33a5d027ee7d8ff Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 14:24:17 +0200 Subject: [PATCH 10/21] more Additions to ECM for candidate elements --- .../python_scripts/empirical_cubature_method.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index da2bea4edeaf..88df2f865ed7 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -127,7 +127,8 @@ def Calculate(self): 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() @@ -153,7 +154,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 isinstance(self.y, np.int64) or isinstance(self.y, np.int32): + 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): @@ -191,6 +198,7 @@ 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: @@ -222,6 +230,7 @@ def _UpdateWeightsInverse(self, A,Aast,a,xold): 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]) From c828703232938ca2feb32ef31edab23fa74ea404 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 14:32:51 +0200 Subject: [PATCH 11/21] contemplating the case where the ECM did not converge due to initial set --- .../RomApplication/python_scripts/hrom_training_utility.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/hrom_training_utility.py b/applications/RomApplication/python_scripts/hrom_training_utility.py index 61a7c8dce94a..0715a8d45d26 100644 --- a/applications/RomApplication/python_scripts/hrom_training_utility.py +++ b/applications/RomApplication/python_scripts/hrom_training_utility.py @@ -125,7 +125,10 @@ def CalculateAndSaveHRomWeights(self): 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, 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: + #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() From 4d0c6d49b490adfe391032721ef7c9084f673a55 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 14:49:08 +0200 Subject: [PATCH 12/21] integers required --- .../RomApplication/python_scripts/hrom_training_utility.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/hrom_training_utility.py b/applications/RomApplication/python_scripts/hrom_training_utility.py index 0715a8d45d26..c218b9b72142 100644 --- a/applications/RomApplication/python_scripts/hrom_training_utility.py +++ b/applications/RomApplication/python_scripts/hrom_training_utility.py @@ -83,7 +83,7 @@ def __init__(self, solver, custom_settings): 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) - 1 # this -1 takes into account the id difference in numpy and Kratos + 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 From b30eb296ec0181e8eb185743f904a6dccf1f25d5 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 15:22:27 +0200 Subject: [PATCH 13/21] RomManager changes to accept initial candidate set --- .../python_scripts/rom_manager.py | 30 +++++++++++++++++-- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/applications/RomApplication/python_scripts/rom_manager.py b/applications/RomApplication/python_scripts/rom_manager.py index b254e81976b9..42dd31162e4e 100644 --- a/applications/RomApplication/python_scripts/rom_manager.py +++ b/applications/RomApplication/python_scripts/rom_manager.py @@ -309,8 +309,12 @@ 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: + #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, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = n_conditions) + self.hyper_reduction_element_selector.Run() simulation.GetHROM_utility().AppendHRomWeightsToRomParameters() simulation.GetHROM_utility().CreateHRomModelParts() @@ -482,6 +486,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() + f["hrom_settings"]["initial_candidate_elements_model_part_list"] = self.general_rom_manager_parameters["HROM"]["initial_candidate_elements_model_part_list"].GetStringArray() + f["hrom_settings"]["initial_candidate_conditions_model_part_list"] = self.general_rom_manager_parameters["HROM"]["initial_candidate_conditions_model_part_list"].GetStringArray() + def _ChangeRomFlags(self, simulation_to_run = 'ROM'): """ This method updates the Flags present in the RomParameters.json file @@ -500,6 +514,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 +678,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 From 696030934d2c5c339e9eff90851101210ec480de Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 16:27:42 +0200 Subject: [PATCH 14/21] bug fix --- .../python_scripts/empirical_cubature_method.py | 6 +++--- applications/RomApplication/python_scripts/rom_manager.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 88df2f865ed7..c29000a1b943 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -134,10 +134,10 @@ def Calculate(self): ExpandedSetFlag = self.expand_candidates_with_complement() #Step 1. Compute new point - if isinstance(self.y, np.int64) or isinstance(self.y, np.int32): + if np.size(self.y)==1: #candidate set consists of a single element indSORT = 0 - i = self.y + i = int(self.y) else: ObjFun = self.G[:,self.y].T @ self.r.T ObjFun = ObjFun.T / self.GnormNOONE[self.y] @@ -156,7 +156,7 @@ def Calculate(self): self.z = np.r_[self.z,i] #self.y = np.delete(self.y,indSORT) - if isinstance(self.y, np.int64) or isinstance(self.y, np.int32): + if np.size(self.y)==1: self.expand_candidates_with_complement() self.y = np.delete(self.y,indSORT) else: diff --git a/applications/RomApplication/python_scripts/rom_manager.py b/applications/RomApplication/python_scripts/rom_manager.py index 42dd31162e4e..e854759c7292 100644 --- a/applications/RomApplication/python_scripts/rom_manager.py +++ b/applications/RomApplication/python_scripts/rom_manager.py @@ -313,7 +313,7 @@ def __LaunchTrainHROM(self, mu_train, store_residuals_projected=False): simulation.GetHROM_utility().hyper_reduction_element_selector.Run() if not simulation.GetHROM_utility().hyper_reduction_element_selector.success: #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, constrain_sum_of_weights=True, constrain_conditions = False, number_of_conditions = n_conditions) + self.hyper_reduction_element_selector.SetUp(u, InitialCandidatesSet = None) self.hyper_reduction_element_selector.Run() simulation.GetHROM_utility().AppendHRomWeightsToRomParameters() simulation.GetHROM_utility().CreateHRomModelParts() From 220b7a04f84f1b2997b183cff549a2edce5ca9e6 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 16:46:46 +0200 Subject: [PATCH 15/21] bug fix in ECM expansion of candidate set --- .../RomApplication/python_scripts/empirical_cubature_method.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index c29000a1b943..b2a9c50e08e8 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -115,7 +115,7 @@ def Run(self): self.Calculate() def expand_candidates_with_complement(self): - self.y = np.union1d(self.y, self.y_complement) + self.y = np.r_[self.y,self.y_complement] print('expanding set to include the complement...') ExpandedSetFlag = True return ExpandedSetFlag From c51ad8536d5f04319837cc904474f8d368a70e7b Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Tue, 19 Sep 2023 16:50:42 +0200 Subject: [PATCH 16/21] updating ECM test --- .../tests/test_empirical_cubature_method.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) 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__") From 52dec05e25c2f63c8a7c1f9a5a688e9990a233da Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 21 Sep 2023 12:38:16 +0200 Subject: [PATCH 17/21] Adding warning. --- .../RomApplication/python_scripts/hrom_training_utility.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/applications/RomApplication/python_scripts/hrom_training_utility.py b/applications/RomApplication/python_scripts/hrom_training_utility.py index c218b9b72142..d112f3ec2234 100644 --- a/applications/RomApplication/python_scripts/hrom_training_utility.py +++ b/applications/RomApplication/python_scripts/hrom_training_utility.py @@ -317,6 +317,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 From 9ab2a944181fd2b3b830e6b634335d1e4fae30de Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 21 Sep 2023 12:38:47 +0200 Subject: [PATCH 18/21] Setting defaults. --- applications/RomApplication/python_scripts/rom_manager.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/RomApplication/python_scripts/rom_manager.py b/applications/RomApplication/python_scripts/rom_manager.py index e854759c7292..af9eeb8596e9 100644 --- a/applications/RomApplication/python_scripts/rom_manager.py +++ b/applications/RomApplication/python_scripts/rom_manager.py @@ -492,9 +492,9 @@ def _AddHromParametersToRomParameters(self,f): 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() - f["hrom_settings"]["initial_candidate_elements_model_part_list"] = self.general_rom_manager_parameters["HROM"]["initial_candidate_elements_model_part_list"].GetStringArray() - f["hrom_settings"]["initial_candidate_conditions_model_part_list"] = self.general_rom_manager_parameters["HROM"]["initial_candidate_conditions_model_part_list"].GetStringArray() + 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'): """ From 06f3b56a744b5291b48e8a3c33d577c71227ed91 Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 21 Sep 2023 12:39:41 +0200 Subject: [PATCH 19/21] Adding warnings in case some candidates entities do not contribute to R. --- .../empirical_cubature_method.py | 27 +++++++++++++++---- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index b2a9c50e08e8..5fa7ac474d79 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -96,12 +96,29 @@ def Initialize(self): #self.y_complement = self.y[rmvpin] self.y = np.delete(self.y,rmvpin) else: - self.y_complement = np.arange(0,M,1) - self.y_complement = np.delete(self.y_complement, self.y)# Set of candidate points (those whose associated column has low norm are removed) + 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 - rmvpin = np.where(self.GnormNOONE[self.y_complement] < TOL_REMOVE) - self.y_complement = np.delete(self.y_complement,rmvpin) + 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 From f6bb99849dc70bfa3178ae8dbbf5a426fe454e01 Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 21 Sep 2023 15:46:46 +0200 Subject: [PATCH 20/21] Adding warning and importing logger. --- .../RomApplication/python_scripts/empirical_cubature_method.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 5fa7ac474d79..42397913035f 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 @@ -148,6 +148,7 @@ def Calculate(self): 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: + Logger.PrintWarning("EmpiricalCubatureMethod", "The Empirical Cubature Method did not converge using the initial set of candidates.") ExpandedSetFlag = self.expand_candidates_with_complement() #Step 1. Compute new point From 057ce8a17dce1076aff337ac81be9a2f48672b35 Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 21 Sep 2023 16:09:14 +0200 Subject: [PATCH 21/21] Updating warnings. --- .../RomApplication/python_scripts/empirical_cubature_method.py | 1 - .../RomApplication/python_scripts/hrom_training_utility.py | 1 + applications/RomApplication/python_scripts/rom_manager.py | 1 + 3 files changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/empirical_cubature_method.py b/applications/RomApplication/python_scripts/empirical_cubature_method.py index 42397913035f..b1573e424ea2 100644 --- a/applications/RomApplication/python_scripts/empirical_cubature_method.py +++ b/applications/RomApplication/python_scripts/empirical_cubature_method.py @@ -148,7 +148,6 @@ def Calculate(self): 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: - Logger.PrintWarning("EmpiricalCubatureMethod", "The Empirical Cubature Method did not converge using the initial set of candidates.") ExpandedSetFlag = self.expand_candidates_with_complement() #Step 1. Compute new point diff --git a/applications/RomApplication/python_scripts/hrom_training_utility.py b/applications/RomApplication/python_scripts/hrom_training_utility.py index d112f3ec2234..9b0a00a5ec1e 100644 --- a/applications/RomApplication/python_scripts/hrom_training_utility.py +++ b/applications/RomApplication/python_scripts/hrom_training_utility.py @@ -126,6 +126,7 @@ def CalculateAndSaveHRomWeights(self): 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() diff --git a/applications/RomApplication/python_scripts/rom_manager.py b/applications/RomApplication/python_scripts/rom_manager.py index af9eeb8596e9..6698289b181c 100644 --- a/applications/RomApplication/python_scripts/rom_manager.py +++ b/applications/RomApplication/python_scripts/rom_manager.py @@ -312,6 +312,7 @@ def __LaunchTrainHROM(self, mu_train, store_residuals_projected=False): 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()