spin-up.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d391c78a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import rasterio\n",
    "import firedrake\n",
    "from firedrake import Constant, sqrt, inner, grad, dx, ds\n",
    "import icepack\n",
    "from icepack.constants import ice_density as ρ_I, gravity as g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0c5f7e3",
   "metadata": {},
   "source": [
    "Load in the work from the previous notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "112b9091",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"modern_state.h5\"\n",
    "with firedrake.CheckpointFile(filename, \"r\") as chk:\n",
    "    mesh = chk.load_mesh()\n",
    "    C = chk.load_function(mesh, \"friction\")\n",
    "    u_0 = chk.load_function(mesh, \"velocity\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2de525e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = mesh.coordinates.dat.data_ro[:]\n",
    "delta = 10e3\n",
    "xmin, xmax = coords[:, 0].min() - delta, coords[:, 0].max() + delta\n",
    "ymin, ymax = coords[:, 1].min() - delta, coords[:, 1].max() + delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f603dc5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "image_filename = icepack.datasets.fetch_mosaic_of_antarctica()\n",
    "with rasterio.open(image_filename, \"r\") as image_file:\n",
    "    height, width = image_file.height, image_file.width\n",
    "    transform = image_file.transform\n",
    "    window = rasterio.windows.from_bounds(\n",
    "        left=xmin,\n",
    "        bottom=ymin,\n",
    "        right=xmax,\n",
    "        top=ymax,\n",
    "        transform=transform,\n",
    "    )\n",
    "    image = image_file.read(indexes=1, window=window, masked=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcadd596",
   "metadata": {},
   "outputs": [],
   "source": [
    "def subplots(*args, **kwargs):\n",
    "    fig, axes = plt.subplots()\n",
    "    axes.set_aspect(\"equal\")\n",
    "    xmin, ymin, xmax, ymax = rasterio.windows.bounds(window, transform)\n",
    "    axes.imshow(\n",
    "        image,\n",
    "        cmap=\"Greys_r\",\n",
    "        vmin=12e3,\n",
    "        vmax=16.38e3,\n",
    "        extent=(xmin, xmax, ymin, ymax),\n",
    "    )\n",
    "\n",
    "    return fig, axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f16ba37",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = firedrake.FunctionSpace(mesh, family=\"CG\", degree=2)\n",
    "V = firedrake.VectorFunctionSpace(mesh, family=\"CG\", degree=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67a2109c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bedmachine_filename = icepack.datasets.fetch_bedmachine_antarctica()\n",
    "thickness_filename = f\"netcdf:{bedmachine_filename}:thickness\"\n",
    "with rasterio.open(thickness_filename, \"r\") as thickness_file:\n",
    "    h_0 = icepack.interpolate(thickness_file, Q)\n",
    "\n",
    "surface_filename = f\"netcdf:{bedmachine_filename}:surface\"\n",
    "with rasterio.open(surface_filename, \"r\") as surface_file:\n",
    "    s_0 = icepack.interpolate(surface_file, Q)\n",
    "    \n",
    "s = s_0.copy(deepcopy=True)\n",
    "h = h_0.copy(deepcopy=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75a72474",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = firedrake.interpolate(s_0 - h_0, Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17eb1c8e",
   "metadata": {},
   "source": [
    "We'll make the accumulation rate a function of the elevation.\n",
    "Remember in the last notebook we found that an accumulation rate of about 10 cm / year overall would keep the system in steady state.\n",
    "First, let's compute the range of surface elevations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03351da5",
   "metadata": {},
   "outputs": [],
   "source": [
    "s_min = Constant(np.floor(s_0.dat.data_ro.min()))\n",
    "s_max = Constant(np.ceil(s_0.dat.data_ro.max()))\n",
    "print(f\"Elevation range: {float(s_min):.0f} - {float(s_max):.0f} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92cccba2",
   "metadata": {},
   "source": [
    "The MATLAB file of the 1D model's output shows min / max accumulation values for the modern as about 12-20 cm / year and 6-10 cm / year for the LGM.\n",
    "What I'm doing below is just for illustrative purposes.\n",
    "I'm making the accumulation scale linearly from 6 cm / year at low elevations to 10 cm / year at the highest ones and clamping them at either end.\n",
    "This doesn't include the effect where accumulation rate goes back down at really high elevations.\n",
    "We can include that (or other things) later if you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e452def",
   "metadata": {},
   "outputs": [],
   "source": [
    "from firedrake import min_value, max_value\n",
    "\n",
    "a_min = Constant(0.06)\n",
    "a_max = Constant(0.10)\n",
    "\n",
    "def accumulation(s):\n",
    "    λ = (s - s_min) / (s_max - s_min)\n",
    "    a = (1 - λ) * a_min + λ * a_max\n",
    "    return min_value(a_max, max_value(a_min, a))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ec1c941",
   "metadata": {},
   "source": [
    "This is just a function that tells us how to calculate the accumulation rate.\n",
    "The code below carries out that calculation and interpolates the results to a field defined over the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4d8693b",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = firedrake.interpolate(accumulation(s_0), Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2886101f",
   "metadata": {},
   "source": [
    "We'll use the same values for the temperature as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a88af56",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = Constant(260.0)\n",
    "A = icepack.rate_factor(T)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc009785",
   "metadata": {},
   "source": [
    "But we'll do something a little tricky with the terminus stress.\n",
    "When terminus BCs are applied, the model assumes that this boundary is adjacent to ocean water.\n",
    "In our case, there's more ice downstream, and we want to apply backpressure consistent with the thickness of this downstream ice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d8d5e7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def terminus(**kwargs):\n",
    "    u = kwargs[\"velocity\"]\n",
    "    h = kwargs[\"thickness\"]\n",
    "    s = kwargs[\"surface\"]\n",
    "    h_D = kwargs[\"thickness_downstream\"]\n",
    "    \n",
    "    τ_I = 0.5 * ρ_I * g * h**2\n",
    "    τ_D = 0.5 * ρ_I * g * h_D**2\n",
    "\n",
    "    ν = firedrake.FacetNormal(mesh)\n",
    "    return (τ_I - τ_D) * inner(u, ν)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "266781d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = icepack.models.IceStream(terminus=terminus)\n",
    "opts = {\n",
    "    \"dirichlet_ids\": [1, 2, 3, 5, 6, 7],\n",
    "    \"diagnostic_solver_type\": \"petsc\",\n",
    "    \"diagnostic_solver_parameters\": {\n",
    "        \"snes_max_it\": 100,\n",
    "        \"snes_type\": \"newtontr\",\n",
    "        \"ksp_type\": \"gmres\",\n",
    "        \"pc_type\": \"lu\",\n",
    "        \"pc_factor_mat_solver_type\": \"mumps\",\n",
    "    },\n",
    "}\n",
    "solver = icepack.solvers.FlowSolver(model, **opts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6812dca",
   "metadata": {},
   "outputs": [],
   "source": [
    "h_D = h_0.copy(deepcopy=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d35bceb",
   "metadata": {},
   "source": [
    "Having modified the terminus pressure calculation, we just pass whatever value we want for the downstream thickness.\n",
    "We're keeping it constant for now but you can make it change in time if you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef7439d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "u = solver.diagnostic_solve(\n",
    "    velocity=u_0,\n",
    "    thickness=h,\n",
    "    surface=s,\n",
    "    fluidity=A,\n",
    "    friction=C,\n",
    "    thickness_downstream=h_D,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3bb2bba",
   "metadata": {},
   "source": [
    "Now let's actually try to simulate things for a while.\n",
    "I started trying this with a timestep of 1 year.\n",
    "Things can explode if you set that value too high, so I worked it down to a timestep of 2 months by trial and error.\n",
    "After the system propagates out some initial transients, we'll dial this back up to a larger value for efficiency's sake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed3a54b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "hs = [h.copy(deepcopy=True)]\n",
    "ss = [s.copy(deepcopy=True)]\n",
    "us = [u.copy(deepcopy=True)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a0df735",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_time = 40.0\n",
    "dt = 1 / 6\n",
    "num_steps = int(final_time / dt)\n",
    "\n",
    "for step in tqdm.trange(num_steps):\n",
    "    a = firedrake.interpolate(accumulation(s), Q)\n",
    "    \n",
    "    h = solver.prognostic_solve(\n",
    "        dt,\n",
    "        thickness=h,\n",
    "        velocity=u,\n",
    "        accumulation=a,\n",
    "        thickness_inflow=h_0,\n",
    "    )\n",
    "    h.interpolate(firedrake.max_value(1.0, h))\n",
    "\n",
    "    s = icepack.compute_surface(thickness=h, bed=b)\n",
    "\n",
    "    u = solver.diagnostic_solve(\n",
    "        velocity=u,\n",
    "        thickness=h,\n",
    "        surface=s,\n",
    "        fluidity=A,\n",
    "        friction=C,\n",
    "        thickness_downstream=h_D,\n",
    "    )\n",
    "    \n",
    "    hs.append(h.copy(deepcopy=True))\n",
    "    ss.append(s.copy(deepcopy=True))\n",
    "    us.append(u.copy(deepcopy=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d08019ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "δh = firedrake.interpolate(h - h_0, Q)\n",
    "\n",
    "fig, axes = subplots()\n",
    "colors = firedrake.tripcolor(δh, vmin=-100, vmax=+100, cmap=\"RdBu_r\", axes=axes)\n",
    "fig.colorbar(colors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de245eed",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = subplots()\n",
    "colors = firedrake.tripcolor(u, axes=axes)\n",
    "fig.colorbar(colors);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38ea1ffc",
   "metadata": {},
   "source": [
    "On average, the ice thickened by about 5.5m over the whole domain, and about 6m at the grounding line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c04f35a",
   "metadata": {},
   "outputs": [],
   "source": [
    "area = firedrake.assemble(Constant(1) * dx(mesh))\n",
    "avg_δh = firedrake.assemble(δh * dx) / area\n",
    "print(f\"Average thickness change: {avg_δh:5.2f} m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49a54579",
   "metadata": {},
   "outputs": [],
   "source": [
    "length = firedrake.assemble(Constant(1) * ds(domain=mesh, subdomain_id=(4,)))\n",
    "terminus_ids = (4,)\n",
    "avg_δh_terminus = firedrake.assemble(δh * ds(terminus_ids)) / length\n",
    "print(f\"Average thickness change at terminus: {avg_δh_terminus:5.2f} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bed5b883",
   "metadata": {},
   "source": [
    "The weird transients have mostly propagated out.\n",
    "There's a blob to grid east of the one nunatak, which I'm guessing is because of an incorrectly specified side wall friction.\n",
    "We can fix that later if it matters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d8da7a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from firedrake import div\n",
    "\n",
    "qs = [firedrake.project(div(h * u), Q) for h, u in zip(hs, us)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7825f1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = subplots()\n",
    "colors = firedrake.tripcolor(qs[-1], vmin=-1.0, vmax=+1.0, cmap=\"RdBu\", axes=axes)\n",
    "fig.colorbar(colors);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1831ee4d",
   "metadata": {},
   "source": [
    "The plot below shows the maximum flux divergence over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25d7b2b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "qmaxs = np.array([abs(q.dat.data_ro).max() for q in qs])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01435594",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots()\n",
    "ts = np.linspace(0.0, final_time, num_steps + 1)\n",
    "axes.set_xlabel(\"years\")\n",
    "axes.set_ylabel(\"meters/year\")\n",
    "axes.plot(ts, qmaxs);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "421249af",
   "metadata": {},
   "source": [
    "Finally, we'll save the results to disk and continue in another notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ead70157",
   "metadata": {},
   "outputs": [],
   "source": [
    "with firedrake.CheckpointFile(\"spin_up_state.h5\", \"w\") as chk:\n",
    "    chk.save_function(C, name=\"friction\")\n",
    "    chk.save_function(h, name=\"thickness\")\n",
    "    chk.save_function(s, name=\"surface\")\n",
    "    chk.save_function(u, name=\"velocity\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "firedrake",
   "language": "python",
   "name": "firedrake"
  },
  "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}