- [TEST] Move geophysics to plugins. NOTE: Probably is not a good idea but we will see once I get into

Author: Leguark

gempy_plugins/assets/geophysics.py

@@ -18,7 +18,7 @@
 import numpy as np
 import aesara
 import aesara.tensor  as T
-from gempy.core.grid_modules.grid_types import CenteredGrid
+from gempy.core.data.grid_types import CenteredGrid
 
 
 class GravityPreprocessing(CenteredGrid):

tests/_geophysics_TO_UPDATE/__init__.py

@@ -0,0 +1,3 @@
+import pytest
+
+pytestmark = pytest.mark.skip(reason="Geophysics module is not yet updated")

tests/_geophysics_TO_UPDATE/test_gravity.py

@@ -0,0 +1,35 @@
+# These two lines are necessary only if GemPy is not installed
+# sys.path.append("../..")
+
+# Importing GemPy
+import gempy as gp
+
+# Importing auxiliary libraries
+import numpy as np
+
+
+def test_gravity(interpolator_gravity):
+    geo_model = gp.create_model('2-layers')
+    gp.init_data(geo_model, extent=[0, 12, -2, 2, 0, 4], resolution=[500, 1, 500])
+    geo_model.add_surfaces('surface 1')
+    geo_model.add_surfaces('surface 2')
+    geo_model.add_surfaces('basement')
+    dz = geo_model._grid.regular_grid.dz
+    geo_model._surfaces.add_surfaces_values([dz, 0, 0], ['dz'])
+    geo_model._surfaces.add_surfaces_values([2.6, 2.4, 3.2], ['density'])
+    geo_model.add_surface_points(3, 0, 3.05, 'surface 1')
+    geo_model.add_surface_points(9, 0, 3.05, 'surface 1')
+
+    geo_model.add_surface_points(3, 0, 1.02, 'surface 2')
+    geo_model.add_surface_points(9, 0, 1.02, 'surface 2')
+
+    geo_model.add_orientations(6, 0, 4, 'surface 1', [0, 0, 1])
+    device_loc = np.array([[6, 0, 4]])
+
+    geo_model.set_centered_grid(device_loc, resolution=[10, 10, 100], radius=16000)
+    geo_model.set_aesara_function(interpolator_gravity)
+    geo_model._interpolator.set_aesara_shared_gravity(pos_density=2)
+    print(geo_model._additional_data)
+    gp.compute_model(geo_model, set_solutions=True, compute_mesh=False)
+    print(geo_model.solutions.fw_gravity)
+    np.testing.assert_almost_equal(geo_model.solutions.fw_gravity, np.array([-1624.1714]), decimal=4)

tests/_geophysics_TO_UPDATE/test_grid_tz.py

@@ -0,0 +1,15 @@
+# These two lines are necessary only if GemPy is not installed
+import sys
+
+sys.path.append("../../../gempy")
+
+# Importing GemPy
+from gempy.core.data.grid_modules import CenteredGrid
+
+# Importing auxiliary libraries
+import numpy as np
+
+
+def test_irregular_grid():
+    g = CenteredGrid()
+    g.set_centered_grid(np.array([0, 0, 0]), resolution=[5, 5, 5], radius = [100,100,100])

tests/_geophysics_TO_UPDATE/test_magnetics.py

@@ -0,0 +1,85 @@
+# Importing GemPy
+import gempy as gp
+from gempy_plugins.assets.geophysics import MagneticsPreprocessing
+
+# Importing auxiliary libraries
+import numpy as np
+import pytest
+
+
+def test_magnetics_api(interpolator_magnetics):
+    # TODO add the check
+
+    geo_model = gp.create_model('test_center_grid_slicing')
+    geo_model.set_default_surfaces()
+    geo_model.add_surface_points(X=-1, Y=0, Z=0, surface='surface1')
+    geo_model.add_surface_points(X=1, Y=0, Z=0, surface='surface1')
+    geo_model.add_orientations(X=0, Y=0, Z=0, surface='surface1', pole_vector=(0, 0, 1))
+    geo_model._surfaces.add_surfaces_values([0.037289, 0.0297], ['susceptibility'])
+
+    # needed constants
+    mu_0 = 4.0 * np.pi * 10e-7  # magnetic permeability in free space [N/A^2]
+    cm = mu_0 / (4.0 * np.pi)  # constant for SI unit
+    incl = 77.0653  # NOAA
+    decl = 6.8116
+    B_ext = 52819.8506939139e-9  # T
+
+    geo_model.set_regular_grid(extent=[-5, 5, -5, 5, -5, 5], resolution=[5, 5, 5])
+    geo_model.set_centered_grid(np.array([[0, 0, 0]]), resolution=[10, 10, 15], radius=5000)
+    geo_model.set_aesara_function(interpolator_magnetics)
+    geo_model._interpolator.set_aesara_shared_magnetics(V='auto', pos_magnetics=1,
+                                                        incl= incl, decl=decl, B_ext=B_ext)
+
+    gp.compute_model(geo_model)
+    print(geo_model._interpolator.aesara_graph.lg0.get_value())
+    print(geo_model.solutions.fw_magnetics)
+    np.testing.assert_almost_equal(geo_model.solutions.fw_magnetics,
+                                   np.array([473.7836]), decimal=4)
+
+
+@pytest.fixture(scope="module")
+def test_magnetics_no_regular_grid(interpolator_magnetics):
+    # TODO add the check
+
+    geo_model = gp.create_model('test_center_grid_slicing')
+    geo_model.set_default_surfaces()
+    geo_model.add_surface_points(X=-1, Y=0, Z=0, surface='surface1')
+    geo_model.add_surface_points(X=1, Y=0, Z=0, surface='surface1')
+    geo_model.add_orientations(X=0, Y=0, Z=0, surface='surface1', pole_vector=(0, 0, 1))
+    geo_model._surfaces.add_surfaces_values([0.037289, 0.0297], ['susceptibility'])
+
+    # needed constants
+    mu_0 = 4.0 * np.pi * 10e-7  # magnetic permeability in free space [N/A^2]
+    cm = mu_0 / (4.0 * np.pi)  # constant for SI unit
+    incl = 77.0653  # NOAA
+    decl = 6.8116
+    B_ext = 52819.8506939139e-9  # T
+
+    geo_model.set_centered_grid(np.array([0, 0, 0]), resolution=[10, 10, 15], radius=5000)
+
+    Vmodel = MagneticsPreprocessing(geo_model._grid.centered_grid).set_Vs_kernel()
+    # gp.set_interpolator(geo_model, output=['magnetics'])
+    geo_model.set_aesara_function(interpolator_magnetics)
+    geo_model._interpolator.set_aesara_shared_magnetics(V= Vmodel, pos_magnetics=1,
+                                                        incl= incl, decl=decl, B_ext=B_ext)
+
+    # geo_model.interpolator.aesara_graph.V.set_value(Vmodel)
+    # geo_model.interpolator.aesara_graph.incl.set_value(incl)
+    # geo_model.interpolator.aesara_graph.decl.set_value(decl)
+    # geo_model.interpolator.aesara_graph.B_ext.set_value(B_ext)
+
+    gp.compute_model(geo_model)
+    np.testing.assert_almost_equal(geo_model.solutions.fw_magnetics,
+                                   np.array([473.7836]), decimal=4)
+    
+    return geo_model
+
+
+def test_center_grid_slicing(test_magnetics_no_regular_grid):
+    geo_model = test_magnetics_no_regular_grid
+
+    geo_model.set_centered_grid(np.array([[0, 0, 0],
+                                [1, 1, 1]]), resolution=[10, 10, 15], radius=5000)
+
+    gp.compute_model(geo_model)
+    print(geo_model._interpolator.aesara_graph.lg0.get_value())