From b029fefc5370cb929facbc12363b3b846dca2af6 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 22 Oct 2021 10:42:38 -0300 Subject: [PATCH 01/17] Add scipy to requirements --- doc/install.rst | 1 + environment.yml | 1 + requirements.txt | 1 + 3 files changed, 3 insertions(+) diff --git a/doc/install.rst b/doc/install.rst index f4a3443d1..e814a77a7 100644 --- a/doc/install.rst +++ b/doc/install.rst @@ -22,6 +22,7 @@ Dependencies * `numpy `__ * `pandas `__ * `numba `__ +* `scipy `__ * `xarray `__ * `scikit-learn `__ * `pooch `__ diff --git a/environment.yml b/environment.yml index 805f8db3f..605121bc4 100644 --- a/environment.yml +++ b/environment.yml @@ -9,6 +9,7 @@ dependencies: - numpy - pandas - numba + - scipy - scikit-learn - pooch>=0.7.0 - verde>=1.5.0 diff --git a/requirements.txt b/requirements.txt index 5efe7d232..1abc0dbfa 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ numpy pandas +scipy scikit-learn numba pooch>=0.7.0 From dd1d9705228079fcfae9bceacb6a34fea1f5e911 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 22 Oct 2021 10:42:56 -0300 Subject: [PATCH 02/17] Draft variable density tesseroids implementation Refactor tesseroid code by moving some functions to a new _tesseroid_utils.py file. Create a new _tesseroid_variable_density.py file with density-based discretization functions and with a gauss_legendre_quadrature function that works with variable density. Add new forward modelling functions for variable density tesseroids only. Add new example that makes use of variable density tesseroids. --- .../forward/tesseroid_variable_density.py | 62 ++ harmonica/forward/_tesseroid_utils.py | 499 ++++++++++++++ .../forward/_tesseroid_variable_density.py | 221 +++++++ harmonica/forward/tesseroid.py | 606 ++++-------------- 4 files changed, 902 insertions(+), 486 deletions(-) create mode 100644 examples/forward/tesseroid_variable_density.py create mode 100644 harmonica/forward/_tesseroid_utils.py create mode 100644 harmonica/forward/_tesseroid_variable_density.py diff --git a/examples/forward/tesseroid_variable_density.py b/examples/forward/tesseroid_variable_density.py new file mode 100644 index 000000000..46ae37763 --- /dev/null +++ b/examples/forward/tesseroid_variable_density.py @@ -0,0 +1,62 @@ +# Copyright (c) 2018 The Harmonica Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +# +# This code is part of the Fatiando a Terra project (https://www.fatiando.org) +# +""" +Tesseroids with variable density +================================ +""" +import matplotlib.pyplot as plt +import cartopy.crs as ccrs +import verde as vd +import boule as bl +import harmonica as hm + + +# Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to +# reference the tesseroid +ellipsoid = bl.WGS84 +mean_radius = ellipsoid.mean_radius + +# Define tesseroid with top surface at the mean Earth radius, a thickness of +# 10km and a linear density function +tesseroids = ( + [-70, -60, -40, -30, mean_radius - 3e3, mean_radius], + [-70, -60, -30, -20, mean_radius - 5e3, mean_radius], + [-60, -50, -40, -30, mean_radius - 7e3, mean_radius], + [-60, -50, -30, -20, mean_radius - 10e3, mean_radius], +) + + +def density(radius): + """Linear density function""" + top = mean_radius + bottom = mean_radius - 10e3 + density_top = 2670 + density_bottom = 3000 + slope = (density_top - density_bottom) / (top - bottom) + return slope * (radius - bottom) + density_bottom + + +# Define computation points on a regular grid at 100km above the mean Earth +# radius +coordinates = vd.grid_coordinates( + region=[-80, -40, -50, -10], + shape=(80, 80), + extra_coords=100e3 + ellipsoid.mean_radius, +) + +# Compute the radial component of the acceleration +gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z") +print(gravity) + +# Plot the gravitational field +fig = plt.figure(figsize=(8, 9)) +ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60)) +img = ax.pcolormesh(*coordinates[:2], gravity, transform=ccrs.PlateCarree()) +plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal") +ax.coastlines() +ax.set_title("Downward component of gravitational acceleration") +plt.show() diff --git a/harmonica/forward/_tesseroid_utils.py b/harmonica/forward/_tesseroid_utils.py new file mode 100644 index 000000000..5f952569e --- /dev/null +++ b/harmonica/forward/_tesseroid_utils.py @@ -0,0 +1,499 @@ +# Copyright (c) 2018 The Harmonica Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +# +# This code is part of the Fatiando a Terra project (https://www.fatiando.org) +# +""" +Utils functions for tesseroid forward modelling +""" +from numba import jit +import numpy as np +from numpy.polynomial.legendre import leggauss + +from .utils import distance_spherical + + +@jit(nopython=True) +def gauss_legendre_quadrature( + longitude, + cosphi, + sinphi, + radius, + tesseroid, + density, + glq_nodes, + glq_weights, + kernel, +): # pylint: disable=too-many-locals + r""" + Compute the effect of a tesseroid on a single observation point through GLQ + + The tesseroid is converted into a set of point masses located on the + scaled nodes of the Gauss-Legendre Quadrature. The number of point masses + created from each tesseroid is equal to the product of the GLQ degrees for + each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass + of each point mass is defined as the product of the tesseroid density + (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`, + :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and + the :math:`\kappa` factor evaluated on the coordinates of the point mass. + + Parameters + ---------- + longitude : float + Longitudinal coordinate of the observation points in radians. + cosphi : float + Cosine of the latitudinal coordinate of the observation point in + radians. + sinphi : float + Sine of the latitudinal coordinate of the observation point in + radians. + radius : float + Radial coordinate of the observation point in meters. + tesseroids : 1d-array + Array containing the boundaries of the tesseroid: + ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + Horizontal boundaries should be in degrees and radial boundaries in + meters. + density : float + Density of the tesseroid in SI units. + glq_nodes : list + Unscaled location of GLQ nodes for each direction. + glq_weights : list + GLQ weigths for each node for each direction. + kernel : func + Kernel function for the gravitational field of point masses. + + """ + # Get tesseroid boundaries + w, e, s, n, bottom, top = tesseroid[:] + # Calculate the A factor for the tesseroid + a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom) + # Unpack nodes and weights + lon_nodes, lat_nodes, rad_nodes = glq_nodes[:] + lon_weights, lat_weights, rad_weights = glq_weights[:] + # Compute effect of the tesseroid on the observation point + # by iterating over the location of the point masses + # (move the iteration along the longitudinal nodes to the bottom for + # optimization: reduce the number of times that the trigonometric functions + # are evaluated) + result = 0.0 + for j, lat_node in enumerate(lat_nodes): + # Get the latitude of the point mass + latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s)) + cosphi_p = np.cos(latitude_p) + sinphi_p = np.sin(latitude_p) + for k, rad_node in enumerate(rad_nodes): + # Get the radius of the point mass + radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom) + # Get kappa constant for the point mass + kappa = radius_p ** 2 * cosphi_p + for i, lon_node in enumerate(lon_nodes): + # Get the longitude of the point mass + longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w)) + # Compute the mass of the point mass + mass = ( + density + * a_factor + * kappa + * lon_weights[i] + * lat_weights[j] + * rad_weights[k] + ) + # Add effect of the current point mass to the result + result += mass * kernel( + longitude, + cosphi, + sinphi, + radius, + longitude_p, + cosphi_p, + sinphi_p, + radius_p, + ) + return result + + +def glq_nodes_weights(glq_degrees): + """ + Calculate GLQ unscaled nodes and weights + + Parameters + ---------- + glq_degrees : list + List of GLQ degrees for each direction: ``longitude``, ``latitude``, + ``radius``. + + Returns + ------- + glq_nodes : list + Unscaled GLQ nodes for each direction: ``longitude``, ``latitude``, + ``radius``. + glq_weights : list + GLQ weights for each node on each direction: ``longitude``, + ``latitude``, ``radius``. + """ + # Unpack GLQ degrees + lon_degree, lat_degree, rad_degree = glq_degrees[:] + # Get nodes coordinates and weights + lon_node, lon_weights = leggauss(lon_degree) + lat_node, lat_weights = leggauss(lat_degree) + rad_node, rad_weights = leggauss(rad_degree) + # Reorder nodes and weights + glq_nodes = (lon_node, lat_node, rad_node) + glq_weights = (lon_weights, lat_weights, rad_weights) + return glq_nodes, glq_weights + + +@jit(nopython=True) +def _adaptive_discretization( + coordinates, + tesseroid, + distance_size_ratio, + stack, + small_tesseroids, + radial_discretization=False, +): + """ + Perform the adaptive discretization algorithm on a tesseroid + + It apply the three or two dimensional adaptive discretization algorithm on + a tesseroid after a single computation point. + + Parameters + ---------- + coordinates : array + Array containing ``longitude``, ``latitude`` and ``radius`` of a single + computation point. + tesseroid : array + Array containing the boundaries of the tesseroid. + distance_size_ratio : float + Value for the distance-size ratio. A greater value will perform more + discretizations. + stack : 2d-array + Array with shape ``(6, stack_size)`` that will temporarly hold the + small tesseroids that are not yet processed. + If too many discretizations will take place, increase the + ``stack_size``. + small_tesseroids : 2d-array + Array with shape ``(6, MAX_DISCRETIZATIONS)`` that will contain every + small tesseroid produced by the adaptive discretization algorithm. + If too many discretizations will take place, increase the + ``MAX_DISCRETIZATIONS``. + radial_discretization : bool (optional) + If ``True`` the three dimensional adaptive discretization will be + applied. + If ``False`` the two dimensional adaptive discretization will be + applied, i.e. the tesseroid will only be split on the ``longitude`` and + ``latitude`` directions. + Default ``False``. + + Returns + ------- + n_splits : int + Total number of small tesseroids generated by the algorithm. + """ + # Create stack of tesseroids + stack[0] = tesseroid + stack_top = 0 + n_splits = 0 + while stack_top >= 0: + # Pop the first tesseroid from the stack + tesseroid = stack[stack_top] + stack_top -= 1 + # Get its dimensions + l_lon, l_lat, l_rad = _tesseroid_dimensions(tesseroid) + # Get distance between computation point and center of tesseroid + distance = _distance_tesseroid_point(coordinates, tesseroid) + # Check inequality + n_lon, n_lat, n_rad = 1, 1, 1 + if distance / l_lon < distance_size_ratio: + n_lon = 2 + if distance / l_lat < distance_size_ratio: + n_lat = 2 + if distance / l_rad < distance_size_ratio and radial_discretization: + n_rad = 2 + # Apply discretization + if n_lon * n_lat * n_rad > 1: + # Raise error if stack overflow + # Number of tesseroids in stack = stack_top + 1 + if (stack_top + 1) + n_lon * n_lat * n_rad > stack.shape[0]: + raise OverflowError("Stack Overflow. Try to increase the stack size.") + stack_top = _split_tesseroid( + tesseroid, n_lon, n_lat, n_rad, stack, stack_top + ) + else: + # Raise error if small_tesseroids overflow + if n_splits + 1 > small_tesseroids.shape[0]: + raise OverflowError( + "Exceeded maximum discretizations." + + " Please increase the MAX_DISCRETIZATIONS." + ) + small_tesseroids[n_splits] = tesseroid + n_splits += 1 + return n_splits + + +@jit(nopython=True) +def _split_tesseroid( + tesseroid, n_lon, n_lat, n_rad, stack, stack_top +): # pylint: disable=too-many-locals + """ + Split tesseroid along each dimension + """ + w, e, s, n, bottom, top = tesseroid[:] + # Compute differential distance + d_lon = (e - w) / n_lon + d_lat = (n - s) / n_lat + d_rad = (top - bottom) / n_rad + for i in range(n_lon): + for j in range(n_lat): + for k in range(n_rad): + stack_top += 1 + stack[stack_top, 0] = w + d_lon * i + stack[stack_top, 1] = w + d_lon * (i + 1) + stack[stack_top, 2] = s + d_lat * j + stack[stack_top, 3] = s + d_lat * (j + 1) + stack[stack_top, 4] = bottom + d_rad * k + stack[stack_top, 5] = bottom + d_rad * (k + 1) + return stack_top + + +@jit(nopython=True) +def _tesseroid_dimensions(tesseroid): + """ + Calculate the dimensions of the tesseroid. + """ + w, e, s, n, bottom, top = tesseroid[:] + w, e, s, n = np.radians(w), np.radians(e), np.radians(s), np.radians(n) + latitude_center = (n + s) / 2 + l_lat = top * np.arccos(np.sin(n) * np.sin(s) + np.cos(n) * np.cos(s)) + l_lon = top * np.arccos( + np.sin(latitude_center) ** 2 + np.cos(latitude_center) ** 2 * np.cos(e - w) + ) + l_rad = top - bottom + return l_lon, l_lat, l_rad + + +@jit(nopython=True) +def _distance_tesseroid_point( + coordinates, tesseroid +): # pylint: disable=too-many-locals + """ + Distance between a computation point and the center of a tesseroid + """ + # Get center of the tesseroid + w, e, s, n, bottom, top = tesseroid[:] + longitude_p = (w + e) / 2 + latitude_p = (s + n) / 2 + radius_p = (bottom + top) / 2 + # Get distance between computation point and tesseroid center + distance = distance_spherical(coordinates, (longitude_p, latitude_p, radius_p)) + return distance + + +def _check_tesseroids(tesseroids): # pylint: disable=too-many-branches + """ + Check if tesseroids boundaries are well defined + + A valid tesseroid should have: + - latitudinal boundaries within the [-90, 90] degrees interval, + - north boundaries greater or equal than the south boundaries, + - radial boundaries positive or zero, + - top boundaries greater or equal than the bottom boundaries, + - longitudinal boundaries within the [-180, 360] degrees interval, + - longitudinal interval must not be greater than one turn around the + globe. + + Some valid tesseroids have its west boundary greater than the east one, + e.g. ``(350, 10, ...)``. On these cases the ``_longitude_continuity`` + function is applied in order to move the longitudinal coordinates to the + [-180, 180) interval. Any valid tesseroid should have east boundaries + greater than the west boundaries before or after applying longitude + continuity. + + Parameters + ---------- + tesseroids : 2d-array + Array containing the boundaries of the tesseroids in the following + order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + Longitudinal and latitudinal boundaries must be in degrees. + The array must have the following shape: (``n_tesseroids``, 6), where + ``n_tesseroids`` is the total number of tesseroids. + + Returns + ------- + tesseroids : 2d-array + Array containing the boundaries of the tesseroids. + If no longitude continuity needs to be applied, the returned array is + the same one as the orignal. + Otherwise, it's copied and its longitudinal boundaries are modified. + """ + west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6)) + err_msg = "Invalid tesseroid or tesseroids. " + # Check if latitudinal boundaries are inside the [-90, 90] interval + invalid = np.logical_or( + np.logical_or(south < -90, south > 90), np.logical_or(north < -90, north > 90) + ) + if (invalid).any(): + err_msg += ( + "The latitudinal boundaries must be inside the [-90, 90] " + + "degrees interval.\n" + ) + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + # Check if south boundary is not greater than the corresponding north + # boundary + invalid = south > north + if (invalid).any(): + err_msg += "The south boundary can't be greater than the north one.\n" + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + # Check if radial boundaries are positive or zero + invalid = np.logical_or(bottom < 0, top < 0) + if (invalid).any(): + err_msg += "The bottom and top radii should be positive or zero.\n" + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + # Check if top boundary is not greater than the corresponding bottom + # boundary + invalid = bottom > top + if (invalid).any(): + err_msg += "The bottom radius boundary can't be greater than the top one.\n" + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + # Check if longitudinal boundaries are inside the [-180, 360] interval + invalid = np.logical_or( + np.logical_or(west < -180, west > 360), np.logical_or(east < -180, east > 360) + ) + if (invalid).any(): + err_msg += ( + "The longitudinal boundaries must be inside the [-180, 360] " + + "degrees interval.\n" + ) + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + # Apply longitude continuity if w > e + if (west > east).any(): + tesseroids = _longitude_continuity(tesseroids) + west, east, south, north, bottom, top = tuple( + tesseroids[:, i] for i in range(6) + ) + # Check if west boundary is not greater than the corresponding east + # boundary, even after applying the longitude continuity + invalid = west > east + if (invalid).any(): + err_msg += "The west boundary can't be greater than the east one.\n" + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + # Check if the longitudinal interval is not grater than one turn around the + # globe + invalid = east - west > 360 + if (invalid).any(): + err_msg += ( + "The difference between east and west boundaries cannot be greater than " + + "one turn around the globe.\n" + ) + for tess in tesseroids[invalid]: + err_msg += "\tInvalid tesseroid: {}\n".format(tess) + raise ValueError(err_msg) + return tesseroids + + +def _check_points_outside_tesseroids( + coordinates, tesseroids +): # pylint: disable=too-many-locals + """ + Check if computation points are not inside the tesseroids + + Parameters + ---------- + coordinates : 2d-array + Array containing the coordinates of the computation points in the + following order: ``longitude``, ``latitude`` and ``radius``. + Both ``longitude`` and ``latitude`` must be in degrees. + The array must have the following shape: (3, ``n_points``), where + ``n_points`` is the total number of computation points. + tesseroids : 2d-array + Array containing the boundaries of the tesseroids in the following + order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + Longitudinal and latitudinal boundaries must be in degrees. + The array must have the following shape: (``n_tesseroids``, 6), where + ``n_tesseroids`` is the total number of tesseroids. + This array of tesseroids must have longitude continuity and valid + boundaries. + Run ``_check_tesseroids`` before. + """ + longitude, latitude, radius = coordinates[:] + west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6)) + # Longitudinal boundaries of the tesseroid must be compared with + # longitudinal coordinates of computation points when moved to + # [0, 360) and [-180, 180). + longitude_360 = longitude % 360 + longitude_180 = ((longitude + 180) % 360) - 180 + inside_longitude = np.logical_or( + np.logical_and( + west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east + ), + np.logical_and( + west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east + ), + ) + inside_latitude = np.logical_and( + south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north + ) + inside_radius = np.logical_and( + bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top + ) + # Build array of booleans. + # The (i, j) element is True if the computation point i is inside the + # tesseroid j. + inside = inside_longitude * inside_latitude * inside_radius + if inside.any(): + err_msg = ( + "Found computation point inside tesseroid. " + + "Computation points must be outside of tesseroids.\n" + ) + for point_i, tess_i in np.argwhere(inside): + err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format( + coordinates[:, point_i], tesseroids[tess_i, :] + ) + raise ValueError(err_msg) + + +def _longitude_continuity(tesseroids): + """ + Modify longitudinal boundaries of tesseroids to ensure longitude continuity + + Longitudinal boundaries of the tesseroids are moved to the ``[-180, 180)`` + degrees interval in case the ``west`` boundary is numerically greater than + the ``east`` one. + + Parameters + ---------- + tesseroids : 2d-array + Longitudinal and latitudinal boundaries must be in degrees. + Array containing the boundaries of each tesseroid: + ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric + spherical coordinate system. + The array must have the following shape: (``n_tesseroids``, 6), where + ``n_tesseroids`` is the total number of tesseroids. + + Returns + ------- + tesseroids : 2d-array + Modified boundaries of the tesseroids. + """ + # Copy the tesseroids to avoid modifying the original tesseroids array + tesseroids = tesseroids.copy() + west, east = tesseroids[:, 0], tesseroids[:, 1] + tess_to_be_changed = west > east + east[tess_to_be_changed] = ((east[tess_to_be_changed] + 180) % 360) - 180 + west[tess_to_be_changed] = ((west[tess_to_be_changed] + 180) % 360) - 180 + return tesseroids diff --git a/harmonica/forward/_tesseroid_variable_density.py b/harmonica/forward/_tesseroid_variable_density.py new file mode 100644 index 000000000..3c210da68 --- /dev/null +++ b/harmonica/forward/_tesseroid_variable_density.py @@ -0,0 +1,221 @@ +# Copyright (c) 2018 The Harmonica Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +# +# This code is part of the Fatiando a Terra project (https://www.fatiando.org) +# +""" +Utils functions for tesseroids with variable density +""" +from numba import jit, prange +import numpy as np +from scipy.optimize import minimize_scalar + +DELTA_RATIO = 0.1 + + +@jit(nopython=True) +def gauss_legendre_quadrature_variable_density( + longitude, + cosphi, + sinphi, + radius, + tesseroid, + density, + glq_nodes, + glq_weights, + kernel, +): # pylint: disable=too-many-locals + r""" + Compute the effect of a tesseroid on a single observation point through GLQ + + The tesseroid is converted into a set of point masses located on the + scaled nodes of the Gauss-Legendre Quadrature. The number of point masses + created from each tesseroid is equal to the product of the GLQ degrees for + each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass + of each point mass is defined as the product of the tesseroid density + (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`, + :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and + the :math:`\kappa` factor evaluated on the coordinates of the point mass. + + Parameters + ---------- + longitude : float + Longitudinal coordinate of the observation points in radians. + cosphi : float + Cosine of the latitudinal coordinate of the observation point in + radians. + sinphi : float + Sine of the latitudinal coordinate of the observation point in + radians. + radius : float + Radial coordinate of the observation point in meters. + tesseroids : 1d-array + Array containing the boundaries of the tesseroid: + ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + Horizontal boundaries should be in degrees and radial boundaries in + meters. + density : func + Density func of the tesseroid in SI units. + glq_nodes : list + Unscaled location of GLQ nodes for each direction. + glq_weights : list + GLQ weigths for each node for each direction. + kernel : func + Kernel function for the gravitational field of point masses. + + """ + # Get tesseroid boundaries + w, e, s, n, bottom, top = tesseroid[:] + # Calculate the A factor for the tesseroid + a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom) + # Unpack nodes and weights + lon_nodes, lat_nodes, rad_nodes = glq_nodes[:] + lon_weights, lat_weights, rad_weights = glq_weights[:] + # Compute effect of the tesseroid on the observation point + # by iterating over the location of the point masses + # (move the iteration along the longitudinal nodes to the bottom for + # optimization: reduce the number of times that the trigonometric functions + # are evaluated) + result = 0.0 + for j, lat_node in enumerate(lat_nodes): + # Get the latitude of the point mass + latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s)) + cosphi_p = np.cos(latitude_p) + sinphi_p = np.sin(latitude_p) + for k, rad_node in enumerate(rad_nodes): + # Get the radius of the point mass + radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom) + density_p = density(radius_p) + # Get kappa constant for the point mass + kappa = radius_p ** 2 * cosphi_p + for i, lon_node in enumerate(lon_nodes): + # Get the longitude of the point mass + longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w)) + # Compute the mass of the point mass + mass = ( + density_p + * a_factor + * kappa + * lon_weights[i] + * lat_weights[j] + * rad_weights[k] + ) + # Add effect of the current point mass to the result + result += mass * kernel( + longitude, + cosphi, + sinphi, + radius, + longitude_p, + cosphi_p, + sinphi_p, + radius_p, + ) + return result + + +# Density-based discretization functions +# -------------------------------------- +def density_based_discretization(tesseroids, density): + """ + Apply density_based discretization to a collection of tesseroids + """ + discretized_tesseroids = [] + for tesseroid in tesseroids: + discretized_tesseroids.extend(_density_based_discretization(tesseroid, density)) + return np.atleast_2d(discretized_tesseroids) + + +def _density_based_discretization(tesseroid, density): + """ + Applies density-based discretization to a single tesseroid + + Splits the tesseroid on the points of maximum density variance + + Parameters + ---------- + tesseroid : tuple + density : func + + Returns + ------- + tesseroids : list + """ + # Define normalized density + def normalized_density(radius): + return (density(radius) - density_min) / (density_max - density_min) + + # Get boundaries of original tesseroid + w, e, s, n, bottom, top = tesseroid[:] + # Get minimum and maximum values of the density + density_min, density_max = density_minmax(density, bottom, top) + # Return the original tesseroid if max and min densities are equal + if np.isclose(density_min, density_max): + return [tesseroid] + # Store the size of the original tesseroid + size_original_tesseroid = top - bottom + # Initialize list of pending and output tesseroids + pending, tesseroids = [tesseroid], [] + # Discretization of the tesseroid + while pending: + tesseroid = pending.pop(0) + bottom, top = tesseroid[-2:] + radius_split, max_diff = maximum_absolute_diff(normalized_density, bottom, top) + size_ratio = (top - bottom) / size_original_tesseroid + if max_diff * size_ratio > DELTA_RATIO: + pending.append([w, e, s, n, radius_split, top]) + pending.append([w, e, s, n, bottom, radius_split]) + else: + tesseroids.append([w, e, s, n, bottom, top]) + return tesseroids + + +def density_minmax(density, bottom, top): + """ + Compute the minimum and maximum value of a bounded density + """ + minimum = minimize_scalar(density, bounds=[bottom, top], method="bounded") + maximum = minimize_scalar( + lambda radius: -density(radius), bounds=[bottom, top], method="bounded" + ) + return minimum.fun, -maximum.fun + + +def maximum_absolute_diff(normalized_density, bottom, top): + """ + Compute maximum abs difference between normalized density and straight line + + The maximum difference is computed within the ``bottom`` and ``top`` + boundaries. + """ + + def neg_absolute_difference(radius): + """ + Define minus absolute diff between normalized density and straight line + """ + return -np.abs( + normalized_density(radius) + - straight_line(radius, normalized_density, bottom, top) + ) + + # Use scipy.optimize.minimize_scalar for maximizing the absolute difference + result = minimize_scalar( + neg_absolute_difference, + bounds=[bottom, top], + method="bounded", + ) + # Get maximum difference and the radius at which it takes place + radius_split = result.x + max_diff = -result.fun + return radius_split, max_diff + + +def straight_line(radius, normalized_density, bottom, top): + """ + Compute the reference straight line that joins points of normalized density + """ + norm_density_bottom = normalized_density(bottom) + norm_density_top = normalized_density(top) + slope = (norm_density_top - norm_density_bottom) / (top - bottom) + return slope * (radius - bottom) + norm_density_bottom diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py index 990ed61fc..4227697c5 100644 --- a/harmonica/forward/tesseroid.py +++ b/harmonica/forward/tesseroid.py @@ -9,14 +9,23 @@ """ from numba import jit, prange import numpy as np -from numpy.polynomial.legendre import leggauss from ..constants import GRAVITATIONAL_CONST -from .utils import distance_spherical from .point_mass import ( kernel_potential_spherical, kernel_g_z_spherical, ) +from ._tesseroid_utils import ( + _check_tesseroids, + _check_points_outside_tesseroids, + _adaptive_discretization, + glq_nodes_weights, + gauss_legendre_quadrature, +) +from ._tesseroid_variable_density import ( + density_based_discretization, + gauss_legendre_quadrature_variable_density, +) STACK_SIZE = 100 MAX_DISCRETIZATIONS = 100000 @@ -124,21 +133,27 @@ def tesseroid_gravity( # Convert coordinates, tesseroids and density to arrays coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3]) tesseroids = np.atleast_2d(tesseroids) - density = np.atleast_1d(density).ravel() # Sanity checks for tesseroids and computation points if not disable_checks: - if density.size != tesseroids.shape[0]: + tesseroids = _check_tesseroids(tesseroids) + _check_points_outside_tesseroids(coordinates, tesseroids) + # Check if density is a function or constant values + if callable(density): + # Compile the density function and run density-based discretization + density = jit(nopython=True)(density) + tesseroids = density_based_discretization(tesseroids, density) + else: + density = np.atleast_1d(density).ravel() + if not disable_checks and density.size != tesseroids.shape[0]: raise ValueError( "Number of elements in density ({}) ".format(density.size) + "mismatch the number of tesseroids ({})".format(tesseroids.shape[0]) ) - tesseroids = _check_tesseroids(tesseroids) - _check_points_outside_tesseroids(coordinates, tesseroids) # Get GLQ unscaled nodes, weights and number of nodes for each small # tesseroid glq_nodes, glq_weights = glq_nodes_weights(GLQ_DEGREES) # Compute gravitational field - dispatcher(parallel)( + dispatcher(parallel, density)( coordinates, tesseroids, density, @@ -157,14 +172,23 @@ def tesseroid_gravity( return result.reshape(cast.shape) -def dispatcher(parallel): +def dispatcher(parallel, density): """ - Return the parallelized or serialized forward modelling function + Return the jitted compiled forward modelling function + + The choice of the forward modelling function is based on whether the + density is a function and if the model should be run in parallel. """ - dispatchers = { - True: jit_tesseroid_gravity_parallel, - False: jit_tesseroid_gravity_serial, - } + if callable(density): + dispatchers = { + True: jit_tesseroid_gravity_variable_density_parallel, + False: jit_tesseroid_gravity_variable_density_serial, + } + else: + dispatchers = { + True: jit_tesseroid_gravity_parallel, + False: jit_tesseroid_gravity_serial, + } return dispatchers[parallel] @@ -207,12 +231,6 @@ def jit_tesseroid_gravity( should be in meters. density : 1d-array Density of each tesseroid in SI units. - stack : 2d-array - Empty array where tesseroids created by adaptive discretization - algorithm will be processed. - small_tesseroids : 2d-array - Empty array where smaller tesseroids created by adaptive discretization - algorithm will be stored. result : 1d-array Array where the gravitational effect of each tesseroid will be added. distance_size_ratio : float @@ -269,489 +287,99 @@ def jit_tesseroid_gravity( ) -@jit(nopython=True) -def gauss_legendre_quadrature( - longitude, - cosphi, - sinphi, - radius, - tesseroid, +def jit_tesseroid_gravity_variable_density( + coordinates, + tesseroids, density, + result, + distance_size_ratio, + radial_adaptive_discretization, glq_nodes, glq_weights, kernel, -): # pylint: disable=too-many-locals - r""" - Compute the effect of a tesseroid on a single observation point through GLQ - - The tesseroid is converted into a set of point masses located on the - scaled nodes of the Gauss-Legendre Quadrature. The number of point masses - created from each tesseroid is equal to the product of the GLQ degrees for - each direction (:math:`N_r`, :math:`N_\lambda`, :math:`N_\phi`). The mass - of each point mass is defined as the product of the tesseroid density - (:math:`\rho`), the GLQ weights for each direction (:math:`W_i^r`, - :math:`W_j^\phi`, :math:`W_k^\lambda`), the scale constant :math:`A` and - the :math:`\kappa` factor evaluated on the coordinates of the point mass. - - Parameters - ---------- - longitude : float - Longitudinal coordinate of the observation points in radians. - cosphi : float - Cosine of the latitudinal coordinate of the observation point in - radians. - sinphi : float - Sine of the latitudinal coordinate of the observation point in - radians. - radius : float - Radial coordinate of the observation point in meters. - tesseroids : 1d-array - Array containing the boundaries of the tesseroid: - ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. - Horizontal boundaries should be in degrees and radial boundaries in - meters. - density : float - Density of the tesseroid in SI units. - glq_nodes : list - Unscaled location of GLQ nodes for each direction. - glq_weights : list - GLQ weigths for each node for each direction. - kernel : func - Kernel function for the gravitational field of point masses. - - """ - # Get tesseroid boundaries - w, e, s, n, bottom, top = tesseroid[:] - # Calculate the A factor for the tesseroid - a_factor = 1 / 8 * np.radians(e - w) * np.radians(n - s) * (top - bottom) - # Unpack nodes and weights - lon_nodes, lat_nodes, rad_nodes = glq_nodes[:] - lon_weights, lat_weights, rad_weights = glq_weights[:] - # Compute effect of the tesseroid on the observation point - # by iterating over the location of the point masses - # (move the iteration along the longitudinal nodes to the bottom for - # optimization: reduce the number of times that the trigonometric functions - # are evaluated) - result = 0.0 - for j, lat_node in enumerate(lat_nodes): - # Get the latitude of the point mass - latitude_p = np.radians(0.5 * (n - s) * lat_node + 0.5 * (n + s)) - cosphi_p = np.cos(latitude_p) - sinphi_p = np.sin(latitude_p) - for k, rad_node in enumerate(rad_nodes): - # Get the radius of the point mass - radius_p = 0.5 * (top - bottom) * rad_node + 0.5 * (top + bottom) - # Get kappa constant for the point mass - kappa = radius_p ** 2 * cosphi_p - for i, lon_node in enumerate(lon_nodes): - # Get the longitude of the point mass - longitude_p = np.radians(0.5 * (e - w) * lon_node + 0.5 * (e + w)) - # Compute the mass of the point mass - mass = ( - density - * a_factor - * kappa - * lon_weights[i] - * lat_weights[j] - * rad_weights[k] - ) - # Add effect of the current point mass to the result - result += mass * kernel( - longitude, - cosphi, - sinphi, - radius, - longitude_p, - cosphi_p, - sinphi_p, - radius_p, - ) - return result - - -def glq_nodes_weights(glq_degrees): - """ - Calculate GLQ unscaled nodes and weights - - Parameters - ---------- - glq_degrees : list - List of GLQ degrees for each direction: ``longitude``, ``latitude``, - ``radius``. - - Returns - ------- - glq_nodes : list - Unscaled GLQ nodes for each direction: ``longitude``, ``latitude``, - ``radius``. - glq_weights : list - GLQ weights for each node on each direction: ``longitude``, - ``latitude``, ``radius``. - """ - # Unpack GLQ degrees - lon_degree, lat_degree, rad_degree = glq_degrees[:] - # Get nodes coordinates and weights - lon_node, lon_weights = leggauss(lon_degree) - lat_node, lat_weights = leggauss(lat_degree) - rad_node, rad_weights = leggauss(rad_degree) - # Reorder nodes and weights - glq_nodes = (lon_node, lat_node, rad_node) - glq_weights = (lon_weights, lat_weights, rad_weights) - return glq_nodes, glq_weights - - -@jit(nopython=True) -def _adaptive_discretization( - coordinates, - tesseroid, - distance_size_ratio, - stack, - small_tesseroids, - radial_discretization=False, -): - """ - Perform the adaptive discretization algorithm on a tesseroid - - It apply the three or two dimensional adaptive discretization algorithm on - a tesseroid after a single computation point. - - Parameters - ---------- - coordinates : array - Array containing ``longitude``, ``latitude`` and ``radius`` of a single - computation point. - tesseroid : array - Array containing the boundaries of the tesseroid. - distance_size_ratio : float - Value for the distance-size ratio. A greater value will perform more - discretizations. - stack : 2d-array - Array with shape ``(6, stack_size)`` that will temporarly hold the - small tesseroids that are not yet processed. - If too many discretizations will take place, increase the - ``stack_size``. - small_tesseroids : 2d-array - Array with shape ``(6, MAX_DISCRETIZATIONS)`` that will contain every - small tesseroid produced by the adaptive discretization algorithm. - If too many discretizations will take place, increase the - ``MAX_DISCRETIZATIONS``. - radial_discretization : bool (optional) - If ``True`` the three dimensional adaptive discretization will be - applied. - If ``False`` the two dimensional adaptive discretization will be - applied, i.e. the tesseroid will only be split on the ``longitude`` and - ``latitude`` directions. - Default ``False``. - - Returns - ------- - n_splits : int - Total number of small tesseroids generated by the algorithm. - """ - # Create stack of tesseroids - stack[0] = tesseroid - stack_top = 0 - n_splits = 0 - while stack_top >= 0: - # Pop the first tesseroid from the stack - tesseroid = stack[stack_top] - stack_top -= 1 - # Get its dimensions - l_lon, l_lat, l_rad = _tesseroid_dimensions(tesseroid) - # Get distance between computation point and center of tesseroid - distance = _distance_tesseroid_point(coordinates, tesseroid) - # Check inequality - n_lon, n_lat, n_rad = 1, 1, 1 - if distance / l_lon < distance_size_ratio: - n_lon = 2 - if distance / l_lat < distance_size_ratio: - n_lat = 2 - if distance / l_rad < distance_size_ratio and radial_discretization: - n_rad = 2 - # Apply discretization - if n_lon * n_lat * n_rad > 1: - # Raise error if stack overflow - # Number of tesseroids in stack = stack_top + 1 - if (stack_top + 1) + n_lon * n_lat * n_rad > stack.shape[0]: - raise OverflowError("Stack Overflow. Try to increase the stack size.") - stack_top = _split_tesseroid( - tesseroid, n_lon, n_lat, n_rad, stack, stack_top - ) - else: - # Raise error if small_tesseroids overflow - if n_splits + 1 > small_tesseroids.shape[0]: - raise OverflowError( - "Exceeded maximum discretizations." - + " Please increase the MAX_DISCRETIZATIONS." - ) - small_tesseroids[n_splits] = tesseroid - n_splits += 1 - return n_splits - - -@jit(nopython=True) -def _split_tesseroid( - tesseroid, n_lon, n_lat, n_rad, stack, stack_top -): # pylint: disable=too-many-locals - """ - Split tesseroid along each dimension - """ - w, e, s, n, bottom, top = tesseroid[:] - # Compute differential distance - d_lon = (e - w) / n_lon - d_lat = (n - s) / n_lat - d_rad = (top - bottom) / n_rad - for i in range(n_lon): - for j in range(n_lat): - for k in range(n_rad): - stack_top += 1 - stack[stack_top, 0] = w + d_lon * i - stack[stack_top, 1] = w + d_lon * (i + 1) - stack[stack_top, 2] = s + d_lat * j - stack[stack_top, 3] = s + d_lat * (j + 1) - stack[stack_top, 4] = bottom + d_rad * k - stack[stack_top, 5] = bottom + d_rad * (k + 1) - return stack_top - - -@jit(nopython=True) -def _tesseroid_dimensions(tesseroid): - """ - Calculate the dimensions of the tesseroid. - """ - w, e, s, n, bottom, top = tesseroid[:] - w, e, s, n = np.radians(w), np.radians(e), np.radians(s), np.radians(n) - latitude_center = (n + s) / 2 - l_lat = top * np.arccos(np.sin(n) * np.sin(s) + np.cos(n) * np.cos(s)) - l_lon = top * np.arccos( - np.sin(latitude_center) ** 2 + np.cos(latitude_center) ** 2 * np.cos(e - w) - ) - l_rad = top - bottom - return l_lon, l_lat, l_rad - - -@jit(nopython=True) -def _distance_tesseroid_point( - coordinates, tesseroid -): # pylint: disable=too-many-locals - """ - Distance between a computation point and the center of a tesseroid - """ - # Get center of the tesseroid - w, e, s, n, bottom, top = tesseroid[:] - longitude_p = (w + e) / 2 - latitude_p = (s + n) / 2 - radius_p = (bottom + top) / 2 - # Get distance between computation point and tesseroid center - distance = distance_spherical(coordinates, (longitude_p, latitude_p, radius_p)) - return distance - - -def _check_tesseroids(tesseroids): # pylint: disable=too-many-branches - """ - Check if tesseroids boundaries are well defined - - A valid tesseroid should have: - - latitudinal boundaries within the [-90, 90] degrees interval, - - north boundaries greater or equal than the south boundaries, - - radial boundaries positive or zero, - - top boundaries greater or equal than the bottom boundaries, - - longitudinal boundaries within the [-180, 360] degrees interval, - - longitudinal interval must not be greater than one turn around the - globe. - - Some valid tesseroids have its west boundary greater than the east one, - e.g. ``(350, 10, ...)``. On these cases the ``_longitude_continuity`` - function is applied in order to move the longitudinal coordinates to the - [-180, 180) interval. Any valid tesseroid should have east boundaries - greater than the west boundaries before or after applying longitude - continuity. - - Parameters - ---------- - tesseroids : 2d-array - Array containing the boundaries of the tesseroids in the following - order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. - Longitudinal and latitudinal boundaries must be in degrees. - The array must have the following shape: (``n_tesseroids``, 6), where - ``n_tesseroids`` is the total number of tesseroids. - - Returns - ------- - tesseroids : 2d-array - Array containing the boundaries of the tesseroids. - If no longitude continuity needs to be applied, the returned array is - the same one as the orignal. - Otherwise, it's copied and its longitudinal boundaries are modified. - """ - west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6)) - err_msg = "Invalid tesseroid or tesseroids. " - # Check if latitudinal boundaries are inside the [-90, 90] interval - invalid = np.logical_or( - np.logical_or(south < -90, south > 90), np.logical_or(north < -90, north > 90) - ) - if (invalid).any(): - err_msg += ( - "The latitudinal boundaries must be inside the [-90, 90] " - + "degrees interval.\n" - ) - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - # Check if south boundary is not greater than the corresponding north - # boundary - invalid = south > north - if (invalid).any(): - err_msg += "The south boundary can't be greater than the north one.\n" - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - # Check if radial boundaries are positive or zero - invalid = np.logical_or(bottom < 0, top < 0) - if (invalid).any(): - err_msg += "The bottom and top radii should be positive or zero.\n" - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - # Check if top boundary is not greater than the corresponding bottom - # boundary - invalid = bottom > top - if (invalid).any(): - err_msg += "The bottom radius boundary can't be greater than the top one.\n" - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - # Check if longitudinal boundaries are inside the [-180, 360] interval - invalid = np.logical_or( - np.logical_or(west < -180, west > 360), np.logical_or(east < -180, east > 360) - ) - if (invalid).any(): - err_msg += ( - "The longitudinal boundaries must be inside the [-180, 360] " - + "degrees interval.\n" - ) - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - # Apply longitude continuity if w > e - if (west > east).any(): - tesseroids = _longitude_continuity(tesseroids) - west, east, south, north, bottom, top = tuple( - tesseroids[:, i] for i in range(6) - ) - # Check if west boundary is not greater than the corresponding east - # boundary, even after applying the longitude continuity - invalid = west > east - if (invalid).any(): - err_msg += "The west boundary can't be greater than the east one.\n" - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - # Check if the longitudinal interval is not grater than one turn around the - # globe - invalid = east - west > 360 - if (invalid).any(): - err_msg += ( - "The difference between east and west boundaries cannot be greater than " - + "one turn around the globe.\n" - ) - for tess in tesseroids[invalid]: - err_msg += "\tInvalid tesseroid: {}\n".format(tess) - raise ValueError(err_msg) - return tesseroids - - -def _check_points_outside_tesseroids( - coordinates, tesseroids -): # pylint: disable=too-many-locals - """ - Check if computation points are not inside the tesseroids - - Parameters - ---------- - coordinates : 2d-array - Array containing the coordinates of the computation points in the - following order: ``longitude``, ``latitude`` and ``radius``. - Both ``longitude`` and ``latitude`` must be in degrees. - The array must have the following shape: (3, ``n_points``), where - ``n_points`` is the total number of computation points. - tesseroids : 2d-array - Array containing the boundaries of the tesseroids in the following - order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. - Longitudinal and latitudinal boundaries must be in degrees. - The array must have the following shape: (``n_tesseroids``, 6), where - ``n_tesseroids`` is the total number of tesseroids. - This array of tesseroids must have longitude continuity and valid - boundaries. - Run ``_check_tesseroids`` before. - """ - longitude, latitude, radius = coordinates[:] - west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6)) - # Longitudinal boundaries of the tesseroid must be compared with - # longitudinal coordinates of computation points when moved to - # [0, 360) and [-180, 180). - longitude_360 = longitude % 360 - longitude_180 = ((longitude + 180) % 360) - 180 - inside_longitude = np.logical_or( - np.logical_and( - west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east - ), - np.logical_and( - west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east - ), - ) - inside_latitude = np.logical_and( - south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north - ) - inside_radius = np.logical_and( - bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top - ) - # Build array of booleans. - # The (i, j) element is True if the computation point i is inside the - # tesseroid j. - inside = inside_longitude * inside_latitude * inside_radius - if inside.any(): - err_msg = ( - "Found computation point inside tesseroid. " - + "Computation points must be outside of tesseroids.\n" - ) - for point_i, tess_i in np.argwhere(inside): - err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format( - coordinates[:, point_i], tesseroids[tess_i, :] - ) - raise ValueError(err_msg) - - -def _longitude_continuity(tesseroids): + dtype, +): # pylint: disable=too-many-locals,invalid-name,not-an-iterable """ - Modify longitudinal boundaries of tesseroids to ensure longitude continuity + Compute gravitational field of tesseroids on computations points - Longitudinal boundaries of the tesseroids are moved to the ``[-180, 180)`` - degrees interval in case the ``west`` boundary is numerically greater than - the ``east`` one. + Perform adaptive discretization, convert each small tesseroid to equivalent + point masses through GLQ and use point masses kernel functions to compute + the gravitational field. Parameters ---------- + coordinates : tuple + Tuple containing the coordinates of the computation points in spherical + geocentric coordinate system in the following order: + ``longitude``, ``latitude``, ``radius``. + Each element of the tuple must be a 1d array. + Both ``longitude`` and ``latitude`` should be in degrees and ``radius`` + in meters. tesseroids : 2d-array - Longitudinal and latitudinal boundaries must be in degrees. Array containing the boundaries of each tesseroid: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top`` under a geocentric spherical coordinate system. The array must have the following shape: (``n_tesseroids``, 6), where ``n_tesseroids`` is the total number of tesseroids. - - Returns - ------- - tesseroids : 2d-array - Modified boundaries of the tesseroids. + All tesseroids must have valid boundary coordinates. + Horizontal boundaries should be in degrees while radial boundaries + should be in meters. + density : func + Density function of every tesseroid in SI units. + result : 1d-array + Array where the gravitational effect of each tesseroid will be added. + distance_size_ratio : float + Value of the distance size ratio. + radial_adaptive_discretization : bool + If ``False``, the adaptive discretization algorithm will split the + tesseroid only on the horizontal direction. + If ``True``, it will perform a three dimensional adaptive + discretization, splitting the tesseroids on every direction. + glq_nodes : list + List containing unscaled GLQ nodes. + glq_weights : list + List containing GLQ weights of the nodes. + kernel : func + Kernel function for the gravitational field of point masses. + dtype : data-type + Data type assigned to the resulting gravitational field. """ - # Copy the tesseroids to avoid modifying the original tesseroids array - tesseroids = tesseroids.copy() - west, east = tesseroids[:, 0], tesseroids[:, 1] - tess_to_be_changed = west > east - east[tess_to_be_changed] = ((east[tess_to_be_changed] + 180) % 360) - 180 - west[tess_to_be_changed] = ((west[tess_to_be_changed] + 180) % 360) - 180 - return tesseroids + # Get coordinates of the observation points + # and precompute trigonometric functions + longitude, latitude, radius = coordinates[:] + longitude_rad = np.radians(longitude) + cosphi = np.cos(np.radians(latitude)) + sinphi = np.sin(np.radians(latitude)) + # Loop over computation points + for l in prange(longitude.size): + # Initialize arrays to perform memory allocation only once + stack = np.empty((STACK_SIZE, 6), dtype=dtype) + small_tesseroids = np.empty((MAX_DISCRETIZATIONS, 6), dtype=dtype) + # Loop over tesseroids + for m in range(tesseroids.shape[0]): + # Apply adaptive discretization on tesseroid + n_splits = _adaptive_discretization( + (longitude[l], latitude[l], radius[l]), + tesseroids[m, :], + distance_size_ratio, + stack, + small_tesseroids, + radial_adaptive_discretization, + ) + # Compute effect of the tesseroid through GLQ + for tess_index in range(n_splits): + tesseroid = small_tesseroids[tess_index, :] + result[l] += gauss_legendre_quadrature_variable_density( + longitude_rad[l], + cosphi[l], + sinphi[l], + radius[l], + tesseroid, + density, + glq_nodes, + glq_weights, + kernel, + ) # Define jitted versions of the forward modelling function @@ -760,3 +388,9 @@ def _longitude_continuity(tesseroids): jit_tesseroid_gravity_parallel = jit(nopython=True, parallel=True)( jit_tesseroid_gravity ) +jit_tesseroid_gravity_variable_density_serial = jit(nopython=True)( + jit_tesseroid_gravity_variable_density +) +jit_tesseroid_gravity_variable_density_parallel = jit(nopython=True, parallel=True)( + jit_tesseroid_gravity_variable_density +) From ea7f731a88def26f625064095d7638fdaee283da Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 22 Oct 2021 13:47:37 -0300 Subject: [PATCH 03/17] Fix import on tests --- harmonica/tests/test_tesseroid.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index 5a4ddf260..1b69ab2a8 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -19,15 +19,18 @@ tesseroid_gravity, _check_tesseroids, _check_points_outside_tesseroids, + _adaptive_discretization, + STACK_SIZE, + MAX_DISCRETIZATIONS, +) +from ..forward._tesseroid_utils import ( _distance_tesseroid_point, _tesseroid_dimensions, _split_tesseroid, - _adaptive_discretization, _longitude_continuity, - STACK_SIZE, - MAX_DISCRETIZATIONS, ) + # Define the accuracy threshold for tesseroids (0.1%) as a # relative error (0.001) ACCURACY_THRESHOLD = 1e-3 From 1abe841a9d821a0d97498480c58745ad7aaf8667 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 22 Oct 2021 18:04:12 -0300 Subject: [PATCH 04/17] Start writing test functions for variable density tesseroids --- .../tests/test_tesseroid_variable_density.py | 261 ++++++++++++++++++ 1 file changed, 261 insertions(+) create mode 100644 harmonica/tests/test_tesseroid_variable_density.py diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py new file mode 100644 index 000000000..1124dc802 --- /dev/null +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -0,0 +1,261 @@ +# Copyright (c) 2018 The Harmonica Developers. +# Distributed under the terms of the BSD 3-Clause License. +# SPDX-License-Identifier: BSD-3-Clause +# +# This code is part of the Fatiando a Terra project (https://www.fatiando.org) +# +""" +Test forward modelling for tesseroids with variable density +""" +import numpy as np +import numpy.testing as npt +import harmonica +import pytest + +from ..forward._tesseroid_variable_density import ( + straight_line, + maximum_absolute_diff, + density_minmax, + _density_based_discretization, +) + + +# Define the accuracy threshold for tesseroids (0.1%) as a +# relative error (0.001) +ACCURACY_THRESHOLD = 1e-3 + + +@pytest.fixture(name="bottom") +def fixture_bottom(): + bottom = 2e3 + return bottom + + +@pytest.fixture(name="top") +def fixture_top(): + top = 5e3 + return top + + +@pytest.fixture(name="quadratic_params") +def fixture_quadratic_params(): + factor = 1e-3 + vertex_radius = 3e3 + vertex_density = 1900.0 + return factor, vertex_radius, vertex_density + + +@pytest.fixture(name="quadratic_density") +def fixture_quadratic_density(quadratic_params): + """ + Return a quadratic density function + + .. math:: + + f(x) = a (x - h)^2 + k + + Where :math:`a` is the ``factor``, :math:`h` is the ``vertex_radius``, + :math:`k` is the ``vertex_density`` and :math:`x` is the ``radius``. + + :: + + +------------------------------------------------+ + | *| + | * | + | * | + | ** | + | * | + | ** | + | * | + | ** | + | ** | + | ** | + |* ** | + | ** *** | + | **** *** | + | ****** ******* | + | ***** | + +------------------------------------------------+ + | | | + bottom vertex_radius top + + """ + factor, vertex_radius, vertex_density = quadratic_params + + def density(radius): + """Quadratic density function""" + return factor * (radius - vertex_radius) ** 2 + vertex_density + + return density + + +@pytest.fixture(name="straight_line_analytic") +def fixture_straight_line_analytic(bottom, top, quadratic_density): + """ + Return the analytic solution for the straigh line of the quadratic density + """ + density_bottom, density_top = quadratic_density(bottom), quadratic_density(top) + slope = (density_top - density_bottom) / (top - bottom) + + def line(radius): + return slope * (radius - bottom) + density_bottom + + return line + + +@pytest.fixture(name="max_abs_diff_analytic") +def fixture_max_abs_diff_analytic( + bottom, top, quadratic_params, quadratic_density, straight_line_analytic +): + r""" + Return analytic solution for maximum absolute difference between quadratic + density and the straight line + + .. math: + + x_m = \frac{m}{2a} + h + + """ + factor, vertex_radius, _ = quadratic_params + density_bottom, density_top = quadratic_density(bottom), quadratic_density(top) + slope = (density_top - density_bottom) / (top - bottom) + radius_split = 0.5 * slope / factor + vertex_radius + max_diff = np.abs( + quadratic_density(radius_split) - straight_line_analytic(radius_split) + ) + return radius_split, max_diff + + +@pytest.fixture(name="quadratic_density_minmax") +def fixture_quadratic_density_minmax(top, quadratic_params, quadratic_density): + """ + Return the analytic maximum and minimum value of the quadratic density + between top and bottom + """ + _, _, vertex_density = quadratic_params + minimum = vertex_density + maximum = quadratic_density(top) + return minimum, maximum + + +def test_straight_line(bottom, top, quadratic_density, straight_line_analytic): + """ + Test the straight_line function + """ + radii = np.linspace(bottom, top, 51) + npt.assert_allclose( + straight_line(radii, quadratic_density, bottom, top), + straight_line_analytic(radii), + ) + + +def test_max_abs_diff(bottom, top, quadratic_density, max_abs_diff_analytic): + """ + Test the maximum absolute difference + + Test against the sine density defined in density_sine_portion fixture. + The solution is analytic. + """ + radius_split_expected, max_diff_expected = max_abs_diff_analytic + radius_split, max_diff = maximum_absolute_diff(quadratic_density, bottom, top) + npt.assert_allclose(radius_split_expected, radius_split) + npt.assert_allclose(max_diff_expected, max_diff) + + +def test_density_minmax(bottom, top, quadratic_density, quadratic_density_minmax): + """ + Test the density_minmax function + """ + density_min, density_max = density_minmax(quadratic_density, bottom, top) + expected_min, expected_max = quadratic_density_minmax + npt.assert_allclose(density_min, expected_min) + npt.assert_allclose(density_max, expected_max) + + +def test_single_density_based_discretization( + bottom, top, quadratic_density, max_abs_diff_analytic +): + """ + Test the density-based discretization algorithm + """ + # Define some dummy horizontal coordinates for the tesseroid + w, e, s, n = -3, 2, -4, 5 + tesseroid = w, e, s, n, bottom, top + tesseroids = _density_based_discretization(tesseroid, quadratic_density) + # With the default value of DELTA=0.1, it should produce only 2 tesseroids + assert len(tesseroids) == 2 + # Check the horizontal coordinates of the tesseroids + for tess in tesseroids: + for coord, original_coord in zip(tess[:4], tesseroid): + npt.assert_allclose(coord, original_coord) + + # Check the radial coordinates + # ---------------------------- + # Take the analytical solution for the radius_split + radius_split, _ = max_abs_diff_analytic + # The first tesseroid in the list should be the one in the top + npt.assert_allclose(tesseroids[0][-2], radius_split) + npt.assert_allclose(tesseroids[0][-1], top) + # The second tesseroid in the list should be the one in the bottom + npt.assert_allclose(tesseroids[1][-2], bottom) + npt.assert_allclose(tesseroids[1][-1], radius_split) + + +def test_density_based_discret_with_delta( + bottom, + top, + quadratic_density, +): + """ + Test the density-based discretization algorithm against values of DELTA + """ + # Define some dummy horizontal coordinates for the tesseroid + w, e, s, n = -3, 2, -4, 5 + tesseroid = w, e, s, n, bottom, top + # Define a collection of values for the delta ratio + deltas = [1e-3, 1e-2, 1e-1, 1e0] + # Define an empty list to count the number of splits for each delta + splits = [] + for delta in deltas: + # Override the DELTA_RATIO global variable + harmonica.forward._tesseroid_variable_density.DELTA_RATIO = delta + # Count the splits generated by density based discretization + splits.append(len(_density_based_discretization(tesseroid, quadratic_density))) + splits = np.array(splits) + # Check if numbers of splits increases for lower deltas + assert (splits[1:] < splits[:-1]).all() + + +def test_density_based_discret_linear_density(): + """ + Test if density-based discretization generates no splits when linear + density is passed + """ + w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50 + tesseroid = [w, e, s, n, bottom, top] + + def linear_density(radius): + """Define a dummy linear density""" + return 3 * radius + 2 + + tesseroids = _density_based_discretization(tesseroid, linear_density) + assert len(tesseroids) == 1 + npt.assert_allclose(tesseroids[0], tesseroid) + + +def test_density_based_discret_constant_density(): + """ + Test if density-based discretization generates no splits when a constant + density function is passed (this should not be done IRL, pass an array of + floats as density instead) + """ + w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50 + tesseroid = [w, e, s, n, bottom, top] + + def stupid_constant_density(radius): + """Define a dummy constant density function""" + return 3 + + tesseroids = _density_based_discretization(tesseroid, stupid_constant_density) + assert len(tesseroids) == 1 + npt.assert_allclose(tesseroids[0], tesseroid) From e1433ff4e97b0bebe04b0ff5eda4f162fc9a36a4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 10:29:37 -0300 Subject: [PATCH 05/17] Add test function for spherical shell w linear density --- .../tests/test_tesseroid_variable_density.py | 101 ++++++++++++++++++ 1 file changed, 101 insertions(+) diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 1124dc802..55c40f92c 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -11,7 +11,11 @@ import numpy.testing as npt import harmonica import pytest +from verde import grid_coordinates +from .utils import require_numba +from ..constants import GRAVITATIONAL_CONST +from .. import tesseroid_gravity from ..forward._tesseroid_variable_density import ( straight_line, maximum_absolute_diff, @@ -25,6 +29,11 @@ ACCURACY_THRESHOLD = 1e-3 +# --------------------------------- +# Test density-based discretization +# --------------------------------- + + @pytest.fixture(name="bottom") def fixture_bottom(): bottom = 2e3 @@ -259,3 +268,95 @@ def stupid_constant_density(radius): tesseroids = _density_based_discretization(tesseroid, stupid_constant_density) assert len(tesseroids) == 1 npt.assert_allclose(tesseroids[0], tesseroid) + + +# ---------------------------- +# Test against spherical shell +# ---------------------------- + + +def spherical_shell_linear( + radius, + bottom, + top, + slope, + constant_term, +): + """ + Analytical solutions of a spherical shell with linear density + """ + constant = np.pi * GRAVITATIONAL_CONST * slope * ( + top ** 4 - bottom ** 4 + ) + 4 / 3.0 * np.pi * GRAVITATIONAL_CONST * constant_term * (top ** 3 - bottom ** 3) + potential = constant / radius + data = { + "potential": potential, + "g_z": 1e5 * (potential / radius), + } + return data + + +def build_spherical_shell(bottom, top, shape=(6, 12)): + """ + Return a set of tesseroids modelling a spherical shell + + Parameters + ---------- + bottom : float + Inner radius of the spherical shell + top : float + Outer radius of the spherical shell + shape : tuple (n_latitude, n_longitude) + """ + region = (-180, 180, -90, 90) + n_lat, n_lon = shape[:] + longitude = np.linspace(*region[:2], n_lon + 1, dtype="float64") + latitude = np.linspace(*region[2:], n_lat + 1, dtype="float64") + west, east = longitude[:-1], longitude[1:] + south, north = latitude[:-1], latitude[1:] + tesseroids = [] + for w, e in zip(west, east): + for s, n in zip(south, north): + tesseroids.append([w, e, s, n, bottom, top]) + return tesseroids + + +@require_numba +@pytest.mark.use_numba +@pytest.mark.parametrize("field", ("potential", "g_z")) +@pytest.mark.parametrize("thickness", (1e2, 1e3, 1e4, 1e5, 1e6)) +def test_spherical_shell_linear_density(field, thickness): + """ + Test numerical results against analytical solution for linear density + """ + # Define a spherical shell made out of tesseroids + top = 6371e3 + bottom = top - thickness + tesseroids = build_spherical_shell(bottom, top) + # Create a set of observation points located on the shell outer radius + region = (-180, 180, -90, 90) + shape = (19, 13) + coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top) + + # Define a linear density + density_outer, density_inner = 2500.0, 3300.0 + slope = (density_outer - density_inner) / (top - bottom) + constant_term = density_outer - slope * top + + def linear_density(radius): + """ + Create a dummy linear density + """ + return slope * radius + constant_term + + # Get analytic solution + analytic = spherical_shell_linear( + coordinates[-1], bottom, top, slope, constant_term + ) + + # Check if numerical results are accurate enough + npt.assert_allclose( + tesseroid_gravity(coordinates, tesseroids, linear_density, field=field), + analytic[field], + rtol=ACCURACY_THRESHOLD, + ) From 38f2374321d96a0bdea3a378085d19884ab02d3b Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 15:02:30 -0300 Subject: [PATCH 06/17] Add tests for shell with exponential density --- .../tests/test_tesseroid_variable_density.py | 127 +++++++++++++++--- 1 file changed, 112 insertions(+), 15 deletions(-) diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 55c40f92c..03c17daa6 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -181,6 +181,32 @@ def test_density_minmax(bottom, top, quadratic_density, quadratic_density_minmax npt.assert_allclose(density_max, expected_max) +def test_density_minmax_exponential_function(bottom, top): + """ + Test density_minmax with an extreme function + + A previous implementation of density_minmax failed this test + """ + # Define an expontential density that has a strong L shape + thickness = top - bottom + b_factor = 50 + + def exponential_density(radius): + """ + Create a dummy exponential density + """ + density_outer, density_inner = 2500.0, 3300.0 + a_factor = (density_inner - density_outer) / (1 - np.exp(-b_factor)) + constant_term = density_inner - a_factor + return ( + a_factor * np.exp(-b_factor * (radius - bottom) / thickness) + constant_term + ) + + density_min, density_max = density_minmax(exponential_density, bottom, top) + npt.assert_allclose(density_min, exponential_density(top)) + npt.assert_allclose(density_max, exponential_density(bottom)) + + def test_single_density_based_discretization( bottom, top, quadratic_density, max_abs_diff_analytic ): @@ -275,13 +301,7 @@ def stupid_constant_density(radius): # ---------------------------- -def spherical_shell_linear( - radius, - bottom, - top, - slope, - constant_term, -): +def analytical_spherical_shell_linear(radius, bottom, top, slope, constant_term): """ Analytical solutions of a spherical shell with linear density """ @@ -296,6 +316,43 @@ def spherical_shell_linear( return data +def analytical_spherical_shell_exponential( + radius, bottom, top, a_factor, b_factor, constant_term +): + r""" + Analytical solutions of a spherical shell with exponential density + + .. math : + \rho(r') = a e^{- b * (r' - R_1) / T} + c + """ + thickness = top - bottom + k = b_factor / thickness + potential = ( + 4 + * np.pi + * GRAVITATIONAL_CONST + * a_factor + / k ** 3 + / radius + * ( + ((bottom * k) ** 2 + 2 * bottom * k + 2) + - ((top * k) ** 2 + 2 * top * k + 2) * np.exp(-k * thickness) + ) + + 4 + / 3 + * np.pi + * GRAVITATIONAL_CONST + * constant_term + * (top ** 3 - bottom ** 3) + / radius + ) + data = { + "potential": potential, + "g_z": 1e5 * (potential / radius), + } + return data + + def build_spherical_shell(bottom, top, shape=(6, 12)): """ Return a set of tesseroids modelling a spherical shell @@ -324,7 +381,7 @@ def build_spherical_shell(bottom, top, shape=(6, 12)): @require_numba @pytest.mark.use_numba @pytest.mark.parametrize("field", ("potential", "g_z")) -@pytest.mark.parametrize("thickness", (1e2, 1e3, 1e4, 1e5, 1e6)) +@pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) def test_spherical_shell_linear_density(field, thickness): """ Test numerical results against analytical solution for linear density @@ -333,11 +390,6 @@ def test_spherical_shell_linear_density(field, thickness): top = 6371e3 bottom = top - thickness tesseroids = build_spherical_shell(bottom, top) - # Create a set of observation points located on the shell outer radius - region = (-180, 180, -90, 90) - shape = (19, 13) - coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top) - # Define a linear density density_outer, density_inner = 2500.0, 3300.0 slope = (density_outer - density_inner) / (top - bottom) @@ -349,14 +401,59 @@ def linear_density(radius): """ return slope * radius + constant_term + # Create a set of observation points located on the shell outer radius + region = (-180, 180, -90, 90) + shape = (19, 13) + coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top) # Get analytic solution - analytic = spherical_shell_linear( + analytic = analytical_spherical_shell_linear( coordinates[-1], bottom, top, slope, constant_term ) - # Check if numerical results are accurate enough npt.assert_allclose( tesseroid_gravity(coordinates, tesseroids, linear_density, field=field), analytic[field], rtol=ACCURACY_THRESHOLD, ) + + +@require_numba +@pytest.mark.use_numba +@pytest.mark.parametrize("field", ("potential", "g_z")) +@pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) +@pytest.mark.parametrize("b_factor", (1, 2, 5, 10, 30, 100)) +def test_spherical_shell_exponential_density(field, thickness, b_factor): + """ + Test numerical results against analytical solution for exponential density + """ + # Define a spherical shell made out of tesseroids + top = 6371e3 + bottom = top - thickness + tesseroids = build_spherical_shell(bottom, top) + # Define a linear density + density_outer, density_inner = 2670.0, 3300.0 + a_factor = (density_inner - density_outer) / (1 - np.exp(-b_factor)) + constant_term = density_inner - a_factor + + def exponential_density(radius): + """ + Create a dummy exponential density + """ + return ( + a_factor * np.exp(-b_factor * (radius - bottom) / thickness) + constant_term + ) + + # Create a set of observation points located on the shell outer radius + region = (-180, 180, -90, 90) + shape = (19, 13) + coordinates = grid_coordinates(region=region, shape=shape, extra_coords=top) + # Get analytic solution + analytic = analytical_spherical_shell_exponential( + coordinates[-1], bottom, top, a_factor, b_factor, constant_term + ) + # Check if numerical results are accurate enough + npt.assert_allclose( + tesseroid_gravity(coordinates, tesseroids, exponential_density, field=field), + analytic[field], + rtol=ACCURACY_THRESHOLD, + ) From 128d35b8a11001dbfb03906123212cfaf86c4ca6 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 15:26:29 -0300 Subject: [PATCH 07/17] Fix wrong behaviour of density_minmax function --- .../forward/_tesseroid_variable_density.py | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/harmonica/forward/_tesseroid_variable_density.py b/harmonica/forward/_tesseroid_variable_density.py index 3c210da68..6b5ca3915 100644 --- a/harmonica/forward/_tesseroid_variable_density.py +++ b/harmonica/forward/_tesseroid_variable_density.py @@ -175,11 +175,24 @@ def density_minmax(density, bottom, top): """ Compute the minimum and maximum value of a bounded density """ - minimum = minimize_scalar(density, bounds=[bottom, top], method="bounded") - maximum = minimize_scalar( - lambda radius: -density(radius), bounds=[bottom, top], method="bounded" + # Calculate min and max density values at the top and bottom boundaries + density_bounds_min, density_bounds_max = np.sort([density(bottom), density(top)]) + # Estimate the minimum value of the density function withing bounds + kwargs = dict(bounds=[bottom, top], method="bounded") + minimum = np.min( + ( + minimize_scalar(density, **kwargs).fun, + density_bounds_min, + ) + ) + # Estimate the maximum value of the density function withing bounds + maximum = np.max( + ( + -minimize_scalar(lambda radius: -density(radius), **kwargs).fun, + density_bounds_max, + ) ) - return minimum.fun, -maximum.fun + return minimum, maximum def maximum_absolute_diff(normalized_density, bottom, top): From 13cc8b4a3712088b70262854ea96b51e070abe45 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 15:57:08 -0300 Subject: [PATCH 08/17] Limit the amount of parameters on exp density --- harmonica/tests/test_tesseroid_variable_density.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 03c17daa6..9aa4c122b 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -421,7 +421,7 @@ def linear_density(radius): @pytest.mark.use_numba @pytest.mark.parametrize("field", ("potential", "g_z")) @pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) -@pytest.mark.parametrize("b_factor", (1, 2, 5, 10, 30, 100)) +@pytest.mark.parametrize("b_factor", (5, 100)) def test_spherical_shell_exponential_density(field, thickness, b_factor): """ Test numerical results against analytical solution for exponential density From d6457c57cd8a3d43d1169c0ae6e07769d58a9be4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 15:57:43 -0300 Subject: [PATCH 09/17] Test a single tesseroid with constant density --- .../tests/test_tesseroid_variable_density.py | 32 +++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 9aa4c122b..450fe10e0 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -296,6 +296,38 @@ def stupid_constant_density(radius): npt.assert_allclose(tesseroids[0], tesseroid) +# ---------------------------------------------- +# Test single tesseroid against constant density +# ---------------------------------------------- + + +@pytest.mark.use_numba +@pytest.mark.parametrize("field", ("potential", "g_z")) +def test_single_tesseroid_against_constant_density(field): + """ + Test the output of a single tesseroid against one with constant density + """ + # Define a single tesseroid with constant density + bottom, top = 5400e3, 6300e3 + tesseroid = [-3, 3, -2, 2, bottom, top] + density = 2900.0 + + # Define a constant density + def constant_density(radius): + return density + + # Define a set of observation points + coordinates = grid_coordinates(region=(-5, 5, -5, 5), spacing=1, extra_coords=top) + + # Compare effects against the constant density implementation + npt.assert_allclose( + tesseroid_gravity(coordinates, [tesseroid], density=[density], field=field), + tesseroid_gravity( + coordinates, [tesseroid], density=constant_density, field=field + ), + ) + + # ---------------------------- # Test against spherical shell # ---------------------------- From fe5a52a5686a0036e13502e437ada0b442004e74 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 16:01:00 -0300 Subject: [PATCH 10/17] Remove use_numba mark where unneeded --- harmonica/tests/test_tesseroid_variable_density.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 450fe10e0..5cd846ecb 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -411,7 +411,6 @@ def build_spherical_shell(bottom, top, shape=(6, 12)): @require_numba -@pytest.mark.use_numba @pytest.mark.parametrize("field", ("potential", "g_z")) @pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) def test_spherical_shell_linear_density(field, thickness): @@ -450,7 +449,6 @@ def linear_density(radius): @require_numba -@pytest.mark.use_numba @pytest.mark.parametrize("field", ("potential", "g_z")) @pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) @pytest.mark.parametrize("b_factor", (5, 100)) From eceb404b0334f7544e5c36bae731fac7c9537544 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 16:03:41 -0300 Subject: [PATCH 11/17] Replace require_numba for run_only_with_numba --- harmonica/tests/test_tesseroid_variable_density.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 5cd846ecb..69db929b5 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -13,7 +13,7 @@ import pytest from verde import grid_coordinates -from .utils import require_numba +from .utils import run_only_with_numba from ..constants import GRAVITATIONAL_CONST from .. import tesseroid_gravity from ..forward._tesseroid_variable_density import ( @@ -410,7 +410,7 @@ def build_spherical_shell(bottom, top, shape=(6, 12)): return tesseroids -@require_numba +@run_only_with_numba @pytest.mark.parametrize("field", ("potential", "g_z")) @pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) def test_spherical_shell_linear_density(field, thickness): @@ -448,7 +448,7 @@ def linear_density(radius): ) -@require_numba +@run_only_with_numba @pytest.mark.parametrize("field", ("potential", "g_z")) @pytest.mark.parametrize("thickness", (1e2, 1e3, 1e6)) @pytest.mark.parametrize("b_factor", (5, 100)) From 0b8a691d5bd5a4f5d16aaed9ebc5077389ede5c8 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Mon, 25 Oct 2021 16:15:30 -0300 Subject: [PATCH 12/17] Fix pylint complains --- .../forward/_tesseroid_variable_density.py | 6 ++++-- .../tests/test_tesseroid_variable_density.py | 19 ++++++++++++------- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/harmonica/forward/_tesseroid_variable_density.py b/harmonica/forward/_tesseroid_variable_density.py index 6b5ca3915..9e3d51232 100644 --- a/harmonica/forward/_tesseroid_variable_density.py +++ b/harmonica/forward/_tesseroid_variable_density.py @@ -7,7 +7,7 @@ """ Utils functions for tesseroids with variable density """ -from numba import jit, prange +from numba import jit import numpy as np from scipy.optimize import minimize_scalar @@ -127,7 +127,9 @@ def density_based_discretization(tesseroids, density): return np.atleast_2d(discretized_tesseroids) -def _density_based_discretization(tesseroid, density): +def _density_based_discretization( + tesseroid, density +): # pylint: disable=too-many-locals """ Applies density-based discretization to a single tesseroid diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index 69db929b5..dbb366652 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -7,11 +7,11 @@ """ Test forward modelling for tesseroids with variable density """ +import pytest import numpy as np import numpy.testing as npt -import harmonica -import pytest from verde import grid_coordinates +import harmonica from .utils import run_only_with_numba from ..constants import GRAVITATIONAL_CONST @@ -36,18 +36,21 @@ @pytest.fixture(name="bottom") def fixture_bottom(): + """Return a bottom boundary""" bottom = 2e3 return bottom @pytest.fixture(name="top") def fixture_top(): + """Return a top boundary""" top = 5e3 return top @pytest.fixture(name="quadratic_params") def fixture_quadratic_params(): + """Return parameters for building a quadratic density function""" factor = 1e-3 vertex_radius = 3e3 vertex_density = 1900.0 @@ -240,7 +243,7 @@ def test_density_based_discret_with_delta( bottom, top, quadratic_density, -): +): # pylint: disable=protected-access """ Test the density-based discretization algorithm against values of DELTA """ @@ -287,7 +290,7 @@ def test_density_based_discret_constant_density(): w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50 tesseroid = [w, e, s, n, bottom, top] - def stupid_constant_density(radius): + def stupid_constant_density(radius): # pylint: disable=unused-argument """Define a dummy constant density function""" return 3 @@ -313,7 +316,7 @@ def test_single_tesseroid_against_constant_density(field): density = 2900.0 # Define a constant density - def constant_density(radius): + def constant_density(radius): # pylint: disable=unused-argument return density # Define a set of observation points @@ -350,7 +353,7 @@ def analytical_spherical_shell_linear(radius, bottom, top, slope, constant_term) def analytical_spherical_shell_exponential( radius, bottom, top, a_factor, b_factor, constant_term -): +): # pylint: disable=too-many-locals r""" Analytical solutions of a spherical shell with exponential density @@ -385,7 +388,9 @@ def analytical_spherical_shell_exponential( return data -def build_spherical_shell(bottom, top, shape=(6, 12)): +def build_spherical_shell( + bottom, top, shape=(6, 12) +): # pylint: disable=too-many-locals """ Return a set of tesseroids modelling a spherical shell From e20831c619eed47d7850f7b647b6facc48e451fe Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 26 Oct 2021 12:43:14 -0300 Subject: [PATCH 13/17] Add Soler2019 reference and improve docstrings --- doc/references.rst | 1 + .../forward/_tesseroid_variable_density.py | 23 +++++++++++++++++++ harmonica/forward/tesseroid.py | 19 ++++++++++++++- 3 files changed, 42 insertions(+), 1 deletion(-) diff --git a/doc/references.rst b/doc/references.rst index 9b1b9badf..21a558d8b 100644 --- a/doc/references.rst +++ b/doc/references.rst @@ -12,5 +12,6 @@ References .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer. .. [Nagy2000] Nagy, D., Papp, G. & Benedek, J.(2000). The gravitational potential and its derivatives for the prism. Journal of Geodesy 74: 552. doi:`10.1007/s001900000116 `__ .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2002). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy. doi:`10.1007/s00190-002-0264-7 `__ +.. [Soler2019] Soler, S. R., Pesce, A., Gimenez, M. E., & Uieda, L. (2019). Gravitational field calculation in spherical coordinates using variable densities in depth, Geophysical Journal International. doi: `10.1093/gji/ggz277 `__ .. [Vajda2004] Vajda, P., Vaníček, P., Novák, P. and Meurers, B. (2004). On evaluation of Newton integrals in geodetic coordinates: Exact formulation and spherical approximation. Contributions to Geophysics and Geodesy, 34(4), 289-314. .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press. diff --git a/harmonica/forward/_tesseroid_variable_density.py b/harmonica/forward/_tesseroid_variable_density.py index 9e3d51232..4b4c99180 100644 --- a/harmonica/forward/_tesseroid_variable_density.py +++ b/harmonica/forward/_tesseroid_variable_density.py @@ -120,6 +120,23 @@ def gauss_legendre_quadrature_variable_density( def density_based_discretization(tesseroids, density): """ Apply density_based discretization to a collection of tesseroids + + Parameters + ---------- + tesseroids : 2d-array + Array containing the coordinates of the tesseroid. Each row of the + array should contain the boundaries of each tesseroid in the following + order: ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + The longitudinal and latitudinal boundaries should be in degrees, while + the radial ones must be in meters. + density : func + Continuous density function of the tesseroid in SI units. + + Returns + ------- + discretized_tesseroids : 2d-array + Array containing the coordinates of radially discretized tesseriods. + Each row of the array will have the boundaries for each new tesseroid. """ discretized_tesseroids = [] for tesseroid in tesseroids: @@ -138,11 +155,17 @@ def _density_based_discretization( Parameters ---------- tesseroid : tuple + Tuple containing the boundaries of the tesseroid: + ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + Horizontal boundaries should be in degrees and radial boundaries in + meters. density : func + Density func of the tesseroid in SI units. Returns ------- tesseroids : list + List containing the boundaries of discretized tesseroids. """ # Define normalized density def normalized_density(radius): diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py index 4227697c5..b004e71d3 100644 --- a/harmonica/forward/tesseroid.py +++ b/harmonica/forward/tesseroid.py @@ -67,8 +67,11 @@ def tesseroid_gravity( spherical coordinate system. The longitudinal and latitudinal boundaries should be in degrees, while the radial ones must be in meters. - density : list or array + density : list, array or func List or array containing the density of each tesseroid in kg/m^3. + Alternatively, it can be a continuous function within the boundaries of + the tesseroids, in which case the variable density tesseroids method + introduced in [Soler2019]_ will be used. field : str Gravitational field that wants to be computed. The available fields are: @@ -102,6 +105,10 @@ def tesseroid_gravity( Gravitational field generated by the tesseroids on the computation points. + References + ---------- + [Soler2019]_ + Examples -------- @@ -123,6 +130,16 @@ def tesseroid_gravity( >>> tesseroid_gravity(coordinates, tesseroid, density, field="g_z") array(112.54539933) + >>> # Define a linear density function for the same tesseroid + >>> def linear_density(radius): + ... density_top = 2670. + ... density_bottom = 3300. + ... slope = (density_top - density_bottom) / (top - bottom) + ... return slope * (radius - bottom) + density_bottom + >>> # Compute the downward acceleration it generates + >>> tesseroid_gravity(coordinates, tesseroid, linear_density, field="g_z") + array(125.80498632) + """ kernels = {"potential": kernel_potential_spherical, "g_z": kernel_g_z_spherical} if field not in kernels: From f0feef27fd2183959a59412bec865243d46095d7 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 26 Oct 2021 15:49:51 -0300 Subject: [PATCH 14/17] Don't jit density function inside tesseroid_gravity The density function should be njit decorated by the user before passing it into `tesseroid_gravity`. --- examples/forward/tesseroid_variable_density.py | 5 +++++ harmonica/forward/tesseroid.py | 10 ++++++---- harmonica/tests/test_tesseroid_variable_density.py | 8 ++++++++ 3 files changed, 19 insertions(+), 4 deletions(-) diff --git a/examples/forward/tesseroid_variable_density.py b/examples/forward/tesseroid_variable_density.py index 46ae37763..2c44f6d90 100644 --- a/examples/forward/tesseroid_variable_density.py +++ b/examples/forward/tesseroid_variable_density.py @@ -13,6 +13,7 @@ import verde as vd import boule as bl import harmonica as hm +from numba import njit # Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to @@ -29,7 +30,11 @@ [-60, -50, -30, -20, mean_radius - 10e3, mean_radius], ) +# Define a linear density function. We should use the njit decorator so Numba +# can run the forward model efficiently. + +@njit def density(radius): """Linear density function""" top = mean_radius diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py index b004e71d3..d40d4a49c 100644 --- a/harmonica/forward/tesseroid.py +++ b/harmonica/forward/tesseroid.py @@ -130,8 +130,11 @@ def tesseroid_gravity( >>> tesseroid_gravity(coordinates, tesseroid, density, field="g_z") array(112.54539933) - >>> # Define a linear density function for the same tesseroid - >>> def linear_density(radius): + >>> # Define a linear density function for the same tesseroid. + >>> # It should be decorated with numba.njit + >>> from numba import njit + >>> @njit + ... def linear_density(radius): ... density_top = 2670. ... density_bottom = 3300. ... slope = (density_top - density_bottom) / (top - bottom) @@ -156,8 +159,7 @@ def tesseroid_gravity( _check_points_outside_tesseroids(coordinates, tesseroids) # Check if density is a function or constant values if callable(density): - # Compile the density function and run density-based discretization - density = jit(nopython=True)(density) + # Run density-based discretization on each tesseroid tesseroids = density_based_discretization(tesseroids, density) else: density = np.atleast_1d(density).ravel() diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index dbb366652..f875545a9 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -11,6 +11,7 @@ import numpy as np import numpy.testing as npt from verde import grid_coordinates +from numba import njit import harmonica from .utils import run_only_with_numba @@ -94,6 +95,7 @@ def fixture_quadratic_density(quadratic_params): """ factor, vertex_radius, vertex_density = quadratic_params + @njit def density(radius): """Quadratic density function""" return factor * (radius - vertex_radius) ** 2 + vertex_density @@ -109,6 +111,7 @@ def fixture_straight_line_analytic(bottom, top, quadratic_density): density_bottom, density_top = quadratic_density(bottom), quadratic_density(top) slope = (density_top - density_bottom) / (top - bottom) + @njit def line(radius): return slope * (radius - bottom) + density_bottom @@ -194,6 +197,7 @@ def test_density_minmax_exponential_function(bottom, top): thickness = top - bottom b_factor = 50 + @njit def exponential_density(radius): """ Create a dummy exponential density @@ -272,6 +276,7 @@ def test_density_based_discret_linear_density(): w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50 tesseroid = [w, e, s, n, bottom, top] + @njit def linear_density(radius): """Define a dummy linear density""" return 3 * radius + 2 @@ -316,6 +321,7 @@ def test_single_tesseroid_against_constant_density(field): density = 2900.0 # Define a constant density + @njit def constant_density(radius): # pylint: disable=unused-argument return density @@ -431,6 +437,7 @@ def test_spherical_shell_linear_density(field, thickness): slope = (density_outer - density_inner) / (top - bottom) constant_term = density_outer - slope * top + @njit def linear_density(radius): """ Create a dummy linear density @@ -470,6 +477,7 @@ def test_spherical_shell_exponential_density(field, thickness, b_factor): a_factor = (density_inner - density_outer) / (1 - np.exp(-b_factor)) constant_term = density_inner - a_factor + @njit def exponential_density(radius): """ Create a dummy exponential density From 0f5de6a433a66f45c178e505a77829f99d820448 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 26 Oct 2021 16:51:15 -0300 Subject: [PATCH 15/17] Ask for a njit decorated density on the docstrings --- harmonica/forward/tesseroid.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py index d40d4a49c..e97bdd911 100644 --- a/harmonica/forward/tesseroid.py +++ b/harmonica/forward/tesseroid.py @@ -67,11 +67,13 @@ def tesseroid_gravity( spherical coordinate system. The longitudinal and latitudinal boundaries should be in degrees, while the radial ones must be in meters. - density : list, array or func + density : list, array or func (decorated with :func:`numba.njit`) List or array containing the density of each tesseroid in kg/m^3. Alternatively, it can be a continuous function within the boundaries of the tesseroids, in which case the variable density tesseroids method introduced in [Soler2019]_ will be used. + If ``density`` is a function, it should be decorated with + :func:`numba.njit`. field : str Gravitational field that wants to be computed. The available fields are: From 92300237a386bcff2bc8a69024976cda088581cc Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 26 Oct 2021 17:09:22 -0300 Subject: [PATCH 16/17] Decorate density functions with jit instead of njit --- doc/conf.py | 1 + examples/forward/tesseroid_variable_density.py | 6 +++--- harmonica/forward/tesseroid.py | 8 ++++---- .../tests/test_tesseroid_variable_density.py | 16 ++++++++-------- 4 files changed, 16 insertions(+), 15 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index 2c6b136a6..a62ee8712 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -42,6 +42,7 @@ intersphinx_mapping = { "python": ("https://docs.python.org/3/", None), "numpy": ("https://docs.scipy.org/doc/numpy/", None), + "numba": ("https://numba.pydata.org/numba-doc/latest/", None), "scipy": ("https://docs.scipy.org/doc/scipy/reference", None), "pandas": ("http://pandas.pydata.org/pandas-docs/stable/", None), "xarray": ("http://xarray.pydata.org/en/stable/", None), diff --git a/examples/forward/tesseroid_variable_density.py b/examples/forward/tesseroid_variable_density.py index 2c44f6d90..ff412bfeb 100644 --- a/examples/forward/tesseroid_variable_density.py +++ b/examples/forward/tesseroid_variable_density.py @@ -13,7 +13,7 @@ import verde as vd import boule as bl import harmonica as hm -from numba import njit +from numba import jit # Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to @@ -30,11 +30,11 @@ [-60, -50, -30, -20, mean_radius - 10e3, mean_radius], ) -# Define a linear density function. We should use the njit decorator so Numba +# Define a linear density function. We should use the jit decorator so Numba # can run the forward model efficiently. -@njit +@jit def density(radius): """Linear density function""" top = mean_radius diff --git a/harmonica/forward/tesseroid.py b/harmonica/forward/tesseroid.py index e97bdd911..9e14a52c5 100644 --- a/harmonica/forward/tesseroid.py +++ b/harmonica/forward/tesseroid.py @@ -67,13 +67,13 @@ def tesseroid_gravity( spherical coordinate system. The longitudinal and latitudinal boundaries should be in degrees, while the radial ones must be in meters. - density : list, array or func (decorated with :func:`numba.njit`) + density : list, array or func decorated with :func:`numba.jit` List or array containing the density of each tesseroid in kg/m^3. Alternatively, it can be a continuous function within the boundaries of the tesseroids, in which case the variable density tesseroids method introduced in [Soler2019]_ will be used. If ``density`` is a function, it should be decorated with - :func:`numba.njit`. + :func:`numba.jit`. field : str Gravitational field that wants to be computed. The available fields are: @@ -134,8 +134,8 @@ def tesseroid_gravity( >>> # Define a linear density function for the same tesseroid. >>> # It should be decorated with numba.njit - >>> from numba import njit - >>> @njit + >>> from numba import jit + >>> @jit ... def linear_density(radius): ... density_top = 2670. ... density_bottom = 3300. diff --git a/harmonica/tests/test_tesseroid_variable_density.py b/harmonica/tests/test_tesseroid_variable_density.py index f875545a9..650911d87 100644 --- a/harmonica/tests/test_tesseroid_variable_density.py +++ b/harmonica/tests/test_tesseroid_variable_density.py @@ -11,7 +11,7 @@ import numpy as np import numpy.testing as npt from verde import grid_coordinates -from numba import njit +from numba import jit import harmonica from .utils import run_only_with_numba @@ -95,7 +95,7 @@ def fixture_quadratic_density(quadratic_params): """ factor, vertex_radius, vertex_density = quadratic_params - @njit + @jit def density(radius): """Quadratic density function""" return factor * (radius - vertex_radius) ** 2 + vertex_density @@ -111,7 +111,7 @@ def fixture_straight_line_analytic(bottom, top, quadratic_density): density_bottom, density_top = quadratic_density(bottom), quadratic_density(top) slope = (density_top - density_bottom) / (top - bottom) - @njit + @jit def line(radius): return slope * (radius - bottom) + density_bottom @@ -197,7 +197,7 @@ def test_density_minmax_exponential_function(bottom, top): thickness = top - bottom b_factor = 50 - @njit + @jit def exponential_density(radius): """ Create a dummy exponential density @@ -276,7 +276,7 @@ def test_density_based_discret_linear_density(): w, e, s, n, bottom, top = -3, 2, -4, 5, 30, 50 tesseroid = [w, e, s, n, bottom, top] - @njit + @jit def linear_density(radius): """Define a dummy linear density""" return 3 * radius + 2 @@ -321,7 +321,7 @@ def test_single_tesseroid_against_constant_density(field): density = 2900.0 # Define a constant density - @njit + @jit def constant_density(radius): # pylint: disable=unused-argument return density @@ -437,7 +437,7 @@ def test_spherical_shell_linear_density(field, thickness): slope = (density_outer - density_inner) / (top - bottom) constant_term = density_outer - slope * top - @njit + @jit def linear_density(radius): """ Create a dummy linear density @@ -477,7 +477,7 @@ def test_spherical_shell_exponential_density(field, thickness, b_factor): a_factor = (density_inner - density_outer) / (1 - np.exp(-b_factor)) constant_term = density_inner - a_factor - @njit + @jit def exponential_density(radius): """ Create a dummy exponential density From bbc8de6ae0d2d9157a9a22e03b365db9a77340d7 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 12 Nov 2021 13:21:30 -0300 Subject: [PATCH 17/17] Add introductory text to example Replace jit for njit decorator on example. --- .../forward/tesseroid_variable_density.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/examples/forward/tesseroid_variable_density.py b/examples/forward/tesseroid_variable_density.py index ff412bfeb..990deb3da 100644 --- a/examples/forward/tesseroid_variable_density.py +++ b/examples/forward/tesseroid_variable_density.py @@ -7,13 +7,28 @@ """ Tesseroids with variable density ================================ + +The :func:`harmonica.tesseroid_gravity` is capable of computing the +gravitational effects of tesseroids whose density is defined through +a continuous function of the radial coordinate. This is achieved by the +application of the method introduced in [Soler2021]_. + +To do so we need to define a regular Python function for the density, which +should have a single argument (the ``radius`` coordinate) and return the +density of the tesseroids at that radial coordinate. +In addition, we need to decorate the density function with +:func:`numba.jit(nopython=True)` or ``numba.njit`` for short. + +On this example we will show how we can compute the gravitational effect of +four tesseroids whose densities are given by a custom linear ``density`` +function. """ import matplotlib.pyplot as plt import cartopy.crs as ccrs import verde as vd import boule as bl import harmonica as hm -from numba import jit +from numba import njit # Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to @@ -34,7 +49,7 @@ # can run the forward model efficiently. -@jit +@njit def density(radius): """Linear density function""" top = mean_radius