{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Escarpment [Py]\n",
                "\n",
                "A simple example simulating the evolution of an escarpment (on a 2D raster grid) under the action of bedrock channel and hillslope erosion.\n",
                "\n",
                "This example is pretty similar to {doc}`mountain_py`, although solving equation {eq}`eq_mountain` without the uplift term {math}`U` on a semi-infinite (periodic) domain. It also starts from a different initial topograhic surface and routes flow using a multiple direction approach."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "import fastscapelib as fs\n",
                "import numpy as np\n",
                "import matplotlib.pyplot as plt"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": { "tags": ["remove-cell"] },
            "outputs": [],
            "source": [
                "# Theme that looks reasonably fine on both dark/light modes\n",
                "import matplotlib\n",
                "matplotlib.style.use('Solarize_Light2')\n",
                "matplotlib.rcParams['axes.grid'] = False"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Setup the Grid, Flow Graph and Eroders\n",
                "\n",
                "Create a {py:class}`~fastscapelib.RasterGrid` of 101x201 nodes with a total length of 10 km in y (rows) and 20 km in x (columns).\n",
                "\n",
                "Set fixed-value boundary at the left border, free (core) boundary at the right border and reflective boundaries for the top-down borders."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "bs = [\n",
                "    fs.NodeStatus.FIXED_VALUE,\n",
                "    fs.NodeStatus.CORE,\n",
                "    fs.NodeStatus.LOOPED,\n",
                "    fs.NodeStatus.LOOPED,\n",
                "]\n",
                "\n",
                "grid = fs.RasterGrid.from_length([101, 201], [1e4, 2e4], bs)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Create a {py:class}`~fastscapelib.FlowGraph` object with multiple direction flow routing and the resolution of closed depressions on the topographic surface. See {ref}`guide-flow-routing-strategies` for more examples on possible flow routing strategies.\n",
                "\n",
                "By default, base level nodes are set from fixed value boundary conditions (all nodes on the left border in this example)."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "flow_graph = fs.FlowGraph(\n",
                "    grid,\n",
                "    [fs.SingleFlowRouter(), fs.MSTSinkResolver(), fs.MultiFlowRouter(1.1)],\n",
                ")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source":
                ["Setup eroder classes (bedrock channel + hillslope) with a given set of parameter values."]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "spl_eroder = fs.SPLEroder(\n",
                "    flow_graph,\n",
                "    k_coef=1e-4,\n",
                "    area_exp=0.4,\n",
                "    slope_exp=1,\n",
                "    tolerance=1e-5,\n",
                ")\n",
                "\n",
                "diffusion_eroder = fs.DiffusionADIEroder(grid, 0.01)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": ["## Setup Initial Conditions and External Forcing"]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source":
                ["Setup initial topography as two nearly flat plateaus (+ small random perturbations) separated by a steep escarpment. Also initialize the array for drainage area."]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "rng = np.random.Generator(np.random.PCG64(1234))\n",
                "\n",
                "init_elevation = rng.uniform(0, 1, size=grid.shape)\n",
                "init_elevation[:, 100:] += 400\n",
                "\n",
                "elevation = init_elevation\n",
                "drainage_area = np.empty_like(elevation)"
            ]
        },
        { "cell_type": "markdown", "metadata": {}, "source": ["## Run the Model"] },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source":
                ["Run the model for a few dozens of time steps (total simulation time: 100k years)."]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "dt = 2e3\n",
                "nsteps = 50\n",
                "\n",
                "for step in range(nsteps):\n",
                "    # flow routing\n",
                "    flow_graph.update_routes(elevation)\n",
                "    \n",
                "    # flow accumulation (drainage area)\n",
                "    flow_graph.accumulate(drainage_area, 1.0)\n",
                "    \n",
                "    # apply channel erosion then hillslope diffusion\n",
                "    spl_erosion = spl_eroder.erode(elevation, drainage_area, dt)\n",
                "    diff_erosion = diffusion_eroder.erode(elevation - spl_erosion, dt)\n",
                "    \n",
                "    # update topography\n",
                "    elevation = elevation - spl_erosion - diff_erosion\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": ["## Plot Outputs and Other Diagnostics\n"]
        },
        { "cell_type": "markdown", "metadata": {}, "source": ["- Topographic elevation"] },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig, ax = plt.subplots(figsize=(12, 6))\n",
                "plt.imshow(elevation)\n",
                "plt.colorbar();"
            ]
        },
        { "cell_type": "markdown", "metadata": {}, "source": ["- Drainage area (log)"] },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fig, ax = plt.subplots(figsize=(12, 6))\n",
                "plt.imshow(np.log(drainage_area), cmap=plt.cm.Blues)\n",
                "plt.colorbar();"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "kernelspec":
            { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" },
        "language_info": {
            "codemirror_mode": { "name": "ipython", "version": 3 },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.9.10"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 4
}