diff --git a/changelog-entries/793.md b/changelog-entries/793.md new file mode 100644 index 000000000..7dc2bfa6d --- /dev/null +++ b/changelog-entries/793.md @@ -0,0 +1 @@ +- Added the `wolf-sheep-soil-creep` tutorial [#793](https://github.com/precice/tutorials/pull/793) \ No newline at end of file diff --git a/wolf-sheep-soil-creep/README.md b/wolf-sheep-soil-creep/README.md new file mode 100644 index 000000000..b084c7285 --- /dev/null +++ b/wolf-sheep-soil-creep/README.md @@ -0,0 +1,58 @@ +--- +title: Wolf-sheep-grass model with soil creep +permalink: tutorials-wolf-sheep-soil-creep.html +keywords: MESA, Landlab, wolf-sheep-grass, soil creep, ABM, agent-based modeling +summary: Example of bi-directional ABM-PDE coupling via preCICE. +--- + +{% note %} +Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/develop/wolf-sheep-soil-creep), as continuously rendered here, or see the [latest released version](https://github.com/precice/tutorials/tree/master/wolf-sheep-soil-creep) (if there is already one). Read how in the [tutorials introduction](https://precice.org/tutorials.html). +{% endnote %} + +## Setup + +This tutorial is based on Landlab's [Wolf-Sheep-Grass Model with Soil Creep](https://landlab.csdms.io/tutorials/agent_based_modeling/wolf_sheep/wolf_sheep_with_soil_creep.html) example, which couples the [Wolf-Sheep model from Mesa](https://mesa.readthedocs.io/latest/examples/advanced/wolf_sheep.html) with a Landlab soil creep model in a single notebook-based workflow. The [source notebook](https://github.com/landlab/landlab/blob/master/docs/source/tutorials/agent_based_modeling/wolf_sheep/wolf_sheep_with_soil_creep.ipynb) is licensed under [MIT](https://landlab.csdms.io/about/license.html). + +The purpose of this tutorial is to demonstrate a simple two-way coupling between an agent-based model (ABM) and a grid-based continuum model via preCICE. The example extends the canonical wolf-sheep-grass ABM with a simple soil-creep feedback on a shared sloping raster grid. When sheep eat grass, the soil beneath the damaged grass becomes more mobile, increasing the local soil transport efficiency and affecting downslope soil transport. The feedback also acts in the other direction: soil thickness affects grass growth by preventing grass from growing where the soil becomes too thin. + +In this implementation, the Mesa wolf-sheep-grass model and the Landlab soil creep model are two separate preCICE participants, exchanging grass cover and soil depth during the simulation. The Mesa participant extracts the current grass map from the wolf-sheep-grass model and sends it to the Landlab participant. The Landlab participant uses this grass map to adapt the soil creep coefficient, evolves the soil thickness, and sends the updated soil depth back to the Mesa participant, where it is used to limit grass growth. This example is intentionally simple and should be understood as a demonstration case for coupling ABM and PDE-style models via preCICE, rather than as a physically calibrated landscape evolution model. + +## Configuration + +preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)): + +![preCICE configuration visualization](images/tutorials-wolf-sheep-soil-creep-precice-config.png) + +## Available solvers + +Soil-Creep PDE participant: + +* Landlab. Numerical modeling of Earth surface dynamics. For more information, have a look at the [Landlab documentation](https://landlab.csdms.io/index.html). + +Wolf-Sheep-Grass ABM participant: + +* MESA. Agent-based modeling (or ABM) framework. For more information have a look at the [MESA documentation](https://mesa.readthedocs.io/latest/index.html) + +## Running the simulation + +Open two separate terminals and start the soil creep and wolf sheep participants by calling the respective run scripts `run.sh` located in each of the participants' directories. For example: + +```bash +cd soil-creep-landlab +./run.sh +``` + +and + +```bash +cd wolf-sheep-grass-mesa +./run.sh +``` + +## Post-processing + +The soil-creep participant generates Python-based visualizations, which are written to `soil-creep-landlab/output`. The following examples were generated by setting `rng=42` in the `WolfSheepScenario` for the wolf-sheep-grass participant. + +![erosion/deposition patterns](images/tutorials-wolf-sheep-soil-creep-erosion_deposition_patterns.png) +![soil thickness](images/tutorials-wolf-sheep-soil-creep-soil_thickness.png) +![grass map](images/tutorials-wolf-sheep-soil-creep-grass_map.png) diff --git a/wolf-sheep-soil-creep/clean-tutorial.sh b/wolf-sheep-soil-creep/clean-tutorial.sh new file mode 120000 index 000000000..4713f5092 --- /dev/null +++ b/wolf-sheep-soil-creep/clean-tutorial.sh @@ -0,0 +1 @@ +../tools/clean-tutorial-base.sh \ No newline at end of file diff --git a/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-erosion_deposition_patterns.png b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-erosion_deposition_patterns.png new file mode 100644 index 000000000..718cabee7 Binary files /dev/null and b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-erosion_deposition_patterns.png differ diff --git a/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-grass_map.png b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-grass_map.png new file mode 100644 index 000000000..4fae85a0c Binary files /dev/null and b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-grass_map.png differ diff --git a/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-precice-config.png b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-precice-config.png new file mode 100644 index 000000000..f39e09e5a Binary files /dev/null and b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-precice-config.png differ diff --git a/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-soil_thickness.png b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-soil_thickness.png new file mode 100644 index 000000000..0bb5d8801 Binary files /dev/null and b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-soil_thickness.png differ diff --git a/wolf-sheep-soil-creep/precice-config.xml b/wolf-sheep-soil-creep/precice-config.xml new file mode 100644 index 000000000..2d988c5ff --- /dev/null +++ b/wolf-sheep-soil-creep/precice-config.xml @@ -0,0 +1,69 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/wolf-sheep-soil-creep/soil-creep-landlab/.gitignore b/wolf-sheep-soil-creep/soil-creep-landlab/.gitignore new file mode 100644 index 000000000..ea1472ec1 --- /dev/null +++ b/wolf-sheep-soil-creep/soil-creep-landlab/.gitignore @@ -0,0 +1 @@ +output/ diff --git a/wolf-sheep-soil-creep/soil-creep-landlab/clean.sh b/wolf-sheep-soil-creep/soil-creep-landlab/clean.sh new file mode 100755 index 000000000..ebce578f2 --- /dev/null +++ b/wolf-sheep-soil-creep/soil-creep-landlab/clean.sh @@ -0,0 +1,9 @@ +#!/usr/bin/env sh +set -e -u + +rm -rfv ./output/ + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/wolf-sheep-soil-creep/soil-creep-landlab/requirements.txt b/wolf-sheep-soil-creep/soil-creep-landlab/requirements.txt new file mode 100644 index 000000000..16efd4ddd --- /dev/null +++ b/wolf-sheep-soil-creep/soil-creep-landlab/requirements.txt @@ -0,0 +1,4 @@ +numpy >1, <2 +matplotlib +landlab +pyprecice~=3.0 diff --git a/wolf-sheep-soil-creep/soil-creep-landlab/run.sh b/wolf-sheep-soil-creep/soil-creep-landlab/run.sh new file mode 100755 index 000000000..14753ac3d --- /dev/null +++ b/wolf-sheep-soil-creep/soil-creep-landlab/run.sh @@ -0,0 +1,23 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh + +exec > >(tee --append "$LOGFILE") 2>&1 + +if [ ! -v PRECICE_TUTORIALS_NO_VENV ]; then + + if [ ! -d ".venv" ]; then + python3 -m venv .venv + source .venv/bin/activate + pip install -r requirements.txt && pip freeze > pip-installed-packages.log + else + source .venv/bin/activate + fi + +fi + +mkdir -p output +python3 soil_creep.py + +close_log diff --git a/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py b/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py new file mode 100644 index 000000000..2f17e3830 --- /dev/null +++ b/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py @@ -0,0 +1,122 @@ +import copy +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +from landlab import RasterModelGrid, imshow_grid +from landlab.components import LinearDiffuser +import precice + +initial_soil_depth = 0.3 + +# Set the soil-thickness scale for limiting creep where little soil is available +hstar = 0.2 + +# Set parameters for two soil creep coefficients: slow (full grass cover) and fast (partial or "eaten" grass cover) +fast_creep = 0.1 +slow_creep = 0.001 + +ground_cover_cmap = copy.copy(mpl.colormaps["YlGn"]) + +# Create a grid the same size as the W-S-G model's grid +width = 20 +height = 20 + +rmg = RasterModelGrid((width, height)) + +# Create elevation field and have it slope down to the south at 10% gradient +elev = rmg.add_zeros("topographic__elevation", at="node") +elev[:] = 0.1 * rmg.y_of_node + +# Have one open boundary on the south side +rmg.set_closed_boundaries_at_grid_edges(True, True, True, False) + +# Remember the starting elevation so we can calculate cumulative erosion/deposition +initial_elev = np.zeros(rmg.number_of_nodes) +initial_elev[:] = elev + +# Also remember the elevation of the prior time step, so we can difference +prior_elev = np.zeros(rmg.number_of_nodes) + +# Create a field for the creep coefficient +creep_coef = rmg.add_zeros("creep_coefficient", at="node") + +# Create a soil-thickness field +soil = rmg.add_zeros("soil__depth", at="node") +soil[:] = initial_soil_depth + +# Instantiate a LinearDiffuser (soil creep) Landlab component +diffuser = LinearDiffuser(rmg, linear_diffusivity=creep_coef) + +# preCICE setup +participant_name = "Soil-Creep" +config_file_name = "../precice-config.xml" +solver_process_index = 0 +solver_process_size = 1 +participant = precice.Participant(participant_name, config_file_name, solver_process_index, solver_process_size) + +positions = [[x, y] for x in range(width) for y in range(height)] +vertex_gm_ids = participant.set_mesh_vertices("Soil-Creep-Mesh", positions) +vertex_soil_ids = participant.set_mesh_vertices("Soil-Depth-Mesh", positions) + +participant.initialize() + +while participant.is_coupling_ongoing(): + solver_dt = 0.2 * rmg.dx * rmg.dx / fast_creep + precice_dt = participant.get_max_time_step_size() + dt = np.minimum(solver_dt, precice_dt) + + gm = participant.read_data("Soil-Creep-Mesh", "Grass", vertex_gm_ids, dt) + + # Assign the higher creep coefficient to cells where the grass has + # been eaten and not yet recovered; the slower value is assigned to + # "fully grown" grass patches. + creep_coef[gm.flatten() == 1] = fast_creep + creep_coef[gm.flatten() == 2] = slow_creep + + # Limit the creep coefficient according to the soil-thickness field, so absent soil cannot move. + creep_coef *= 1.0 - np.exp(-soil / hstar) + + # Remember the current elevation before LinearDiffuser updates the grid's + # topographic__elevation field in place. + prior_elev[:] = elev + + # Run the soil-creep model + diffuser.run_one_step(dt) + + # Update the soil cover + soil += elev - prior_elev + + participant.write_data("Soil-Depth-Mesh", "Soil", vertex_soil_ids, soil) + + participant.advance(dt) + +participant.finalize() + +# Calculate and plot the erosion/deposition patterns +plt.figure() +ero_dep = elev - initial_elev +maxchange = np.amax(np.abs(ero_dep)) +imshow_grid( + rmg, + ero_dep, + vmin=-maxchange, + vmax=maxchange, + cmap="coolwarm_r", + colorbar_label="Depth of soil accumulation (+) or loss (-), m", +) +plt.savefig("output/erosion_deposition_patterns.png") +plt.close() + +# Soil thickness +plt.figure() +imshow_grid(rmg, soil, colorbar_label="Soil thickness, m") +plt.savefig("output/soil_thickness.png") +plt.close() + +# Ground cover +plt.figure() +imshow_grid( + rmg, gm, cmap=ground_cover_cmap, colorbar_label="Ground cover (1 = bare, 2 = grass)" +) +plt.savefig("output/grass_map.png") +plt.close() diff --git a/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/clean.sh b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/clean.sh new file mode 100755 index 000000000..494c80414 --- /dev/null +++ b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/clean.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/requirements.txt b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/requirements.txt new file mode 100644 index 000000000..f9ccca730 --- /dev/null +++ b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/requirements.txt @@ -0,0 +1,5 @@ +numpy >1, <2 +matplotlib +mesa>=3 +pyprecice~=3.0 +networkx diff --git a/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/run.sh b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/run.sh new file mode 100755 index 000000000..43c20ae0d --- /dev/null +++ b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/run.sh @@ -0,0 +1,22 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh + +exec > >(tee --append "$LOGFILE") 2>&1 + +if [ ! -v PRECICE_TUTORIALS_NO_VENV ]; then + + if [ ! -d ".venv" ]; then + python3 -m venv .venv + source .venv/bin/activate + pip install -r requirements.txt && pip freeze > pip-installed-packages.log + else + source .venv/bin/activate + fi + +fi + +python3 wolf_sheep_grass.py + +close_log diff --git a/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/wolf_sheep_grass.py b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/wolf_sheep_grass.py new file mode 100644 index 000000000..54b267e4e --- /dev/null +++ b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/wolf_sheep_grass.py @@ -0,0 +1,97 @@ +from mesa.examples.advanced.wolf_sheep.agents import GrassPatch +from mesa.examples.advanced.wolf_sheep.model import WolfSheep +from mesa.examples.advanced.wolf_sheep.model import WolfSheepScenario +import copy +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +import precice + +min_depth_for_grass = 0.2 + +# Use the same time step size as the Soil-Creep participant +solver_dt = 2.0 + +# Initialize Wolf-Sheep-Grass model +wss = WolfSheepScenario( + initial_sheep=20, + initial_wolves=10, + grass=True, + grass_regrowth_time=15, # give grass a fighting chance... + rng=42, # fixed seed for reproducible results +) +ws = WolfSheep(wss) + +ground_cover_cmap = copy.copy(mpl.colormaps["YlGn"]) + + +def generate_grass_map(model): + grass_map = np.zeros((model.grid.width, model.grid.height)) + for cell in model.grid: + (x, y) = cell.coordinate + cell_content = cell.agents + for agent in cell_content: + if type(agent) is GrassPatch: + if agent.fully_grown: + grass_map[x][y] = 2 + else: + grass_map[x][y] = 1 + return grass_map + + +def plot_grass_map(grass_map): + plt.imshow(grass_map, interpolation="nearest", cmap=ground_cover_cmap) + plt.colorbar() + + +def limit_grass_by_soil(wsg_model, soil, min_soil_depth): + soilmatrix = soil.reshape((wsg_model.width, wsg_model.height)) + for cell in wsg_model.grid: + (x, y) = cell.coordinate + cell_content = cell.agents + if soilmatrix[x][y] < min_soil_depth: + for agent in cell_content: + if type(agent) is GrassPatch: + agent.fully_grown = False + + +width = ws.grid.width +height = ws.grid.height + +# preCICE setup +participant_name = "Wolf-Sheep-Grass" +config_file_name = "../precice-config.xml" +solver_process_index = 0 +solver_process_size = 1 +participant = precice.Participant(participant_name, config_file_name, solver_process_index, solver_process_size) + +positions = [[x, y] for x in range(width) for y in range(height)] +vertex_gm_ids = participant.set_mesh_vertices("Wolf-Sheep-Grass-Mesh", positions) +vertex_soil_ids = participant.set_mesh_vertices("Soil-Grass-Mesh", positions) + +soil = np.zeros([width * height]) + +if participant.requires_initial_data(): + gm = generate_grass_map(ws) + participant.write_data("Wolf-Sheep-Grass-Mesh", "Grass", vertex_gm_ids, gm.flatten()) + +participant.initialize() + +while participant.is_coupling_ongoing(): + precice_dt = participant.get_max_time_step_size() + dt = np.minimum(solver_dt, precice_dt) + + soil = participant.read_data("Soil-Grass-Mesh", "Soil", vertex_soil_ids, dt) + + # Update the grass cover + limit_grass_by_soil(ws, soil, min_depth_for_grass) + + # Run the W-S-G model + ws.step() + + gm = generate_grass_map(ws) + participant.write_data("Wolf-Sheep-Grass-Mesh", "Grass", vertex_gm_ids, gm.flatten()) + + participant.advance(dt) + +participant.finalize()