Forward models of prisms gravity fields with Choclo (#400)

Replace the prisms gravity kernels for Choclo functions. Extend the
capabilities of the `prism_gravity` function allowing to compute all the
gravitational acceleration and tensor components. Add tests that compare
Harmonica results against dumb Choclo calls and a Laplacian equation
test for the diagonal components of the tensor. Ditch symmetry tests
that were already covered by Choclo. Improve docstrings. Raise warning
if any computation point falls in a singular point for the tensor
components. Add tests for these new checks and warnings.

---------

Co-authored-by: Lu Li <[email protected]>

doc/user_guide/forward_modelling/prism.rst

@@ -48,6 +48,97 @@ computation point:
  print(potential, "J/kg")
 
 
+Gravitational fields
+^^^^^^^^^^^^^^^^^^^^
+
+The :func:`harmonica.prism_gravity` is able to compute the gravitational
+potential (``"potential"``), the acceleration components (``"g_e"``, ``"g_n"``,
+``"g_z"``), and tensor components (``"g_ee"``, ``"g_nn"``, ``"g_zz"``,
+``"g_en"``, ``"g_ez"``, ``"g_nz"``).
+
+
+Build a regular grid of computation points located a 10m above the zero height:
+
+.. jupyter-execute::
+
+ import verde as vd
+
+ region = (-10e3, 10e3, -10e3, 10e3)
+ shape = (51, 51)
+ height = 10
+ coordinates = vd.grid_coordinates(region, shape=shape, extra_coords=height)
+
+Define a single prism:
+
+.. jupyter-execute::
+
+ prism = [-2e3, 2e3, -2e3, 2e3, -1.6e3, -900]
+ density = 3300
+
+Compute the gravitational fields that this prism generate on each observation
+point:
+
+.. jupyter-execute::
+
+ fields = (
+ "potential",
+ "g_e", "g_n", "g_z",
+ "g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"
+ )
+
+ results = {}
+ for field in fields:
+ results[field] = hm.prism_gravity(coordinates, prism, density, field=field)
+
+Plot the results:
+
+.. jupyter-execute::
+
+ import matplotlib.pyplot as plt
+
+ plt.pcolormesh(coordinates[0], coordinates[1], results["potential"])
+ plt.gca().set_aspect("equal")
+ plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
+ plt.colorbar(label="J/kg")
+ plt.show()
+
+
+.. jupyter-execute::
+
+ fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+ for field, ax in zip(("g_e", "g_n", "g_z"), axes):
+ tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+ ax.set_aspect("equal")
+ ax.set_title(field)
+ ax.ticklabel_format(style="sci", scilimits=(0, 0))
+ plt.colorbar(tmp, ax=ax, label="mGal", orientation="horizontal", pad=0.08)
+ plt.show()
+
+.. jupyter-execute::
+
+ fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+ for field, ax in zip(("g_ee", "g_nn", "g_zz"), axes):
+ tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+ ax.set_aspect("equal")
+ ax.set_title(field)
+ ax.ticklabel_format(style="sci", scilimits=(0, 0))
+ plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
+ plt.show()
+
+.. jupyter-execute::
+
+ fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+ for field, ax in zip(("g_en", "g_ez", "g_nz"), axes):
+ tmp = ax.pcolormesh(coordinates[0], coordinates[1], results[field])
+ ax.set_aspect("equal")
+ ax.set_title(field)
+ ax.ticklabel_format(style="sci", scilimits=(0, 0))
+ plt.colorbar(tmp, ax=ax, label="Eotvos", orientation="horizontal", pad=0.08)
+ plt.show()
+
 Passing multiple prisms
 ^^^^^^^^^^^^^^^^^^^^^^^
 

environment.yml

@@ -18,6 +18,7 @@ dependencies:
  - verde>=1.7.0
  - xarray>=0.16
  - xrft>=1.0
+ - choclo>=0.1
  # Optional requirements
  - pyvista>=0.27
  - vtk>=9