Add fenicsx and simplified Landlab solvers

wolf-sheep-soil-creep/precice-config.xml

@@ -1,69 +1,61 @@
 <?xml version="1.0" encoding="UTF-8" ?>
 <precice-configuration>
+
   <log>
     <sink
       filter="%Severity% > debug"
       format="---[precice] %ColorizedSeverity% %Message%"
       enabled="true" />
   </log>
 
+
   <data:scalar name="Grass" />
   <data:scalar name="Soil" />
 
   <mesh name="Wolf-Sheep-Grass-Mesh" dimensions="2">
     <use-data name="Grass" />
+    <use-data name="Soil" />
   </mesh>
 
   <mesh name="Soil-Creep-Mesh" dimensions="2">
     <use-data name="Grass" />
-  </mesh>
-
-  <mesh name="Soil-Grass-Mesh" dimensions="2">
-    <use-data name="Soil" />
-  </mesh>
-
-  <mesh name="Soil-Depth-Mesh" dimensions="2">
     <use-data name="Soil" />
   </mesh>
 
   <participant name="Wolf-Sheep-Grass">
-    <receive-mesh name="Soil-Depth-Mesh" from="Soil-Creep" />
+    <receive-mesh name="Soil-Creep-Mesh" from="Soil-Creep" />
     <provide-mesh name="Wolf-Sheep-Grass-Mesh" />
-    <provide-mesh name="Soil-Grass-Mesh" />
     <write-data name="Grass" mesh="Wolf-Sheep-Grass-Mesh" />
-    <read-data name="Soil" mesh="Soil-Grass-Mesh" />
+    <read-data name="Soil" mesh="Wolf-Sheep-Grass-Mesh" />
     <mapping:nearest-neighbor
       direction="read"
-      from="Soil-Depth-Mesh"
-      to="Soil-Grass-Mesh"
+      from="Soil-Creep-Mesh"
+      to="Wolf-Sheep-Grass-Mesh"
       constraint="consistent" />
   </participant>
 
   <participant name="Soil-Creep">
     <receive-mesh name="Wolf-Sheep-Grass-Mesh" from="Wolf-Sheep-Grass" />
     <provide-mesh name="Soil-Creep-Mesh" />
-    <provide-mesh name="Soil-Depth-Mesh" />
     <read-data name="Grass" mesh="Soil-Creep-Mesh" />
-    <write-data name="Soil" mesh="Soil-Depth-Mesh" />
+    <write-data name="Soil" mesh="Soil-Creep-Mesh" />
     <mapping:nearest-neighbor
       direction="read"
       from="Wolf-Sheep-Grass-Mesh"
       to="Soil-Creep-Mesh"
       constraint="consistent" />
   </participant>
 
+
   <m2n:sockets acceptor="Wolf-Sheep-Grass" connector="Soil-Creep" exchange-directory=".." />
 
+
   <coupling-scheme:serial-explicit>
     <participants first="Soil-Creep" second="Wolf-Sheep-Grass" />
     <time-window-size value="2" />
     <max-time value="100" />
-    <exchange
-      data="Grass"
-      mesh="Wolf-Sheep-Grass-Mesh"
-      from="Wolf-Sheep-Grass"
-      to="Soil-Creep"
-      initialize="true" />
-    <exchange data="Soil" mesh="Soil-Depth-Mesh" from="Soil-Creep" to="Wolf-Sheep-Grass" />
+    <exchange data="Grass" mesh="Wolf-Sheep-Grass-Mesh" from="Wolf-Sheep-Grass" to="Soil-Creep" initialize="true"/>
+    <exchange data="Soil" mesh="Soil-Creep-Mesh" from="Soil-Creep" to="Wolf-Sheep-Grass" />
   </coupling-scheme:serial-explicit>
+
 </precice-configuration>

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 .

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"
+        }
+      ]
+    }
+  ]
+}

wolf-sheep-soil-creep/soil-creep-fenicsx/requirements.txt

@@ -0,0 +1,3 @@
+numpy >1, <2
+fenicsxprecice~=1.0
+matplotlib>=3,<4

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

wolf-sheep-soil-creep/soil-creep-fenicsx/soil_creep.py

@@ -0,0 +1,241 @@
+"""
+Linear diffusion equation of soil with Dirichlet boundary.
+  z' = div(D*grad(z))    in the square [0,19] x [0,19]
+  z = 0                  on the bottom boundary
+"""
+
+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()