Merge pull request #12140 from KratosMultiphysics/RomApp_FixBugIndexes

[RomApp] Fixing Ids issues in ROM files.

applications/RomApplication/python_scripts/calculate_rom_basis_output_process.py

@@ -183,8 +183,12 @@ def _PrintRomBasis(self, snapshots_matrix):
 
         elif self.rom_basis_output_format == "numpy":
             # Storing modes in Numpy format
+            node_ids = []
+            for node in self.model_part.Nodes:
+                node_ids.append(node.Id)
+            node_ids = numpy.array(node_ids)
             numpy.save(self.rom_basis_output_folder / "RightBasisMatrix.npy", u)
-            numpy.save(self.rom_basis_output_folder / "NodeIds.npy", numpy.arange(1,((u.shape[0]+1)/n_nodal_unknowns), 1, dtype=int))
+            numpy.save(self.rom_basis_output_folder / "NodeIds.npy", node_ids)
         else:
             err_msg = "Unsupported output format {}.".format(self.rom_basis_output_format)
             raise Exception(err_msg)

applications/RomApplication/python_scripts/hrom_training_utility.py

@@ -64,10 +64,13 @@ def __init__(self, solver, custom_settings):
                 raise Exception('The model part named "' + model_part_name + '" does not exist in the model')
 
 
+        #TODO use cpp mappings
+        root_model_part = self.solver.GetComputingModelPart().GetRootModelPart()
+        number_of_elements = root_model_part.NumberOfElements()
+        self._setup_mappings(root_model_part, number_of_elements)
+
         # 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):
@@ -80,16 +83,19 @@ def __init__(self, solver, custom_settings):
             else:
                 this_modelpart_element_ids = KratosROM.RomAuxiliaryUtilities.GetElementIdsInModelPart(candidate_elements_model_part)
             if len(this_modelpart_element_ids)>0:
-                candidate_ids = np.r_[candidate_ids, np.array(this_modelpart_element_ids)]
+                this_modelpart_element_ids = self.map_element_ids_to_numpy_indexes(this_modelpart_element_ids)
+                candidate_ids = np.r_[candidate_ids, this_modelpart_element_ids]
+
+
+        initial_candidate_conditions_model_part_list = settings["initial_candidate_conditions_model_part_list"].GetStringArray()
 
-        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]
-
+                this_modelpart_condition_ids = self.map_condition_ids_to_numpy_indexes(this_modelpart_condition_ids)
+                candidate_ids = np.r_[candidate_ids, this_modelpart_condition_ids]
         if np.size(candidate_ids)>0:
             self.candidate_ids = np.unique(candidate_ids).astype(int)
         else:
@@ -99,6 +105,47 @@ def __init__(self, solver, custom_settings):
         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())
 
+    def _setup_mappings(self, root_model_part, number_of_elements):
+        self.element_id_to_numpy_index_mapping = {}
+        self.numpy_index_to_element_id_mapping = {}
+        for index, element in enumerate(root_model_part.Elements):
+            self.numpy_index_to_element_id_mapping[index] = element.Id-1 #FIXME -1
+            self.element_id_to_numpy_index_mapping[element.Id-1] = index #FIXME -1
+
+        self.condition_id_to_numpy_index_mapping = {}
+        self.numpy_index_to_condition_id_mapping = {}
+        for index, condition in enumerate(root_model_part.Conditions):
+            self.numpy_index_to_condition_id_mapping[index+number_of_elements] = condition.Id -1 +number_of_elements  #FIXME -1 #FIXME +number_of_elements We should remove redundant fixes
+            self.condition_id_to_numpy_index_mapping[condition.Id-1] = index+number_of_elements #FIXME -1
+
+
+    def map_condition_ids_to_numpy_indexes(self, this_modelpart_condition_ids):
+        this_modelpart_indexes_numpy = []
+        for cond_id in this_modelpart_condition_ids:
+            this_modelpart_indexes_numpy.append(self.condition_id_to_numpy_index_mapping[cond_id])
+        return np.array(this_modelpart_indexes_numpy)
+
+
+    def map_element_ids_to_numpy_indexes(self, this_modelpart_element_ids):
+        this_modelpart_indexes_numpy = []
+        for elem_id in this_modelpart_element_ids:
+            this_modelpart_indexes_numpy.append(self.condition_id_to_numpy_index_mapping[elem_id])
+        return np.array(this_modelpart_indexes_numpy)
+
+
+    def map_numpy_indexes_to_element_and_conditions_ids(self,indexes,number_of_elements):
+
+        kratos_indexes = []
+        for i in range(np.size(indexes)):
+            if indexes[i]<=number_of_elements-1:
+                kratos_indexes.append(self.numpy_index_to_element_id_mapping[indexes[i]])
+            else:
+                kratos_indexes.append(self.numpy_index_to_condition_id_mapping[indexes[i]])
+
+        return np.array(kratos_indexes) #FIXME -1
+
+
+
     def AppendCurrentStepResiduals(self):
         # Get the computing model part from the solver implementing the problem physics
         computing_model_part = self.solver.GetComputingModelPart()
@@ -268,6 +315,7 @@ def AppendHRomWeightsToRomParameters(self):
         number_of_elements = self.solver.GetComputingModelPart().GetRootModelPart().NumberOfElements()
         weights = np.squeeze(self.hyper_reduction_element_selector.w)
         indexes = self.hyper_reduction_element_selector.z
+        indexes = self.map_numpy_indexes_to_element_and_conditions_ids(indexes,number_of_elements)
 
         # Create dictionary with HROM weights (Only used for the expansion of the selected Conditions to include their parent Elements)
         hrom_weights = self.__CreateDictionaryWithRomElementsAndWeights(weights,indexes,number_of_elements)

applications/RomApplication/tests/calculate_rom_basis_output_process_test_files/RomParameters_test_Results.json

@@ -1,42 +1,52 @@
 {
+    "rom_manager": false,
+    "train_hrom": false,
+    "run_hrom": false,
+    "projection_strategy": "galerkin",
+    "assembling_strategy": "global",
+    "rom_format": "json",
     "rom_settings": {
+        "rom_bns_settings": {},
         "nodal_unknowns": [
             "TEMPERATURE"
         ],
         "number_of_rom_dofs": 2,
-        "petrov_galerkin_number_of_rom_dofs": 0,
-        "rom_bns_settings":{}
+        "petrov_galerkin_number_of_rom_dofs": 0
+    },
+    "hrom_settings": {
+        "hrom_format": "json"
     },
     "nodal_modes": {
-        "1": [
+        "2": [
             [
-                -0.17574638455134867,
-                0.7543959227866553
+                -0.01983405679273213,
+                0.5993762115792658
             ]
         ],
-        "2": [
+        "5": [
             [
-                -0.2949770806015417,
-                0.4615067950960123
+                -0.04391666176273041,
+                0.5763832650962621
             ]
         ],
-        "3": [
+        "12": [
             [
-                -0.4142077766517345,
-                0.16861766740536666
+                -0.10010940669272668,
+                0.5227330566359412
             ]
         ],
-        "4": [
+        "69": [
             [
-                -0.5334384727019273,
-                -0.12427146028527797
+                -0.5576789011226962,
+                0.08586707345904437
             ]
         ],
-        "5": [
+        "102": [
             [
-                -0.6526691687521202,
-                -0.4171605879759215
+                -0.8225875557926786,
+                -0.16705533785390972
             ]
         ]
-    }
+    },
+    "elements_and_weights": {}
 }

applications/RomApplication/tests/test_RomApplication.py

@@ -10,7 +10,8 @@
 from test_structural_rom import TestStructuralRom
 from test_randomized_singular_value_decomposition import TestRandomizedSVD
 from test_empirical_cubature_method import TestEmpiricalCubatureMethod
-from test_calculate_rom_basis_output_process import TestCalculateRomBasisOutputProcess
+from test_calculate_rom_basis_output_process_json import TestCalculateRomBasisOutputProcessJSON
+from test_calculate_rom_basis_output_process_numpy import TestCalculateRomBasisOutputProcessNumpy
 from test_compressible_potiential_rom import TestCompressiblePotentialRom
 from test_fluid_lspg_rom import TestFluidLSPGRom
 from test_thermal_lspg_rom import TestThermalLSPGRom
@@ -42,7 +43,8 @@ def AssembleTestSuites():
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestFluidRom]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestThermalRom]))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestStructuralRom]))
-    smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestCalculateRomBasisOutputProcess]))
+    smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestCalculateRomBasisOutputProcessJSON]))
+    smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestCalculateRomBasisOutputProcessNumpy]))
     smallSuite.addTest(TestRandomizedSVD('test_radomized_svd'))
     smallSuite.addTest(TestEmpiricalCubatureMethod('test_empirical_cubature_method'))
     smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestCompressiblePotentialRom]))

applications/RomApplication/tests/test_calculate_rom_basis_output_process.py → applications/RomApplication/tests/test_calculate_rom_basis_output_process_json.py

@@ -7,7 +7,7 @@
 from KratosMultiphysics.compare_two_files_check_process import CompareTwoFilesCheckProcess
 from KratosMultiphysics.RomApplication.calculate_rom_basis_output_process import CalculateRomBasisOutputProcess
 
-class TestCalculateRomBasisOutputProcess(KratosUnittest.TestCase):
+class TestCalculateRomBasisOutputProcessJSON(KratosUnittest.TestCase):
 
     def setUp(self):
         # Test data
@@ -21,9 +21,10 @@ def setUp(self):
         self.model = KratosMultiphysics.Model()
         model_part = self.model.CreateModelPart("MainModelPart")
         model_part.AddNodalSolutionStepVariable(KratosMultiphysics.TEMPERATURE)
-        n_nodes = 5
+        node_ids_non_ordered = [2,5,12,69,102]
+        n_nodes = len(node_ids_non_ordered)
         for i in range(n_nodes):
-            model_part.CreateNewNode(i+1,float(i),0.0,0.0)
+            model_part.CreateNewNode(node_ids_non_ordered[i],float(i),0.0,0.0)
 
     def testCalculateRomBasisOutputProcess(self):
         # Create a CalculateROMBasisOutputProcess instance

applications/RomApplication/tests/test_calculate_rom_basis_output_process_numpy.py

@@ -0,0 +1,124 @@
+import os
+
+import KratosMultiphysics
+import KratosMultiphysics.KratosUnittest as KratosUnittest
+import KratosMultiphysics.kratos_utilities as kratos_utilities
+from KratosMultiphysics.RomApplication.calculate_rom_basis_output_process import CalculateRomBasisOutputProcess
+import numpy as np
+
+class TestCalculateRomBasisOutputProcessNumpy(KratosUnittest.TestCase):
+
+    def setUp(self):
+        # Test data
+        self.print_output = False
+        self.work_folder = "calculate_rom_basis_output_process_test_files"
+
+        # List containing the fake total times
+        self.time_steps_list = [1.0,2.0,4.0]
+
+        # Create a model part to perform the ROM basis calculation
+        self.model = KratosMultiphysics.Model()
+        model_part = self.model.CreateModelPart("MainModelPart")
+        model_part.AddNodalSolutionStepVariable(KratosMultiphysics.TEMPERATURE)
+        node_ids_non_ordered = [2,5,12,69,102]
+        n_nodes = len(node_ids_non_ordered)
+        for i in range(n_nodes):
+            model_part.CreateNewNode(node_ids_non_ordered[i],float(i),0.0,0.0)
+
+    def testCalculateRomBasisOutputProcess(self):
+        # Create a CalculateROMBasisOutputProcess instance
+        self.process_settings = KratosMultiphysics.Parameters("""{
+            "model_part_name": "MainModelPart",
+            "snapshots_control_type": "step",
+            "snapshots_interval": 1.0,
+            "nodal_unknowns": ["TEMPERATURE"],
+            "rom_basis_output_format": "numpy",
+            "rom_basis_output_name": "RomParameters_test",
+            "rom_basis_output_folder": "rom_data_test",
+            "svd_truncation_tolerance": 1.0e-6
+        }""")
+
+        # Run a "fake" simulation to calculate the ROM basis from its results
+        self.__ExecuteTest(self.process_settings)
+
+        # Check results
+        self.__CheckResults()
+
+    def tearDown(self):
+        with KratosUnittest.WorkFolderScope(self.work_folder, __file__):
+            if not self.print_output:
+                for path in os.listdir():
+                    full_path = os.path.join(os.getcwd(), path)
+                    if not "_Results" in path:
+                        if os.path.isfile(full_path):
+                            kratos_utilities.DeleteFileIfExisting(full_path)
+                        elif os.path.isdir(full_path):
+                            kratos_utilities.DeleteDirectoryIfExisting(full_path)
+
+
+    def __ExecuteTest(self, process_settings):
+        # Emulate a simulation to get the ROM basis output
+        with KratosUnittest.WorkFolderScope(self.work_folder, __file__):
+            # Create a calculate ROM basis process instance
+            rom_basis_process = CalculateRomBasisOutputProcess(self.model, process_settings)
+
+            # Emulate simulation to fill the database and calculate the ROM basis
+            rom_basis_process.ExecuteInitialize()
+            rom_basis_process.Check()
+            rom_basis_process.ExecuteBeforeSolutionLoop()
+
+            model_part = self.model.GetModelPart("MainModelPart")
+            for time,step in zip(self.time_steps_list, range(1,len(self.time_steps_list)+1,1)):
+                # Fake time advance
+                model_part.CloneTimeStep(time)
+                model_part.ProcessInfo[KratosMultiphysics.STEP] = step
+
+                # Set nodal values
+                if step % 2 == 0:
+                    get_value = lambda node_id, time : node_id*time**2 + step
+                else:
+                    get_value = lambda node_id, time : -node_id/time**step
+                for node in model_part.Nodes:
+                    node.SetSolutionStepValue(KratosMultiphysics.TEMPERATURE, get_value(node.Id, time))
+
+                # Call the process instance methods
+                rom_basis_process.ExecuteInitializeSolutionStep()
+                if rom_basis_process.IsOutputStep():
+                    rom_basis_process.ExecuteBeforeOutputStep()
+                    rom_basis_process.PrintOutput()
+                    rom_basis_process.ExecuteAfterOutputStep()
+                rom_basis_process.ExecuteFinalizeSolutionStep()
+
+            rom_basis_process.ExecuteFinalize()
+
+    def __CheckResults(self):
+        with KratosUnittest.WorkFolderScope(self.work_folder, __file__):
+            # Load ROM basis output file
+            output_folder = self.process_settings["rom_basis_output_folder"].GetString()
+            RightBasisOutput_name = os.path.join(output_folder, "{}.{}".format("RightBasisMatrix", "npy")) #TODO allow name customization in Numpy
+            NodeIdsOutput_name = os.path.join(output_folder, "{}.{}".format("NodeIds", "npy")) #TODO allow name customization in Numpy
+            RightBasisOutput = np.load(RightBasisOutput_name)
+            NodeIdsOutput = np.load(NodeIdsOutput_name)
+
+            # Load reference file
+            RightBasisReference_name = "{}_Results.{}".format("RightBasisMatrix", "npy")
+            NodeIdsReference_name = "{}_Results.{}".format("NodeIds", "npy")
+            RightBasisReference = np.load(RightBasisReference_name)
+            NodeIdsReference = np.load(NodeIdsReference_name)
+
+
+
+            for node_output, node_reference in zip(NodeIdsOutput,NodeIdsReference):
+                self.assertEqual(node_output, node_reference)
+            for i in range(RightBasisOutput.shape[0]):
+                for j in range(RightBasisOutput.shape[1]):
+                    self.assertAlmostEqual(RightBasisOutput[i,j], RightBasisReference[i,j])
+
+            #JSON part tested in the JSON test already
+
+
+##########################################################################################
+
+if __name__ == '__main__':
+    KratosMultiphysics.Logger.GetDefaultOutput().SetSeverity(KratosMultiphysics.Logger.Severity.WARNING)
+    KratosUnittest.main()