diff --git a/docs/source/conf.py b/docs/source/conf.py index 7df3b1a..567807f 100644 --- a/docs/source/conf.py +++ b/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()}. """ diff --git a/plate_simulation/driver.py b/plate_simulation/driver.py index 139355d..06c75b3 100644 --- a/plate_simulation/driver.py +++ b/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): diff --git a/plate_simulation/models/series.py b/plate_simulation/models/series.py index 5d642b2..ae2d336 100644 --- a/plate_simulation/models/series.py +++ b/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): diff --git a/tests/models/series_test.py b/tests/models/series_test.py index 58e06bb..dc77cbd 100644 --- a/tests/models/series_test.py +++ b/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,7 +104,6 @@ def test_scenario(tmp_path): mesh=octree, background=0.0, history=[overburden, lithology], - name="model", ) scenario = Scenario( @@ -113,20 +111,19 @@ def test_scenario(tmp_path): 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) diff --git a/tests/runtest/gravity_test.py b/tests/runtest/gravity_test.py index a270435..551c72f 100644 --- a/tests/runtest/gravity_test.py +++ b/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