Merge pull request #47 from MiraGeoscience/GEOPY-1882

GEOPY-1882: Reciprocal applied to Magnetic and density

docs/source/conf.py

@@ -7,8 +7,9 @@
 #  (see LICENSE file at the root of this source code package).                         '
 # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
 
-from importlib.metadata import version
 from datetime import datetime
+from importlib.metadata import version
+
 
 # Configuration file for the Sphinx documentation builder.
 #
@@ -51,13 +52,15 @@
 
 html_theme = "alabaster"
 html_theme_options = {
-    'description': f"version {release}",
+    "description": f"version {release}",
 }
 html_static_path = ["_static"]
 
+
 def get_copyright_notice():
     return f"Copyright {datetime.now().strftime(project_copyright)}"
 
+
 rst_epilog = f"""
 .. |copyright_notice| replace:: {get_copyright_notice()}.
 """

plate_simulation/driver.py

@@ -237,10 +237,19 @@ def make_model(self) -> FloatData:
             mesh=self.mesh,
             background=self.params.model.background,
             history=[dikes, overburden, erosion],
-            name=self.params.model.name,
         )
 
-        return scenario.geologize()
+        geology = scenario.geologize()
+
+        if self.simulation_parameters.physical_property == "conductivity":
+            geology **= -1.0
+
+        with fetch_active_workspace(self.params.geoh5, mode="r+"):
+            model: FloatData = self.mesh.add_data(  # type: ignore
+                {self.params.model.name: {"values": geology}}
+            )
+
+        return model
 
     @staticmethod
     def start(ifile: str | Path | InputFile):

plate_simulation/models/series.py

@@ -14,7 +14,6 @@
 
 import numpy as np
 from geoh5py import Workspace
-from geoh5py.data import FloatData
 from geoh5py.objects import Octree
 from geoh5py.shared.utils import fetch_active_workspace
 
@@ -96,7 +95,6 @@ class Scenario(Series):
     :param background: Initial value that will fill any areas of the model
         not covered by event realizations.
     :param history: Geological events that form the model.
-    :param name: Name of the model that will be added to the mesh object.
     """
 
     def __init__(
@@ -106,14 +104,12 @@ def __init__(
         mesh: Octree,
         background: float,
         history: Sequence[Event | Series],
-        name: str = "model",
     ):
         super().__init__(history)
         self.workspace = workspace
         self.mesh = mesh
         self.background = background
         self.history = Geology(history)
-        self.name = name
 
     @property
     def mesh(self) -> Octree:
@@ -126,20 +122,15 @@ def mesh(self, val: Octree):
             raise ValueError("Mesh must have n_cells.")
         self._mesh = val
 
-    def geologize(self) -> FloatData:
+    def geologize(self) -> np.ndarray:
         """Realize the geological events in the scenario"""
         with fetch_active_workspace(self.workspace, mode="r+"):
             if self.mesh.n_cells is None:
                 raise ValueError("Mesh must have n_cells.")
-            geology = (
-                super().realize(self.mesh, np.ones(self.mesh.n_cells) * self.background)
-                ** -1.0
+            geology = super().realize(
+                self.mesh, np.ones(self.mesh.n_cells) * self.background
             )
-            model: FloatData = self.mesh.add_data(  # type: ignore
-                {self.name: {"values": geology}}
-            )
-
-        return model
+        return geology
 
 
 class GeologyViolationError(Exception):

tests/models/series_test.py

@@ -93,7 +93,6 @@ def test_scenario(tmp_path):
                 mesh=octree,
                 background=0.0,
                 history=[lithology, erosion, overburden],
-                name="model",
             )
 
         with pytest.raises(
@@ -105,28 +104,26 @@ def test_scenario(tmp_path):
                 mesh=octree,
                 background=0.0,
                 history=[overburden, lithology],
-                name="model",
             )
 
         scenario = Scenario(
             workspace=ws,
             mesh=octree,
             background=100.0,
             history=[lithology, overburden, erosion],
-            name="model",
         )
         model = scenario.geologize()
-        assert model.values is not None
+        assert model is not None
 
         ind = octree.centroids[:, 2] > 0.0
-        assert all(np.isnan(model.values[ind]))
+        assert all(np.isnan(model[ind]))
         ind = (octree.centroids[:, 2] < 0.0) & (octree.centroids[:, 2] > -1.0)
-        assert all(model.values[ind] == 0.1)
+        assert all(model[ind] == 10.0)
         ind = (octree.centroids[:, 2] < -1.0) & (octree.centroids[:, 2] > -2.0)
-        assert all(model.values[ind] == 0.01)
+        assert all(model[ind] == 100.0)
         ind = (octree.centroids[:, 2] < -2.0) & (octree.centroids[:, 2] > -5.0)
-        assert all(model.values[ind] == 1.0)
+        assert all(model[ind] == 1.0)
         ind = (octree.centroids[:, 2] < -5.0) & (octree.centroids[:, 2] > -10.0)
-        assert all(model.values[ind] == 0.5)
+        assert all(model[ind] == 2.0)
         ind = octree.centroids[:, 2] < -10.0
-        np.testing.assert_allclose(model.values[ind], 1 / 3.0)
+        np.testing.assert_allclose(model[ind], 3.0)

tests/runtest/gravity_test.py

@@ -9,6 +9,7 @@
 
 from copy import deepcopy
 
+import numpy as np
 from geoh5py import Workspace
 from geoh5py.groups import SimPEGGroup
 from simpeg_drivers.constants import default_ui_json
@@ -35,11 +36,11 @@ def test_gravity_plate_simulation(tmp_path):
             max_distance=200.0,
         )
 
-        overburden_params = OverburdenParams(thickness=50.0, overburden=5.0)
+        overburden_params = OverburdenParams(thickness=50.0, overburden=0.2)
 
         plate_params = PlateParams(
             name="plate",
-            plate=2.0,
+            plate=0.5,
             elevation=-250.0,
             width=100.0,
             strike_length=100.0,
@@ -51,7 +52,7 @@ def test_gravity_plate_simulation(tmp_path):
 
         model_params = ModelParams(
             name="density",
-            background=1000.0,
+            background=0.0,
             overburden=overburden_params,
             plate=plate_params,
         )
@@ -77,3 +78,5 @@ def test_gravity_plate_simulation(tmp_path):
         )
         driver = PlateSimulationDriver(params)
         driver.run()
+
+        assert np.nanmax(driver.model.values) == 0.5