GEOPY-2686: Improve topography augmentation for DC and MT

simpeg_drivers/components/topography.py

@@ -28,6 +28,7 @@
 from geoh5py.data import NumericData
 from geoh5py.objects.surveys.electromagnetics.base import LargeLoopGroundEMSurvey
 from geoh5py.shared import Entity
+from scipy.spatial import cKDTree
 
 from simpeg_drivers.components.data import InversionData
 from simpeg_drivers.components.locations import InversionLocations
@@ -37,7 +38,6 @@
     active_from_xyz,
     floating_active,
     get_containing_cells,
-    get_neighbouring_cells,
     mask_vertices_and_cells,
 )
 
@@ -121,7 +121,7 @@ def active_cells(self, mesh: InversionMesh, data: InversionData) -> np.ndarray:
             active_cells = active_from_xyz(
                 mesh.entity,
                 vertices,
-                grid_reference="bottom" if forced_to_surface else "center",
+                grid_reference="center",
                 triangulation=cells,
             )
 
@@ -181,22 +181,14 @@ def expand_actives(
         containing_cells = get_containing_cells(mesh.mesh, data)
         active_cells[containing_cells] = True
 
-        # Apply extra active cells to ensure connectivity for tree meshes
-        if isinstance(mesh.mesh, TreeMesh):
-            neighbours = get_neighbouring_cells(mesh.mesh, containing_cells)
-            neighbours_xy = np.r_[neighbours[0] + neighbours[1]]
-
-            neighbours_xy = neighbours_xy[neighbours_xy != -1]
-            # Make sure the new actives are connected to the old actives
-            new_actives = ~active_cells[neighbours_xy]
-            if np.any(new_actives):
-                neighbours = get_neighbouring_cells(
-                    mesh.mesh, neighbours_xy[new_actives]
-                )
-                neighbours_z = np.r_[neighbours[2][0]]
-                neighbours_z = neighbours_z[neighbours_z != -1]
-                active_cells[neighbours_z] = True  # z-axis neighbours
-
-            active_cells[neighbours_xy] = True  # xy-axis neighbours
+        # Apply extra active cells to ensure connectivity for neighbours
+        tree = cKDTree(mesh.mesh.cell_centers[containing_cells])
+        rad, ind = tree.query(mesh.mesh.cell_centers)
+        neighbours_xy = rad < (3 * mesh.mesh.h[0].min())
+        neighbours_xy &= (
+            mesh.mesh.cell_centers[containing_cells, :][ind, -1]
+            >= mesh.mesh.cell_centers[:, -1]
+        )
+        active_cells[neighbours_xy] = True  # xy-axis neighbours
 
         return active_cells

simpeg_drivers/utils/utils.py

@@ -505,7 +505,11 @@ def get_containing_cells(
             inds = np.r_[inds, np.hstack(line_ind)]
 
     elif isinstance(mesh, TensorMesh):
-        locations = data.drape_locations(np.unique(data.locations, axis=0))
+        potentials = data.entity.vertices
+        currents = data.entity.current_electrodes.vertices
+        locations = np.unique(np.r_[potentials, currents], axis=0)
+
+        locations = data.drape_locations(np.unique(locations, axis=0))
         xi = np.searchsorted(mesh.nodes_x, locations[:, 0]) - 1
         yi = np.searchsorted(mesh.nodes_y, locations[:, -1]) - 1
         inds = xi + yi * mesh.shape_cells[0]

tests/run_tests/driver_dc_2d_test.py

@@ -43,7 +43,7 @@
 # To test the full run and validate the inversion.
 # Move this file out of the test directory and run.
 
-target_run = {"data_norm": 11.136935742296085, "phi_d": 5570, "phi_m": 314}
+target_run = {"data_norm": 10.945383968572745, "phi_d": 5210, "phi_m": 340}
 
 
 def test_dc_2d_fwr_run(

tests/run_tests/driver_dc_test.py

@@ -38,7 +38,7 @@
 # To test the full run and validate the inversion.
 # Move this file out of the test directory and run.
 
-target_run = {"data_norm": 0.15320935486917722, "phi_d": 25.7, "phi_m": 3580}
+target_run = {"data_norm": 0.14272756694409652, "phi_d": 12.6, "phi_m": 3580}
 
 
 def test_dc_3d_fwr_run(

tests/run_tests/driver_ip_2d_test.py

@@ -39,7 +39,7 @@
 # To test the full run and validate the inversion.
 # Move this file out of the test directory and run.
 
-target_run = {"data_norm": 0.1244717397585979, "phi_d": 15500, "phi_m": 0.0002845}
+target_run = {"data_norm": 0.12345083692791135, "phi_d": 15700, "phi_m": 0.000206}
 
 
 def test_ip_2d_fwr_run(

tests/run_tests/driver_ip_test.py

@@ -37,7 +37,7 @@
 # To test the full run and validate the inversion.
 # Move this file out of the test directory and run.
 
-target_run = {"data_norm": 0.005219568872697125, "phi_d": 680, "phi_m": 1.01e-6}
+target_run = {"data_norm": 0.004747013060439004, "phi_d": 562, "phi_m": 1.43e-06}
 
 
 def test_ip_3d_fwr_run(

tests/run_tests/driver_joint_cross_gradient_test.py

@@ -325,13 +325,13 @@ def test_joint_cross_gradient_inv_run(
     # the scaling from its total misfit.
     np.testing.assert_allclose(
         driver.directives.scale_misfits.scalings,
-        [1, 0.7558, 0.5, 0.5, 0.6710],
+        [1.0, 0.755269, 0.5, 0.5, 0.675174],
         atol=1e-3,
     )
     # Check that scaling * chi factor is reflected in data misfit multipliers
     np.testing.assert_allclose(
         driver.data_misfit.multipliers,
-        [0.8, 0.6046, 0.5, 0.5, 0.6710],
+        [0.8, 0.604215, 0.5, 0.5, 0.675174],
         atol=1e-3,
     )
 

tests/run_tests/driver_mt_test.py

@@ -42,7 +42,7 @@
 # To test the full run and validate the inversion.
 # Move this file out of the test directory and run.
 
-target_run = {"data_norm": 0.025087759238073448, "phi_d": 0.385, "phi_m": 3}
+target_run = {"data_norm": 0.003584600661140228, "phi_d": 4.45, "phi_m": 5.56}
 
 
 def setup_data(workspace, survey):