diff --git a/data/example_config.toml b/data/example_config.toml deleted file mode 100644 index 7e122077f..000000000 --- a/data/example_config.toml +++ /dev/null @@ -1,11 +0,0 @@ -[core] -modules = ["soil"] - -[core.grid] -cell_nx = 10 -cell_ny = 10 - -[core.timing] -start_date = "2020-01-01" -run_length = "50 years" -update_interval = "1 day" diff --git a/docs/source/conf.py b/docs/source/conf.py index ca4d34325..ff919dbc4 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -67,6 +67,7 @@ "sphinx.ext.autosummary", # "sphinx.ext.autosectionlabel", # Generates hard to trace exception "sphinxcontrib.bibtex", + "sphinxcontrib.mermaid", "myst_nb", "sphinx_rtd_theme", ] diff --git a/docs/source/data_recipes/ERA5_preprocessing_example.md b/docs/source/data_recipes/ERA5_preprocessing_example.md index fb22a1f84..61cd06d23 100644 --- a/docs/source/data_recipes/ERA5_preprocessing_example.md +++ b/docs/source/data_recipes/ERA5_preprocessing_example.md @@ -163,7 +163,7 @@ The dummy model iterates over time indices rather than real datetime. Therefore, a `time_index` coordinate to the dataset: ```{code-cell} ipython3 -dataset_xy_timeindex = dataset_xy_dummy.rename_dims({'time':'time_index'}).assign_coords({'time_index':np.arange(0,24,1)}) +dataset_xy_timeindex = dataset_xy_dummy.assign_coords({'time_index':np.arange(0,24,1)}) dataset_xy_timeindex.coords ``` diff --git a/docs/source/virtual_rainforest/core/data.md b/docs/source/virtual_rainforest/core/data.md index 39297f75a..7c9c3e1d5 100644 --- a/docs/source/virtual_rainforest/core/data.md +++ b/docs/source/virtual_rainforest/core/data.md @@ -194,6 +194,11 @@ configuration section. This can be used to automatically load multiple datasets a Data object. The configuration file is TOML formatted and should contain an entry like the example below for each variable to be loaded. +**NOTE**: At the moment, +`core.data.variable` tags cannot be used across multiple toml config files without +causing `ConfigurationError: Duplicated entries in config files: core.data.variable` to +be raised. This means that all variables need to be combined in one `config` file. + ```toml [[core.data.variable]] file="'../../data/xy_dim.nc'" diff --git a/docs/source/virtual_rainforest/main_simulation.md b/docs/source/virtual_rainforest/main_simulation.md index bd46fad67..03f073c0d 100644 --- a/docs/source/virtual_rainforest/main_simulation.md +++ b/docs/source/virtual_rainforest/main_simulation.md @@ -1,71 +1,159 @@ # Virtual Rainforest simulation flow -This document describes the main simulation flow of the Virtual Rainforest model. This -flow will now be set out step by step. +This document describes the main simulation flow of the Virtual Rainforest model. The +main stages are: + +* setup of the **simulation core** that provides shared resources and functions +* setup of the individual **science models** that simulate the behaviour of different +components of the Virtual Rainforest, and +* iteration of the simulation over the configured timescale. + +```{mermaid} +flowchart TD + A[vr_run] --> B + B --> F + C --> F + D --> F + D --> G + E --> F + H --> H2 + subgraph Core + direction LR + B[Load configuration] --> C[Create grid] + C --> D[Load data] + D --> E[Validate timescale] + end + subgraph Setup science models + direction LR + F[Model configuration] --> G[Model setup] + G --> H[Model spinup] + end + subgraph Simulation + H2[Save initial state] --> I[Start time] + I --> J[Update interval] + J --> K[Update science models] + K --> J + K --> L[End time] + L --> M[Save final state] + end +``` + +## Core setup + +The first stage in the simulation is the configuration and initialisation of the core +resources and functionality. + +### Loading configuration files + +First, a set of user-provided configuration files in `TOML` format for a simulation are +loaded. These files are then validated to ensure: + +* that they are valid `TOML`, +* and that all the required settings are present and not duplicated. + +Some settings will be filled automatically from defaults settings and so can be omitted, +but validation will fail if mandatory settings are omitted. Further details can be found +in the [configuration documentation](./core/config.md). + +### Grid creation + +Next, the spatial structure of the simulation is configured as a [`Grid` +object](./core/grid.md) that defines the area, coordinate system and geometry of the +individual cells that will be used in the simulation. + +### Loading and validation of input data + +All of the data required to initialise and run the simulation is then loaded into an +internal [`Data` object](./core/data.md). The model configuration sets the locations of +files containing required variables and this configuration is passed into the +{meth}`~virtual_rainforest.core.data.Data.load_data_config` method, which ensures that: + +* the input files are valid and can be read, and +* that the data in files is congruent with the rest of the configuration, such as + checking the dimensionality and shape of [core axes](./core/axes.md) like the spatial + grid. + +### Simulation timescale + +The simulation runs between two dates with an update interval at which each science +model is recalculated. These values are defined in the `core` configuration and are +now validated to ensure that the start date, end date and update interval are sensible. + +```{note} +The simulation uses 12 equal length months (30.4375 days) and equal length years (365.25 +days), ignoring leap years. +``` + +## Science models -## Loading configuration files +The Virtual Rainforest is implemented as model objects, each of which is responsible for +simulating a particular aspect of the rainforest ecosystem. The models used for the +specific simulation run can be set in the configuration and will typically include the +four standard models: -As a first step, configuration files are loaded in and validated. These files are of -`toml` format and are found based on paths provided by the user. When these files are -loaded in they are validated to ensure that they are valid `toml`, and that all the -required settings are present and not duplicated. The Virtual Rainforest code provides -default settings for some configuration options and these will be used automatically if -that setting is not found in the configuration, but if no default is set then the -validation will fail. Further details can be found in the [configuration -documentation](./core/config.md). +* the [`AbioticModel`](../api/abiotic.md) or + [`AbioticSimpleModel`](../api/abiotic_simple.md), +* the `AnimalModel`, +* the `PlantModel` and the +* [`SoilModel`](../api/soil.md) -## Grid creation +but this can be [extended to include new models](../development/defining_new_models.md) +or adopt different combinations of models. -The next step is creating a grid, which defines the spatial structure of the simulation: -the area, coordinate system and geometry of the individual cells that will be used in -the simulation. The [grid documentation](./core/grid.md) describes this object further. -This grid is then used to create an empty `Data` object of the correct size to describe -the simulation. +Once a list of models to configure has been extracted from the configuration, all +science models run through a set of steps to prepare for the simulation to start. Each +step is represented using a set of standard model methods that are run in the following +order. -## Loading and validation of input data +### Model configuration -This `Data` object now needs to be populated with data. The data used in the simulation -is stored in a set of data files which are specified in the configuration. This data is -then loaded using the `data.load_data_config` method, which has built in validation -procedures to ensure that the loaded data can be mapped onto the spatial grid and also -other core dimensions for the configured simulation. Further details can be found in the -[data system](./core/data.md) and [core axes](./core/axes.md) documentation. +The loaded configuration should include the configuration details for each individual +science model. These are now used to initialise each requested model using the +{meth}`~virtual_rainforest.core.base_model.BaseModel.from_config` method defined +for each model. This method checks that the configuration is valid for the science +model. -## Configuration of specific models +### Model setup -The Virtual Rainforest is implemented as model objects, each of which is responsible for -simulating a particular aspect of the rainforest ecosystem. The models used for the -specific simulation run can be set in the configuration. This will typically be the four -standard models ([`AbioticModel`](../api/abiotic.md), `AnimalModel`, `PlantModel` and -[`SoilModel`](../api/soil.md)), but this can be extended to include new models or -different combinations of models. For more information about implementing new models, -see [this page](../development/defining_new_models.md) about the required module -structure and the model API. Once a list of models to configure has been extracted from -the configuration, they are then all configured. +Some models require an additional setup step to calculate values for internal variables +from the initial loaded data or to set up further structures within the model, such as +representations of plant or animal communities. Each model will run the +{meth}`~virtual_rainforest.core.base_model.BaseModel.setup` method defined for the +specific model. In simple science models, this method may not actually need to do +anything. + +### Model spinup + +Some models may then require a spin up step to allow initial variables to reach an +equilibrium before running the main simulation. Again, each model will run the +{meth}`~virtual_rainforest.core.base_model.BaseModel.spinup` method defined for the +specific model, and again this may not need to do anything for simple models. + +### Model update + +At this point, the model instance is now ready for simulation. The +{meth}`~virtual_rainforest.core.base_model.BaseModel.update` method for each science +model is run as part of the simulation process described below. -## Extracting simulation timing details +## Simulation process -The configuration contains a start time, an end time and a time interval for the update -of all models. These details are extracted from configuration, with a check performed -to ensure that the simulation will update at least once between the start and end time -of the simulation. It is important to note that because months and years are of -inconsistent length they are currently averaged over. This means that 1 month becomes -30.4375 days, and 1 year becomes 365.25 days. +Now that the simulation core and science models have been configure and initialised, +along with any setup or spinup steps, the simulation itself starts. -## Saving the initial state +### Saving the initial state The `data` object has now been populated with all of the configured data required to run -the model. It can now optionally be saved to a file set in the configuration to generate -a single data file of the initial model state. +the model. The simulation configuration can optionally provide a filepath that will be +used to output a single data file of the initial simulation state. -## Simulating over time +### Simulation -The previously extracted timing details are used to setup a timing loop, which runs from -the start time to the end time with a time step set by the update interval. At each +The science models are now iterated over the configured simulation timescale, running +from the start time to the end time with a time step set by the update interval. At each step all models are updated. -## Saving the final state +### Saving the final state -After the full simulation loop has been completed the final simulation state can be -saved. Whether to save this state and where to save it to are again configuration -options, but in this case the default is for the final state to be saved. +After the full simulation loop has been completed, the final simulation state held in +the `Data` object can be optionally be saved to a path provided in the configuration, +defaulting to saving the data. diff --git a/docs/source/virtual_rainforest/usage.md b/docs/source/virtual_rainforest/usage.md index c872affda..69ce13d7d 100644 --- a/docs/source/virtual_rainforest/usage.md +++ b/docs/source/virtual_rainforest/usage.md @@ -21,10 +21,19 @@ access to all relevant dependencies. ## Running an example Virtual Rainforest simulation -An example configuration for the Virtual Rainforest model is provided in -`data/example_config.toml`. A simulation using this configuration can be run through the -command line by calling: +An example configuration file and the dummy data required to run the `abiotic_simple` +and `soil` models are stored in the Imperial College Research Data Store (RDS). To +access this data you must have been added to the the Virtual Rainforest RDS project. +Details of how to access RDS can be found +[here](https://wiki.imperial.ac.uk/pages/viewpage.action?spaceKey=HPC&title=Research+Data+Store) +(Imperial College login required). + +Once you are connected to RDS a simulation using this configuration and data set can be +run through the command line by calling: ```shell -vr_run data/example_config.toml +vr_run /Volumes/virtual_rainforest/live/dummy_data/ ``` + +**Note:** The directory path above is specific to macOS and will be different for +Windows and Linux machines. diff --git a/poetry.lock b/poetry.lock index 016ccc3d5..ded19a739 100644 --- a/poetry.lock +++ b/poetry.lock @@ -1577,6 +1577,14 @@ python-versions = ">=3.5" [package.extras] test = ["flake8", "mypy", "pytest"] +[[package]] +name = "sphinxcontrib-mermaid" +version = "0.9.2" +description = "Mermaid diagrams in yours Sphinx powered docs" +category = "dev" +optional = false +python-versions = ">=3.7" + [[package]] name = "sphinxcontrib-qthelp" version = "1.0.3" @@ -1797,7 +1805,7 @@ testing = ["flake8 (<5)", "func-timeout", "jaraco.functools", "jaraco.itertools" [metadata] lock-version = "1.1" python-versions = ">=3.9,<3.12" -content-hash = "f88a906d7dbec3b88d3d669b340642db285cb228a0c6e0d897bf0e7b5dfc81f3" +content-hash = "df797fc1106331773d681494de5e21b73014f846a52844909e520dc1a7826570" [metadata.files] alabaster = [ @@ -2145,6 +2153,7 @@ debugpy = [ {file = "debugpy-1.6.6-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:9b5d1b13d7c7bf5d7cf700e33c0b8ddb7baf030fcf502f76fc061ddd9405d16c"}, {file = "debugpy-1.6.6-cp38-cp38-win32.whl", hash = "sha256:70ab53918fd907a3ade01909b3ed783287ede362c80c75f41e79596d5ccacd32"}, {file = "debugpy-1.6.6-cp38-cp38-win_amd64.whl", hash = "sha256:c05349890804d846eca32ce0623ab66c06f8800db881af7a876dc073ac1c2225"}, + {file = "debugpy-1.6.6-cp39-cp39-macosx_11_0_x86_64.whl", hash = "sha256:11a0f3a106f69901e4a9a5683ce943a7a5605696024134b522aa1bfda25b5fec"}, {file = "debugpy-1.6.6-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a771739902b1ae22a120dbbb6bd91b2cae6696c0e318b5007c5348519a4211c6"}, {file = "debugpy-1.6.6-cp39-cp39-win32.whl", hash = "sha256:549ae0cb2d34fc09d1675f9b01942499751d174381b6082279cf19cdb3c47cbe"}, {file = "debugpy-1.6.6-cp39-cp39-win_amd64.whl", hash = "sha256:de4a045fbf388e120bb6ec66501458d3134f4729faed26ff95de52a754abddb1"}, @@ -2760,6 +2769,8 @@ psutil = [ {file = "psutil-5.9.4-cp36-abi3-macosx_10_9_x86_64.whl", hash = "sha256:54d5b184728298f2ca8567bf83c422b706200bcbbfafdc06718264f9393cfeb7"}, {file = "psutil-5.9.4-cp36-abi3-manylinux_2_12_i686.manylinux2010_i686.manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:16653106f3b59386ffe10e0bad3bb6299e169d5327d3f187614b1cb8f24cf2e1"}, {file = "psutil-5.9.4-cp36-abi3-manylinux_2_12_x86_64.manylinux2010_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:54c0d3d8e0078b7666984e11b12b88af2db11d11249a8ac8920dd5ef68a66e08"}, + {file = "psutil-5.9.4-cp36-abi3-win32.whl", hash = "sha256:149555f59a69b33f056ba1c4eb22bb7bf24332ce631c44a319cec09f876aaeff"}, + {file = "psutil-5.9.4-cp36-abi3-win_amd64.whl", hash = "sha256:fd8522436a6ada7b4aad6638662966de0d61d241cb821239b2ae7013d41a43d4"}, {file = "psutil-5.9.4-cp38-abi3-macosx_11_0_arm64.whl", hash = "sha256:6001c809253a29599bc0dfd5179d9f8a5779f9dffea1da0f13c53ee568115e1e"}, {file = "psutil-5.9.4.tar.gz", hash = "sha256:3d7f9739eb435d4b1338944abe23f49584bde5395f27487d2ee25ad9a8774a62"}, ] @@ -3164,6 +3175,10 @@ sphinxcontrib-jsmath = [ {file = "sphinxcontrib-jsmath-1.0.1.tar.gz", hash = "sha256:a9925e4a4587247ed2191a22df5f6970656cb8ca2bd6284309578f2153e0c4b8"}, {file = "sphinxcontrib_jsmath-1.0.1-py2.py3-none-any.whl", hash = "sha256:2ec2eaebfb78f3f2078e73666b1415417a116cc848b72e5172e596c871103178"}, ] +sphinxcontrib-mermaid = [ + {file = "sphinxcontrib-mermaid-0.9.2.tar.gz", hash = "sha256:252ef13dd23164b28f16d8b0205cf184b9d8e2b714a302274d9f59eb708e77af"}, + {file = "sphinxcontrib_mermaid-0.9.2-py3-none-any.whl", hash = "sha256:6795a72037ca55e65663d2a2c1a043d636dc3d30d418e56dd6087d1459d98a5d"}, +] sphinxcontrib-qthelp = [ {file = "sphinxcontrib-qthelp-1.0.3.tar.gz", hash = "sha256:4c33767ee058b70dba89a6fc5c1892c0d57a54be67ddd3e7875a18d14cba5a72"}, {file = "sphinxcontrib_qthelp-1.0.3-py2.py3-none-any.whl", hash = "sha256:bd9fc24bcb748a8d51fd4ecaade681350aa63009a347a8c14e637895444dfab6"}, @@ -3265,7 +3280,6 @@ virtualenv = [ {file = "virtualenv-20.19.0.tar.gz", hash = "sha256:37a640ba82ed40b226599c522d411e4be5edb339a0c0de030c0dc7b646d61590"}, ] wcwidth = [ - {file = "wcwidth-0.2.6-py2.py3-none-any.whl", hash = "sha256:795b138f6875577cd91bba52baf9e445cd5118fd32723b460e30a0af30ea230e"}, {file = "wcwidth-0.2.6.tar.gz", hash = "sha256:a5220780a404dbe3353789870978e472cfe477761f06ee55077256e509b156d0"}, ] wheel = [ diff --git a/pyproject.toml b/pyproject.toml index 6f074b94c..909beca57 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -59,6 +59,7 @@ myst-nb = "^0.16.0" sphinx-rtd-theme = "^1.0.0" autodocsumm = "^0.2.8" pydocstyle = "^6.1.1" +sphinxcontrib-mermaid = "^0.9.2" [build-system] requires = ["poetry-core>=1.2.0"] diff --git a/tests/conftest.py b/tests/conftest.py index 1896a934e..207dd035a 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -239,7 +239,7 @@ def dummy_climate_data(layer_roles_fixture): data["air_temperature_ref"] = DataArray( np.full((3, 3), 30), - dims=["cell_id", "time"], + dims=["cell_id", "time_index"], ) data["mean_annual_temperature"] = DataArray( np.full((3), 20), @@ -247,19 +247,19 @@ def dummy_climate_data(layer_roles_fixture): ) data["relative_humidity_ref"] = DataArray( np.full((3, 3), 90), - dims=["cell_id", "time"], + dims=["cell_id", "time_index"], ) data["vapour_pressure_deficit_ref"] = DataArray( np.full((3, 3), 0.14), - dims=["cell_id", "time"], + dims=["cell_id", "time_index"], ) data["atmospheric_pressure_ref"] = DataArray( np.full((3, 3), 1.5), - dims=["cell_id", "time"], + dims=["cell_id", "time_index"], ) data["atmospheric_co2_ref"] = DataArray( np.full((3, 3), 400), - dims=["cell_id", "time"], + dims=["cell_id", "time_index"], ) leaf_area_index = np.repeat(a=[np.nan, 1.0, np.nan], repeats=[1, 3, 11]) @@ -290,7 +290,7 @@ def dummy_climate_data(layer_roles_fixture): ) data["precipitation"] = DataArray( - [[20, 30, 200], [20, 30, 200], [20, 30, 200]], dims=["time", "cell_id"] + [[20, 30, 200], [20, 30, 200], [20, 30, 200]], dims=["time_index", "cell_id"] ) data["soil_moisture"] = xr.concat( [ diff --git a/tests/core/test_utils.py b/tests/core/test_utils.py index 38a964add..3e4f2e40a 100644 --- a/tests/core/test_utils.py +++ b/tests/core/test_utils.py @@ -1,12 +1,13 @@ """Testing the utility functions.""" -from logging import CRITICAL +from contextlib import nullcontext as does_not_raise +from logging import CRITICAL, ERROR from pathlib import Path import pytest from tests.conftest import log_check -from virtual_rainforest.core.exceptions import ConfigurationError +from virtual_rainforest.core.exceptions import ConfigurationError, InitialisationError @pytest.mark.parametrize( @@ -58,3 +59,75 @@ def test_check_outfile(caplog, mocker, out_path, expected_log_entries): check_outfile(Path(out_path)) log_check(caplog, expected_log_entries) + + +@pytest.mark.parametrize( + "soil_layers, canopy_layers, raises, exp_log", + [ + ( + 2, + 10, + does_not_raise(), + (), + ), + ( + -2, + 10, + pytest.raises(InitialisationError), + ( + ( + ERROR, + "There has to be at least one soil layer in the Virtual " + "Rainforest!", + ), + ), + ), + ( + 2, + -3, + pytest.raises(InitialisationError), + ( + ( + ERROR, + "There has to be at least one canopy layer in the Virtual " + "Rainforest!", + ), + ), + ), + ( + 2.5, + 10, + pytest.raises(InitialisationError), + ( + ( + ERROR, + "The number of soil layers must be an integer!", + ), + ), + ), + ( + 2, + 3.4, + pytest.raises(InitialisationError), + ( + ( + ERROR, + "The number of canopy layers must be an integer!", + ), + ), + ), + ], +) +def test_set_layer_roles(soil_layers, canopy_layers, raises, caplog, exp_log): + """Test correct order of layers.""" + from virtual_rainforest.core.utils import set_layer_roles + + with raises: + result = set_layer_roles(canopy_layers, soil_layers) + + assert result == ( + ["above"] + ["canopy"] * 10 + ["subcanopy"] + ["surface"] + ["soil"] * 2 + ) + + # Final check that expected logging entries are produced + log_check(caplog, exp_log) diff --git a/tests/models/abiotic_simple/test_abiotic_simple_model.py b/tests/models/abiotic_simple/test_abiotic_simple_model.py index bcdeb624a..d21819d0b 100644 --- a/tests/models/abiotic_simple/test_abiotic_simple_model.py +++ b/tests/models/abiotic_simple/test_abiotic_simple_model.py @@ -16,24 +16,11 @@ ) -def test_set_layer_roles(): - """Test correct order of layers.""" - from virtual_rainforest.models.abiotic_simple.abiotic_simple_model import ( - set_layer_roles, - ) - - assert set_layer_roles(10, 2) == ( - ["above"] + ["canopy"] * 10 + ["subcanopy"] + ["surface"] + ["soil"] * 2 - ) - - @pytest.mark.parametrize( - "soil_layers,canopy_layers,ini_soil_moisture,raises,expected_log_entries", + "ini_soil_moisture,raises,expected_log_entries", [ ( - 2, - 10, - 50.0, + 0.5, does_not_raise(), ( ( @@ -74,68 +61,16 @@ def test_set_layer_roles(): ), ), ( - -2, - 10, - 50.0, - pytest.raises(InitialisationError), - ( - ( - ERROR, - "There has to be at least one soil layer in the abiotic model!", - ), - ), - ), - ( - 2, - -3, - 50.0, - pytest.raises(InitialisationError), - ( - ( - ERROR, - "There has to be at least one canopy layer in the abiotic model!", - ), - ), - ), - ( - 2.5, - 10, - 50.0, - pytest.raises(InitialisationError), - ( - ( - ERROR, - "The number of soil layers must be an integer!", - ), - ), - ), - ( - 2, - 3.4, - 50.0, - pytest.raises(InitialisationError), - ( - ( - ERROR, - "The number of canopy layers must be an integer!", - ), - ), - ), - ( - 2, - 10, - -50.0, + -0.5, pytest.raises(InitialisationError), ( ( ERROR, - "The initial soil moisture has to be between 0 and 100!", + "The initial soil moisture has to be between 0 and 1!", ), ), ), ( - 2, - 10, DataArray([50, 30, 20]), pytest.raises(InitialisationError), ( @@ -150,12 +85,12 @@ def test_set_layer_roles(): def test_abiotic_simple_model_initialization( caplog, dummy_climate_data, - soil_layers, - canopy_layers, ini_soil_moisture, raises, expected_log_entries, layer_roles_fixture, + soil_layers=2, + canopy_layers=10, ): """Test `AbioticSimpleModel` initialization.""" @@ -202,12 +137,14 @@ def test_abiotic_simple_model_initialization( "timing": { "start_date": "2020-01-01", "update_interval": "1 week", - } + }, + "layers": { + "soil_layers": 2, + "canopy_layers": 10, + }, }, "abiotic_simple": { - "soil_layers": 2, - "canopy_layers": 10, - "initial_soil_moisture": 50.0, + "initial_soil_moisture": 0.5, }, }, pint.Quantity("1 week"), @@ -291,12 +228,14 @@ def test_generate_abiotic_simple_model( "timing": { "start_date": "2020-01-01", "update_interval": "1 week", - } + }, + "layers": { + "soil_layers": 2, + "canopy_layers": 10, + }, }, "abiotic_simple": { - "soil_layers": 2, - "canopy_layers": 10, - "initial_soil_moisture": 50.0, + "initial_soil_moisture": 0.5, }, }, pint.Quantity("1 week"), @@ -322,7 +261,7 @@ def test_setup( model.setup() - soil_moisture_values = np.repeat(a=[np.nan, 50], repeats=[13, 2]) + soil_moisture_values = np.repeat(a=[np.nan, 0.5], repeats=[13, 2]) xr.testing.assert_allclose( dummy_climate_data["soil_moisture"], @@ -342,7 +281,7 @@ def test_setup( dummy_climate_data["vapour_pressure_deficit_ref"], DataArray( np.full((3, 3), 0.141727), - dims=["cell_id", "time"], + dims=["cell_id", "time_index"], coords={ "cell_id": [0, 1, 2], }, diff --git a/tests/models/abiotic_simple/test_simple_regression.py b/tests/models/abiotic_simple/test_simple_regression.py index 0abf0e147..c0e404799 100644 --- a/tests/models/abiotic_simple/test_simple_regression.py +++ b/tests/models/abiotic_simple/test_simple_regression.py @@ -19,7 +19,7 @@ def test_log_interpolation(dummy_climate_data, layer_roles_fixture): # temperature result = log_interpolation( data=data, - reference_data=data["air_temperature_ref"].isel(time=0), + reference_data=data["air_temperature_ref"].isel(time_index=0), leaf_area_index_sum=leaf_area_index_sum, layer_roles=layer_roles_fixture, layer_heights=data["layer_heights"], @@ -62,7 +62,7 @@ def test_log_interpolation(dummy_climate_data, layer_roles_fixture): # relative humidity result_hum = log_interpolation( data=data, - reference_data=data["relative_humidity_ref"].isel(time=0), + reference_data=data["relative_humidity_ref"].isel(time_index=0), leaf_area_index_sum=leaf_area_index_sum, layer_roles=layer_roles_fixture, layer_heights=data["layer_heights"], @@ -112,7 +112,7 @@ def test_calculate_saturation_vapour_pressure(dummy_climate_data): data = dummy_climate_data result = calculate_saturation_vapour_pressure( - data["air_temperature_ref"].isel(time=0) + data["air_temperature_ref"].isel(time_index=0) ) exp_output = DataArray( @@ -324,7 +324,7 @@ def test_calculate_soil_moisture(dummy_climate_data, layer_roles_fixture): ) data = dummy_climate_data - precipitation_surface = data["precipitation"].isel(time=0) * ( + precipitation_surface = data["precipitation"].isel(time_index=0) * ( 1 - 0.1 * data["leaf_area_index"].sum(dim="layers") ) diff --git a/tests/models/soil/test_soil_model.py b/tests/models/soil/test_soil_model.py index aae8a497d..2dc5a3559 100644 --- a/tests/models/soil/test_soil_model.py +++ b/tests/models/soil/test_soil_model.py @@ -22,7 +22,10 @@ def soil_model_fixture(dummy_carbon_data): from virtual_rainforest.models.soil.soil_model import SoilModel config = { - "core": {"timing": {"start_date": "2020-01-01", "update_interval": "12 hours"}}, + "core": { + "timing": {"start_date": "2020-01-01", "update_interval": "12 hours"}, + "layers": {"soil_layers": 2, "canopy_layers": 10}, + }, } return SoilModel.from_config(dummy_carbon_data, config, pint.Quantity("12 hours")) @@ -54,14 +57,6 @@ def soil_model_fixture(dummy_carbon_data): DEBUG, "soil model: required var 'bulk_density' checked", ), - ( - DEBUG, - "soil model: required var 'soil_moisture' checked", - ), - ( - DEBUG, - "soil model: required var 'soil_temperature' checked", - ), ( DEBUG, "soil model: required var 'percent_clay' checked", @@ -93,14 +88,6 @@ def soil_model_fixture(dummy_carbon_data): ERROR, "soil model: init data missing required var 'bulk_density'", ), - ( - ERROR, - "soil model: init data missing required var 'soil_moisture'", - ), - ( - ERROR, - "soil model: init data missing required var 'soil_temperature'", - ), ( ERROR, "soil model: init data missing required var 'percent_clay'", @@ -139,14 +126,6 @@ def soil_model_fixture(dummy_carbon_data): DEBUG, "soil model: required var 'bulk_density' checked", ), - ( - DEBUG, - "soil model: required var 'soil_moisture' checked", - ), - ( - DEBUG, - "soil model: required var 'soil_temperature' checked", - ), ( DEBUG, "soil model: required var 'percent_clay' checked", @@ -181,9 +160,9 @@ def test_soil_model_initialization( [0.05, 0.02, 0.1, -0.005], dims=["cell_id"] ) # Initialise model with bad data object - model = SoilModel(carbon_data, pint.Quantity("1 week")) + model = SoilModel(carbon_data, pint.Quantity("1 week"), 2, 10) else: - model = SoilModel(dummy_carbon_data, pint.Quantity("1 week")) + model = SoilModel(dummy_carbon_data, pint.Quantity("1 week"), 2, 10) # In cases where it passes then checks that the object has the right properties assert set( @@ -219,7 +198,8 @@ def test_soil_model_initialization( "timing": { "start_date": "2020-01-01", "update_interval": "12 hours", - } + }, + "layers": {"soil_layers": 2, "canopy_layers": 10}, }, }, pint.Quantity("12 hours"), @@ -250,14 +230,6 @@ def test_soil_model_initialization( DEBUG, "soil model: required var 'bulk_density' checked", ), - ( - DEBUG, - "soil model: required var 'soil_moisture' checked", - ), - ( - DEBUG, - "soil model: required var 'soil_temperature' checked", - ), ( DEBUG, "soil model: required var 'percent_clay' checked", @@ -440,7 +412,10 @@ def test_order_independance(dummy_carbon_data, soil_model_fixture): # Use this new data to make a new soil model object config = { - "core": {"timing": {"start_date": "2020-01-01", "update_interval": "12 hours"}}, + "core": { + "timing": {"start_date": "2020-01-01", "update_interval": "12 hours"}, + "layers": {"soil_layers": 2, "canopy_layers": 10}, + }, } new_soil_model = SoilModel.from_config(new_data, config, pint.Quantity("12 hours")) diff --git a/tests/test_main.py b/tests/test_main.py index 045f62d86..b345a8073 100644 --- a/tests/test_main.py +++ b/tests/test_main.py @@ -82,7 +82,7 @@ def test_select_models(caplog, model_list, no_models, raises, expected_log_entri "config,update_interval,output,raises,expected_log_entries", [ pytest.param( - {}, + {"core": {"layers": {"soil_layers": 2, "canopy_layers": 10}}}, pint.Quantity("7 days"), "SoilModel(update_interval = 7 day)", does_not_raise(), @@ -113,14 +113,6 @@ def test_select_models(caplog, model_list, no_models, raises, expected_log_entri DEBUG, "soil model: required var 'bulk_density' checked", ), - ( - DEBUG, - "soil model: required var 'soil_moisture' checked", - ), - ( - DEBUG, - "soil model: required var 'soil_temperature' checked", - ), ( DEBUG, "soil model: required var 'percent_clay' checked", @@ -129,7 +121,7 @@ def test_select_models(caplog, model_list, no_models, raises, expected_log_entri id="valid config", ), pytest.param( - {}, + {"core": {"layers": {"soil_layers": 2, "canopy_layers": 10}}}, pint.Quantity("1 minute"), None, pytest.raises(InitialisationError), @@ -153,7 +145,7 @@ def test_select_models(caplog, model_list, no_models, raises, expected_log_entri id="update interval too short", ), pytest.param( - {}, + {"core": {"layers": {"soil_layers": 2, "canopy_layers": 10}}}, pint.Quantity("1 year"), None, pytest.raises(InitialisationError), diff --git a/virtual_rainforest/core/core_schema.json b/virtual_rainforest/core/core_schema.json index b166962f7..813df8ae5 100644 --- a/virtual_rainforest/core/core_schema.json +++ b/virtual_rainforest/core/core_schema.json @@ -145,6 +145,29 @@ "out_path_initial", "out_path_final" ] + }, + "layers": { + "description": "Layers to create vertical structure", + "type": "object", + "properties": { + "soil_layers": { + "description": "Number of soil layers to simulate", + "type": "integer", + "exclusiveMinimum": 0, + "default": 2 + }, + "canopy_layers": { + "description": "Number of canopy layers to simulate", + "type": "integer", + "exclusiveMinimum": 0, + "default": 10 + } + }, + "default": {}, + "required": [ + "soil_layers", + "canopy_layers" + ] } }, "default": {}, @@ -153,7 +176,8 @@ "data_output_options", "grid", "timing", - "modules" + "modules", + "layers" ] } }, diff --git a/virtual_rainforest/core/data.py b/virtual_rainforest/core/data.py index 398fedab0..9f1ea89e2 100644 --- a/virtual_rainforest/core/data.py +++ b/virtual_rainforest/core/data.py @@ -107,7 +107,13 @@ [[core.data.variable]] var_name="elev" -Data configurations must not contain repeated data variable names. +Data configurations must not contain repeated data variable names. NOTE: At the moment, +```core.data.variable``` tags cannot be used across multiple toml config files without +causing ```ConfigurationError: Duplicated entries in config files: core.data.variable``` +to be raised. This means that all variables need to be combined in one ```config``` +file. + + .. code-block:: python diff --git a/virtual_rainforest/core/utils.py b/virtual_rainforest/core/utils.py index f0e8add2e..e8ecb95bf 100644 --- a/virtual_rainforest/core/utils.py +++ b/virtual_rainforest/core/utils.py @@ -6,8 +6,9 @@ """ # noqa: D205, D415 from pathlib import Path +from typing import List -from virtual_rainforest.core.exceptions import ConfigurationError +from virtual_rainforest.core.exceptions import ConfigurationError, InitialisationError from virtual_rainforest.core.logger import LOGGER @@ -57,3 +58,64 @@ def check_outfile(merge_file_path: Path) -> None: raise to_raise return None + + +def set_layer_roles(canopy_layers: int, soil_layers: int) -> List[str]: + """Create a list of layer roles. + + This function creates a list of layer roles for the vertical dimension of the + Virtual Rainforest. The layer above the canopy is defined as 0 (canopy height + 2m) + and the index increases towards the bottom of the soil column. The canopy includes a + maximum number of canopy layers (defined in config) which are filled from the top + with canopy node heights from the plant module (the rest is set to NaN). Below the + canopy, we currently set one subcanopy layer (around 1.5m above ground) and one + surface layer (0.1 m above ground). Below ground, we include a maximum number of + soil layers (defined in config); the deepest layer is currently set to 1 m as the + temperature there is fairly constant and equals the mean annual temperature. + + Args: + canopy_layers: number of canopy layers soil_layers: number of soil layers + + Raises: + InitialisationError: If the number soil or canopy layers are not both positive + integers + + Returns: + List of canopy layer roles + """ + + # sanity checks for soil and canopy layers + if soil_layers < 1: + to_raise = InitialisationError( + "There has to be at least one soil layer in the Virtual Rainforest!" + ) + LOGGER.error(to_raise) + raise to_raise + + if not isinstance(soil_layers, int): + to_raise = InitialisationError("The number of soil layers must be an integer!") + LOGGER.error(to_raise) + raise to_raise + + if canopy_layers < 1: + to_raise = InitialisationError( + "There has to be at least one canopy layer in the Virtual Rainforest!" + ) + LOGGER.error(to_raise) + raise to_raise + + if canopy_layers != int(canopy_layers): + to_raise = InitialisationError( + "The number of canopy layers must be an integer!" + ) + LOGGER.error(to_raise) + raise to_raise + + layer_roles = ( + ["above"] + + ["canopy"] * canopy_layers + + ["subcanopy"] + + ["surface"] + + ["soil"] * soil_layers + ) + return layer_roles diff --git a/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py b/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py index 712b05287..bf9181f65 100644 --- a/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py +++ b/virtual_rainforest/models/abiotic_simple/abiotic_simple_model.py @@ -18,7 +18,7 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c from __future__ import annotations -from typing import Any, List +from typing import Any import numpy as np from pint import Quantity @@ -28,6 +28,7 @@ class as a child of the :class:`~virtual_rainforest.core.base_model.BaseModel` c from virtual_rainforest.core.data import Data from virtual_rainforest.core.exceptions import InitialisationError from virtual_rainforest.core.logger import LOGGER +from virtual_rainforest.core.utils import set_layer_roles from virtual_rainforest.models.abiotic_simple import simple_regression @@ -39,14 +40,15 @@ class AbioticSimpleModel(BaseModel): update_interval: Time to wait between updates of the model state. soil_layers: The number of soil layers to be modelled. canopy_layers: The initial number of canopy layers to be modelled. - initial_soil_moisture: The initial soil moisture for all layers. + initial_soil_moisture: The initial soil moisture for all layers, + [relative water content] (between 0.0 and 1.0). """ model_name = "abiotic_simple" """The model name for use in registering the model and logging.""" lower_bound_on_time_scale = "1 day" """Shortest time scale that abiotic simple model can sensibly capture.""" - upper_bound_on_time_scale = "30 day" + upper_bound_on_time_scale = "1 month" """Longest time scale that abiotic simple model can sensibly capture.""" required_init_vars = ( # TODO add temporal axis ("air_temperature_ref", ("spatial",)), @@ -69,44 +71,15 @@ def __init__( initial_soil_moisture: float, **kwargs: Any, ): - # sanity checks for soil and canopy layers - if soil_layers < 1: - to_raise = InitialisationError( - "There has to be at least one soil layer in the abiotic model!" - ) - LOGGER.error(to_raise) - raise to_raise - - if not isinstance(soil_layers, int): - to_raise = InitialisationError( - "The number of soil layers must be an integer!" - ) - LOGGER.error(to_raise) - raise to_raise - - if canopy_layers < 1: - to_raise = InitialisationError( - "There has to be at least one canopy layer in the abiotic model!" - ) - LOGGER.error(to_raise) - raise to_raise - - if canopy_layers != int(canopy_layers): - to_raise = InitialisationError( - "The number of canopy layers must be an integer!" - ) - LOGGER.error(to_raise) - raise to_raise - # sanity checks for initial soil moisture if type(initial_soil_moisture) is not float: to_raise = InitialisationError("The initial soil moisture must be a float!") LOGGER.error(to_raise) raise to_raise - if initial_soil_moisture < 0 or initial_soil_moisture > 100: + if initial_soil_moisture < 0 or initial_soil_moisture > 1: to_raise = InitialisationError( - "The initial soil moisture has to be between 0 and 100!" + "The initial soil moisture has to be between 0 and 1!" ) LOGGER.error(to_raise) raise to_raise @@ -144,8 +117,8 @@ def from_config( """ # Find number of soil and canopy layers - soil_layers = config["abiotic_simple"]["soil_layers"] - canopy_layers = config["abiotic_simple"]["canopy_layers"] + soil_layers = config["core"]["layers"]["soil_layers"] + canopy_layers = config["core"]["layers"]["canopy_layers"] initial_soil_moisture = config["abiotic_simple"]["initial_soil_moisture"] LOGGER.info( @@ -164,7 +137,7 @@ def setup(self) -> None: steps. Both variables are added directly to the self.data object. """ - # Create 1-dimenaional numpy array filled with initial soil moisture values for + # Create 1-dimensional numpy array filled with initial soil moisture values for # all soil layers and np.nan for atmosphere layers soil_moisture_values = np.repeat( a=[np.nan, self.initial_soil_moisture], @@ -188,6 +161,18 @@ def setup(self) -> None: name="soil_moisture", ) + # create soil temperature array + self.data["soil_temperature"] = DataArray( + np.full((len(self.layer_roles), len(self.data.grid.cell_id)), np.nan), + dims=["layers", "cell_id"], + coords={ + "layers": np.arange(0, len(self.layer_roles)), + "layer_roles": ("layers", self.layer_roles), + "cell_id": self.data.grid.cell_id, + }, + name="soil_temperature", + ) + # calculate vapour pressure deficit at reference height for all time steps self.data[ "vapour_pressure_deficit_ref" @@ -217,32 +202,3 @@ def update(self) -> None: def cleanup(self) -> None: """Placeholder function for abiotic model cleanup.""" - - -def set_layer_roles(canopy_layers: int, soil_layers: int) -> List[str]: - """Create a list of layer roles. - - This function creates a list of layer roles for the vertical dimension of the - Virtual Rainforest. The layer above the canopy is defined as 0 (canopy height + 2m) - and the index increases towards the bottom of the soil column. The canopy includes - a maximum number of canopy layers (defined in config) which are filled from the top - with canopy node heights from the plant module (the rest is set to NaN). Below the - canopy, we currently set one subcanopy layer (around 1.5m above ground) and one - surface layer (0.1 m above ground). Below ground, we include a maximum number of - soil layers (defined in config); the deepest layer is currently set to 1 m as the - temperature there is fairly constant and equals the mean annual temperature. - - Args: - canopy_layers: number of canopy layers - soil_layers: number of soil layers - - Returns: - List of canopy layer roles - """ - return ( - ["above"] - + ["canopy"] * canopy_layers - + ["subcanopy"] - + ["surface"] - + ["soil"] * soil_layers - ) diff --git a/virtual_rainforest/models/abiotic_simple/abiotic_simple_schema.json b/virtual_rainforest/models/abiotic_simple/abiotic_simple_schema.json index 5076d4e22..42b063c26 100644 --- a/virtual_rainforest/models/abiotic_simple/abiotic_simple_schema.json +++ b/virtual_rainforest/models/abiotic_simple/abiotic_simple_schema.json @@ -5,29 +5,16 @@ "description": "Configuration settings for the simple abiotic module", "type": "object", "properties": { - "soil_layers": { - "description": "Number of soil layers to simulate", - "type": "integer", - "exclusiveMinimum": 0, - "default": 2 - }, - "canopy_layers": { - "description": "Number of canopy layers to simulate", - "type": "integer", - "exclusiveMinimum": 0, - "default": 10 - }, "initial_soil_moisture": { "description": "Initial soil moisture for all layers", "type": "number", "exclusiveMinimum": 0, - "default": 50 + "default": 0.5, + "units": "Relative water content (between 0.0 and 1.0)" } }, "default": {}, "required": [ - "soil_layers", - "canopy_layers", "initial_soil_moisture" ] } diff --git a/virtual_rainforest/models/abiotic_simple/simple_regression.py b/virtual_rainforest/models/abiotic_simple/simple_regression.py index 5ee1df532..ade56c62b 100644 --- a/virtual_rainforest/models/abiotic_simple/simple_regression.py +++ b/virtual_rainforest/models/abiotic_simple/simple_regression.py @@ -31,7 +31,7 @@ """ MicroclimateParameters: Dict[str, float] = { - "soil_moisture_capacity": 90, + "soil_moisture_capacity": 0.9, "water_interception_factor": 0.1, "saturation_vapour_pressure_factor1": 0.61078, "saturation_vapour_pressure_factor2": 7.5, @@ -133,12 +133,14 @@ def run_simple_regression( log interpolation. Note that currently no conservation of water and energy! water_interception_factor: Factor that determines how much rainfall is intercepted by stem and canopy - soil_moisture_capacity: soil moisture capacity for water + soil_moisture_capacity: soil moisture capacity for water, relative water content + (between 0.0 and 1.0) Returns: Dict of DataArrays for air temperature [C], relative humidity [-], vapour pressure deficit [kPa], soil temperature [C], atmospheric pressure [kPa], - atmospheric :math:`\ce{CO2}` [ppm], soil moisture [-], and surface runoff [mm] + atmospheric :math:`\ce{CO2}` [ppm], soil moisture [relative water content], and + surface runoff [mm] """ # TODO correct gap between 1.5 m and 2m reference height for LAI = 0 @@ -152,7 +154,7 @@ def run_simple_regression( for var in ["air_temperature", "relative_humidity", "vapour_pressure_deficit"]: output[var] = log_interpolation( data=data, - reference_data=data[var + "_ref"].isel(time=time_index), + reference_data=data[var + "_ref"].isel(time_index=time_index), leaf_area_index_sum=leaf_area_index_sum, layer_roles=layer_roles, layer_heights=data["layer_heights"], @@ -164,7 +166,7 @@ def run_simple_regression( # Mean atmospheric pressure profile, [kPa] output["atmospheric_pressure"] = ( data["atmospheric_pressure_ref"] - .isel(time=time_index) + .isel(time_index=time_index) .where(output["air_temperature"].coords["layer_roles"] != "soil") .rename("atmospheric_pressure") .T @@ -173,7 +175,7 @@ def run_simple_regression( # Mean atmospheric C02 profile, [ppm] output["atmospheric_co2"] = ( data["atmospheric_co2_ref"] - .isel(time=0) + .isel(time_index=0) .where(output["air_temperature"].coords["layer_roles"] != "soil") .rename("atmospheric_co2") .T @@ -202,7 +204,7 @@ def run_simple_regression( ) # Precipitation at the surface is reduced as a function of leaf area index - precipitation_surface = data["precipitation"].isel(time=time_index) * ( + precipitation_surface = data["precipitation"].isel(time_index=time_index) * ( 1 - water_interception_factor * data["leaf_area_index"].sum(dim="layers") ) @@ -399,7 +401,7 @@ def calculate_soil_moisture( "soil_moisture_capacity" ], ) -> Tuple[DataArray, DataArray]: - """Calculate surface runoff and update soil mositure content. + """Calculate surface runoff and update soil moisture. Soil moisture and surface runoff are calculated with a simple bucket model: if precipitation exceeds soil moisture capacity (see MicroclimateParameters), the @@ -412,10 +414,12 @@ def calculate_soil_moisture( surface, soil) precipitation_surface: precipitation that reaches surface, [mm], current_soil_moisture: current soil moisture at upper layer, [mm], - soil_moisture_capacity: soil moisture capacity (optional) + soil_moisture_capacity: soil moisture capacity (optional), [relative water + content] Returns: - current soil moisture for one layer, [mm], surface runoff, [mm] + current soil moisture for one layer, [relative water content], surface runoff, + [mm] """ # calculate how much water can be added to soil before capacity is reached available_capacity = soil_moisture_capacity - current_soil_moisture.mean( diff --git a/virtual_rainforest/models/soil/soil_model.py b/virtual_rainforest/models/soil/soil_model.py index 53b39fd63..f22c84737 100644 --- a/virtual_rainforest/models/soil/soil_model.py +++ b/virtual_rainforest/models/soil/soil_model.py @@ -30,6 +30,7 @@ from virtual_rainforest.core.data import Data from virtual_rainforest.core.exceptions import InitialisationError from virtual_rainforest.core.logger import LOGGER +from virtual_rainforest.core.utils import set_layer_roles from virtual_rainforest.models.soil.carbon import calculate_soil_carbon_updates @@ -64,8 +65,6 @@ class SoilModel(BaseModel): ("soil_c_pool_microbe", ("spatial",)), ("pH", ("spatial",)), ("bulk_density", ("spatial",)), - ("soil_moisture", ("spatial",)), - ("soil_temperature", ("spatial",)), ("percent_clay", ("spatial",)), ) """Required initialisation variables for the soil model. @@ -78,6 +77,8 @@ def __init__( self, data: Data, update_interval: Quantity, + soil_layers: int, + canopy_layers: int, **kwargs: Any, ): super().__init__(data, update_interval, **kwargs) @@ -98,11 +99,11 @@ def __init__( """A Data instance providing access to the shared simulation data.""" self.update_interval """The time interval between model updates.""" - # Find first soil layer from the soil temperature in data object + # create a list of layer roles + layer_roles = set_layer_roles(canopy_layers, soil_layers) + # Find first soil layer from the list of layer roles self.top_soil_layer_index = next( - i - for i, v in enumerate(data["soil_moisture"].coords["layer_roles"]) - if v == "soil" + i for i, v in enumerate(layer_roles) if v == "soil" ) """The layer in the data object representing the first soil layer.""" # TODO - At the moment the soil model only cares about the very top layer. As @@ -124,11 +125,15 @@ def from_config( update_interval: Frequency with which all models are updated """ + # Find number of soil and canopy layers + soil_layers = config["core"]["layers"]["soil_layers"] + canopy_layers = config["core"]["layers"]["canopy_layers"] + LOGGER.info( "Information required to initialise the soil model successfully " "extracted." ) - return cls(data, update_interval) + return cls(data, update_interval, soil_layers, canopy_layers) def setup(self) -> None: """Placeholder function to setup up the soil model."""