Merge pull request #494 from gprMax/mpi

Improve MPI fractal surface generation performance

gprMax/fractals/fractal_surface.py

@@ -298,16 +298,31 @@ class MPIFractalSurface(FractalSurface):
         A_shape = np.array(A.shape)
         A_substart = np.array(A.substart)
 
-        # 3D array of random numbers to be convolved with the fractal function
+        # 2D array of random numbers to be convolved with the fractal function
         rng = np.random.default_rng(seed=self.seed)
 
-        for index in np.ndindex(*A.global_shape):
-            index = np.array(index)
-            if any(index < A_substart) or any(index >= A_substart + A_shape):
-                rng.standard_normal()
-            else:
-                index -= A_substart
-                A[index[0], index[1]] = rng.standard_normal()
+        # We need to generate random numbers for the whole domain in the
+        # correct order (and throw away ones we don't need) to ensure
+        # reproducibility when running with MPI domain decomposition
+
+        # We use the following terms:
+        # x - number of rows
+        # y - number of cells
+        cells_per_row = A.global_shape[Dim.Y]
+
+        skip_to_next_row = cells_per_row - A_shape[Dim.Y]
+
+        # Skip to the start of the fractal surface
+        rng.standard_normal(size=A_substart[Dim.X] * cells_per_row)
+        rng.standard_normal(size=A_substart[Dim.Y])
+
+        # Generate numbers for the first row
+        A[0, :] = rng.standard_normal(size=A_shape[Dim.Y])
+
+        # Generate numbers for the remaining rows
+        for row in range(1, A_shape[Dim.X]):
+            rng.standard_normal(size=skip_to_next_row)
+            A[row, :] = rng.standard_normal(size=A_shape[Dim.Y])
 
         A_hat = newDistArray(fft)
         assert isinstance(A_hat, DistArray)

gprMax/fractals/fractal_volume.py

@@ -369,34 +369,37 @@ class MPIFractalVolume(FractalVolume):
         # We need to generate random numbers for the whole domain in the
         # correct order (and throw away ones we don't need) to ensure
         # reproducibility when running with MPI domain decomposition
+
+        # We use the following terms:
+        # x - number of planes
+        # y - number of rows
+        # z - number of cells
         cells_per_row = A.global_shape[Dim.Z]
         cells_per_plane = A.global_shape[Dim.Y] * cells_per_row
 
-        # Skip forward in the x dimension
-        planes_to_skip = A_substart[Dim.X]
-        rng.standard_normal(size=planes_to_skip * cells_per_plane)
+        skip_to_next_row = cells_per_row - A_shape[Dim.Z]
+        skip_to_next_plane = cells_per_plane - (A_shape[Dim.Y] * cells_per_row)
 
-        for plane in range(A_shape[Dim.X]):
-            # Skip forward in the y dimension
-            rows_to_skip = A_substart[Dim.Y]
-            rng.standard_normal(size=rows_to_skip * cells_per_row)
+        # Skip to the start of the fractal volume
+        rng.standard_normal(size=A_substart[Dim.X] * cells_per_plane)
+        rng.standard_normal(size=A_substart[Dim.Y] * cells_per_row)
+        rng.standard_normal(size=A_substart[Dim.Z])
 
+        # Generate numbers for the first row of the first plane
+        A[0, 0, :] = rng.standard_normal(size=A_shape[Dim.Z])
+
+        # Generate remaining numbers for the first plane
+        for row in range(1, A_shape[Dim.Y]):
+            rng.standard_normal(size=skip_to_next_row)
+            A[0, row, :] = rng.standard_normal(size=A_shape[Dim.Z])
+
+        # Generate numbers for the remaining planes
+        for plane in range(1, A_shape[Dim.X]):
+            rng.standard_normal(size=skip_to_next_plane)
             for row in range(A_shape[Dim.Y]):
-                # Skip forward in the z dimension
-                columns_to_skip = A_substart[Dim.Z]
-                rng.standard_normal(size=columns_to_skip)
-
-                # Generate column of numbers in the z dimension
+                rng.standard_normal(size=skip_to_next_row)
                 A[plane, row, :] = rng.standard_normal(size=A_shape[Dim.Z])
 
-                # Skip rest of the z dimension
-                columns_to_skip = A.global_shape[Dim.Z] - columns_to_skip - A_shape[Dim.Z]
-                rng.standard_normal(size=columns_to_skip)
-
-            # Skip rest of the y dimension
-            rows_to_skip = A.global_shape[Dim.Y] - rows_to_skip - A_shape[Dim.Y]
-            rng.standard_normal(size=rows_to_skip * cells_per_row)
-
         A_hat = newDistArray(fft)
         assert isinstance(A_hat, DistArray)
 

reframe_tests/tests/test_geometry.py

@@ -238,6 +238,9 @@ class TestBoxGeometryNoPmlMpi(MpiMixin, TestBoxGeometryNoPml):
     test_dependency = TestBoxGeometryNoPml
 
 
+# Fails for the 'non_axis_aligned_cone' model due to slight differences
+# in the model built by the MPI and non-MPI implementations. This is
+# caused by floating point errors when building the geometry.
 @rfm.simple_test
 class TestConeGeometryMpi(MpiMixin, TestConeGeometry):
     tags = {"test", "mpi", "geometery", "cone"}
@@ -245,6 +248,9 @@ class TestConeGeometryMpi(MpiMixin, TestConeGeometry):
     test_dependency = TestConeGeometry
 
 
+# Fails for the 'non_axis_aligned_cylinder' model due to slight
+# differences in the model built by the MPI and non-MPI implementations.
+# This is caused by floating point errors when building the geometry.
 @rfm.simple_test
 class TestCylinderGeometryMpi(MpiMixin, TestCylinderGeometry):
     tags = {"test", "mpi", "geometery", "cylindrical", "sector", "cylindrical_sector"}
@@ -319,6 +325,9 @@ class TestSphereGeometryMpi(MpiMixin, TestSphereGeometry):
     test_dependency = TestSphereGeometry
 
 
+# Fails for the 'triangle_z_rigid' model due to slight differences in
+# the model built by the MPI and non-MPI implementations. This is
+# caused by floating point errors when building the geometry.
 @rfm.simple_test
 class TestTriangleGeometryMpi(MpiMixin, TestTriangleGeometry):
     tags = {"test", "mpi", "geometery", "triangle"}

reframe_tests/tests/test_geometry_views.py

@@ -28,6 +28,10 @@ class TestGeometryViewVoxelMPI(MpiMixin, TestGeometryViewVoxel):
     test_dependency = TestGeometryViewVoxel
 
 
+# Fails as the VTKHDF file format for unstructured grids uses seperate
+# partitions for each MPI rank. This means the internal structure of the
+# file is dependant on the number MPI ranks and just directly comparing
+# the files does not check correctness.
 @rfm.simple_test
 class TestGeometryViewFineMPI(MpiMixin, TestGeometryViewFine):
     tags = {"test", "mpi", "geometry only", "geometry view"}