Test g_z of prisms against infinite slab

Add test that compares the analytical solution of the g_z generated by
an infinite slab with the one generated by a large enough prism. Also
tests convergence: checks if the result gets closer to the analytical
solution if the horizontal dimensions of the prisms are increased.

harmonica/tests/test_prism.py

@@ -5,6 +5,7 @@
 import numpy as np
 import numpy.testing as npt
 
+from ..constants import GRAVITATIONAL_CONST
 from ..forward.prism import (
  prism_gravity,
  _check_prisms,
@@ -138,3 +139,55 @@ def test_custom_log():
  x = np.linspace(1, 100, 101)
  for x_i in x:
  npt.assert_allclose(log(x_i), np.log(x_i))
+
+
+def gravity_infinite_slab(thickness, density):
+ """
+ Compute the downward gravity acceleration generated by a infinite slab
+
+ Parameters
+ ----------
+ thickness : float
+ Thickness of the infinite slab in meters.
+ density : float
+ Density of the infinite slab in kg/m^3.
+
+ Returns
+ -------
+ result : float
+ Downward component of the gravity acceleration generated by the infinite slab
+ in mGal.
+ """
+ result = 2 * np.pi * GRAVITATIONAL_CONST * density * thickness
+ # Convert the result for g_z to mGal
+ result *= 1e5
+ return result
+
+
+@pytest.mark.use_numba
+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
+ thickness = 10.5
+ sizes = np.logspace(3, 9, 7)
+ bottom, top = -thickness, 0
+ density = 2670
+ # Compute the gravity fields generated by each prism
+ results = np.zeros_like(sizes)
+ for i, size in enumerate(sizes):
+ west, east = -size / 2, size / 2
+ south, north = west, east
+ prism = [west, east, south, north, bottom, top]
+ results[i] = prism_gravity(coordinates, prism, density, field="g_z")
+ # Check convergence: assert if as the prism size increases, the result gets closer
+ # to the analytical solution for an infinite slab
+ analytical = gravity_infinite_slab(thickness, density)
+ diffs = abs((results - analytical) / analytical)
+ assert (diffs[1:] < diffs[:-1]).all()
+ # Check if the largest size is close enough to the analytical solution
+ npt.assert_allclose(analytical, results[-1])