Merge pull request #340 from MiraGeoscience/GEOPY-2519

GEOPY-2519: Synthetic topography for runtests is overly sampled

Author: domfournier

simpeg_drivers/utils/synthetics/driver.py

@@ -50,7 +50,7 @@ def topography(self):
             assert self.options is not None
             entity = get_topography_surface(
                 geoh5=self.geoh5,
-                options=self.options.survey,
+                options=self.options,
             )
         self._topography = entity
         return self._topography

simpeg_drivers/utils/synthetics/meshes/octrees.py

@@ -75,7 +75,7 @@ def get_octree_mesh(
     :return mesh: The discretize TreeMesh object for computations.
     """
 
-    mesh = get_base_octree(survey, topography, cell_size, (0, 0, 2), padding_distance)
+    mesh = get_base_octree(survey, topography, cell_size, (0, 1), padding_distance)
 
     mesh = OctreeDriver.refine_tree_from_points(
         mesh, survey.vertices, levels=refinement, finalize=False

simpeg_drivers/utils/synthetics/surveys/factory.py

@@ -26,21 +26,21 @@
 
 def grid_layout(
     limits: list[float],
-    station_spacing: int,
-    line_spacing: int,
+    n_stations: int,
+    n_lines: int,
     topography: Callable,
 ):
     """
     Generates grid locations based on limits and spacing.
 
     :param limits: Tuple of (xmin, xmax, ymin, ymax).
-    :param station_spacing: Number of stations along each line.
-    :param line_spacing: Number of lines in the grid.
+    :param n_stations: Number of stations along each line.
+    :param n_lines: Number of lines in the grid.
     :param topography: Callable that generates the topography (z values).
     """
 
-    x = np.linspace(limits[0], limits[1], station_spacing)
-    y = np.linspace(limits[2], limits[3], line_spacing)
+    x = np.linspace(limits[0], limits[1], n_stations)
+    y = np.linspace(limits[2], limits[3], n_lines)
     X, Y = np.meshgrid(x, y)
     Z = topography(X, Y)
 
@@ -62,8 +62,8 @@ def get_survey(
 
     X, Y, Z = grid_layout(
         limits=options.limits,
-        station_spacing=options.n_stations,
-        line_spacing=options.n_lines,
+        n_stations=options.n_stations,
+        n_lines=options.n_lines,
         topography=options.topography,
     )
     Z += options.drape

simpeg_drivers/utils/synthetics/topography.py

@@ -13,24 +13,36 @@
 from geoh5py.objects import DrapeModel, Octree, Surface
 from scipy.spatial import Delaunay
 
-from simpeg_drivers.utils.synthetics.options import SurveyOptions
+from simpeg_drivers.utils.synthetics.options import (
+    MeshOptions,
+    SurveyOptions,
+    SyntheticsComponentsOptions,
+)
 from simpeg_drivers.utils.synthetics.surveys.factory import grid_layout
 from simpeg_drivers.utils.utils import active_from_xyz
 
 
-def get_topography_surface(geoh5: Workspace, options: SurveyOptions) -> Surface:
+def get_topography_surface(
+    geoh5: Workspace, options: SyntheticsComponentsOptions
+) -> Surface:
     """
-    Returns a topography surface with 2x the limits of the survey.
+    Returns topography with same limits as the mesh and 2x resolution of the survey.
 
     :param geoh5: Geoh5 workspace.
-    :param options: Survey options. Extents will be 2x the survey extents.
+    :param options: Options containing survey and mesh specifications.
     """
 
+    width, height = compute_mesh_extents(options.survey, options.mesh)
     X, Y, Z = grid_layout(
-        limits=[4 * k for k in options.limits],  # type: ignore
-        station_spacing=int(np.ceil((options.limits[1] - options.limits[0]) / 16)),
-        line_spacing=int(np.ceil((options.limits[3] - options.limits[2]) / 16)),
-        topography=options.topography,
+        limits=[
+            options.survey.center[0] - width / 2,
+            options.survey.center[0] + width / 2,
+            options.survey.center[1] - height / 2,
+            options.survey.center[1] + height / 2,
+        ],
+        n_stations=int(2 * width / (options.survey.width / options.survey.n_stations)),
+        n_lines=int(2 * height / (options.survey.height / options.survey.n_lines)),
+        topography=options.survey.topography,
     )
 
     vertices = np.column_stack(
@@ -45,6 +57,32 @@ def get_topography_surface(geoh5: Workspace, options: SurveyOptions) -> Surface:
     )
 
 
+def compute_mesh_extents(
+    survey_options: SurveyOptions, mesh_options: MeshOptions
+) -> tuple[float, float]:
+    """
+    Estimates the extent of the mesh from survey and mesh options.
+
+    :param survey_options: Survey options.
+    :param mesh_options: Mesh options.
+
+    :return: mesh width including padding.
+    :return: mesh height including padding.
+    """
+    width = survey_options.width
+    height = survey_options.height
+    cell_size = mesh_options.cell_size
+    padding = mesh_options.padding_distance
+
+    def next_pow2_cells(span, cell_size, padding):
+        return 2 ** np.ceil(np.log2((span + 2 * padding) / cell_size))
+
+    return (
+        cell_size[0] * next_pow2_cells(width, cell_size[0], padding),
+        cell_size[1] * next_pow2_cells(height, cell_size[1], padding),
+    )
+
+
 def get_active(
     mesh: Octree | DrapeModel, topography: Surface, name: str = "active_cells"
 ):

tests/plate_simulation/runtest/sweep_test.py

@@ -16,17 +16,16 @@
 from simpeg_drivers import assets_path
 from simpeg_drivers.plate_simulation.options import PlateSimulationOptions
 from simpeg_drivers.plate_simulation.sweep.driver import PlateSweepDriver
-from simpeg_drivers.plate_simulation.sweep.options import SweepOptions
 from simpeg_drivers.potential_fields.gravity.options import GravityForwardOptions
-from simpeg_drivers.utils.synthetics.options import SurveyOptions
+from simpeg_drivers.utils.synthetics.options import SyntheticsComponentsOptions
 from simpeg_drivers.utils.synthetics.surveys.factory import get_survey
 from simpeg_drivers.utils.synthetics.topography import get_topography_surface
 
 
 def setup_plate_sweep(workspace) -> SimPEGGroup:
-    survey_options = SurveyOptions()
-    data = get_survey(workspace, method="gravity", options=survey_options)
-    topo = get_topography_surface(workspace, survey_options)
+    options = SyntheticsComponentsOptions()
+    data = get_survey(workspace, method="gravity", options=options.survey)
+    topo = get_topography_surface(workspace, options)
 
     gravity = SimPEGGroup.create(workspace, name="gravity fwd")
     options = GravityForwardOptions.model_construct()

tests/run_tests/driver_2d_rotated_gradients_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.688327405060795, "phi_d": 190, "phi_m": 70.5}
+target_run = {"data_norm": 0.644727102942656, "phi_d": 182, "phi_m": 66.3}
 
 
 def test_dc2d_rotated_grad_fwr_run(

tests/run_tests/driver_airborne_fem_1d_test.py

@@ -41,7 +41,7 @@
 # To test the full run and validate the inversion.
 # Move this file out of the test directory and run.
 
-target_run = {"data_norm": 654.2956399855235, "phi_d": 58300, "phi_m": 396}
+target_run = {"data_norm": 652.4097142104288, "phi_d": 52800, "phi_m": 341}
 
 
 def test_fem_fwr_1d_run(

tests/run_tests/driver_airborne_tem_1d_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": 5.08855682e-10, "phi_d": 20.2, "phi_m": 169000}
+target_run = {"data_norm": 5.043446361560613e-10, "phi_d": 23, "phi_m": 194000}
 
 
 def test_airborne_tem_1d_fwr_run(

tests/run_tests/driver_airborne_tem_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": 4.0476883635425e-08, "phi_d": 22700000, "phi_m": 1030}
+target_run = {"data_norm": 4.047707815081421e-08, "phi_d": 22100000, "phi_m": 2920}
 
 
 def test_bad_waveform(tmp_path: Path):

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": 0.6441011110410102, "phi_d": 644, "phi_m": 107}
+target_run = {"data_norm": 0.606697124141973, "phi_d": 693, "phi_m": 88.5}
 
 
 def test_dc_2d_fwr_run(