diff --git a/src/harmonica/__init__.py b/src/harmonica/__init__.py
index 6881b9cc0..191a5eca1 100644
--- a/src/harmonica/__init__.py
+++ b/src/harmonica/__init__.py
@@ -18,9 +18,12 @@
 from ._forward.dipole import dipole_magnetic
 from ._forward.ellipsoids import Ellipsoid, ellipsoid_gravity, ellipsoid_magnetic
 from ._forward.point import point_gravity
-from ._forward.prism_gravity import prism_gravity
-from ._forward.prism_layer import DatasetAccessorPrismLayer, prism_layer
-from ._forward.prism_magnetic import prism_magnetic
+from ._forward.prisms import (
+    DatasetAccessorPrismLayer,
+    prism_gravity,
+    prism_layer,
+    prism_magnetic,
+)
 from ._forward.tesseroid_gravity import tesseroid_gravity
 from ._forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer
 from ._gravity_corrections import bouguer_correction
diff --git a/src/harmonica/_forward/prisms/__init__.py b/src/harmonica/_forward/prisms/__init__.py
new file mode 100644
index 000000000..af872175b
--- /dev/null
+++ b/src/harmonica/_forward/prisms/__init__.py
@@ -0,0 +1,13 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Forward modelling of rectangular prisms.
+"""
+
+from .gravity import prism_gravity
+from .layer import DatasetAccessorPrismLayer, prism_layer
+from .magnetic import prism_magnetic
diff --git a/src/harmonica/_forward/prism_gravity.py b/src/harmonica/_forward/prisms/gravity.py
similarity index 93%
rename from src/harmonica/_forward/prism_gravity.py
rename to src/harmonica/_forward/prisms/gravity.py
index fe98794e9..ee64393d8 100644
--- a/src/harmonica/_forward/prism_gravity.py
+++ b/src/harmonica/_forward/prisms/gravity.py
@@ -30,7 +30,8 @@
 )
 from numba import jit, prange
 
-from .utils import initialize_progressbar
+from ..utils import initialize_progressbar
+from .utils import check_prisms
 
 # Define dictionary with available gravity fields for prisms
 FIELDS = {
@@ -209,7 +210,7 @@ def prism_gravity(
                 f"Number of elements in density ({density.size}) "
                 + f"mismatch the number of prisms ({prisms.shape[0]})"
             )
-        _check_prisms(prisms)
+        check_prisms(prisms)
         _check_singular_points(coordinates, prisms, field)
     # Discard null prisms (zero volume or zero density)
     prisms, density = _discard_null_prisms(prisms, density)
@@ -448,40 +449,6 @@ def _any_singular_point_g_nz(coordinates, prisms):
     return False
 
 
-def _check_prisms(prisms):
-    """
-    Check if prisms boundaries are well defined.
-
-    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.
-    """
-    west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6))
-    err_msg = "Invalid prism or prisms. "
-    bad_we = west > east
-    bad_sn = south > north
-    bad_bt = bottom > top
-    if bad_we.any():
-        err_msg += "The west boundary can't be greater than the east one.\n"
-        for prism in prisms[bad_we]:
-            err_msg += f"\tInvalid prism: {prism}\n"
-        raise ValueError(err_msg)
-    if bad_sn.any():
-        err_msg += "The south boundary can't be greater than the north one.\n"
-        for prism in prisms[bad_sn]:
-            err_msg += f"\tInvalid prism: {prism}\n"
-        raise ValueError(err_msg)
-    if bad_bt.any():
-        err_msg += "The bottom boundary can't be greater than the top one.\n"
-        for prism in prisms[bad_bt]:
-            err_msg += f"\tInvalid tesseroid: {prism}\n"
-        raise ValueError(err_msg)
-
-
 def _discard_null_prisms(prisms, density):
     """
     Discard prisms with zero volume or zero density.
@@ -494,7 +461,7 @@ def _discard_null_prisms(prisms, density):
         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.
+        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.
diff --git a/src/harmonica/_forward/prism_layer.py b/src/harmonica/_forward/prisms/layer.py
similarity index 72%
rename from src/harmonica/_forward/prism_layer.py
rename to src/harmonica/_forward/prisms/layer.py
index 975960b50..fa72ab0c2 100644
--- a/src/harmonica/_forward/prism_layer.py
+++ b/src/harmonica/_forward/prisms/layer.py
@@ -10,12 +10,14 @@
 
 import warnings
 
+import numba
 import numpy as np
 import verde as vd
 import xarray as xr
 
-from ..visualization import prism_to_pyvista
-from .prism_gravity import prism_gravity
+from ...visualization import prism_to_pyvista
+from ..utils import initialize_progressbar
+from .gravity import FIELDS
 
 
 def prism_layer(
@@ -313,23 +315,23 @@ def gravity(
         self,
         coordinates,
         field,
-        progressbar=False,
+        *,
         density_name="density",
         thickness_threshold=None,
-        **kwargs,
+        parallel=True,
+        progressbar=False,
     ):
         r"""
         Compute the gravity generated by the layer of prisms.
 
-        Uses :func:`harmonica.prism_gravity` for computing the gravity field
-        generated by the prisms of the layer.
+        Compute the gravity field generated by the layer of prisms on a set of
+        observation points.
         The density of the prisms will be assigned from the ``data_var`` chosen
         through the ``density_name`` argument.
         Ignores the prisms which ``top`` or ``bottom`` boundaries are
         ``np.nan``'s.
         Prisms thinner than a given threshold can be optionally ignored through
         the ``thickness_threshold`` argument.
-        All ``kwargs`` will be passed to :func:`harmonica.prism_gravity`.
 
         Parameters
         ----------
@@ -346,18 +348,20 @@ def gravity(
             - Downward acceleration: ``g_z``
             - Diagonal tensor components: ``g_ee``, ``g_nn``, ``g_zz``
             - Non-diagonal tensor components: ``g_en``, ``g_ez``, ``g_nz``
-        progressbar : bool (optional)
-            If True, a progress bar of the computation will be printed to
-            standard error (stderr). Requires :mod:`numba_progress` to be
-            installed. Default to ``False``.
-        density_name : str (optional)
+        density_name : str, optional
             Name of the property layer (or ``data_var`` of the
             :class:`xarray.Dataset`) that will be used for the density of each
             prism in the layer. Default to ``"density"``
-        thickness_threshold : float or None
+        thickness_threshold : float or None, optional
             Prisms thinner than this threshold will be ignored in the
             forward gravity calculation. If None, every prism with non-zero
             volume will be considered. Default to None.
+        parallel : bool, optional
+            Whether to run the computation of gravity fields in parallel or in serial.
+        progressbar : bool, optional
+            If True, a progress bar of the computation will be printed to
+            standard error (stderr). Requires :mod:`numba_progress` to be
+            installed. Default to ``False``.
 
         Returns
         -------
@@ -369,68 +373,59 @@ def gravity(
         --------
         harmonica.prism_gravity
         """
-        # Get boundaries and density of the prisms
-        boundaries = self._to_prisms()
+        # Sanity check the field
+        if field not in FIELDS:
+            msg = f"Gravitational field '{field}' not recognized."
+            raise ValueError(msg)
+
+        # Check if prism layer defines a regular grid in horizontal coords
+        _check_regular_grid(self._obj.easting.values, self._obj.northing.values)
+
+        # Ravel coordinates to 1D
+        cast = np.broadcast(*coordinates[:3])
+        coordinates = tuple(np.atleast_1d(c).ravel() for c in coordinates[:3])
+
+        # Determine parallel or serial forward modelling function
+        numba_function = (
+            _forward_gravity_prism_layer_parallel
+            if parallel
+            else _forward_gravity_prism_layer_serial
+        )
+
+        # Forward model the prism layer
+        thickness_threshold = (
+            0.0 if thickness_threshold is None else thickness_threshold
+        )
         density = self._obj[density_name].values
-        # Get the mask for selecting only the prisms whose top boundary, bottom
-        # boundary and density have no nans
-        mask = self._get_nonans_mask(property_name=density_name)
-        # Select only the boundaries and density elements for masked prisms
-        boundaries = boundaries[mask.ravel()]
-        density = density[mask]
-        # Discard thin prisms and their densities
-        if thickness_threshold is not None:
-            boundaries, density = _discard_thin_prisms(
-                boundaries,
+        if np.isnan(density).any():
+            msg = (
+                "Found NaN values in 'density' property of the prisms layer. "
+                "Their respective prisms will be ignored."
+            )
+            warnings.warn(msg, stacklevel=2)
+        with initialize_progressbar(coordinates[0].size, progressbar) as progress_proxy:
+            result = numba_function(
+                coordinates,
+                self._obj.easting.values,
+                self._obj.northing.values,
+                self._obj.bottom.values,
+                self._obj.top.values,
                 density,
+                FIELDS[field],
                 thickness_threshold,
+                progress_proxy,
             )
-        # Return gravity field of prisms
-        return prism_gravity(
-            coordinates,
-            prisms=boundaries,
-            density=density,
-            field=field,
-            progressbar=progressbar,
-            **kwargs,
-        )
-
-    def _get_nonans_mask(self, property_name=None):
-        """
-        Build a mask for prisms with no nans on top, bottom or a property.
 
-        Parameters
-        ----------
-        property_name : str (optional)
-            Name of the property layer (or ``data_var`` of the
-            :class:`xarray.Dataset`) that will be used for masking the prisms
-            in the layer.
-
-        Returns
-        -------
-        mask : 2d-array
-            Array of bools that can be used as a mask for selecting prisms with
-            no nans on top boundaries, bottom boundaries and the passed
-            property.
-        """
-        # Mask the prisms that contains no nans on top and bottom boundaries
-        mask = np.logical_and(
-            np.logical_not(np.isnan(self._obj.top.values)),
-            np.logical_not(np.isnan(self._obj.bottom.values)),
-        )
-        # Mask the prisms that contains nans on the selected property
-        if property_name is not None:
-            mask_property = np.logical_not(np.isnan(self._obj[property_name].values))
-            # Warn if a nan is found within the masked property
-            if not mask_property[mask].all():
-                warnings.warn(
-                    f"Found missing values in '{property_name}' property "
-                    + "of the prisms layer. The prisms with a nan as "
-                    + f"'{property_name}' will be ignored.",
-                    stacklevel=1,
-                )
-            mask = np.logical_and(mask, mask_property)
-        return mask
+        # 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 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 _to_prisms(self):
         """
@@ -519,40 +514,113 @@ def to_pyvista(self, drop_null_prisms=True):
         return prism_to_pyvista(prisms, properties=properties)
 
 
-def _discard_thin_prisms(
-    prisms,
-    density,
-    thickness_threshold,
+def _forward_gravity_prism_layer(
+    coordinates,
+    prisms_easting,
+    prisms_northing,
+    prisms_bottom,
+    prisms_top,
+    densities,
+    forward_func,
+    thickness_threshold=0.0,
+    progress_proxy=None,
 ):
     """
-    Discard prisms with a thickness below a threshold.
+    Forward model the gravity fields of prisms in a prism layer.
+
+    Builds the boundaries of each one the prisms in the layer on the fly, iterating over
+    the easting and northing dimensional coordinates. This function is intended to avoid
+    allocating all prism boundaries in memory before computation.
+
+    The function ignores prisms with zero density, prisms thinner than the
+    ``thickness_threshold``, and prisms with NaN as top or bottom boundaries.
+
+    .. important::
+
+        This function is intended to be decorated with Numba before use.
 
     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.
-    density : 1d-array
-        Array containing the density of each prism in kg/m^3. Must have the
-        same size as the number of prisms.
+    coordinates : tuple of (n,) array
+        Tuple with the easting, northing, and upward coordinates of the observation
+        points.
+    prisms_easting : (n_e,) array
+        1D array with the easting coordinates of the center of the prisms.
+        Use the easting dimensional coordinate of the prism layer.
+    prisms_northing : (n_n,) array
+        1D array with the northing coordinates of the center of the prisms.
+        Use the northing dimensional coordinate of the prism layer.
+    prisms_bottom, prisms_top : (n_n, n_e) array
+        2D arrays with the bottom and top boundaries of the prisms.
+    densities : (n_n, n_e) array
+        2D array with the densities of the prisms in kg/m3.
+    forward_func : callable
+        Choclo function to forward model a gravity field of prisms.
     thickness_threshold : float
-        Prisms thinner than this threshold will be discarded.
+        Prisms thinner than this threshold will be ignored in the
+        forward gravity calculation. If None, every prism with non-zero
+        volume will be considered.
+    progress_proxy : :class:`numba_progress.ProgressBar` or None
+        Instance of :class:`numba_progress.ProgressBar` that gets updated after
+        each iteration on the observation points. Use None if no progress bar
+        is should be used.
 
     Returns
     -------
-    prisms : 2d-array
-        A copy of the ``prisms`` array that doesn't include the thin prisms.
-    density : 1d-array
-        A copy of the ``density`` array that doesn't include the density values
-        for thin prisms.
+    result : array
+        Gravity field in SI units.
     """
-    bottom, top = prisms[:, -2], prisms[:, -1]
-    # Mark prisms with thickness < threshold  as null prisms
-    thickness = top - bottom
-    null_prisms = thickness < thickness_threshold
-    # Keep only thick prisms and their densities
-    prisms = prisms[np.logical_not(null_prisms), :]
-    density = density[np.logical_not(null_prisms)]
-    return prisms, density
+    easting, northing, upward = coordinates
+    n_coords = easting.size
+    half_spacing_east = (prisms_easting[1] - prisms_easting[0]) / 2
+    half_spacing_north = (prisms_northing[1] - prisms_northing[0]) / 2
+
+    # Check if we need to update the progressbar on each iteration
+    update_progressbar = progress_proxy is not None
+
+    result = np.zeros(n_coords, dtype=np.float64)
+    for i in numba.prange(n_coords):
+        for j, easting_center in enumerate(prisms_easting):
+            west = easting_center - half_spacing_east
+            east = easting_center + half_spacing_east
+            for k, northing_center in enumerate(prisms_northing):
+                # Ignore prisms with zero or NaN density
+                density = densities[k, j]
+                if density == 0.0 or np.isnan(density):
+                    continue
+
+                # Ignore thin prisms
+                bottom, top = prisms_bottom[k, j], prisms_top[k, j]
+                if top - bottom < thickness_threshold:
+                    continue
+
+                # Ignore prisms with nan top or bottom
+                if np.isnan(top) or np.isnan(bottom):
+                    continue
+
+                south = northing_center - half_spacing_north
+                north = northing_center + half_spacing_north
+                result[i] += forward_func(
+                    easting[i],
+                    northing[i],
+                    upward[i],
+                    west,
+                    east,
+                    south,
+                    north,
+                    bottom,
+                    top,
+                    density,
+                )
+        # Update progress bar if called
+        if update_progressbar:
+            progress_proxy.update(1)
+    return result
+
+
+_forward_gravity_prism_layer_parallel = numba.jit(nopython=True, parallel=True)(
+    _forward_gravity_prism_layer
+)
+_forward_gravity_prism_layer_serial = numba.jit(nopython=True, parallel=False)(
+    _forward_gravity_prism_layer
+)
diff --git a/src/harmonica/_forward/prism_magnetic.py b/src/harmonica/_forward/prisms/magnetic.py
similarity index 99%
rename from src/harmonica/_forward/prism_magnetic.py
rename to src/harmonica/_forward/prisms/magnetic.py
index 65b139389..43fc5e344 100644
--- a/src/harmonica/_forward/prism_magnetic.py
+++ b/src/harmonica/_forward/prisms/magnetic.py
@@ -14,8 +14,8 @@
 from choclo.prism import magnetic_e, magnetic_field, magnetic_n, magnetic_u
 from numba import jit, prange
 
-from .prism_gravity import _check_prisms
-from .utils import initialize_progressbar
+from ..utils import initialize_progressbar
+from .utils import check_prisms
 
 VALID_FIELDS = ("b", "b_e", "b_n", "b_u")
 FORWARD_FUNCTIONS = {
@@ -412,7 +412,7 @@ def _discard_null_prisms(prisms, magnetization):
         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.
+        Run ``check_prisms`` before.
     magnetization : 2d-array
         Array containing the magnetization vector of each prism in
         :math:`Am^{-1}`. Each vector will be a row in the 2d-array.
@@ -454,7 +454,7 @@ def _run_sanity_checks(prisms, magnetization):
             f"Number of magnetization vectors ({magnetization[0].size}) "
             + f"mismatch the number of prisms ({prisms.shape[0]})"
         )
-    _check_prisms(prisms)
+    check_prisms(prisms)
 
 
 # Define jitted versions of the forward modelling function
diff --git a/src/harmonica/_forward/prisms/utils.py b/src/harmonica/_forward/prisms/utils.py
new file mode 100644
index 000000000..9a2468295
--- /dev/null
+++ b/src/harmonica/_forward/prisms/utils.py
@@ -0,0 +1,43 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Utilities for forward modelling of rectangular prisms.
+"""
+
+
+def check_prisms(prisms):
+    """
+    Check if prisms boundaries are well defined.
+
+    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.
+    """
+    west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6))
+    err_msg = "Invalid prism or prisms. "
+    bad_we = west > east
+    bad_sn = south > north
+    bad_bt = bottom > top
+    if bad_we.any():
+        err_msg += "The west boundary can't be greater than the east one.\n"
+        for prism in prisms[bad_we]:
+            err_msg += f"\tInvalid prism: {prism}\n"
+        raise ValueError(err_msg)
+    if bad_sn.any():
+        err_msg += "The south boundary can't be greater than the north one.\n"
+        for prism in prisms[bad_sn]:
+            err_msg += f"\tInvalid prism: {prism}\n"
+        raise ValueError(err_msg)
+    if bad_bt.any():
+        err_msg += "The bottom boundary can't be greater than the top one.\n"
+        for prism in prisms[bad_bt]:
+            err_msg += f"\tInvalid tesseroid: {prism}\n"
+        raise ValueError(err_msg)
diff --git a/src/harmonica/_forward/utils.py b/src/harmonica/_forward/utils.py
index 1831a5f63..aa1b46264 100644
--- a/src/harmonica/_forward/utils.py
+++ b/src/harmonica/_forward/utils.py
@@ -349,7 +349,7 @@ def initialize_progressbar(total, use_progressbar):
         Number of iterations that will be assumed as the total for the
         :class:`numba_progress.ProgressBar`.
     use_progressbar : bool
-        Weather to initialize a progressbar or not. If True, then the function
+        Whether to initialize a progressbar or not. If True, then the function
         will return None.
 
     Returns
diff --git a/test/test_prism.py b/test/test_prism.py
index 077e1470f..daeb950ac 100644
--- a/test/test_prism.py
+++ b/test/test_prism.py
@@ -33,12 +33,9 @@
 except ImportError:  # pragma: no cover
     ProgressBar = None
 
-from harmonica import bouguer_correction
-from harmonica._forward.prism_gravity import (
-    _check_prisms,
-    _discard_null_prisms,
-    prism_gravity,
-)
+from harmonica import bouguer_correction, prism_gravity
+from harmonica._forward.prisms.gravity import _discard_null_prisms
+from harmonica._forward.prisms.utils import check_prisms
 
 from .utils import run_only_with_numba
 
@@ -75,24 +72,24 @@ def test_invalid_density_array():
 
 
 def test_invalid_prisms():
-    """Check if invalid prism boundaries are caught by _check_prisms."""
+    """Check if invalid prism boundaries are caught by check_prisms."""
     w, e, s, n, bottom, top = -100, 100, -100, 100, -200, -100
     # Check if it works properly on valid prisms
-    _check_prisms(np.atleast_2d([w, e, s, n, bottom, top]))
+    check_prisms(np.atleast_2d([w, e, s, n, bottom, top]))
     # Check if it works properly on valid prisms with zero volume
-    _check_prisms(np.atleast_2d([w, w, s, n, bottom, top]))
-    _check_prisms(np.atleast_2d([w, e, s, s, bottom, top]))
-    _check_prisms(np.atleast_2d([w, e, s, n, bottom, bottom]))
+    check_prisms(np.atleast_2d([w, w, s, n, bottom, top]))
+    check_prisms(np.atleast_2d([w, e, s, s, bottom, top]))
+    check_prisms(np.atleast_2d([w, e, s, n, bottom, bottom]))
     # Test invalid boundaries
     msg = "The west boundary can't be greater than the east one"
     with pytest.raises(ValueError, match=msg):
-        _check_prisms(np.atleast_2d([e, w, s, n, bottom, top]))
+        check_prisms(np.atleast_2d([e, w, s, n, bottom, top]))
     msg = "The south boundary can't be greater than the north one"
     with pytest.raises(ValueError, match=msg):
-        _check_prisms(np.atleast_2d([w, e, n, s, bottom, top]))
+        check_prisms(np.atleast_2d([w, e, n, s, bottom, top]))
     msg = "The bottom boundary can't be greater than the top one"
     with pytest.raises(ValueError, match=msg):
-        _check_prisms(np.atleast_2d([w, e, s, n, top, bottom]))
+        check_prisms(np.atleast_2d([w, e, s, n, top, bottom]))
 
 
 def test_discard_null_prisms():
diff --git a/test/test_prism_layer.py b/test/test_prism_layer.py
index acd560dd0..e6f2ab939 100644
--- a/test/test_prism_layer.py
+++ b/test/test_prism_layer.py
@@ -8,9 +8,10 @@
 Test prisms layer.
 """
 
-import warnings
+import re
 from unittest.mock import patch
 
+import bordado as bd
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -18,7 +19,7 @@
 import xarray as xr
 
 from harmonica import prism_gravity, prism_layer
-from harmonica._forward.prism_layer import _discard_thin_prisms
+from harmonica._forward.prisms.gravity import FIELDS
 
 try:
     import pyvista
@@ -238,109 +239,8 @@ def test_prism_layer_get_prism_by_index():
             )
 
 
-def test_nonans_prisms_mask(dummy_layer):
-    """
-    Check if the mask for nonans prism is correctly created.
-    """
-    (easting, northing), surface, reference, _ = dummy_layer
-    shape = (northing.size, easting.size)
-
-    # No nan in top nor bottom
-    # ------------------------
-    layer = prism_layer((easting, northing), surface, reference)
-    expected_mask = np.ones(shape, dtype=bool)
-    mask = layer.prism_layer._get_nonans_mask()
-    npt.assert_allclose(mask, expected_mask)
-
-    # Nans in top only
-    # ----------------
-    layer = prism_layer((easting, northing), surface, reference)
-    expected_mask = np.ones(shape, dtype=bool)
-    for index in ((1, 2), (2, 3)):
-        layer.top[index] = np.nan
-        expected_mask[index] = False
-    mask = layer.prism_layer._get_nonans_mask()
-    npt.assert_allclose(mask, expected_mask)
-
-    # Nans in bottom only
-    # -------------------
-    layer = prism_layer((easting, northing), surface, reference)
-    expected_mask = np.ones(shape, dtype=bool)
-    for index in ((2, 1), (3, 2)):
-        layer.bottom[index] = np.nan
-        expected_mask[index] = False
-    mask = layer.prism_layer._get_nonans_mask()
-    npt.assert_allclose(mask, expected_mask)
-
-    # Nans in top and bottom
-    # ----------------------
-    layer = prism_layer((easting, northing), surface, reference)
-    expected_mask = np.ones(shape, dtype=bool)
-    for index in ((1, 2), (2, 3)):
-        layer.top[index] = np.nan
-        expected_mask[index] = False
-    for index in ((1, 2), (2, 1), (3, 2)):
-        layer.bottom[index] = np.nan
-        expected_mask[index] = False
-    mask = layer.prism_layer._get_nonans_mask()
-    npt.assert_allclose(mask, expected_mask)
-
-
-def test_nonans_prisms_mask_property(
-    dummy_layer,
-):
-    """
-    Check if the method masks the property and raises a warning.
-    """
-    (easting, northing), surface, reference, density = dummy_layer
-    shape = (northing.size, easting.size)
-
-    # Nans in top and property (on the same prisms)
-    # ---------------------------------------------
-    expected_mask = np.ones_like(surface, dtype=bool)
-    indices = ((1, 2), (2, 3))
-    # Set some elements of surface and density as nans
-    for index in indices:
-        surface[index] = np.nan
-        density[index] = np.nan
-        expected_mask[index] = False
-    layer = prism_layer(
-        (easting, northing), surface, reference, properties={"density": density}
-    )
-    # Check if no warning is raised
-    with warnings.catch_warnings(record=True) as warn:
-        mask = layer.prism_layer._get_nonans_mask(property_name="density")
-        assert len(warn) == 0
-    npt.assert_allclose(mask, expected_mask)
-
-    # Nans in top and property (not precisely on the same prisms)
-    # -----------------------------------------------------------
-    surface = np.arange(20, dtype=float).reshape(shape)
-    density = 2670 * np.ones_like(surface)
-    expected_mask = np.ones_like(surface, dtype=bool)
-    # Set some elements of surface as nans
-    indices = ((1, 2), (2, 3))
-    for index in indices:
-        surface[index] = np.nan
-        expected_mask[index] = False
-    # Set a different set of elements of density as nans
-    indices = ((2, 2), (0, 1))
-    for index in indices:
-        density[index] = np.nan
-        expected_mask[index] = False
-    layer = prism_layer(
-        (easting, northing), surface, reference, properties={"density": density}
-    )
-    # Check if warning is raised
-    with warnings.catch_warnings(record=True) as warn:
-        mask = layer.prism_layer._get_nonans_mask(property_name="density")
-        assert len(warn) == 1
-        assert issubclass(warn[-1].category, UserWarning)
-    npt.assert_allclose(mask, expected_mask)
-
-
 @pytest.mark.use_numba
-@pytest.mark.parametrize("field", ["potential", "g_z"])
+@pytest.mark.parametrize("field", FIELDS)
 def test_prism_layer_gravity(field, dummy_layer):
     """
     Check if gravity method works as expected.
@@ -361,6 +261,24 @@ def test_prism_layer_gravity(field, dummy_layer):
     )
 
 
+def test_prism_layer_invalid_field(dummy_layer):
+    """
+    Test error after passing invalid field.
+    """
+    coordinates = bd.grid_coordinates(
+        (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0
+    )
+    (easting, northing), surface, reference, density = dummy_layer
+    layer = prism_layer(
+        (easting, northing), surface, reference, properties={"density": density}
+    )
+
+    invalid_field = "invalid field"
+    msg = re.escape(f"Gravitational field '{invalid_field}' not recognized.")
+    with pytest.raises(ValueError, match=msg):
+        layer.prism_layer.gravity(coordinates, field=invalid_field)
+
+
 @pytest.mark.use_numba
 @pytest.mark.parametrize("field", ["potential", "g_z"])
 def test_prism_layer_gravity_surface_nans(field, dummy_layer, prism_layer_with_holes):
@@ -401,9 +319,9 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h
         prisms_coords, surface, reference, properties={"density": density}
     )
     # Check if warning is raised after passing density with nans
-    with warnings.catch_warnings(record=True) as warn:
+    msg = re.escape("Found NaN values in 'density' property of the prisms layer.")
+    with pytest.warns(match=msg):
         result = layer.prism_layer.gravity(coordinates, field=field)
-        assert len(warn) == 1
     # Check if it generates the expected gravity field
     prisms, rho = prism_layer_with_holes
     npt.assert_allclose(
@@ -544,42 +462,3 @@ def test_gravity_discarded_thin_prisms(dummy_layer):
         coordinates, field="g_z", thickness_threshold=5
     )
     npt.assert_allclose(gravity_manually_removed, gravity_threshold_removed)
-
-
-def test_discard_thin_prisms():
-    """
-    Check if thin prisms are properly discarded.
-    """
-    # create set of 4 prisms (west, east, south, north, bottom, top)
-    prism_boundaries = np.array(
-        [
-            [-5000.0, 5000.0, -5000.0, 5000.0, 0.0, 55.1],
-            [5000.0, 15000.0, -5000.0, 5000.0, 0.0, 55.01],
-            [-5000.0, 5000.0, 5000.0, 15000.0, 0.0, 35.0],
-            [5000.0, 15000.0, 5000.0, 15000.0, 0.0, 84.0],
-        ]
-    )
-
-    # assign densities to each prism
-    densities = np.array([2306, 2122, 2190, 2069])
-
-    # drop prisms and respective densities thinner than 55.05
-    # (2nd and 3rd prisms)
-    thick_prisms, thick_densities = _discard_thin_prisms(
-        prism_boundaries,
-        densities,
-        thickness_threshold=55.05,
-    )
-
-    # manually remove prisms and densities of thin prisms
-    expected_prisms = np.array(
-        [
-            [-5000.0, 5000.0, -5000.0, 5000.0, 0.0, 55.1],
-            [5000.0, 15000.0, 5000.0, 15000.0, 0.0, 84.0],
-        ]
-    )
-    expected_densities = np.array([2306, 2069])
-
-    # check the correct prisms and densities were discarded
-    npt.assert_allclose(expected_prisms, thick_prisms)
-    npt.assert_allclose(expected_densities, thick_densities)
diff --git a/test/test_prism_magnetic.py b/test/test_prism_magnetic.py
index 89f847ead..505ccf9e0 100644
--- a/test/test_prism_magnetic.py
+++ b/test/test_prism_magnetic.py
@@ -21,7 +21,7 @@
     ProgressBar = None
 
 from harmonica import prism_magnetic
-from harmonica._forward.prism_magnetic import VALID_FIELDS
+from harmonica._forward.prisms.magnetic import VALID_FIELDS
 
 from .utils import run_only_with_numba