diff --git a/doc/user_guide/forward_modelling/prism.rst b/doc/user_guide/forward_modelling/prism.rst index 1f6dfbc6b..9093b146c 100644 --- a/doc/user_guide/forward_modelling/prism.rst +++ b/doc/user_guide/forward_modelling/prism.rst @@ -48,6 +48,97 @@ computation point: print(potential, "J/kg") +Gravitational fields +^^^^^^^^^^^^^^^^^^^^ + +The :func:`harmonica.prism_gravity` is able to compute the gravitational +potential (``"potential"``), the acceleration components (``"g_e"``, ``"g_n"``, +``"g_z"``), and tensor components (``"g_ee"``, ``"g_nn"``, ``"g_zz"``, +``"g_en"``, ``"g_ez"``, ``"g_nz"``). + + +Build a regular grid of computation points located a 10m above the zero height: + +.. jupyter-execute:: + + import verde as vd + + region = (-10e3, 10e3, -10e3, 10e3) + shape = (51, 51) + height = 10 + coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height) + +Define a single prism: + +.. jupyter-execute:: + + prism = [-2e3, 2e3, -2e3, 2e3, -1.6e3, -900] + density = 3300 + +Compute the gravitational fields that this prism generate on each observation +point: + +.. jupyter-execute:: + + fields = ( + "potential", + "g_e", "g_n", "g_z", + "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz" + ) + + results = {} + for field in fields: + results[field] = hm.prism_gravity(coordinates, prism, density, field=field) + +Plot the results: + +.. jupyter-execute:: + + import matplotlib.pyplot as plt + + plt.pcolormesh(coordinates[0], coordinates[1], results["potential"]) + plt.gca().set_aspect("equal") + plt.gca().ticklabel_format(style="sci", scilimits=(0, 0)) + plt.colorbar(label="J/kg") + plt.show() + + +.. jupyter-execute:: + + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) + + for field, ax in zip(("g_e", "g_n", "g_z"), axes): + tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field]) + ax.set_aspect("equal") + ax.set_title(field) + ax.ticklabel_format(style="sci", scilimits=(0, 0)) + plt.colorbar(tmp, ax=ax, label="mGal", orientation="horizontal", pad=0.08) + plt.show() + +.. jupyter-execute:: + + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) + + for field, ax in zip(("g_ee", "g_nn", "g_zz"), axes): + tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field]) + ax.set_aspect("equal") + ax.set_title(field) + ax.ticklabel_format(style="sci", scilimits=(0, 0)) + plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08) + plt.show() + +.. jupyter-execute:: + + fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8)) + + for field, ax in zip(("g_en", "g_ez", "g_nz"), axes): + tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field]) + ax.set_aspect("equal") + ax.set_title(field) + ax.ticklabel_format(style="sci", scilimits=(0, 0)) + plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08) + plt.show() + Passing multiple prisms ^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/environment.yml b/environment.yml index e7c399cd9..79deff76c 100644 --- a/environment.yml +++ b/environment.yml @@ -18,6 +18,7 @@ dependencies: - verde>=1.7.0 - xarray>=0.16 - xrft>=1.0 + - choclo>=0.1 # Optional requirements - pyvista>=0.27 - vtk>=9 diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 7a4a3a712..3ac4013ec 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -7,17 +7,48 @@ """ Forward modelling for prisms """ +import warnings + import numpy as np +from choclo.prism import ( + gravity_e, + gravity_ee, + gravity_en, + gravity_eu, + gravity_n, + gravity_nn, + gravity_nu, + gravity_pot, + gravity_u, + gravity_uu, +) +from choclo.prism._utils import ( + is_point_on_easting_edge, + is_point_on_northing_edge, + is_point_on_upward_edge, +) from numba import jit, prange +# Define dictionary with available gravity fields for prisms +FIELDS = { + "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, +} + # Attempt to import numba_progress try: from numba_progress import ProgressBar except ImportError: ProgressBar = None -from ..constants import GRAVITATIONAL_CONST - def prism_gravity( coordinates, @@ -32,23 +63,18 @@ def prism_gravity( """ Gravitational fields of right-rectangular prisms in Cartesian coordinates - The gravitational fields are computed through the analytical solutions - given by [Nagy2000]_ and [Nagy2002]_, which are valid on the entire domain. - This means that the computation can be done at any point, either outside or - inside the prism. - - This implementation makes use of the modified arctangent function proposed - by [Fukushima2020]_ (eq. 12) so that the potential field to satisfies - Poisson's equation in the entire domain. Moreover, the logarithm function - was also modified in order to solve the singularities that the analytical - solution has on some points (see [Nagy2000]_). + Compute the gravitational potential, gravitational acceleration and tensor + components generated by a collection of prisms on a set of observation + points. .. warning:: - The **z direction points upwards**, i.e. positive and negative values - of ``upward`` represent points above and below the surface, - respectively. But remember that the ``g_z`` field returns the downward - component of the gravitational acceleration so that positive density - contrasts produce positive anomalies. + 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. Parameters ---------- @@ -58,11 +84,11 @@ def prism_gravity( system. All coordinates should be in meters. prisms : list, 1d-array, or 2d-array List or array containing the coordinates of the prism(s) in the - following order: - west, east, south, north, bottom, top in a Cartesian coordinate system. + following order: west, east, south, north, bottom, top, in a Cartesian + coordinate system. All coordinates should be in meters. Coordinates for more than one prism can be provided. In this case, *prisms* should be a list of lists - or 2d-array (with one prism per line). + or 2d-array (with one prism per row). density : list or array List or array containing the density of each prism in kg/m^3. field : str @@ -70,7 +96,11 @@ def prism_gravity( The available fields are: - Gravitational potential: ``potential`` + - Eastward acceleration: ``g_e`` + - Northward acceleration: ``g_n`` - Downward acceleration: ``g_z`` + - Diagonal tensor components: ``g_ee``, ``g_nn``, ``g_zz`` + - Non-diagonal tensor components: ``g_en``, ``g_ez``, ``g_nz`` parallel : bool (optional) If True the computations will run in parallel using Numba built-in @@ -95,6 +125,36 @@ def prism_gravity( result : array Gravitational field generated by the prisms on the computation points. + Notes + ----- + This function makes use of :mod:`choclo.prism` forward modelling functions + to compute each gravitational field. + + The gravitational potential (``"potential"``) and the acceleration + components (``"g_e"``, ``"g_n"`` and ``"g_z"``) are well defined on the + entire domain. + Tensor components aren't defined on prism vertices. + Diagonal tensor components aren't defined on edges normal to the direction + of the tensor (e.g. `"g_ee"` is not defined on edges parallel to northing + and upward directions). + Non-diagonal tensor components aren't defined on edges normal to the + remaining direction of the tensor (e.g. `"g_en"` is not defined on edges + parallel to the upward direction). + The function returns :func:`numpy.nan` on every singular point. + + The diagonal tensor components aren't defined on observation points that + belong to the faces normal to their direction (e.g. ``"g_zz"`` is not + define on horizontal faces): two different limits exist when approaching + from the inside and from the outside of the prism. + This function returns the limit of these components while approaching from + the outside. + + References + ---------- + * [Nagy2000]_ + * [Nagy2002]_ + * [Fukushima2020]_ + Examples -------- @@ -104,9 +164,12 @@ def prism_gravity( >>> prism = [-34, 5, -18, 14, -345, -146] >>> # Set prism density to 2670 kg/m³ >>> density = 2670 - >>> # Define three computation points along the easting axe at 30m above - >>> # the surface - >>> coordinates = ([-40, 0, 40], [0, 0, 0], [30, 30, 30]) + >>> # Define three computation points along the easting direction at 30m + >>> # above the surface + >>> easting = [-40, 0, 40] + >>> northing = [0, 0, 0] + >>> upward = [30, 30, 30] + >>> coordinates = (easting, northing, upward) >>> # Compute the downward component of the gravitational acceleration that >>> # the prism generates on the computation points >>> gz = prism_gravity(coordinates, prism, density, field="g_z") @@ -123,8 +186,7 @@ def prism_gravity( (-0.05380, 0.02908, 0.11237) """ - kernels = {"potential": kernel_potential, "g_z": kernel_g_z} - if field not in kernels: + if field not in FIELDS: raise ValueError("Gravitational field {} not recognized".format(field)) # Figure out the shape and size of the output array cast = np.broadcast(*coordinates[:3]) @@ -141,40 +203,252 @@ def prism_gravity( + "mismatch the number of prisms ({})".format(prisms.shape[0]) ) _check_prisms(prisms) + _check_singular_points(coordinates, prisms, field) # Discard null prisms (zero volume or zero density) prisms, density = _discard_null_prisms(prisms, density) # Show progress bar for 'jit_prism_gravity' function if progressbar: if ProgressBar is None: raise ImportError( - "Missing optional dependency 'numba_progress' required if progressbar=True" + "Missing optional dependency 'numba_progress' required " + "if progressbar=True" ) progress_proxy = ProgressBar(total=coordinates[0].size) else: progress_proxy = None + # Choose parallelized or serialized forward function + if parallel: + gravity_prism_func = jit_prism_gravity_parallel + else: + gravity_prism_func = jit_prism_gravity_serial # Compute gravitational field - dispatcher(parallel)( - coordinates, prisms, density, kernels[field], result, progress_proxy + gravity_prism_func( + coordinates, prisms, density, FIELDS[field], result, progress_proxy ) - result *= GRAVITATIONAL_CONST # Close previously created progress bars - if progressbar: + if progress_proxy: progress_proxy.close() + # Invert sign of gravity_u, gravity_eu, gravity_nu + if field in ("g_z", "g_ez", "g_nz"): + result *= -1 # Convert to more convenient units - if field == "g_z": + if field in ("g_e", "g_n", "g_z"): result *= 1e5 # SI to mGal + # Convert to more convenient units + if field in ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"): + result *= 1e9 # SI to Eotvos return result.reshape(cast.shape) -def dispatcher(parallel): +def _check_singular_points(coordinates, prisms, field): """ - Return the parallelized or serialized forward modelling function + Check if any observation point is a singular point for tensor components + + The analytic solutions for the tensor components of the prism have some + singular points: + + - All prism vertices are singular points for every tensor component. + - Diagonal components aren't defined on edges perpendicular to the + component direction (e.g. ``g_ee`` is not defined on edges parallel to + northing and upward directions). + - Non-diagonal components aren't defined on edges perpendicular to the + two directions of the component (e.g. ``g_en`` is not defined on edges + parallel to the upward direction). """ - dispatchers = { - True: jit_prism_gravity_parallel, - False: jit_prism_gravity_serial, + functions = { + "g_ee": _any_singular_point_g_ee, + "g_nn": _any_singular_point_g_nn, + "g_zz": _any_singular_point_g_zz, + "g_en": _any_singular_point_g_en, + "g_ez": _any_singular_point_g_ez, + "g_nz": _any_singular_point_g_nz, } - return dispatchers[parallel] + if field not in functions: + return None + if functions[field](coordinates, prisms): + warnings.warn( + "Found observation point on singular point of a prism.", UserWarning + ) + + +@jit(nopython=True) +def _any_singular_point_g_ee(coordinates, prisms): + """ + Check observation points as singular points of g_ee + """ + easting, northing, upward = coordinates + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if is_point_on_northing_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ) or is_point_on_upward_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ): + return True + return False + + +@jit(nopython=True) +def _any_singular_point_g_nn(coordinates, prisms): + """ + Check observation points as singular points of g_nn + """ + easting, northing, upward = coordinates + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if is_point_on_easting_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ) or is_point_on_upward_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ): + return True + return False + + +@jit(nopython=True) +def _any_singular_point_g_zz(coordinates, prisms): + """ + Check observation points as singular points of g_zz + """ + easting, northing, upward = coordinates + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if is_point_on_easting_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ) or is_point_on_northing_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ): + return True + return False + + +@jit(nopython=True) +def _any_singular_point_g_en(coordinates, prisms): + """ + Check observation points as singular points of g_en + """ + easting, northing, upward = coordinates + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if is_point_on_upward_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ): + return True + return False + + +@jit(nopython=True) +def _any_singular_point_g_ez(coordinates, prisms): + """ + Check observation points as singular points of g_ez + """ + easting, northing, upward = coordinates + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if is_point_on_northing_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ): + return True + return False + + +@jit(nopython=True) +def _any_singular_point_g_nz(coordinates, prisms): + """ + Check observation points as singular points of g_nz + """ + easting, northing, upward = coordinates + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if is_point_on_easting_edge( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + ): + return True + return False def _check_prisms(prisms): @@ -248,7 +522,9 @@ def _discard_null_prisms(prisms, density): return prisms, density -def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy=None): +def jit_prism_gravity( + coordinates, prisms, density, forward_func, out, progress_proxy=None +): """ Compute gravitational field of prisms on computations points @@ -266,109 +542,38 @@ def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy= density : 1d-array Array containing the density of each prism in kg/m^3. Must have the same size as the number of prisms. - kernel : func - Kernel function that will be used to compute the desired field. + forward_func : func + Forward modelling function that will be used to compute the desired + field. It could be one of the forward modelling functions in + :mod:`choclo.prism`. out : 1d-array Array where the resulting field values will be stored. Must have the same size as the arrays contained on ``coordinates``. """ + # Unpack coordinates + easting, northing, upward = coordinates # Check if we need to update the progressbar on each iteration update_progressbar = progress_proxy is not None # Iterate over computation points and prisms - for l in prange(coordinates[0].size): + for l in prange(easting.size): for m in range(prisms.shape[0]): - # Iterate over the prism vertices to compute the result of the - # integration (see Nagy et al., 2000) - for i in range(2): - # Compute shifted easting coordinate - shift_east = prisms[m, 1 - i] - coordinates[0][l] - for j in range(2): - # Compute shifted northing coordinate - shift_north = prisms[m, 3 - j] - coordinates[1][l] - for k in range(2): - # Compute shifted upward coordinate - shift_upward = prisms[m, 5 - k] - coordinates[2][l] - # If i, j or k is 1, the corresponding shifted - # coordinate will refer to the lower boundary, - # meaning the corresponding term should have a minus - # sign. - out[l] += ( - density[m] - * (-1) ** (i + j + k) - * kernel(shift_east, shift_north, shift_upward) - ) + out[l] += forward_func( + easting[l], + northing[l], + upward[l], + prisms[m, 0], + prisms[m, 1], + prisms[m, 2], + prisms[m, 3], + prisms[m, 4], + prisms[m, 5], + density[m], + ) # Update progress bar if called if update_progressbar: progress_proxy.update(1) -@jit(nopython=True) -def kernel_potential(easting, northing, upward): - """ - Kernel function for potential gravitational field generated by a prism - """ - radius = np.sqrt(easting**2 + northing**2 + upward**2) - kernel = ( - easting * northing * safe_log(upward + radius) - + northing * upward * safe_log(easting + radius) - + easting * upward * safe_log(northing + radius) - - 0.5 * easting**2 * safe_atan2(upward * northing, easting * radius) - - 0.5 * northing**2 * safe_atan2(upward * easting, northing * radius) - - 0.5 * upward**2 * safe_atan2(easting * northing, upward * radius) - ) - return kernel - - -@jit(nopython=True) -def kernel_g_z(easting, northing, upward): - """ - Kernel for downward component of gravitational acceleration of a prism - """ - radius = np.sqrt(easting**2 + northing**2 + upward**2) - kernel = ( - easting * safe_log(northing + radius) - + northing * safe_log(easting + radius) - - upward * safe_atan2(easting * northing, upward * radius) - ) - return kernel - - -@jit(nopython=True) -def safe_atan2(y, x): - """ - Principal value of the arctangent expressed as a two variable function - - This modification has to be made to the arctangent function so the - gravitational field of the prism satisfies the Poisson's equation. - Therefore, it guarantees that the fields satisfies the symmetry properties - of the prism. This modified function has been defined according to - [Fukushima2020]_. - """ - if x != 0: - result = np.arctan(y / x) - else: - if y > 0: - result = np.pi / 2 - elif y < 0: - result = -np.pi / 2 - else: - result = 0 - return result - - -@jit(nopython=True) -def safe_log(x): - """ - Modified log to return 0 for log(0). - The limits in the formula terms tend to 0 (see [Nagy2000]_). - """ - if np.abs(x) < 1e-10: - result = 0 - else: - result = np.log(x) - return result - - # Define jitted versions of the forward modelling function jit_prism_gravity_serial = jit(nopython=True)(jit_prism_gravity) jit_prism_gravity_parallel = jit(nopython=True, parallel=True)(jit_prism_gravity) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 214d159e3..8a13c94d4 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -13,6 +13,18 @@ import numpy.testing as npt import pytest import verde as vd +from choclo.prism import ( + gravity_e, + gravity_ee, + gravity_en, + gravity_eu, + gravity_n, + gravity_nn, + gravity_nu, + gravity_pot, + gravity_u, + gravity_uu, +) try: from numba_progress import ProgressBar @@ -20,13 +32,7 @@ ProgressBar = None from .. import bouguer_correction -from .._forward.prism import ( - _check_prisms, - _discard_null_prisms, - prism_gravity, - safe_atan2, - safe_log, -) +from .._forward.prism import _check_prisms, _discard_null_prisms, prism_gravity from .utils import run_only_with_numba @@ -41,7 +47,7 @@ def test_invalid_field(): def test_invalid_density_array(): "Check if error is raised when density shape does not match prisms shape" - # Create a set of 4 tesseroids + # Create a set of 4 prisms prisms = [ [-100, 0, -100, 0, -200, -100], [-100, 0, 0, 100, -200, -100], @@ -90,7 +96,7 @@ def test_discard_null_prisms(): ) density = np.array([2670, 0, 3300, 3000, 2800, 2700]) prisms, density = _discard_null_prisms(prisms, density) - npt.assert_allclose(prisms, np.array([[0, 10, -50, 33, -2e3, 150]])) + npt.assert_allclose(prisms, np.array([[0.0, 10.0, -50.0, 33.0, -2e3, 150.0]])) npt.assert_allclose(density, np.array([2670])) @@ -139,254 +145,112 @@ def test_disable_checks(): npt.assert_allclose(invalid_result, -valid_result) -@pytest.mark.use_numba -def test_potential_field_symmetry(): - "Test if the potential field satisfies symmetry" - prism = [-100, 100, -100, 100, -100, 100] - density = 2670 - # Create six outside computation points located on the normal directions to - # the prism faces and at the same distance from its center - coordinates = ( - [-200, 200, 0, 0, 0, 0], - [0, 0, -200, 200, 0, 0], - [0, 0, 0, 0, -200, 200], - ) - result = prism_gravity(coordinates, prism, density, field="potential") - npt.assert_allclose(result[0], result) - # Create six inside computation points located on the normal directions to - # the prism faces and at the same distance from its center - coordinates = ([-50, 50, 0, 0, 0, 0], [0, 0, -50, 50, 0, 0], [0, 0, 0, 0, -50, 50]) - result = prism_gravity(coordinates, prism, density, field="potential") - npt.assert_allclose(result[0], result) - # Create twelve outside computation points located on the diagonal - # directions to the prism faces and at the same distance from its center. - # They can be divided into three sets: one made by those points that live - # on the horizontal plane that passes through the prism center, and the - # other two that live on the pair of vertical and perpendicular planes that - # also passes through the center of the prism. - coordinates = ( - [-200, -200, 200, 200, -200, -200, 200, 200, 0, 0, 0, 0], - [-200, 200, -200, 200, 0, 0, 0, 0, -200, -200, 200, 200], - [0, 0, 0, 0, -200, 200, -200, 200, -200, 200, -200, 200], - ) - result = prism_gravity(coordinates, prism, density, field="potential") - npt.assert_allclose(result[0], result) - # Create the same twelve points as before, but now all points fall inside - # the prism - coordinates = ( - [-50, -50, 50, 50, -50, -50, 50, 50, 0, 0, 0, 0], - [-50, 50, -50, 50, 0, 0, 0, 0, -50, -50, 50, 50], - [0, 0, 0, 0, -50, 50, -50, 50, -50, 50, -50, 50], +class TestAgainstChoclo: + """ + Test forward modelling functions against dumb Choclo runs + """ + + @pytest.fixture() + def sample_prisms(self): + """ + Return three sample prisms + """ + prisms = np.array( + [ + [-10, 10, -10, 0, -10, 0], + [-10, 0, -10, 10, -20, -10], + [5, 15, 5, 15, -15, -5], + ], + dtype=float, + ) + densities = np.array([400, -200, 300], dtype=float) + return prisms, densities + + @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), + ], ) - result = prism_gravity(coordinates, prism, density, field="potential") - npt.assert_allclose(result[0], result) + def test_against_choclo( + self, + field, + choclo_func, + sample_coordinates, + sample_prisms, + ): + """ + Tests forward functions against dumb runs on Choclo + """ + easting, northing, upward = sample_coordinates + prisms, densities = sample_prisms + # Compute expected results with dumb choclo calls + expected_result = np.zeros_like(easting) + for i in range(easting.size): + for j in range(densities.size): + expected_result[i] += choclo_func( + easting[i], + northing[i], + upward[i], + prisms[j, 0], + prisms[j, 1], + prisms[j, 2], + prisms[j, 3], + prisms[j, 4], + prisms[j, 5], + densities[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 = prism_gravity(sample_coordinates, prisms, densities, field=field) + npt.assert_allclose(result, expected_result) -@pytest.mark.use_numba -def test_g_z_symmetry_outside(): - """ - Test if the g_z field satisfies symmetry - - In order to test if the computed g_z satisfies the symmetry of a square - prism we will define several set of computation points: - - A. Two points located on the vertical axis of the prism (``easting == 0`` - and ``northing == 0``), one above and one below the prism at the same - distance from its center. - B. Four points located on the ``upward == 0`` plane around the prism - distributed normally to its faces , i.e. only one of the horizontal - coordinates will be nonzero. - C. Same as points defined in B, but located on a plane above the prism. - D. Same as points defined in B, but located on a plane below the prism. - E. Same as points defined in B, but located on a plane slightly above the - ``upward == 0`` plane. - F. Same as points defined in B, but located on a plane slightly below the - ``upward == 0`` plane. - G. Four points located on the ``upward == 0`` plane around the prism - distributed on the diagonal directions , i.e. both horizontal - coordinates will be equal and nonzero. - H. Same as points defined in G, but located on an plane above the prism. - I. Same as points defined in G, but located on an plane below the prism. - J. Same as points defined in G, but located on a plane slightly above the - ``upward == 0`` plane. - K. Same as points defined in G, but located on a plane slightly below the - ``upward == 0`` plane. - - All computation points defined on the previous groups fall outside of the - prism. - - The g_z field for a square prism (the horizontal dimensions of the prism - are equal) must satisfy the following symmetry rules: - - - The g_z values on points A must be opposite. - - The g_z values on points B must be all zero. - - The g_z values on points C must be all equal. - - The g_z values on points D must be all equal. - - The g_z values on points E must be all equal. - - The g_z values on points F must be all equal. - - The g_z values on points C and D must be opposite. - - The g_z values on points E and F must be opposite. - - The g_z values on points G must be all zero. - - The g_z values on points H must be all equal. - - The g_z values on points I must be all equal. - - The g_z values on points J must be all equal. - - The g_z values on points K must be all equal. - - The g_z values on points H and I must be opposite. - - The g_z values on points J and K must be opposite. - """ - prism = [-100, 100, -100, 100, -150, 150] - density = 2670 - computation_points = { - "A": ([0, 0], [0, 0], [-200, 200]), - "B": ([-200, 200, 0, 0], [0, 0, -200, 200], [0, 0, 0, 0]), - "C": ([-200, 200, 0, 0], [0, 0, -200, 200], [200, 200, 200, 200]), - "D": ([-200, 200, 0, 0], [0, 0, -200, 200], [-200, -200, -200, -200]), - "E": ([-200, 200, 0, 0], [0, 0, -200, 200], [1, 1, 1, 1]), - "F": ([-200, 200, 0, 0], [0, 0, -200, 200], [-1, -1, -1, -1]), - "G": ([-200, 200, 0, 0], [0, 0, -200, 200], [0, 0, 0, 0]), - "H": ([-200, -200, 200, 200], [-200, 200, -200, 200], [200, 200, 200, 200]), - "I": ([-200, -200, 200, 200], [-200, 200, -200, 200], [-200, -200, -200, -200]), - "J": ([-200, -200, 200, 200], [-200, 200, -200, 200], [1, 1, 1, 1]), - "K": ([-200, -200, 200, 200], [-200, 200, -200, 200], [-1, -1, -1, -1]), - } - # Compute g_z on each set of points - results = {} - for group, coords in computation_points.items(): - results[group] = prism_gravity(coords, prism, density, field="g_z") - # Check symmetries - # Values on A must be opposite, and the value of g_z at the point above the - # prism must have the same sign as the density, while the one below should - # have the opposite - npt.assert_allclose(results["A"][0], -results["A"][1]) - npt.assert_allclose(np.sign(results["A"][0]), -np.sign(density)) - npt.assert_allclose(np.sign(results["A"][1]), np.sign(density)) - # Values on C, D, E, F, H, I, J, K must be all equal within each set - for group in ["C", "D", "E", "F", "H", "I", "J", "K"]: - npt.assert_allclose(results[group][0], results[group]) - # Values on B and G must be zero - for group in ["B", "G"]: - npt.assert_allclose(0, results[group]) - # Values on C and D, E and F, H and I, J and K must be opposite - # Moreover, the set of points that are above the prism must have the same - # sign as the density, while the ones below should have the opposite - for above, below in (("C", "D"), ("E", "F"), ("H", "I"), ("J", "K")): - npt.assert_allclose(results[above], -results[below]) - npt.assert_allclose(np.sign(results[above]), np.sign(density)) - npt.assert_allclose(np.sign(results[below]), -np.sign(density)) - - -def test_g_z_symmetry_inside(): +@run_only_with_numba +def test_laplace(): """ - Test g_z symmetry on computation points that fall inside the prism - - In order to test if the computed g_z satisfies the symmetry of a square - prism on computation points that fall inside the prism, we will define - several set of computation points: - - A. Two points located on the vertical axis of the prism (``easting == 0`` - and ``northing == 0``), one above and one below the center of prism, - but at the same distance from it. - B. Four points located on the ``upward == 0`` plane around the prism - distributed normally to its faces , i.e. only one of the horizontal - coordinates will be nonzero. - C. Same as points defined in B, but located on a plane above the - ``upward == 0`` plane. - D. Same as points defined in B, but located on a plane below the - ``upward == 0`` plane. - E. Four points located on the ``upward == 0`` plane around the prism - distributed on the diagonal directions , i.e. both horizontal - coordinates will be equal and nonzero. - F. Same as points defined in E, but located on a plane above the - ``upward == 0`` plane. - G. Same as points defined in E, but located on a plane below the - ``upward == 0`` plane. - - All computation points defined on the previous groups fall outside of the - prism. - - The g_z field for a square prism (the horizontal dimensions of the prism - are equal) must satisfy the following symmetry rules: - - - The g_z values on points A must be opposite. - - The g_z values on points B must be all zero. - - The g_z values on points C must be all equal. - - The g_z values on points D must be all equal. - - The g_z values on points C and D must be opposite. - - The g_z values on points E must be all zero. - - The g_z values on points F must be all equal. - - The g_z values on points G must be all equal. - - The g_z values on points F and G must be opposite. + Test if the diagonal components satisfy Laplace equation """ - prism = [-100, 100, -100, 100, -150, 150] - density = 2670 - computation_points = { - "A": ([0, 0], [0, 0], [-50, 50]), - "B": ([-50, 50, 0, 0], [0, 0, -50, 50], [0, 0, 0, 0]), - "C": ([-50, 50, 0, 0], [0, 0, -50, 50], [50, 50, 50, 50]), - "D": ([-50, 50, 0, 0], [0, 0, -50, 50], [-50, -50, -50, -50]), - "E": ([-50, -50, 50, 50], [-50, 50, -50, 50], [0, 0, 0, 0]), - "F": ([-50, -50, 50, 50], [-50, 50, -50, 50], [50, 50, 50, 50]), - "G": ([-50, -50, 50, 50], [-50, 50, -50, 50], [-50, -50, -50, -50]), + region = (-10e3, 10e3, -10e3, 10e3) + coords = vd.grid_coordinates(region, shape=(10, 10), extra_coords=300) + prisms = [ + [1e3, 7e3, -5e3, 2e3, -1e3, -500], + [-4e3, 1e3, 4e3, 10e3, -2e3, 200], + ] + densities = [2670.0, 2900.0] + diagonal_components = { + field: prism_gravity(coords, prisms, densities, field=field) + for field in ("g_ee", "g_nn", "g_zz") } - # Compute g_z on each set of points - results = {} - for group, coords in computation_points.items(): - results[group] = prism_gravity(coords, prism, density, field="g_z") - # Check symmetries - # Values on A must be opposite, and the value of g_z at the point above the - # center of the prism must have the same sign as the density, while the one - # below should have the opposite - npt.assert_allclose(results["A"][0], -results["A"][1]) - npt.assert_allclose(np.sign(results["A"][0]), -np.sign(density)) - npt.assert_allclose(np.sign(results["A"][1]), np.sign(density)) - # Values on C, D, F, G must be all equal within each set - for group in ["C", "D", "F", "G"]: - npt.assert_allclose(results[group][0], results[group]) - # Values on B and E must be zero - for group in ["B", "E"]: - npt.assert_allclose(0, results[group]) - # Values on C and D, F and G must be opposite - # Moreover, the set of points that are above the center of the prism must - # have the same sign as the density, while the ones below should have the - # opposite - for above, below in (("C", "D"), ("F", "G")): - npt.assert_allclose(results[above], -results[below]) - npt.assert_allclose(np.sign(results[above]), np.sign(density)) - npt.assert_allclose(np.sign(results[below]), -np.sign(density)) - - -@pytest.mark.use_numba -def test_safe_atan2(): - "Test the safe_atan2 function" - # Test safe_atan2 for one point per quadrant - # First quadrant - x, y = 1, 1 - npt.assert_allclose(safe_atan2(y, x), np.pi / 4) - # Second quadrant - x, y = -1, 1 - npt.assert_allclose(safe_atan2(y, x), -np.pi / 4) - # Third quadrant - x, y = -1, -1 - npt.assert_allclose(safe_atan2(y, x), np.pi / 4) - # Forth quadrant - x, y = 1, -1 - npt.assert_allclose(safe_atan2(y, x), -np.pi / 4) - # Test safe_atan2 if the denominator is equal to zero - npt.assert_allclose(safe_atan2(1, 0), np.pi / 2) - npt.assert_allclose(safe_atan2(-1, 0), -np.pi / 2) - # Test safe_atan2 if both numerator and denominator are equal to zero - assert safe_atan2(0, 0) == 0 - - -@pytest.mark.use_numba -def test_safe_log(): - "Test the safe_log function" - # Check if safe_log function satisfies safe_log(0) == 0 - assert safe_log(0) == 0 - # Check if safe_log behaves like the natural logarithm in case that x != 0 - x = np.linspace(1, 100, 101) - for x_i in x: - npt.assert_allclose(safe_log(x_i), np.log(x_i)) + npt.assert_allclose( + diagonal_components["g_ee"] + diagonal_components["g_nn"], + -diagonal_components["g_zz"], + ) @pytest.mark.use_numba @@ -394,15 +258,19 @@ 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 + # Define an observation point at 1.5m above zero + height = 1.5 + coordinates = (0, 0, height) + # Define prisms with a top surface located at the same height as the + # observation point. Each prisms will have a thickness of 10.5m, and + # horizontal dimensions from 1e3 to 1e9m, and density of 2670 kg/m^3. + top = height thickness = 10.5 + bottom = top - thickness sizes = np.logspace(3, 9, 7) - bottom, top = -thickness, 0 density = 2670 - # Compute the gravity fields generated by each prism + # Compute the gravity field generated by each prism + # (from smaller to larger) results = np.zeros_like(sizes) for i, size in enumerate(sizes): prism = [-size / 2, size / 2, -size / 2, size / 2, bottom, top] @@ -410,8 +278,8 @@ def test_prism_against_infinite_slab(): # Check convergence: assert if as the prism size increases, the result gets # closer to the analytical solution for an infinite slab analytical = bouguer_correction(np.array(thickness), density) - diffs = abs((results - analytical) / analytical) - assert (diffs[1:] < diffs[:-1]).all() + errors = abs(analytical - results) + assert (errors[1:] < errors[:-1]).all() # Check if the largest size is close enough to the analytical solution npt.assert_allclose(analytical, results[-1]) @@ -486,3 +354,128 @@ def test_numba_progress_missing_error(): prism_gravity( coordinates, prisms, densities, field="potential", progressbar=True ) + + +class TestSingularPoints: + """ + Tests tensor components on singular points of the prism + """ + + @pytest.fixture + def sample_prism(self): + """Return a sample prism""" + return np.array([-10.3, 5.4, 8.6, 14.3, -30.3, 2.4]) + + @pytest.fixture + def sample_density(self): + """Return a sample density for the sample prism""" + return np.array([2900.0]) + + def get_vertices(self, prism): + """ + Return the vertices of the prism as points + """ + easting, northing, upward = tuple( + c.ravel() for c in np.meshgrid(prism[:2], prism[2:4], prism[4:6]) + ) + return easting, northing, upward + + def get_easting_edges_center(self, prism): + """ + Return points on the center of prism edges parallel to easting + """ + easting_c = (prism[0] + prism[1]) / 2 + northing, upward = tuple(c.ravel() for c in np.meshgrid(prism[2:4], prism[4:6])) + easting = np.full_like(northing, easting_c) + return easting, northing, upward + + def get_northing_edges_center(self, prism): + """ + Return points on the center of prism edges parallel to northing + """ + northing_c = (prism[2] + prism[3]) / 2 + easting, upward = tuple(c.ravel() for c in np.meshgrid(prism[0:2], prism[4:6])) + northing = np.full_like(easting, northing_c) + return easting, northing, upward + + def get_upward_edges_center(self, prism): + """ + Return points on the center of prism edges parallel to upward + """ + upward_c = (prism[4] + prism[5]) / 2 + easting, northing = tuple( + c.ravel() for c in np.meshgrid(prism[0:2], prism[2:4]) + ) + upward = np.full_like(easting, upward_c) + return easting, northing, upward + + @pytest.mark.use_numba + @pytest.mark.parametrize("field", ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz")) + def test_on_vertices(self, sample_prism, sample_density, field): + """ + Test tensor components when observation points fall on prism vertices + """ + easting, northing, upward = self.get_vertices(sample_prism) + for i in range(easting.size): + msg = "Found observation point" + with pytest.warns(UserWarning, match=msg): + prism_gravity( + (easting[i], northing[i], upward[i]), + sample_prism, + sample_density, + field=field, + ) + + @pytest.mark.use_numba + @pytest.mark.parametrize("field", ("g_nn", "g_zz", "g_nz")) + def test_on_easting_edges(self, sample_prism, sample_density, field): + """ + Test tensor components that have singular points on edges parallel to + easting direction + """ + easting, northing, upward = self.get_easting_edges_center(sample_prism) + for i in range(easting.size): + msg = "Found observation point" + with pytest.warns(UserWarning, match=msg): + prism_gravity( + (easting[i], northing[i], upward[i]), + sample_prism, + sample_density, + field=field, + ) + + @pytest.mark.use_numba + @pytest.mark.parametrize("field", ("g_ee", "g_zz", "g_ez")) + def test_on_northing_edges(self, sample_prism, sample_density, field): + """ + Test tensor components that have singular points on edges parallel to + easting direction + """ + easting, northing, upward = self.get_northing_edges_center(sample_prism) + for i in range(easting.size): + msg = "Found observation point" + with pytest.warns(UserWarning, match=msg): + prism_gravity( + (easting[i], northing[i], upward[i]), + sample_prism, + sample_density, + field=field, + ) + + @pytest.mark.use_numba + @pytest.mark.parametrize("field", ("g_ee", "g_nn", "g_en")) + def test_on_upward_edges(self, sample_prism, sample_density, field): + """ + Test tensor components that have singular points on edges parallel to + easting direction + """ + easting, northing, upward = self.get_upward_edges_center(sample_prism) + for i in range(easting.size): + msg = "Found observation point" + with pytest.warns(UserWarning, match=msg): + prism_gravity( + (easting[i], northing[i], upward[i]), + sample_prism, + sample_density, + field=field, + ) diff --git a/setup.cfg b/setup.cfg index a9d43db09..5a5101047 100644 --- a/setup.cfg +++ b/setup.cfg @@ -48,6 +48,7 @@ install_requires = xarray>=0.16 verde>=1.7 xrft>=1.0 + choclo>=0.1 [options.extras_require] visualizations =