Repository: fatiando/harmonica
Pull request #400 (2023-04-17 23:02:42)

Forward models of prisms gravity fields with Choclo
Replace the prisms gravity kernels for Choclo functions. Extend the capabilities of the `prism_gravity` function allowing to compute all the gravitational acceleration and tensor components. Add tests that compare Harmonica results against dumb Choclo calls and a Laplacian equation test for the diagonal components of the tensor. Ditch symmetry tests that were already covered by Choclo. Improve docstrings. Raise warning if any computation point falls in a singular point for the tensor components. Add tests for these new checks and warnings.


santisoler: Thanks @LL-Geo! I just pushed a few fixes. We can ignore codecov failing: we have less code now, since was moved to Choclo, so codecov sees that we are running less coverage percentage, but we are actually testing everything we want to test. I'm merging this!


Patch:
diff --git a/doc/user_guide/forward_modelling/prism.rst b/doc/user_guide/forward_modelling/prism.rst
index 1f6dfbc..9093b14 100644
--- a/doc/user_guide/forward_modelling/prism.rst
+++ b/doc/user_guide/forward_modelling/prism.rst
@@ -48,6 +48,97 @@ computation point:
    print(potential, "J/kg")
 
 
+Gravitational fields
+^^^^^^^^^^^^^^^^^^^^
+
+The :func:`harmonica.prism_gravity` is able to compute the gravitational
+potential (``"potential"``), the acceleration components (``"g_e"``, ``"g_n"``,
+``"g_z"``), and tensor components (``"g_ee"``, ``"g_nn"``, ``"g_zz"``,
+``"g_en"``, ``"g_ez"``, ``"g_nz"``).
+
+
+Build a regular grid of computation points located a 10m above the zero height:
+
+.. jupyter-execute::
+
+   import verde as vd
+
+   region = (-10e3, 10e3, -10e3, 10e3)
+   shape = (51, 51)
+   height = 10
+   coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height)
+
+Define a single prism:
+
+.. jupyter-execute::
+
+   prism = [-2e3, 2e3, -2e3, 2e3, -1.6e3, -900]
+   density = 3300
+
+Compute the gravitational fields that this prism generate on each observation
+point:
+
+.. jupyter-execute::
+
+   fields = (
+      "potential",
+      "g_e", "g_n", "g_z",
+      "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"
+   )
+
+   results = {}
+   for field in fields:
+      results[field] = hm.prism_gravity(coordinates, prism, density, field=field)
+
+Plot the results:
+
+.. jupyter-execute::
+
+   import matplotlib.pyplot as plt
+
+   plt.pcolormesh(coordinates[0], coordinates[1], results["potential"])
+   plt.gca().set_aspect("equal")
+   plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
+   plt.colorbar(label="J/kg")
+   plt.show()
+
+
+.. jupyter-execute::
+
+   fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+   for field, ax in zip(("g_e", "g_n", "g_z"), axes):
+      tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+      ax.set_aspect("equal")
+      ax.set_title(field)
+      ax.ticklabel_format(style="sci", scilimits=(0, 0))
+      plt.colorbar(tmp, ax=ax, label="mGal", orientation="horizontal", pad=0.08)
+   plt.show()
+
+.. jupyter-execute::
+
+   fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+   for field, ax in zip(("g_ee", "g_nn", "g_zz"), axes):
+      tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+      ax.set_aspect("equal")
+      ax.set_title(field)
+      ax.ticklabel_format(style="sci", scilimits=(0, 0))
+      plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
+   plt.show()
+
+.. jupyter-execute::
+
+   fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+   for field, ax in zip(("g_en", "g_ez", "g_nz"), axes):
+      tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+      ax.set_aspect("equal")
+      ax.set_title(field)
+      ax.ticklabel_format(style="sci", scilimits=(0, 0))
+      plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
+   plt.show()
+
 Passing multiple prisms
 ^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/environment.yml b/environment.yml
index e7c399c..79deff7 100644
--- a/environment.yml
+++ b/environment.yml
@@ -18,6 +18,7 @@ dependencies:
   - verde>=1.7.0
   - xarray>=0.16
   - xrft>=1.0
+  - choclo>=0.1
   # Optional requirements
   - pyvista>=0.27
   - vtk>=9
diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py
index 7a4a3a7..3ac4013 100644
--- a/harmonica/_forward/prism.py
+++ b/harmonica/_forward/prism.py
@@ -7,17 +7,48 @@
 """
 Forward modelling for prisms
 """
+import warnings
+
 import numpy as np
+from choclo.prism import (
+    gravity_e,
+    gravity_ee,
+    gravity_en,
+    gravity_eu,
+    gravity_n,
+    gravity_nn,
+    gravity_nu,
+    gravity_pot,
+    gravity_u,
+    gravity_uu,
+)
+from choclo.prism._utils import (
+    is_point_on_easting_edge,
+    is_point_on_northing_edge,
+    is_point_on_upward_edge,
+)
 from numba import jit, prange
 
+# Define dictionary with available gravity fields for prisms
+FIELDS = {
+    "potential": gravity_pot,
+    "g_e": gravity_e,
+    "g_n": gravity_n,
+    "g_z": gravity_u,
+    "g_ee": gravity_ee,
+    "g_nn": gravity_nn,
+    "g_zz": gravity_uu,
+    "g_en": gravity_en,
+    "g_ez": gravity_eu,
+    "g_nz": gravity_nu,
+}
+
 # Attempt to import numba_progress
 try:
     from numba_progress import ProgressBar
 except ImportError:
     ProgressBar = None
 
-from ..constants import GRAVITATIONAL_CONST
-
 
 def prism_gravity(
     coordinates,
@@ -32,23 +63,18 @@ def prism_gravity(
     """
     Gravitational fields of right-rectangular prisms in Cartesian coordinates
 
-    The gravitational fields are computed through the analytical solutions
-    given by [Nagy2000]_ and [Nagy2002]_, which are valid on the entire domain.
-    This means that the computation can be done at any point, either outside or
-    inside the prism.
-
-    This implementation makes use of the modified arctangent function proposed
-    by [Fukushima2020]_ (eq. 12) so that the potential field to satisfies
-    Poisson's equation in the entire domain. Moreover, the logarithm function
-    was also modified in order to solve the singularities that the analytical
-    solution has on some points (see [Nagy2000]_).
+    Compute the gravitational potential, gravitational acceleration and tensor
+    components generated by a collection of prisms on a set of observation
+    points.
 
     .. warning::
-        The **z direction points upwards**, i.e. positive and negative values
-        of ``upward`` represent points above and below the surface,
-        respectively. But remember that the ``g_z`` field returns the downward
-        component of the gravitational acceleration so that positive density
-        contrasts produce positive anomalies.
+        The **vertical direction points upwards**, i.e. positive and negative
+        values of ``upward`` represent points above and below the surface,
+        respectively. But ``g_z`` field returns the **downward component** of
+        the gravitational acceleration so that positive density contrasts
+        produce positive anomalies. The same applies to the tensor components,
+        i.e. the ``g_ez`` is the non-diagonal **easting-downward** tensor
+        component.
 
     Parameters
     ----------
@@ -58,11 +84,11 @@ def prism_gravity(
         system. All coordinates should be in meters.
     prisms : list, 1d-array, or 2d-array
         List or array containing the coordinates of the prism(s) in the
-        following order:
-        west, east, south, north, bottom, top in a Cartesian coordinate system.
+        following order: west, east, south, north, bottom, top, in a Cartesian
+        coordinate system.
         All coordinates should be in meters. Coordinates for more than one
         prism can be provided. In this case, *prisms* should be a list of lists
-        or 2d-array (with one prism per line).
+        or 2d-array (with one prism per row).
     density : list or array
         List or array containing the density of each prism in kg/m^3.
     field : str
@@ -70,7 +96,11 @@ def prism_gravity(
         The available fields are:
 
         - Gravitational potential: ``potential``
+        - Eastward acceleration: ``g_e``
+        - Northward acceleration: ``g_n``
         - Downward acceleration: ``g_z``
+        - Diagonal tensor components: ``g_ee``, ``g_nn``, ``g_zz``
+        - Non-diagonal tensor components: ``g_en``, ``g_ez``, ``g_nz``
 
     parallel : bool (optional)
         If True the computations will run in parallel using Numba built-in
@@ -95,6 +125,36 @@ def prism_gravity(
     result : array
         Gravitational field generated by the prisms on the computation points.
 
+    Notes
+    -----
+    This function makes use of :mod:`choclo.prism` forward modelling functions
+    to compute each gravitational field.
+
+    The gravitational potential (``"potential"``) and the acceleration
+    components (``"g_e"``, ``"g_n"`` and ``"g_z"``) are well defined on the
+    entire domain.
+    Tensor components aren't defined on prism vertices.
+    Diagonal tensor components aren't defined on edges normal to the direction
+    of the tensor (e.g. `"g_ee"` is not defined on edges parallel to northing
+    and upward directions).
+    Non-diagonal tensor components aren't defined on edges normal to the
+    remaining direction of the tensor (e.g. `"g_en"` is not defined on edges
+    parallel to the upward direction).
+    The function returns :func:`numpy.nan` on every singular point.
+
+    The diagonal tensor components aren't defined on observation points that
+    belong to the faces normal to their direction (e.g. ``"g_zz"`` is not
+    define on horizontal faces): two different limits exist when approaching
+    from the inside and from the outside of the prism.
+    This function returns the limit of these components while approaching from
+    the outside.
+
+    References
+    ----------
+    * [Nagy2000]_
+    * [Nagy2002]_
+    * [Fukushima2020]_
+
     Examples
     --------
 
@@ -104,9 +164,12 @@ def prism_gravity(
     >>> prism = [-34, 5, -18, 14, -345, -146]
     >>> # Set prism density to 2670 kg/m³
     >>> density = 2670
-    >>> # Define three computation points along the easting axe at 30m above
-    >>> # the surface
-    >>> coordinates = ([-40, 0, 40], [0, 0, 0], [30, 30, 30])
+    >>> # Define three computation points along the easting direction at 30m
+    >>> # above the surface
+    >>> easting = [-40, 0, 40]
+    >>> northing = [0, 0, 0]
+    >>> upward = [30, 30, 30]
+    >>> coordinates = (easting, northing, upward)
     >>> # Compute the downward component of the gravitational acceleration that
     >>> # the prism generates on the computation points
     >>> gz = prism_gravity(coordinates, prism, density, field="g_z")
@@ -123,8 +186,7 @@ def prism_gravity(
     (-0.05380, 0.02908, 0.11237)
 
     """
-    kernels = {"potential": kernel_potential, "g_z": kernel_g_z}
-    if field not in kernels:
+    if field not in FIELDS:
         raise ValueError("Gravitational field {} not recognized".format(field))
     # Figure out the shape and size of the output array
     cast = np.broadcast(*coordinates[:3])
@@ -141,40 +203,252 @@ def prism_gravity(
                 + "mismatch the number of prisms ({})".format(prisms.shape[0])
             )
         _check_prisms(prisms)
+        _check_singular_points(coordinates, prisms, field)
     # 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:
             raise ImportError(
-                "Missing optional dependency 'numba_progress' required if progressbar=True"
+                "Missing optional dependency 'numba_progress' required "
+                "if progressbar=True"
             )
         progress_proxy = ProgressBar(total=coordinates[0].size)
     else:
         progress_proxy = None
+    # Choose parallelized or serialized forward function
+    if parallel:
+        gravity_prism_func = jit_prism_gravity_parallel
+    else:
+        gravity_prism_func = jit_prism_gravity_serial
     # Compute gravitational field
-    dispatcher(parallel)(
-        coordinates, prisms, density, kernels[field], result, progress_proxy
+    gravity_prism_func(
+        coordinates, prisms, density, FIELDS[field], result, progress_proxy
     )
-    result *= GRAVITATIONAL_CONST
     # Close previously created progress bars
-    if progressbar:
+    if progress_proxy:
         progress_proxy.close()
+    # Invert sign of gravity_u, gravity_eu, gravity_nu
+    if field in ("g_z", "g_ez", "g_nz"):
+        result *= -1
     # Convert to more convenient units
-    if field == "g_z":
+    if field in ("g_e", "g_n", "g_z"):
         result *= 1e5  # SI to mGal
+    # Convert to more convenient units
+    if field in ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"):
+        result *= 1e9  # SI to Eotvos
     return result.reshape(cast.shape)
 
 
-def dispatcher(parallel):
+def _check_singular_points(coordinates, prisms, field):
     """
-    Return the parallelized or serialized forward modelling function
+    Check if any observation point is a singular point for tensor components
+
+    The analytic solutions for the tensor components of the prism have some
+    singular points:
+
+      - All prism vertices are singular points for every tensor component.
+      - Diagonal components aren't defined on edges perpendicular to the
+        component direction (e.g. ``g_ee`` is not defined on edges parallel to
+        northing and upward directions).
+      - Non-diagonal components aren't defined on edges perpendicular to the
+        two directions of the component (e.g. ``g_en`` is not defined on edges
+        parallel to the upward direction).
     """
-    dispatchers = {
-        True: jit_prism_gravity_parallel,
-        False: jit_prism_gravity_serial,
+    functions = {
+        "g_ee": _any_singular_point_g_ee,
+        "g_nn": _any_singular_point_g_nn,
+        "g_zz": _any_singular_point_g_zz,
+        "g_en": _any_singular_point_g_en,
+        "g_ez": _any_singular_point_g_ez,
+        "g_nz": _any_singular_point_g_nz,
     }
-    return dispatchers[parallel]
+    if field not in functions:
+        return None
+    if functions[field](coordinates, prisms):
+        warnings.warn(
+            "Found observation point on singular point of a prism.", UserWarning
+        )
+
+
+@jit(nopython=True)
+def _any_singular_point_g_ee(coordinates, prisms):
+    """
+    Check observation points as singular points of g_ee
+    """
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    n_prisms = prisms.shape[0]
+    for l in range(n_coords):
+        for m in range(n_prisms):
+            if is_point_on_northing_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ) or is_point_on_upward_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ):
+                return True
+    return False
+
+
+@jit(nopython=True)
+def _any_singular_point_g_nn(coordinates, prisms):
+    """
+    Check observation points as singular points of g_nn
+    """
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    n_prisms = prisms.shape[0]
+    for l in range(n_coords):
+        for m in range(n_prisms):
+            if is_point_on_easting_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ) or is_point_on_upward_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ):
+                return True
+    return False
+
+
+@jit(nopython=True)
+def _any_singular_point_g_zz(coordinates, prisms):
+    """
+    Check observation points as singular points of g_zz
+    """
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    n_prisms = prisms.shape[0]
+    for l in range(n_coords):
+        for m in range(n_prisms):
+            if is_point_on_easting_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ) or is_point_on_northing_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ):
+                return True
+    return False
+
+
+@jit(nopython=True)
+def _any_singular_point_g_en(coordinates, prisms):
+    """
+    Check observation points as singular points of g_en
+    """
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    n_prisms = prisms.shape[0]
+    for l in range(n_coords):
+        for m in range(n_prisms):
+            if is_point_on_upward_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ):
+                return True
+    return False
+
+
+@jit(nopython=True)
+def _any_singular_point_g_ez(coordinates, prisms):
+    """
+    Check observation points as singular points of g_ez
+    """
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    n_prisms = prisms.shape[0]
+    for l in range(n_coords):
+        for m in range(n_prisms):
+            if is_point_on_northing_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ):
+                return True
+    return False
+
+
+@jit(nopython=True)
+def _any_singular_point_g_nz(coordinates, prisms):
+    """
+    Check observation points as singular points of g_nz
+    """
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    n_prisms = prisms.shape[0]
+    for l in range(n_coords):
+        for m in range(n_prisms):
+            if is_point_on_easting_edge(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+            ):
+                return True
+    return False
 
 
 def _check_prisms(prisms):
@@ -248,7 +522,9 @@ def _discard_null_prisms(prisms, density):
     return prisms, density
 
 
-def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy=None):
+def jit_prism_gravity(
+    coordinates, prisms, density, forward_func, out, progress_proxy=None
+):
     """
     Compute gravitational field of prisms on computations points
 
@@ -266,109 +542,38 @@ def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy=
     density : 1d-array
         Array containing the density of each prism in kg/m^3. Must have the
         same size as the number of prisms.
-    kernel : func
-        Kernel function that will be used to compute the desired field.
+    forward_func : func
+        Forward modelling function that will be used to compute the desired
+        field. It could be one of the forward modelling functions in
+        :mod:`choclo.prism`.
     out : 1d-array
         Array where the resulting field values will be stored.
         Must have the same size as the arrays contained on ``coordinates``.
     """
+    # Unpack coordinates
+    easting, northing, upward = coordinates
     # Check if we need to update the progressbar on each iteration
     update_progressbar = progress_proxy is not None
     # Iterate over computation points and prisms
-    for l in prange(coordinates[0].size):
+    for l in prange(easting.size):
         for m in range(prisms.shape[0]):
-            # Iterate over the prism vertices to compute the result of the
-            # integration (see Nagy et al., 2000)
-            for i in range(2):
-                # Compute shifted easting coordinate
-                shift_east = prisms[m, 1 - i] - coordinates[0][l]
-                for j in range(2):
-                    # Compute shifted northing coordinate
-                    shift_north = prisms[m, 3 - j] - coordinates[1][l]
-                    for k in range(2):
-                        # Compute shifted upward coordinate
-                        shift_upward = prisms[m, 5 - k] - coordinates[2][l]
-                        # If i, j or k is 1, the corresponding shifted
-                        # coordinate will refer to the lower boundary,
-                        # meaning the corresponding term should have a minus
-                        # sign.
-                        out[l] += (
-                            density[m]
-                            * (-1) ** (i + j + k)
-                            * kernel(shift_east, shift_north, shift_upward)
-                        )
+            out[l] += forward_func(
+                easting[l],
+                northing[l],
+                upward[l],
+                prisms[m, 0],
+                prisms[m, 1],
+                prisms[m, 2],
+                prisms[m, 3],
+                prisms[m, 4],
+                prisms[m, 5],
+                density[m],
+            )
         # Update progress bar if called
         if update_progressbar:
             progress_proxy.update(1)
 
 
-@jit(nopython=True)
-def kernel_potential(easting, northing, upward):
-    """
-    Kernel function for potential gravitational field generated by a prism
-    """
-    radius = np.sqrt(easting**2 + northing**2 + upward**2)
-    kernel = (
-        easting * northing * safe_log(upward + radius)
-        + northing * upward * safe_log(easting + radius)
-        + easting * upward * safe_log(northing + radius)
-        - 0.5 * easting**2 * safe_atan2(upward * northing, easting * radius)
-        - 0.5 * northing**2 * safe_atan2(upward * easting, northing * radius)
-        - 0.5 * upward**2 * safe_atan2(easting * northing, upward * radius)
-    )
-    return kernel
-
-
-@jit(nopython=True)
-def kernel_g_z(easting, northing, upward):
-    """
-    Kernel for downward component of gravitational acceleration of a prism
-    """
-    radius = np.sqrt(easting**2 + northing**2 + upward**2)
-    kernel = (
-        easting * safe_log(northing + radius)
-        + northing * safe_log(easting + radius)
-        - upward * safe_atan2(easting * northing, upward * radius)
-    )
-    return kernel
-
-
-@jit(nopython=True)
-def safe_atan2(y, x):
-    """
-    Principal value of the arctangent expressed as a two variable function
-
-    This modification has to be made to the arctangent function so the
-    gravitational field of the prism satisfies the Poisson's equation.
-    Therefore, it guarantees that the fields satisfies the symmetry properties
-    of the prism. This modified function has been defined according to
-    [Fukushima2020]_.
-    """
-    if x != 0:
-        result = np.arctan(y / x)
-    else:
-        if y > 0:
-            result = np.pi / 2
-        elif y < 0:
-            result = -np.pi / 2
-        else:
-            result = 0
-    return result
-
-
-@jit(nopython=True)
-def safe_log(x):
-    """
-    Modified log to return 0 for log(0).
-    The limits in the formula terms tend to 0 (see [Nagy2000]_).
-    """
-    if np.abs(x) < 1e-10:
-        result = 0
-    else:
-        result = np.log(x)
-    return result
-
-
 # Define jitted versions of the forward modelling function
 jit_prism_gravity_serial = jit(nopython=True)(jit_prism_gravity)
 jit_prism_gravity_parallel = jit(nopython=True, parallel=True)(jit_prism_gravity)
diff --git a/setup.cfg b/setup.cfg
index a9d43db..5a51010 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -48,6 +48,7 @@ install_requires =
     xarray>=0.16
     verde>=1.7
     xrft>=1.0
+    choclo>=0.1
 
 [options.extras_require]
 visualizations =


Test patch:
diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py
index 214d159..8a13c94 100644
--- a/harmonica/tests/test_prism.py
+++ b/harmonica/tests/test_prism.py
@@ -13,6 +13,18 @@ import numpy as np
 import numpy.testing as npt
 import pytest
 import verde as vd
+from choclo.prism import (
+    gravity_e,
+    gravity_ee,
+    gravity_en,
+    gravity_eu,
+    gravity_n,
+    gravity_nn,
+    gravity_nu,
+    gravity_pot,
+    gravity_u,
+    gravity_uu,
+)
 
 try:
     from numba_progress import ProgressBar
@@ -20,13 +32,7 @@ except ImportError:
     ProgressBar = None
 
 from .. import bouguer_correction
-from .._forward.prism import (
-    _check_prisms,
-    _discard_null_prisms,
-    prism_gravity,
-    safe_atan2,
-    safe_log,
-)
+from .._forward.prism import _check_prisms, _discard_null_prisms, prism_gravity
 from .utils import run_only_with_numba
 
 
@@ -41,7 +47,7 @@ def test_invalid_field():
 
 def test_invalid_density_array():
     "Check if error is raised when density shape does not match prisms shape"
-    # Create a set of 4 tesseroids
+    # Create a set of 4 prisms
     prisms = [
         [-100, 0, -100, 0, -200, -100],
         [-100, 0, 0, 100, -200, -100],
@@ -90,7 +96,7 @@ def test_discard_null_prisms():
     )
     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(prisms, np.array([[0.0, 10.0, -50.0, 33.0, -2e3, 150.0]]))
     npt.assert_allclose(density, np.array([2670]))
 
 
@@ -139,254 +145,112 @@ def test_disable_checks():
     npt.assert_allclose(invalid_result, -valid_result)
 
 
-@pytest.mark.use_numba
-def test_potential_field_symmetry():
-    "Test if the potential field satisfies symmetry"
-    prism = [-100, 100, -100, 100, -100, 100]
-    density = 2670
-    # Create six outside computation points located on the normal directions to
-    # the prism faces and at the same distance from its center
-    coordinates = (
-        [-200, 200, 0, 0, 0, 0],
-        [0, 0, -200, 200, 0, 0],
-        [0, 0, 0, 0, -200, 200],
-    )
-    result = prism_gravity(coordinates, prism, density, field="potential")
-    npt.assert_allclose(result[0], result)
-    # Create six inside computation points located on the normal directions to
-    # the prism faces and at the same distance from its center
-    coordinates = ([-50, 50, 0, 0, 0, 0], [0, 0, -50, 50, 0, 0], [0, 0, 0, 0, -50, 50])
-    result = prism_gravity(coordinates, prism, density, field="potential")
-    npt.assert_allclose(result[0], result)
-    # Create twelve outside computation points located on the diagonal
-    # directions to the prism faces and at the same distance from its center.
-    # They can be divided into three sets: one made by those points that live
-    # on the horizontal plane that passes through the prism center, and the
-    # other two that live on the pair of vertical and perpendicular planes that
-    # also passes through the center of the prism.
-    coordinates = (
-        [-200, -200, 200, 200, -200, -200, 200, 200, 0, 0, 0, 0],
-        [-200, 200, -200, 200, 0, 0, 0, 0, -200, -200, 200, 200],
-        [0, 0, 0, 0, -200, 200, -200, 200, -200, 200, -200, 200],
-    )
-    result = prism_gravity(coordinates, prism, density, field="potential")
-    npt.assert_allclose(result[0], result)
-    # Create the same twelve points as before, but now all points fall inside
-    # the prism
-    coordinates = (
-        [-50, -50, 50, 50, -50, -50, 50, 50, 0, 0, 0, 0],
-        [-50, 50, -50, 50, 0, 0, 0, 0, -50, -50, 50, 50],
-        [0, 0, 0, 0, -50, 50, -50, 50, -50, 50, -50, 50],
+class TestAgainstChoclo:
+    """
+    Test forward modelling functions against dumb Choclo runs
+    """
+
+    @pytest.fixture()
+    def sample_prisms(self):
+        """
+        Return three sample prisms
+        """
+        prisms = np.array(
+            [
+                [-10, 10, -10, 0, -10, 0],
+                [-10, 0, -10, 10, -20, -10],
+                [5, 15, 5, 15, -15, -5],
+            ],
+            dtype=float,
+        )
+        densities = np.array([400, -200, 300], dtype=float)
+        return prisms, densities
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return four sample observation points
+        """
+        easting = np.array([-5, 10, 0, 15], dtype=float)
+        northing = np.array([14, -4, 11, 0], dtype=float)
+        upward = np.array([9, 6, 6, 12], dtype=float)
+        return (easting, northing, upward)
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize(
+        "field, choclo_func",
+        [
+            ("potential", gravity_pot),
+            ("g_e", gravity_e),
+            ("g_n", gravity_n),
+            ("g_z", gravity_u),
+            ("g_ee", gravity_ee),
+            ("g_nn", gravity_nn),
+            ("g_zz", gravity_uu),
+            ("g_en", gravity_en),
+            ("g_ez", gravity_eu),
+            ("g_nz", gravity_nu),
+        ],
     )
-    result = prism_gravity(coordinates, prism, density, field="potential")
-    npt.assert_allclose(result[0], result)
+    def test_against_choclo(
+        self,
+        field,
+        choclo_func,
+        sample_coordinates,
+        sample_prisms,
+    ):
+        """
+        Tests forward functions against dumb runs on Choclo
+        """
+        easting, northing, upward = sample_coordinates
+        prisms, densities = sample_prisms
+        # Compute expected results with dumb choclo calls
+        expected_result = np.zeros_like(easting)
+        for i in range(easting.size):
+            for j in range(densities.size):
+                expected_result[i] += choclo_func(
+                    easting[i],
+                    northing[i],
+                    upward[i],
+                    prisms[j, 0],
+                    prisms[j, 1],
+                    prisms[j, 2],
+                    prisms[j, 3],
+                    prisms[j, 4],
+                    prisms[j, 5],
+                    densities[j],
+                )
+        if field in ("g_z", "g_ez", "g_nz"):
+            expected_result *= -1  # invert sign
+        if field in ("g_e", "g_n", "g_z"):
+            expected_result *= 1e5  # convert to mGal
+        if field in ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"):
+            expected_result *= 1e9  # convert to Eotvos
+        # Compare with Harmonica results
+        result = prism_gravity(sample_coordinates, prisms, densities, field=field)
+        npt.assert_allclose(result, expected_result)
 
 
-@pytest.mark.use_numba
-def test_g_z_symmetry_outside():
-    """
-    Test if the g_z field satisfies symmetry
-
-    In order to test if the computed g_z satisfies the symmetry of a square
-    prism we will define several set of computation points:
-
-    A. Two points located on the vertical axis of the prism (``easting == 0``
-       and ``northing == 0``), one above and one below the prism at the same
-       distance from its center.
-    B. Four points located on the ``upward == 0`` plane around the prism
-       distributed normally to its faces , i.e. only one of the horizontal
-       coordinates will be nonzero.
-    C. Same as points defined in B, but located on a plane above the prism.
-    D. Same as points defined in B, but located on a plane below the prism.
-    E. Same as points defined in B, but located on a plane slightly above the
-       ``upward == 0`` plane.
-    F. Same as points defined in B, but located on a plane slightly below the
-       ``upward == 0`` plane.
-    G. Four points located on the ``upward == 0`` plane around the prism
-       distributed on the diagonal directions , i.e. both horizontal
-       coordinates will be equal and nonzero.
-    H. Same as points defined in G, but located on an plane above the prism.
-    I. Same as points defined in G, but located on an plane below the prism.
-    J. Same as points defined in G, but located on a plane slightly above the
-       ``upward == 0`` plane.
-    K. Same as points defined in G, but located on a plane slightly below the
-       ``upward == 0`` plane.
-
-    All computation points defined on the previous groups fall outside of the
-    prism.
-
-    The g_z field for a square prism (the horizontal dimensions of the prism
-    are equal) must satisfy the following symmetry rules:
-
-    - The g_z values on points A must be opposite.
-    - The g_z values on points B must be all zero.
-    - The g_z values on points C must be all equal.
-    - The g_z values on points D must be all equal.
-    - The g_z values on points E must be all equal.
-    - The g_z values on points F must be all equal.
-    - The g_z values on points C and D must be opposite.
-    - The g_z values on points E and F must be opposite.
-    - The g_z values on points G must be all zero.
-    - The g_z values on points H must be all equal.
-    - The g_z values on points I must be all equal.
-    - The g_z values on points J must be all equal.
-    - The g_z values on points K must be all equal.
-    - The g_z values on points H and I must be opposite.
-    - The g_z values on points J and K must be opposite.
-    """
-    prism = [-100, 100, -100, 100, -150, 150]
-    density = 2670
-    computation_points = {
-        "A": ([0, 0], [0, 0], [-200, 200]),
-        "B": ([-200, 200, 0, 0], [0, 0, -200, 200], [0, 0, 0, 0]),
-        "C": ([-200, 200, 0, 0], [0, 0, -200, 200], [200, 200, 200, 200]),
-        "D": ([-200, 200, 0, 0], [0, 0, -200, 200], [-200, -200, -200, -200]),
-        "E": ([-200, 200, 0, 0], [0, 0, -200, 200], [1, 1, 1, 1]),
-        "F": ([-200, 200, 0, 0], [0, 0, -200, 200], [-1, -1, -1, -1]),
-        "G": ([-200, 200, 0, 0], [0, 0, -200, 200], [0, 0, 0, 0]),
-        "H": ([-200, -200, 200, 200], [-200, 200, -200, 200], [200, 200, 200, 200]),
-        "I": ([-200, -200, 200, 200], [-200, 200, -200, 200], [-200, -200, -200, -200]),
-        "J": ([-200, -200, 200, 200], [-200, 200, -200, 200], [1, 1, 1, 1]),
-        "K": ([-200, -200, 200, 200], [-200, 200, -200, 200], [-1, -1, -1, -1]),
-    }
-    # Compute g_z on each set of points
-    results = {}
-    for group, coords in computation_points.items():
-        results[group] = prism_gravity(coords, prism, density, field="g_z")
-    # Check symmetries
-    # Values on A must be opposite, and the value of g_z at the point above the
-    # prism must have the same sign as the density, while the one below should
-    # have the opposite
-    npt.assert_allclose(results["A"][0], -results["A"][1])
-    npt.assert_allclose(np.sign(results["A"][0]), -np.sign(density))
-    npt.assert_allclose(np.sign(results["A"][1]), np.sign(density))
-    # Values on C, D, E, F, H, I, J, K must be all equal within each set
-    for group in ["C", "D", "E", "F", "H", "I", "J", "K"]:
-        npt.assert_allclose(results[group][0], results[group])
-    # Values on B and G must be zero
-    for group in ["B", "G"]:
-        npt.assert_allclose(0, results[group])
-    # Values on C and D, E and F, H and I, J and K must be opposite
-    # Moreover, the set of points that are above the prism must have the same
-    # sign as the density, while the ones below should have the opposite
-    for above, below in (("C", "D"), ("E", "F"), ("H", "I"), ("J", "K")):
-        npt.assert_allclose(results[above], -results[below])
-        npt.assert_allclose(np.sign(results[above]), np.sign(density))
-        npt.assert_allclose(np.sign(results[below]), -np.sign(density))
-
-
-def test_g_z_symmetry_inside():
+@run_only_with_numba
+def test_laplace():
     """
-    Test g_z symmetry on computation points that fall inside the prism
-
-    In order to test if the computed g_z satisfies the symmetry of a square
-    prism on computation points that fall inside the prism, we will define
-    several set of computation points:
-
-    A. Two points located on the vertical axis of the prism (``easting == 0``
-       and ``northing == 0``), one above and one below the center of prism,
-       but at the same distance from it.
-    B. Four points located on the ``upward == 0`` plane around the prism
-       distributed normally to its faces , i.e. only one of the horizontal
-       coordinates will be nonzero.
-    C. Same as points defined in B, but located on a plane above the
-       ``upward == 0`` plane.
-    D. Same as points defined in B, but located on a plane below the
-       ``upward == 0`` plane.
-    E. Four points located on the ``upward == 0`` plane around the prism
-       distributed on the diagonal directions , i.e. both horizontal
-       coordinates will be equal and nonzero.
-    F. Same as points defined in E, but located on a plane above the
-       ``upward == 0`` plane.
-    G. Same as points defined in E, but located on a plane below the
-       ``upward == 0`` plane.
-
-    All computation points defined on the previous groups fall outside of the
-    prism.
-
-    The g_z field for a square prism (the horizontal dimensions of the prism
-    are equal) must satisfy the following symmetry rules:
-
-    - The g_z values on points A must be opposite.
-    - The g_z values on points B must be all zero.
-    - The g_z values on points C must be all equal.
-    - The g_z values on points D must be all equal.
-    - The g_z values on points C and D must be opposite.
-    - The g_z values on points E must be all zero.
-    - The g_z values on points F must be all equal.
-    - The g_z values on points G must be all equal.
-    - The g_z values on points F and G must be opposite.
+    Test if the diagonal components satisfy Laplace equation
     """
-    prism = [-100, 100, -100, 100, -150, 150]
-    density = 2670
-    computation_points = {
-        "A": ([0, 0], [0, 0], [-50, 50]),
-        "B": ([-50, 50, 0, 0], [0, 0, -50, 50], [0, 0, 0, 0]),
-        "C": ([-50, 50, 0, 0], [0, 0, -50, 50], [50, 50, 50, 50]),
-        "D": ([-50, 50, 0, 0], [0, 0, -50, 50], [-50, -50, -50, -50]),
-        "E": ([-50, -50, 50, 50], [-50, 50, -50, 50], [0, 0, 0, 0]),
-        "F": ([-50, -50, 50, 50], [-50, 50, -50, 50], [50, 50, 50, 50]),
-        "G": ([-50, -50, 50, 50], [-50, 50, -50, 50], [-50, -50, -50, -50]),
+    region = (-10e3, 10e3, -10e3, 10e3)
+    coords = vd.grid_coordinates(region, shape=(10, 10), extra_coords=300)
+    prisms = [
+        [1e3, 7e3, -5e3, 2e3, -1e3, -500],
+        [-4e3, 1e3, 4e3, 10e3, -2e3, 200],
+    ]
+    densities = [2670.0, 2900.0]
+    diagonal_components = {
+        field: prism_gravity(coords, prisms, densities, field=field)
+        for field in ("g_ee", "g_nn", "g_zz")
     }
-    # Compute g_z on each set of points
-    results = {}
-    for group, coords in computation_points.items():
-        results[group] = prism_gravity(coords, prism, density, field="g_z")
-    # Check symmetries
-    # Values on A must be opposite, and the value of g_z at the point above the
-    # center of the prism must have the same sign as the density, while the one
-    # below should have the opposite
-    npt.assert_allclose(results["A"][0], -results["A"][1])
-    npt.assert_allclose(np.sign(results["A"][0]), -np.sign(density))
-    npt.assert_allclose(np.sign(results["A"][1]), np.sign(density))
-    # Values on C, D, F, G must be all equal within each set
-    for group in ["C", "D", "F", "G"]:
-        npt.assert_allclose(results[group][0], results[group])
-    # Values on B and E must be zero
-    for group in ["B", "E"]:
-        npt.assert_allclose(0, results[group])
-    # Values on C and D, F and G must be opposite
-    # Moreover, the set of points that are above the center of the prism must
-    # have the same sign as the density, while the ones below should have the
-    # opposite
-    for above, below in (("C", "D"), ("F", "G")):
-        npt.assert_allclose(results[above], -results[below])
-        npt.assert_allclose(np.sign(results[above]), np.sign(density))
-        npt.assert_allclose(np.sign(results[below]), -np.sign(density))
-
-
-@pytest.mark.use_numba
-def test_safe_atan2():
-    "Test the safe_atan2 function"
-    # Test safe_atan2 for one point per quadrant
-    # First quadrant
-    x, y = 1, 1
-    npt.assert_allclose(safe_atan2(y, x), np.pi / 4)
-    # Second quadrant
-    x, y = -1, 1
-    npt.assert_allclose(safe_atan2(y, x), -np.pi / 4)
-    # Third quadrant
-    x, y = -1, -1
-    npt.assert_allclose(safe_atan2(y, x), np.pi / 4)
-    # Forth quadrant
-    x, y = 1, -1
-    npt.assert_allclose(safe_atan2(y, x), -np.pi / 4)
-    # Test safe_atan2 if the denominator is equal to zero
-    npt.assert_allclose(safe_atan2(1, 0), np.pi / 2)
-    npt.assert_allclose(safe_atan2(-1, 0), -np.pi / 2)
-    # Test safe_atan2 if both numerator and denominator are equal to zero
-    assert safe_atan2(0, 0) == 0
-
-
-@pytest.mark.use_numba
-def test_safe_log():
-    "Test the safe_log function"
-    # Check if safe_log function satisfies safe_log(0) == 0
-    assert safe_log(0) == 0
-    # Check if safe_log behaves like the natural logarithm in case that x != 0
-    x = np.linspace(1, 100, 101)
-    for x_i in x:
-        npt.assert_allclose(safe_log(x_i), np.log(x_i))
+    npt.assert_allclose(
+        diagonal_components["g_ee"] + diagonal_components["g_nn"],
+        -diagonal_components["g_zz"],
+    )
 
 
 @pytest.mark.use_numba
@@ -394,15 +258,19 @@ def test_prism_against_infinite_slab():
     """
     Test if g_z of a large prism matches the solution for an infinite slab
     """
-    # Define an observation point 1.5 m above the prism
-    coordinates = (0, 0, 1.5)
-    # Define prisms with thickness of 10.5 m and horizontal dimensions from 1e3
-    # to 1e9m and density of 2670 kg/m^3
+    # Define an observation point at 1.5m above zero
+    height = 1.5
+    coordinates = (0, 0, height)
+    # Define prisms with a top surface located at the same height as the
+    # observation point. Each prisms will have a thickness of 10.5m, and
+    # horizontal dimensions from 1e3 to 1e9m, and density of 2670 kg/m^3.
+    top = height
     thickness = 10.5
+    bottom = top - thickness
     sizes = np.logspace(3, 9, 7)
-    bottom, top = -thickness, 0
     density = 2670
-    # Compute the gravity fields generated by each prism
+    # Compute the gravity field generated by each prism
+    # (from smaller to larger)
     results = np.zeros_like(sizes)
     for i, size in enumerate(sizes):
         prism = [-size / 2, size / 2, -size / 2, size / 2, bottom, top]
@@ -410,8 +278,8 @@ def test_prism_against_infinite_slab():
     # Check convergence: assert if as the prism size increases, the result gets
     # closer to the analytical solution for an infinite slab
     analytical = bouguer_correction(np.array(thickness), density)
-    diffs = abs((results - analytical) / analytical)
-    assert (diffs[1:] < diffs[:-1]).all()
+    errors = abs(analytical - results)
+    assert (errors[1:] < errors[:-1]).all()
     # Check if the largest size is close enough to the analytical solution
     npt.assert_allclose(analytical, results[-1])
 
@@ -486,3 +354,128 @@ def test_numba_progress_missing_error():
         prism_gravity(
             coordinates, prisms, densities, field="potential", progressbar=True
         )
+
+
+class TestSingularPoints:
+    """
+    Tests tensor components on singular points of the prism
+    """
+
+    @pytest.fixture
+    def sample_prism(self):
+        """Return a sample prism"""
+        return np.array([-10.3, 5.4, 8.6, 14.3, -30.3, 2.4])
+
+    @pytest.fixture
+    def sample_density(self):
+        """Return a sample density for the sample prism"""
+        return np.array([2900.0])
+
+    def get_vertices(self, prism):
+        """
+        Return the vertices of the prism as points
+        """
+        easting, northing, upward = tuple(
+            c.ravel() for c in np.meshgrid(prism[:2], prism[2:4], prism[4:6])
+        )
+        return easting, northing, upward
+
+    def get_easting_edges_center(self, prism):
+        """
+        Return points on the center of prism edges parallel to easting
+        """
+        easting_c = (prism[0] + prism[1]) / 2
+        northing, upward = tuple(c.ravel() for c in np.meshgrid(prism[2:4], prism[4:6]))
+        easting = np.full_like(northing, easting_c)
+        return easting, northing, upward
+
+    def get_northing_edges_center(self, prism):
+        """
+        Return points on the center of prism edges parallel to northing
+        """
+        northing_c = (prism[2] + prism[3]) / 2
+        easting, upward = tuple(c.ravel() for c in np.meshgrid(prism[0:2], prism[4:6]))
+        northing = np.full_like(easting, northing_c)
+        return easting, northing, upward
+
+    def get_upward_edges_center(self, prism):
+        """
+        Return points on the center of prism edges parallel to upward
+        """
+        upward_c = (prism[4] + prism[5]) / 2
+        easting, northing = tuple(
+            c.ravel() for c in np.meshgrid(prism[0:2], prism[2:4])
+        )
+        upward = np.full_like(easting, upward_c)
+        return easting, northing, upward
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"))
+    def test_on_vertices(self, sample_prism, sample_density, field):
+        """
+        Test tensor components when observation points fall on prism vertices
+        """
+        easting, northing, upward = self.get_vertices(sample_prism)
+        for i in range(easting.size):
+            msg = "Found observation point"
+            with pytest.warns(UserWarning, match=msg):
+                prism_gravity(
+                    (easting[i], northing[i], upward[i]),
+                    sample_prism,
+                    sample_density,
+                    field=field,
+                )
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", ("g_nn", "g_zz", "g_nz"))
+    def test_on_easting_edges(self, sample_prism, sample_density, field):
+        """
+        Test tensor components that have singular points on edges parallel to
+        easting direction
+        """
+        easting, northing, upward = self.get_easting_edges_center(sample_prism)
+        for i in range(easting.size):
+            msg = "Found observation point"
+            with pytest.warns(UserWarning, match=msg):
+                prism_gravity(
+                    (easting[i], northing[i], upward[i]),
+                    sample_prism,
+                    sample_density,
+                    field=field,
+                )
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", ("g_ee", "g_zz", "g_ez"))
+    def test_on_northing_edges(self, sample_prism, sample_density, field):
+        """
+        Test tensor components that have singular points on edges parallel to
+        easting direction
+        """
+        easting, northing, upward = self.get_northing_edges_center(sample_prism)
+        for i in range(easting.size):
+            msg = "Found observation point"
+            with pytest.warns(UserWarning, match=msg):
+                prism_gravity(
+                    (easting[i], northing[i], upward[i]),
+                    sample_prism,
+                    sample_density,
+                    field=field,
+                )
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", ("g_ee", "g_nn", "g_en"))
+    def test_on_upward_edges(self, sample_prism, sample_density, field):
+        """
+        Test tensor components that have singular points on edges parallel to
+        easting direction
+        """
+        easting, northing, upward = self.get_upward_edges_center(sample_prism)
+        for i in range(easting.size):
+            msg = "Found observation point"
+            with pytest.warns(UserWarning, match=msg):
+                prism_gravity(
+                    (easting[i], northing[i], upward[i]),
+                    sample_prism,
+                    sample_density,
+                    field=field,
+                )