Used prange for some geo build methods

2025-08-07 15:10:13 +08:00 · 2022-11-28 21:55:14 +00:00
--- a/gprMax/cmds_geometry/box.py
+++ b/gprMax/cmds_geometry/box.py
@@ -20,6 +20,8 @@ import logging

 import numpy as np

+import gprMax.config as config
+
 from ..cython.geometry_primitives import build_box
 from ..materials import Material
 from .cmds_geometry import UserObjectGeometry, rotate_2point_object
@@ -49,9 +51,9 @@ class Box(UserObjectGeometry):
        self.axis = axis
        self.angle = angle
        self.origin = origin
-        self.dorotate = True
+        self.do_rotate = True

-    def __dorotate(self):
+    def _do_rotate(self):
        """Perform rotation."""
        pts = np.array([self.kwargs['p1'], self.kwargs['p2']])
        rot_pts = rotate_2point_object(pts, self.axis, self.angle, self.origin)
@@ -66,8 +68,8 @@ class Box(UserObjectGeometry):
            logger.exception(self.__str__() + ' Please specify two points.')
            raise

-        if self.dorotate:
-            self.__dorotate()
+        if self.do_rotate:
+            self._do_rotate()

        # Check materials have been specified
        # Isotropic case
@@ -136,8 +138,9 @@ class Box(UserObjectGeometry):
                # Append the new material object to the materials list
                grid.materials.append(m)

-        build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, 
-                  averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
+        build_box(xs, xf, ys, yf, zs, zf, config.get_model_config().ompthreads, 
+                  numID, numIDx, numIDy, numIDz, averaging, grid.solid, 
+                  grid.rigidE, grid.rigidH, grid.ID)

        dielectricsmoothing = 'on' if averaging else 'off'

--- a/gprMax/cmds_geometry/fractal_box_builder.py
+++ b/gprMax/cmds_geometry/fractal_box_builder.py
@@ -46,7 +46,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                volume.originalzs = volume.zs
                volume.originalzf = volume.zf

-                # Extend the volume to accomodate any rough surfaces, grass, or roots
+                # Extend the volume to accomodate any rough surfaces, grass, 
+                # or roots
                for surface in volume.fractalsurfaces:
                    if surface.surfaceID == 'xminus':
                        if surface.fractalrange[0] < volume.xs:
@@ -73,7 +74,8 @@ class FractalBoxBuilder(UserObjectGeometry):
                            volume.nz += surface.fractalrange[1] - volume.zf
                            volume.zf = surface.fractalrange[1]

-                # If there is only 1 bin then a normal material is being used, otherwise a mixing model
+                # If there is only 1 bin then a normal material is being used, 
+                # otherwise a mixing model
                if volume.nbins == 1:
                    volume.fractalvolume = np.ones((volume.nx, volume.ny, volume.nz),
                                                    dtype=config.sim_config.dtypes['float_or_double'])
@@ -85,7 +87,8 @@ class FractalBoxBuilder(UserObjectGeometry):

                volume.generate_volume_mask()

-                # Apply any rough surfaces and add any surface water to the 3D mask array
+                # Apply any rough surfaces and add any surface water to the 
+                # 3D mask array
                for surface in volume.fractalsurfaces:
                    if surface.surfaceID == 'xminus':
                        for i in range(surface.fractalrange[0], surface.fractalrange[1]):
@@ -303,6 +306,7 @@ class FractalBoxBuilder(UserObjectGeometry):
                data = volume.fractalvolume.astype('int16', order='C')
                mask = volume.mask.copy(order='C')
                build_voxels_from_array_mask(volume.xs, volume.ys, volume.zs,
+                                             config.get_model_config().ompthreads,
                                             waternumID, grassnumID, volume.averaging,
                                             mask, data, grid.solid, grid.rigidE,
                                             grid.rigidH, grid.ID)
@@ -318,6 +322,7 @@ class FractalBoxBuilder(UserObjectGeometry):
                    volume.fractalvolume += volume.mixingmodel.startmaterialnum

                data = volume.fractalvolume.astype('int16', order='C')
-                build_voxels_from_array(volume.xs, volume.ys, volume.zs, 0,
+                build_voxels_from_array(volume.xs, volume.ys, volume.zs, 
+                                        config.get_model_config().ompthreads, 0,
                                        volume.averaging, data, grid.solid,
                                        grid.rigidE, grid.rigidH, grid.ID)
--- a/gprMax/cmds_geometry/geometry_objects_read.py
+++ b/gprMax/cmds_geometry/geometry_objects_read.py
@@ -128,9 +128,9 @@ class GeometryObjectsRead(UserObjectGeometry):
                        f'{matfile}.')
        except KeyError:
            averaging = False
-            build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, 
-                                    data, grid.solid, grid.rigidE, grid.rigidH, 
-                                    grid.ID)
+            build_voxels_from_array(xs, ys, zs, config.get_model_config().ompthreads,
+                                    numexistmaterials, averaging, data, 
+                                    grid.solid, grid.rigidE, grid.rigidH, grid.ID)
            logger.info(self.grid_name(grid) + f'Geometry objects from file ' +
                        f'(voxels only){geofile} inserted at {xs * grid.dx:g}m, ' +
                        f'{ys * grid.dy:g}m, {zs * grid.dz:g}m, with corresponding ' +
--- a/gprMax/cython/geometry_primitives.pyx
+++ b/gprMax/cython/geometry_primitives.pyx
@@ -19,6 +19,7 @@
 import numpy as np
 cimport numpy as np
 np.seterr(divide='raise')
+from cython.parallel import prange

 from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_Ex
 from gprMax.cython.yee_cell_setget_rigid cimport set_rigid_Ey
@@ -323,7 +324,7 @@ cpdef void build_voxel(
    np.int8_t[:, :, :, ::1] rigidE,
    np.int8_t[:, :, :, ::1] rigidH,
    np.uint32_t[:, :, :, ::1] ID
-):
+) nogil:
    """Set values in the solid, rigid and ID arrays for a Yee voxel.

    Args:
@@ -643,6 +644,7 @@ cpdef void build_box(
    int yf,
    int zs,
    int zf,
+    int nthreads,
    int numID,
    int numIDx,
    int numIDy,
@@ -657,6 +659,7 @@ cpdef void build_box(

    Args:
        xs, xf, ys, yf, zs, zf: ints for cell coordinates of entire box.
+        nthreads: int for number of threads to use
        numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
        averaging: bint for whether material property averaging will occur for 
                    the object.
@@ -666,50 +669,56 @@ cpdef void build_box(
    cdef Py_ssize_t i, j, k

    if averaging:
-        for i in range(xs, xf):
+        for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
            for j in range(ys, yf):
                for k in range(zs, zf):
                    solid[i, j, k] = numID
                    unset_rigid_E(i, j, k, rigidE)
                    unset_rigid_H(i, j, k, rigidH)
    else:
-        for i in range(xs, xf):
+        for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
            for j in range(ys, yf):
                for k in range(zs, zf):
                    solid[i, j, k] = numID
                    set_rigid_E(i, j, k, rigidE)
                    set_rigid_H(i, j, k, rigidH)
-
-        for i in range(xs, xf):
-            for j in range(ys, yf + 1):
-                for k in range(zs, zf + 1):
                    ID[0, i, j, k] = numIDx
-
-        for i in range(xs, xf + 1):
-            for j in range(ys, yf):
-                for k in range(zs, zf + 1):
                    ID[1, i, j, k] = numIDy
-
-        for i in range(xs, xf + 1):
-            for j in range(ys, yf + 1):
-                for k in range(zs, zf):
                    ID[2, i, j, k] = numIDz
-
-        for i in range(xs, xf + 1):
-            for j in range(ys, yf):
-                for k in range(zs, zf):
                    ID[3, i, j, k] = numIDx
-
-        for i in range(xs, xf):
-            for j in range(ys, yf + 1):
-                for k in range(zs, zf):
                    ID[4, i, j, k] = numIDy
-
-        for i in range(xs, xf):
-            for j in range(ys, yf):
-                for k in range(zs, zf + 1):
                    ID[5, i, j, k] = numIDz

+        for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
+            j = yf
+            k = zf
+            ID[0, i, j, k] = numIDx
+
+        i = xf
+        for j in prange(ys, yf, nogil=True, schedule='static', num_threads=nthreads):
+            for k in range(zf, zf + 1):
+                ID[1, i, j, k] = numIDy
+
+        i = xf
+        j = yf
+        for k in prange(zs, zf, nogil=True, schedule='static', num_threads=nthreads):
+            ID[2, i, j, k] = numIDz
+
+        i = xf
+        for j in prange(ys, yf, nogil=True, schedule='static', num_threads=nthreads):
+            for k in range(zs, zf):
+                ID[3, i, j, k] = numIDx
+
+        for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
+            j = yf
+            for k in range(zs, zf):
+                ID[4, i, j, k] = numIDy
+
+        for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
+            for j in range(ys, yf):
+                k = zf
+                ID[5, i, j, k] = numIDz
+

 cpdef void build_cylinder(
    float x1,
@@ -748,7 +757,8 @@ cpdef void build_cylinder(

    cdef Py_ssize_t i, j, k
    cdef int xs, xf, ys, yf, zs, zf, xc, yc, zc
-    cdef float f1f2mag, f2f1mag, f1ptmag, f2ptmag, dot1, dot2, factor1, factor2, theta1, theta2, distance1, distance2
+    cdef float f1f2mag, f2f1mag, f1ptmag, f2ptmag, dot1, dot2, factor1, factor2 
+    cdef float theta1, theta2, distance1, distance2
    cdef bint build, x_align, y_align, z_align
    cdef np.ndarray f1f2, f2f1, f1pt, f2pt

@@ -977,6 +987,7 @@ cpdef void build_voxels_from_array(
    int xs,
    int ys,
    int zs,
+    int nthreads,
    int numexistmaterials,
    bint averaging,
    np.int16_t[:, :, ::1] data,
@@ -990,6 +1001,7 @@ cpdef void build_voxels_from_array(
    Args:
        xs, ys, zs: ints for cell coordinates of position of start of array in 
                    domain.
+        nthreads: int for number of threads to use
        numexistmaterials: int for number of existing materials in model prior 
                            to building voxels.
        averaging: bint for whether material property averaging will occur for 
@@ -1023,12 +1035,12 @@ cpdef void build_voxels_from_array(
    else:
        zf = zs + data.shape[2]

-    for i in range(xs, xf):
+    for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
        for j in range(ys, yf):
            for k in range(zs, zf):
                numID = data[i - xs, j - ys, k - zs]
                if numID >= 0:
-                    numID += numexistmaterials
+                    numID = numID + numexistmaterials
                    build_voxel(i, j, k, numID, numID, numID, numID, averaging, solid, rigidE, rigidH, ID)


@@ -1036,6 +1048,7 @@ cpdef void build_voxels_from_array_mask(
    int xs,
    int ys,
    int zs,
+    int nthreads,
    int waternumID,
    int grassnumID,
    bint averaging,
@@ -1050,6 +1063,7 @@ cpdef void build_voxels_from_array_mask(

    Args:
        xs, ys, zs: ints for cell coordinates of position of start of array in domain.
+        nthreads: int for number of threads to use
        waternumID, grassnumID: ints for numeric ID of water and grass materials.
        averaging: bint for whether material property averaging will occur for 
                the object.
@@ -1066,7 +1080,7 @@ cpdef void build_voxels_from_array_mask(
    yf = ys + data.shape[1]
    zf = zs + data.shape[2]

-    for i in range(xs, xf):
+    for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
        for j in range(ys, yf):
            for k in range(zs, zf):
                if mask[i - xs, j - ys, k - zs] == 1:
--- a/gprMax/cython/yee_cell_setget_rigid.pxd
+++ b/gprMax/cython/yee_cell_setget_rigid.pxd
@@ -23,23 +23,23 @@ cimport numpy as np
 # Get and set functions for the rigid electric component array. The rigid array
 # is 4D with the 1st dimension holding the 12 electric edge components of a
 # cell - Ex1, Ex2, Ex3, Ex4, Ey1, Ey2, Ey3, Ey4, Ez1, Ez2, Ez3, Ez4
-cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef void set_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef void set_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef void set_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
-cdef void unset_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE)
+cdef bint get_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef bint get_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef bint get_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef void set_rigid_Ex(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef void set_rigid_Ey(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef void set_rigid_Ez(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef void set_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil
+cdef void unset_rigid_E(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidE) nogil

 # Get and set functions for the rigid magnetic component array. The rigid array
 # is 4D with the 1st dimension holding the 6 magnetic edge components - Hx1,
 # Hx2, Hy1, Hy2, Hz1, Hz2
-cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef void set_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef void set_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef void set_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef void set_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
-cdef void unset_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH)
+cdef bint get_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef bint get_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef bint get_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef void set_rigid_Hx(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef void set_rigid_Hy(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef void set_rigid_Hz(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef void set_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
+cdef void unset_rigid_H(int i, int j, int k, np.int8_t[:, :, :, ::1] rigidH) nogil
--- a/gprMax/cython/yee_cell_setget_rigid.pyx
+++ b/gprMax/cython/yee_cell_setget_rigid.pyx
@@ -28,7 +28,7 @@ cdef bint get_rigid_Ex(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:

    cdef bint result
    result = False
@@ -51,7 +51,7 @@ cdef bint get_rigid_Ey(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    cdef bint result
    result = False
    if rigidE[4, i, j, k]:
@@ -73,7 +73,7 @@ cdef bint get_rigid_Ez(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    cdef bint result
    result = False
    if rigidE[8, i, j, k]:
@@ -95,7 +95,7 @@ cdef void set_rigid_Ex(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    rigidE[0, i, j, k] = True
    if j != 0:
        rigidE[1, i, j - 1, k] = True
@@ -110,7 +110,7 @@ cdef void set_rigid_Ey(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    rigidE[4, i, j, k] = True
    if i != 0:
        rigidE[7, i - 1, j, k] = True
@@ -125,7 +125,7 @@ cdef void set_rigid_Ez(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    rigidE[8, i, j, k] = True
    if i != 0:
        rigidE[9, i - 1, j, k] = True
@@ -140,7 +140,7 @@ cdef void set_rigid_E(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    rigidE[:, i, j, k] = True


@@ -149,7 +149,7 @@ cdef void unset_rigid_E(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidE
-):
+) nogil:
    rigidE[:, i, j, k] = False

 # Get and set functions for the rigid magnetic component array. The rigid array
@@ -160,7 +160,7 @@ cdef bint get_rigid_Hx(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    cdef bint result
    result = False
    if rigidH[0, i, j, k]:
@@ -176,7 +176,7 @@ cdef bint get_rigid_Hy(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    cdef bint result
    result = False
    if rigidH[2, i, j, k]:
@@ -192,7 +192,7 @@ cdef bint get_rigid_Hz(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    cdef bint result
    result = False
    if rigidH[4, i, j, k]:
@@ -208,7 +208,7 @@ cdef void set_rigid_Hx(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    rigidH[0, i, j, k] = True
    if i != 0:
        rigidH[1, i - 1, j, k] = True
@@ -219,7 +219,7 @@ cdef void set_rigid_Hy(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    rigidH[2, i, j, k] = True
    if j != 0:
        rigidH[3, i, j - 1, k] = True
@@ -230,7 +230,7 @@ cdef void set_rigid_Hz(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    rigidH[4, i, j, k] = True
    if k != 0:
        rigidH[5, i, j, k - 1] = True
@@ -241,7 +241,7 @@ cdef void set_rigid_H(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    rigidH[:, i, j, k] = True


@@ -250,5 +250,5 @@ cdef void unset_rigid_H(
    int j, 
    int k, 
    np.int8_t[:, :, :, ::1] rigidH
-):
+) nogil:
    rigidH[:, i, j, k] = False
--- a/gprMax/fractals.py
+++ b/gprMax/fractals.py
@@ -16,10 +16,11 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.

-import gprMax.config as config
 import numpy as np
 from scipy import fftpack

+import gprMax.config as config
+
 from .cython.fractals_generate import generate_fractal2D, generate_fractal3D
 from .utilities.utilities import round_value