From 8cea0f917e6515db095c19a9d184890a36d8bb6b Mon Sep 17 00:00:00 2001 From: Lu Li Date: Fri, 7 Jul 2023 23:29:34 +0800 Subject: [PATCH] Revert "Update point source use Choclo" This reverts commit 9bc2e0405ccdf1026d64c922369e4d09d1bd9053. --- harmonica/_forward/point.py | 235 ++++++++++++++++++++------ harmonica/_forward/tesseroid.py | 38 +++-- harmonica/tests/test_point_gravity.py | 121 ++----------- 3 files changed, 217 insertions(+), 177 deletions(-) diff --git a/harmonica/_forward/point.py b/harmonica/_forward/point.py index 2f9bf7bca..57876d956 100644 --- a/harmonica/_forward/point.py +++ b/harmonica/_forward/point.py @@ -9,22 +9,10 @@ """ import numpy as np -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_spherical_core +from .utils import check_coordinate_system, distance_cartesian, distance_spherical_core def point_gravity( @@ -228,11 +216,9 @@ def point_gravity( dispatcher(coordinate_system, parallel)( *coordinates, *points, masses, result, kernel ) - # Invert sign of gravity_u, gravity_eu, gravity_nu - if field in ("g_z", "g_ez", "g_nz"): - result *= -1 + result *= GRAVITATIONAL_CONST # Convert to more convenient units - if field in ("g_e", "g_n", "g_z"): + if field in ("g_easting", "g_northing", "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: @@ -263,27 +249,27 @@ def get_kernel(coordinate_system, field): """ kernels = { "cartesian": { - "potential": gravity_pot, - "g_z": gravity_u, - "g_n": gravity_n, - "g_e": gravity_e, + "potential": kernel_potential_cartesian, + "g_z": kernel_g_z_cartesian, + "g_northing": kernel_g_northing_cartesian, + "g_easting": kernel_g_easting_cartesian, # diagonal tensor components - "g_ee": gravity_ee, - "g_nn": gravity_nn, - "g_zz": gravity_uu, + "g_ee": kernel_g_ee_cartesian, + "g_nn": kernel_g_nn_cartesian, + "g_zz": kernel_g_zz_cartesian, # non-diagonal tensor components - "g_en": gravity_en, - "g_ez": gravity_eu, - "g_nz": gravity_nu, - "g_ne": gravity_en, - "g_ze": gravity_eu, - "g_zn": gravity_nu, + "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, }, "spherical": { "potential": kernel_potential_spherical, - "g_z": kernel_gravity_u_spherical, - "g_n": None, - "g_e": None, + "g_z": kernel_g_z_spherical, + "g_northing": None, + "g_easting": None, }, } if field not in kernels[coordinate_system]: @@ -295,7 +281,160 @@ def get_kernel(coordinate_system, field): # ------------------------------------------ -# Kernel functions for Spherical coordinates +# Kernel functions for Cartesian 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 # ------------------------------------------ @@ -309,7 +448,7 @@ def kernel_potential_spherical( distance, _, _ = distance_spherical_core( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ) - return 1 / distance * GRAVITATIONAL_CONST + return 1 / distance # Acceleration components @@ -317,11 +456,11 @@ def kernel_potential_spherical( @jit(nopython=True) -def kernel_gravity_u_spherical( +def kernel_g_z_spherical( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ): """ - Kernel for upward component of gravitational acceleration + Kernel for downward component of gravitational acceleration Use spherical coordinates """ @@ -329,19 +468,11 @@ def kernel_gravity_u_spherical( longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p ) delta_z = radius - radius_p * cospsi - return -GRAVITATIONAL_CONST * delta_z / distance**3 + return delta_z / distance**3 def point_mass_cartesian( - easting, - northing, - upward, - easting_p, - northing_p, - upward_p, - masses, - out, - forward_func, + easting, northing, upward, easting_p, northing_p, upward_p, masses, out, kernel ): """ Compute gravitational field of point masses in Cartesian coordinates @@ -358,21 +489,19 @@ 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``. - 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`. + kernel : func + Kernel function that will be used to compute the gravitational field on + the computation points. """ for l in prange(easting.size): for m in range(easting_p.size): - out[l] += forward_func( + out[l] += masses[m] * kernel( 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 ddfa175bc..4154f2e29 100644 --- a/harmonica/_forward/tesseroid.py +++ b/harmonica/_forward/tesseroid.py @@ -10,6 +10,7 @@ import numpy as np from numba import jit, prange +from ..constants import GRAVITATIONAL_CONST from ._tesseroid_utils import ( _adaptive_discretization, _check_points_outside_tesseroids, @@ -22,7 +23,7 @@ density_based_discretization, gauss_legendre_quadrature_variable_density, ) -from .point import kernel_gravity_u_spherical, kernel_potential_spherical +from .point import kernel_g_z_spherical, kernel_potential_spherical STACK_SIZE = 100 MAX_DISCRETIZATIONS = 100000 @@ -40,7 +41,7 @@ def tesseroid_gravity( dtype=np.float64, disable_checks=False, ): - """ + r""" Compute gravitational field of tesseroids on computation points. .. warning:: @@ -50,15 +51,21 @@ def tesseroid_gravity( It is equivalent to the opposite of the radial component, therefore it's positive if the acceleration vector points inside the spheroid. + .. 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`). + Parameters ---------- - coordinates : list or 1d-array - List or array containing ``longitude``, ``latitude`` and ``radius`` of - the computation points defined on a spherical geocentric coordinate - system. - Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` - in meters. - tesseroids : list or 1d-array + coordinates : list of arrays + List of arrays containing the ``longitude``, ``latitude`` and + ``radius`` coordinates of the computation points, defined on + a spherical geocentric coordinate system. Both ``longitude`` and + ``latitude`` should be in degrees and ``radius`` in meters. + tesseroids : list or 2d-array List or array containing the coordinates of the tesseroid: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric spherical coordinate system. @@ -102,7 +109,8 @@ def tesseroid_gravity( ------- result : array Gravitational field generated by the tesseroids on the computation - points. + points. Gravitational potential is returned in + :math:`\text{J}/\text{kg}` and acceleration components in mgal. References ---------- @@ -143,10 +151,7 @@ def tesseroid_gravity( ... ) """ - kernels = { - "potential": kernel_potential_spherical, - "g_z": kernel_gravity_u_spherical, - } + kernels = {"potential": kernel_potential_spherical, "g_z": kernel_g_z_spherical} if field not in kernels: raise ValueError("Gravitational field {} not recognized".format(field)) # Figure out the shape and size of the output array @@ -188,9 +193,10 @@ def tesseroid_gravity( kernels[field], dtype, ) - # Convert to more convenient units and Invert sign + result *= GRAVITATIONAL_CONST + # Convert to more convenient units if field == "g_z": - result *= -1e5 # SI to mGal + 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 e7913ecdc..9f710daea 100644 --- a/harmonica/tests/test_point_gravity.py +++ b/harmonica/tests/test_point_gravity.py @@ -14,18 +14,6 @@ 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 @@ -59,7 +47,7 @@ def test_not_implemented_field(): point_mass = [0.0, 0.0, 0.0] mass = 1.0 coordinate_system = "spherical" - for field in ("g_n", "g_e"): + for field in ("g_northing", "g_easting"): with pytest.raises(NotImplementedError): point_gravity( coordinates, @@ -171,7 +159,7 @@ def test_potential_symmetry_cartesian(): @pytest.mark.use_numba -@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z")) +@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z")) def test_acceleration_symmetry_cartesian(field): """ Test if the acceleration components verify the expected symmetry @@ -187,10 +175,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_n": + if field == "g_northing": northing[0] += distance northing[1] -= distance - elif field == "g_e": + elif field == "g_easting": easting[0] += distance easting[1] -= distance elif field == "g_z": @@ -216,8 +204,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_e", - "g_n", "g_z"). + Acceleration component that needs to be approximated ("g_easting", + "g_northing", "g_z"). delta : float Distance use to compute the finite difference in meters. @@ -231,9 +219,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_e": + if field == "g_easting": index = 0 - elif field == "g_n": + elif field == "g_northing": index = 1 elif field == "g_z": index = 2 @@ -259,7 +247,7 @@ def acceleration_finite_differences(coordinates, point, mass, field, delta=0.05) @pytest.mark.use_numba -@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z")) +@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z")) @pytest.mark.parametrize( "coordinates, point, mass", ( @@ -283,7 +271,7 @@ def test_acceleration_finite_diff_cartesian(coordinates, point, mass, field): @pytest.mark.use_numba -@pytest.mark.parametrize("field", ("g_n", "g_e", "g_z")) +@pytest.mark.parametrize("field", ("g_northing", "g_easting", "g_z")) def test_acceleration_sign(field): """ Test if acceleration components have the correct sign @@ -293,9 +281,9 @@ def test_acceleration_sign(field): mass = [2670] # Define computation points coordinates = [np.zeros(3) for i in range(3)] - if field == "g_e": + if field == "g_easting": coordinates[0] = np.array([-150.7, -10, 79]) - elif field == "g_n": + elif field == "g_northing": coordinates[1] = np.array([0, 100.2, 210.7]) elif field == "g_z": coordinates[2] = np.array([100.11, -300.7, -400]) @@ -312,7 +300,7 @@ def test_acceleration_sign(field): "field", # fmt: off ( - "potential", "g_z", "g_n", "g_e", + "potential", "g_z", "g_northing", "g_easting", "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz" ), # fmt: on @@ -776,86 +764,3 @@ 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[j, 0], - points[j, 1], - points[j, 2], - masses[j], - ) - if field in ("g_z", "g_ez", "g_nz"): - 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)