diff --git a/wolf-sheep-soil-creep/README.md b/wolf-sheep-soil-creep/README.md
index 29985368d..c784a540a 100644
--- a/wolf-sheep-soil-creep/README.md
+++ b/wolf-sheep-soil-creep/README.md
@@ -3,7 +3,7 @@ title: Wolf-sheep-grass model with soil creep
permalink: tutorials-wolf-sheep-soil-creep.html
aliases:
- /tutorials-wolf-sheep-soil-creep.html
-keywords: MESA, Landlab, wolf-sheep-grass, soil creep, ABM, agent-based modeling
+keywords: MESA, Landlab, FEniCSx, wolf-sheep-grass, soil creep, ABM, agent-based modeling
summary: Example of bi-directional ABM-PDE coupling via preCICE.
---
@@ -31,6 +31,10 @@ 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).
+* Standalone reimplementation of the Landlab participant, tailored to this example.
+
+* FEniCSx. For more information, have a look at the [FEniCSx adapter](https://precice.org/adapter-fenicsx.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)
@@ -44,6 +48,20 @@ cd soil-creep-landlab
./run.sh
```
+or
+
+```bash
+cd soil-creep-simplified-landlab
+./run.sh
+```
+
+or
+
+```bash
+cd soil-creep-fenicsx
+./run.sh
+```
+
and
```bash
@@ -53,7 +71,7 @@ cd wolf-sheep-grass-mesa
## 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.
+The soil-creep participant generates Python-based visualizations, which are written to `soil-creep-landlab/output` or `soil-creep-fenicsx/output`, depending on which participant is used. The following examples were generated by Landlab with `rng=42` in the `WolfSheepScenario` for the wolf-sheep-grass participant.
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
index f39e09e5a..aa0fed213 100644
Binary files a/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-precice-config.png and b/wolf-sheep-soil-creep/images/tutorials-wolf-sheep-soil-creep-precice-config.png differ
diff --git a/wolf-sheep-soil-creep/precice-config.xml b/wolf-sheep-soil-creep/precice-config.xml
index 2d988c5ff..5edbeca2b 100644
--- a/wolf-sheep-soil-creep/precice-config.xml
+++ b/wolf-sheep-soil-creep/precice-config.xml
@@ -12,39 +12,31 @@
+
-
-
-
-
-
-
-
-
+
-
-
+
-
-
+
-
+
diff --git a/wolf-sheep-soil-creep/soil-creep-fenicsx/clean.sh b/wolf-sheep-soil-creep/soil-creep-fenicsx/clean.sh
new file mode 100755
index 000000000..a2e82c73e
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-fenicsx/clean.sh
@@ -0,0 +1,11 @@
+#!/usr/bin/env sh
+set -e -u
+
+SCRIPT_DIR="$(cd "$(dirname "$0")" && pwd)"
+cd "$SCRIPT_DIR"
+
+rm -rfv ./output/
+
+. ../../tools/cleaning-tools.sh
+
+clean_fenicsx .
diff --git a/wolf-sheep-soil-creep/soil-creep-fenicsx/precice-adapter-config.json b/wolf-sheep-soil-creep/soil-creep-fenicsx/precice-adapter-config.json
new file mode 100644
index 000000000..3948d55d6
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-fenicsx/precice-adapter-config.json
@@ -0,0 +1,19 @@
+{
+ "participant_name": "Soil-Creep",
+ "precice_config_file_path": "../precice-config.xml",
+ "interfaces": [
+ {
+ "mesh_name": "Soil-Creep-Mesh",
+ "write_data": [
+ {
+ "name": "Soil"
+ }
+ ],
+ "read_data": [
+ {
+ "name": "Grass"
+ }
+ ]
+ }
+ ]
+}
diff --git a/wolf-sheep-soil-creep/soil-creep-fenicsx/requirements.txt b/wolf-sheep-soil-creep/soil-creep-fenicsx/requirements.txt
new file mode 100644
index 000000000..7cef380a1
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-fenicsx/requirements.txt
@@ -0,0 +1,4 @@
+numpy >1, <2
+fenicsxprecice~=1.0
+matplotlib>=3,<4
+pyprecice~=3.0
diff --git a/wolf-sheep-soil-creep/soil-creep-fenicsx/run.sh b/wolf-sheep-soil-creep/soil-creep-fenicsx/run.sh
new file mode 100755
index 000000000..16c90cc8d
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-fenicsx/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 --system-site-packages .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-fenicsx/soil_creep.py b/wolf-sheep-soil-creep/soil-creep-fenicsx/soil_creep.py
new file mode 100644
index 000000000..2ecf34b23
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-fenicsx/soil_creep.py
@@ -0,0 +1,244 @@
+"""
+Linear diffusion equation of soil with Dirichlet boundary.
+ u' = div(D*grad(u)) in the square [0,19] x [0,19]
+ u = 0 on the bottom boundary at y=0
+"""
+
+import copy
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import ufl
+import numpy as np
+from petsc4py import PETSc
+from mpi4py import MPI
+from dolfinx import default_scalar_type, fem, mesh
+from dolfinx.fem.petsc import (
+ assemble_vector,
+ assemble_matrix,
+ create_vector,
+ apply_lifting,
+ set_bc,
+)
+from fenicsxprecice import Adapter, CouplingMesh
+
+
+def bottom_boundary(x):
+ return np.isclose(x[1], 0)
+
+# Volume coupling
+
+
+def coupling_boundary(x):
+ return np.ones(x.shape[1], dtype=bool)
+
+
+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
+
+# Set grid size (same as W-S-G model's grid)
+# Adjust refinement by changing number of nodes ny/nx (e.g. double refinement by setting ny = nx = 40).
+nx = ny = 20
+length_x = length_y = 19.0
+dx = length_x / (nx - 1)
+
+fenics_dt = 0.2 * dx * dx / fast_creep
+
+comm = MPI.COMM_WORLD
+
+# Create domain and function space
+domain = mesh.create_rectangle(
+ comm,
+ [np.array([0.0, 0.0]), np.array([length_x, length_y])],
+ [nx - 1, ny - 1],
+ mesh.CellType.triangle,
+)
+V = fem.functionspace(domain, ("Lagrange", 1))
+Q = fem.functionspace(domain, ("DG", 0))
+
+# u is elev function
+u_n = fem.Function(V)
+u_n.interpolate(lambda x: 0.1 * x[1])
+
+initial_elev = fem.Function(V)
+initial_elev.x.array[:] = u_n.x.array
+
+# D is diffusivity function / creep coefficient
+D = fem.Function(Q)
+
+# Functions for coupling
+soil = fem.Function(V)
+soil.x.array[:] = initial_soil_depth
+
+grass = fem.Function(V)
+
+# Auxiliary functions
+grass_diffusivity = fem.Function(V)
+soil_scaling = fem.Function(V)
+
+# Dirichlet boundary condition
+bottom_dofs = fem.locate_dofs_geometrical(V, bottom_boundary)
+bc = fem.dirichletbc(default_scalar_type(0), bottom_dofs, V)
+
+# preCICE setup
+participant = Adapter(adapter_config_filename="precice-adapter-config.json", mpi_comm=comm)
+
+coupling_mesh = CouplingMesh("Soil-Creep-Mesh", coupling_boundary, {"Grass": V}, {"Soil": soil})
+participant.initialize([coupling_mesh])
+
+# Define the variational formualation
+u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
+dt_constant = fem.Constant(domain, default_scalar_type(fenics_dt))
+
+a = u * v * ufl.dx + dt_constant * D * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
+L = u_n * v * ufl.dx
+
+bilinear_form = fem.form(a)
+linear_form = fem.form(L)
+
+# Create the matrix and vector for the linear problem
+b = create_vector(fem.extract_function_spaces(linear_form))
+uh = fem.Function(V)
+
+# Define a linear variational solver
+solver = PETSc.KSP().create(domain.comm)
+solver.setType(PETSc.KSP.Type.PREONLY)
+solver.getPC().setType(PETSc.PC.Type.LU)
+
+
+while participant.is_coupling_ongoing():
+ precice_dt = participant.get_max_time_step_size()
+ dt = np.min([fenics_dt, precice_dt])
+ dt_constant.value = default_scalar_type(dt)
+
+ participant.read_data("Soil-Creep-Mesh", "Grass", dt, grass)
+
+ grass_diffusivity.x.array[:] = np.where(
+ grass.x.array == 1,
+ fast_creep,
+ slow_creep,
+ )
+ grass_diffusivity.x.scatter_forward()
+
+ soil_scaling.x.array[:] = 1.0 - np.exp(-soil.x.array / hstar)
+ soil_scaling.x.scatter_forward()
+
+ D_expr = fem.Expression(
+ grass_diffusivity * soil_scaling,
+ Q.element.interpolation_points,
+ )
+ D.interpolate(D_expr)
+
+ # Reassemble problem because D is updated
+ A = assemble_matrix(bilinear_form, bcs=[bc])
+ A.assemble()
+ solver.setOperators(A)
+
+ # Update the right hand side reusing the initial vector
+ with b.localForm() as loc_b:
+ loc_b.set(0)
+ assemble_vector(b, linear_form)
+
+ apply_lifting(b, [bilinear_form], [[bc]])
+ b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
+
+ set_bc(b, [bc])
+
+ # Solve linear problem
+ solver.solve(b, uh.x.petsc_vec)
+ uh.x.scatter_forward()
+
+ soil.x.array[:] += uh.x.array - u_n.x.array
+ soil.x.scatter_forward()
+
+ # Update solution at previous time step (u_n)
+ u_n.x.array[:] = uh.x.array
+
+ participant.write_data("Soil-Creep-Mesh", "Soil", soil)
+
+ participant.advance(dt)
+
+
+participant.finalize()
+
+
+dofs_coupling_coordinates = V.tabulate_dof_coordinates()[:, :2].copy()
+
+
+def values_to_grid(values):
+ grid = np.full((ny, nx), np.nan)
+
+ for coord, value in zip(dofs_coupling_coordinates, values):
+ i = int(round(coord[0] / length_x * (nx - 1)))
+ j = int(round(coord[1] / length_y * (ny - 1)))
+ grid[j, i] = value
+
+ return grid
+
+
+# Mask for closed boundaries (mimic Landlab)
+closed = np.zeros((ny, nx), dtype=bool)
+
+closed[1:-1, 0] = True
+closed[1:-1, -1] = True
+closed[-1, :] = True
+
+# Calculate and plot the erosion/deposition patterns
+ero_dep = u_n.x.array - initial_elev.x.array
+ero_dep_grid = values_to_grid(ero_dep)
+
+maxchange = np.max(np.abs(ero_dep_grid))
+
+cmap = mpl.colormaps["coolwarm_r"].copy()
+cmap.set_bad("black")
+
+plt.figure()
+plt.imshow(
+ np.ma.array(ero_dep_grid, mask=closed),
+ origin="lower",
+ vmin=-maxchange,
+ vmax=maxchange,
+ cmap=cmap
+)
+plt.colorbar(label="Depth of soil accumulation (+) or loss (-), m")
+plt.savefig("output/erosion_deposition_patterns.png")
+plt.close()
+
+# Soil thickness
+soil_grid = values_to_grid(soil.x.array)
+
+cmap = mpl.colormaps["pink"].copy()
+cmap.set_bad("black")
+
+plt.figure()
+plt.imshow(
+ np.ma.array(soil_grid, mask=closed),
+ origin="lower",
+ cmap=cmap,
+)
+plt.colorbar(label="Soil thickness, m")
+plt.savefig("output/soil_thickness.png")
+plt.close()
+
+# Ground cover
+gm_grid = values_to_grid(grass.x.array)
+
+cmap = mpl.colormaps["YlGn"].copy()
+cmap.set_bad("black")
+
+plt.figure()
+plt.imshow(
+ np.ma.array(gm_grid, mask=closed),
+ origin="lower",
+ cmap=cmap,
+ vmin=1,
+ vmax=2,
+)
+plt.colorbar(label="Ground cover (1 = bare, 2 = grass)")
+plt.savefig("output/grass_map.png")
+plt.close()
diff --git a/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py b/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py
index 2f17e3830..2debf84ee 100644
--- a/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py
+++ b/wolf-sheep-soil-creep/soil-creep-landlab/soil_creep.py
@@ -23,6 +23,8 @@
rmg = RasterModelGrid((width, height))
+solver_dt = 0.2 * rmg.dx * rmg.dx / fast_creep
+
# 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
@@ -54,24 +56,22 @@
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)
+positions = [[x, y] for y in range(height) for x in range(width)]
+vertex_ids = participant.set_mesh_vertices("Soil-Creep-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)
+ gm = participant.read_data("Soil-Creep-Mesh", "Grass", vertex_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
+ creep_coef[gm == 1] = fast_creep
+ creep_coef[gm == 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)
@@ -86,7 +86,7 @@
# Update the soil cover
soil += elev - prior_elev
- participant.write_data("Soil-Depth-Mesh", "Soil", vertex_soil_ids, soil)
+ participant.write_data("Soil-Creep-Mesh", "Soil", vertex_ids, soil)
participant.advance(dt)
diff --git a/wolf-sheep-soil-creep/soil-creep-simplified-landlab/clean.sh b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/clean.sh
new file mode 100755
index 000000000..ebce578f2
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-simplified-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-simplified-landlab/requirements.txt b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/requirements.txt
new file mode 100644
index 000000000..29f7bd20c
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/requirements.txt
@@ -0,0 +1,3 @@
+numpy >1, <2
+matplotlib>=3,<4
+pyprecice~=3.0
diff --git a/wolf-sheep-soil-creep/soil-creep-simplified-landlab/run.sh b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/run.sh
new file mode 100755
index 000000000..6cfe51b1c
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/run.sh
@@ -0,0 +1,21 @@
+#!/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-simplified-landlab/soil_creep.py b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/soil_creep.py
new file mode 100644
index 000000000..d0613c524
--- /dev/null
+++ b/wolf-sheep-soil-creep/soil-creep-simplified-landlab/soil_creep.py
@@ -0,0 +1,252 @@
+import copy
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+import precice
+
+
+class LandlabLinearDiffuserClone:
+ CORE = 0 # regular node
+ FIXED_VALUE = 1 # Dirichlet fixed elevation
+ CLOSED = 4 # Neumann zero-flux
+
+ ALPHA = 0.15 # time-step stability factor
+
+ def __init__(self, z, D, dx):
+ self.z = z # grid
+ self.D = D # linear diffusivity
+ self.dx = dx # node spacing
+ self.ny, self.nx = z.shape # grid height and width
+ self.id = np.arange(self.nx * self.ny).reshape(self.ny, self.nx)
+
+ self.status = np.full((self.ny, self.nx), self.FIXED_VALUE, dtype=np.uint8)
+ self.status[1:-1, 1:-1] = self.CORE
+
+ # Leave the south boundary open and close the others
+ # Corresponds to Landlab's rmg.set_closed_boundaries_at_grid_edges(True, True, True, False)
+ # with the boundary order: right, top, left, bottom
+ self.status[:, -1] = self.CLOSED
+ self.status[-1, :] = self.CLOSED
+ self.status[:, 0] = self.CLOSED
+
+ self._build_links()
+
+ def _build_links(self): # Initalize list of active links: between two core nodes or between core and fixed
+ ny, nx = self.ny, self.nx
+ links = []
+
+ for j in range(ny):
+ for i in range(nx - 1):
+ links.append((self.id[j, i], self.id[j, i + 1]))
+
+ for j in range(ny - 1):
+ for i in range(nx):
+ links.append((self.id[j, i], self.id[j + 1, i]))
+
+ status = self.status.ravel()
+ active_links = []
+
+ for n1, n2 in links:
+ s1 = status[n1]
+ s2 = status[n2]
+
+ if s1 == self.CLOSED or s2 == self.CLOSED:
+ continue
+
+ if s1 != self.CORE and s2 != self.CORE:
+ continue
+
+ active_links.append((n1, n2))
+
+ self.links = np.array(active_links, dtype=int)
+
+ def step(self, dt):
+ D = self.D.ravel()
+
+ active_D = np.array([
+ np.maximum(D[n1], D[n2])
+ for n1, n2 in self.links
+ ])
+
+ # Divide the timestep dt passed to the diffuser into internal stability-limited substeps.
+ # Perform as many full substeps as possible (repeats), followed by one smaller substep for the remainder.
+ internal_dt = self.ALPHA * self.dx * self.dx / np.nanmax(active_D)
+
+ repeats = int(dt // internal_dt)
+ remainder = dt - repeats * internal_dt
+
+ for _ in range(repeats):
+ self._substep(internal_dt)
+
+ if remainder > 0:
+ self._substep(remainder)
+
+ def _substep(self, dt):
+ z = self.z.ravel()
+ D = self.D.ravel()
+ status = self.status.ravel()
+
+ dzdt = np.zeros_like(z)
+
+ for n1, n2 in self.links: # Iterate over all active links (tuples)
+ Dk = np.maximum(D[n1], D[n2]) # max link diffusivity
+ q = -Dk * (z[n2] - z[n1]) / self.dx
+
+ if status[n1] == self.CORE:
+ dzdt[n1] -= q / self.dx
+ if status[n2] == self.CORE:
+ dzdt[n2] += q / self.dx
+
+ z[status == self.CORE] += dt * dzdt[status == self.CORE]
+ self.z = z.reshape(self.ny, self.nx)
+
+
+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
+
+# Set grid size (same as W-S-G model's grid)
+# Adjust refinement by changing number of nodes ny/nx (e.g. double refinement by setting ny = nx = 40).
+ny = nx = 20
+length_x = 19
+length_y = 19
+dx = length_x / (nx - 1)
+
+solver_dt = 0.2 * dx * dx / fast_creep
+
+# Grass map field for plotting
+gm = np.zeros((ny, nx))
+
+# Create elevation field and have it slope down to the south at 10% gradient
+y = np.arange(ny)[:, None] * np.ones((ny, nx)) * dx
+elev = 0.1 * y
+
+# Remember the starting elevation so we can calculate cumulative erosion/deposition
+initial_elev = np.zeros((ny, nx))
+initial_elev[:] = elev
+
+# Also remember the elevation of the prior time step, so we can difference
+prior_elev = np.zeros((ny, nx))
+
+# Create a field for the creep coefficient
+creep_coef = np.zeros((ny, nx))
+
+# Create a soil-thickness field
+soil = np.full((ny, nx), initial_soil_depth)
+
+# Instantiate a LinearDiffuser (soil creep) Landlab component
+diffuser = LandlabLinearDiffuserClone(
+ z=elev,
+ D=creep_coef,
+ dx=dx,
+)
+
+# 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)
+
+x = np.linspace(0, length_x, nx)
+y = np.linspace(0, length_y, ny)
+
+positions = [[x0, y0] for y0 in y for x0 in x]
+vertex_ids = participant.set_mesh_vertices("Soil-Creep-Mesh", positions)
+
+participant.initialize()
+
+while participant.is_coupling_ongoing():
+ precice_dt = participant.get_max_time_step_size()
+ dt = np.minimum(solver_dt, precice_dt)
+
+ gm_flat = participant.read_data("Soil-Creep-Mesh", "Grass", vertex_ids, dt)
+ gm[:, :] = np.asarray(gm_flat).reshape(ny, nx)
+
+ # 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 == 1] = fast_creep
+ creep_coef[gm == 2] = slow_creep
+
+ # Adjust the creep coefficient to account for soil depth
+ creep_coef[:, :] *= 1.0 - np.exp(-soil / hstar)
+
+ # Remember the current elevation
+ prior_elev[:, :] = elev
+
+ # Run the soil-creep model
+ diffuser.step(dt)
+
+ # Update elevation
+ elev = diffuser.z
+
+ # Update the soil cover
+ soil += elev - prior_elev
+
+ participant.write_data("Soil-Creep-Mesh", "Soil", vertex_ids, soil.ravel())
+
+ participant.advance(dt)
+
+participant.finalize()
+
+
+# Mask for closed boundaries (mimic Landlab)
+closed = np.zeros((ny, nx), dtype=bool)
+
+closed[1:-1, 0] = True
+closed[1:-1, -1] = True
+closed[-1, :] = True
+
+# Calculate and plot the erosion/deposition patterns
+plt.figure()
+ero_dep = elev - initial_elev
+maxchange = np.max(np.abs(ero_dep))
+
+cmap = mpl.colormaps["coolwarm_r"].copy()
+cmap.set_bad("black")
+
+plt.figure()
+plt.imshow(
+ np.ma.array(ero_dep, mask=closed),
+ origin="lower",
+ vmin=-maxchange,
+ vmax=maxchange,
+ cmap=cmap
+)
+plt.colorbar(label="Depth of soil accumulation (+) or loss (-), m")
+plt.savefig("output/erosion_deposition_patterns.png")
+plt.close()
+
+# Soil thickness
+cmap = mpl.colormaps["pink"].copy()
+cmap.set_bad("black")
+
+plt.figure()
+plt.imshow(
+ np.ma.array(soil, mask=closed),
+ origin="lower",
+ cmap=cmap
+)
+plt.colorbar(label="Soil thickness, m")
+plt.savefig("output/soil_thickness.png")
+plt.close()
+
+# Ground cover
+cmap = mpl.colormaps["YlGn"].copy()
+cmap.set_bad("black")
+
+plt.figure()
+plt.imshow(
+ np.ma.array(gm, mask=closed),
+ origin="lower",
+ cmap=cmap
+)
+plt.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/wolf_sheep_grass.py b/wolf-sheep-soil-creep/wolf-sheep-grass-mesa/wolf_sheep_grass.py
index 54b267e4e..21e4f27fd 100644
--- 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
@@ -65,23 +65,22 @@ def limit_grass_by_soil(wsg_model, soil, min_soil_depth):
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)
+positions = [[x, y] for y in range(height) for x in range(width)]
-soil = np.zeros([width * height])
+vertex_ids = participant.set_mesh_vertices("Wolf-Sheep-Grass-Mesh", positions)
if participant.requires_initial_data():
gm = generate_grass_map(ws)
- participant.write_data("Wolf-Sheep-Grass-Mesh", "Grass", vertex_gm_ids, gm.flatten())
+ participant.write_data("Wolf-Sheep-Grass-Mesh", "Grass", vertex_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)
+ soil = participant.read_data("Wolf-Sheep-Grass-Mesh", "Soil", vertex_ids, dt)
# Update the grass cover
limit_grass_by_soil(ws, soil, min_depth_for_grass)
@@ -90,8 +89,9 @@ def limit_grass_by_soil(wsg_model, soil, min_soil_depth):
ws.step()
gm = generate_grass_map(ws)
- participant.write_data("Wolf-Sheep-Grass-Mesh", "Grass", vertex_gm_ids, gm.flatten())
+ participant.write_data("Wolf-Sheep-Grass-Mesh", "Grass", vertex_ids, gm.flatten())
participant.advance(dt)
+
participant.finalize()