diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 70008ab5..d5bf96c9 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -91,7 +91,7 @@ repos: - id: mixed-line-ending exclude: ^\.idea/.*\.xml$ - id: name-tests-test - exclude: testing_utils.py + exclude: targets.py - id: pretty-format-json args: - --autofix diff --git a/simpeg_drivers/utils/synthetics/__init__.py b/simpeg_drivers/utils/synthetics/__init__.py new file mode 100644 index 00000000..4d06f672 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/utils/synthetics/driver.py b/simpeg_drivers/utils/synthetics/driver.py new file mode 100644 index 00000000..abd7f714 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/driver.py @@ -0,0 +1,109 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.data import FloatData +from geoh5py.objects import DrapeModel, ObjectBase, Octree, Surface + +from simpeg_drivers.utils.synthetics.meshes.factory import get_mesh +from simpeg_drivers.utils.synthetics.models import get_model +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_active, + get_topography_surface, +) + + +class SyntheticsComponents: + """Creates a workspace populated with objects for simulation and subsequent inversion.""" + + def __init__( + self, + geoh5: Workspace, + options: SyntheticsComponentsOptions | None = None, + ): + if options is None: + options = SyntheticsComponentsOptions() + + self.geoh5 = geoh5 + self.options = options + self._topography: Surface | None = None + self._survey: ObjectBase | None = None + self._mesh: Octree | DrapeModel | None = None + self._active: FloatData | None = None + self._model: FloatData | None = None + + @property + def topography(self): + entity = self.geoh5.get_entity("topography")[0] + assert isinstance(entity, Surface | type(None)) + if entity is None: + assert self.options is not None + entity = get_topography_surface( + geoh5=self.geoh5, + options=self.options.survey, + ) + self._topography = entity + return self._topography + + @property + def survey(self): + entity = self.geoh5.get_entity(self.options.survey.name)[0] + assert isinstance(entity, ObjectBase | type(None)) + if entity is None: + assert self.options is not None + entity = get_survey( + geoh5=self.geoh5, + method=self.options.method, + options=self.options.survey, + ) + self._survey = entity + return self._survey + + @property + def mesh(self): + entity = self.geoh5.get_entity(self.options.mesh.name)[0] + assert isinstance(entity, Octree | DrapeModel | type(None)) + if entity is None: + assert self.options is not None + entity = get_mesh( + self.options.method, + survey=self.survey, + topography=self.topography, + options=self.options.mesh, + ) + self._mesh = entity + return self._mesh + + @property + def active(self): + entity = self.geoh5.get_entity(self.options.active.name)[0] + assert isinstance(entity, FloatData | type(None)) + if entity is None: + entity = get_active(self.mesh, self.topography) + self._active = entity + return self._active + + @property + def model(self): + entity = self.geoh5.get_entity(self.options.model.name)[0] + assert isinstance(entity, FloatData | type(None)) + if entity is None: + assert self.options is not None + entity = get_model( + method=self.options.method, + mesh=self.mesh, + active=self.active.values, + options=self.options.model, + ) + self._model = entity + return self._model diff --git a/simpeg_drivers/utils/synthetics/meshes/__init__.py b/simpeg_drivers/utils/synthetics/meshes/__init__.py new file mode 100644 index 00000000..4d06f672 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/meshes/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/utils/synthetics/meshes/factory.py b/simpeg_drivers/utils/synthetics/meshes/factory.py new file mode 100644 index 00000000..19fc12a5 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/meshes/factory.py @@ -0,0 +1,45 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from discretize import TensorMesh, TreeMesh +from geoh5py.objects import CellObject, DrapeModel, Octree, Points, Surface + +from simpeg_drivers.utils.synthetics.options import MeshOptions + +from .octrees import get_octree_mesh +from .tensors import get_tensor_mesh + + +def get_mesh( + method: str, + survey: Points, + topography: Surface, + options: MeshOptions, +) -> DrapeModel | Octree: + """Factory for mesh creation with behaviour modified by the provided method.""" + + if "2d" in method: + assert isinstance(survey, CellObject) + return get_tensor_mesh( + survey=survey, + cell_size=options.cell_size, + padding_distance=options.padding_distance, + name=options.name, + ) + + return get_octree_mesh( + survey=survey, + topography=topography, + cell_size=options.cell_size, + refinement=options.refinement, + padding_distance=options.padding_distance, + refine_on_receivers=method in ["fdem", "airborne tdem"], + name=options.name, + ) diff --git a/simpeg_drivers/utils/synthetics/meshes/octrees.py b/simpeg_drivers/utils/synthetics/meshes/octrees.py new file mode 100644 index 00000000..7c3c9043 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/meshes/octrees.py @@ -0,0 +1,90 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from discretize import TreeMesh +from discretize.utils import mesh_builder_xyz +from geoh5py.objects import Octree, Points, Surface +from octree_creation_app.driver import OctreeDriver +from octree_creation_app.utils import treemesh_2_octree + + +def get_base_octree( + survey: Points, + topography: Surface, + cell_size: tuple[float, float, float], + refinement: tuple, + padding: float, +) -> TreeMesh: + """ + Generate a survey centered TreeMesh object with topography refinement. + + :param survey: Survey object with vertices that define the core of the + tensor mesh. + :param topography: Surface used to refine the topography. + :param cell_size: Tuple defining the cell size in all directions. + :param refinement: Tuple containing the number of cells to refine at each + level around the topography. + :param padding: Distance to pad the mesh in all directions. + + :return mesh: The discretize TreeMesh object for computations. + """ + padding_distance = np.ones((3, 2)) * padding + mesh = mesh_builder_xyz( + survey.vertices - np.r_[cell_size] / 2.0, + cell_size, + depth_core=100.0, + padding_distance=padding_distance, + mesh_type="TREE", + tree_diagonal_balance=False, + ) + mesh = OctreeDriver.refine_tree_from_surface( + mesh, topography, levels=refinement, finalize=False + ) + + return mesh + + +def get_octree_mesh( + survey: Points, + topography: Surface, + cell_size: tuple[float, float, float], + refinement: tuple, + padding_distance: float, + refine_on_receivers: bool, + name: str = "mesh", +) -> Octree: + """Generate a survey centered mesh with topography and survey refinement. + + :param survey: Survey object with vertices that define the core of the + tensor mesh and the source refinement for EM methods. + :param topography: Surface used to refine the topography. + :param cell_size: Tuple defining the cell size in all directions. + :param refinement: Tuple containing the number of cells to refine at each + level around the topography. + :param padding: Distance to pad the mesh in all directions. + :param refine_on_receivers: Refine on the survey locations or not. + + :return entity: The geoh5py Octree object to store the results of + computation in the shared cells of the computational mesh. + :return mesh: The discretize TreeMesh object for computations. + """ + + mesh = get_base_octree(survey, topography, cell_size, refinement, padding_distance) + + if refine_on_receivers: + mesh = OctreeDriver.refine_tree_from_points( + mesh, survey.vertices, levels=[2], finalize=False + ) + + mesh.finalize() + entity = treemesh_2_octree(survey.workspace, mesh, name=name) + + return entity diff --git a/simpeg_drivers/utils/synthetics/meshes/tensors.py b/simpeg_drivers/utils/synthetics/meshes/tensors.py new file mode 100644 index 00000000..145fb944 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/meshes/tensors.py @@ -0,0 +1,55 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py.data import IntegerData +from geoh5py.objects import CellObject, DrapeModel + +from simpeg_drivers.utils.utils import get_drape_model + + +def get_tensor_mesh( + survey: CellObject, + cell_size: tuple[float, float, float], + padding_distance: float, + line_id: int = 101, + name: str = "mesh", +) -> DrapeModel: + """ + Generate a tensor mesh and the colocated DrapeModel. + + :param survey: Survey object with vertices that define the core of the + tensor mesh. + :param cell_size: Tuple defining the cell size in all directions. + :param padding_distance: Distance to pad the mesh in all directions. + :param line_id: Chooses line from the survey to define the drape model. + + :return entity: The DrapeModel object that shares cells with the discretize + tensor mesh and which stores the results of computations. + :return mesh: The discretize tensor mesh object for computations. + """ + + line_data = survey.get_entity("line_ids")[0] + assert isinstance(line_data, IntegerData) + lines = line_data.values + entity, mesh, _ = get_drape_model( # pylint: disable=unbalanced-tuple-unpacking + survey.workspace, + name, + survey.vertices[np.unique(survey.cells[lines == line_id, :]), :], + [cell_size[0], cell_size[2]], + 100.0, + [padding_distance] * 2 + [padding_distance] * 2, + 1.1, + parent=None, + return_colocated_mesh=True, + return_sorting=True, + ) + + return entity diff --git a/simpeg_drivers/utils/synthetics/models.py b/simpeg_drivers/utils/synthetics/models.py new file mode 100644 index 00000000..0ba1b219 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/models.py @@ -0,0 +1,53 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoapps_utils.modelling.plates import make_plate +from geoh5py.data import FloatData +from geoh5py.objects import DrapeModel, Octree + +from simpeg_drivers.utils.synthetics.options import ModelOptions + + +def get_model( + method: str, + mesh: Octree | DrapeModel, + active: np.ndarray, + options: ModelOptions, +) -> FloatData: + """ + Create a halfspace and plate model in the cell centers of the provided mesh. + + :param method: The geophysical method controlling the factory behaviour + :param mesh: The mesh whose cell centers the model will be defined on. + :param plate: The plate object defining the location and orientation of the + plate anomaly. + :param background: Value given to the halfspace. + :param anomaly: Value given to the plate. + """ + + cell_centers = mesh.centroids.copy() + + model = make_plate( + points=cell_centers, + plate=options.plate, + background=options.background, + anomaly=options.anomaly, + ) + + if "1d" in method: + model = options.background * np.ones(mesh.n_cells) + inside_anomaly = (mesh.centroids[:, 2] < 0) & (mesh.centroids[:, 2] > -20) + model[inside_anomaly] = options.anomaly + + model[~active] = np.nan + model = mesh.add_data({options.name: {"values": model}}) + assert isinstance(model, FloatData) + return model diff --git a/simpeg_drivers/utils/synthetics/options.py b/simpeg_drivers/utils/synthetics/options.py new file mode 100644 index 00000000..de3cfdd5 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/options.py @@ -0,0 +1,68 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from collections.abc import Callable + +from geoapps_utils.modelling.plates import PlateModel +from geoapps_utils.utils.locations import gaussian +from pydantic import BaseModel, ConfigDict + + +class SurveyOptions(BaseModel): + center: tuple[float, float] = (0.0, 0.0) + width: float = 200.0 + height: float = 200.0 + drape: float = 0.0 + n_stations: int = 20 + n_lines: int = 5 + topography: Callable = lambda x, y: gaussian(x, y, amplitude=50.0, width=100.0) + name: str = "survey" + + @property + def limits(self) -> list[float]: + return [ + self.center[0] - self.width / 2, + self.center[0] + self.width / 2, + self.center[1] - self.height / 2, + self.center[1] + self.height / 2, + ] + + +class MeshOptions(BaseModel): + cell_size: tuple[float, float, float] = (5.0, 5.0, 5.0) + refinement: tuple = (4, 6) + padding_distance: float = 100.0 + name: str = "mesh" + + +class ModelOptions(BaseModel): + background: float = 0.0 + anomaly: float = 1.0 + plate: PlateModel = PlateModel( + strike_length=40.0, + dip_length=40.0, + width=40.0, + origin=(0.0, 0.0, 10.0), + ) + name: str = "model" + + +class ActiveCellsOptions(BaseModel): + name: str = "active_cells" + + +class SyntheticsComponentsOptions(BaseModel): + model_config = ConfigDict(arbitrary_types_allowed=True) + + method: str = "gravity" + survey: SurveyOptions = SurveyOptions() + mesh: MeshOptions = MeshOptions() + model: ModelOptions = ModelOptions() + active: ActiveCellsOptions = ActiveCellsOptions() diff --git a/simpeg_drivers/utils/synthetics/surveys/__init__.py b/simpeg_drivers/utils/synthetics/surveys/__init__.py new file mode 100644 index 00000000..76e3205e --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/utils/synthetics/surveys/dcip.py b/simpeg_drivers/utils/synthetics/surveys/dcip.py new file mode 100644 index 00000000..385b1cc0 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/dcip.py @@ -0,0 +1,79 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import CurrentElectrode, PotentialElectrode + + +def generate_dc_survey( + workspace: Workspace, + X: np.ndarray, + Y: np.ndarray, + Z: np.ndarray | None = None, + name: str = "survey", +) -> PotentialElectrode: + """Generate a DC survey object from survey grid locations.""" + + if Z is None: + Z = np.zeros_like(X) + + vertices = np.c_[X.ravel(), Y.ravel(), Z.ravel()] + parts = np.kron(np.arange(X.shape[0]), np.ones(X.shape[1])).astype("int") + currents = CurrentElectrode.create( + workspace, name=f"{name}_tx", vertices=vertices, parts=parts + ) + currents.add_default_ab_cell_id() + n_dipoles = 9 + dipoles = [] + current_id = [] + + for val in currents.ab_cell_id.values: + cell_id = int(currents.ab_map[val]) - 1 + + for dipole in range(n_dipoles): + dipole_ids = currents.cells[cell_id, :] + 2 + dipole + + if ( + any(dipole_ids > (currents.n_vertices - 1)) + or len( + np.unique(parts[np.r_[currents.cells[cell_id, 0], dipole_ids[1]]]) + ) + > 1 + ): + continue + + dipoles += [dipole_ids] + current_id += [val] + + potentials = PotentialElectrode.create( + workspace, + name=name, + vertices=vertices, + cells=np.vstack(dipoles).astype("uint32"), + ) + line_id = potentials.vertices[potentials.cells[:, 0], 1] + line_id = (line_id - np.min(line_id) + 1).astype(np.int32) + line_reference = {0: "Unknown"} + line_reference.update({k: str(k) for k in np.unique(line_id)}) + potentials.add_data( + { + "line_ids": { + "values": line_id, + "type": "REFERENCED", + "value_map": line_reference, + } + } + ) + potentials.ab_cell_id = np.hstack(current_id).astype("int32") + potentials.current_electrodes = currents + currents.potential_electrodes = potentials + + return potentials diff --git a/simpeg_drivers/utils/synthetics/surveys/factory.py b/simpeg_drivers/utils/synthetics/surveys/factory.py new file mode 100644 index 00000000..23c1082e --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/factory.py @@ -0,0 +1,93 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from collections.abc import Callable + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects import ObjectBase, Points + +from simpeg_drivers.utils.synthetics.options import SurveyOptions + +from .dcip import generate_dc_survey +from .frequency_domain.fdem import generate_fdem_survey +from .natural_sources.magnetotellurics import generate_magnetotellurics_survey +from .natural_sources.tipper import generate_tipper_survey +from .time_domain.airborne_tdem import generate_airborne_tdem_survey +from .time_domain.ground_tdem import generate_tdem_survey + + +def grid_layout( + limits: list[float], + station_spacing: int, + line_spacing: 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 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, Y = np.meshgrid(x, y) + Z = topography(X, Y) + + return X, Y, Z + + +def get_survey( + geoh5: Workspace, + method: str, + options: SurveyOptions, +) -> ObjectBase: + """ + Factory for survey creation with behaviour modified by the provided method. + + :param geoh5: The workspace to create the survey in. + :param method: The geophysical method controlling the factory behaviour. + :param options: Survey options. + """ + + X, Y, Z = grid_layout( + limits=options.limits, + station_spacing=options.n_stations, + line_spacing=options.n_lines, + topography=options.topography, + ) + Z += options.drape + + if "current" in method or "polarization" in method: + return generate_dc_survey(geoh5, X, Y, Z, name=options.name) + + if "magnetotellurics" in method: + return generate_magnetotellurics_survey(geoh5, X, Y, Z, name=options.name) + + if "tipper" in method: + return generate_tipper_survey(geoh5, X, Y, Z, name=options.name) + + if method in ["fdem", "fem", "fdem 1d"]: + return generate_fdem_survey(geoh5, X, Y, Z, name=options.name) + + if "tdem" in method: + if "airborne" in method: + return generate_airborne_tdem_survey(geoh5, X, Y, Z, name=options.name) + else: + return generate_tdem_survey(geoh5, X, Y, Z, name=options.name) + + return Points.create( + geoh5, + vertices=np.column_stack([X.flatten(), Y.flatten(), Z.flatten()]), + name=options.name, + ) diff --git a/simpeg_drivers/utils/synthetics/surveys/frequency_domain/__init__.py b/simpeg_drivers/utils/synthetics/surveys/frequency_domain/__init__.py new file mode 100644 index 00000000..76e3205e --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/frequency_domain/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/utils/synthetics/surveys/frequency_domain/fdem.py b/simpeg_drivers/utils/synthetics/surveys/frequency_domain/fdem.py new file mode 100644 index 00000000..13a62705 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/frequency_domain/fdem.py @@ -0,0 +1,64 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoapps_utils.utils.locations import mask_large_connections +from geoh5py import Workspace +from geoh5py.objects import AirborneFEMReceivers, AirborneFEMTransmitters + + +frequency_config = [ + {"Coaxial data": False, "Frequency": 900, "Offset": 7.86}, + {"Coaxial data": False, "Frequency": 7200, "Offset": 7.86}, + {"Coaxial data": False, "Frequency": 56000, "Offset": 6.3}, +] + + +def generate_fdem_survey( + geoh5: Workspace, X: np.ndarray, Y: np.ndarray, Z: np.ndarray, name: str = "survey" +) -> AirborneFEMReceivers: + """Create an FDEM survey object from survey grid locations.""" + + vertices = np.column_stack([X.flatten(), Y.flatten(), Z.flatten()]) + survey = AirborneFEMReceivers.create(geoh5, vertices=vertices, name=name) + + survey.metadata["EM Dataset"]["Frequency configurations"] = frequency_config + + tx_locs_list = [] + frequency_list = [] + for config in frequency_config: + tx_vertices = vertices.copy() + tx_vertices[:, 0] -= config["Offset"] + tx_locs_list.append(tx_vertices) + frequency_list.append([[config["Frequency"]] * len(vertices)]) + tx_locs = np.vstack(tx_locs_list) + freqs = np.hstack(frequency_list).flatten() + + transmitters = AirborneFEMTransmitters.create( + geoh5, vertices=tx_locs, name=f"{name}_tx" + ) + survey.transmitters = transmitters + survey.channels = [float(k["Frequency"]) for k in frequency_config] + + transmitters.add_data( + { + "Tx frequency": { + "values": freqs, + "association": "VERTEX", + "primitive_type": "REFERENCED", + "value_map": {k: str(k) for k in freqs}, + } + } + ) + + survey.remove_cells(mask_large_connections(survey, 200.0)) + transmitters.remove_cells(mask_large_connections(transmitters, 200.0)) + + return survey diff --git a/simpeg_drivers/utils/synthetics/surveys/layout.py b/simpeg_drivers/utils/synthetics/surveys/layout.py new file mode 100644 index 00000000..c8aa1763 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/layout.py @@ -0,0 +1,13 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +from collections.abc import Callable + +import numpy as np diff --git a/simpeg_drivers/utils/synthetics/surveys/natural_sources/__init__.py b/simpeg_drivers/utils/synthetics/surveys/natural_sources/__init__.py new file mode 100644 index 00000000..4d06f672 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/natural_sources/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/utils/synthetics/surveys/natural_sources/magnetotellurics.py b/simpeg_drivers/utils/synthetics/surveys/natural_sources/magnetotellurics.py new file mode 100644 index 00000000..3268f5d8 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/natural_sources/magnetotellurics.py @@ -0,0 +1,42 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +from geoh5py.objects.surveys.electromagnetics.magnetotellurics import MTReceivers + + +def generate_magnetotellurics_survey( + geoh5: Workspace, + X: np.ndarray, + Y: np.ndarray, + Z: np.ndarray, + channels: tuple = (10.0, 100.0, 1000.0), + name: str = "survey", +) -> MTReceivers: + """Create a Magnetotellurics survey object from survey grid locations.""" + + survey = MTReceivers.create( + geoh5, + vertices=np.column_stack([X.flatten(), Y.flatten(), Z.flatten()]), + name=name, + components=[ + "Zxx (real)", + "Zxx (imag)", + "Zxy (real)", + "Zxy (imag)", + "Zyx (real)", + "Zyx (imag)", + "Zyy (real)", + "Zyy (imag)", + ], + channels=list(channels), + ) + return survey diff --git a/simpeg_drivers/utils/synthetics/surveys/natural_sources/tipper.py b/simpeg_drivers/utils/synthetics/surveys/natural_sources/tipper.py new file mode 100644 index 00000000..5ba6ec01 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/natural_sources/tipper.py @@ -0,0 +1,47 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoapps_utils.utils.locations import mask_large_connections +from geoh5py import Workspace +from geoh5py.objects.surveys.electromagnetics.tipper import ( + TipperBaseStations, + TipperReceivers, +) + + +def generate_tipper_survey( + geoh5: Workspace, + X: np.ndarray, + Y: np.ndarray, + Z: np.ndarray, + channels: tuple = (10.0, 100.0, 1000.0), + name: str = "survey", +) -> TipperReceivers: + """Create a Tipper survey object from survey grid locations.""" + vertices = np.column_stack([X.flatten(), Y.flatten(), Z.flatten()]) + survey = TipperReceivers.create( + geoh5, + vertices=vertices, + name=name, + components=[ + "Txz (real)", + "Txz (imag)", + "Tyz (real)", + "Tyz (imag)", + ], + channels=list(channels), + ) + survey.base_stations = TipperBaseStations.create( + geoh5, vertices=np.c_[vertices[0, :]].T + ) + survey.remove_cells(mask_large_connections(survey, 200.0)) + + return survey diff --git a/simpeg_drivers/utils/synthetics/surveys/time_domain/__init__.py b/simpeg_drivers/utils/synthetics/surveys/time_domain/__init__.py new file mode 100644 index 00000000..1f6de18c --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/time_domain/__init__.py @@ -0,0 +1,24 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np + + +CHANNELS = np.r_[3e-04, 6e-04, 1.2e-03] * 1e3 +WAVEFORM = np.c_[ + np.r_[ + np.arange(-0.002, -0.0001, 5e-4), + np.arange(-0.0004, 0.0, 1e-4), + np.arange(0.0, 0.002, 5e-4), + ] + * 1e3 + + 2.0, + np.r_[np.linspace(0, 1, 4), np.linspace(0.9, 0.0, 4), np.zeros(4)], +] diff --git a/simpeg_drivers/utils/synthetics/surveys/time_domain/airborne_tdem.py b/simpeg_drivers/utils/synthetics/surveys/time_domain/airborne_tdem.py new file mode 100644 index 00000000..1fa96adb --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/time_domain/airborne_tdem.py @@ -0,0 +1,47 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoapps_utils.utils.locations import mask_large_connections +from geoh5py import Workspace +from geoh5py.objects import ( + AirborneTEMReceivers, + AirborneTEMTransmitters, +) + +from simpeg_drivers.utils.synthetics.surveys.time_domain import CHANNELS, WAVEFORM + + +def generate_airborne_tdem_survey( + geoh5: Workspace, + X: np.ndarray, + Y: np.ndarray, + Z: np.ndarray, + channels: np.ndarray = CHANNELS, + waveform: np.ndarray = WAVEFORM, + name: str = "survey", +) -> AirborneTEMReceivers: + """Create an Airborne TDEM survey object from survey grid locations""" + vertices = np.column_stack([X.flatten(), Y.flatten(), Z.flatten()]) + survey = AirborneTEMReceivers.create(geoh5, vertices=vertices, name=name) + transmitters = AirborneTEMTransmitters.create( + geoh5, vertices=vertices, name=f"{name}_tx" + ) + mask = mask_large_connections(survey, 200.0) + survey.remove_cells(mask) + transmitters.remove_cells(mask) + + survey.transmitters = transmitters + survey.channels = channels + survey.waveform = waveform + survey.timing_mark = 2.0 + survey.unit = "Milliseconds (ms)" + + return survey diff --git a/simpeg_drivers/utils/synthetics/surveys/time_domain/ground_tdem.py b/simpeg_drivers/utils/synthetics/surveys/time_domain/ground_tdem.py new file mode 100644 index 00000000..146346cd --- /dev/null +++ b/simpeg_drivers/utils/synthetics/surveys/time_domain/ground_tdem.py @@ -0,0 +1,108 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of geoapps-utils package. ' +# ' +# geoapps-utils is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoapps_utils.utils.locations import gaussian +from geoh5py import Workspace +from geoh5py.objects import ( + LargeLoopGroundTEMReceivers, + LargeLoopGroundTEMTransmitters, +) + +from simpeg_drivers.utils.synthetics.surveys.time_domain import CHANNELS, WAVEFORM + + +def generate_tdem_survey( + geoh5: Workspace, + X: np.ndarray, + Y: np.ndarray, + Z: np.ndarray, + channels: np.ndarray = CHANNELS, + waveform: np.ndarray = WAVEFORM, + name: str = "survey", +) -> LargeLoopGroundTEMReceivers: + """Create a large loop TDEM survey object from survey grid locations.""" + + flatten = len(np.unique(Z)) == 1 + vertices = np.column_stack([X.flatten(), Y.flatten(), Z.flatten()]) + center = np.mean(vertices, axis=0) + if flatten: + center[2] -= np.unique(Z) + n_lines = X.shape[0] + arrays = [ + np.c_[ + X[: int(n_lines / 2), :].flatten(), + Y[: int(n_lines / 2), :].flatten(), + Z[: int(n_lines / 2), :].flatten(), + ], + np.c_[ + X[int(n_lines / 2) :, :].flatten(), + Y[int(n_lines / 2) :, :].flatten(), + Z[int(n_lines / 2) :, :].flatten(), + ], + ] + loops = [] + loop_cells = [] + loop_id = [] + count = 0 + for ind, array in enumerate(arrays): + loop_id += [np.ones(array.shape[0]) * (ind + 1)] + min_loc = np.min(array, axis=0) + max_loc = np.max(array, axis=0) + loop = np.vstack( + [ + np.c_[ + np.ones(5) * min_loc[0], + np.linspace(min_loc[1], max_loc[1], 5), + ], + np.c_[ + np.linspace(min_loc[0], max_loc[0], 5)[1:], + np.ones(4) * max_loc[1], + ], + np.c_[ + np.ones(4) * max_loc[0], + np.linspace(max_loc[1], min_loc[1], 5)[1:], + ], + np.c_[ + np.linspace(max_loc[0], min_loc[0], 5)[1:-1], + np.ones(3) * min_loc[1], + ], + ] + ) + loop = (loop - np.mean(loop, axis=0)) * 1.5 + np.mean(loop, axis=0) + elevation = ( + np.ones(len(loop)) * np.unique(Z) + if flatten + else gaussian(loop[:, 0], loop[:, 1], amplitude=50.0, width=100.0) + ) + loop = np.c_[loop, elevation] + loops += [loop + np.asarray(center)] + loop_cells += [np.c_[np.arange(15) + count, np.arange(15) + count + 1]] + loop_cells += [np.c_[count + 15, count]] + count += 16 + + transmitters = LargeLoopGroundTEMTransmitters.create( + geoh5, + vertices=np.vstack(loops), + cells=np.vstack(loop_cells), + name=f"{name}_tx", + ) + transmitters.tx_id_property = transmitters.parts + 1 + survey = LargeLoopGroundTEMReceivers.create(geoh5, name=name, vertices=vertices) + survey.transmitters = transmitters + survey.tx_id_property = np.hstack(loop_id) + + survey.channels = channels + + survey.waveform = waveform + survey.timing_mark = 2.0 + survey.unit = "Milliseconds (ms)" + + return survey diff --git a/simpeg_drivers/utils/synthetics/topography.py b/simpeg_drivers/utils/synthetics/topography.py new file mode 100644 index 00000000..98ba9236 --- /dev/null +++ b/simpeg_drivers/utils/synthetics/topography.py @@ -0,0 +1,52 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import numpy as np +from geoh5py import Workspace +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.surveys.factory import grid_layout +from simpeg_drivers.utils.utils import active_from_xyz + + +def get_topography_surface(geoh5: Workspace, options: SurveyOptions) -> Surface: + """ + Returns a topography surface with 2x the limits of the survey. + + :param geoh5: Geoh5 workspace. + :param options: Survey options. Extents will be 2x the survey extents. + """ + + X, Y, Z = grid_layout( + limits=[2 * k for k in options.limits], # type: ignore + station_spacing=int(np.ceil((options.limits[1] - options.limits[0]) / 4)), + line_spacing=int(np.ceil((options.limits[3] - options.limits[2]) / 4)), + topography=options.topography, + ) + + vertices = np.column_stack( + [X.flatten(order="F"), Y.flatten(order="F"), Z.flatten(order="F")] + ) + + return Surface.create( + geoh5, + vertices=vertices, + cells=Delaunay(vertices[:, :2]).simplices, # pylint: disable=no-member + name="topography", + ) + + +def get_active( + mesh: Octree | DrapeModel, topography: Surface, name: str = "active_cells" +): + active = active_from_xyz(mesh, topography.vertices, grid_reference="top") + return mesh.add_data({name: {"values": active}}) diff --git a/simpeg_drivers/utils/utils.py b/simpeg_drivers/utils/utils.py index 425a2ae4..56f22591 100644 --- a/simpeg_drivers/utils/utils.py +++ b/simpeg_drivers/utils/utils.py @@ -14,7 +14,6 @@ import warnings from copy import deepcopy from typing import TYPE_CHECKING -from uuid import UUID import numpy as np from discretize import TensorMesh, TreeMesh @@ -436,31 +435,6 @@ def get_drape_model( return val -def get_inversion_output(h5file: str | Workspace, inversion_group: str | UUID): - """ - Recover inversion iterations from a ContainerGroup comments. - """ - if isinstance(h5file, Workspace): - workspace = h5file - else: - workspace = Workspace(h5file) - - try: - group = workspace.get_entity(inversion_group)[0] - except IndexError as exc: - raise IndexError( - f"BaseInversion group {inversion_group} could not be found in the target geoh5 {h5file}" - ) from exc - - outfile = group.get_entity("SimPEG.out")[0] - out = list(outfile.file_bytes.decode("utf-8").replace("\r", "").split("\n"))[:-1] - cols = out.pop(0).split(" ") - out = [[string_to_numeric(k) for k in elem.split(" ")] for elem in out] - out = dict(zip(cols, list(map(list, zip(*out, strict=True))), strict=True)) - - return out - - def xyz_2_drape_model( workspace, locations, depths, name=None, parent=None ) -> DrapeModel: diff --git a/tests/data_test.py b/tests/data_test.py index 851cdad8..d3aaa81e 100644 --- a/tests/data_test.py +++ b/tests/data_test.py @@ -31,36 +31,42 @@ from simpeg_drivers.potential_fields.magnetic_vector.options import ( MVIInversionOptions, ) -from tests.testing_utils import setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) def get_mvi_params(tmp_path: Path, **kwargs) -> MVIInversionOptions: - geoh5, entity, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=(2,), - n_electrodes=2, - n_lines=2, - inversion_type="magnetic_vector", - ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions(n_stations=2, n_lines=2), + mesh=MeshOptions(refinement=(2,)), + model=ModelOptions(anomaly=0.05), + ) + components = SyntheticsComponents(geoh5=geoh5, options=opts) + with geoh5.open(): - tmi_channel = survey.add_data( - {"tmi": {"values": np.random.rand(survey.n_vertices)}} + tmi_channel = components.survey.add_data( + {"tmi": {"values": np.random.rand(components.survey.n_vertices)}} + ) + params = MVIInversionOptions.build( + geoh5=geoh5, + data_object=components.survey, + tmi_channel=tmi_channel, + tmi_uncertainty=1.0, + topography_object=components.topography, + mesh=components.model.parent, + starting_model=components.model, + inducing_field_strength=50000.0, + inducing_field_inclination=60.0, + inducing_field_declination=30.0, + **kwargs, ) - params = MVIInversionOptions.build( - geoh5=geoh5, - data_object=survey, - tmi_channel=tmi_channel, - tmi_uncertainty=1.0, - topography_object=topography, - mesh=model.parent, - starting_model=model, - inducing_field_strength=50000.0, - inducing_field_inclination=60.0, - inducing_field_declination=30.0, - **kwargs, - ) return params @@ -242,23 +248,20 @@ def test_get_survey(tmp_path: Path): def test_data_parts(tmp_path: Path): n_lines = 8 - geoh5, entity, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.01, - anomaly=10, - n_electrodes=10, - n_lines=n_lines, - drape_height=0.0, - inversion_type="direct current 3d", - flatten=False, + opts = SyntheticsComponentsOptions( + method="direct current 3d", + survey=SurveyOptions(n_stations=10, n_lines=n_lines), + mesh=MeshOptions(), + model=ModelOptions(background=0.01, anomaly=10.0), ) - with geoh5.open(): + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) params = DC3DForwardOptions.build( geoh5=geoh5, - data_object=survey, - topography_object=topography, - mesh=model.parent, - starting_model=model, + data_object=components.survey, + topography_object=components.topography, + mesh=components.model.parent, + starting_model=components.model, ) data = InversionData(geoh5, params) diff --git a/tests/driver_test.py b/tests/driver_test.py index 72ca8629..60c413a1 100644 --- a/tests/driver_test.py +++ b/tests/driver_test.py @@ -11,35 +11,41 @@ from pathlib import Path import numpy as np +from geoh5py import Workspace -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import GravityInversionOptions from simpeg_drivers.potential_fields.gravity.driver import GravityInversionDriver -from tests.testing_utils import setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) def test_smallness_terms(tmp_path: Path): n_grid_points = 2 refinement = (2,) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - flatten=False, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions(n_stations=n_grid_points, n_lines=n_grid_points), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.75), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - with geoh5.open(): - gz = survey.add_data({"gz": {"values": np.ones(survey.n_vertices)}}) - mesh = model.parent + gz = components.survey.add_data( + {"gz": {"values": np.ones(components.survey.n_vertices)}} + ) + mesh = components.model.parent params = GravityInversionOptions.build( geoh5=geoh5, mesh=mesh, - topography_object=topography, + topography_object=components.topography, data_object=gz.parent, starting_model=1e-4, reference_model=0.0, @@ -61,23 +67,22 @@ def test_target_chi(tmp_path: Path, caplog): n_grid_points = 2 refinement = (2,) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - flatten=False, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions(n_station=n_grid_points, n_lines=n_grid_points), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.75), ) - - with geoh5.open(): - gz = survey.add_data({"gz": {"values": np.ones(survey.n_vertices)}}) - mesh = model.parent + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + gz = components.survey.add_data( + {"gz": {"values": np.ones(components.survey.n_vertices)}} + ) + mesh = components.model.parent params = GravityInversionOptions.build( geoh5=geoh5, mesh=mesh, - topography_object=topography, + topography_object=components.topography, data_object=gz.parent, gz_channel=gz, gz_uncertainty=2e-3, diff --git a/tests/locations_test.py b/tests/locations_test.py index 948e8825..efdd6230 100644 --- a/tests/locations_test.py +++ b/tests/locations_test.py @@ -16,42 +16,44 @@ import pytest from geoh5py import Workspace from geoh5py.objects import Curve, Grid2D, Points -from scipy.spatial import ConvexHull from simpeg_drivers.components.locations import InversionLocations -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import MVIInversionOptions +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) from simpeg_drivers.utils.utils import tile_locations -from tests.testing_utils import Geoh5Tester, setup_inversion_workspace def get_mvi_params(tmp_path: Path) -> MVIInversionOptions: - geoh5, enitiy, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=(2,), - n_electrodes=2, - n_lines=2, - inversion_type="magnetic_vector", + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions(n_lines=2, n_stations=2), + mesh=MeshOptions(refinement=(2,)), + model=ModelOptions(background=0.0, anomaly=0.05), ) - with geoh5.open(): - tmi_channel = survey.add_data( - {"tmi": {"values": np.random.rand(survey.n_vertices)}} + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + tmi_channel = components.survey.add_data( + {"tmi": {"values": np.random.rand(components.survey.n_vertices)}} ) - params = MVIInversionOptions.build( - geoh5=geoh5, - data_object=survey, - tmi_channel=tmi_channel, - tmi_uncertainty=1.0, - topography_object=topography, - mesh=model.parent, - starting_model=model, - inducing_field_strength=50000.0, - inducing_field_inclination=60.0, - inducing_field_declination=30.0, - ) - return params + params = MVIInversionOptions.build( + geoh5=geoh5, + data_object=components.survey, + tmi_channel=tmi_channel, + tmi_uncertainty=1.0, + topography_object=components.topography, + mesh=components.model.parent, + starting_model=components.model, + inducing_field_strength=50000.0, + inducing_field_inclination=60.0, + inducing_field_declination=30.0, + ) + return params def test_mask(tmp_path: Path): diff --git a/tests/meshes_test.py b/tests/meshes_test.py index d492c1c8..abcbb37f 100644 --- a/tests/meshes_test.py +++ b/tests/meshes_test.py @@ -22,45 +22,49 @@ from simpeg_drivers.components import InversionMesh from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import MVIInversionOptions -from tests.testing_utils import Geoh5Tester, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) def get_mvi_params(tmp_path: Path) -> MVIInversionOptions: - geoh5, entity, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=(2,), - n_electrodes=4, - n_lines=4, - inversion_type="magnetic_vector", + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions(n_stations=4, n_lines=4), + mesh=MeshOptions(refinement=(2,)), + model=ModelOptions(anomaly=0.05), ) - with geoh5.open(): - mesh = model.parent - tmi_channel, gyz_channel = survey.add_data( + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + mesh = components.model.parent + tmi_channel, gyz_channel = components.survey.add_data( { - "tmi": {"values": np.random.rand(survey.n_vertices)}, - "gyz": {"values": np.random.rand(survey.n_vertices)}, + "tmi": {"values": np.random.rand(components.survey.n_vertices)}, + "gyz": {"values": np.random.rand(components.survey.n_vertices)}, } ) - elevation = topography.add_data( - {"elevation": {"values": topography.vertices[:, 2]}} + elevation = components.topography.add_data( + {"elevation": {"values": components.topography.vertices[:, 2]}} ) - params = MVIInversionOptions.build( - geoh5=geoh5, - data_object=survey, - tmi_channel=tmi_channel, - tmi_uncertainty=0.01, - active_cells=ActiveCellsOptions( - topography_object=topography, topography=elevation - ), - inducing_field_strength=50000.0, - inducing_field_inclination=60.0, - inducing_field_declination=30.0, - mesh=mesh, - starting_model=model, - ) + params = MVIInversionOptions.build( + geoh5=geoh5, + data_object=components.survey, + tmi_channel=tmi_channel, + tmi_uncertainty=0.01, + active_cells=ActiveCellsOptions( + topography_object=components.topography, topography=elevation + ), + inducing_field_strength=50000.0, + inducing_field_inclination=60.0, + inducing_field_declination=30.0, + mesh=mesh, + starting_model=components.model, + ) return params diff --git a/tests/models_test.py b/tests/models_test.py index 0352c22c..1a6a03d6 100644 --- a/tests/models_test.py +++ b/tests/models_test.py @@ -13,6 +13,7 @@ from pathlib import Path import numpy as np +from geoh5py import Workspace from geoh5py.objects import Points from simpeg_drivers.components import ( @@ -25,52 +26,56 @@ from simpeg_drivers.potential_fields.magnetic_vector.driver import ( MVIInversionDriver, ) -from tests.testing_utils import Geoh5Tester, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) def get_mvi_params(tmp_path: Path) -> MVIInversionOptions: - geoh5, entity, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=(2,), - n_electrodes=2, - n_lines=2, - inversion_type="magnetic_vector", + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions(n_stations=2, n_lines=2), + mesh=MeshOptions(refinement=(2,)), + model=ModelOptions(anomaly=0.05), ) - with geoh5.open(): - mesh = model.parent + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + mesh = components.model.parent ref_inducing = mesh.add_data( { "reference_inclination": {"values": np.ones(mesh.n_cells) * 79.0}, "reference_declination": {"values": np.ones(mesh.n_cells) * 11.0}, } ) - tmi_channel = survey.add_data( + tmi_channel = components.survey.add_data( { - "tmi": {"values": np.random.rand(survey.n_vertices)}, + "tmi": {"values": np.random.rand(components.survey.n_vertices)}, } ) - elevation = topography.add_data( - {"elevation": {"values": topography.vertices[:, 2]}} + elevation = components.topography.add_data( + {"elevation": {"values": components.topography.vertices[:, 2]}} + ) + params = MVIInversionOptions.build( + geoh5=geoh5, + data_object=components.survey, + tmi_channel=tmi_channel, + tmi_uncertainty=1.0, + mesh=mesh, + active_cells=ActiveCellsOptions( + topography_object=components.topography, topography=elevation + ), + starting_model=1e-04, + inducing_field_inclination=79.0, + inducing_field_declination=11.0, + reference_model=0.0, + inducing_field_strength=50000.0, + reference_inclination=ref_inducing[0], + reference_declination=ref_inducing[1], ) - params = MVIInversionOptions.build( - geoh5=geoh5, - data_object=survey, - tmi_channel=tmi_channel, - tmi_uncertainty=1.0, - mesh=mesh, - active_cells=ActiveCellsOptions( - topography_object=topography, topography=elevation - ), - starting_model=1e-04, - inducing_field_inclination=79.0, - inducing_field_declination=11.0, - reference_model=0.0, - inducing_field_strength=50000.0, - reference_inclination=ref_inducing[0], - reference_declination=ref_inducing[1], - ) return params diff --git a/tests/plate_simulation/runtest/driver_test.py b/tests/plate_simulation/runtest/driver_test.py index b85d7736..351503f7 100644 --- a/tests/plate_simulation/runtest/driver_test.py +++ b/tests/plate_simulation/runtest/driver_test.py @@ -8,54 +8,58 @@ # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -import numpy as np from geoh5py import Workspace from geoh5py.groups import SimPEGGroup -from geoh5py.objects import AirborneTEMReceivers, ObjectBase, Octree, Surface from geoh5py.ui_json import InputFile from simpeg_drivers import assets_path -from simpeg_drivers.electromagnetics.time_domain.options import TDEMForwardOptions from simpeg_drivers.plate_simulation.driver import ( - PlateSimulationDriver, PlateSimulationOptions, ) from simpeg_drivers.plate_simulation.models.options import ModelOptions from simpeg_drivers.plate_simulation.options import MeshOptions from simpeg_drivers.potential_fields.gravity.options import GravityForwardOptions -from tests.testing_utils import setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions as SyntheticsMeshOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + ModelOptions as SyntheticsModelOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + SurveyOptions, + SyntheticsComponentsOptions, +) # pylint: disable=too-many-statements def test_plate_simulation_params_from_input_file(tmp_path): - geoh5, mesh, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.0, - n_electrodes=8, - n_lines=8, - inversion_type="gravity", - flatten=False, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions(n_stations=8, n_lines=8), + mesh=SyntheticsMeshOptions(), + model=SyntheticsModelOptions(anomaly=0.0), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - with geoh5.open() as ws: ifile = InputFile.read_ui_json( assets_path() / "uijson" / "plate_simulation.ui.json", validate=False ) ifile.data["name"] = "test_gravity_plate_simulation" - ifile.data["geoh5"] = ws + ifile.data["geoh5"] = geoh5 # Add simulation parameter - gravity_inversion = SimPEGGroup.create(ws) + gravity_inversion = SimPEGGroup.create(geoh5) options = GravityForwardOptions.model_construct() fwr_ifile = InputFile.read_ui_json(options.default_ui_json) options_dict = fwr_ifile.ui_json options_dict["inversion_type"] = "gravity" options_dict["forward_only"] = True - options_dict["geoh5"] = str(ws.h5file) - options_dict["topography_object"]["value"] = str(topography.uid) - options_dict["data_object"]["value"] = str(survey.uid) + options_dict["geoh5"] = str(geoh5.h5file) + options_dict["topography_object"]["value"] = str(components.topography.uid) + options_dict["data_object"]["value"] = str(components.survey.uid) gravity_inversion.options = options_dict ifile.data["simulation"] = gravity_inversion @@ -95,9 +99,12 @@ def test_plate_simulation_params_from_input_file(tmp_path): assert simulation_parameters.inversion_type == "gravity" assert simulation_parameters.forward_only - assert simulation_parameters.geoh5.h5file == ws.h5file - assert simulation_parameters.active_cells.topography_object.uid == topography.uid - assert simulation_parameters.data_object.uid == survey.uid + assert simulation_parameters.geoh5.h5file == geoh5.h5file + assert ( + simulation_parameters.active_cells.topography_object.uid + == components.topography.uid + ) + assert simulation_parameters.data_object.uid == components.survey.uid assert isinstance(params.mesh, MeshOptions) assert params.mesh.u_cell_size == 10.0 diff --git a/tests/plate_simulation/runtest/gravity_test.py b/tests/plate_simulation/runtest/gravity_test.py index 40d796ab..95ded860 100644 --- a/tests/plate_simulation/runtest/gravity_test.py +++ b/tests/plate_simulation/runtest/gravity_test.py @@ -20,21 +20,28 @@ ) from simpeg_drivers.plate_simulation.options import MeshOptions, PlateSimulationOptions from simpeg_drivers.potential_fields.gravity.options import GravityForwardOptions -from tests.testing_utils import setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions as SyntheticsMeshOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + ModelOptions as SyntheticsModelOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + SurveyOptions, + SyntheticsComponentsOptions, +) def test_gravity_plate_simulation(tmp_path): - geoh5, mesh, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.0, - n_electrodes=8, - n_lines=8, - inversion_type="gravity", - flatten=False, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions(n_stations=8, n_lines=8, drape=5.0), + mesh=SyntheticsMeshOptions(), + model=SyntheticsModelOptions(anomaly=0.0), ) - - with geoh5.open() as ws: + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) mesh_params = MeshOptions( u_cell_size=10.0, v_cell_size=10.0, @@ -66,19 +73,19 @@ def test_gravity_plate_simulation(tmp_path): ) options = GravityForwardOptions.build( - topography_object=topography, - data_object=survey, - geoh5=ws, + topography_object=components.topography, + data_object=components.survey, + geoh5=geoh5, starting_model=0.1, ) - gravity_inversion = SimPEGGroup.create(ws) + gravity_inversion = SimPEGGroup.create(geoh5) gravity_inversion.options = options.serialize() params = PlateSimulationOptions( title="test", run_command="run", - geoh5=ws, + geoh5=geoh5, mesh=mesh_params, model=model_params, simulation=gravity_inversion, diff --git a/tests/run_tests/driver_2d_rotated_gradients_test.py b/tests/run_tests/driver_2d_rotated_gradients_test.py index 2b3a41b1..affd14d7 100644 --- a/tests/run_tests/driver_2d_rotated_gradients_test.py +++ b/tests/run_tests/driver_2d_rotated_gradients_test.py @@ -11,14 +11,11 @@ from __future__ import annotations from pathlib import Path -from unittest.mock import patch import numpy as np from geoapps_utils.modelling.plates import PlateModel -from geoapps_utils.utils.importing import GeoAppsError from geoh5py.groups.property_group import PropertyGroup from geoh5py.workspace import Workspace -from pytest import raises from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( DC2DForwardDriver, @@ -29,12 +26,20 @@ DC2DInversionOptions, ) from simpeg_drivers.options import ( - ActiveCellsOptions, DrapeModelOptions, LineSelectionOptions, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -49,47 +54,49 @@ def test_dc2d_rotated_grad_fwr_run( n_lines=3, refinement=(4, 6), ): - plate_model = PlateModel( - strike_length=1000.0, - dip_length=150.0, - width=20.0, - origin=(50.0, 0.0, -30), - direction=90, - dip=45, - ) + filepath = Path(tmp_path) / "inversion_test.ui.geoh5" + with Workspace.create(filepath) as geoh5: + # Run the forward + components = SyntheticsComponents( + geoh5=geoh5, + options=SyntheticsComponentsOptions( + method="direct current 2d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions( + background=0.01, + anomaly=10.0, + plate=PlateModel( + strike_length=1000.0, + dip_length=150.0, + width=20.0, + origin=(50.0, 0.0, -30), + direction=90, + dip=45, + ), + ), + ), + ) - # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model=plate_model, - background=0.01, - anomaly=10.0, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - inversion_type="direct current 2d", - drape_height=0.0, - flatten=False, - ) - line_selection = LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0], - line_id=101, - ) - params = DC2DForwardOptions.build( - geoh5=geoh5, - data_object=survey, - line_selection=line_selection, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - horizontal_padding=100.0, - vertical_padding=100.0, - expansion_factor=1.1, - ), - starting_model=model, - topography_object=topography, - ) + line_selection = LineSelectionOptions( + line_object=components.survey.get_data("line_ids")[0], + line_id=101, + ) + params = DC2DForwardOptions.build( + geoh5=geoh5, + data_object=components.survey, + line_selection=line_selection, + drape_model=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=100.0, + horizontal_padding=100.0, + vertical_padding=100.0, + expansion_factor=1.1, + ), + starting_model=components.model, + topography_object=components.topography, + ) fwr_driver = DC2DForwardDriver(params) fwr_driver.run() @@ -109,7 +116,7 @@ def test_dc2d_rotated_grad_run( with Workspace(workpath) as geoh5: potential = geoh5.get_entity("Iteration_0_dc")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5) orig_potential = potential.values.copy() mesh = geoh5.get_entity("Models")[0] @@ -140,7 +147,7 @@ def test_dc2d_rotated_grad_run( vertical_padding=100.0, expansion_factor=1.1, ), - topography_object=topography, + topography_object=components.topography, line_selection=LineSelectionOptions( line_object=geoh5.get_entity("line_ids")[0], line_id=101, diff --git a/tests/run_tests/driver_airborne_fem_1d_test.py b/tests/run_tests/driver_airborne_fem_1d_test.py index 1f676489..9bf6b052 100644 --- a/tests/run_tests/driver_airborne_fem_1d_test.py +++ b/tests/run_tests/driver_airborne_fem_1d_test.py @@ -26,9 +26,19 @@ FDEM1DForwardOptions, FDEM1DInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -44,28 +54,27 @@ def test_fem_fwr_1d_run( cell_size=(20.0, 20.0, 20.0), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=1e-4, - anomaly=0.1, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - drape_height=10.0, - cell_size=cell_size, - padding_distance=400, - inversion_type="fdem 1d", - flatten=False, - ) - params = FDEM1DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - z_real_channel_bool=True, - z_imag_channel_bool=True, + opts = SyntheticsComponentsOptions( + method="fdem 1d", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=10.0 + ), + mesh=MeshOptions( + cell_size=cell_size, refinement=refinement, padding_distance=400.0 + ), + model=ModelOptions(background=1e-4, anomaly=0.1), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5=geoh5, options=opts) + params = FDEM1DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + z_real_channel_bool=True, + z_imag_channel_bool=True, + ) fwr_driver = FDEM1DForwardDriver(params) fwr_driver.run() @@ -77,32 +86,26 @@ def test_fem_1d_run(tmp_path: Path, max_iterations=1, pytest=True): workpath = tmp_path.parent / "test_fem_fwr_1d_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: - survey = next( - child - for child in geoh5.get_entity("Airborne_rx") - if not isinstance(child.parent, SimPEGGroup) - ) - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5) data = {} uncertainties = {} - components = { + channels = { "z_real": "z_real", "z_imag": "z_imag", } - for comp, cname in components.items(): + for chan, cname in channels.items(): data[cname] = [] uncertainties[f"{cname} uncertainties"] = [] - for ind, freq in enumerate(survey.channels): - data_entity = geoh5.get_entity(f"Iteration_0_{comp}_[{ind}]")[0].copy( - parent=survey + for ind, freq in enumerate(components.survey.channels): + data_entity = geoh5.get_entity(f"Iteration_0_{chan}_[{ind}]")[0].copy( + parent=components.survey ) data[cname].append(data_entity) abs_val = np.abs(data_entity.values) - uncert = survey.add_data( + uncert = components.survey.add_data( { - f"uncertainty_{comp}_[{ind}]": { + f"uncertainty_{chan}_[{ind}]": { "values": np.ones_like(abs_val) * freq / 200.0 # * 2**(np.abs(ind-1)) @@ -110,27 +113,27 @@ def test_fem_1d_run(tmp_path: Path, max_iterations=1, pytest=True): } ) uncertainties[f"{cname} uncertainties"].append( - uncert.copy(parent=survey) + uncert.copy(parent=components.survey) ) - data_groups = survey.add_components_data(data) - uncert_groups = survey.add_components_data(uncertainties) + data_groups = components.survey.add_components_data(data) + uncert_groups = components.survey.add_components_data(uncertainties) data_kwargs = {} - for comp, data_group, uncert_group in zip( - components, data_groups, uncert_groups, strict=True + for chan, data_group, uncert_group in zip( + channels, data_groups, uncert_groups, strict=True ): - data_kwargs[f"{comp}_channel"] = data_group - data_kwargs[f"{comp}_uncertainty"] = uncert_group + data_kwargs[f"{chan}_channel"] = data_group + data_kwargs[f"{chan}_uncertainty"] = uncert_group orig_z_real_1 = geoh5.get_entity("Iteration_0_z_real_[0]")[0].values # Run the inverse params = FDEM1DInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, starting_model=1e-3, reference_model=1e-3, alpha_s=0.0, diff --git a/tests/run_tests/driver_airborne_tem_1d_test.py b/tests/run_tests/driver_airborne_tem_1d_test.py index 3f706843..9eb0244d 100644 --- a/tests/run_tests/driver_airborne_tem_1d_test.py +++ b/tests/run_tests/driver_airborne_tem_1d_test.py @@ -15,7 +15,6 @@ import numpy as np from geoh5py.groups import SimPEGGroup from geoh5py.workspace import Workspace -from pytest import raises from simpeg_drivers.electromagnetics.time_domain_1d.driver import ( TDEM1DForwardDriver, @@ -25,9 +24,19 @@ TDEM1DForwardOptions, TDEM1DInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -42,28 +51,30 @@ def test_airborne_tem_1d_fwr_run( cell_size=(20.0, 20.0, 20.0), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.1, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - cell_size=cell_size, - refinement=refinement, - inversion_type="airborne tdem 1d", - drape_height=10.0, - padding_distance=400.0, - flatten=False, - ) - params = TDEM1DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - z_channel_bool=True, - solver_type="Mumps", + opts = SyntheticsComponentsOptions( + method="airborne tdem 1d", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=10.0 + ), + mesh=MeshOptions( + cell_size=cell_size, refinement=refinement, padding_distance=400.0 + ), + model=ModelOptions(background=0.1), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents( + geoh5, + options=opts, + ) + params = TDEM1DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + z_channel_bool=True, + solver_type="Mumps", + ) fwr_driver = TDEM1DForwardDriver(params) @@ -80,31 +91,25 @@ def test_airborne_tem_1d_run(tmp_path: Path, max_iterations=1, pytest=True): ) with Workspace(workpath) as geoh5: - survey = next( - child - for child in geoh5.get_entity("Airborne_rx") - if not isinstance(child.parent, SimPEGGroup) - ) - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5=geoh5) data = {} uncertainties = {} - components = { + channels = { "z": "dBzdt", } - for comp, cname in components.items(): + for chan, cname in channels.items(): data[cname] = [] uncertainties[f"{cname} uncertainties"] = [] - for ii, _ in enumerate(survey.channels): - data_entity = geoh5.get_entity(f"Iteration_0_{comp}_[{ii}]")[0].copy( - parent=survey + for ii, _ in enumerate(components.survey.channels): + data_entity = geoh5.get_entity(f"Iteration_0_{chan}_[{ii}]")[0].copy( + parent=components.survey ) data[cname].append(data_entity) - uncert = survey.add_data( + uncert = components.survey.add_data( { - f"uncertainty_{comp}_[{ii}]": { + f"uncertainty_{chan}_[{ii}]": { "values": np.ones_like(data_entity.values) * (np.percentile(np.abs(data_entity.values), 10) / 2.0) } @@ -112,16 +117,16 @@ def test_airborne_tem_1d_run(tmp_path: Path, max_iterations=1, pytest=True): ) uncertainties[f"{cname} uncertainties"].append(uncert) - survey.add_components_data(data) - survey.add_components_data(uncertainties) + components.survey.add_components_data(data) + components.survey.add_components_data(uncertainties) data_kwargs = {} - for comp in components: - data_kwargs[f"{comp}_channel"] = survey.fetch_property_group( - name=f"dB{comp}dt" + for chan in channels: + data_kwargs[f"{chan}_channel"] = components.survey.fetch_property_group( + name=f"dB{chan}dt" ) - data_kwargs[f"{comp}_uncertainty"] = survey.fetch_property_group( - name=f"dB{comp}dt uncertainties" + data_kwargs[f"{chan}_uncertainty"] = components.survey.fetch_property_group( + name=f"dB{chan}dt uncertainties" ) orig_dBzdt = geoh5.get_entity("Iteration_0_z_[0]")[0].values @@ -129,9 +134,9 @@ def test_airborne_tem_1d_run(tmp_path: Path, max_iterations=1, pytest=True): # Run the inverse params = TDEM1DInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, starting_model=5e-1, reference_model=1e-1, s_norm=0.0, diff --git a/tests/run_tests/driver_airborne_tem_test.py b/tests/run_tests/driver_airborne_tem_test.py index 1fa77241..f29d0c4e 100644 --- a/tests/run_tests/driver_airborne_tem_test.py +++ b/tests/run_tests/driver_airborne_tem_test.py @@ -13,7 +13,6 @@ from pathlib import Path import numpy as np -from geoapps_utils.modelling.plates import PlateModel from geoh5py.groups import SimPEGGroup from geoh5py.workspace import Workspace from pytest import raises @@ -26,9 +25,19 @@ TDEMForwardOptions, TDEMInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -40,32 +49,30 @@ def test_bad_waveform(tmp_path: Path): n_grid_points = 3 refinement = (2,) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.001, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - inversion_type="airborne tdem", - drape_height=10.0, - padding_distance=400.0, - flatten=False, - ) - params = TDEMForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - x_channel_bool=True, - y_channel_bool=True, - z_channel_bool=True, + opts = SyntheticsComponentsOptions( + method="airborne tdem", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=10.0 + ), + mesh=MeshOptions(refinement=refinement, padding_distance=400.0), + model=ModelOptions(background=0.001), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = TDEMForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + x_channel_bool=True, + y_channel_bool=True, + z_channel_bool=True, + ) - fwr_driver = TDEMForwardDriver(params) + fwr_driver = TDEMForwardDriver(params) - survey.channels[-1] = 1000.0 + params.data_object.channels[-1] = 1000.0 with raises(ValueError, match="The latest time"): _ = fwr_driver.inversion_data.survey @@ -78,30 +85,29 @@ def test_airborne_tem_fwr_run( cell_size=(20.0, 20.0, 20.0), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.001, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - cell_size=cell_size, - refinement=refinement, - inversion_type="airborne tdem", - drape_height=10.0, - padding_distance=400.0, - flatten=False, - ) - params = TDEMForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - x_channel_bool=True, - y_channel_bool=True, - z_channel_bool=True, - solver_type="Mumps", + opts = SyntheticsComponentsOptions( + method="airborne tdem", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=10.0 + ), + mesh=MeshOptions( + cell_size=cell_size, refinement=refinement, padding_distance=400.0 + ), + model=ModelOptions(background=0.001), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = TDEMForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + x_channel_bool=True, + y_channel_bool=True, + z_channel_bool=True, + solver_type="Mumps", + ) fwr_driver = TDEMForwardDriver(params) @@ -116,32 +122,25 @@ def test_airborne_tem_run(tmp_path: Path, max_iterations=1, pytest=True): ) with Workspace(workpath) as geoh5: - survey = next( - child - for child in geoh5.get_entity("Airborne_rx") - if not isinstance(child.parent, SimPEGGroup) - ) - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] - + components = SyntheticsComponents(geoh5) data = {} uncertainties = {} - components = { + channels = { "z": "dBzdt", } - for comp, cname in components.items(): + for chan, cname in channels.items(): data[cname] = [] uncertainties[f"{cname} uncertainties"] = [] - for ii, _ in enumerate(survey.channels): - data_entity = geoh5.get_entity(f"Iteration_0_{comp}_[{ii}]")[0].copy( - parent=survey + for ii, _ in enumerate(components.survey.channels): + data_entity = geoh5.get_entity(f"Iteration_0_{chan}_[{ii}]")[0].copy( + parent=components.survey ) data[cname].append(data_entity) - uncert = survey.add_data( + uncert = components.survey.add_data( { - f"uncertainty_{comp}_[{ii}]": { + f"uncertainty_{chan}_[{ii}]": { "values": np.ones_like(data_entity.values) * (np.median(np.abs(data_entity.values))) } @@ -149,16 +148,16 @@ def test_airborne_tem_run(tmp_path: Path, max_iterations=1, pytest=True): ) uncertainties[f"{cname} uncertainties"].append(uncert) - survey.add_components_data(data) - survey.add_components_data(uncertainties) + components.survey.add_components_data(data) + components.survey.add_components_data(uncertainties) data_kwargs = {} - for comp in components: - data_kwargs[f"{comp}_channel"] = survey.fetch_property_group( - name=f"dB{comp}dt" + for chan in channels: + data_kwargs[f"{chan}_channel"] = components.survey.fetch_property_group( + name=f"dB{chan}dt" ) - data_kwargs[f"{comp}_uncertainty"] = survey.fetch_property_group( - name=f"dB{comp}dt uncertainties" + data_kwargs[f"{chan}_uncertainty"] = components.survey.fetch_property_group( + name=f"dB{chan}dt uncertainties" ) orig_dBzdt = geoh5.get_entity("Iteration_0_z_[0]")[0].values @@ -166,9 +165,9 @@ def test_airborne_tem_run(tmp_path: Path, max_iterations=1, pytest=True): # Run the inverse params = TDEMInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, starting_model=1e-3, reference_model=1e-3, chi_factor=1.0, diff --git a/tests/run_tests/driver_dc_2d_test.py b/tests/run_tests/driver_dc_2d_test.py index 66a15975..2a4dc370 100644 --- a/tests/run_tests/driver_dc_2d_test.py +++ b/tests/run_tests/driver_dc_2d_test.py @@ -14,7 +14,6 @@ from pathlib import Path import numpy as np -from geoh5py.data import ReferencedData from geoh5py.workspace import Workspace from simpeg_drivers.electricals.direct_current.two_dimensions.driver import ( @@ -26,12 +25,22 @@ DC2DInversionOptions, ) from simpeg_drivers.options import ( - ActiveCellsOptions, DrapeModelOptions, LineSelectionOptions, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -47,36 +56,34 @@ def test_dc_2d_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.01, - anomaly=10, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - inversion_type="direct current 2d", - drape_height=0.0, - flatten=False, - ) - line_selection = LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0], - line_id=101, - ) - params = DC2DForwardOptions.build( - geoh5=geoh5, - data_object=survey, - line_selection=line_selection, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - horizontal_padding=100.0, - vertical_padding=100.0, - expansion_factor=1.1, - ), - starting_model=model, - topography_object=topography, + opts = SyntheticsComponentsOptions( + method="direct current 2d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=0.01, anomaly=10), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + + line_selection = LineSelectionOptions( + line_object=components.survey.get_data("line_ids")[0], + line_id=101, + ) + params = DC2DForwardOptions.build( + geoh5=geoh5, + data_object=components.survey, + line_selection=line_selection, + drape_model=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=100.0, + horizontal_padding=100.0, + vertical_padding=100.0, + expansion_factor=1.1, + ), + starting_model=components.model, + topography_object=components.topography, + ) fwr_driver = DC2DForwardDriver(params) fwr_driver.run() @@ -88,7 +95,7 @@ def test_dc_2d_run(tmp_path: Path, max_iterations=1, pytest=True): with Workspace(workpath) as geoh5: potential = geoh5.get_entity("Iteration_0_dc")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5) # Run the inverse params = DC2DInversionOptions.build( @@ -101,7 +108,7 @@ def test_dc_2d_run(tmp_path: Path, max_iterations=1, pytest=True): vertical_padding=100.0, expansion_factor=1.1, ), - topography_object=topography, + topography_object=components.topography, line_selection=LineSelectionOptions( line_object=geoh5.get_entity("line_ids")[0], line_id=101, diff --git a/tests/run_tests/driver_dc_b2d_rotated_gradients_test.py b/tests/run_tests/driver_dc_b2d_rotated_gradients_test.py index 9484e2df..259d54ea 100644 --- a/tests/run_tests/driver_dc_b2d_rotated_gradients_test.py +++ b/tests/run_tests/driver_dc_b2d_rotated_gradients_test.py @@ -31,12 +31,22 @@ FileControlOptions, ) from simpeg_drivers.options import ( - ActiveCellsOptions, DrapeModelOptions, LineSelectionOptions, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -48,45 +58,43 @@ def test_dc_rotated_p3d_fwr_run( tmp_path: Path, n_electrodes=10, n_lines=3, refinement=(4, 6) ): - plate_model = PlateModel( - strike_length=1000.0, - dip_length=150.0, - width=20.0, - origin=(0.0, 0.0, -50), - direction=90, - dip=45, - ) - - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model=plate_model, - background=0.01, - anomaly=10.0, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - inversion_type="direct current pseudo 3d", - drape_height=0.0, - flatten=False, - ) - params = DCBatch2DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), - topography_object=topography, - data_object=survey, - starting_model=model, - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0] + opts = SyntheticsComponentsOptions( + method="direct current pseudo 3d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions( + background=0.01, + anomaly=10.0, + plate=PlateModel( + strike_length=1000.0, + dip_length=150.0, + width=20.0, + origin=(0.0, 0.0, -50), + direction=90, + dip=45, + ), ), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = DCBatch2DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + drape_model=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=100.0, + expansion_factor=1.1, + horizontal_padding=1000.0, + vertical_padding=1000.0, + ), + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + line_selection=LineSelectionOptions( + line_object=components.survey.get_data("line_ids")[0] + ), + ) fwr_driver = DCBatch2DForwardDriver(params) fwr_driver.run() @@ -103,29 +111,29 @@ def test_dc_rotated_gradient_p3d_run( ) with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) potential = geoh5.get_entity("Iteration_0_dc")[0] - out_group = geoh5.get_entity("Line 1")[0].parent - mesh = out_group.get_entity("mesh")[0] # Finds the octree mesh - topography = geoh5.get_entity("topography")[0] # Create property group with orientation - dip = np.ones(mesh.n_cells) * 45 - azimuth = np.ones(mesh.n_cells) * 90 + dip = np.ones(components.mesh.n_cells) * 45 + azimuth = np.ones(components.mesh.n_cells) * 90 - data_list = mesh.add_data( + data_list = components.mesh.add_data( { "azimuth": {"values": azimuth}, "dip": {"values": dip}, } ) pg = PropertyGroup( - mesh, properties=data_list, property_group_type="Dip direction & dip" + components.mesh, + properties=data_list, + property_group_type="Dip direction & dip", ) # Run the inverse params = DCBatch2DInversionOptions.build( geoh5=geoh5, - mesh=mesh, + mesh=components.mesh, drape_model=DrapeModelOptions( u_cell_size=5.0, v_cell_size=5.0, @@ -134,7 +142,7 @@ def test_dc_rotated_gradient_p3d_run( horizontal_padding=1000.0, vertical_padding=1000.0, ), - topography_object=topography, + topography_object=components.topography, data_object=potential.parent, gradient_rotation=pg, potential_channel=potential, diff --git a/tests/run_tests/driver_dc_b2d_test.py b/tests/run_tests/driver_dc_b2d_test.py index d6bba018..9e7ffe83 100644 --- a/tests/run_tests/driver_dc_b2d_test.py +++ b/tests/run_tests/driver_dc_b2d_test.py @@ -29,12 +29,22 @@ FileControlOptions, ) from simpeg_drivers.options import ( - ActiveCellsOptions, DrapeModelOptions, LineSelectionOptions, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -50,35 +60,32 @@ def test_dc_p3d_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.01, - anomaly=10.0, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - inversion_type="direct current pseudo 3d", - drape_height=0.0, - flatten=False, - ) - params = DCBatch2DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=1000.0, - vertical_padding=1000.0, - ), - topography_object=topography, - data_object=survey, - starting_model=model, - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0] - ), + opts = SyntheticsComponentsOptions( + method="direct current pseudo 3d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=0.01, anomaly=10.0), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5=geoh5, options=opts) + params = DCBatch2DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + drape_model=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=100.0, + expansion_factor=1.1, + horizontal_padding=1000.0, + vertical_padding=1000.0, + ), + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + line_selection=LineSelectionOptions( + line_object=components.survey.get_data("line_ids")[0] + ), + ) fwr_driver = DCBatch2DForwardDriver(params) fwr_driver.run() @@ -93,15 +100,13 @@ def test_dc_p3d_run( workpath = tmp_path.parent / "test_dc_p3d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) potential = geoh5.get_entity("Iteration_0_dc")[0] - out_group = geoh5.get_entity("Line 1")[0].parent - mesh = out_group.get_entity("mesh")[0] # Finds the octree mesh - topography = geoh5.get_entity("topography")[0] # Run the inverse params = DCBatch2DInversionOptions.build( geoh5=geoh5, - mesh=mesh, + mesh=components.mesh, drape_model=DrapeModelOptions( u_cell_size=5.0, v_cell_size=5.0, @@ -110,7 +115,7 @@ def test_dc_p3d_run( horizontal_padding=1000.0, vertical_padding=1000.0, ), - topography_object=topography, + topography_object=components.topography, data_object=potential.parent, potential_channel=potential, potential_uncertainty=1e-3, diff --git a/tests/run_tests/driver_dc_test.py b/tests/run_tests/driver_dc_test.py index dc891abd..a4ec28a3 100644 --- a/tests/run_tests/driver_dc_test.py +++ b/tests/run_tests/driver_dc_test.py @@ -23,9 +23,19 @@ DC3DForwardOptions, DC3DInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -41,37 +51,38 @@ def test_dc_3d_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.01, - anomaly=10, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - drape_height=0.0, - inversion_type="direct current 3d", - flatten=False, + opts = SyntheticsComponentsOptions( + method="direct current 3d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=0.01, anomaly=10.0), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + + # Randomly flip order of receivers + old = np.random.randint(0, components.survey.cells.shape[0], n_electrodes) + indices = np.ones(components.survey.cells.shape[0], dtype=bool) + indices[old] = False + + tx_id = np.r_[ + components.survey.ab_cell_id.values[indices], + components.survey.ab_cell_id.values[~indices], + ] + cells = np.vstack( + [components.survey.cells[indices, :], components.survey.cells[~indices, :]] + ) - # Randomly flip order of receivers - old = np.random.randint(0, survey.cells.shape[0], n_electrodes) - indices = np.ones(survey.cells.shape[0], dtype=bool) - indices[old] = False - - tx_id = np.r_[survey.ab_cell_id.values[indices], survey.ab_cell_id.values[~indices]] - cells = np.vstack([survey.cells[indices, :], survey.cells[~indices, :]]) - - with survey.workspace.open(): - survey.ab_cell_id = tx_id - survey.cells = cells + components.survey.ab_cell_id = tx_id + components.survey.cells = cells - params = DC3DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - ) + params = DC3DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + ) fwr_driver = DC3DForwardDriver(params) fwr_driver.run() @@ -87,15 +98,14 @@ def test_dc_3d_run( workpath = tmp_path.parent / "test_dc_3d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) potential = geoh5.get_entity("Iteration_0_dc")[0] - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] # Run the inverse params = DC3DInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, + mesh=components.mesh, + topography_object=components.topography, data_object=potential.parent, starting_model=1e-2, reference_model=1e-2, @@ -140,24 +150,21 @@ def test_dc_single_line_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.01, - anomaly=10, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - drape_height=0.0, - inversion_type="direct current 3d", - flatten=False, - ) - params = DC3DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, + opts = SyntheticsComponentsOptions( + method="direct current 3d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=0.01, anomaly=10.0), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = DC3DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + ) fwr_driver = DC3DForwardDriver(params) assert np.all(fwr_driver.window.window["size"] > 0) diff --git a/tests/run_tests/driver_fem_test.py b/tests/run_tests/driver_fem_test.py index 695e6ecb..f4fb0918 100644 --- a/tests/run_tests/driver_fem_test.py +++ b/tests/run_tests/driver_fem_test.py @@ -28,9 +28,19 @@ FDEMForwardOptions, FDEMInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -41,29 +51,26 @@ def test_fem_name_change(tmp_path, caplog): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=1e-3, - anomaly=1.0, - n_electrodes=2, - n_lines=2, - refinement=(2,), - drape_height=15.0, - padding_distance=400, - inversion_type="fdem", + opts = SyntheticsComponentsOptions( + method="fdem", + survey=SurveyOptions(n_stations=2, n_lines=2, drape=15.0), + mesh=MeshOptions(refinement=(2,), padding_distance=400.0), + model=ModelOptions(background=1e-3), ) - with caplog.at_level(logging.WARNING): - FDEMForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - z_real_channel_bool=True, - z_imag_channel_bool=True, - inversion_type="fem", - ) - assert "fdem" in caplog.text + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + with caplog.at_level(logging.WARNING): + FDEMForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + z_real_channel_bool=True, + z_imag_channel_bool=True, + inversion_type="fem", + ) + assert "fdem" in caplog.text def test_fem_fwr_run( @@ -73,35 +80,38 @@ def test_fem_fwr_run( cell_size=(20.0, 20.0, 20.0), ): # Run the forward - plate_model = PlateModel( - strike_length=40.0, - dip_length=40.0, - width=40.0, - origin=(0.0, 0.0, -50.0), - ) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model=plate_model, - background=1e-3, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - drape_height=15.0, - cell_size=cell_size, - padding_distance=400, - inversion_type="fdem", - flatten=True, - ) - params = FDEMForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - z_real_channel_bool=True, - z_imag_channel_bool=True, + opts = SyntheticsComponentsOptions( + method="fdem", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + drape=15.0, + topography=lambda x, y: np.zeros(x.shape), + ), + mesh=MeshOptions( + cell_size=cell_size, refinement=refinement, padding_distance=400.0 + ), + model=ModelOptions( + background=1e-3, + plate=PlateModel( + strike_length=40.0, + dip_length=40.0, + width=40.0, + origin=(0.0, 0.0, -50.0), + ), + ), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = FDEMForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + z_real_channel_bool=True, + z_imag_channel_bool=True, + ) fwr_driver = FDEMForwardDriver(params) fwr_driver.run() @@ -113,32 +123,26 @@ def test_fem_run(tmp_path: Path, max_iterations=1, pytest=True): workpath = tmp_path.parent / "test_fem_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: - survey = next( - child - for child in geoh5.get_entity("Airborne_rx") - if not isinstance(child.parent, SimPEGGroup) - ) - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5) data = {} uncertainties = {} - components = { + channels = { "z_real": "z_real", "z_imag": "z_imag", } - for comp, cname in components.items(): + for chan, cname in channels.items(): data[cname] = [] uncertainties[f"{cname} uncertainties"] = [] - for ind, freq in enumerate(survey.channels): - data_entity = geoh5.get_entity(f"Iteration_0_{comp}_[{ind}]")[0].copy( - parent=survey + for ind, freq in enumerate(components.survey.channels): + data_entity = geoh5.get_entity(f"Iteration_0_{chan}_[{ind}]")[0].copy( + parent=components.survey ) data[cname].append(data_entity) abs_val = np.abs(data_entity.values) - uncert = survey.add_data( + uncert = components.survey.add_data( { - f"uncertainty_{comp}_[{ind}]": { + f"uncertainty_{chan}_[{ind}]": { "values": np.ones_like(abs_val) * freq / 200.0 # * 2**(np.abs(ind-1)) @@ -146,27 +150,27 @@ def test_fem_run(tmp_path: Path, max_iterations=1, pytest=True): } ) uncertainties[f"{cname} uncertainties"].append( - uncert.copy(parent=survey) + uncert.copy(parent=components.survey) ) - data_groups = survey.add_components_data(data) - uncert_groups = survey.add_components_data(uncertainties) + data_groups = components.survey.add_components_data(data) + uncert_groups = components.survey.add_components_data(uncertainties) data_kwargs = {} - for comp, data_group, uncert_group in zip( - components, data_groups, uncert_groups, strict=True + for chan, data_group, uncert_group in zip( + channels, data_groups, uncert_groups, strict=True ): - data_kwargs[f"{comp}_channel"] = data_group - data_kwargs[f"{comp}_uncertainty"] = uncert_group + data_kwargs[f"{chan}_channel"] = data_group + data_kwargs[f"{chan}_uncertainty"] = uncert_group orig_z_real_1 = geoh5.get_entity("Iteration_0_z_real_[0]")[0].values # Run the inverse params = FDEMInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, starting_model=1e-3, reference_model=1e-3, alpha_s=0.0, diff --git a/tests/run_tests/driver_grav_test.py b/tests/run_tests/driver_grav_test.py index 9c120092..8c7917b7 100644 --- a/tests/run_tests/driver_grav_test.py +++ b/tests/run_tests/driver_grav_test.py @@ -14,12 +14,11 @@ from unittest.mock import patch import numpy as np -from geoapps_utils.modelling.plates import PlateModel from geoapps_utils.utils.importing import GeoAppsError +from geoapps_utils.utils.locations import gaussian from geoh5py.workspace import Workspace from pytest import raises -from simpeg_drivers.options import ActiveCellsOptions, ModelOptions from simpeg_drivers.potential_fields import ( GravityForwardOptions, GravityInversionOptions, @@ -28,8 +27,17 @@ GravityForwardDriver, GravityInversionDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -42,34 +50,30 @@ def test_gravity_fwr_run( n_grid_points=2, refinement=(2,), ): - plate_model = PlateModel( - strike_length=40.0, - dip_length=40.0, - width=40.0, - origin=(0.0, 0.0, 10.0), - ) + filepath = Path(tmp_path) / "inversion_test.ui.geoh5" + with Workspace.create(filepath) as geoh5: + # Run the forward + components = SyntheticsComponents( + geoh5=geoh5, + options=SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0 + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.75), + ), + ) - # Run the forward - geoh5, mesh, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model, - background=0.0, - anomaly=0.75, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - inversion_type="gravity", - flatten=False, - ) + params = GravityForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + gz_channel_bool=True, + ) - params = GravityForwardOptions.build( - geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, - starting_model=model, - gz_channel_bool=True, - ) fwr_driver = GravityForwardDriver(params) fwr_driver.run() @@ -80,15 +84,14 @@ def test_array_too_large_run( workpath = tmp_path.parent / "test_gravity_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) gz = geoh5.get_entity("Iteration_0_gz")[0] - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] # Run the inverse params = GravityInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, + mesh=components.mesh, + topography_object=components.topography, data_object=gz.parent, gz_channel=gz, gz_uncertainty=1e-4, @@ -114,25 +117,23 @@ def test_gravity_run( workpath = tmp_path.parent / "test_gravity_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: - group = geoh5.get_entity("Gravity Forward")[0] gz = geoh5.get_entity("Iteration_0_gz")[0] orig_gz = gz.values.copy() - mesh = group.get_entity("mesh")[0] - model = mesh.get_entity("starting_model")[0] + components = SyntheticsComponents(geoh5) - inds = (mesh.centroids[:, 0] > -35) & (mesh.centroids[:, 0] < 35) - norms = np.ones(mesh.n_cells) * 2 + inds = (components.mesh.centroids[:, 0] > -35) & ( + components.mesh.centroids[:, 0] < 35 + ) + norms = np.ones(components.mesh.n_cells) * 2 norms[inds] = 0 - gradient_norms = mesh.add_data({"norms": {"values": norms}}) + gradient_norms = components.mesh.add_data({"norms": {"values": norms}}) # Test mesh UBC ordered - ind = np.argsort(mesh.octree_cells, order=["K", "J", "I"]) - mesh.octree_cells = mesh.octree_cells[ind] - model.values = model.values[ind] + ind = np.argsort(components.mesh.octree_cells, order=["K", "J", "I"]) + components.mesh.octree_cells = components.mesh.octree_cells[ind] + components.model.values = components.model.values[ind] gradient_norms.values = gradient_norms.values[ind] - topography = geoh5.get_entity("topography")[0] - # Turn some values to nan values = gz.values.copy() values[0] = np.nan @@ -141,7 +142,7 @@ def test_gravity_run( # Run the inverse params = GravityInversionOptions.build( geoh5=geoh5, - mesh=mesh, + mesh=components.mesh, data_object=gz.parent, s_norm=0.0, x_norm=gradient_norms, @@ -154,7 +155,7 @@ def test_gravity_run( initial_beta_ratio=1e-2, percentile=100, starting_model=1e-4, - topography_object=topography, + topography_object=components.topography, reference_model=0.0, save_sensitivities=True, ) diff --git a/tests/run_tests/driver_ground_tem_test.py b/tests/run_tests/driver_ground_tem_test.py index 33783acc..9eedbf08 100644 --- a/tests/run_tests/driver_ground_tem_test.py +++ b/tests/run_tests/driver_ground_tem_test.py @@ -26,9 +26,19 @@ TDEMForwardOptions, TDEMInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) logger = getLogger(__name__) @@ -48,37 +58,39 @@ def test_tiling_ground_tem( **_, ): # Run the forward - plate_model = PlateModel( - strike_length=40.0, - dip_length=40.0, - width=40.0, - origin=(0.0, 0.0, -50.0), - ) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model=plate_model, - background=0.001, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - inversion_type="ground tdem", - drape_height=5.0, - padding_distance=1000.0, - flatten=True, - ) - - params = TDEMForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - x_channel_bool=True, - y_channel_bool=True, - z_channel_bool=True, - tile_spatial=4, + opts = SyntheticsComponentsOptions( + method="ground tdem", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + drape=5.0, + topography=lambda x, y: np.zeros(x.shape), + name="ground_tdem_survey", + ), + mesh=MeshOptions(refinement=refinement, padding_distance=1000.0), + model=ModelOptions( + background=0.001, + plate=PlateModel( + strike_length=40.0, + dip_length=40.0, + width=40.0, + origin=(0.0, 0.0, -50.0), + ), + ), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = TDEMForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + x_channel_bool=True, + y_channel_bool=True, + z_channel_bool=True, + tile_spatial=4, + ) fwr_driver = TDEMForwardDriver(params) tiles = fwr_driver.get_tiles() @@ -86,7 +98,7 @@ def test_tiling_ground_tem( assert len(tiles) == 4 for tile in tiles: - assert len(np.unique(survey.tx_id_property.values[tile])) == 1 + assert len(np.unique(components.survey.tx_id_property.values[tile])) == 1 fwr_driver.run() @@ -102,43 +114,46 @@ def test_ground_tem_fwr_run( if pytest: caplog.set_level(INFO) # Run the forward - plate_model = PlateModel( - strike_length=40.0, - dip_length=40.0, - width=40.0, - origin=(0.0, 0.0, -50.0), - ) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model=plate_model, - background=0.001, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - inversion_type="ground tdem", - drape_height=5.0, - cell_size=cell_size, - padding_distance=1000.0, - flatten=True, - ) - params = TDEMForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - x_channel_bool=True, - y_channel_bool=True, - z_channel_bool=True, - solver_type="Mumps", + opts = SyntheticsComponentsOptions( + method="ground tdem", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + drape=5.0, + topography=lambda x, y: np.zeros(x.shape), + ), + mesh=MeshOptions( + cell_size=cell_size, refinement=refinement, padding_distance=1000.0 + ), + model=ModelOptions( + background=0.001, + plate=PlateModel( + strike_length=40.0, + dip_length=40.0, + width=40.0, + origin=(0.0, 0.0, -50.0), + ), + ), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + components.survey.transmitters.remove_cells([15]) + params = TDEMForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + x_channel_bool=True, + y_channel_bool=True, + z_channel_bool=True, + solver_type="Mumps", + ) fwr_driver = TDEMForwardDriver(params) - with survey.workspace.open(): - survey.transmitters.remove_cells([15]) - survey.tx_id_property.name = "tx_id" + with components.survey.workspace.open(): + components.survey.tx_id_property.name = "tx_id" assert fwr_driver.inversion_data.survey.source_list[0].n_segments == 16 if pytest: @@ -161,29 +176,25 @@ def test_ground_tem_run(tmp_path: Path, max_iterations=1, pytest=True): ) with Workspace(workpath) as geoh5: - simpeg_group = geoh5.get_entity("Time-domain EM (TEM) Forward")[0] - survey = simpeg_group.get_entity("Ground TEM Rx")[0] - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] - + components = SyntheticsComponents(geoh5) data = {} uncertainties = {} - components = { + channels = { "z": "dBzdt", } - for comp, cname in components.items(): + for chan, cname in channels.items(): data[cname] = [] uncertainties[f"{cname} uncertainties"] = [] - for ii, _ in enumerate(survey.channels): - data_entity = geoh5.get_entity(f"Iteration_0_{comp}_[{ii}]")[0].copy( - parent=survey + for ii, _ in enumerate(components.survey.channels): + data_entity = geoh5.get_entity(f"Iteration_0_{chan}_[{ii}]")[0].copy( + parent=components.survey ) data[cname].append(data_entity) - uncert = survey.add_data( + uncert = components.survey.add_data( { - f"uncertainty_{comp}_[{ii}]": { + f"uncertainty_{chan}_[{ii}]": { "values": np.ones_like(data_entity.values) * np.median(np.abs(data_entity.values)) / 2.0 @@ -192,16 +203,16 @@ def test_ground_tem_run(tmp_path: Path, max_iterations=1, pytest=True): ) uncertainties[f"{cname} uncertainties"].append(uncert) - survey.add_components_data(data) - survey.add_components_data(uncertainties) + components.survey.add_components_data(data) + components.survey.add_components_data(uncertainties) data_kwargs = {} - for comp in components: - data_kwargs[f"{comp}_channel"] = survey.fetch_property_group( - name=f"Iteration_0_{comp}" + for chan in channels: + data_kwargs[f"{chan}_channel"] = components.survey.fetch_property_group( + name=f"Iteration_0_{chan}" ) - data_kwargs[f"{comp}_uncertainty"] = survey.fetch_property_group( - name=f"dB{comp}dt uncertainties" + data_kwargs[f"{chan}_uncertainty"] = components.survey.fetch_property_group( + name=f"dB{chan}dt uncertainties" ) orig_dBzdt = geoh5.get_entity("Iteration_0_z_[0]")[0].values @@ -209,9 +220,9 @@ def test_ground_tem_run(tmp_path: Path, max_iterations=1, pytest=True): # Run the inverse params = TDEMInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, starting_model=1e-3, reference_model=1e-3, chi_factor=0.1, @@ -239,7 +250,7 @@ def test_ground_tem_run(tmp_path: Path, max_iterations=1, pytest=True): output = get_inversion_output( driver.params.geoh5.h5file, driver.params.out_group.uid ) - assert driver.inversion_data.entity.tx_id_property.name == "tx_id" + assert driver.inversion_data.entity.tx_id_property.name == "Transmitter ID" output["data"] = orig_dBzdt if pytest: check_target(output, target_run) diff --git a/tests/run_tests/driver_ip_2d_test.py b/tests/run_tests/driver_ip_2d_test.py index 5c2248f5..6ca51ffb 100644 --- a/tests/run_tests/driver_ip_2d_test.py +++ b/tests/run_tests/driver_ip_2d_test.py @@ -23,9 +23,20 @@ IP2DForwardDriver, IP2DInversionDriver, ) -from simpeg_drivers.options import ActiveCellsOptions, LineSelectionOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.options import LineSelectionOptions +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -41,30 +52,27 @@ def test_ip_2d_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=1e-6, - anomaly=1e-1, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - inversion_type="induced polarization 2d", - flatten=False, - drape_height=0.0, - ) - params = IP2DForwardOptions.build( - geoh5=geoh5, - data_object=survey, - mesh=model.parent, - topography_object=topography, - starting_model=model, - conductivity_model=1e2, - model_type="Resistivity (Ohm-m)", - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0], - line_id=101, - ), + opts = SyntheticsComponentsOptions( + method="induced polarization 2d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=1e-6, anomaly=1e-1), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = IP2DForwardOptions.build( + geoh5=geoh5, + data_object=components.survey, + mesh=components.mesh, + topography_object=components.topography, + starting_model=components.model, + conductivity_model=1e2, + model_type="Resistivity (Ohm-m)", + line_selection=LineSelectionOptions( + line_object=geoh5.get_entity("line_ids")[0], + line_id=101, + ), + ) fwr_driver = IP2DForwardDriver(params) fwr_driver.run() @@ -80,15 +88,14 @@ def test_ip_2d_run( workpath = tmp_path.parent / "test_ip_2d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) chargeability = geoh5.get_entity("Iteration_0_ip")[0] - mesh = geoh5.get_entity("Models")[0] - topography = geoh5.get_entity("topography")[0] # Run the inverse params = IP2DInversionOptions.build( geoh5=geoh5, - mesh=mesh, - topography_object=topography, + mesh=components.mesh, + topography_object=components.topography, data_object=chargeability.parent, chargeability_channel=chargeability, chargeability_uncertainty=2e-4, diff --git a/tests/run_tests/driver_ip_b2d_test.py b/tests/run_tests/driver_ip_b2d_test.py index 7fd3d686..52082182 100644 --- a/tests/run_tests/driver_ip_b2d_test.py +++ b/tests/run_tests/driver_ip_b2d_test.py @@ -32,8 +32,19 @@ DrapeModelOptions, LineSelectionOptions, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -49,39 +60,36 @@ def test_ip_p3d_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=1e-6, - anomaly=1e-1, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - inversion_type="induced polarization pseudo 3d", - drape_height=0.0, - flatten=False, + opts = SyntheticsComponentsOptions( + method="induced polarization pseudo 3d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=1e-6, anomaly=1e-1), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - params = IPBatch2DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - drape_model=DrapeModelOptions( - u_cell_size=5.0, - v_cell_size=5.0, - depth_core=100.0, - expansion_factor=1.1, - horizontal_padding=100.0, - vertical_padding=100.0, - ), - active_cells=ActiveCellsOptions( - topography_object=topography, - ), - data_object=survey, - conductivity_model=1e-2, - starting_model=model, - line_selection=LineSelectionOptions( - line_object=geoh5.get_entity("line_ids")[0] - ), - ) + params = IPBatch2DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + drape_model=DrapeModelOptions( + u_cell_size=5.0, + v_cell_size=5.0, + depth_core=100.0, + expansion_factor=1.1, + horizontal_padding=100.0, + vertical_padding=100.0, + ), + active_cells=ActiveCellsOptions( + topography_object=components.topography, + ), + data_object=components.survey, + conductivity_model=1e-2, + starting_model=components.model, + line_selection=LineSelectionOptions( + line_object=geoh5.get_entity("line_ids")[0] + ), + ) fwr_driver = IPBatch2DForwardDriver(params) fwr_driver.run() @@ -97,15 +105,13 @@ def test_ip_p3d_run( workpath = tmp_path.parent / "test_ip_p3d_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: + components = SyntheticsComponents(geoh5) chargeability = geoh5.get_entity("Iteration_0_ip")[0] - out_group = geoh5.get_entity("Line 1")[0].parent - mesh = out_group.get_entity("mesh")[0] # Finds the octree mesh - topography = geoh5.get_entity("topography")[0] # Run the inverse params = IPBatch2DInversionOptions.build( geoh5=geoh5, - mesh=mesh, + mesh=components.mesh, drape_model=DrapeModelOptions( u_cell_size=5.0, v_cell_size=5.0, @@ -114,7 +120,7 @@ def test_ip_p3d_run( horizontal_padding=1000.0, vertical_padding=1000.0, ), - topography_object=topography, + topography_object=components.topography, data_object=chargeability.parent, chargeability_channel=chargeability, chargeability_uncertainty=2e-4, diff --git a/tests/run_tests/driver_ip_test.py b/tests/run_tests/driver_ip_test.py index b17377dd..b3d444d8 100644 --- a/tests/run_tests/driver_ip_test.py +++ b/tests/run_tests/driver_ip_test.py @@ -22,9 +22,19 @@ IP3DForwardDriver, IP3DInversionDriver, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -40,25 +50,22 @@ def test_ip_3d_fwr_run( refinement=(4, 6), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=1e-6, - anomaly=1e-1, - n_electrodes=n_electrodes, - n_lines=n_lines, - refinement=refinement, - drape_height=0.0, - inversion_type="induced polarization 3d", - flatten=False, - ) - params = IP3DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - conductivity_model=1e-2, + opts = SyntheticsComponentsOptions( + method="induced polarization 3d", + survey=SurveyOptions(n_stations=n_electrodes, n_lines=n_lines), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(background=1e-6, anomaly=1e-1), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = IP3DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + conductivity_model=1e-2, + ) fwr_driver = IP3DForwardDriver(params) fwr_driver.run() diff --git a/tests/run_tests/driver_joint_cross_gradient_test.py b/tests/run_tests/driver_joint_cross_gradient_test.py index e343e4cc..c7d10530 100644 --- a/tests/run_tests/driver_joint_cross_gradient_test.py +++ b/tests/run_tests/driver_joint_cross_gradient_test.py @@ -25,7 +25,6 @@ ) from simpeg_drivers.joint.joint_cross_gradient import JointCrossGradientOptions from simpeg_drivers.joint.joint_cross_gradient.driver import JointCrossGradientDriver -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import ( GravityForwardOptions, GravityInversionOptions, @@ -40,14 +39,29 @@ MVIForwardDriver, MVIInversionDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + ActiveCellsOptions as SyntheticsActiveCellsOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. # Move this file out of the test directory and run. target_run = {"data_norm": 53.29585552088844, "phi_d": 7970, "phi_m": 26.7} +INDUCING_FIELD = (50000.0, 90.0, 0.0) def test_joint_cross_gradient_fwr_run( @@ -57,74 +71,75 @@ def test_joint_cross_gradient_fwr_run( refinement=(2,), ): # Create local problem A - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - drape_height=15.0, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - ) - params = GravityForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=15.0, name="survey A" + ), + mesh=MeshOptions(refinement=refinement, name="mesh A"), + model=ModelOptions(anomaly=0.75, name="model A"), + active=SyntheticsActiveCellsOptions(name="active A"), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = GravityForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + ) fwr_driver_a = GravityForwardDriver(params) with geoh5.open(): - _, _, model, survey, _ = setup_inversion_workspace( - tmp_path, + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + drape=15.0, + name="survey B", + ), + mesh=MeshOptions(refinement=refinement, name="mesh B"), + model=ModelOptions(anomaly=0.05, name="model B"), + active=SyntheticsActiveCellsOptions(name="active B"), + ) + components = SyntheticsComponents(geoh5, options=opts) + params = MVIForwardOptions.build( geoh5=geoh5, - background=0.0, - anomaly=0.05, - drape_height=15.0, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - flatten=False, + mesh=components.mesh, + topography_object=components.topography, + inducing_field_strength=INDUCING_FIELD[0], + inducing_field_inclination=INDUCING_FIELD[1], + inducing_field_declination=INDUCING_FIELD[2], + data_object=components.survey, + starting_model=components.model, ) - inducing_field = (50000.0, 90.0, 0.0) - params = MVIForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - inducing_field_strength=inducing_field[0], - inducing_field_inclination=inducing_field[1], - inducing_field_declination=inducing_field[2], - data_object=survey, - starting_model=model, - ) fwr_driver_b = MVIForwardDriver(params) with geoh5.open(): - _, _, model, survey, _ = setup_inversion_workspace( - tmp_path, - geoh5=geoh5, - background=0.01, - anomaly=10, - n_electrodes=n_grid_points, - n_lines=n_lines, - refinement=refinement, - drape_height=0.0, - inversion_type="direct current 3d", - flatten=False, + opts = SyntheticsComponentsOptions( + method="direct current 3d", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_lines, name="survey C" + ), + mesh=MeshOptions(refinement=refinement, name="mesh C"), + model=ModelOptions(background=0.01, anomaly=10, name="model C"), + active=SyntheticsActiveCellsOptions(name="active C"), ) + components = SyntheticsComponents(geoh5, options=opts) - params = DC3DForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - ) + params = DC3DForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + ) fwr_driver_c = DC3DForwardDriver(params) with geoh5.open(): - fwr_driver_c.inversion_data.entity.name = "survey" + fwr_driver_c.inversion_data.entity.name = "survey C" # Force co-location of meshes for driver in [fwr_driver_b, fwr_driver_c]: @@ -159,30 +174,24 @@ def test_joint_cross_gradient_inv_run( drivers = [] orig_data = [] - for group_name in [ - "Gravity Forward", - "Magnetic Vector Forward", - "Direct Current 3D Forward", - ]: - group = geoh5.get_entity(group_name)[0] - - if not isinstance(group, SimPEGGroup): - continue - - mesh = group.get_entity("mesh")[0] - survey = group.get_entity("survey")[0] - - data = None - for child in survey.children: - if isinstance(child, FloatData): - data = child - - if data is None: - raise ValueError("No data found in survey") + for suffix in "ABC": + components = SyntheticsComponents( + geoh5=geoh5, + options=SyntheticsComponentsOptions( + method="joint", + survey=SurveyOptions(name=f"survey {suffix}"), + mesh=MeshOptions(name=f"mesh {suffix}"), + model=ModelOptions(name=f"model {suffix}"), + active=SyntheticsActiveCellsOptions(name=f"active {suffix}"), + ), + ) + mesh = components.mesh + survey = components.survey + data = next(k for k in survey.children if "Iteration_0" in k.name) orig_data.append(data.values) - if group.options["inversion_type"] == "gravity": + if suffix == "A": params = GravityInversionOptions.build( geoh5=geoh5, mesh=mesh, @@ -195,7 +204,7 @@ def test_joint_cross_gradient_inv_run( reference_model=0.0, ) drivers.append(GravityInversionDriver(params)) - elif group.options["inversion_type"] == "direct current 3d": + elif suffix == "C": params = DC3DInversionOptions.build( geoh5=geoh5, mesh=mesh, @@ -218,15 +227,9 @@ def test_joint_cross_gradient_inv_run( mesh=mesh, alpha_s=1.0, topography_object=topography, - inducing_field_strength=group.options["inducing_field_strength"][ - "value" - ], - inducing_field_inclination=group.options[ - "inducing_field_inclination" - ]["value"], - inducing_field_declination=group.options[ - "inducing_field_declination" - ]["value"], + inducing_field_strength=INDUCING_FIELD[0], + inducing_field_inclination=INDUCING_FIELD[1], + inducing_field_declination=INDUCING_FIELD[2], data_object=survey, starting_model=1e-4, reference_model=0.0, diff --git a/tests/run_tests/driver_joint_pgi_homogeneous_test.py b/tests/run_tests/driver_joint_pgi_homogeneous_test.py index 6ccbefcc..24424995 100644 --- a/tests/run_tests/driver_joint_pgi_homogeneous_test.py +++ b/tests/run_tests/driver_joint_pgi_homogeneous_test.py @@ -11,10 +11,8 @@ from __future__ import annotations from pathlib import Path -from unittest.mock import patch import numpy as np -from geoapps_utils.utils.importing import GeoAppsError from geoh5py.data import FloatData from geoh5py.groups.property_group import GroupTypeEnum, PropertyGroup from geoh5py.groups.simpeg import SimPEGGroup @@ -38,14 +36,29 @@ from simpeg_drivers.potential_fields.magnetic_vector.driver import ( MVIForwardDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + ActiveCellsOptions as SyntheticsActiveCellsOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. # Move this file out of the test directory and run. target_run = {"data_norm": 390.70695894864303, "phi_d": 2020, "phi_m": 8.98} +INDUCING_FIELD = (50000.0, 90.0, 0.0) def test_homogeneous_fwr_run( @@ -54,56 +67,62 @@ def test_homogeneous_fwr_run( refinement=(2,), ): # Create local problem A - geoh5, mesh, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - drape_height=15.0, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + drape=15.0, + name="survey A", + ), + mesh=MeshOptions(refinement=refinement, name="mesh A"), + model=ModelOptions(anomaly=0.75, name="model A"), + active=SyntheticsActiveCellsOptions(name="active A"), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - # Change half the model - ind = mesh.centroids[:, 0] > 0 - model.values[ind] = 0.05 + # Change half the model + ind = components.mesh.centroids[:, 0] > 0 + components.model.values[ind] = 0.05 - params = GravityForwardOptions.build( - geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=survey, - starting_model=model, - ) + params = GravityForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + ) fwr_driver_a = GravityForwardDriver(params) with geoh5.open(): - _, mesh, model, survey, _ = setup_inversion_workspace( - tmp_path, - geoh5=geoh5, - background=0.0, - anomaly=0.05, - drape_height=15.0, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - flatten=False, + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + drape=15.0, + name="survey B", + ), + mesh=MeshOptions(refinement=refinement, name="mesh B"), + model=ModelOptions(anomaly=0.05, name="model B"), + active=SyntheticsActiveCellsOptions(name="active B"), ) - inducing_field = (50000.0, 90.0, 0.0) - # Change half the model - ind = mesh.centroids[:, 0] > 0 - model.values[ind] = 0.01 + components = SyntheticsComponents(geoh5, options=opts) + # Change half the model + ind = components.mesh.centroids[:, 0] > 0 + components.model.values[ind] = 0.01 - params = MVIForwardOptions.build( - geoh5=geoh5, - mesh=mesh, - topography_object=topography, - inducing_field_strength=inducing_field[0], - inducing_field_inclination=inducing_field[1], - inducing_field_declination=inducing_field[2], - data_object=survey, - starting_model=model, - ) + params = MVIForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + inducing_field_strength=INDUCING_FIELD[0], + inducing_field_inclination=INDUCING_FIELD[1], + inducing_field_declination=INDUCING_FIELD[2], + data_object=components.survey, + starting_model=components.model, + ) fwr_driver_b = MVIForwardDriver(params) fwr_driver_a.run() @@ -128,19 +147,21 @@ def test_homogeneous_run( petrophysics = None gradient_rotation = None mesh = None - for group_name in [ - "Gravity Forward", - "Magnetic Vector Forward", - ]: - group = geoh5.get_entity(group_name)[0] - - if not isinstance(group, SimPEGGroup): - continue - - mesh = group.get_entity("mesh")[0] - survey = group.get_entity("survey")[0] + for suffix in "AB": + components = SyntheticsComponents( + geoh5=geoh5, + options=SyntheticsComponentsOptions( + method="joint", + survey=SurveyOptions(name=f"survey {suffix}"), + mesh=MeshOptions(name=f"mesh {suffix}"), + model=ModelOptions(name=f"model {suffix}"), + active=SyntheticsActiveCellsOptions(name=f"active {suffix}"), + ), + ) + mesh = components.mesh + survey = components.survey - if group_name == "Gravity Forward": + if suffix == "A": global_mesh = mesh.copy(parent=geoh5) model = global_mesh.get_entity("starting_model")[0] @@ -173,19 +194,13 @@ def test_homogeneous_run( parent=global_mesh, ) - data = None - for child in survey.children: - if isinstance(child, FloatData): - data = child - - if data is None: - raise ValueError("No data found in survey") + data = next(k for k in survey.children if "Iteration_0" in k.name) orig_data.append(data.values) ref_model = mesh.get_entity("starting_model")[0].copy(name="ref_model") ref_model.values = ref_model.values / 2.0 - if group.options["inversion_type"] == "gravity": + if suffix == "A": params = GravityInversionOptions.build( geoh5=geoh5, mesh=mesh, @@ -202,15 +217,9 @@ def test_homogeneous_run( geoh5=geoh5, mesh=mesh, topography_object=topography, - inducing_field_strength=group.options["inducing_field_strength"][ - "value" - ], - inducing_field_inclination=group.options[ - "inducing_field_inclination" - ]["value"], - inducing_field_declination=group.options[ - "inducing_field_declination" - ]["value"], + inducing_field_strength=INDUCING_FIELD[0], + inducing_field_inclination=INDUCING_FIELD[1], + inducing_field_declination=INDUCING_FIELD[2], data_object=survey, starting_model=ref_model, reference_model=ref_model, @@ -248,7 +257,7 @@ def test_homogeneous_run( check_target(output, target_run) out_group = run_ws.get_entity(driver.params.out_group.uid)[0] - mesh = out_group.get_entity("mesh")[0] + mesh = out_group.get_entity("mesh A")[0] petro_model = mesh.get_entity("petrophysical_model")[0] assert len(np.unique(petro_model.values)) == 4 diff --git a/tests/run_tests/driver_joint_surveys_test.py b/tests/run_tests/driver_joint_surveys_test.py index 793f5bc1..1aa1fb8b 100644 --- a/tests/run_tests/driver_joint_surveys_test.py +++ b/tests/run_tests/driver_joint_surveys_test.py @@ -22,8 +22,22 @@ GravityInversionOptions, ) from simpeg_drivers.potential_fields.gravity.driver import GravityInversionDriver -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + ActiveCellsOptions as SyntheticsActiveCellsOptions, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -38,21 +52,24 @@ def test_joint_surveys_fwr_run( refinement=(2,), ): # Create local problem A - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - ) - params = GravityForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0, name="survey A" + ), + mesh=MeshOptions(refinement=refinement, name="mesh A"), + model=ModelOptions(anomaly=0.75, name="model A"), + active=SyntheticsActiveCellsOptions(name="active A"), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = GravityForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + ) fwr_driver_a = GravityInversionDriver(params) with fwr_driver_a.out_group.workspace.open(): @@ -60,24 +77,26 @@ def test_joint_surveys_fwr_run( # Create local problem B with geoh5.open(): - _, _, model, survey, _ = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - refinement=[0, 2], - n_electrodes=int(n_grid_points / 2), - n_lines=int(n_grid_points / 2), - flatten=False, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=int(n_grid_points / 2), + n_lines=int(n_grid_points / 2), + drape=10.0, + name="survey B", + ), + mesh=MeshOptions(refinement=(0, 2), name="mesh B"), + model=ModelOptions(anomaly=0.75, name="model B"), + active=SyntheticsActiveCellsOptions(name="active B"), + ) + components = SyntheticsComponents(geoh5, options=opts) + params = GravityForwardOptions.build( geoh5=geoh5, - drape_height=10.0, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, ) - params = GravityForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - ) fwr_driver_b = GravityInversionDriver(params) with fwr_driver_b.out_group.workspace.open(): diff --git a/tests/run_tests/driver_mag_automesh_test.py b/tests/run_tests/driver_mag_automesh_test.py index 1ff96236..6b88e266 100644 --- a/tests/run_tests/driver_mag_automesh_test.py +++ b/tests/run_tests/driver_mag_automesh_test.py @@ -13,20 +13,21 @@ from pathlib import Path import numpy as np -from dask.distributed import LocalCluster, performance_report -from geoh5py.workspace import Workspace +from geoh5py import Workspace -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import ( MagneticForwardOptions, - MagneticInversionOptions, ) from simpeg_drivers.potential_fields.magnetic_scalar.driver import ( MagneticForwardDriver, - MagneticInversionDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) TARGET = 1132.1998 @@ -38,27 +39,28 @@ def test_automesh( refinement=(4, 8), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - flatten=False, - ) - inducing_field = (49999.8, 90.0, 0.0) - params = MagneticForwardOptions.build( - forward_only=True, - geoh5=geoh5, - mesh=None, - topography_object=topography, - inducing_field_strength=inducing_field[0], - inducing_field_inclination=inducing_field[1], - inducing_field_declination=inducing_field[2], - data_object=survey, - starting_model=model, + opts = SyntheticsComponentsOptions( + method="magnetic_scalar", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0 + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.05), ) + with Workspace.create(tmp_path / "forward_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + inducing_field = (49999.8, 90.0, 0.0) + params = MagneticForwardOptions.build( + forward_only=True, + geoh5=geoh5, + mesh=None, + topography_object=components.topography, + inducing_field_strength=inducing_field[0], + inducing_field_inclination=inducing_field[1], + inducing_field_declination=inducing_field[2], + data_object=components.survey, + starting_model=components.model, + ) fwr_driver = MagneticForwardDriver(params) fwr_driver.run() diff --git a/tests/run_tests/driver_mag_test.py b/tests/run_tests/driver_mag_test.py index 11e70f45..a6d95cf0 100644 --- a/tests/run_tests/driver_mag_test.py +++ b/tests/run_tests/driver_mag_test.py @@ -16,7 +16,6 @@ from dask.distributed import LocalCluster, performance_report from geoh5py.workspace import Workspace -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import ( MagneticForwardOptions, MagneticInversionOptions, @@ -25,8 +24,19 @@ MagneticForwardDriver, MagneticInversionDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -41,29 +51,29 @@ def test_susceptibility_fwr_run( refinement=(2,), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - inversion_type="magnetic", - flatten=False, - ) - inducing_field = (49999.8, 90.0, 0.0) - - params = MagneticForwardOptions.build( - forward_only=True, - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - inducing_field_strength=inducing_field[0], - inducing_field_inclination=inducing_field[1], - inducing_field_declination=inducing_field[2], - data_object=survey, - starting_model=model, + opts = SyntheticsComponentsOptions( + method="magnetic", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0 + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.05), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + inducing_field = (49999.8, 90.0, 0.0) + + params = MagneticForwardOptions.build( + forward_only=True, + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + inducing_field_strength=inducing_field[0], + inducing_field_inclination=inducing_field[1], + inducing_field_declination=inducing_field[2], + data_object=components.survey, + starting_model=components.model, + ) # params.workpath = tmp_path fwr_driver = MagneticForwardDriver(params) fwr_driver.run() @@ -88,7 +98,8 @@ def test_susceptibility_run( with Workspace(workpath) as geoh5: tmi = geoh5.get_entity("Iteration_0_tmi")[0] orig_tmi = tmi.values.copy() - mesh = geoh5.get_entity("mesh")[0] + components = SyntheticsComponents(geoh5=geoh5) + mesh = components.mesh active_cells = mesh.get_entity("active_cells")[0] inducing_field = (50000.0, 90.0, 0.0) diff --git a/tests/run_tests/driver_mt_test.py b/tests/run_tests/driver_mt_test.py index 5933f942..51f85961 100644 --- a/tests/run_tests/driver_mt_test.py +++ b/tests/run_tests/driver_mt_test.py @@ -15,7 +15,6 @@ from pathlib import Path import numpy as np -from dask.distributed import LocalCluster, performance_report from geoh5py.groups import SimPEGGroup from geoh5py.workspace import Workspace @@ -27,9 +26,19 @@ MTForwardOptions, MTInversionOptions, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -92,35 +101,31 @@ def test_magnetotellurics_fwr_run( cell_size=(20.0, 20.0, 20.0), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.01, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - cell_size=cell_size, - drape_height=0.0, - inversion_type="magnetotellurics", - flatten=False, - ) - params = MTForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - background_conductivity=1e-2, - zxx_real_channel_bool=True, - zxx_imag_channel_bool=True, - zxy_real_channel_bool=True, - zxy_imag_channel_bool=True, - zyx_real_channel_bool=True, - zyx_imag_channel_bool=True, - zyy_real_channel_bool=True, - zyy_imag_channel_bool=True, - solver_type="Mumps", + opts = SyntheticsComponentsOptions( + method="magnetotellurics", + survey=SurveyOptions(n_stations=n_grid_points, n_lines=n_grid_points), + mesh=MeshOptions(cell_size=cell_size, refinement=refinement), + model=ModelOptions(background=0.01), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + params = MTForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + background_conductivity=1e-2, + zxx_real_channel_bool=True, + zxx_imag_channel_bool=True, + zxy_real_channel_bool=True, + zxy_imag_channel_bool=True, + zyx_real_channel_bool=True, + zyx_imag_channel_bool=True, + zyy_real_channel_bool=True, + zyy_imag_channel_bool=True, + solver_type="Mumps", + ) fwr_driver = MTForwardDriver(params) fwr_driver.run() @@ -137,13 +142,10 @@ def test_magnetotellurics_run(tmp_path: Path, max_iterations=1, pytest=True): ) with Workspace(workpath) as geoh5: - survey = next( - child - for child in geoh5.get_entity("survey") - if not isinstance(child.parent, SimPEGGroup) - ) - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5) + survey = components.survey + mesh = components.mesh + topography = components.topography data_kwargs = setup_data(geoh5, survey) orig_zyy_real_1 = geoh5.get_entity("Iteration_0_zyy_real_[0]")[0].values diff --git a/tests/run_tests/driver_mvi_test.py b/tests/run_tests/driver_mvi_test.py index 73a13ee8..48fc2cec 100644 --- a/tests/run_tests/driver_mvi_test.py +++ b/tests/run_tests/driver_mvi_test.py @@ -19,9 +19,7 @@ from geoh5py.groups.property_group import GroupTypeEnum from geoh5py.objects import Curve from geoh5py.workspace import Workspace -from pytest import warns -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import ( MVIForwardOptions, MVIInversionOptions, @@ -30,8 +28,19 @@ MVIForwardDriver, MVIInversionDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -46,33 +55,37 @@ def test_magnetic_vector_fwr_run( refinement=(2,), ): # Run the forward - geoh5, _, model, points, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0 + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.05), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - # Unitest dealing with Curve - with geoh5.open(): - survey = Curve.create(geoh5, name=points.name, vertices=points.vertices) - geoh5.remove_entity(points) - inducing_field = (50000.0, 90.0, 0.0) - params = MVIForwardOptions.build( - forward_only=True, - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - inducing_field_strength=inducing_field[0], - inducing_field_inclination=inducing_field[1], - inducing_field_declination=inducing_field[2], - data_object=survey, - starting_model=model, - starting_inclination=45, - starting_declination=270, - ) + # Unitest dealing with Curve + + _ = Curve.create( + geoh5, name=components.survey.name, vertices=components.survey.vertices + ) + geoh5.remove_entity(components.survey) + inducing_field = (50000.0, 90.0, 0.0) + params = MVIForwardOptions.build( + forward_only=True, + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + inducing_field_strength=inducing_field[0], + inducing_field_inclination=inducing_field[1], + inducing_field_declination=inducing_field[2], + data_object=components.survey, + starting_model=components.model, + starting_inclination=45, + starting_declination=270, + ) fwr_driver = MVIForwardDriver(params) fwr_driver.run() @@ -95,8 +108,9 @@ def test_magnetic_vector_run( with Workspace(workpath) as geoh5: tmi = geoh5.get_entity("Iteration_0_tmi")[0] orig_tmi = tmi.values.copy() - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5=geoh5) + mesh = components.mesh + topography = components.topography inducing_field = (50000.0, 90.0, 0.0) dip, direction = mesh.add_data( { diff --git a/tests/run_tests/driver_rotated_gradients_test.py b/tests/run_tests/driver_rotated_gradients_test.py index 958752bb..8b29ded0 100644 --- a/tests/run_tests/driver_rotated_gradients_test.py +++ b/tests/run_tests/driver_rotated_gradients_test.py @@ -11,16 +11,13 @@ from __future__ import annotations from pathlib import Path -from unittest.mock import patch import numpy as np from geoapps_utils.modelling.plates import PlateModel -from geoapps_utils.utils.importing import GeoAppsError +from geoapps_utils.utils.locations import gaussian from geoh5py.groups.property_group import PropertyGroup from geoh5py.workspace import Workspace -from pytest import raises -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import ( GravityForwardOptions, GravityInversionOptions, @@ -29,8 +26,19 @@ GravityForwardDriver, GravityInversionDriver, ) -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -45,34 +53,40 @@ def test_gravity_rotated_grad_fwr_run( refinement=(2,), ): # Run the forward - plate_model = PlateModel( - strike_length=500.0, - dip_length=150.0, - width=20.0, - origin=(0.0, 0.0, -10.0), - direction=60.0, - dip=70.0, - ) - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - plate_model=plate_model, - background=0.0, - anomaly=0.75, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - center=(0.0, 0.0, 15.0), - flatten=False, - ) - params = GravityForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - gz_channel_bool=True, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=n_grid_points, + n_lines=n_grid_points, + center=(0.0, 0.0), + drape=5.0, + topography=lambda x, y: gaussian(x, y, amplitude=50.0, width=100.0) + 15, + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions( + anomaly=0.75, + plate=PlateModel( + strike_length=500.0, + dip_length=150.0, + width=20.0, + origin=(0.0, 0.0, -10.0), + direction=60.0, + dip=70.0, + ), + ), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + + params = GravityForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + gz_channel_bool=True, + ) fwr_driver = GravityForwardDriver(params) fwr_driver.run() @@ -93,8 +107,9 @@ def test_rotated_grad_run( with Workspace(workpath) as geoh5: gz = geoh5.get_entity("Iteration_0_gz")[0] orig_gz = gz.values.copy() - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5=geoh5) + mesh = components.mesh + topography = components.topography # Create property group with orientation dip = np.ones(mesh.n_cells) * 70 diff --git a/tests/run_tests/driver_tile_estimator_test.py b/tests/run_tests/driver_tile_estimator_test.py index 0c74ef21..ea3edcc6 100644 --- a/tests/run_tests/driver_tile_estimator_test.py +++ b/tests/run_tests/driver_tile_estimator_test.py @@ -13,15 +13,21 @@ from pathlib import Path import numpy as np +from geoh5py import Workspace -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import MagneticInversionOptions from simpeg_drivers.potential_fields.magnetic_scalar.driver import ( MagneticInversionDriver, ) +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) from simpeg_drivers.utils.tile_estimate import TileEstimator, TileParameters from simpeg_drivers.utils.utils import simpeg_group_to_driver -from tests.testing_utils import setup_inversion_workspace def test_tile_estimator_run( @@ -31,32 +37,31 @@ def test_tile_estimator_run( ): inducing_field = (49999.8, 90.0, 0.0) # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=refinement, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - flatten=False, + + opts = SyntheticsComponentsOptions( + method="magnetic_scalar", + survey=SurveyOptions(n_stations=n_grid_points, n_lines=n_grid_points), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.05), ) - with geoh5.open(): - tmi_channel = survey.add_data( + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + tmi_channel = components.survey.add_data( { - "tmi": {"values": np.random.rand(survey.n_vertices)}, + "tmi": {"values": np.random.rand(components.survey.n_vertices)}, } ) - params = MagneticInversionOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - inducing_field_strength=inducing_field[0], - inducing_field_inclination=inducing_field[1], - inducing_field_declination=inducing_field[2], - data_object=survey, - tmi_channel=tmi_channel, - starting_model=model, - ) + params = MagneticInversionOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + inducing_field_strength=inducing_field[0], + inducing_field_inclination=inducing_field[1], + inducing_field_declination=inducing_field[2], + data_object=components.survey, + tmi_channel=tmi_channel, + starting_model=components.model, + ) driver = MagneticInversionDriver(params) tile_params = TileParameters(geoh5=geoh5, simulation=driver.out_group) diff --git a/tests/run_tests/driver_tipper_test.py b/tests/run_tests/driver_tipper_test.py index 12dcc2a0..e865f7cb 100644 --- a/tests/run_tests/driver_tipper_test.py +++ b/tests/run_tests/driver_tipper_test.py @@ -24,9 +24,19 @@ TipperForwardDriver, TipperInversionDriver, ) -from simpeg_drivers.options import ActiveCellsOptions -from simpeg_drivers.utils.utils import get_inversion_output -from tests.testing_utils import check_target, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import ( + SyntheticsComponents, +) +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) +from tests.utils.targets import ( + check_target, + get_inversion_output, +) # To test the full run and validate the inversion. @@ -42,32 +52,30 @@ def test_tipper_fwr_run( cell_size=(20.0, 20.0, 20.0), ): # Run the forward - geoh5, _, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=100.0, - anomaly=1.0, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - cell_size=cell_size, - inversion_type="tipper", - drape_height=15.0, - flatten=False, + opts = SyntheticsComponentsOptions( + method="tipper", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=15.0 + ), + mesh=MeshOptions(cell_size=cell_size, refinement=refinement), + model=ModelOptions(background=100.0), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - params = TipperForwardOptions.build( - geoh5=geoh5, - mesh=model.parent, - topography_object=topography, - data_object=survey, - starting_model=model, - model_type="Resistivity (Ohm-m)", - background_conductivity=1e2, - txz_real_channel_bool=True, - txz_imag_channel_bool=True, - tyz_real_channel_bool=True, - tyz_imag_channel_bool=True, - ) + params = TipperForwardOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=components.survey, + starting_model=components.model, + model_type="Resistivity (Ohm-m)", + background_conductivity=1e2, + txz_real_channel_bool=True, + txz_imag_channel_bool=True, + tyz_real_channel_bool=True, + tyz_imag_channel_bool=True, + ) fwr_driver = TipperForwardDriver(params) @@ -82,13 +90,10 @@ def test_tipper_run(tmp_path: Path, max_iterations=1, pytest=True): workpath = tmp_path.parent / "test_tipper_fwr_run0" / "inversion_test.ui.geoh5" with Workspace(workpath) as geoh5: - survey = next( - child - for child in geoh5.get_entity("survey") - if not isinstance(child.parent, SimPEGGroup) - ) - mesh = geoh5.get_entity("mesh")[0] - topography = geoh5.get_entity("topography")[0] + components = SyntheticsComponents(geoh5=geoh5) + survey = components.survey + mesh = components.mesh + topography = components.topography data = {} uncertainties = {} diff --git a/tests/run_tests/sensitivity_cutoff_test.py b/tests/run_tests/sensitivity_cutoff_test.py index 73c16292..945f5e3b 100644 --- a/tests/run_tests/sensitivity_cutoff_test.py +++ b/tests/run_tests/sensitivity_cutoff_test.py @@ -21,10 +21,15 @@ from simpeg_drivers.depth_of_investigation.sensitivity_cutoff.options import ( SensitivityCutoffOptions, ) -from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import GravityInversionOptions from simpeg_drivers.potential_fields.gravity.driver import GravityInversionDriver -from tests.testing_utils import setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) def setup_inversion_results( @@ -32,36 +37,38 @@ def setup_inversion_results( n_grid_points=2, refinement=(2,), ): - geoh5, mesh, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.75, - n_electrodes=n_grid_points, - n_lines=n_grid_points, - refinement=refinement, - flatten=False, + opts = SyntheticsComponentsOptions( + method="gravity", + survey=SurveyOptions( + n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0 + ), + mesh=MeshOptions(refinement=refinement), + model=ModelOptions(anomaly=0.75), ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) - # Run the inverse with save_sensitivities=True - with geoh5.open(): - gz = survey.add_data({"gz": {"values": np.random.randn(len(survey.vertices))}}) - - params = GravityInversionOptions.build( - geoh5=geoh5, - mesh=mesh, - topography_object=topography, - data_object=gz.parent, - starting_model=1e-4, - reference_model=0.0, - s_norm=0.0, - gz_channel=gz, - gz_uncertainty=2e-3, - lower_bound=0.0, - max_global_iterations=1, - initial_beta_ratio=1e-2, - percentile=100, - save_sensitivities=True, - ) + # Run the inverse with save_sensitivities=True + gz = components.survey.add_data( + {"gz": {"values": np.random.randn(len(components.survey.vertices))}} + ) + + params = GravityInversionOptions.build( + geoh5=geoh5, + mesh=components.mesh, + topography_object=components.topography, + data_object=gz.parent, + starting_model=1e-4, + reference_model=0.0, + s_norm=0.0, + gz_channel=gz, + gz_uncertainty=2e-3, + lower_bound=0.0, + max_global_iterations=1, + initial_beta_ratio=1e-2, + percentile=100, + save_sensitivities=True, + ) params.write_ui_json(path=tmp_path / "Inv_run.ui.json") GravityInversionDriver.start(str(tmp_path / "Inv_run.ui.json")) @@ -74,8 +81,9 @@ def test_sensitivity_percent_cutoff_run(tmp_path): ) with Workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5) sensitivity = geoh5.get_entity("Iteration_1_sensitivities")[0] - mesh = sensitivity.parent + mesh = components.mesh params = SensitivityCutoffOptions( geoh5=geoh5, mesh=mesh, @@ -100,8 +108,9 @@ def test_sensitivity_cutoff_percentile_run(tmp_path): ) with Workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5) sensitivity = geoh5.get_entity("Iteration_1_sensitivities")[0] - mesh = sensitivity.parent + mesh = components.mesh params = SensitivityCutoffOptions( geoh5=geoh5, mesh=mesh, @@ -128,8 +137,9 @@ def test_sensitivity_cutoff_log_percent_run(tmp_path): ) with Workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5) sensitivity = geoh5.get_entity("Iteration_1_sensitivities")[0] - mesh = sensitivity.parent + mesh = components.mesh params = SensitivityCutoffOptions( geoh5=geoh5, mesh=mesh, diff --git a/tests/testing_utils.py b/tests/testing_utils.py deleted file mode 100644 index da542e98..00000000 --- a/tests/testing_utils.py +++ /dev/null @@ -1,532 +0,0 @@ -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' -# Copyright (c) 2025 Mira Geoscience Ltd. ' -# ' -# This file is part of simpeg-drivers package. ' -# ' -# simpeg-drivers is distributed under the terms and conditions of the MIT License ' -# (see LICENSE file at the root of this source code package). ' -# ' -# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -import warnings -from pathlib import Path -from uuid import UUID - -import numpy as np -from discretize.utils import mesh_builder_xyz -from geoapps_utils.modelling.plates import PlateModel, make_plate -from geoh5py import Workspace -from geoh5py.objects import ( - AirborneFEMReceivers, - AirborneFEMTransmitters, - AirborneTEMReceivers, - AirborneTEMTransmitters, - CurrentElectrode, - DrapeModel, - LargeLoopGroundTEMReceivers, - LargeLoopGroundTEMTransmitters, - MTReceivers, - Points, - PotentialElectrode, - Surface, - TipperBaseStations, - TipperReceivers, -) -from octree_creation_app.driver import OctreeDriver -from octree_creation_app.utils import treemesh_2_octree -from scipy.spatial import Delaunay -from simpeg import utils - -from simpeg_drivers.utils.utils import active_from_xyz, get_drape_model - - -class Geoh5Tester: - """Create temp workspace, copy entities, and setup params class.""" - - def __init__(self, geoh5, path, name, params_class=None): - self.geoh5 = geoh5 - self.tmp_path = Path(path) / name - self.ws = Workspace.create(self.tmp_path) - - if params_class is not None: - self.params = params_class(validate=False, geoh5=self.ws) - self.has_params = True - else: - self.has_params = False - - def copy_entity(self, uid): - entity = self.ws.get_entity(uid) - if not entity or entity[0] is None: - return self.geoh5.get_entity(uid)[0].copy(parent=self.ws) - return entity[0] - - def set_param(self, param, value): - if self.has_params: - try: - uid = UUID(value) - entity = self.copy_entity(uid) - setattr(self.params, param, entity) - except (AttributeError, ValueError): - setattr(self.params, param, value) - else: - msg = "No params class has been initialized." - raise (ValueError(msg)) - - def make(self): - if self.has_params: - return self.ws, self.params - else: - return self.ws - - -def check_target(output: dict, target: dict, tolerance=0.05): - """ - Check inversion output metrics against hard-valued target. - :param output: Dictionary containing keys for 'data', 'phi_d' and 'phi_m'. - :param target: Dictionary containing keys for 'data_norm', 'phi_d' and 'phi_m'.\ - :param tolerance: Tolerance between output and target measured as: |a-b|/b - """ - print( - f"Output: 'data_norm': {np.linalg.norm(output['data'])}, 'phi_d': {output['phi_d'][1]}, 'phi_m': {output['phi_m'][1]}" - ) - print(f"Target: {target}") - - if any(np.isnan(output["data"])): - warnings.warn( - "Skipping data norm comparison due to nan (used to bypass lone faulty test run in GH actions)." - ) - else: - np.testing.assert_array_less( - np.abs(np.linalg.norm(output["data"]) - target["data_norm"]) - / target["data_norm"], - tolerance, - ) - - np.testing.assert_array_less( - np.abs(output["phi_m"][1] - target["phi_m"]) / target["phi_m"], tolerance - ) - np.testing.assert_array_less( - np.abs(output["phi_d"][1] - target["phi_d"]) / target["phi_d"], tolerance - ) - - -def gaussian_topo_drape(x, y, flatten=False): - """Topography Gaussian function""" - b = 100 - A = 50 - if flatten: - return np.zeros_like(x) - else: - return A * np.exp(-0.5 * ((x / b) ** 2.0 + (y / b) ** 2.0)) - - -def generate_dc_survey( - workspace: Workspace, - x_loc: np.ndarray, - y_loc: np.ndarray, - z_loc: np.ndarray | None = None, -) -> PotentialElectrode: - """ - Utility function to generate a DC survey. - """ - # Create sources along line - if z_loc is None: - z_loc = np.zeros_like(x_loc) - - vertices = np.c_[x_loc.ravel(), y_loc.ravel(), z_loc.ravel()] - parts = np.kron(np.arange(x_loc.shape[0]), np.ones(x_loc.shape[1])).astype("int") - currents = CurrentElectrode.create(workspace, vertices=vertices, parts=parts) - currents.add_default_ab_cell_id() - n_dipoles = 9 - dipoles = [] - current_id = [] - - for val in currents.ab_cell_id.values: - cell_id = int(currents.ab_map[val]) - 1 - - for dipole in range(n_dipoles): - dipole_ids = currents.cells[cell_id, :] + 2 + dipole - - if ( - any(dipole_ids > (currents.n_vertices - 1)) - or len( - np.unique(parts[np.r_[currents.cells[cell_id, 0], dipole_ids[1]]]) - ) - > 1 - ): - continue - - dipoles += [dipole_ids] - current_id += [val] - - potentials = PotentialElectrode.create( - workspace, vertices=vertices, cells=np.vstack(dipoles).astype("uint32") - ) - line_id = potentials.vertices[potentials.cells[:, 0], 1] - line_id = (line_id - np.min(line_id) + 1).astype(np.int32) - line_reference = {0: "Unknown"} - line_reference.update({k: str(k) for k in np.unique(line_id)}) - potentials.add_data( - { - "line_ids": { - "values": line_id, - "type": "REFERENCED", - "value_map": line_reference, - } - } - ) - potentials.ab_cell_id = np.hstack(current_id).astype("int32") - potentials.current_electrodes = currents - currents.potential_electrodes = potentials - - return potentials - - -def generate_fdem_survey(geoh5, vertices): - survey = AirborneFEMReceivers.create(geoh5, vertices=vertices, name="Airborne_rx") - freq_metadata = [ - {"Coaxial data": False, "Frequency": 900, "Offset": 7.86}, - {"Coaxial data": False, "Frequency": 7200, "Offset": 7.86}, - {"Coaxial data": False, "Frequency": 56000, "Offset": 6.3}, - ] - survey.metadata["EM Dataset"]["Frequency configurations"] = freq_metadata - - tx_locs = [] - freqs = [] - for meta in freq_metadata: - tx_vertices = vertices.copy() - tx_vertices[:, 0] -= meta["Offset"] - tx_locs.append(tx_vertices) - freqs.append([[meta["Frequency"]] * len(vertices)]) - tx_locs = np.vstack(tx_locs) - freqs = np.hstack(freqs).flatten() - - transmitters = AirborneFEMTransmitters.create( - geoh5, vertices=tx_locs, name="Airborne_tx" - ) - survey.transmitters = transmitters - survey.channels = [900.0, 7200.0, 56000.0] - - transmitters.add_data( - { - "Tx frequency": { - "values": freqs, - "association": "VERTEX", - "primitive_type": "REFERENCED", - "value_map": {k: str(k) for k in freqs}, - } - } - ) - - dist = np.linalg.norm( - survey.vertices[survey.cells[:, 0], :] - survey.vertices[survey.cells[:, 1], :], - axis=1, - ) - survey.remove_cells(np.where(dist > 200.0)[0]) - dist = np.linalg.norm( - transmitters.vertices[transmitters.cells[:, 0], :] - - transmitters.vertices[transmitters.cells[:, 1], :], - axis=1, - ) - transmitters.remove_cells(np.where(dist > 200.0)[0]) - - return survey - - -def generate_tdem_survey(geoh5, vertices, n_lines, flatten=False, airborne=False): - if airborne: - survey = AirborneTEMReceivers.create( - geoh5, vertices=vertices, name="Airborne_rx" - ) - transmitters = AirborneTEMTransmitters.create( - geoh5, vertices=vertices, name="Airborne_tx" - ) - - dist = np.linalg.norm( - survey.vertices[survey.cells[:, 0], :] - - survey.vertices[survey.cells[:, 1], :], - axis=1, - ) - survey.remove_cells(np.where(dist > 200.0)[0]) - transmitters.remove_cells(np.where(dist > 200.0)[0]) - survey.transmitters = transmitters - else: - X = vertices[:, 0].reshape((n_lines, -1)) - Y = vertices[:, 1].reshape((n_lines, -1)) - Z = vertices[:, 2].reshape((n_lines, -1)) - center = np.mean(vertices, axis=0) - arrays = [ - np.c_[ - X[: int(n_lines / 2), :].flatten(), - Y[: int(n_lines / 2), :].flatten(), - Z[: int(n_lines / 2), :].flatten(), - ], - np.c_[ - X[int(n_lines / 2) :, :].flatten(), - Y[int(n_lines / 2) :, :].flatten(), - Z[int(n_lines / 2) :, :].flatten(), - ], - ] - loops = [] - loop_cells = [] - loop_id = [] - count = 0 - for ind, array in enumerate(arrays): - loop_id += [np.ones(array.shape[0]) * (ind + 1)] - min_loc = np.min(array, axis=0) - max_loc = np.max(array, axis=0) - loop = np.vstack( - [ - np.c_[ - np.ones(5) * min_loc[0], - np.linspace(min_loc[1], max_loc[1], 5), - ], - np.c_[ - np.linspace(min_loc[0], max_loc[0], 5)[1:], - np.ones(4) * max_loc[1], - ], - np.c_[ - np.ones(4) * max_loc[0], - np.linspace(max_loc[1], min_loc[1], 5)[1:], - ], - np.c_[ - np.linspace(max_loc[0], min_loc[0], 5)[1:-1], - np.ones(3) * min_loc[1], - ], - ] - ) - loop = (loop - np.mean(loop, axis=0)) * 1.5 + np.mean(loop, axis=0) - loop = np.c_[ - loop, gaussian_topo_drape(loop[:, 0], loop[:, 1], flatten=flatten) - ] - loops += [loop + np.asarray(center)] - loop_cells += [np.c_[np.arange(15) + count, np.arange(15) + count + 1]] - loop_cells += [np.c_[count + 15, count]] - count += 16 - - transmitters = LargeLoopGroundTEMTransmitters.create( - geoh5, - vertices=np.vstack(loops), - cells=np.vstack(loop_cells), - ) - transmitters.tx_id_property = transmitters.parts + 1 - survey = LargeLoopGroundTEMReceivers.create(geoh5, vertices=np.vstack(vertices)) - survey.transmitters = transmitters - survey.tx_id_property = np.hstack(loop_id) - # survey.parts = np.repeat(np.arange(n_lines), n_electrodes) - - survey.channels = np.r_[3e-04, 6e-04, 1.2e-03] * 1e3 - waveform = np.c_[ - np.r_[ - np.arange(-0.002, -0.0001, 5e-4), - np.arange(-0.0004, 0.0, 1e-4), - np.arange(0.0, 0.002, 5e-4), - ] - * 1e3 - + 2.0, - np.r_[np.linspace(0, 1, 4), np.linspace(0.9, 0.0, 4), np.zeros(4)], - ] - survey.waveform = waveform - survey.timing_mark = 2.0 - survey.unit = "Milliseconds (ms)" - - return survey - - -plate_model_default = PlateModel( - strike_length=40.0, - dip_length=40.0, - width=40.0, - origin=(0.0, 0.0, 10.0), -) - - -def setup_inversion_workspace( - work_dir, - plate_model: PlateModel = plate_model_default, - background=None, - anomaly=None, - cell_size=(5.0, 5.0, 5.0), - center=(0.0, 0.0, 0.0), - n_electrodes=20, - n_lines=5, - refinement=(4, 6), - x_limits=(-100.0, 100.0), - y_limits=(-100.0, 100.0), - padding_distance=100, - drape_height=5.0, - inversion_type="other", - flatten=False, - geoh5=None, -): - """ - Creates a workspace populated with objects to simulate/invert a simple model. - - rot_xyz: Rotation angles in degrees about the x, y, and z axes. - """ - filepath = Path(work_dir) / "inversion_test.ui.geoh5" - if geoh5 is None: - if filepath.is_file(): - filepath.unlink() - geoh5 = Workspace.create(filepath) - # Topography - xx, yy = np.meshgrid(np.linspace(-200.0, 200.0, 50), np.linspace(-200.0, 200.0, 50)) - - zz = gaussian_topo_drape(xx, yy, flatten=flatten) - topo = np.c_[ - utils.mkvc(xx) + center[0], - utils.mkvc(yy) + center[1], - utils.mkvc(zz) + center[2], - ] - triang = Delaunay(topo[:, :2]) - topography = Surface.create( - geoh5, - vertices=topo, - cells=triang.simplices, # pylint: disable=E1101 - name="topography", - ) - - # Observation points - n_electrodes = ( - 4 - if ("polarization" in inversion_type or "current" in inversion_type) - & (n_electrodes < 4) - else n_electrodes - ) - xr = np.linspace(x_limits[0], x_limits[1], int(n_electrodes)) - yr = np.linspace(y_limits[0], y_limits[1], int(n_lines)) - X, Y = np.meshgrid(xr, yr) - Z = gaussian_topo_drape(X, Y, flatten=flatten) + drape_height - - vertices = np.c_[ - utils.mkvc(X.T) + center[0], - utils.mkvc(Y.T) + center[1], - utils.mkvc(Z.T) + center[2], - ] - - if "polarization" in inversion_type or "current" in inversion_type: - survey = generate_dc_survey(geoh5, X, Y, Z) - - elif "magnetotellurics" in inversion_type: - survey = MTReceivers.create( - geoh5, - vertices=vertices, - name="survey", - components=[ - "Zxx (real)", - "Zxx (imag)", - "Zxy (real)", - "Zxy (imag)", - "Zyx (real)", - "Zyx (imag)", - "Zyy (real)", - "Zyy (imag)", - ], - channels=[10.0, 100.0, 1000.0], - ) - - elif "tipper" in inversion_type: - survey = TipperReceivers.create( - geoh5, - vertices=vertices, - name="survey", - components=[ - "Txz (real)", - "Txz (imag)", - "Tyz (real)", - "Tyz (imag)", - ], - ) - survey.base_stations = TipperBaseStations.create( - geoh5, vertices=np.c_[vertices[0, :]].T - ) - survey.channels = [10.0, 100.0, 1000.0] - dist = np.linalg.norm( - survey.vertices[survey.cells[:, 0], :] - - survey.vertices[survey.cells[:, 1], :], - axis=1, - ) - # survey.cells = survey.cells[dist < 100.0, :] - survey.remove_cells(np.where(dist > 200)[0]) - - elif inversion_type in ["fdem", "fem", "fdem 1d"]: - survey = generate_fdem_survey(geoh5, vertices) - - elif "tdem" in inversion_type: - survey = generate_tdem_survey( - geoh5, - vertices, - n_lines, - flatten=flatten, - airborne="airborne" in inversion_type, - ) - - else: - survey = Points.create( - geoh5, - vertices=vertices, - name="survey", - ) - - # Create a mesh - if "2d" in inversion_type: - lines = survey.get_entity("line_ids")[0].values - entity, mesh, _ = get_drape_model( # pylint: disable=unbalanced-tuple-unpacking - geoh5, - "Models", - survey.vertices[np.unique(survey.cells[lines == 101, :]), :], - [cell_size[0], cell_size[2]], - 100.0, - [padding_distance] * 2 + [padding_distance] * 2, - 1.1, - parent=None, - return_colocated_mesh=True, - return_sorting=True, - ) - active = active_from_xyz(entity, topography.vertices, grid_reference="top") - - else: - padDist = np.ones((3, 2)) * padding_distance - mesh = mesh_builder_xyz( - vertices - np.r_[cell_size] / 2.0, - cell_size, - depth_core=100.0, - padding_distance=padDist, - mesh_type="TREE", - tree_diagonal_balance=False, - ) - mesh = OctreeDriver.refine_tree_from_surface( - mesh, - topography, - levels=refinement, - finalize=False, - ) - - if inversion_type in ["fdem", "airborne tdem"]: - mesh = OctreeDriver.refine_tree_from_points( - mesh, - vertices, - levels=[2], - finalize=False, - ) - - mesh.finalize() - entity = treemesh_2_octree(geoh5, mesh, name="mesh") - active = active_from_xyz(entity, topography.vertices, grid_reference="top") - - # Create the Model - cell_centers = entity.centroids.copy() - - model = make_plate( - cell_centers, plate_model, background=background, anomaly=anomaly - ) - - if "1d" in inversion_type: - model = background * np.ones(mesh.nC) - model[(mesh.cell_centers[:, 2] < 0) & (mesh.cell_centers[:, 2] > -20)] = anomaly - - model[~active] = np.nan - model = entity.add_data({"model": {"values": model}}) - geoh5.close() - return geoh5, entity, model, survey, topography diff --git a/tests/topography_test.py b/tests/topography_test.py index 26fc148f..c3c0b2bb 100644 --- a/tests/topography_test.py +++ b/tests/topography_test.py @@ -13,47 +13,52 @@ from pathlib import Path import numpy as np +from geoh5py import Workspace from simpeg_drivers.components import InversionData, InversionMesh, InversionTopography from simpeg_drivers.options import ActiveCellsOptions from simpeg_drivers.potential_fields import MVIInversionOptions -from tests.testing_utils import Geoh5Tester, setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) def test_get_locations(tmp_path: Path): - geoh5, mesh, model, survey, topography = setup_inversion_workspace( - tmp_path, - background=0.0, - anomaly=0.05, - refinement=(2,), - n_electrodes=2, - n_lines=2, - inversion_type="magnetic_vector", + opts = SyntheticsComponentsOptions( + method="magnetic_vector", + survey=SurveyOptions(n_stations=2, n_lines=2), + mesh=MeshOptions(refinement=(2,)), + model=ModelOptions(anomaly=0.05), ) - with geoh5.open(): - tmi_channel, gyz_channel = survey.add_data( + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + tmi_channel, gyz_channel = components.survey.add_data( { - "tmi": {"values": np.random.rand(survey.n_vertices)}, - "gyz": {"values": np.random.rand(survey.n_vertices)}, + "tmi": {"values": np.random.rand(components.survey.n_vertices)}, + "gyz": {"values": np.random.rand(components.survey.n_vertices)}, } ) - elevation = topography.add_data( - {"elevation": {"values": topography.vertices[:, 2]}} + elevation = components.topography.add_data( + {"elevation": {"values": components.topography.vertices[:, 2]}} + ) + params = MVIInversionOptions.build( + geoh5=geoh5, + mesh=components.mesh, + data_object=components.survey, + tmi_channel=tmi_channel, + tmi_uncertainty=0.01, + active_cells=ActiveCellsOptions( + topography_object=components.topography, topography=elevation + ), + inducing_field_strength=50000.0, + inducing_field_inclination=60.0, + inducing_field_declination=0.0, + starting_model=1.0, ) - params = MVIInversionOptions.build( - geoh5=geoh5, - mesh=mesh, - data_object=survey, - tmi_channel=tmi_channel, - tmi_uncertainty=0.01, - active_cells=ActiveCellsOptions( - topography_object=topography, topography=elevation - ), - inducing_field_strength=50000.0, - inducing_field_inclination=60.0, - inducing_field_declination=0.0, - starting_model=1.0, - ) geoh5 = params.geoh5 with geoh5.open(): @@ -66,14 +71,14 @@ def test_get_locations(tmp_path: Path): # Check that boundary cells are handled properly # Shift one of the survey vertices to the corner - survey.vertices[0, :] = mesh.centroids[0, :] - geoh5.update_attribute(survey, "vertices") + components.survey.vertices[0, :] = components.mesh.centroids[0, :] + geoh5.update_attribute(components.survey, "vertices") inversion_mesh = InversionMesh(geoh5, params) inversion_data = InversionData(geoh5, params) params.inversion_type = "magnetotellurics" active_cells = topo.active_cells(inversion_mesh, inversion_data) - mesh.add_data( + components.mesh.add_data( { "active_cells": { "values": active_cells, diff --git a/tests/uijson_test.py b/tests/uijson_test.py index 1a63769a..dc1a4089 100644 --- a/tests/uijson_test.py +++ b/tests/uijson_test.py @@ -25,11 +25,17 @@ import simpeg_drivers from simpeg_drivers.driver import InversionDriver from simpeg_drivers.line_sweep.driver import LineSweepDriver -from simpeg_drivers.options import ActiveCellsOptions, Deprecations, IRLSOptions +from simpeg_drivers.options import Deprecations, IRLSOptions from simpeg_drivers.potential_fields.gravity.options import GravityInversionOptions from simpeg_drivers.potential_fields.gravity.uijson import GravityInversionUIJson from simpeg_drivers.uijson import SimPEGDriversUIJson -from tests.testing_utils import setup_inversion_workspace +from simpeg_drivers.utils.synthetics.driver import SyntheticsComponents +from simpeg_drivers.utils.synthetics.options import ( + MeshOptions, + ModelOptions, + SurveyOptions, + SyntheticsComponentsOptions, +) logger = logging.getLogger(__name__) @@ -235,23 +241,28 @@ def test_gravity_uijson(tmp_path): import warnings warnings.filterwarnings("error") - geoh5, _, starting_model, survey, topography = setup_inversion_workspace( - tmp_path, background=0.0, anomaly=0.75, inversion_type="gravity" - ) - with geoh5.open(): - gz_channel = survey.add_data({"gz": {"values": np.ones(survey.n_vertices)}}) - gz_uncerts = survey.add_data({"gz_unc": {"values": np.ones(survey.n_vertices)}}) - - opts = GravityInversionOptions.build( - version="old news", - geoh5=geoh5, - data_object=survey, - gz_channel=gz_channel, - gz_uncertainty=gz_uncerts, - mesh=starting_model.parent, - starting_model=starting_model, - topography_object=topography, + opts = SyntheticsComponentsOptions( + method="gravity", model=ModelOptions(anomaly=0.75) ) + with Workspace.create(tmp_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) + gz_channel = components.survey.add_data( + {"gz": {"values": np.ones(components.survey.n_vertices)}} + ) + gz_uncerts = components.survey.add_data( + {"gz_unc": {"values": np.ones(components.survey.n_vertices)}} + ) + + opts = GravityInversionOptions.build( + version="old news", + geoh5=geoh5, + data_object=components.survey, + gz_channel=gz_channel, + gz_uncertainty=gz_uncerts, + mesh=components.mesh, + starting_model=components.model, + topography_object=components.topography, + ) params_uijson_path = tmp_path / "from_params.ui.json" opts.write_ui_json(params_uijson_path) @@ -326,21 +337,25 @@ def test_legacy_uijson(tmp_path: Path): ) work_path.mkdir(parents=True) - geoh5, mesh, model, survey, topo = setup_inversion_workspace( - work_path, - background=1.0, - anomaly=2.0, - n_electrodes=10, - n_lines=3, - inversion_type=inversion_type, + opts = SyntheticsComponentsOptions( + method=inversion_type, + survey=SurveyOptions( + n_stations=10, + n_lines=3, + ), + mesh=MeshOptions(), + model=ModelOptions( + background=1.0, + anomaly=2.0, + ), ) - - with geoh5.open(mode="r+"): + with Workspace.create(work_path / "inversion_test.ui.geoh5") as geoh5: + components = SyntheticsComponents(geoh5, options=opts) ifile.data["geoh5"] = geoh5 - ifile.data["mesh"] = mesh - ifile.data["starting_model"] = model - ifile.data["data_object"] = survey - ifile.data["topography_object"] = topo + ifile.data["mesh"] = components.mesh + ifile.data["starting_model"] = components.model + ifile.data["data_object"] = components.survey + ifile.data["topography_object"] = components.topography # Test deprecated name ifile.data["coolingFactor"] = 4.0 @@ -350,19 +365,19 @@ def test_legacy_uijson(tmp_path: Path): ifile.data["line_object"] = line_id if not forward: - n_vals = survey.n_vertices + n_vals = components.survey.n_vertices if ( "direct current" in inversion_type or "induced polarization" in inversion_type ): - n_vals = survey.n_cells + n_vals = components.survey.n_cells - channels = getattr(survey, "channels", [1]) + channels = getattr(components.survey, "channels", [1]) data = [] for channel in channels: data.append( - survey.add_data( + components.survey.add_data( { CHANNEL_NAME[inversion_type] + f"[{channel}]": { "values": np.ones(n_vals) @@ -372,7 +387,7 @@ def test_legacy_uijson(tmp_path: Path): ) if len(data) > 1: - channel = survey.add_data_to_group(data, "Group") + channel = components.survey.add_data_to_group(data, "Group") else: channel = data[0] diff --git a/tests/utils/__init__.py b/tests/utils/__init__.py new file mode 100644 index 00000000..4d06f672 --- /dev/null +++ b/tests/utils/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/tests/utils/targets.py b/tests/utils/targets.py new file mode 100644 index 00000000..ca12a074 --- /dev/null +++ b/tests/utils/targets.py @@ -0,0 +1,72 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +import warnings +from uuid import UUID + +import numpy as np +from geoapps_utils.utils.conversions import string_to_numeric +from geoh5py import Workspace + + +def get_inversion_output(h5file: str | Workspace, inversion_group: str | UUID): + """ + Recover inversion iterations from a ContainerGroup comments. + """ + if isinstance(h5file, Workspace): + workspace = h5file + else: + workspace = Workspace(h5file) + + try: + group = workspace.get_entity(inversion_group)[0] + except IndexError as exc: + raise IndexError( + f"BaseInversion group {inversion_group} could not be found in the target geoh5 {h5file}" + ) from exc + + outfile = group.get_entity("SimPEG.out")[0] + out = list(outfile.file_bytes.decode("utf-8").replace("\r", "").split("\n"))[:-1] + cols = out.pop(0).split(" ") + out = [[string_to_numeric(k) for k in elem.split(" ")] for elem in out] + out = dict(zip(cols, list(map(list, zip(*out, strict=True))), strict=True)) + + return out + + +def check_target(output: dict, target: dict, tolerance=0.05): + """ + Check inversion output metrics against hard-valued target. + :param output: Dictionary containing keys for 'data', 'phi_d' and 'phi_m'. + :param target: Dictionary containing keys for 'data_norm', 'phi_d' and 'phi_m'.\ + :param tolerance: Tolerance between output and target measured as: |a-b|/b + """ + print( + f"Output: 'data_norm': {np.linalg.norm(output['data'])}, 'phi_d': {output['phi_d'][1]}, 'phi_m': {output['phi_m'][1]}" + ) + print(f"Target: {target}") + + if any(np.isnan(output["data"])): + warnings.warn( + "Skipping data norm comparison due to nan (used to bypass lone faulty test run in GH actions)." + ) + else: + np.testing.assert_array_less( + np.abs(np.linalg.norm(output["data"]) - target["data_norm"]) + / target["data_norm"], + tolerance, + ) + + np.testing.assert_array_less( + np.abs(output["phi_m"][1] - target["phi_m"]) / target["phi_m"], tolerance + ) + np.testing.assert_array_less( + np.abs(output["phi_d"][1] - target["phi_d"]) / target["phi_d"], tolerance + )