Repository: fatiando/harmonica
Pull request #334 (2022-08-12 19:00:04)

Discard prisms with no volume or zero density
Discard null prisms before passing them to the forward modelling
functions. A null prism is the one that has either zero volume or zero
density. Since they don't contribute to the generated field, we can save
some computation time by discarding them.

**Relevant issues/PRs:**

Related to #65 (@aguspesce noticed some zero division errors when passing zero volume tesseroids to the forward model computation).



Patch:
diff --git a/harmonica/forward/prism.py b/harmonica/forward/prism.py
index 4d22c29..3e4d202 100644
--- a/harmonica/forward/prism.py
+++ b/harmonica/forward/prism.py
@@ -141,6 +141,8 @@ def prism_gravity(
                 + "mismatch the number of prisms ({})".format(prisms.shape[0])
             )
         _check_prisms(prisms)
+    # Discard null prisms (zero volume or zero density)
+    prisms, density = _discard_null_prisms(prisms, density)
     # Show progress bar for 'jit_prism_gravity' function
     if progressbar:
         if ProgressBar is None:
@@ -186,8 +188,6 @@ def _check_prisms(prisms):
         ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
         The array must have the following shape: (``n_prisms``, 6), where
         ``n_prisms`` is the total number of prisms.
-        This array of prisms must have valid boundaries.
-        Run ``_check_prisms`` before.
     """
     west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6))
     err_msg = "Invalid prism or prisms. "
@@ -211,6 +211,43 @@ def _check_prisms(prisms):
         raise ValueError(err_msg)
 
 
+def _discard_null_prisms(prisms, density):
+    """
+    Discard prisms with zero volume or zero density
+
+    Parameters
+    ----------
+    prisms : 2d-array
+        Array containing the boundaries of the prisms in the following order:
+        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        The array must have the following shape: (``n_prisms``, 6), where
+        ``n_prisms`` is the total number of prisms.
+        This array of prisms must have valid boundaries.
+        Run ``_check_prisms`` before.
+    density : 1d-array
+        Array containing the density of each prism in kg/m^3. Must have the
+        same size as the number of prisms.
+
+    Returns
+    -------
+    prisms : 2d-array
+        A copy of the ``prisms`` array that doesn't include the null prisms
+        (prisms with zero volume or zero density).
+    density : 1d-array
+        A copy of the ``density`` array that doesn't include the density values
+        for null prisms (prisms with zero volume or zero density).
+    """
+    west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6))
+    # Mark prisms with zero volume as null prisms
+    null_prisms = (west == east) | (south == north) | (bottom == top)
+    # Mark prisms with zero density as null prisms
+    null_prisms[density == 0] = True
+    # Keep only non null prisms
+    prisms = prisms[np.logical_not(null_prisms), :]
+    density = density[np.logical_not(null_prisms)]
+    return prisms, density
+
+
 def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy=None):
     """
     Compute gravitational field of prisms on computations points


Test patch:
diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py
index d3c5ad5..b254ae7 100644
--- a/harmonica/tests/test_prism.py
+++ b/harmonica/tests/test_prism.py
@@ -19,7 +19,13 @@ try:
 except ImportError:
     ProgressBar = None
 
-from ..forward.prism import _check_prisms, prism_gravity, safe_atan2, safe_log
+from ..forward.prism import (
+    _check_prisms,
+    _discard_null_prisms,
+    prism_gravity,
+    safe_atan2,
+    safe_log,
+)
 from ..gravity_corrections import bouguer_correction
 from .utils import run_only_with_numba
 
@@ -67,6 +73,51 @@ def test_invalid_prisms():
         _check_prisms(np.atleast_2d([w, e, s, n, top, bottom]))
 
 
+def test_discard_null_prisms():
+    """
+    Test if discarding invalid prisms works as expected
+    """
+    # Define a set of sample prisms, including invalid ones
+    prisms = np.array(
+        [
+            [0, 10, -50, 33, -2e3, 150],  # ok prism
+            [-10, 0, 33, 66, -1e3, -100],  # ok prism (will set zero density)
+            [7, 7, -50, 50, -3e3, -1e3],  # no volume due to easting bounds
+            [-50, 50, 7, 7, -3e3, -1e3],  # no volume due to northing bounds
+            [-50, 50, -50, 50, -3e3, -3e3],  # no volume due to upward bounds
+            [7, 7, 7, 7, -200, -200],  # no volume due to multiple bounds
+        ]
+    )
+    density = np.array([2670, 0, 3300, 3000, 2800, 2700])
+    prisms, density = _discard_null_prisms(prisms, density)
+    npt.assert_allclose(prisms, np.array([[0, 10, -50, 33, -2e3, 150]]))
+    npt.assert_allclose(density, np.array([2670]))
+
+
+@pytest.mark.use_numba
+def test_forward_with_null_prisms():
+    """
+    Test if the forward model with null prisms gives sensible results
+    """
+    # Create a set of observation points
+    coordinates = vd.grid_coordinates(
+        region=(-50, 50, -50, 50), shape=(3, 3), extra_coords=0
+    )
+    # Build a set of prisms that includes null ones (no volume or zero density)
+    prisms = [
+        [-100, 100, -200, 200, -10e3, -5e3],  # ok prism
+        [100, 200, 200, 300, -10e3, -4e3],  # no density prism
+        [30, 30, -200, 200, -10e3, -5e3],  # no volume (easting)
+        [100, 200, 30, 30, -10e3, -5e3],  # no volume (northing)
+        [-100, 100, -200, 200, -10e3, -10e3],  # no volume (upward)
+    ]
+    density = [2600, 0, 3000, 3100, 3200]
+    npt.assert_allclose(
+        prism_gravity(coordinates, prisms, density, field="g_z"),
+        prism_gravity(coordinates, [prisms[0]], [density[0]], field="g_z"),
+    )
+
+
 @pytest.mark.use_numba
 def test_disable_checks():
     "Check if the disable_checks flag works properly"