Solve numerical instabilities of almost spheroids (#611)

Solve the numerical instabilities in gravity and magnetic fields due to
prolate, oblate or triaxial ellipsoids that are close enough to a
sphere. Fallback to sphere analytic solutions when the ellipsoid's
semiaxes lengths are close enough within a relative tolerance.

**Relevant issues/PRs:**

Part of #594

Author: santisoler

harmonica/_forward/ellipsoids/gravity.py

@@ -17,7 +17,12 @@
 
 from ...errors import NoPhysicalPropertyWarning
 from ...typing import Coordinates, Ellipsoid
-from .utils import calculate_lambda, get_elliptical_integrals, is_internal
+from .utils import (
+    calculate_lambda,
+    get_elliptical_integrals,
+    is_almost_a_sphere,
+    is_internal,
+)
 
 
 def ellipsoid_gravity(
@@ -152,7 +157,7 @@ def _compute_gravity_ellipsoid(
     # Mask internal points
     internal = is_internal(x, y, z, a, b, c)
 
-    if a == b == c:
+    if is_almost_a_sphere(a, b, c):
         # Fallback to sphere equations which are simpler
         factor = -4 / 3 * np.pi * a**3 * GRAVITATIONAL_CONST * density
         gx, gy, gz = tuple(np.zeros_like(x) for _ in range(3))

harmonica/_forward/ellipsoids/magnetic.py

@@ -25,6 +25,7 @@
     calculate_lambda,
     get_derivatives_of_elliptical_integrals,
     get_elliptical_integrals,
+    is_almost_a_sphere,
     is_internal,
 )
 
@@ -359,14 +360,14 @@ def get_demagnetization_tensor_internal(a: float, b: float, c: float):
 
     where :math:`\mathbf{M}` is the magnetization vector of the ellipsoid.
     """
-    if a > b > c:
+    if is_almost_a_sphere(a, b, c):
+        n_diagonal = 1 / 3 * np.ones(3)
+    elif a > b > c:
         n_diagonal = _demag_tensor_triaxial_internal(a, b, c)
     elif a > b and b == c:
         n_diagonal = _demag_tensor_prolate_internal(a, b)
     elif a < b and b == c:
         n_diagonal = _demag_tensor_oblate_internal(a, b)
-    elif a == b == c:
-        n_diagonal = 1 / 3 * np.ones(3)
     else:
         msg = "Could not determine ellipsoid type for values given."
         raise ValueError(msg)
@@ -532,7 +533,7 @@ def get_demagnetization_tensor_external(
 
     coords = (x, y, z)
 
-    if a == b == c:
+    if is_almost_a_sphere(a, b, c):
         r = np.sqrt(x**2 + y**2 + z**2)
         r5, r3 = r**5, r**3
         factor = -(a**3 / 3)

harmonica/_forward/ellipsoids/utils.py

@@ -8,6 +8,11 @@
 import numpy.typing as npt
 from scipy.special import ellipeinc, ellipkinc
 
+# Relative tolerance for two ellipsoid semiaxes to be considered almost equal.
+# E.g.: two semiaxes a and b are considered almost equal if:
+# | a - b | <  max(a, b) * SEMIAXES_RTOL
+SEMIAXES_RTOL = 1e-5
+
 
 def is_internal(x, y, z, a, b, c):
     """
@@ -27,6 +32,41 @@ def is_internal(x, y, z, a, b, c):
     return ((x**2) / (a**2) + (y**2) / (b**2) + (z**2) / (c**2)) < 1
 
 
+def is_almost_a_sphere(a: float, b: float, c: float) -> bool:
+    """
+    Check if a given ellipsoid approximates a sphere.
+
+    Returns True if ellipsoid's semiaxes lenghts are close enough to each other to be
+    approximated by a sphere.
+
+    Parameters
+    ----------
+    a, b, c: float
+        Ellipsoid's semiaxes lenghts.
+
+    Returns
+    -------
+    bool
+    """
+    # Exactly a sphere
+    if a == b == c:
+        return True
+
+    # Prolate or oblate that is almost a sphere
+    if b == c and np.abs(a - b) < SEMIAXES_RTOL * max(a, b):
+        return True
+
+    # Triaxial that is almost a sphere
+    if (  # noqa: SIM103
+        a != b != c
+        and np.abs(a - b) < SEMIAXES_RTOL * max(a, b)
+        and np.abs(b - c) < SEMIAXES_RTOL * max(b, c)
+    ):
+        return True
+
+    return False
+
+
 def calculate_lambda(x, y, z, a, b, c):
     """
     Get the lambda ellipsoidal coordinate for a given ellipsoid and observation points.

harmonica/tests/ellipsoids/test_gravity.py

@@ -663,3 +663,109 @@ def test_warning(self, ellipsoid_class, ellipsoid_args):
         # Check the gravity acceleration components are zero
         for g_component in (gx, gy, gz):
             assert g_component == 0.0
+
+
+class TestNumericalInstability:
+    """
+    Test fix for numerical instabilities when semiaxes lengths are very similar.
+
+    When the semiaxes are almost equal to each other (in the order of machine
+    precision), analytic solutions for prolate, oblate, and triaxial ellipsoids might
+    fall under singularities, triggering numerical instabilities.
+
+    Given two semiaxes a and b, define a ratio as:
+
+    .. code::
+
+        ratio  = | a - b | / max(a, b)
+
+    These test functions will build the ellipsoids using very small values of this
+    ratio.
+    """
+
+    # Properties of the sphere
+    radius = 50.0
+    center = (0, 0, 0)
+    density = 200.0
+
+    def get_coordinates(self, azimuth=45, polar=45):
+        """
+        Generate coordinates of observation points along a certain direction.
+        """
+        r = np.linspace(self.radius, 5 * self.radius, 501)
+        azimuth, polar = np.rad2deg(azimuth), np.rad2deg(polar)
+        easting = r * np.cos(azimuth) * np.cos(polar)
+        northing = r * np.sin(azimuth) * np.cos(polar)
+        upward = r * np.sin(polar)
+        return (easting, northing, upward)
+
+    @pytest.fixture
+    def sphere(self):
+        return Sphere(self.radius, center=self.center, density=self.density)
+
+    def test_gravity_prolate(self, sphere):
+        """
+        Test gravity field of prolate ellipsoid that is almost a sphere.
+        """
+        coordinates = self.get_coordinates()
+
+        # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities
+        ratio = 1e-9
+
+        a = self.radius
+        b = (1 - ratio) * a
+        ellipsoid = ProlateEllipsoid(
+            a, b, yaw=0, pitch=0, center=self.center, density=self.density
+        )
+        ge_sphere, gn_sphere, gu_sphere = ellipsoid_gravity(coordinates, sphere)
+        ge_ell, gn_ell, gu_ell = ellipsoid_gravity(coordinates, ellipsoid)
+
+        rtol, atol = 1e-7, 1e-8
+        np.testing.assert_allclose(ge_ell, ge_sphere, rtol=rtol, atol=atol)
+        np.testing.assert_allclose(gn_ell, gn_sphere, rtol=rtol, atol=atol)
+        np.testing.assert_allclose(gu_ell, gu_sphere, rtol=rtol, atol=atol)
+
+    def test_gravity_oblate(self, sphere):
+        """
+        Test gravity field of oblate ellipsoid that is almost a sphere.
+        """
+        coordinates = self.get_coordinates()
+
+        # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities
+        ratio = 1e-14
+
+        a = self.radius
+        b = (1 + ratio) * a
+        ellipsoid = OblateEllipsoid(
+            a, b, yaw=0, pitch=0, center=self.center, density=self.density
+        )
+        ge_sphere, gn_sphere, gu_sphere = ellipsoid_gravity(coordinates, sphere)
+        ge_ell, gn_ell, gu_ell = ellipsoid_gravity(coordinates, ellipsoid)
+
+        rtol, atol = 1e-7, 1e-8
+        np.testing.assert_allclose(ge_ell, ge_sphere, rtol=rtol, atol=atol)
+        np.testing.assert_allclose(gn_ell, gn_sphere, rtol=rtol, atol=atol)
+        np.testing.assert_allclose(gu_ell, gu_sphere, rtol=rtol, atol=atol)
+
+    def test_gravity_triaxial(self, sphere):
+        """
+        Test gravity field of triaxial ellipsoid that is almost a sphere.
+        """
+        coordinates = self.get_coordinates()
+
+        # Semiaxes ratio. Sufficiently small to trigger the numerical instabilities
+        ratio = 1e-14
+
+        a = self.radius
+        b = (1 - ratio) * a
+        c = (1 - 2 * ratio) * a
+        ellipsoid = TriaxialEllipsoid(
+            a, b, c, yaw=0, pitch=0, roll=0, center=self.center, density=self.density
+        )
+        ge_sphere, gn_sphere, gu_sphere = ellipsoid_gravity(coordinates, sphere)
+        ge_ell, gn_ell, gu_ell = ellipsoid_gravity(coordinates, ellipsoid)
+
+        rtol, atol = 1e-7, 1e-8
+        np.testing.assert_allclose(ge_ell, ge_sphere, rtol=rtol, atol=atol)
+        np.testing.assert_allclose(gn_ell, gn_sphere, rtol=rtol, atol=atol)
+        np.testing.assert_allclose(gu_ell, gu_sphere, rtol=rtol, atol=atol)