From a202db2dd734c2800b058fc739809a33aca5c23c Mon Sep 17 00:00:00 2001 From: "Jean Pierre Mandujano G." Date: Thu, 23 Jul 2026 14:53:26 -0500 Subject: [PATCH] Add linearity tests for forward models Test that the field generated by a set of sources equals the sum of the fields of each source computed separately. Covers point_gravity (Cartesian and spherical), dipole_magnetic, prism_gravity, prism_magnetic and tesseroid_gravity, parametrized by field. --- test/test_linearity.py | 372 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 372 insertions(+) create mode 100644 test/test_linearity.py diff --git a/test/test_linearity.py b/test/test_linearity.py new file mode 100644 index 000000000..a91eed988 --- /dev/null +++ b/test/test_linearity.py @@ -0,0 +1,372 @@ +# 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) +# +""" +Test linearity of the forward models. + +The forward models are linear on the sources: the field generated by a set of +sources should be equal to the sum of the fields generated by each source +separately. +""" + +import bordado as bd +import boule +import numpy as np +import numpy.testing as npt +import pytest + +from harmonica import ( + dipole_magnetic, + point_gravity, + prism_gravity, + prism_magnetic, + tesseroid_gravity, +) +from harmonica._forward.dipole import VALID_FIELDS as DIPOLE_FIELDS +from harmonica._forward.prisms.magnetic import VALID_FIELDS as PRISM_MAGNETIC_FIELDS + +GRAVITY_FIELDS = ( + "potential", + "g_e", + "g_n", + "g_z", + "g_ee", + "g_nn", + "g_zz", + "g_en", + "g_ez", + "g_nz", +) + + +class TestPointGravityLinearity: + """ + Test linearity of ``point_gravity``. + """ + + @pytest.fixture() + def sample_points(self): + """ + Return four sample point masses in Cartesian coordinates. + """ + points = ( + np.array([-50.0, 0.0, 55.0, 10.0]), + np.array([20.0, 0.0, -35.0, 15.0]), + np.array([-100.0, -150.0, -50.0, -20.0]), + ) + return points + + @pytest.fixture() + def sample_masses(self): + """ + Return masses for the sample points, including negative ones. + """ + return np.array([2.0e3, -3.0e3, 5.0e3, -1.5e3]) + + @pytest.fixture() + def sample_coordinates(self): + """ + Return a grid of observation points above the sources. + """ + return bd.grid_coordinates( + region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10 + ) + + @pytest.mark.parametrize("field", GRAVITY_FIELDS) + def test_linearity(self, sample_coordinates, sample_points, sample_masses, field): + """ + Compare the field of all points vs. the sum of their individual fields. + """ + result = point_gravity(sample_coordinates, sample_points, sample_masses, field) + expected = sum( + point_gravity( + sample_coordinates, + tuple(p[i : i + 1] for p in sample_points), + sample_masses[i : i + 1], + field, + ) + for i in range(sample_masses.size) + ) + npt.assert_allclose(result, expected) + + @pytest.mark.parametrize("field", ["potential", "g_z"]) + def test_linearity_spherical(self, sample_masses, field): + """ + Compare joint vs. summed individual fields in spherical coordinates. + """ + radius = boule.WGS84.mean_radius + points = ( + np.array([-50.0, 0.0, 55.0, 10.0]), + np.array([20.0, 0.0, -35.0, 15.0]), + radius - np.array([1e3, 2e3, 5e3, 10e3]), + ) + coordinates = bd.grid_coordinates( + region=(-100, 100, -70, 70), spacing=20, non_dimensional_coords=radius + ) + result = point_gravity( + coordinates, points, sample_masses, field, coordinate_system="spherical" + ) + expected = sum( + point_gravity( + coordinates, + tuple(p[i : i + 1] for p in points), + sample_masses[i : i + 1], + field, + coordinate_system="spherical", + ) + for i in range(sample_masses.size) + ) + npt.assert_allclose(result, expected) + + +class TestDipoleMagneticLinearity: + """ + Test linearity of ``dipole_magnetic``. + """ + + @pytest.fixture() + def sample_dipoles(self): + """ + Return four sample dipoles. + """ + dipoles = ( + np.array([-10.0, 0.0, 10.0, 5.0]), + np.array([2.0, 0.0, -1.0, 1.5]), + np.array([-10.0, -15.0, -5.0, -2.0]), + ) + return dipoles + + @pytest.fixture() + def sample_magnetic_moments(self): + """ + Return sample magnetic moment vectors for the dipoles. + """ + magnetic_moments = ( + np.array([1.0, -1.0, 3.0, -1.5]), + np.array([1.0, -0.5, -1.0, 2.0]), + np.array([1.0, 2.0, -3.0, -2.0]), + ) + return magnetic_moments + + @pytest.fixture() + def sample_coordinates(self): + """ + Return a grid of observation points above the dipoles. + """ + return bd.grid_coordinates( + region=(-20, 20, -20, 20), spacing=4, non_dimensional_coords=10 + ) + + @pytest.mark.parametrize("field", DIPOLE_FIELDS) + def test_linearity( + self, sample_coordinates, sample_dipoles, sample_magnetic_moments, field + ): + """ + Compare the field of all dipoles vs. the sum of their individual fields. + """ + result = dipole_magnetic( + sample_coordinates, sample_dipoles, sample_magnetic_moments, field=field + ) + n_dipoles = sample_dipoles[0].size + expected = sum( + np.asarray( + dipole_magnetic( + sample_coordinates, + tuple(d[i : i + 1] for d in sample_dipoles), + tuple(m[i : i + 1] for m in sample_magnetic_moments), + field=field, + ) + ) + for i in range(n_dipoles) + ) + npt.assert_allclose(result, expected) + + +class TestPrismGravityLinearity: + """ + Test linearity of ``prism_gravity``. + """ + + @pytest.fixture() + def sample_prisms(self): + """ + Return four non-overlapping sample prisms. + """ + prisms = np.array( + [ + [-100.0, -50.0, -100.0, -50.0, -30.0, -10.0], + [-50.0, 50.0, 50.0, 100.0, -50.0, -20.0], + [50.0, 100.0, -100.0, 50.0, -20.0, -5.0], + [-100.0, 50.0, -50.0, 50.0, -60.0, -40.0], + ] + ) + return prisms + + @pytest.fixture() + def sample_densities(self): + """ + Return densities for the sample prisms, including negative ones. + """ + return np.array([2670.0, -1500.0, 3300.0, -2900.0]) + + @pytest.fixture() + def sample_coordinates(self): + """ + Return a grid of observation points above the prisms. + """ + return bd.grid_coordinates( + region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10 + ) + + @pytest.mark.parametrize("field", GRAVITY_FIELDS) + def test_linearity( + self, sample_coordinates, sample_prisms, sample_densities, field + ): + """ + Compare the field of all prisms vs. the sum of their individual fields. + """ + result = prism_gravity( + sample_coordinates, sample_prisms, sample_densities, field + ) + expected = sum( + prism_gravity( + sample_coordinates, + sample_prisms[i : i + 1], + sample_densities[i : i + 1], + field, + ) + for i in range(sample_densities.size) + ) + npt.assert_allclose(result, expected) + + +class TestPrismMagneticLinearity: + """ + Test linearity of ``prism_magnetic``. + """ + + @pytest.fixture() + def sample_prisms(self): + """ + Return four non-overlapping sample prisms. + """ + prisms = np.array( + [ + [-100.0, -50.0, -100.0, -50.0, -30.0, -10.0], + [-50.0, 50.0, 50.0, 100.0, -50.0, -20.0], + [50.0, 100.0, -100.0, 50.0, -20.0, -5.0], + [-100.0, 50.0, -50.0, 50.0, -60.0, -40.0], + ] + ) + return prisms + + @pytest.fixture() + def sample_magnetizations(self): + """ + Return sample magnetization vectors for the prisms. + """ + magnetizations = ( + np.array([1.0, -1.0, 3.0, -1.5]), + np.array([1.0, -0.5, -1.0, 2.0]), + np.array([1.0, 2.0, -3.0, -2.0]), + ) + return magnetizations + + @pytest.fixture() + def sample_coordinates(self): + """ + Return a grid of observation points above the prisms. + """ + return bd.grid_coordinates( + region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10 + ) + + @pytest.mark.parametrize("field", PRISM_MAGNETIC_FIELDS) + def test_linearity( + self, sample_coordinates, sample_prisms, sample_magnetizations, field + ): + """ + Compare the field of all prisms vs. the sum of their individual fields. + """ + result = prism_magnetic( + sample_coordinates, sample_prisms, sample_magnetizations, field=field + ) + expected = sum( + np.asarray( + prism_magnetic( + sample_coordinates, + sample_prisms[i : i + 1], + tuple(m[i : i + 1] for m in sample_magnetizations), + field=field, + ) + ) + for i in range(sample_prisms.shape[0]) + ) + npt.assert_allclose(result, expected) + + +class TestTesseroidGravityLinearity: + """ + Test linearity of ``tesseroid_gravity``. + """ + + @pytest.fixture() + def sample_tesseroids(self): + """ + Return three non-overlapping sample tesseroids. + """ + ellipsoid = boule.WGS84 + top = ellipsoid.mean_radius + bottom = top - 1e3 + tesseroids = np.array( + [ + [-10.0, 10.0, -10.0, 10.0, bottom, top], + [20.0, 30.0, 20.0, 30.0, bottom - 1e3, top - 500], + [-50.0, -40.0, -30.0, -20.0, bottom - 2e3, top], + ] + ) + return tesseroids + + @pytest.fixture() + def sample_densities(self): + """ + Return densities for the sample tesseroids, including a negative one. + """ + return np.array([2670.0, -1500.0, 3300.0]) + + @pytest.fixture() + def sample_coordinates(self): + """ + Return a grid of observation points above the tesseroids. + """ + return bd.grid_coordinates( + region=(-60, 40, -40, 40), + spacing=20, + non_dimensional_coords=boule.WGS84.mean_radius + 100, + ) + + @pytest.mark.use_numba + @pytest.mark.parametrize("field", ["potential", "g_z"]) + def test_linearity( + self, sample_coordinates, sample_tesseroids, sample_densities, field + ): + """ + Compare the field of all tesseroids vs. the sum of their individual + fields. + """ + result = tesseroid_gravity( + sample_coordinates, sample_tesseroids, sample_densities, field=field + ) + expected = sum( + tesseroid_gravity( + sample_coordinates, + sample_tesseroids[i : i + 1], + sample_densities[i : i + 1], + field=field, + ) + for i in range(sample_densities.size) + ) + npt.assert_allclose(result, expected)