diff --git a/harmonica/_forward/point.py b/harmonica/_forward/point.py index 57876d956..861df1327 100644 --- a/harmonica/_forward/point.py +++ b/harmonica/_forward/point.py @@ -9,10 +9,22 @@ """ import numpy as np +from choclo.constants import GRAVITATIONAL_CONST +from choclo.point import ( + gravity_e, + gravity_ee, + gravity_en, + gravity_eu, + gravity_n, + gravity_nn, + gravity_nu, + gravity_pot, + gravity_u, + gravity_uu, +) from numba import jit, prange -from ..constants import GRAVITATIONAL_CONST -from .utils import check_coordinate_system, distance_cartesian, distance_spherical_core +from .utils import check_coordinate_system, distance_spherical_core def point_gravity( @@ -27,10 +39,91 @@ def point_gravity( r""" Compute gravitational fields of point masses. - It can compute the gravitational fields of point masses on a set of - computation points defined either in Cartesian or geocentric spherical + Compute the gravitational potential, gravitational acceleration and tensor + components generated by a collection of point masses on a set of + observation points defined either in Cartesian or geocentric spherical coordinates. + .. warning:: + + The **vertical direction points upwards**, i.e. positive and negative + values of ``upward`` represent points above and below the surface, + respectively. But ``g_z`` field returns the **downward component** of + the gravitational acceleration so that positive density contrasts + produce positive anomalies. The same applies to the tensor components, + i.e. the ``g_ez`` is the non-diagonal **easting-downward** tensor + component. + + .. important:: + + - The gravitational potential is returned in + :math:`\text{J}/\text{kg}`. + - The gravity acceleration components are returned in mGal + (:math:`\text{m}/\text{s}^2`). + - The tensor components are returned in Eotvos (:math:`\text{s}^{-2}`). + + Parameters + ---------- + coordinates : list of arrays + List of arrays containing the coordinates of computation points in the + following order: ``easting``, ``northing`` and ``upward`` (if + coordinates given in Cartesian coordinates), or ``longitude``, + ``latitude`` and ``radius`` (if given on a spherical geocentric + coordinate system). + All ``easting``, ``northing`` and ``upward`` should be in meters. + Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` + in meters. + points : list or array + List or array containing the coordinates of the point masses in the + following order: ``easting``, ``northing`` and ``upward`` (if + coordinates given in Cartesian coordinates), or ``longitude``, + ``latitude`` and ``radius`` (if given on a spherical geocentric + coordinate system). + All ``easting``, ``northing`` and ``upward`` should be in meters. + Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` + in meters. + masses : list or array + List or array containing the mass of each point mass in kg. + field : str + Gravitational field that wants to be computed. + The available fields coordinates are: + + - Gravitational potential: ``potential`` + - Easting acceleration: ``g_e`` + - Northing acceleration: ``g_n`` + - Downward acceleration: ``g_z`` + - Tensor components: + - ``g_ee`` + - ``g_nn`` + - ``g_zz`` + - ``g_en`` + - ``g_ez`` + - ``g_nz`` + + coordinate_system : str (optional) + Coordinate system of the coordinates of the computation points and the + 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``. + + Returns + ------- + result : array + Gravitational field generated by the ``point_mass`` on the computation + points defined in ``coordinates``. + The potential is given in SI units, the accelerations in mGal and the + Marussi tensor components in Eotvos. + + Notes + ----- The gravitational potential field generated by a point mass with mass :math:`m` located at a point :math:`Q` on a computation point :math:`P` can be computed as: @@ -133,66 +226,6 @@ def point_gravity( positive if the acceleration vector points inside the spheroid. - Parameters - ---------- - coordinates : list of arrays - List of arrays containing the coordinates of computation points in the - following order: ``easting``, ``northing`` and ``upward`` (if - coordinates given in Cartesian coordinates), or ``longitude``, - ``latitude`` and ``radius`` (if given on a spherical geocentric - coordinate system). - All ``easting``, ``northing`` and ``upward`` should be in meters. - Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` - in meters. - points : list or array - List or array containing the coordinates of the point masses in the - following order: ``easting``, ``northing`` and ``upward`` (if - coordinates given in Cartesian coordinates), or ``longitude``, - ``latitude`` and ``radius`` (if given on a spherical geocentric - coordinate system). - All ``easting``, ``northing`` and ``upward`` should be in meters. - Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` - in meters. - masses : list or array - List or array containing the mass of each point mass in kg. - field : str - Gravitational field that wants to be computed. - The available fields coordinates are: - - - Gravitational potential: ``potential`` - - Downward acceleration: ``g_z`` - - Northing acceleration: ``g_northing`` - - Easting acceleration: ``g_easting`` - - Tensor components: - - ``g_ee`` - - ``g_nn`` - - ``g_zz`` - - ``g_en`` - - ``g_ez`` - - ``g_nz`` - - coordinate_system : str (optional) - Coordinate system of the coordinates of the computation points and the - 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``. - - - Returns - ------- - result : array - Gravitational field generated by the ``point_mass`` on the computation - points defined in ``coordinates``. - The potential is given in SI units, the accelerations in mGal and the - Marussi tensor components in Eotvos. """ # Sanity checks for coordinate_system check_coordinate_system( @@ -216,9 +249,11 @@ def point_gravity( dispatcher(coordinate_system, parallel)( *coordinates, *points, masses, result, kernel ) - result *= GRAVITATIONAL_CONST + # Invert sign of gravity_u, gravity_eu, gravity_nu + if field in ("g_z", "g_ez", "g_ze", "g_nz", "g_zn"): + result *= -1 # Convert to more convenient units - if field in ("g_easting", "g_northing", "g_z"): + if field in ("g_e", "g_n", "g_z"): result *= 1e5 # SI to mGal tensors = ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz", "g_ne", "g_ze", "g_zn") if field in tensors: @@ -249,27 +284,27 @@ def get_kernel(coordinate_system, field): """ kernels = { "cartesian": { - "potential": kernel_potential_cartesian, - "g_z": kernel_g_z_cartesian, - "g_northing": kernel_g_northing_cartesian, - "g_easting": kernel_g_easting_cartesian, + "potential": gravity_pot, + "g_e": gravity_e, + "g_n": gravity_n, + "g_z": gravity_u, # diagonal tensor components - "g_ee": kernel_g_ee_cartesian, - "g_nn": kernel_g_nn_cartesian, - "g_zz": kernel_g_zz_cartesian, + "g_ee": gravity_ee, + "g_nn": gravity_nn, + "g_zz": gravity_uu, # non-diagonal tensor components - "g_en": kernel_g_en_cartesian, - "g_ez": kernel_g_ez_cartesian, - "g_nz": kernel_g_nz_cartesian, - "g_ne": kernel_g_en_cartesian, - "g_ze": kernel_g_ez_cartesian, - "g_zn": kernel_g_nz_cartesian, + "g_en": gravity_en, + "g_ez": gravity_eu, + "g_nz": gravity_nu, + "g_ne": gravity_en, + "g_ze": gravity_eu, + "g_zn": gravity_nu, }, "spherical": { - "potential": kernel_potential_spherical, - "g_z": kernel_g_z_spherical, - "g_northing": None, - "g_easting": None, + "potential": potential_spherical, + "g_z": gravity_u_spherical, + "g_n": None, + "g_e": None, }, } if field not in kernels[coordinate_system]: @@ -281,165 +316,12 @@ def get_kernel(coordinate_system, field): # ------------------------------------------ -# Kernel functions for Cartesian coordinates +# Kernel functions for Spherical coordinates # ------------------------------------------ @jit(nopython=True) -def kernel_potential_cartesian( - easting, northing, upward, easting_p, northing_p, upward_p -): - """ - Kernel function for gravitational potential field in Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return 1 / distance - - -# Acceleration components -# ----------------------- - - -@jit(nopython=True) -def kernel_g_z_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel for downward component of gravitational acceleration - - Use Cartesian coords. - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - # Remember that the ``g_z`` field returns the downward component of the - # gravitational acceleration. As a consequence, it is multiplied by -1. - # Notice that the ``g_z`` does not have the minus signal observed at the - # compoents ``g_northing`` and ``g_easting``. - return (upward - upward_p) / distance**3 - - -@jit(nopython=True) -def kernel_g_northing_cartesian( - easting, northing, upward, easting_p, northing_p, upward_p -): - """ - Kernel function for northing component of gravitational acceleration - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return -(northing - northing_p) / distance**3 - - -@jit(nopython=True) -def kernel_g_easting_cartesian( - easting, northing, upward, easting_p, northing_p, upward_p -): - """ - Kernel function for easting component of gravitational acceleration - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return -(easting - easting_p) / distance**3 - - -# Tensor components -# ----------------- - - -@jit(nopython=True) -def kernel_g_ee_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel function for g_ee component of gravitational tensor - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return 3 * (easting - easting_p) ** 2 / distance**5 - 1 / distance**3 - - -@jit(nopython=True) -def kernel_g_nn_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel function for g_nn component of gravitational tensor - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return 3 * (northing - northing_p) ** 2 / distance**5 - 1 / distance**3 - - -@jit(nopython=True) -def kernel_g_zz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel function for g_zz component of gravitational tensor - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return 3 * (upward - upward_p) ** 2 / distance**5 - 1 / distance**3 - - -@jit(nopython=True) -def kernel_g_en_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel function for g_en component of gravitational tensor - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - return 3 * (easting - easting_p) * (northing - northing_p) / distance**5 - - -@jit(nopython=True) -def kernel_g_ez_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel function for g_ez component of gravitational tensor - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - # Add a minus sign to account that the z axis points downwards. - return -3 * (easting - easting_p) * (upward - upward_p) / distance**5 - - -@jit(nopython=True) -def kernel_g_nz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p): - """ - Kernel function for g_nz component of gravitational tensor - - Use Cartesian coordinates - """ - distance = distance_cartesian( - (easting, northing, upward), (easting_p, northing_p, upward_p) - ) - # Add a minus sign to account that the z axis points downwards. - return -3 * (northing - northing_p) * (upward - upward_p) / distance**5 - - -# ------------------------------------------ -# Kernel functions for Cartesian coordinates -# ------------------------------------------ - - -@jit(nopython=True) -def kernel_potential_spherical( +def potential_spherical( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ): """ @@ -448,7 +330,7 @@ def kernel_potential_spherical( distance, _, _ = distance_spherical_core( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ) - return 1 / distance + return 1 / distance * GRAVITATIONAL_CONST # Acceleration components @@ -456,11 +338,11 @@ def kernel_potential_spherical( @jit(nopython=True) -def kernel_g_z_spherical( +def gravity_u_spherical( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ): """ - Kernel for downward component of gravitational acceleration + Kernel for upward component of gravitational acceleration Use spherical coordinates """ @@ -468,11 +350,19 @@ def kernel_g_z_spherical( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ) delta_z = radius - radius_p * cospsi - return delta_z / distance**3 + return -GRAVITATIONAL_CONST * delta_z / distance**3 def point_mass_cartesian( - easting, northing, upward, easting_p, northing_p, upward_p, masses, out, kernel + easting, + northing, + upward, + easting_p, + northing_p, + upward_p, + masses, + out, + forward_func, ): """ Compute gravitational field of point masses in Cartesian coordinates @@ -489,19 +379,21 @@ def point_mass_cartesian( 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. + forward_func : func + forward_func function that will be used to compute the gravitational + field on the computation points. It could be one of the forward + modelling functions in :mod:`choclo.point`. """ for l in prange(easting.size): for m in range(easting_p.size): - out[l] += masses[m] * kernel( + out[l] += forward_func( easting[l], northing[l], upward[l], easting_p[m], northing_p[m], upward_p[m], + masses[m], ) diff --git a/harmonica/_forward/tesseroid.py b/harmonica/_forward/tesseroid.py index eb6a7de00..c86d82e53 100644 --- a/harmonica/_forward/tesseroid.py +++ b/harmonica/_forward/tesseroid.py @@ -10,7 +10,6 @@ import numpy as np from numba import jit, prange -from ..constants import GRAVITATIONAL_CONST from ._tesseroid_utils import ( _adaptive_discretization, _check_tesseroids, @@ -23,7 +22,7 @@ density_based_discretization, gauss_legendre_quadrature_variable_density, ) -from .point import kernel_g_z_spherical, kernel_potential_spherical +from .point import gravity_u_spherical, potential_spherical from .utils import initialize_progressbar STACK_SIZE = 100 @@ -157,7 +156,10 @@ def tesseroid_gravity( ... ) """ - kernels = {"potential": kernel_potential_spherical, "g_z": kernel_g_z_spherical} + kernels = { + "potential": potential_spherical, + "g_z": gravity_u_spherical, + } if field not in kernels: raise ValueError("Gravitational field {} not recognized".format(field)) # Figure out the shape and size of the output array @@ -201,9 +203,11 @@ def tesseroid_gravity( dtype, progress_proxy, ) - result *= GRAVITATIONAL_CONST + # Invert sign + if field in ("g_z",): + result *= -1 # Convert to more convenient units - if field == "g_z": + if field in ("g_z",): result *= 1e5 # SI to mGal return result.reshape(cast.shape) diff --git a/harmonica/tests/test_point_gravity.py b/harmonica/tests/test_point_gravity.py index 560b61d83..172e84448 100644 --- a/harmonica/tests/test_point_gravity.py +++ b/harmonica/tests/test_point_gravity.py @@ -14,6 +14,18 @@ import numpy.testing as npt import pytest import verde as vd +from choclo.point import ( + gravity_e, + gravity_ee, + gravity_en, + gravity_eu, + gravity_n, + gravity_nn, + gravity_nu, + gravity_pot, + gravity_u, + gravity_uu, +) from .._forward.point import point_gravity from .._forward.utils import distance_cartesian @@ -47,7 +59,7 @@ def test_not_implemented_field(): point_mass = [0.0, 0.0, 0.0] mass = 1.0 coordinate_system = "spherical" - for field in ("g_northing", "g_easting"): + for field in ("g_n", "g_e"): with pytest.raises(NotImplementedError): point_gravity( coordinates, @@ -159,7 +171,7 @@ def test_potential_symmetry_cartesian(): @pytest.mark.use_numba -@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z")) +@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z")) def test_acceleration_symmetry_cartesian(field): """ Test if the acceleration components verify the expected symmetry @@ -175,10 +187,10 @@ def test_acceleration_symmetry_cartesian(field): easting = point_mass[0] * np.ones(2) northing = point_mass[1] * np.ones(2) upward = point_mass[2] * np.ones(2) - if field == "g_northing": + if field == "g_n": northing[0] += distance northing[1] -= distance - elif field == "g_easting": + elif field == "g_e": easting[0] += distance easting[1] -= distance elif field == "g_z": @@ -204,8 +216,8 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05) mass : float Mass of the point source. field : str - Acceleration component that needs to be approximated ("g_easting", - "g_northing", "g_z"). + Acceleration component that needs to be approximated ("g_e", + "g_n", "g_z"). delta : float Distance use to compute the finite difference in meters. @@ -219,9 +231,9 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05) # Build a two computation points slightly shifted from the original # computation point by a small delta coordinates_pair = tuple([coord, coord] for coord in coordinates) - if field == "g_easting": + if field == "g_e": index = 0 - elif field == "g_northing": + elif field == "g_n": index = 1 elif field == "g_z": index = 2 @@ -247,7 +259,7 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05) @pytest.mark.use_numba -@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z")) +@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z")) @pytest.mark.parametrize( "coordinates, point, mass", ( @@ -271,7 +283,7 @@ def test_acceleration_finite_diff_cartesian(coordinates, point, mass, field): @pytest.mark.use_numba -@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z")) +@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z")) def test_acceleration_sign(field): """ Test if acceleration components have the correct sign @@ -281,9 +293,9 @@ def test_acceleration_sign(field): mass = [2670] # Define computation points coordinates = [np.zeros(3) for i in range(3)] - if field == "g_easting": + if field == "g_e": coordinates[0] = np.array([-150.7, -10, 79]) - elif field == "g_northing": + elif field == "g_n": coordinates[1] = np.array([0, 100.2, 210.7]) elif field == "g_z": coordinates[2] = np.array([100.11, -300.7, -400]) @@ -300,7 +312,7 @@ def test_acceleration_sign(field): "field", # fmt: off ( - "potential", "g_z", "g_northing", "g_easting", + "potential", "g_z", "g_n", "g_e", "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz" ), # fmt: on @@ -549,20 +561,12 @@ def tensor_finite_differences(coordinates, point, mass, field, delta=0.05): # Determine the direction along which the finite difference will be # computed direction_i, direction_j = field[-2], field[-1] - if direction_i == "e": - direction_i = "easting" - if direction_i == "n": - direction_i = "northing" - if direction_j == "e": - direction_j = "easting" - if direction_j == "n": - direction_j = "northing" # Build a two computation points slightly shifted from the original # computation point by a small delta coordinates_pair = tuple([coord, coord] for coord in coordinates) - if direction_j == "easting": + if direction_j == "e": index = 0 - elif direction_j == "northing": + elif direction_j == "n": index = 1 elif direction_j == "z": index = 2 @@ -763,3 +767,86 @@ def test_point_mass_spherical_parallel(): parallel=True, ) npt.assert_allclose(result_serial, result_parallel) + + +class TestAgainstChoclo: + """ + Test forward modelling functions against dumb Choclo runs + """ + + @pytest.fixture() + def sample_points(self): + """ + Return three sample prisms + """ + points = np.array( + [ + [-5, 5, -5, 5], + [-5, -5, 5, 5], + [-10, -10, -10, -5], + ], + dtype=float, + ) + masses = np.array([100.0, -100.0, 200.0, -300.0], dtype=float) + return points, masses + + @pytest.fixture() + def sample_coordinates(self): + """ + Return four sample observation points + """ + easting = np.array([-5, 10, 0, 15], dtype=float) + northing = np.array([14, -4, 11, 0], dtype=float) + upward = np.array([9, 6, 6, 12], dtype=float) + return (easting, northing, upward) + + @pytest.mark.use_numba + @pytest.mark.parametrize( + "field, choclo_func", + [ + ("potential", gravity_pot), + ("g_e", gravity_e), + ("g_n", gravity_n), + ("g_z", gravity_u), + ("g_ee", gravity_ee), + ("g_nn", gravity_nn), + ("g_zz", gravity_uu), + ("g_en", gravity_en), + ("g_ez", gravity_eu), + ("g_nz", gravity_nu), + ], + ) + def test_against_choclo( + self, + field, + choclo_func, + sample_coordinates, + sample_points, + ): + """ + Tests forward functions against dumb runs on Choclo + """ + easting, northing, upward = sample_coordinates + points, masses = sample_points + # Compute expected results with dumb choclo calls + expected_result = np.zeros_like(easting) + for i in range(easting.size): + for j in range(masses.size): + expected_result[i] += choclo_func( + easting[i], + northing[i], + upward[i], + points[0, j], + points[1, j], + points[2, j], + masses[j], + ) + if field in ("g_z", "g_ez", "g_ze", "g_nz", "g_zn"): + expected_result *= -1 # invert sign + if field in ("g_e", "g_n", "g_z"): + expected_result *= 1e5 # convert to mGal + if field in ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"): + expected_result *= 1e9 # convert to Eotvos + # Compare with Harmonica results + result = point_gravity(sample_coordinates, points, masses, field=field) + npt.assert_allclose(result, expected_result) diff --git a/harmonica/tests/test_transformations.py b/harmonica/tests/test_transformations.py index 2906e05ad..16c5b845b 100644 --- a/harmonica/tests/test_transformations.py +++ b/harmonica/tests/test_transformations.py @@ -124,7 +124,7 @@ def fixture_sample_g_n(sample_grid_coords, sample_sources): Return g_n field of sample points on sample grid coords """ points, masses = sample_sources - g_n = point_gravity(sample_grid_coords, points, masses, field="g_northing") + g_n = point_gravity(sample_grid_coords, points, masses, field="g_n") g_n = vd.make_xarray_grid( sample_grid_coords, g_n, @@ -140,7 +140,7 @@ def fixture_sample_g_e(sample_grid_coords, sample_sources): Return g_e field of sample points on sample grid coords """ points, masses = sample_sources - g_e = point_gravity(sample_grid_coords, points, masses, field="g_easting") + g_e = point_gravity(sample_grid_coords, points, masses, field="g_e") g_e = vd.make_xarray_grid( sample_grid_coords, g_e,