From 5a34467c528ed66176bd77e7f1cdcbd8efdce3ac Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 17 Apr 2023 14:07:27 -0700 Subject: [PATCH 01/21] Use Choclo kernels for prism gravity forwards Add Choclo as a required dependency of Harmonica and replace the current kernels for prisms gravity with the ones in Choclo. --- environment.yml | 1 + harmonica/_forward/prism.py | 105 +++++----------------------------- harmonica/tests/test_prism.py | 42 +------------- setup.cfg | 1 + 4 files changed, 16 insertions(+), 133 deletions(-) diff --git a/environment.yml b/environment.yml index af140524f..83106c9ed 100644 --- a/environment.yml +++ b/environment.yml @@ -18,6 +18,7 @@ dependencies: - verde>=1.7.0 - xarray>=0.16 - xrft>=1.0 + - choclo>=0.0.1 # Optional requirements - pyvista>=0.27 - vtk>=9 diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index f2c69e5b2..c5c759e8a 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -8,8 +8,12 @@ Forward modelling for prisms """ import numpy as np +from choclo.prism import gravity_pot, gravity_u from numba import jit, prange +# Define dictionary with available gravity fields for prisms +FIELDS = {"potential": gravity_pot, "g_z": gravity_u} + # Attempt to import numba_progress try: from numba_progress import ProgressBar @@ -123,8 +127,7 @@ def prism_gravity( (-0.05379, 0.02908, 0.11235) """ - 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]) @@ -154,15 +157,15 @@ def prism_gravity( progress_proxy = None # Compute gravitational field dispatcher(parallel)( - coordinates, prisms, density, kernels[field], result, progress_proxy + coordinates, prisms, density, FIELDS[field], result, progress_proxy ) - result *= GRAVITATIONAL_CONST # Close previously created progress bars if progressbar: progress_proxy.close() # Convert to more convenient units if field == "g_z": result *= 1e5 # SI to mGal + result *= -1 # invert sign (choclo computes upward component) return result.reshape(cast.shape) @@ -272,103 +275,21 @@ def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy= 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] += kernel( + easting[l], northing[l], upward[l], prisms[m, :], 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..170935292 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -20,13 +20,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 @@ -355,40 +349,6 @@ def test_g_z_symmetry_inside(): 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)) - - @pytest.mark.use_numba def test_prism_against_infinite_slab(): """ diff --git a/setup.cfg b/setup.cfg index a9d43db09..8b9faca67 100644 --- a/setup.cfg +++ b/setup.cfg @@ -48,6 +48,7 @@ install_requires = xarray>=0.16 verde>=1.7 xrft>=1.0 + choclo>=0.0.1 [options.extras_require] visualizations = From 0e5faecb15e230bccc63d7a34a22e90dcdba72d1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 17 Apr 2023 14:48:49 -0700 Subject: [PATCH 02/21] Remove unused GRAVITATIONAL_CONST global --- harmonica/_forward/prism.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index c5c759e8a..4b0f45fe2 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -20,8 +20,6 @@ except ImportError: ProgressBar = None -from ..constants import GRAVITATIONAL_CONST - def prism_gravity( coordinates, From bb51d128413a42db87a7163dfc7bb4d3a466b171 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 17 Apr 2023 15:14:46 -0700 Subject: [PATCH 03/21] Add easting and northing components of acceleration Allow the `prism_gravity` function to compute the easting and northing components of the gravitational acceleration of prisms, using Choclo kernels for it. Add tests that compare with bare Choclo results. Remove the dispatcher and replace it with and `if`/`else` statement for selecting the parallelized or serialized jitted function. --- harmonica/_forward/prism.py | 39 +++++++++---------- harmonica/tests/test_prism.py | 70 +++++++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+), 19 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 4b0f45fe2..39eed3230 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -8,11 +8,16 @@ Forward modelling for prisms """ import numpy as np -from choclo.prism import gravity_pot, gravity_u +from choclo.prism import gravity_pot, gravity_e, gravity_n, gravity_u from numba import jit, prange # Define dictionary with available gravity fields for prisms -FIELDS = {"potential": gravity_pot, "g_z": gravity_u} +FIELDS = { + "potential": gravity_pot, + "g_e": gravity_e, + "g_n": gravity_n, + "g_z": gravity_u, +} # Attempt to import numba_progress try: @@ -72,6 +77,8 @@ def prism_gravity( The available fields are: - Gravitational potential: ``potential`` + - Eastward acceleration: ``g_e`` + - Northward acceleration: ``g_n`` - Downward acceleration: ``g_z`` parallel : bool (optional) @@ -153,31 +160,25 @@ def prism_gravity( progress_proxy = ProgressBar(total=coordinates[0].size) else: progress_proxy = None + # Choose parallelized or serialized forward function + if parallel: + forward_func = jit_prism_gravity_parallel + else: + forward_func = jit_prism_gravity_serial # Compute gravitational field - dispatcher(parallel)( - coordinates, prisms, density, FIELDS[field], result, progress_proxy - ) + forward_func(coordinates, prisms, density, FIELDS[field], result, progress_proxy) # Close previously created progress bars - if progressbar: + if progress_proxy: progress_proxy.close() - # Convert to more convenient units + # Invert sign of gravity_u (upward component) if field == "g_z": + result *= -1 + # Convert to more convenient units + if field in ("g_e", "g_n", "g_z"): result *= 1e5 # SI to mGal - result *= -1 # invert sign (choclo computes upward component) return result.reshape(cast.shape) -def dispatcher(parallel): - """ - Return the parallelized or serialized forward modelling function - """ - dispatchers = { - True: jit_prism_gravity_parallel, - False: jit_prism_gravity_serial, - } - return dispatchers[parallel] - - def _check_prisms(prisms): """ Check if prisms boundaries are well defined diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 170935292..5fc12360d 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -13,6 +13,7 @@ import numpy.testing as npt import pytest import verde as vd +from choclo.prism import gravity_pot, gravity_e, gravity_n, gravity_u try: from numba_progress import ProgressBar @@ -176,6 +177,75 @@ def test_potential_field_symmetry(): npt.assert_allclose(result[0], result) +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), + ], + ) + 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, :], densities[j] + ) + if field in ("g_e", "g_n", "g_z"): + expected_result *= 1e5 + if field == "g_z": + expected_result *= -1 + # 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(): """ From 930efe0c5b8a258d90a0ee14f270f4003b53a711 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 17 Apr 2023 15:54:38 -0700 Subject: [PATCH 04/21] Minor style improvement in test --- harmonica/tests/test_prism.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 5fc12360d..7ba62fda1 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -85,7 +85,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])) From 0e960b463f7cfb3ee68afa05782dd40d333fc39f Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 17 Apr 2023 15:57:50 -0700 Subject: [PATCH 05/21] Allow to compute gravity tensor components Add tests against dumb Choclo calls and a Laplace equation test for the diagonal tensor components. Improve the warning on the docstring explaining the direction of `z`. --- harmonica/_forward/prism.py | 36 ++++++++++++++++++++----- harmonica/tests/test_prism.py | 49 ++++++++++++++++++++++++++++++++--- 2 files changed, 75 insertions(+), 10 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 39eed3230..8f445d8fc 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -8,7 +8,18 @@ Forward modelling for prisms """ import numpy as np -from choclo.prism import gravity_pot, gravity_e, gravity_n, gravity_u +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 numba import jit, prange # Define dictionary with available gravity fields for prisms @@ -17,6 +28,12 @@ "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 @@ -51,11 +68,13 @@ def prism_gravity( solution has on some points (see [Nagy2000]_). .. warning:: - The **z direction points upwards**, i.e. positive and negative values - of ``upward`` represent points above and below the surface, + The **vertical 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. + 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 ---------- @@ -80,6 +99,8 @@ def prism_gravity( - 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 @@ -170,12 +191,15 @@ def prism_gravity( # Close previously created progress bars if progress_proxy: progress_proxy.close() - # Invert sign of gravity_u (upward component) - if field == "g_z": + # 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 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) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 7ba62fda1..fb1e8f23b 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -13,7 +13,18 @@ import numpy.testing as npt import pytest import verde as vd -from choclo.prism import gravity_pot, gravity_e, gravity_n, gravity_u +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 @@ -216,6 +227,12 @@ def sample_coordinates(self): ("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( @@ -237,15 +254,39 @@ def test_against_choclo( expected_result[i] += choclo_func( easting[i], northing[i], upward[i], prisms[j, :], 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 - if field == "g_z": - expected_result *= -1 + 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) +@run_only_with_numba +def test_laplace(): + """ + Test if the diagonal components satisfy Laplace equation + """ + 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") + } + npt.assert_allclose( + diagonal_components["g_ee"] + diagonal_components["g_nn"], + -diagonal_components["g_zz"], + ) + + @pytest.mark.use_numba def test_g_z_symmetry_outside(): """ From c7ec09639bcb2fe3aec304a8720a77fa88cfbb39 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 17 Apr 2023 15:59:35 -0700 Subject: [PATCH 06/21] Ditch symmetry tests The tests in Choclo already cover those and they are more extensive. --- harmonica/tests/test_prism.py | 216 ---------------------------------- 1 file changed, 216 deletions(-) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index fb1e8f23b..20ae43f12 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -145,49 +145,6 @@ 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], - ) - result = prism_gravity(coordinates, prism, density, field="potential") - npt.assert_allclose(result[0], result) - - class TestAgainstChoclo: """ Test forward modelling functions against dumb Choclo runs @@ -287,179 +244,6 @@ def test_laplace(): ) -@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(): - """ - 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. - """ - 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]), - } - # 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_prism_against_infinite_slab(): """ From 2498439157559e4b0a6ad80c074827eb4570f7c4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 18 Apr 2023 12:13:30 -0700 Subject: [PATCH 07/21] Warn for prism gravity tensor on singular points Check if the gravity tensor components of a prism are being computed at any singular point. Choclo will return a ``np.nan`` on those points, but it's worth warning the user about it before they get the results. Add tests that check if the warnings are being raised on those cases. --- harmonica/_forward/prism.py | 138 ++++++++++++++++++++++++++++++++++ harmonica/tests/test_prism.py | 125 ++++++++++++++++++++++++++++++ 2 files changed, 263 insertions(+) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 8f445d8fc..dd87265f5 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -7,6 +7,8 @@ """ Forward modelling for prisms """ +import warnings + import numpy as np from choclo.prism import ( gravity_e, @@ -20,6 +22,11 @@ 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 @@ -170,6 +177,7 @@ 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 @@ -203,6 +211,136 @@ def prism_gravity( return result.reshape(cast.shape) +def _check_singular_points(coordinates, prisms, field): + """ + 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). + """ + any_singular_point = False + easting, northing, upward = coordinates + if field == "g_ee": + any_singular_point = _any_singular_point_diagonal( + easting, + northing, + upward, + prisms, + is_point_on_northing_edge, + is_point_on_upward_edge, + ) + elif field == "g_nn": + any_singular_point = _any_singular_point_diagonal( + easting, + northing, + upward, + prisms, + is_point_on_easting_edge, + is_point_on_upward_edge, + ) + elif field == "g_zz": + any_singular_point = _any_singular_point_diagonal( + easting, + northing, + upward, + prisms, + is_point_on_easting_edge, + is_point_on_northing_edge, + ) + elif field == "g_en": + any_singular_point = _any_singular_point_non_diagonal( + easting, northing, upward, prisms, is_point_on_upward_edge + ) + elif field == "g_ez": + any_singular_point = _any_singular_point_non_diagonal( + easting, northing, upward, prisms, is_point_on_northing_edge + ) + elif field == "g_nz": + any_singular_point = _any_singular_point_non_diagonal( + easting, northing, upward, prisms, is_point_on_easting_edge + ) + if any_singular_point: + warnings.warn( + "Found observation point on singular point of a prism.", UserWarning + ) + + +@jit(nopython=True) +def _any_singular_point_diagonal( + easting, northing, upward, prisms, check_func1, check_func2 +): + """ + Check singular points for diagonal tensor components + """ + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if check_func1( + easting[l], northing[l], upward[l], prisms[m, :] + ) or check_func2(easting[l], northing[l], upward[l], prisms[m, :]): + return True + return False + + +@jit(nopython=True) +def _any_singular_point_non_diagonal(easting, northing, upward, prisms, check_func): + """ + Check singular points for non-diagonal tensor components + """ + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in range(n_coords): + for m in range(n_prisms): + if check_func(easting[l], northing[l], upward[l], prisms[m, :]): + return True + return False + + +@jit(nopython=True, parallel=True) +def _check_singular_points_diagonal(easting, northing, upward, prisms, check_functions): + """ + Check singular points for diagonal tensor components + """ + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in prange(n_coords): + for m in prange(n_prisms): + if check_func1( + easting[l], northing[l], upward[l], prisms[:, m] + ) or check_func2(easting[l], northing[l], upward[l], prisms[:, m]): + raise ValueError( + "Found observation point " + f"'({easting[l]}, {northing[l]}, {upward[l]})' located at a " + f"singular point for prism '{prisms[:, m]}'." + ) + + +@jit(nopython=True, parallel=True) +def _check_singular_points_non_diagonal(easting, northing, upward, prisms, check_func): + """ + Check singular points for non-diagonal tensor components + """ + n_coords = easting.size + n_prisms = prisms.shape[0] + for l in prange(n_coords): + for m in prange(n_prisms): + if check_func(easting[l], northing[l], upward[l], prisms[:, m]): + raise ValueError( + "Found observation point " + f"'({easting[l]}, {northing[l]}, {upward[l]})' located at a " + f"singular point for prism '{prisms[:, m]}'." + ) + + def _check_prisms(prisms): """ Check if prisms boundaries are well defined diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 20ae43f12..0b0007077 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -341,3 +341,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, + ) From d0f0d2e05ff029887dc7749f0679d84035833d72 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 24 Apr 2023 12:03:56 -0700 Subject: [PATCH 08/21] Fix wrong word in test comment --- harmonica/tests/test_prism.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 20ae43f12..aa58a66ff 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -47,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], From 32ecfc97886d24a8825996ec937122d9d2b30f1c Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 5 May 2023 10:47:47 -0700 Subject: [PATCH 09/21] Update to new Choclo functions that take only scalars --- harmonica/_forward/prism.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 8f445d8fc..d802cdb10 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -306,7 +306,16 @@ def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy= for l in prange(easting.size): for m in range(prisms.shape[0]): out[l] += kernel( - easting[l], northing[l], upward[l], prisms[m, :], density[m] + 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: From d397f3728534a356fa842d80476efd883391257f Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 5 May 2023 10:48:20 -0700 Subject: [PATCH 10/21] Rename "kernel" argument to "forward_func" --- harmonica/_forward/prism.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index d802cdb10..023d78ebe 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -274,7 +274,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 @@ -292,8 +294,9 @@ 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. out : 1d-array Array where the resulting field values will be stored. Must have the same size as the arrays contained on ``coordinates``. @@ -305,7 +308,7 @@ def jit_prism_gravity(coordinates, prisms, density, kernel, out, progress_proxy= # Iterate over computation points and prisms for l in prange(easting.size): for m in range(prisms.shape[0]): - out[l] += kernel( + out[l] += forward_func( easting[l], northing[l], upward[l], From f205e3d670886df6939bb74e48c3bac20e463255 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 5 May 2023 10:59:36 -0700 Subject: [PATCH 11/21] Update tests to use the new Choclo functions --- environment.yml | 2 +- harmonica/tests/test_prism.py | 11 ++++++++++- setup.cfg | 2 +- 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/environment.yml b/environment.yml index 83106c9ed..9ab8d5697 100644 --- a/environment.yml +++ b/environment.yml @@ -18,7 +18,7 @@ dependencies: - verde>=1.7.0 - xarray>=0.16 - xrft>=1.0 - - choclo>=0.0.1 + - choclo>=0.1 # Optional requirements - pyvista>=0.27 - vtk>=9 diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index aa58a66ff..65a6cf8e8 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -209,7 +209,16 @@ def test_against_choclo( 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, :], densities[j] + 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 diff --git a/setup.cfg b/setup.cfg index 8b9faca67..5a5101047 100644 --- a/setup.cfg +++ b/setup.cfg @@ -48,7 +48,7 @@ install_requires = xarray>=0.16 verde>=1.7 xrft>=1.0 - choclo>=0.0.1 + choclo>=0.1 [options.extras_require] visualizations = From 83fa73184e5a894b52412c040ae0dcd2b52fda1a Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 12 May 2023 14:36:24 -0700 Subject: [PATCH 12/21] Update expected values in example Update the expected values after we updated the gravitational constant in Choclo. --- harmonica/_forward/prism.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 023d78ebe..56a0494ca 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -141,7 +141,7 @@ def prism_gravity( >>> # the prism generates on the computation points >>> gz = prism_gravity(coordinates, prism, density, field="g_z") >>> print("({:.5f}, {:.5f}, {:.5f})".format(*gz)) - (0.06551, 0.06628, 0.06173) + (0.06552, 0.06629, 0.06174) Define two prisms with positive and negative density contrasts @@ -150,7 +150,7 @@ def prism_gravity( >>> # Compute the g_z that the prisms generate on the computation points >>> gz = prism_gravity(coordinates, prisms, densities, field="g_z") >>> print("({:.5f}, {:.5f}, {:.5f})".format(*gz)) - (-0.05379, 0.02908, 0.11235) + (-0.05380, 0.02908, 0.11237) """ if field not in FIELDS: From fa732f66f425f716c28a0e9d1738e060f69a0c04 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 12 May 2023 16:54:29 -0700 Subject: [PATCH 13/21] Fix test function against infinite slab Modify the test function that compares prism forward model with the analytic solution for an infinite slab: locate the observation point on top of the slab. --- harmonica/tests/test_prism.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 65a6cf8e8..77d7f6238 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -258,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 + # Define an observation point at 1.5m above zero + height = 1.5 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 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] @@ -274,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]) From 16e317229b18686f8a336d864567bbc35bd2cdcf Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 12 May 2023 17:01:34 -0700 Subject: [PATCH 14/21] =?UTF-8?q?Actually,=20use=20the=20height=20variable?= =?UTF-8?q?=20=F0=9F=A4=A6?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- harmonica/tests/test_prism.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/test_prism.py b/harmonica/tests/test_prism.py index 77d7f6238..9751177c6 100644 --- a/harmonica/tests/test_prism.py +++ b/harmonica/tests/test_prism.py @@ -260,7 +260,7 @@ def test_prism_against_infinite_slab(): """ # Define an observation point at 1.5m above zero height = 1.5 - coordinates = (0, 0, 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. From 7305fc305b31be2d2523de92d30948c4b5ce9aa1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 17 May 2023 15:09:10 -0700 Subject: [PATCH 15/21] Rename forward_func variable to avoid confusion Rename forward_func variable to avoid confusion with the forward_func argument of the jit functions. --- harmonica/_forward/prism.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 56a0494ca..062f1aece 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -183,11 +183,13 @@ def prism_gravity( progress_proxy = None # Choose parallelized or serialized forward function if parallel: - forward_func = jit_prism_gravity_parallel + gravity_prism_func = jit_prism_gravity_parallel else: - forward_func = jit_prism_gravity_serial + gravity_prism_func = jit_prism_gravity_serial # Compute gravitational field - forward_func(coordinates, prisms, density, FIELDS[field], result, progress_proxy) + gravity_prism_func( + coordinates, prisms, density, FIELDS[field], result, progress_proxy + ) # Close previously created progress bars if progress_proxy: progress_proxy.close() @@ -296,7 +298,8 @@ def jit_prism_gravity( same size as the number of prisms. forward_func : func Forward modelling function that will be used to compute the desired - field. + 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``. From ba115ce8658362a1651bf1e1f8d9b1c49937efa2 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 17 May 2023 15:37:45 -0700 Subject: [PATCH 16/21] Improve docstrings of prism_gravity Move some details to a new Notes section, improve some definition of the parameters, improve the warning regarding the upward direction, add a References section. --- harmonica/_forward/prism.py | 74 +++++++++++++++++++++++++------------ 1 file changed, 50 insertions(+), 24 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 062f1aece..47aca3c95 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -56,39 +56,32 @@ 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 **vertical 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 same applies to the tensor - components, i.e. the ``g_ez`` is the non-diagonal easting-downward - tensor component. + 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 ---------- - coordinates : list or 1d-array - List or array containing ``easting``, ``northing`` and ``upward`` of - the computation points defined on a Cartesian coordinate system. - All coordinates should be in meters. + coordinates : list of arrays + List of arrays containing the ``easting``, ``northing`` and ``upward`` + coordinates of the computation points defined on a Cartesian coordinate + 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. + 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 @@ -125,6 +118,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 functions returns the limit of these components while approaching from + the outside. + + References + ---------- + * [Nagy2000]_ + * [Nagy2002]_ + * [Fukushima2020]_ + Examples -------- @@ -134,9 +157,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") From 8694500240a8e4dad7eb87913d853b1bf8645d1f Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 17 May 2023 16:29:40 -0700 Subject: [PATCH 17/21] Fix check for singular points Redesign the checker for singular points to make it simpler and fix the issue with Numba that made some tests to fail. --- harmonica/_forward/prism.py | 166 ++++++++++++++++++------------------ 1 file changed, 85 insertions(+), 81 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index c43352c03..e8720534a 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -210,7 +210,8 @@ def prism_gravity( 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: @@ -254,119 +255,122 @@ def _check_singular_points(coordinates, prisms, field): two directions of the component (e.g. ``g_en`` is not defined on edges parallel to the upward direction). """ - any_singular_point = False - easting, northing, upward = coordinates - if field == "g_ee": - any_singular_point = _any_singular_point_diagonal( - easting, - northing, - upward, - prisms, - is_point_on_northing_edge, - is_point_on_upward_edge, - ) - elif field == "g_nn": - any_singular_point = _any_singular_point_diagonal( - easting, - northing, - upward, - prisms, - is_point_on_easting_edge, - is_point_on_upward_edge, - ) - elif field == "g_zz": - any_singular_point = _any_singular_point_diagonal( - easting, - northing, - upward, - prisms, - is_point_on_easting_edge, - is_point_on_northing_edge, - ) - elif field == "g_en": - any_singular_point = _any_singular_point_non_diagonal( - easting, northing, upward, prisms, is_point_on_upward_edge - ) - elif field == "g_ez": - any_singular_point = _any_singular_point_non_diagonal( - easting, northing, upward, prisms, is_point_on_northing_edge - ) - elif field == "g_nz": - any_singular_point = _any_singular_point_non_diagonal( - easting, northing, upward, prisms, is_point_on_easting_edge - ) - if any_singular_point: + 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, + } + 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_diagonal( - easting, northing, upward, prisms, check_func1, check_func2 -): +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, :] + ) or is_point_on_upward_edge( + easting[l], northing[l], upward[l], *prisms[m, :] + ): + return True + return False + + +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, :] + ) or is_point_on_upward_edge( + easting[l], northing[l], upward[l], *prisms[m, :] + ): + return True + return False + + +def _any_singular_point_g_zz(coordinates, prisms): """ - Check singular points for diagonal tensor components + 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 check_func1( - easting[l], northing[l], upward[l], prisms[m, :] - ) or check_func2(easting[l], northing[l], upward[l], prisms[m, :]): + if is_point_on_easting_edge( + easting[l], northing[l], upward[l], *prisms[m, :] + ) or is_point_on_northing_edge( + easting[l], northing[l], upward[l], *prisms[m, :] + ): return True return False -@jit(nopython=True) -def _any_singular_point_non_diagonal(easting, northing, upward, prisms, check_func): +def _any_singular_point_g_en(coordinates, prisms): """ - Check singular points for non-diagonal tensor components + 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 check_func(easting[l], northing[l], upward[l], prisms[m, :]): + if is_point_on_upward_edge( + easting[l], northing[l], upward[l], *prisms[m, :] + ): return True return False -@jit(nopython=True, parallel=True) -def _check_singular_points_diagonal(easting, northing, upward, prisms, check_functions): +def _any_singular_point_g_ez(coordinates, prisms): """ - Check singular points for diagonal tensor components + 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 prange(n_coords): - for m in prange(n_prisms): - if check_func1( - easting[l], northing[l], upward[l], prisms[:, m] - ) or check_func2(easting[l], northing[l], upward[l], prisms[:, m]): - raise ValueError( - "Found observation point " - f"'({easting[l]}, {northing[l]}, {upward[l]})' located at a " - f"singular point for prism '{prisms[:, m]}'." - ) - - -@jit(nopython=True, parallel=True) -def _check_singular_points_non_diagonal(easting, northing, upward, prisms, check_func): + 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, :] + ): + return True + return False + + +def _any_singular_point_g_nz(coordinates, prisms): """ - Check singular points for non-diagonal tensor components + 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 prange(n_coords): - for m in prange(n_prisms): - if check_func(easting[l], northing[l], upward[l], prisms[:, m]): - raise ValueError( - "Found observation point " - f"'({easting[l]}, {northing[l]}, {upward[l]})' located at a " - f"singular point for prism '{prisms[:, m]}'." - ) + 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, :] + ): + return True + return False def _check_prisms(prisms): From 70250f7fced55e1547ef3f7db565d33024549cd3 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 17 May 2023 16:30:28 -0700 Subject: [PATCH 18/21] Fix line too long --- harmonica/_forward/prism.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index e8720534a..dc0874d89 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -84,8 +84,8 @@ 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 row). From fa139292688853c659ef8020f02d9b491a747c5d Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 23 May 2023 11:26:05 -0700 Subject: [PATCH 19/21] Fix typo Co-authored-by: Lu Li <54405391+LL-Geo@users.noreply.github.com> --- harmonica/_forward/prism.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index dc0874d89..83859791e 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -146,7 +146,7 @@ def prism_gravity( 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 functions returns the limit of these components while approaching from + This function returns the limit of these components while approaching from the outside. References From 59a2dff7184aa11635b5b27adc7bd48f80af58d9 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 23 May 2023 12:18:37 -0700 Subject: [PATCH 20/21] Run checks for singular points with Numba Run checks for singular points using Numba jitted functions to speed up the for loops. --- harmonica/_forward/prism.py | 96 +++++++++++++++++++++++++++++++++---- 1 file changed, 87 insertions(+), 9 deletions(-) diff --git a/harmonica/_forward/prism.py b/harmonica/_forward/prism.py index 83859791e..3ac4013ec 100644 --- a/harmonica/_forward/prism.py +++ b/harmonica/_forward/prism.py @@ -271,6 +271,7 @@ def _check_singular_points(coordinates, prisms, field): ) +@jit(nopython=True) def _any_singular_point_g_ee(coordinates, prisms): """ Check observation points as singular points of g_ee @@ -281,14 +282,31 @@ def _any_singular_point_g_ee(coordinates, prisms): 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, :] + 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, :] + 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 @@ -299,14 +317,31 @@ def _any_singular_point_g_nn(coordinates, prisms): 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, :] + 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, :] + 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 @@ -317,14 +352,31 @@ def _any_singular_point_g_zz(coordinates, prisms): 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, :] + 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, :] + 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 @@ -335,12 +387,21 @@ def _any_singular_point_g_en(coordinates, prisms): 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, :] + 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 @@ -351,12 +412,21 @@ def _any_singular_point_g_ez(coordinates, prisms): 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, :] + 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 @@ -367,7 +437,15 @@ def _any_singular_point_g_nz(coordinates, prisms): 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, :] + 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 From 13da7a5576755bf032879a1fee051beef32253e1 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Wed, 24 May 2023 11:39:41 -0700 Subject: [PATCH 21/21] Add user guide section for other fields for prisms Add a section to the prism forward modelling page in the user guide showing how we can compute all gravitational acceleration and tensor components of a single prism. --- doc/user_guide/forward_modelling/prism.rst | 91 ++++++++++++++++++++++ 1 file changed, 91 insertions(+) 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 ^^^^^^^^^^^^^^^^^^^^^^^