Repository: fatiando/harmonica (python)

Issue:
Parallelize forward models
**Description of the desired feature**

Since #203 I found out that the EQL predictions were using a single CPU, achieving very high and annoying computation times.
After parallelizing the `predict` method, the computation time of the gridding significantly dropped.
I think we should also parallelize the forward models of **points**, **prisms** and **tesseroids** and benchmark them against their non-parallel versions.

One possible option is to use `prange` on jitted functions (enabling `parallel=True`), although maybe a Dask related solution would create more robust solutions that could be directly applied to distributed systems.

**Are you willing to help implement and maintain this feature?** Yes

Pull request description:
Optimize forward models
Parallelize the outer loops of `prism_gravity` and `point_mass_gravity`. Add an
optional `parallel` argument to `prism_gravity` and `point_mass_gravity` to run
them in parallel or in serial at will. Refactor and optimize tesseroids forward
model: avoid defining a set of point masses to compute their effect, instead
compute results in runtime. Simplify its code, reduce memory consumption and
prevent tesseroid forward model from using `jit_point_mass_spherical` function
which has been parallelized, but we don't want to use its parallel version when
computing GLQ on tesseroids forward model.

Fixes #204 

**Reminders**:

- [ ] Run `make format` and `make check` to make sure the code follows the style guide.
- [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
- [ ] Add new public functions/methods/classes to `doc/api/index.rst` and the base `__init__.py` file for the package.
- [ ] Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
- [ ] Add your full name, affiliation, and [ORCID](https://orcid.org) (optional) to the `AUTHORS.md` file (if you haven't already) in case you'd like to be listed as an author on the [Zenodo](https://zenodo.org/communities/fatiando) archive of the next release.


Patch:
diff --git a/harmonica/forward/point_mass.py b/harmonica/forward/point_mass.py
index 4829c4f..d5b9501 100644
--- a/harmonica/forward/point_mass.py
+++ b/harmonica/forward/point_mass.py
@@ -8,14 +8,20 @@
 Forward modelling for point masses
 """
 import numpy as np
-from numba import jit
+from numba import jit, prange
 
 from ..constants import GRAVITATIONAL_CONST
 from .utils import check_coordinate_system, distance_cartesian, distance_spherical_core
 
 
 def point_mass_gravity(
-    coordinates, points, masses, field, coordinate_system="cartesian", dtype="float64"
+    coordinates,
+    points,
+    masses,
+    field,
+    coordinate_system="cartesian",
+    parallel=True,
+    dtype="float64",
 ):
     r"""
     Compute gravitational fields of point masses.
@@ -162,6 +168,11 @@ def point_mass_gravity(
         point masses.
         Available coordinates systems: ``cartesian``, ``spherical``.
         Default ``cartesian``.
+    parallel : bool (optional)
+        If True the computations will run in parallel using Numba built-in
+        parallelization. If False, the forward model will run on a single core.
+        Might be useful to disable parallelization if the forward model is run
+        by an already parallelized workflow. Default to True.
     dtype : data-type (optional)
         Data type assigned to resulting gravitational field. Default to
         ``np.float64``.
@@ -175,27 +186,10 @@ def point_mass_gravity(
         The potential is given in SI units, the accelerations in mGal and the
         Marussi tensor components in Eotvos.
     """
-    # Organize dispatchers and kernel functions inside dictionaries
-    dispatchers = {
-        "cartesian": jit_point_mass_cartesian,
-        "spherical": jit_point_mass_spherical,
-    }
-    kernels = {
-        "cartesian": {
-            "potential": kernel_potential_cartesian,
-            "g_z": kernel_g_z_cartesian,
-            "g_northing": kernel_g_northing_cartesian,
-            "g_easting": kernel_g_easting_cartesian,
-        },
-        "spherical": {
-            "potential": kernel_potential_spherical,
-            "g_z": kernel_g_z_spherical,
-        },
-    }
-    # Sanity checks for coordinate_system and field
-    check_coordinate_system(coordinate_system)
-    if field not in kernels[coordinate_system]:
-        raise ValueError("Gravitational field {} not recognized".format(field))
+    # Sanity checks for coordinate_system
+    check_coordinate_system(
+        coordinate_system, valid_coord_systems=("cartesian", "spherical")
+    )
     # Figure out the shape and size of the output array
     cast = np.broadcast(*coordinates[:3])
     result = np.zeros(cast.size, dtype=dtype)
@@ -210,8 +204,9 @@ def point_mass_gravity(
             + "mismatch the number of points ({})".format(points[0].size)
         )
     # Compute gravitational field
-    dispatchers[coordinate_system](
-        *coordinates, *points, masses, result, kernels[coordinate_system][field]
+    kernel = get_kernel(coordinate_system, field)
+    dispatcher(coordinate_system, parallel)(
+        *coordinates, *points, masses, result, kernel
     )
     result *= GRAVITATIONAL_CONST
     # Convert to more convenient units
@@ -220,39 +215,47 @@ def point_mass_gravity(
     return result.reshape(cast.shape)
 
 
-@jit(nopython=True)
-def jit_point_mass_cartesian(
-    easting, northing, upward, easting_p, northing_p, upward_p, masses, out, kernel
-):  # pylint: disable=invalid-name
+def dispatcher(coordinate_system, parallel):
     """
-    Compute gravitational field of point masses in Cartesian coordinates
+    Return the appropriate forward model function
+    """
+    dispatchers = {
+        "cartesian": {
+            True: point_mass_cartesian_parallel,
+            False: point_mass_cartesian_serial,
+        },
+        "spherical": {
+            True: point_mass_spherical_parallel,
+            False: point_mass_spherical_serial,
+        },
+    }
+    return dispatchers[coordinate_system][parallel]
 
-    Parameters
-    ----------
-    easting, northing, upward : 1d-arrays
-        Coordinates of computation points in Cartesian coordinate system.
-    easting_p, northing_p, upward_p : 1d-arrays
-        Coordinates of point masses in Cartesian coordinate system.
-    masses : 1d-array
-        Mass of each point mass in SI units.
-    out : 1d-array
-        Array where the gravitational field on each computation point will be
-        appended.
-        It must have the same size of ``easting``, ``northing`` and ``upward``.
-    kernel : func
-        Kernel function that will be used to compute the gravitational field on
-        the computation points.
+
+def get_kernel(coordinate_system, field):
     """
-    for l in range(easting.size):
-        for m in range(easting_p.size):
-            out[l] += masses[m] * kernel(
-                easting[l],
-                northing[l],
-                upward[l],
-                easting_p[m],
-                northing_p[m],
-                upward_p[m],
-            )
+    Return the appropriate kernel
+    """
+    kernels = {
+        "cartesian": {
+            "potential": kernel_potential_cartesian,
+            "g_z": kernel_g_z_cartesian,
+            "g_northing": kernel_g_northing_cartesian,
+            "g_easting": kernel_g_easting_cartesian,
+        },
+        "spherical": {
+            "potential": kernel_potential_spherical,
+            "g_z": kernel_g_z_spherical,
+            "g_northing": None,
+            "g_easting": None,
+        },
+    }
+    if field not in kernels[coordinate_system]:
+        raise ValueError("Gravitational field '{}' not recognized".format(field))
+    kernel = kernels[coordinate_system][field]
+    if kernel is None:
+        raise NotImplementedError
+    return kernel
 
 
 @jit(nopython=True)
@@ -312,9 +315,69 @@ def kernel_g_easting_cartesian(
 
 
 @jit(nopython=True)
-def jit_point_mass_spherical(
+def kernel_potential_spherical(
+    longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
+):
+    """
+    Kernel function for potential gravitational field in spherical coordinates
+    """
+    distance, _, _ = distance_spherical_core(
+        longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
+    )
+    return 1 / distance
+
+
+@jit(nopython=True)
+def kernel_g_z_spherical(
+    longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
+):
+    """
+    Kernel for downward component of gravitational gradient in spherical coords
+    """
+    distance, cospsi, _ = distance_spherical_core(
+        longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
+    )
+    delta_z = radius - radius_p * cospsi
+    return delta_z / distance ** 3
+
+
+def point_mass_cartesian(
+    easting, northing, upward, easting_p, northing_p, upward_p, masses, out, kernel
+):  # pylint: disable=invalid-name,not-an-iterable
+    """
+    Compute gravitational field of point masses in Cartesian coordinates
+
+    Parameters
+    ----------
+    easting, northing, upward : 1d-arrays
+        Coordinates of computation points in Cartesian coordinate system.
+    easting_p, northing_p, upward_p : 1d-arrays
+        Coordinates of point masses in Cartesian coordinate system.
+    masses : 1d-array
+        Mass of each point mass in SI units.
+    out : 1d-array
+        Array where the gravitational field on each computation point will be
+        appended.
+        It must have the same size of ``easting``, ``northing`` and ``upward``.
+    kernel : func
+        Kernel function that will be used to compute the gravitational field on
+        the computation points.
+    """
+    for l in prange(easting.size):
+        for m in range(easting_p.size):
+            out[l] += masses[m] * kernel(
+                easting[l],
+                northing[l],
+                upward[l],
+                easting_p[m],
+                northing_p[m],
+                upward_p[m],
+            )
+
+
+def point_mass_spherical(
     longitude, latitude, radius, longitude_p, latitude_p, radius_p, masses, out, kernel
-):  # pylint: disable=invalid-name
+):  # pylint: disable=invalid-name,not-an-iterable
     """
     Compute gravitational field of point masses in spherical coordinates
 
@@ -347,7 +410,7 @@ def jit_point_mass_spherical(
     cosphi_p = np.cos(latitude_p)
     sinphi_p = np.sin(latitude_p)
     # Compute gravitational field
-    for l in range(longitude.size):
+    for l in prange(longitude.size):
         for m in range(longitude_p.size):
             out[l] += masses[m] * kernel(
                 longitude[l],
@@ -361,28 +424,9 @@ def jit_point_mass_spherical(
             )
 
 
-@jit(nopython=True)
-def kernel_potential_spherical(
-    longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
-):
-    """
-    Kernel function for potential gravitational field in spherical coordinates
-    """
-    distance, _, _ = distance_spherical_core(
-        longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
-    )
-    return 1 / distance
-
-
-@jit(nopython=True)
-def kernel_g_z_spherical(
-    longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
-):
-    """
-    Kernel for downward component of gravitational gradient in spherical coords
-    """
-    distance, cospsi, _ = distance_spherical_core(
-        longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
-    )
-    delta_z = radius - radius_p * cospsi
-    return delta_z / distance ** 3
+# Define jitted versions of the forward modelling functions
+# pylint: disable=invalid-name
+point_mass_cartesian_serial = jit(nopython=True)(point_mass_cartesian)
+point_mass_cartesian_parallel = jit(nopython=True, parallel=True)(point_mass_cartesian)
+point_mass_spherical_serial = jit(nopython=True)(point_mass_spherical)
+point_mass_spherical_parallel = jit(nopython=True, parallel=True)(point_mass_spherical)
diff --git a/harmonica/forward/prism.py b/harmonica/forward/prism.py
index be86a68..7b4e82b 100644
--- a/harmonica/forward/prism.py
+++ b/harmonica/forward/prism.py
@@ -8,13 +8,19 @@
 Forward modelling for prisms
 """
 import numpy as np
-from numba import jit
+from numba import jit, prange
 
 from ..constants import GRAVITATIONAL_CONST
 
 
 def prism_gravity(
-    coordinates, prisms, density, field, dtype="float64", disable_checks=False
+    coordinates,
+    prisms,
+    density,
+    field,
+    parallel=True,
+    dtype="float64",
+    disable_checks=False,
 ):
     """
     Gravitational fields of right-rectangular prisms in Cartesian coordinates
@@ -59,6 +65,11 @@ def prism_gravity(
         - Gravitational potential: ``potential``
         - Downward acceleration: ``g_z``
 
+    parallel : bool (optional)
+        If True the computations will run in parallel using Numba built-in
+        parallelization. If False, the forward model will run on a single core.
+        Might be useful to disable parallelization if the forward model is run
+        by an already parallelized workflow. Default to True.
     dtype : data-type (optional)
         Data type assigned to the resulting gravitational field. Default to
         ``np.float64``.
@@ -123,7 +134,7 @@ def prism_gravity(
             )
         _check_prisms(prisms)
     # Compute gravitational field
-    jit_prism_gravity(coordinates, prisms, density, kernels[field], result)
+    dispatcher(parallel)(coordinates, prisms, density, kernels[field], result)
     result *= GRAVITATIONAL_CONST
     # Convert to more convenient units
     if field == "g_z":
@@ -131,6 +142,17 @@ def prism_gravity(
     return result.reshape(cast.shape)
 
 
+def dispatcher(parallel):
+    """
+    Return the parallelized or serialized forward modelling function
+    """
+    dispatchers = {
+        True: jit_prism_gravity_parallel,
+        False: jit_prism_gravity_serial,
+    }
+    return dispatchers[parallel]
+
+
 def _check_prisms(prisms):
     """
     Check if prisms boundaries are well defined
@@ -167,10 +189,9 @@ def _check_prisms(prisms):
         raise ValueError(err_msg)
 
 
-@jit(nopython=True)
 def jit_prism_gravity(
     coordinates, prisms, density, kernel, out
-):  # pylint: disable=invalid-name
+):  # pylint: disable=invalid-name,not-an-iterable
     """
     Compute gravitational field of prisms on computations points
 
@@ -195,7 +216,7 @@ def jit_prism_gravity(
         Must have the same size as the arrays contained on ``coordinates``.
     """
     # Iterate over computation points and prisms
-    for l in range(coordinates[0].size):
+    for l in prange(coordinates[0].size):
         for m in range(prisms.shape[0]):
             # Iterate over the prism boundaries to compute the result of the
             # integration (see Nagy et al., 2000)
@@ -284,3 +305,9 @@ def safe_log(x):
     else:
         result = np.log(x)
     return result
+
+
+# Define jitted versions of the forward modelling function
+# pylint: disable=invalid-name
+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/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py
index bd3d86f..9cb7578 100644
--- a/harmonica/forward/tesseroid.py
+++ b/harmonica/forward/tesseroid.py
@@ -14,7 +14,6 @@ from numpy.polynomial.legendre import leggauss
 from ..constants import GRAVITATIONAL_CONST
 from .utils import distance_spherical
 from .point_mass import (
-    jit_point_mass_spherical,
     kernel_potential_spherical,
     kernel_g_z_spherical,
 )
@@ -174,12 +173,10 @@ def tesseroid_gravity(
     distance_size_ratio = distance_size_ratii[field]
     # Get GLQ unscaled nodes, weights and number of nodes for each small
     # tesseroid
-    n_nodes, glq_nodes, glq_weights = glq_nodes_weights(glq_degrees)
+    glq_nodes, glq_weights = glq_nodes_weights(glq_degrees)
     # Initialize arrays to perform memory allocation only once
     stack = np.empty((stack_size, 6), dtype=dtype)
     small_tesseroids = np.empty((max_discretizations, 6), dtype=dtype)
-    point_masses = np.empty((3, n_nodes * max_discretizations), dtype=dtype)
-    weights = np.empty(n_nodes * max_discretizations, dtype=dtype)
     # Compute gravitational field
     jit_tesseroid_gravity(
         coordinates,
@@ -187,12 +184,9 @@ def tesseroid_gravity(
         density,
         stack,
         small_tesseroids,
-        point_masses,
-        weights,
         result,
         distance_size_ratio,
         radial_adaptive_discretization,
-        n_nodes,
         glq_nodes,
         glq_weights,
         kernels[field],
@@ -204,19 +198,16 @@ def tesseroid_gravity(
     return result.reshape(cast.shape)
 
 
-@jit(nopython=True, parallel=True)
+@jit(nopython=True)
 def jit_tesseroid_gravity(
     coordinates,
     tesseroids,
     density,
     stack,
     small_tesseroids,
-    point_masses,
-    weights,
     result,
     distance_size_ratio,
     radial_discretization,
-    n_nodes,
     glq_nodes,
     glq_weights,
     kernel,
@@ -235,6 +226,8 @@ def jit_tesseroid_gravity(
         geocentric coordinate system in the following order:
         ``longitude``, ``latitude``, ``radius``.
         Each element of the tuple must be a 1d array.
+        Both ``longitude`` and ``latitude`` should be in degrees and ``radius``
+        in meters.
     tesseroids : 2d-array
         Array containing the boundaries of each tesseroid:
         ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric
@@ -242,6 +235,8 @@ def jit_tesseroid_gravity(
         The array must have the following shape: (``n_tesseroids``, 6), where
         ``n_tesseroids`` is the total number of tesseroids.
         All tesseroids must have valid boundary coordinates.
+        Horizontal boundaries should be in degrees while radial boundaries
+        should be in meters.
     density : 1d-array
         Density of each tesseroid in SI units.
     stack : 2d-array
@@ -263,9 +258,6 @@ def jit_tesseroid_gravity(
         tesseroid only on the horizontal direction.
         If ``True``, it will perform a three dimensional adaptive
         discretization, splitting the tesseroids on every direction.
-    n_nodes : int
-        Total number of equivalent point masses that will be generated for each
-        split tesseroid.
     glq_nodes : list
         List containing unscaled GLQ nodes.
     glq_weights : list
@@ -273,121 +265,142 @@ def jit_tesseroid_gravity(
     kernel : func
         Kernel function for the gravitational field of point masses.
     """
-    for l in range(tesseroids.shape[0]):
-        tesseroid = tesseroids[l, :]
-        for m in range(coordinates[0].size):
+    # Get coordinates of the observation points
+    # and precompute trigonometric functions
+    longitude, latitude, radius = coordinates[:]
+    longitude_rad = np.radians(longitude)
+    cosphi = np.cos(np.radians(latitude))
+    sinphi = np.sin(np.radians(latitude))
+    for l in range(longitude.size):
+        for m in range(tesseroids.shape[0]):
             # Apply adaptive discretization on tesseroid
             n_splits = _adaptive_discretization(
-                (coordinates[0][m], coordinates[1][m], coordinates[2][m]),
-                tesseroid,
+                (longitude[l], latitude[l], radius[l]),
+                tesseroids[m, :],
                 distance_size_ratio,
                 stack,
                 small_tesseroids,
                 radial_discretization,
             )
-            # Get total number of point masses, their coordinates and weights
-            n_point_masses = n_nodes * n_splits
-            tesseroids_to_point_masses(
-                small_tesseroids[:n_splits],
-                glq_nodes,
-                glq_weights,
-                point_masses,
-                weights,
-            )
-            # Compute gravitational fields
-            jit_point_mass_spherical(
-                coordinates[0][m : m + 1],  # slice lon to pass a single element array
-                coordinates[1][m : m + 1],  # slice lat to pass a single element array
-                coordinates[2][m : m + 1],  # slice rad to pass a single element array
-                point_masses[0, :n_point_masses],
-                point_masses[1, :n_point_masses],
-                point_masses[2, :n_point_masses],
-                density[l] * weights[:n_point_masses],
-                result[m : m + 1],
-                kernel,
-            )
+            # Compute effect of the tesseroid through GLQ
+            for tess_index in range(n_splits):
+                tesseroid = small_tesseroids[tess_index, :]
+                result[l] += gauss_legendre_quadrature(
+                    longitude_rad[l],
+                    cosphi[l],
+                    sinphi[l],
+                    radius[l],
+                    tesseroid,
+                    density[m],
+                    glq_nodes,
+                    glq_weights,
+                    kernel,
+                )
 
 
 @jit(nopython=True)
-def tesseroids_to_point_masses(
-    tesseroids, glq_nodes, glq_weights, point_masses, weights
-):  # pylint: disable=too-many-locals,invalid-name
+def gauss_legendre_quadrature(
+    longitude,
+    cosphi,
+    sinphi,
+    radius,
+    tesseroid,
+    density,
+    glq_nodes,
+    glq_weights,
+    kernel,
+):  # pylint: disable=too-many-locals
     r"""
-    Convert tesseroids to equivalent point masses on nodes of GLQ
+    Compute the effect of a tesseroid on a single observation point through GLQ
 
-    Each tesseroid is converted into a set of point masses located on the
+    The tesseroid is converted into a set of point masses located on the
     scaled nodes of the Gauss-Legendre Quadrature. The number of point masses
     created from each tesseroid is equal to the product of the GLQ degrees for
-    each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). It also
-    compute a weight value for each point mass defined as the product of the
-    GLQ weights for each direction (:math:`W_i^r`, :math:`W_j^\phi`,
-    :math:`W_k^\lambda`), the scale constant :math:`A` and the :math:`\kappa`
-    factor evaluated on the coordinates of the point mass.
+    each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass
+    of each point mass is defined as the product of the tesseroid density
+    (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`,
+    :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and
+    the :math:`\kappa` factor evaluated on the coordinates of the point mass.
 
     Parameters
     ----------
-    tesseroids : 2d-array
-        Array containing the boundaries of each tesseroid:
-        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric
-        spherical coordinate system.
-        The array must have the following shape: (``n_tesseroids``, 6), where
-        ``n_tesseroids`` is the total number of tesseroids.
-        All tesseroids must have valid boundary coordinates.
+    longitude : float
+        Longitudinal coordinate of the observation points in radians.
+    cosphi : float
+        Cosine of the latitudinal coordinate of the observation point in
+        radians.
+    sinphi : float
+        Sine of the latitudinal coordinate of the observation point in
+        radians.
+    radius : float
+        Radial coordinate of the observation point in meters.
+    tesseroids : 1d-array
+        Array containing the boundaries of the tesseroid:
+        ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+        Horizontal boundaries should be in degrees and radial boundaries in
+        meters.
+    density : float
+        Density of the tesseroid in SI units.
     glq_nodes : list
         Unscaled location of GLQ nodes for each direction.
     glq_weights : list
         GLQ weigths for each node for each direction.
-    point_masses : 2d-array
-        Empty array with shape ``(3, n)``, where ``n`` is the total number of
-        point masses computed as the product of number of tesseroids and the
-        GLQ degrees for each direction.
-        The location of the point masses will be located inside this array.
-    weights : 1d-array
-        Empty array with ``n`` elements.
-        It will contain the weight constant for each point mass.
+    kernel : func
+        Kernel function for the gravitational field of point masses.
 
     """
+    # Get tesseroid boundaries
+    w, e, s, n, bottom, top = tesseroid[:]
+    # Calculate the A factor for the tesseroid
+    a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom)
     # Unpack nodes and weights
     lon_nodes, lat_nodes, rad_nodes = glq_nodes[:]
     lon_weights, lat_weights, rad_weights = glq_weights[:]
-    # Recover GLQ degrees from nodes
-    lon_glq_degree = len(lon_nodes)
-    lat_glq_degree = len(lat_nodes)
-    rad_glq_degree = len(rad_nodes)
-    # Convert each tesseroid to a point mass
-    mass_index = 0
-    for tesseroid_index in range(len(tesseroids)):
-        w = tesseroids[tesseroid_index, 0]
-        e = tesseroids[tesseroid_index, 1]
-        s = tesseroids[tesseroid_index, 2]
-        n = tesseroids[tesseroid_index, 3]
-        bottom = tesseroids[tesseroid_index, 4]
-        top = tesseroids[tesseroid_index, 5]
-        A_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom)
-        for i in range(lon_glq_degree):
-            for j in range(lat_glq_degree):
-                for k in range(rad_glq_degree):
-                    # Compute coordinates of each point mass
-                    longitude = 0.5 * (e - w) * lon_nodes[i] + 0.5 * (e + w)
-                    latitude = 0.5 * (n - s) * lat_nodes[j] + 0.5 * (n + s)
-                    radius = 0.5 * (top - bottom) * rad_nodes[k] + 0.5 * (top + bottom)
-                    kappa = radius ** 2 * np.cos(np.radians(latitude))
-                    point_masses[0, mass_index] = longitude
-                    point_masses[1, mass_index] = latitude
-                    point_masses[2, mass_index] = radius
-                    weights[mass_index] = (
-                        A_factor
-                        * kappa
-                        * lon_weights[i]
-                        * lat_weights[j]
-                        * rad_weights[k]
-                    )
-                    mass_index += 1
+    # Compute effect of the tesseroid on the observation point
+    # by iterating over the location of the point masses
+    # (move the iteration along the longitudinal nodes to the bottom for
+    # optimization: reduce the number of times that the trigonometric functions
+    # are evaluated)
+    result = 0.0
+    for j, lat_node in enumerate(lat_nodes):
+        # Get the latitude of the point mass
+        latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s))
+        cosphi_p = np.cos(latitude_p)
+        sinphi_p = np.sin(latitude_p)
+        for k, rad_node in enumerate(rad_nodes):
+            # Get the radius of the point mass
+            radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom)
+            # Get kappa constant for the point mass
+            kappa = radius_p ** 2 * cosphi_p
+            for i, lon_node in enumerate(lon_nodes):
+                # Get the longitude of the point mass
+                longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w))
+                # Compute the mass of the point mass
+                mass = (
+                    density
+                    * a_factor
+                    * kappa
+                    * lon_weights[i]
+                    * lat_weights[j]
+                    * rad_weights[k]
+                )
+                # Add effect of the current point mass to the result
+                result += mass * kernel(
+                    longitude,
+                    cosphi,
+                    sinphi,
+                    radius,
+                    longitude_p,
+                    cosphi_p,
+                    sinphi_p,
+                    radius_p,
+                )
+    return result
 
 
 def glq_nodes_weights(glq_degrees):
     """
-    Calculate GLQ unscaled nodes, weights and total number of nodes
+    Calculate GLQ unscaled nodes and weights
 
     Parameters
     ----------
@@ -397,8 +410,6 @@ def glq_nodes_weights(glq_degrees):
 
     Returns
     -------
-    n_nodes : int
-        Total number of nodes computed as the product of the GLQ degrees.
     glq_nodes : list
         Unscaled GLQ nodes for each direction: ``longitude``, ``latitude``,
         ``radius``.
@@ -408,8 +419,6 @@ def glq_nodes_weights(glq_degrees):
     """
     # Unpack GLQ degrees
     lon_degree, lat_degree, rad_degree = glq_degrees[:]
-    # Get number of point masses
-    n_nodes = np.prod(glq_degrees)
     # Get nodes coordinates and weights
     lon_node, lon_weights = leggauss(lon_degree)
     lat_node, lat_weights = leggauss(lat_degree)
@@ -417,7 +426,7 @@ def glq_nodes_weights(glq_degrees):
     # Reorder nodes and weights
     glq_nodes = (lon_node, lat_node, rad_node)
     glq_weights = (lon_weights, lat_weights, rad_weights)
-    return n_nodes, glq_nodes, glq_weights
+    return glq_nodes, glq_weights
 
 
 @jit(nopython=True)


Test patch:
diff --git a/harmonica/tests/test_point_mass.py b/harmonica/tests/test_point_mass.py
index 8a1d3fe..05795de 100644
--- a/harmonica/tests/test_point_mass.py
+++ b/harmonica/tests/test_point_mass.py
@@ -10,11 +10,14 @@ Test forward modelling for point masses.
 import numpy as np
 import numpy.testing as npt
 import pytest
+import verde as vd
 
 from ..constants import GRAVITATIONAL_CONST
 from ..forward.point_mass import point_mass_gravity
 from ..forward.utils import distance_cartesian
 
+from .utils import require_numba
+
 
 def test_invalid_coordinate_system():
     "Check if invalid coordinate system is passed"
@@ -31,6 +34,25 @@ def test_invalid_coordinate_system():
         )
 
 
+def test_not_implemented_field():
+    """
+    Check if NotImplementedError is raised after asking a non-implemented field
+    """
+    coordinates = [0.0, 0.0, 0.0]
+    point_mass = [0.0, 0.0, 0.0]
+    mass = 1.0
+    coordinate_system = "spherical"
+    for field in ("g_northing", "g_easting"):
+        with pytest.raises(NotImplementedError):
+            point_mass_gravity(
+                coordinates,
+                point_mass,
+                mass,
+                field,
+                coordinate_system,
+            )
+
+
 def test_invalid_field():
     "Check if an invalid gravitational field is passed as argument"
     coordinates = [0.0, 0.0, 0.0]
@@ -574,6 +596,25 @@ def test_g_easting_sign():
     assert np.sign(mass) == -np.sign(results[2])
 
 
+@require_numba
+def test_point_mass_cartesian_parallel():
+    """
+    Check if parallel and serial runs return the same result
+    """
+    region = (2e3, 10e3, -3e3, 5e3)
+    points = vd.scatter_points(region, size=30, extra_coords=-1e3, random_state=0)
+    masses = np.arange(points[0].size)
+    coordinates = vd.grid_coordinates(region=region, spacing=1e3, extra_coords=0)
+    for field in ("potential", "g_z", "g_northing", "g_easting"):
+        result_serial = point_mass_gravity(
+            coordinates, points, masses, field=field, parallel=False
+        )
+        result_parallel = point_mass_gravity(
+            coordinates, points, masses, field=field, parallel=True
+        )
+        npt.assert_allclose(result_serial, result_parallel)
+
+
 # ---------------------------
 # Spherical coordinates tests
 # ---------------------------
@@ -692,3 +733,35 @@ def test_point_mass_potential_on_same_meridian():
                     ),
                     analytical["potential"],
                 )
+
+
+@require_numba
+def test_point_mass_spherical_parallel():
+    """
+    Check if parallel and serial runs return the same result
+    """
+    region = (2, 10, -3, 5)
+    radius = 6400e3
+    points = vd.scatter_points(
+        region, size=30, extra_coords=radius - 10e3, random_state=0
+    )
+    masses = np.arange(points[0].size)
+    coordinates = vd.grid_coordinates(region=region, spacing=1, extra_coords=radius)
+    for field in ("potential", "g_z"):
+        result_serial = point_mass_gravity(
+            coordinates,
+            points,
+            masses,
+            field=field,
+            coordinate_system="spherical",
+            parallel=False,
+        )
+        result_parallel = point_mass_gravity(
+            coordinates,
+            points,
+            masses,
+            field=field,
+            coordinate_system="spherical",
+            parallel=True,
+        )
+        npt.assert_allclose(result_serial, result_parallel)
diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py
index 3e01f9d..81e6f56 100644
--- a/harmonica/tests/test_prism.py
+++ b/harmonica/tests/test_prism.py
@@ -10,10 +10,13 @@ Test forward modelling for prisms.
 import pytest
 import numpy as np
 import numpy.testing as npt
+import verde as vd
 
 from ..gravity_corrections import bouguer_correction
 from ..forward.prism import prism_gravity, _check_prisms, safe_atan2, safe_log
 
+from .utils import require_numba
+
 
 def test_invalid_field():
     "Check if passing an invalid field raises an error"
@@ -358,3 +361,28 @@ def test_prism_against_infinite_slab():
     assert (diffs[1:] < diffs[:-1]).all()
     # Check if the largest size is close enough to the analytical solution
     npt.assert_allclose(analytical, results[-1])
+
+
+@require_numba
+def test_prisms_parallel_vs_serial():
+    """
+    Check if the parallelized run returns the same results as the serial one
+    """
+    prisms = [
+        [-100, 0, -100, 0, -10, 0],
+        [0, 100, -100, 0, -10, 0],
+        [-100, 0, 0, 100, -10, 0],
+        [0, 100, 0, 100, -10, 0],
+    ]
+    densities = [2000, 3000, 4000, 5000]
+    coordinates = vd.grid_coordinates(
+        region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+    )
+    for field in ("potential", "g_z"):
+        result_parallel = prism_gravity(
+            coordinates, prisms, densities, field=field, parallel=True
+        )
+        result_serial = prism_gravity(
+            coordinates, prisms, densities, field=field, parallel=False
+        )
+        npt.assert_allclose(result_parallel, result_serial)