Skip to content

Commit

Permalink
Update point source use Choclo
Browse files Browse the repository at this point in the history 
1. Update point source use Choclo
2. Rename the kernel from kernel_g_z_spherical to kernel_gravity_u_spherical
3. Update tesseroid associate with kernel name change
4. Update Test add Choclo
  • Loading branch information
@LL-Geo
LL-Geo committed Jul 7, 2023
1 parent a8b4359 commit 9bc2e04
Show file tree
Hide file tree
Showing 3 changed files with 177 additions and 217 deletions.
235 changes: 53 additions & 182 deletions harmonica/_forward/point.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,22 @@
"""

import numpy as np
from choclo.point import (
gravity_e,
gravity_ee,
gravity_en,
gravity_eu,
gravity_n,
gravity_nn,
gravity_nu,
gravity_pot,
gravity_u,
gravity_uu,
)
from numba import jit, prange

from ..constants import GRAVITATIONAL_CONST
from .utils import check_coordinate_system, distance_cartesian, distance_spherical_core
from .utils import check_coordinate_system, distance_spherical_core


def point_gravity(
Expand Down Expand Up @@ -216,9 +228,11 @@ def point_gravity(
dispatcher(coordinate_system, parallel)(
*coordinates, *points, masses, result, kernel
)
result *= GRAVITATIONAL_CONST
# Invert sign of gravity_u, gravity_eu, gravity_nu
if field in ("g_z", "g_ez", "g_nz"):
result *= -1
# Convert to more convenient units
if field in ("g_easting", "g_northing", "g_z"):
if field in ("g_e", "g_n", "g_z"):
result *= 1e5 # SI to mGal
tensors = ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz", "g_ne", "g_ze", "g_zn")
if field in tensors:
Expand Down Expand Up @@ -249,27 +263,27 @@ def get_kernel(coordinate_system, field):
"""
kernels = {
"cartesian": {
"potential": kernel_potential_cartesian,
"g_z": kernel_g_z_cartesian,
"g_northing": kernel_g_northing_cartesian,
"g_easting": kernel_g_easting_cartesian,
"potential": gravity_pot,
"g_z": gravity_u,
"g_n": gravity_n,
"g_e": gravity_e,
# diagonal tensor components
"g_ee": kernel_g_ee_cartesian,
"g_nn": kernel_g_nn_cartesian,
"g_zz": kernel_g_zz_cartesian,
"g_ee": gravity_ee,
"g_nn": gravity_nn,
"g_zz": gravity_uu,
# non-diagonal tensor components
"g_en": kernel_g_en_cartesian,
"g_ez": kernel_g_ez_cartesian,
"g_nz": kernel_g_nz_cartesian,
"g_ne": kernel_g_en_cartesian,
"g_ze": kernel_g_ez_cartesian,
"g_zn": kernel_g_nz_cartesian,
"g_en": gravity_en,
"g_ez": gravity_eu,
"g_nz": gravity_nu,
"g_ne": gravity_en,
"g_ze": gravity_eu,
"g_zn": gravity_nu,
},
"spherical": {
"potential": kernel_potential_spherical,
"g_z": kernel_g_z_spherical,
"g_northing": None,
"g_easting": None,
"g_z": kernel_gravity_u_spherical,
"g_n": None,
"g_e": None,
},
}
if field not in kernels[coordinate_system]:
Expand All @@ -281,160 +295,7 @@ def get_kernel(coordinate_system, field):


# ------------------------------------------
# Kernel functions for Cartesian coordinates
# ------------------------------------------


@jit(nopython=True)
def kernel_potential_cartesian(
easting, northing, upward, easting_p, northing_p, upward_p
):
"""
Kernel function for gravitational potential field in Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return 1 / distance


# Acceleration components
# -----------------------


@jit(nopython=True)
def kernel_g_z_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel for downward component of gravitational acceleration
Use Cartesian coords.
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
# Remember that the ``g_z`` field returns the downward component of the
# gravitational acceleration. As a consequence, it is multiplied by -1.
# Notice that the ``g_z`` does not have the minus signal observed at the
# compoents ``g_northing`` and ``g_easting``.
return (upward - upward_p) / distance**3


@jit(nopython=True)
def kernel_g_northing_cartesian(
easting, northing, upward, easting_p, northing_p, upward_p
):
"""
Kernel function for northing component of gravitational acceleration
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return -(northing - northing_p) / distance**3


@jit(nopython=True)
def kernel_g_easting_cartesian(
easting, northing, upward, easting_p, northing_p, upward_p
):
"""
Kernel function for easting component of gravitational acceleration
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return -(easting - easting_p) / distance**3


# Tensor components
# -----------------


@jit(nopython=True)
def kernel_g_ee_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel function for g_ee component of gravitational tensor
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return 3 * (easting - easting_p) ** 2 / distance**5 - 1 / distance**3


@jit(nopython=True)
def kernel_g_nn_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel function for g_nn component of gravitational tensor
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return 3 * (northing - northing_p) ** 2 / distance**5 - 1 / distance**3


@jit(nopython=True)
def kernel_g_zz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel function for g_zz component of gravitational tensor
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return 3 * (upward - upward_p) ** 2 / distance**5 - 1 / distance**3


@jit(nopython=True)
def kernel_g_en_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel function for g_en component of gravitational tensor
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
return 3 * (easting - easting_p) * (northing - northing_p) / distance**5


@jit(nopython=True)
def kernel_g_ez_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel function for g_ez component of gravitational tensor
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
# Add a minus sign to account that the z axis points downwards.
return -3 * (easting - easting_p) * (upward - upward_p) / distance**5


@jit(nopython=True)
def kernel_g_nz_cartesian(easting, northing, upward, easting_p, northing_p, upward_p):
"""
Kernel function for g_nz component of gravitational tensor
Use Cartesian coordinates
"""
distance = distance_cartesian(
(easting, northing, upward), (easting_p, northing_p, upward_p)
)
# Add a minus sign to account that the z axis points downwards.
return -3 * (northing - northing_p) * (upward - upward_p) / distance**5


# ------------------------------------------
# Kernel functions for Cartesian coordinates
# Kernel functions for Spherical coordinates
# ------------------------------------------


Expand All @@ -448,31 +309,39 @@ def kernel_potential_spherical(
distance, _, _ = distance_spherical_core(
longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
)
return 1 / distance
return 1 / distance * GRAVITATIONAL_CONST


# Acceleration components
# -------------------


@jit(nopython=True)
def kernel_g_z_spherical(
def kernel_gravity_u_spherical(
longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
):
"""
Kernel for downward component of gravitational acceleration
Kernel for upward component of gravitational acceleration
Use spherical coordinates
"""
distance, cospsi, _ = distance_spherical_core(
longitude, cosphi, sinphi, radius, longitude_p, cosphi_p, sinphi_p, radius_p
)
delta_z = radius - radius_p * cospsi
return delta_z / distance**3
return -GRAVITATIONAL_CONST * delta_z / distance**3


def point_mass_cartesian(
easting, northing, upward, easting_p, northing_p, upward_p, masses, out, kernel
easting,
northing,
upward,
easting_p,
northing_p,
upward_p,
masses,
out,
forward_func,
):
"""
Compute gravitational field of point masses in Cartesian coordinates
Expand All @@ -489,19 +358,21 @@ def point_mass_cartesian(
Array where the gravitational field on each computation point will be
appended.
It must have the same size of ``easting``, ``northing`` and ``upward``.
kernel : func
Kernel function that will be used to compute the gravitational field on
the computation points.
forward_func : func
forward_func function that will be used to compute the gravitational
field on the computation points. It could be one of the forward
modelling functions in :mod:`choclo.point`.
"""
for l in prange(easting.size):
for m in range(easting_p.size):
out[l] += masses[m] * kernel(
out[l] += forward_func(
easting[l],
northing[l],
upward[l],
easting_p[m],
northing_p[m],
upward_p[m],
masses[m],
)


Expand Down
Loading

0 comments on commit 9bc2e04

Please sign in to comment.