{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "slide" } }, "outputs": [], "source": [ "using AIBECS\n", "using PyPlot, PyCall\n", "using LinearAlgebra\n", "using GR, Interact" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "

AIBECS.jl

\n", "\n", "*The ideal tool for exploring global marine biogeochemical cycles*\n", "\n", "**A**lgebraic **I**mplicit **B**iogeochemical **E**lemental **C**ycling **S**ystem\n", "\n", "Check it on GitHub (look for [AIBECS.jl](https://github.com/briochemc/AIBECS.jl))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", "\n", "
A Julia package" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Outline\n", "\n", "1. Motivation and concept\n", "1. Example: Phosphorus cycle in OCIM1" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Motivation: Starting from the AWESOME OCIM\n", "\n", "Features: \n", "\n", "\n", "- **GUI**\n", "- **simple to use**\n", "- **fast**\n", "- **good circulation**
(thanks to the OCIM1 by **Tim DeVries** (UCSB))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Motivation: comes the AIBECS\n", "\n", "\n", "\n", "Features (at present)\n", "\n", "\\- **no GUI** but **simple to use** and **fast**
\n", "\\- **nonlinear** systems
\n", "\\- **multiple tracers**
\n", "\\- **swappable circulations** (not just the OCIM)
\n", "\\- made for **parameter estimation/optimization** and **sensitivity analysis**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "# AIBECS Concept: a simple interface\n", "\n", "To build a BGC model with the AIBECS, you just need to\n", "\n", "**1.** Specify the **local sources and sinks**\n", "\n", "**2.** Chose the **ocean circulation**
\n", "\n", "(**3.** Specify the **particle transport**, if any)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "\n", "\n", "# AIBECS concept: Vectorization\n", "\n", "The **3D ocean grid** is rearranged
into a **1D column vector**.\n", "\n", "And **linear operators** are represented by **matrices**, allowing use of powerful linear algebra tools." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "\n", "

Example: A phosphorus-cycling model

\n", "\n", "Dissolved inorganic phosphrous (**DIP**)
\n", "(transported by the **ocean circulation**)\n", "\n", "$$\\left[\\frac{\\partial}{\\partial t} + \\color{RoyalBlue}{\\nabla \\cdot (\\boldsymbol{u} + \\mathbf{K}\\cdot\\nabla )}\\right] x_\\mathsf{DIP} = \\color{forestgreen}{-U(x_\\mathsf{DIP}) + R(x_\\mathsf{POP})},$$\n", "\n", "and particulate organic phosphrous (**POP**)
\n", "(transported by **sinking particles**)\n", "\n", "\n", "$$\\left[\\frac{\\partial}{\\partial t} + \\color{DarkOrchid}{\\nabla \\cdot \\boldsymbol{w}}\\right] x_\\mathsf{POP} = \\color{forestgreen}{U(x_\\mathsf{DIP}) - R(x_\\mathsf{POP})}.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Translating to AIBECS Code is “easy”\n", "\n", "To use AIBECS, we must recast each tracer equation for $\\boldsymbol{x}_k$ into the generic form:\n", "\n", "$$\\frac{\\partial \\boldsymbol{x}_k}{\\partial t} + \\color{RoyalBlue}{\\underbrace{\\mathbf{T}(\\boldsymbol{p})}_\\textrm{transport}} \\, \\boldsymbol{x} = \\color{ForestGreen}{\\underbrace{\\boldsymbol{G}(\\boldsymbol{x}_1, \\ldots, \\boldsymbol{x}_N, \\boldsymbol{p})}_\\textrm{local sources and sinks}}$$\n", "\n", "where $\\boldsymbol{p} =$ vector of model parameters\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "For **DIP**\n", "\n", "$$\\left[\\frac{\\partial}{\\partial t} + \\color{RoyalBlue}{\\nabla \\cdot (\\boldsymbol{u} + \\mathbf{K}\\cdot\\nabla )}\\right] x_\\mathsf{DIP} = \\color{forestgreen}{-U(x_\\mathsf{DIP}) + R(x_\\mathsf{POP})},$$\n", "\n", "becomes\n", "\n", "$$\\left[\\frac{\\partial}{\\partial t} + \\color{RoyalBlue}{\\mathbf{T}_{\\mathsf{DIP}}(\\boldsymbol{p})}\\right] \\boldsymbol{x}_\\mathsf{DIP} = \\color{forestgreen}{G_\\mathsf{DIP}(\\boldsymbol{x}_\\mathsf{DIP}, \\boldsymbol{x}_\\mathsf{POP}, \\boldsymbol{p})}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "For **POP**\n", "\n", "$$\\left[\\frac{\\partial}{\\partial t} + \\color{DarkOrchid}{\\nabla \\cdot \\boldsymbol{w}}\\right] x_\\mathsf{POP} = \\color{forestgreen}{U(x_\\mathsf{DIP}) - R(x_\\mathsf{POP})}.$$\n", "\n", "becomes\n", "\n", "$$\\left[\\frac{\\partial}{\\partial t} + \\color{RoyalBlue}{\\mathbf{T}_{\\mathsf{POP}}(\\boldsymbol{p})}\\right] \\boldsymbol{x}_\\mathsf{POP} = \\color{forestgreen}{G_\\mathsf{POP}(\\boldsymbol{x}_\\mathsf{DIP}, \\boldsymbol{x}_\\mathsf{POP}, \\boldsymbol{p})}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## The transport for DIP is the ocean circulation\n", "\n", "For DIP, the advective–eddy-diffusive transport operator, $\\nabla \\cdot (\\boldsymbol{u} + \\mathbf{K}\\cdot\\nabla)$, is converted into the **OCIM1 transport matrix** `T`, which can be loaded via the `OCIM1.load` function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "wet3D, grd, T = OCIM1.load();" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We assign `T` as the transport for DIP this way" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "T_DIP(p) = T" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "\n", "\n", "### Small diversion: What is this matrix doing?\n", "\n", "The matrix $\\mathbf{T}$ acts on the column vector of the tracers.\n", "\n", "What does `T` look like?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "T" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## The transport for sinking particles\n", "\n", "For POP, the **particle flux divergence (PFD)** operator, $\\nabla \\cdot \\boldsymbol{w}$, is created via `buildPFD`. We assign the transport for POP via" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "T_POP(p) = buildPFD(grd, wet3D, sinking_speed = w(p))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The settling velocity, `w(p)`, is assumed linearly increasing with depth `z` to yield a \"Martin curve profile\"" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "w(p) = p.w₀ .+ p.w′ * (z .|> ustrip)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "\n", "\n", "The vector of depths, `z`, is defined through `iwet`, which converts from 3D to 1D." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "iwet = findall(wet3D)\n", "z = grd.depth_3D[iwet] " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### Uptake:\n", "$$\\frac{1}{τ} \\, \\frac{\\boldsymbol{x}_\\mathrm{DIP}^2}{\\boldsymbol{x}_\\mathrm{DIP} + k}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "relu(x) = (x .≥ 0) .* x\n", "zₑ = 80u\"m\" # depth of the euphotic zone\n", "function uptake(DIP, p)\n", " τ, k = p.τ, p.k\n", " DIP⁺ = relu(DIP)\n", " return @. 1/τ * DIP⁺^2 / (DIP⁺ + k) * (z ≤ zₑ)\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$$\\frac{1}{τ} \\, \\frac{\\boldsymbol{x}_\\mathrm{DIP}^2}{\\boldsymbol{x}_\\mathrm{DIP} + k} = \\underbrace{\\frac{1}{τ} \\,\\frac{\\boldsymbol{x}_\\mathrm{DIP}}{\\boldsymbol{x}_\\mathrm{DIP} + k}}_{\\textrm{Michaelis-Menten}} \\cdot \\boldsymbol{x}_\\mathrm{DIP}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "function plot_uptake()\n", " PyPlot.figure(figsize=(8,4))\n", " k = 0.5\n", " τ = 5.0\n", " x = 0:0.01:1.5\n", " y = [1/τ*x^2/(x+k) for x in x]\n", " y2 = [1/τ*x for x in x]\n", " PyPlot.plot(x,y)\n", " PyPlot.plot(x,y2)\n", " PyPlot.plot(k*one.(x),y, \":\")\n", " PyPlot.xlabel(\"DIP (mol m⁻³)\")\n", " PyPlot.ylabel(\"uptake\")\n", " PyPlot.legend((\"uptake (mol m⁻³ s⁻¹)\", \"1/τ DIP\", \"k (=$k mol m⁻³)\"))\n", " PyPlot.title(\"Uptake rate\")\n", "end\n", "plot_uptake()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### Remineralization:\n", "$$\\kappa \\, \\boldsymbol{x}_\\mathrm{POP}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "remineralization(POP, p) = p.κ * POP" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### “Geological” restoring\n", "\n", "$$\\frac{x_\\mathrm{geo} - \\boldsymbol{x}_\\mathrm{DIP}}{\\tau_\\mathrm{geo}}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "geores(x, p) = (p.xgeo .- x) / p.τgeo" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Local sources and sinks\n", "\n", "### Net sum of local sources and sinks" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "G_DIP(DIP, POP, p) = -uptake(DIP, p) + remineralization(POP, p) + geores(DIP, p)\n", "G_POP(DIP, POP, p) = uptake(DIP, p) - remineralization(POP, p)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Model parameters\n", "\n", "First, create a table of parameters using the AIBECS interface" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "t = empty_parameter_table() # empty table of parameters\n", "add_parameter!(t, :xgeo, 2.12u\"mmol/m^3\")\n", "add_parameter!(t, :τgeo, 1.0u\"Myr\")\n", "add_parameter!(t, :k, 6.62u\"μmol/m^3\")\n", "add_parameter!(t, :w₀, 0.64u\"m/d\")\n", "add_parameter!(t, :w′, 0.13u\"1/d\")\n", "add_parameter!(t, :κ, 0.19u\"1/d\")\n", "add_parameter!(t, :τ, 236.52u\"d\")\n", "t" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We then create the parameters vector `p` via" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "initialize_Parameters_type(t, \"Pcycle_Parameters\") # Generate the parameter type\n", "p = Pcycle_Parameters()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# State function (and Jacobian)\n", "\n", "The \"state function\" $\\boldsymbol{F}$ is the **total** rate of change of the whole system (i.e., both DIP and POP) due to local sources and sinks **and** transport:\n", "\n", "$$\\frac{\\partial \\boldsymbol{x}}{\\partial t} = \\color{Brown}{\\boldsymbol{F}(\\boldsymbol{x}, \\boldsymbol{p})}$$\n", "\n", "For this coupled DIP–POP system, \n", "\n", "$$\\boldsymbol{F}\\left(\\begin{bmatrix}\\boldsymbol{x}_\\mathsf{DIP} \\\\ \\boldsymbol{x}_\\mathsf{POP}\\end{bmatrix}, \\boldsymbol{p}\\right) = \\begin{bmatrix}\\color{forestgreen}{\\boldsymbol{G}_\\mathsf{DIP}(\\boldsymbol{x}_\\mathsf{DIP}, \\boldsymbol{x}_\\mathsf{POP}, \\boldsymbol{p})} - \\color{royalblue}{\\boldsymbol{T}_\\mathsf{DIP}} \\, \\boldsymbol{x}_\\mathsf{DIP} \\\\ \\color{forestgreen}{\\boldsymbol{G}_\\mathsf{POP}(\\boldsymbol{x}_\\mathsf{DIP}, \\boldsymbol{x}_\\mathsf{POP}, \\boldsymbol{p})} - \\color{royalblue}{\\boldsymbol{T}_\\mathsf{POP}} \\, \\boldsymbol{x}_\\mathsf{POP} \\end{bmatrix}$$\n", "\n", "We generate `F` and `∇ₓF` via" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "nb = length(iwet)\n", "F, ∇ₓF = state_function_and_Jacobian((T_DIP,T_POP), (G_DIP,G_POP), nb) ;" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "skip" } }, "source": [ "## Why the Jacobian `∇ₓF`?\n", "\n", "The Jacobian matrix is $\\nabla_{\\boldsymbol{x}}\\boldsymbol{F}(\\boldsymbol{x},\\boldsymbol{p}) = \\left[\\frac{\\partial F_i}{\\partial x_j}\\right]_{i,j}$, is useful for \n", "- **implicit** time-steps\n", "- **solving** the **steady-state** system\n", "- **optimization** / **uncertainty analysis**\n", "\n", "With AIBECS, the **Jacobian** is **computed automatically** for you under the hood... using **dual numbers**!
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Solve directly for the steady state\n", "\n", "We can directly solve for the **steady-state**, $\\boldsymbol{s}$,
\n", "(using `CTKAlg()`, a quasi-Newton root-finding algorithm from C. T. Kelley)\n", "\n", "i.e., such that $\\boldsymbol{F}(\\boldsymbol{s},\\boldsymbol{p}) = 0$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "x = ones(2nb)\n", "prob = SteadyStateProblem(F, ∇ₓF, x, p)\n", "s = solve(prob, CTKAlg(), preprint=\" \").u" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Plot DIP" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "DIP, POP = state_to_tracers(s, nb, 2)\n", "DIP_unit = u\"mmol/m^3\"\n", "lon, lat = grd.lon |> ustrip, grd.lat |> ustrip\n", "function horizontal_slice(x, levels, title)\n", " x_3D = fill(NaN, size(grd))\n", " x_3D[iwet] .= x .|> ustrip\n", " GR.figure(figsize=(10,5))\n", " @manipulate for iz in 1:size(grd)[3]\n", " GR.xlim([0,360])\n", " GR.title(string(title, \" at $(AIBECS.round(grd.depth[iz])) depth\"))\n", " GR.contourf(lon, lat, x_3D[:,:,iz]', levels=levels)\n", " end\n", "end\n", "horizontal_slice(DIP * u\"mol/m^3\" .|> DIP_unit, 0:0.3:3, \"P-cycle model: DIP [mmol m^{-3}]\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Plot POP" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "function zonal_slice(x, levels, title)\n", " x_3D = fill(NaN, size(grd))\n", " x_3D[iwet] .= x .|> ustrip\n", " GR.figure(figsize=(10,5))\n", " @manipulate for longitude in (grd.lon .|> ustrip)\n", " GR.title(string(title, \" at $(round(longitude))°\"))\n", " ilon = findfirst(ustrip.(grd.lon) .== longitude)\n", " GR.contourf(lat, reverse(-grd.depth .|> ustrip), reverse(x_3D[:,ilon,:], dims=2), levels=levels)\n", " end\n", "end\n", "zonal_slice(POP .|> relu .|> log10, -10:1:10, \"P-cycle model: POP [log\\\\(mol m^{-3}\\\\)]\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "We can also make publication-quality plots using Python's Cartopy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "ccrs = pyimport(\"cartopy.crs\")\n", "cfeature = pyimport(\"cartopy.feature\")\n", "lon_cyc = [lon; 360+lon[1]]\n", "DIP_3D = fill(NaN, size(grd)); DIP_3D[iwet] = DIP * u\"mol/m^3\" .|> u\"mmol/m^3\" .|> ustrip\n", "DIP_3D_cyc = cat(DIP_3D, DIP_3D[:,[1],:], dims=2)\n", "f1 = PyPlot.figure(figsize=(12,5))\n", "@manipulate for iz in 1:24\n", " withfig(f1, clear=true) do\n", " ax = PyPlot.subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))\n", " ax.add_feature(cfeature.COASTLINE, edgecolor=\"#000000\") # black coast lines\n", " ax.add_feature(cfeature.LAND, facecolor=\"#CCCCCC\") # gray land\n", " plt = PyPlot.contourf(lon_cyc, lat, DIP_3D_cyc[:,:,iz], levels=0:0.1:3.5, \n", " transform=ccrs.PlateCarree(), zorder=-1, extend=\"both\")\n", " plt2 = PyPlot.contour(plt, lon_cyc, lat, DIP_3D_cyc[:,:,iz], levels=plt.levels[1:5:end], \n", " transform=ccrs.PlateCarree(), zorder=0, extend=\"both\", colors=\"k\")\n", " ax.clabel(plt2, fmt=\"%2.1f\", colors=\"k\", fontsize=14)\n", " cbar = PyPlot.colorbar(plt, orientation=\"vertical\", extend=\"both\", ticks=plt2.levels)\n", " cbar.add_lines(plt2)\n", " cbar.ax.tick_params(labelsize=14)\n", " cbar.set_label(label=\"mmol / m³\", fontsize=16)\n", " PyPlot.title(\"Publication-quality DIP with Cartopy; depth = $(string(round(typeof(1u\"m\"),grd.depth[iz])))\", fontsize=20)\n", " end\n", "end" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.2.0", "language": "julia", "name": "julia-1.2" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.2.0" }, "rise": { "scroll": true } }, "nbformat": 4, "nbformat_minor": 4 }