From 9c87a1a441c0f6e120a243fdbc7c84d46d54a86b Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Fri, 22 Jul 2022 17:21:48 +0200 Subject: [PATCH 01/10] Adding HROM to fluid analysis ROM --- .../fluid_dynamics_analysis_rom.py | 45 ++++++++++++++++++- 1 file changed, 43 insertions(+), 2 deletions(-) diff --git a/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py b/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py index 68805c9a8da7..1229b6a56908 100644 --- a/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py +++ b/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py @@ -2,21 +2,34 @@ import KratosMultiphysics.FluidDynamicsApplication import KratosMultiphysics.RomApplication as romapp from KratosMultiphysics.RomApplication import python_solvers_wrapper_rom as solver_wrapper +from KratosMultiphysics.RomApplication.empirical_cubature_method import EmpiricalCubatureMethod from KratosMultiphysics.FluidDynamicsApplication.fluid_dynamics_analysis import FluidDynamicsAnalysis import json +import numpy as np class FluidDynamicsAnalysisROM(FluidDynamicsAnalysis): - def __init__(self,model,project_parameters): + def __init__(self,model,project_parameters, hyper_reduction_element_selector = None): KratosMultiphysics.Logger.PrintWarning('\x1b[1;31m[DEPRECATED CLASS] \x1b[0m',"\'FluidDynamicsAnalysisROM\'", "class is deprecated. Use the \'RomAnalysis\' one instead.") super().__init__(model,project_parameters) + if hyper_reduction_element_selector != None : + if hyper_reduction_element_selector == "EmpiricalCubature": + self.hyper_reduction_element_selector = EmpiricalCubatureMethod() + self.time_step_residual_matrix_container = [] + else: + err_msg = "The requested element selection method \"" + hyper_reduction_element_selector + "\" is not in the rom application\n" + err_msg += "Available options are: \"EmpiricalCubature\"" + raise Exception(err_msg) + else: + self.hyper_reduction_element_selector = None + #### Internal functions #### def _CreateSolver(self): """ Create the Solver (and create and import the ModelPart if it is not alread in the model) """ ## Solver construction - with open('RomParameters.json') as rom_parameters: + with open(self.path +'/RomParameters.json') as rom_parameters: rom_settings = KratosMultiphysics.Parameters(rom_parameters.read()) self.project_parameters["solver_settings"].AddValue("rom_settings", rom_settings["rom_settings"]) return solver_wrapper.CreateSolverByParameters(self.model, self.project_parameters["solver_settings"],self.project_parameters["problem_data"]["parallel_type"].GetString()) @@ -42,3 +55,31 @@ def ModifyAfterSolverInitialize(self): aux[j,i] = nodal_modes[Counter][j][i] node.SetValue(romapp.ROM_BASIS, aux ) # ROM basis counter+=1 + if self.hyper_reduction_element_selector != None: + if self.hyper_reduction_element_selector.Name == "EmpiricalCubature": + self.ResidualUtilityObject = romapp.RomResidualsUtility(self._GetSolver().GetComputingModelPart(), self.project_parameters["solver_settings"]["rom_settings"], self._GetSolver()._GetScheme()) + + + def FinalizeSolutionStep(self): + if self.hyper_reduction_element_selector != None: + if self.hyper_reduction_element_selector.Name == "EmpiricalCubature": + print('\n\n\n\nGenerating matrix of residuals') + ResMat = self.ResidualUtilityObject.GetResiduals() + NP_ResMat = np.array(ResMat, copy=False) + self.time_step_residual_matrix_container.append(NP_ResMat) + super().FinalizeSolutionStep() + + def Finalize(self): + super().Finalize() + if self.hyper_reduction_element_selector != None: + if self.hyper_reduction_element_selector.Name == "EmpiricalCubature": + self.residuals_snapshots = self._ObtainBasis() + + def _ObtainBasis(self): + ### Building the Snapshot matrix #### + for i in range (len(self.time_step_residual_matrix_container)): + if i == 0: + SnapshotMatrix = self.time_step_residual_matrix_container[i] + else: + SnapshotMatrix = np.c_[SnapshotMatrix, self.time_step_residual_matrix_container[i]] + return SnapshotMatrix \ No newline at end of file From 6f3f695b5507ebd975280a0257ff22f781931a6a Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Fri, 22 Jul 2022 17:24:34 +0200 Subject: [PATCH 02/10] Adding path argument for simulations in cluster that require it --- .../python_scripts/fluid_dynamics_analysis_rom.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py b/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py index 1229b6a56908..161df65d9023 100644 --- a/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py +++ b/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py @@ -10,7 +10,7 @@ class FluidDynamicsAnalysisROM(FluidDynamicsAnalysis): - def __init__(self,model,project_parameters, hyper_reduction_element_selector = None): + def __init__(self,model,project_parameters, path = '.', hyper_reduction_element_selector = None): KratosMultiphysics.Logger.PrintWarning('\x1b[1;31m[DEPRECATED CLASS] \x1b[0m',"\'FluidDynamicsAnalysisROM\'", "class is deprecated. Use the \'RomAnalysis\' one instead.") super().__init__(model,project_parameters) if hyper_reduction_element_selector != None : @@ -23,6 +23,7 @@ def __init__(self,model,project_parameters, hyper_reduction_element_selector = N raise Exception(err_msg) else: self.hyper_reduction_element_selector = None + self.path = path #### Internal functions #### @@ -41,7 +42,7 @@ def ModifyAfterSolverInitialize(self): """Here is where the ROM_BASIS is imposed to each node""" super().ModifyAfterSolverInitialize() computing_model_part = self._solver.GetComputingModelPart() - with open('RomParameters.json') as f: + with open(self.path + '/RomParameters.json') as f: data = json.load(f) nodal_dofs = len(data["rom_settings"]["nodal_unknowns"]) nodal_modes = data["nodal_modes"] From 7d803134163fd94d87d58dcdd8ccf3dc5a21c7e5 Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Wed, 27 Jul 2022 12:28:05 +0200 Subject: [PATCH 03/10] placing the path attribute in the right place --- .../python_scripts/fluid_dynamics_analysis_rom.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py b/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py index 161df65d9023..260da6a2fc7a 100644 --- a/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py +++ b/applications/RomApplication/python_scripts/fluid_dynamics_analysis_rom.py @@ -12,6 +12,7 @@ class FluidDynamicsAnalysisROM(FluidDynamicsAnalysis): def __init__(self,model,project_parameters, path = '.', hyper_reduction_element_selector = None): KratosMultiphysics.Logger.PrintWarning('\x1b[1;31m[DEPRECATED CLASS] \x1b[0m',"\'FluidDynamicsAnalysisROM\'", "class is deprecated. Use the \'RomAnalysis\' one instead.") + self.path = path super().__init__(model,project_parameters) if hyper_reduction_element_selector != None : if hyper_reduction_element_selector == "EmpiricalCubature": @@ -23,7 +24,7 @@ def __init__(self,model,project_parameters, path = '.', hyper_reduction_element_ raise Exception(err_msg) else: self.hyper_reduction_element_selector = None - self.path = path + #### Internal functions #### From bf293db7a4459f80f0ae680002a0a1abb1a4b91a Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Wed, 27 Jul 2022 12:28:41 +0200 Subject: [PATCH 04/10] adding auxiliary functions for workflow --- .../auxiliary_functions_workflow.py | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 applications/RomApplication/python_scripts/auxiliary_functions_workflow.py diff --git a/applications/RomApplication/python_scripts/auxiliary_functions_workflow.py b/applications/RomApplication/python_scripts/auxiliary_functions_workflow.py new file mode 100644 index 000000000000..c4fb0f164ba9 --- /dev/null +++ b/applications/RomApplication/python_scripts/auxiliary_functions_workflow.py @@ -0,0 +1,20 @@ +import dislib as ds +from dislib.data.array import Array + +def load_blocks_array(blocks, shape, block_size): + if shape[0] < block_size[0] or shape[1] < block_size[1]: + raise ValueError("The block size is greater than the ds-array") + return Array(blocks, shape=shape, top_left_shape=block_size, + reg_shape=block_size, sparse=False) + +def load_blocks_rechunk(blocks, shape, block_size, new_block_size): + if shape[0] < new_block_size[0] or shape[1] < new_block_size[1]: + raise ValueError("The block size requested for rechunk" + "is greater than the ds-array") + final_blocks = [[]] + # Este bucle lo puse por si los Future objects se guardan en una lista, en caso de que la forma de guardarlos cambie, también cambiará un poco este bucle. + # Si blocks se pasa ya como (p. ej) [[Future_object, Future_object]] no hace falta. + for block in blocks: + final_blocks[0].append(block) + arr = load_blocks_array(final_blocks, shape, block_size) + return arr.rechunk(new_block_size) From b4e5dda74b4ef53df74bf59be8e785d519e152ed Mon Sep 17 00:00:00 2001 From: Raul Bravo Date: Wed, 27 Jul 2022 12:29:43 +0200 Subject: [PATCH 05/10] Adding randomized SVD using COMPSs paralelization --- .../python_scripts/parallel_svd.py | 88 +++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 applications/RomApplication/python_scripts/parallel_svd.py diff --git a/applications/RomApplication/python_scripts/parallel_svd.py b/applications/RomApplication/python_scripts/parallel_svd.py new file mode 100644 index 000000000000..dd8fe374a2e3 --- /dev/null +++ b/applications/RomApplication/python_scripts/parallel_svd.py @@ -0,0 +1,88 @@ +import dislib as ds +import numpy as np +from pycompss.api.task import task +from pycompss.api.constraint import constraint +from pycompss.api.api import compss_wait_on +from pycompss.api.parameter import Type, COLLECTION_IN, Depth +from dislib.data.array import Array + +from KratosMultiphysics.RomApplication.auxiliary_functions_workflow import load_blocks_array, load_blocks_rechunk + + +@constraint(computingUnits=2) +@task(Y_blocks={Type: COLLECTION_IN, Depth: 2}, returns=2) +def my_qr(Y_blocks): + Y = np.block(Y_blocks) + Q,R = np.linalg.qr(Y, mode='reduced') + return Q,R + + +@constraint(computingUnits=1) +@task(B_blocks={Type: COLLECTION_IN, Depth: 2}, returns=2) +def my_svd(B_blocks): + B = np.block(B_blocks) + U_hat, s, _ = np.linalg.svd(B, full_matrices=False) + return U_hat, s + + + + +def rsvd(A,desired_rank,oversampling=10,row_splits=10,column_splits=1): + +#-----Dimensions-------- + k = desired_rank + p = k + oversampling + n = A.shape[0] + m = A.shape[1] + A_row_chunk_size = int( n / row_splits) + A_column_chunk_size = int( m / column_splits) + A = A.rechunk((A_row_chunk_size,A_column_chunk_size)) +#-----Matrix Omega Initialization-------- + omega_column_chunk_size = p + omega_row_chunk_size = A_column_chunk_size + Omega = ds.random_array(shape=(m, p), block_size=(omega_row_chunk_size, omega_column_chunk_size) ) # Create a random projection matrix Omega of size mxp, for this test, p is of 110. +#---------------------------------- + +# STEP 1 (DISTRIBUTED): Sample the column space of A. Y results into an nxp matrix. + Y = A @ Omega + +# STEP 2 (SERIAL): Serial QR, to be done in distributed. Q results into a matrix of nxp + Q,_ = my_qr(Y._blocks) + Q=load_blocks_rechunk([Q], shape = (n, p), block_size= (n, p), new_block_size=(A_row_chunk_size, omega_column_chunk_size)) + +#"SERIAL STEPS FINISH" + B = Q.T @ A +# STEP 3 (DISTRIBUTED): Project A into the orthonormal basis. B results into a matrix of pxm. + +#""""SERIAL STEPS START""" + U_hat, s = my_svd(B._blocks) + U_hat = load_blocks_rechunk([U_hat], shape = (p, m), block_size = (p, m), new_block_size=(omega_column_chunk_size, omega_column_chunk_size)) + +#""""SERIAL STEPS FINISH""" + U_hat = U_hat[:,:k] #U_hat results into a matrix of pxk. + U = Q @ U_hat #STEP 5 (DISTRIBUTED): Project the reduced basis into Q. U results into a matrix of nxk. + U = U.collect() #STEP 6 (collecting the required basis): + s = compss_wait_on(s) + s = s[:k] + return U, s + +def truncated_svd(Matrix, epsilon=0): + Matrix = compss_wait_on(Matrix) + M,N=np.shape(Matrix) + dimMATRIX = max(M,N) + U, s, V = np.linalg.svd(Matrix, full_matrices=False) #U --> M xN, V --> N x N + V = V.T + tol = dimMATRIX*np.finfo(float).eps*max(s)/2 + R = np.sum(s > tol) # Definition of numerical rank + if epsilon == 0: + K = R + else: + SingVsq = np.multiply(s,s) + SingVsq.sort() + normEf2 = np.sqrt(np.cumsum(SingVsq)) + epsilon = epsilon*normEf2[-1] #relative tolerance + T = (sum(normEf2 Date: Wed, 27 Jul 2022 17:09:33 +0200 Subject: [PATCH 06/10] Adding tsqr(svd) for MPI. For eFlows4HPC deliverable M20. --- .../RomApplication/python_scripts/tsqr.py | 251 ++++++++++++++++++ .../tests/test_RomApplication.py | 3 +- ...omized_singular_value_decomposition_mpi.py | 81 ++++++ 3 files changed, 334 insertions(+), 1 deletion(-) create mode 100644 applications/RomApplication/python_scripts/tsqr.py create mode 100644 applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py diff --git a/applications/RomApplication/python_scripts/tsqr.py b/applications/RomApplication/python_scripts/tsqr.py new file mode 100644 index 000000000000..437942f551aa --- /dev/null +++ b/applications/RomApplication/python_scripts/tsqr.py @@ -0,0 +1,251 @@ +import os +# os.environ['MKL_NUM_THREADS'] = '1' +# os.environ['OPENBLAS_NUM_THREADS']='1' +from mpi4py import MPI +import numpy as np +import numpy.linalg + + +def next_power_of_2(n): + ''' + Find the next power of 2 of n + ''' + p = 1 + if (n and not(n & (n - 1))): + return n + while (p < n): + p <<= 1 + return p + +# Ai and Bi are the local parts of A and B +# this function implements Ai.T@Bi and gathers it to root + + +def distributed_transpose_mult(Ai, Bi, comm, root): + + AitBi = Ai.T @ Bi # multiplication of local blocks + + if(comm.rank == root): + AtB = np.zeros(AitBi.shape, dtype=np.double) + else: + AtB = 0 + + comm.Reduce([AitBi, MPI.DOUBLE], [AtB, MPI.DOUBLE], op=MPI.SUM, root=0) + return AtB + +# A -= Q Qt A + +def OrthogonalProjector(Ai, Qi, comm): + QtA = distributed_transpose_mult(Qi, Ai, comm, 0) + QtA = comm.bcast(QtA, 0) + Ai -= Qi @ QtA + +# apply randomization + +def Randomize(Ai, rand_size, comm): + omega = 0 + np.random.seed(0) + if(comm.Get_rank() == 0): + omega = np.random.random_sample((Ai.shape[1], rand_size)) + omega = comm.bcast(omega, 0) + return Ai @ omega + +def TotalNorm(Ai, comm): + norm = np.array([0.0]) + local_norm2 = np.array([np.linalg.norm(Ai)**2],dtype=np.double) + #comm.Allreduce([local_norm2, MPI.DOUBLE], [norm, MPI.DOUBLE], op=MPI.SUM) + comm.Allreduce(local_norm2, norm, op=MPI.SUM) + norm = np.sqrt(norm[0]) + return norm + +def tsqr(Ai): + ''' + Parallel QR factorization using Lapack + Q(m,n) is the Q matrix + R(n,n) is the R matrix + ''' + + # Recover rank and size + MPI_COMM = MPI.COMM_WORLD # Communications macro + MPI_RANK = MPI_COMM.Get_rank() # Who are you? who? who? + MPI_SIZE = MPI_COMM.Get_size() # Total number of processors used (workers) + m, n = Ai.shape + # Algorithm 1 from Demmel et al (2012) + # 1: QR Factorization on Ai to obtain Q1i and Ri + Q1i, R = numpy.linalg.qr(Ai) + nextPower = next_power_of_2(MPI_SIZE) + nlevels = int(np.log2(nextPower)) + QW = np.eye(n, dtype=np.double) + C = np.zeros((2 * n, n), np.double) + Q2l = np.zeros((2 * n * nlevels, n), np.double) + blevel = 1 + for ilevel in range(nlevels): + # Store R in the upper part of the C matrix + C[:n, :] = R + # Decide who sends and who recieves, use R as buffer + prank = MPI_RANK ^ blevel + if MPI_RANK & blevel: + if prank < MPI_SIZE: + MPI_COMM.send(R, prank, tag=0) + else: + if prank < MPI_SIZE: + R = MPI_COMM.recv(source=prank, tag=0) + # Store R in the lower part of the C matrix + C[n:, :] = R + # 2: QR from the C matrix, reuse C and R + Q2i, R = np.linalg.qr(C) + # Store Q2i from this level + Q2l[2 * n * ilevel:2 * n * ilevel + 2 * n, :] = Q2i + blevel <<= 1 + # At this point R is correct on processor 0 + # Broadcast R and its part of the Q matrix + if MPI_SIZE > 1: + blevel = 1 << (nlevels - 1) + mask = blevel - 1 + for ilevel in reversed(range(nlevels)): + if MPI_RANK & mask == 0: + # Obtain Q2i for this level - use C as buffer + C = Q2l[2 * n * ilevel:2 * n * ilevel + 2 * n, :] + # Multiply by QW either set to identity or allocated to a value + # Store into Q2i + Q2i = C @ QW # matmul(C, QW) + # Communications scheme + prank = MPI_RANK ^ blevel + if MPI_RANK & blevel: + if prank < MPI_SIZE: + C = MPI_COMM.recv(source=prank, tag=0) + # Recover R from the upper part of C and QW from the lower + # part + R = C[:n, :] + QW = C[n:, :] + else: + if prank < MPI_SIZE: + # Set up C matrix for sending + # Store R in the upper part and Q2i on the lower part + # Store Q2i of this rank to QW + C[:n, :] = R + C[n:, :] = Q2i[n:, :] + QW = Q2i[:n, :] + MPI_COMM.send(C, prank, tag=0) + blevel >>= 1 + mask >>= 1 + # Multiply Q1i and QW to obtain Qi + Qi = Q1i @ QW # matmul(Q1i, QW) + return Qi + + + +def randomized_orthogonalization(Ai, comm, mu=0, R=0): + Mlocal, N = Ai.shape # C has dimensions M by N + + #obtain M as the total size + mlocal = np.array([Mlocal]) + mtotal = np.array([0]) + comm.Allreduce([mlocal, MPI.INT], [mtotal, MPI.INT], op=MPI.SUM) + M = mtotal[0] + + c = nC = TotalNorm(Ai, comm) # Norm of the initial residual + + if mu == 0: + mu = max(M, N) * np.finfo(float).eps * nC / 2 # Machine precision parameter + else: + mu = nC*mu #it is relative to the initial norm + + dRmax = np.ceil(0.25 * (min(M, N))) + dRmin = min(1, min(M, N)) + dRmin = max(dRmin, np.ceil(0.05 * min(M, N))) + + if(R==0): + R = np.ceil(0.005 * (min(M, N))) # Initial guess for the rank of C + dR = int(R) + TypeRankEstimate = 1 # Exponential + i = 1 # iteration counter + nC_old = c + R_old = 0 + + Q = None + B = None + + while nC > mu: + Qi = tsqr(Randomize(Ai, dR, comm)) + if Q is not None: + for j in range(4): # repeat application of projector to avoid numeric cancellation + OrthogonalProjector(Qi, Q, comm) # Qi=qr(Qi - Q @(Q.T @ Qi)) + Qi = tsqr(Qi) + Bi = distributed_transpose_mult(Qi, Ai, comm, 0) + Bi = comm.bcast(Bi, 0) # broadcast Bi to all nodes + + Ai -= Qi @ Bi + if i == 1: + Q = Qi + if comm.Get_rank()==0: + B = Bi + else: + Q = np.c_[Q, Qi] # da.concatenate([Q,Qi], axis=1) #how to make this concatenation??? + if comm.Get_rank()==0: + B = np.r_[B, Bi] # da.concatenate([B,Bi], axis=0) + nC = TotalNorm(Ai, comm) + if(comm.Get_rank() == 0): + print('iter = ', i, ' nC = ', nC, ' dR = ', dR, ' R = ', Q.shape[1]) + R_new = Q.shape[1] + if TypeRankEstimate == 0: + Rest = R_old + (R_new - R_old) / (nC - nC_old) * \ + (mu - nC_old) # Estimated rank (linear) + elif TypeRankEstimate == 1: + Rest = R_old + (R_new - R_old) / (np.log(nC) - np.log(nC_old)) * ( + np.log(mu) - np.log(nC_old)) # Logarithmic Rank Estimation + else: + Rest = R_old + dRmin + dR = np.ceil(Rest - R_new) + dR = min(dR, dRmax) + dR = max(dR, dRmin) + Rest = R_new + dR + dR = int(dR) + if Rest >= N: + print("FULL") + Q = 'FULL' + B = np.array([]) + # break + i += 1 + R_old = R + nC_old = nC + R = Rest + + return Q, B + + + +def svd_parallel(Q, B, comm, epsilon=0): + + if(comm.rank == 0): + # u, s, _ = np.linalg.svd(B, full_matrices=False) + M,N=np.shape(B) + dimMATRIX = max(M,N) + u, s, v = np.linalg.svd(B, full_matrices=False) #U --> M xN, V --> N x N + v = v.T + tol = dimMATRIX*np.finfo(float).eps*max(s)/2 + R = np.sum(s > tol) # Definition of numerical rank + if epsilon == 0: + K = R + else: + SingVsq = np.multiply(s,s) + SingVsq.sort() + normEf2 = np.sqrt(np.cumsum(SingVsq)) + epsilon = epsilon*normEf2[-1] #relative tolerance + T = (sum(normEf2 Date: Thu, 28 Jul 2022 10:50:56 +0200 Subject: [PATCH 07/10] Adding test to mpi. --- .../tests/test_RomApplication.py | 2 - .../tests/test_RomApplication_mpi.py | 42 +++++++++++++++++++ ...omized_singular_value_decomposition_mpi.py | 5 +-- .../tests/tests_RomApplication_mpi.py | 42 +++++++++++++++++++ 4 files changed, 86 insertions(+), 5 deletions(-) create mode 100644 applications/RomApplication/tests/test_RomApplication_mpi.py create mode 100644 applications/RomApplication/tests/tests_RomApplication_mpi.py diff --git a/applications/RomApplication/tests/test_RomApplication.py b/applications/RomApplication/tests/test_RomApplication.py index f809183c63e7..612def957898 100644 --- a/applications/RomApplication/tests/test_RomApplication.py +++ b/applications/RomApplication/tests/test_RomApplication.py @@ -12,7 +12,6 @@ from test_empirical_cubature_method import TestEmpiricalCubatureMethod from test_calculate_rom_basis_output_process import TestCalculateRomBasisOutputProcess from test_compressible_potiential_rom import TestCompressiblePotentialRom -from test_randomized_singular_value_decomposition_mpi import TestRandomizedSVDMPI def AssembleTestSuites(): ''' Populates the test suites to run. @@ -40,7 +39,6 @@ def AssembleTestSuites(): smallSuite.addTest(TestRandomizedSVD('test_radomized_svd')) smallSuite.addTest(TestEmpiricalCubatureMethod('test_empirical_cubature_method')) smallSuite.addTests(KratosUnittest.TestLoader().loadTestsFromTestCases([TestCompressiblePotentialRom])) - smallSuite.addTest(TestRandomizedSVDMPI('test_radomized_svd_mpi')) # Create a test suit that contains all the tests from every testCase # in the list: diff --git a/applications/RomApplication/tests/test_RomApplication_mpi.py b/applications/RomApplication/tests/test_RomApplication_mpi.py new file mode 100644 index 000000000000..0b3e477e7e8b --- /dev/null +++ b/applications/RomApplication/tests/test_RomApplication_mpi.py @@ -0,0 +1,42 @@ +# Importing the Kratos Library +import KratosMultiphysics as KM + +if not KM.IsDistributedRun(): + raise Exception("This test script can only be executed in MPI!") + +# Import Kratos "wrapper" for unittests +import KratosMultiphysics.KratosUnittest as KratosUnittest + +# Import the tests or test_classes to create the suits + +# Shell tests +from test_randomized_singular_value_decomposition_mpi import TestRandomizedSVDMPI + +def AssembleTestSuites(): + ''' Populates the test suites to run. + Populates the test suites to run. At least, it should pupulate the suites: + "small", "nighlty" and "all" + Return + ------ + suites: A dictionary of suites + The set of suites with its test_cases added. + ''' + suites = KratosUnittest.KratosSuites + + ### Small MPI tests ######################################################## + smallMPISuite = suites['mpi_small'] + + ### Nightly MPI tests ###################################################### + nightlyMPISuite = suites['mpi_nightly'] + nightlyMPISuite.addTests(smallMPISuite) + nightlyMPISuite.addTest(TestRandomizedSVDMPI('test_radomized_svd_mpi')) + + ### Full MPI set ########################################################### + allMPISuite = suites['mpi_all'] + allMPISuite.addTests(nightlyMPISuite) # already contains the smallMPISuite + + return suites + + +if __name__ == '__main__': + KratosUnittest.runTests(AssembleTestSuites()) diff --git a/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py b/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py index 688ce76f12ec..312f11404f5e 100644 --- a/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py +++ b/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py @@ -13,7 +13,7 @@ except: numpy_and_mpi4py_available = False -def synthetic_matrix(degree, rows = 40,repetitions=4): +def synthetic_matrix(degree, rows = 36,repetitions=4): TestMatrix = np.zeros((rows,degree)) x = np.linspace(0,1,rows) for i in range(degree): @@ -23,7 +23,7 @@ def synthetic_matrix(degree, rows = 40,repetitions=4): class TestRandomizedSVDMPI(KratosUnittest.TestCase): @KratosUnittest.skipUnless(numpy_and_mpi4py_available, "numpy and mpi4py are required for TSQR MPI") - def test_radomized_svd(self): + def test_radomized_svd_mpi(self): comm = MPI.COMM_WORLD # Communications macro svd_truncation_tolerance = 1e-6 rank = comm.Get_rank() @@ -70,7 +70,6 @@ def test_radomized_svd(self): Randomized_Reconstruction = U_local@np.diag(s)@v.T #reconstruct matrix #check that the difference of the reconstruction is below tolerance - print(np.linalg.norm(Randomized_Reconstruction - TestMatrix)) self.assertLess( np.linalg.norm(Randomized_Reconstruction - TestMatrix), svd_truncation_tolerance*np.linalg.norm(TestMatrix)) # Cleaning diff --git a/applications/RomApplication/tests/tests_RomApplication_mpi.py b/applications/RomApplication/tests/tests_RomApplication_mpi.py new file mode 100644 index 000000000000..2ef1f6ea268c --- /dev/null +++ b/applications/RomApplication/tests/tests_RomApplication_mpi.py @@ -0,0 +1,42 @@ +# Importing the Kratos Library +import KratosMultiphysics as KM + +if not KM.IsDistributedRun(): + raise Exception("This test script can only be executed in MPI!") + +# Import Kratos "wrapper" for unittests +import KratosMultiphysics.KratosUnittest as KratosUnittest + +# Import the tests or test_classes to create the suits + +# Shell tests +from adjoint_fluid_test import AdjointFluidTest + +def AssembleTestSuites(): + ''' Populates the test suites to run. + Populates the test suites to run. At least, it should pupulate the suites: + "small", "nighlty" and "all" + Return + ------ + suites: A dictionary of suites + The set of suites with its test_cases added. + ''' + suites = KratosUnittest.KratosSuites + + ### Small MPI tests ######################################################## + smallMPISuite = suites['mpi_small'] + + ### Nightly MPI tests ###################################################### + nightlyMPISuite = suites['mpi_nightly'] + nightlyMPISuite.addTests(smallMPISuite) + nightlyMPISuite.addTest(AdjointFluidTest('testCylinder')) + + ### Full MPI set ########################################################### + allMPISuite = suites['mpi_all'] + allMPISuite.addTests(nightlyMPISuite) # already contains the smallMPISuite + + return suites + + +if __name__ == '__main__': + KratosUnittest.runTests(AssembleTestSuites()) From e04608fac80ca1b0d363714569f389b2e1d09208 Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 28 Jul 2022 11:15:42 +0200 Subject: [PATCH 08/10] Additional comments --- .../test_randomized_singular_value_decomposition_mpi.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py b/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py index 312f11404f5e..d06f96bf096a 100644 --- a/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py +++ b/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py @@ -24,6 +24,13 @@ class TestRandomizedSVDMPI(KratosUnittest.TestCase): @KratosUnittest.skipUnless(numpy_and_mpi4py_available, "numpy and mpi4py are required for TSQR MPI") def test_radomized_svd_mpi(self): + """This algorithm it's being tested for a given tolerance which must be always fullfilled according to random theory + Inside MPI algorithm, random matrices are being seeded ("0") since the mpi algorithm requires it + for all cases. In addition, note that it's being compared by assertLess to a given tolerance, + this must be satisfied for random theory in all cases. + Referring to Brunton's book "There is an extensive literature on random matrix theory, + where the above stereotypes (this case) are almost certainly true, meaning that they are + true with high probability""" comm = MPI.COMM_WORLD # Communications macro svd_truncation_tolerance = 1e-6 rank = comm.Get_rank() From e6fab34cfdd99ea867af342b203f85e6824c21fb Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 28 Jul 2022 11:18:30 +0200 Subject: [PATCH 09/10] Corrected comments --- .../test_randomized_singular_value_decomposition_mpi.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py b/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py index d06f96bf096a..da1c91d97407 100644 --- a/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py +++ b/applications/RomApplication/tests/test_randomized_singular_value_decomposition_mpi.py @@ -24,10 +24,9 @@ class TestRandomizedSVDMPI(KratosUnittest.TestCase): @KratosUnittest.skipUnless(numpy_and_mpi4py_available, "numpy and mpi4py are required for TSQR MPI") def test_radomized_svd_mpi(self): - """This algorithm it's being tested for a given tolerance which must be always fullfilled according to random theory - Inside MPI algorithm, random matrices are being seeded ("0") since the mpi algorithm requires it - for all cases. In addition, note that it's being compared by assertLess to a given tolerance, - this must be satisfied for random theory in all cases. + """According to random theory, this algorithm is being tested for a given tolerance that must always be met. + Random matrices are seeded ("0") within the MPI algorithm because the mpi algorithm requires it in all cases. + Furthermore, note that assertLess compares it to a given tolerance, which must be satisfied in all cases for random theory. Referring to Brunton's book "There is an extensive literature on random matrix theory, where the above stereotypes (this case) are almost certainly true, meaning that they are true with high probability""" From fb0e67b38f61787c395973ad611b9b7493f58436 Mon Sep 17 00:00:00 2001 From: SADPR Date: Thu, 28 Jul 2022 12:17:17 +0200 Subject: [PATCH 10/10] Deleting duplicated file --- .../tests/tests_RomApplication_mpi.py | 42 ------------------- 1 file changed, 42 deletions(-) delete mode 100644 applications/RomApplication/tests/tests_RomApplication_mpi.py diff --git a/applications/RomApplication/tests/tests_RomApplication_mpi.py b/applications/RomApplication/tests/tests_RomApplication_mpi.py deleted file mode 100644 index 2ef1f6ea268c..000000000000 --- a/applications/RomApplication/tests/tests_RomApplication_mpi.py +++ /dev/null @@ -1,42 +0,0 @@ -# Importing the Kratos Library -import KratosMultiphysics as KM - -if not KM.IsDistributedRun(): - raise Exception("This test script can only be executed in MPI!") - -# Import Kratos "wrapper" for unittests -import KratosMultiphysics.KratosUnittest as KratosUnittest - -# Import the tests or test_classes to create the suits - -# Shell tests -from adjoint_fluid_test import AdjointFluidTest - -def AssembleTestSuites(): - ''' Populates the test suites to run. - Populates the test suites to run. At least, it should pupulate the suites: - "small", "nighlty" and "all" - Return - ------ - suites: A dictionary of suites - The set of suites with its test_cases added. - ''' - suites = KratosUnittest.KratosSuites - - ### Small MPI tests ######################################################## - smallMPISuite = suites['mpi_small'] - - ### Nightly MPI tests ###################################################### - nightlyMPISuite = suites['mpi_nightly'] - nightlyMPISuite.addTests(smallMPISuite) - nightlyMPISuite.addTest(AdjointFluidTest('testCylinder')) - - ### Full MPI set ########################################################### - allMPISuite = suites['mpi_all'] - allMPISuite.addTests(nightlyMPISuite) # already contains the smallMPISuite - - return suites - - -if __name__ == '__main__': - KratosUnittest.runTests(AssembleTestSuites())