Add docs

.github/workflows/Documentation.yml

@@ -0,0 +1,23 @@
+name: Documentation
+on:
+  push:
+    branches:
+      - main
+    tags: '*'
+  pull_request:
+    types: [opened, synchronize, reopened]
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: '1'
+      - name: Install dependencies
+        run: julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+        run: julia --project=docs --color=yes docs/make.jl

docs/Project.toml

@@ -0,0 +1,7 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+
+[compat]
+Documenter = "0.27"
+Literate = "2"

docs/assets/swm_equations.jl

@@ -0,0 +1,110 @@
+# # ShallowWaters.jl - The equations
+#
+# The shallow water equations for the prognostic variables velocity $\mathbf{u} = (u,v)$ and sea surface elevation $\eta$ over the 2-dimensional domain $\Omega$ in $x,y$ are
+#
+# ```math
+# \begin{align}
+# \partial_t u &+ u\partial_xu + v\partial_yu - fv = -g\partial_x\eta + D_x(u,v,\eta) + F_x, \\
+# \partial_t v &+ u\partial_xv + v\partial_yv + fu = -g\partial_y\eta + D_y(u,v,\eta) + F_y, \\
+# \partial_t \eta &+ \partial_x(uh) + \partial_y(vh) = 0.
+# \end{align}
+# ```
+# where the first two are the momentum equations for $u,v$ and the latter is the continuity equation. The layer thickness is $h = \eta + H$ with $H=H(x,y)$ being the bottom topography. The gravitational acceleration is $g$, the coriolis parameter $f = f(y)$ depends only (i.e. latitude) only and the beta-plane approximation is used. $(D_x,D_y) = \mathbf{D}$ are the dissipative terms
+# 
+# ```math
+# \mathbf{D} = -\frac{c_D}{h}\vert \mathbf{u} \vert \mathbf{u} - u \nabla^4\mathbf{u}
+# ```
+
+# which are a sum of a quadratic bottom drag with dimensionless coefficient $c_D$ and a biharmonic diffusion with viscosity coefficient $u$. $\mathbf{F} = (F_x,F_y)$ is the wind forcing which can depend on time and space, i.e. $F_x = F_x(x,y,t)$ and $F_y = F_y(x,y,t)$.
+
+# ## The vector invariant formulation
+
+# The Bernoulli potential $p$ is introduced as
+
+# ```math
+# p = \frac{1}{2}(u^2 + v^2) + gh
+# ```
+
+# The relative vorticity $\zeta = \partial_xv + \partial_yu$ lets us define the potential vorticity $q$ as
+
+# ```math
+# q = \frac{f + \zeta}{h}
+# ```
+
+# such that we can rewrite the shallow water equations as
+
+# ```math
+# \begin{align}
+# \partial_t u &= qhv -g\partial_xp + D_x(u,v,\eta) + F_x, \\
+# \partial_t v &= -qhu -g\partial_yp + D_y(u,v,\eta) + F_y, \\ 
+# \partial_t \eta &= -\partial_x(uh) -\partial_y(vh).
+# \end{align}
+# ```
+
+# ## Runge-Kutta time discretisation
+
+# Let
+# ```math
+# R(u,v,\eta) = \begin{pmatrix} 
+#                 qhv-g\partial_xp+F_x 
+#                 -qhu-g\partial_yp+F_y 
+#                 -\partial_x(uh) -\partial_y(vh)
+#                 \end{pmatrix}
+# ```
+# 
+# be the non-dissipative right-hand side, i.e. excluding the dissipative terms $\mathbf{D}$. Then we dicretise the time derivative with 4th order Runge-Kutta with $\mathbf{k}_n = (u_n,v_n,\eta_n)$ by
+# 
+# ```math
+# \begin{align}
+# \mathbf{d}_1 &= R(\mathbf{k}_n) , \\
+# \mathbf{d}_2 &= R(\mathbf{k}_n + \frac{1}{2}\Delta t \mathbf{d}_1), \\
+# \mathbf{d}_3 &= R(\mathbf{k}_n + \frac{1}{2}\Delta t \mathbf{d}_2), \\
+# \mathbf{d}_4 &= R(\mathbf{k}_n + \Delta t \mathbf{d}_3), \\
+# u_{n+1}^*,v_{n+1}^*,\eta_{n+1} = \mathbf{k}_{n+1} &= \mathbf{k}_n + \frac{1}{6}\Delta t(\mathbf{d}_1 + 2\mathbf{d}_2 + 2\mathbf{d}_3 + \mathbf{d}_1),
+# \end{align}
+# ```
+
+# and the dissipative terms are then added semi-implictly.
+
+# ```math
+# \begin{align}
+# u_{n+1} = u_{n+1}^* + \Delta t D_x(u_{n+1}^*,v_{n+1}^*,\eta_{n+1}) , \\
+# v_{n+1} = v_{n+1}^* + \Delta t D_y(u_{n+1}^*,v_{n+1}^*,\eta_{n+1}).
+# \end{align}
+# ```
+
+
+# Consequently, the dissipative terms only have to be evaluated once per time step, which reduces the computational cost of the right=hand side drastically. This is motivated as the Courant number $C = \sqrt{gH_0}$ is restricted by gravity waves, which are caused by $\partial_t\mathbf{u} = -g\nabla\eta$ and $\partial_t\eta = -H_0\nabla \cdot \mathbf{u}$. The other terms have longer time scales, and it is therefore sufficient to solve those with larger time steps. With this scheme, the shallow water model runs stable at $C=1$."
+
+# ## Strong stability preserving Runge-Kutta with semi-implicit continuity equation
+
+# We split the right-hand side into the momentum equations and the continuity equation
+
+# ```math
+# R_m(u,v,\eta) = \begin{pmatrix} 
+#                 qhv-g\partial_xp+F_x 
+#                 -qhu-g\partial_yp+F_y 
+#                 \end{pmatrix}, \quad R_\eta(u,v,\eta) = -\partial_x(uh) -\partial_y(vh)
+# ```
+
+# The 4-stage strong stability preserving Runge-Kutta scheme, with a semi-implicit treatment of the continuity equation then reads as
+
+# ```math
+# \begin{align}
+# \mathbf{u}_1 &= \mathbf{u}_n + \frac{1}{2} \Delta t R_m(u_n,v_n,\eta_n), \quad \text{then}
+# \quad \eta_1 = \eta_n + \frac{1}{2} \Delta t R_\eta(u_1,v_1,\eta_n), \\
+# \mathbf{u}_2 &= \mathbf{u}_1 + \frac{1}{2} \Delta t R_m(u_1,v_1,\eta_1), \quad \text{then}
+# \quad \eta_2 = \eta_1 + \frac{1}{2} \Delta t R_\eta(u_2,v_2,\eta_1) , \\
+# \mathbf{u}_3 &= \frac{2}{3}\mathbf{u}_n + \frac{1}{3}\mathbf{u}_2 + \frac{1}{6} \Delta t R_m(u_2,v_2,\eta_2), \quad \text{then}
+# \quad \eta_3 = \frac{2}{3}\eta_n + \frac{1}{3}\eta_2 + \frac{1}{6} \Delta t R_\eta(u_3,v_3,\eta_2), \\
+# \mathbf{u}_{n+1} &= \mathbf{u}_3 + \frac{1}{2} \Delta t R_m(u_3,v_3,\eta_3), \quad \text{then}
+# \quad \eta_{n+1} = \eta_3 + \frac{1}{2} \Delta t R_\eta(u_{n+1},v_{n+1},\eta_3).
+# \end{align}
+# ```
+
+# ## Splitting the continuity equation
+
+# From
+# ```math
+# \partial_t = -\partial_x(uh) - \partial_y(vh)
+# ```

docs/make.jl

@@ -0,0 +1,34 @@
+using Documenter
+using Literate
+using ShallowWaters
+
+EXAMPLE = joinpath(@__DIR__, "assets", "swm_equations.jl")
+OUTPUT = joinpath(@__DIR__, "src")
+
+# Generate markdown
+binder_badge = "# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/milankl/ShallowWaters.jl/gh-pages?labpath=dev%2Fswm_equations.ipynb)"
+function preprocess_docs(content)
+  return string(binder_badge, "\n\n", content)
+end
+
+Literate.markdown(EXAMPLE, OUTPUT; preprocess=preprocess_docs, codefence="```julia" => "```")
+Literate.notebook(EXAMPLE, OUTPUT)
+
+pages = [
+  "Introduction" => "index.md",
+  "Tutorial" => "swm_equations.md",
+  "Reference" => "reference.md",
+]
+
+makedocs(
+  sitename = "ShallowWaters.jl",
+  format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
+  modules = [ShallowWaters],
+  pages = pages,
+)
+
+deploydocs(
+  repo = "github.com/milankl/ShallowWaters.jl.git",
+  push_preview = true,
+  devbranch = "main",
+)

docs/src/index.md

@@ -0,0 +1,128 @@
+# ShallowWaters.jl - A type-flexible 16-bit shallow water model
+[![CI](https://github.com/milankl/ShallowWaters.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/milankl/ShallowWaters.jl/actions/workflows/CI.yml)
+[![DOI](https://zenodo.org/badge/132787050.svg)](https://zenodo.org/badge/latestdoi/132787050)
+![sst](https://github.com/milankl/ShallowWaters.jl/tree/main/figs/isambard_float16.png?raw=true "Float16 simulation with ShallowWaters.jl on Isambard's A64FX")
+
+A shallow water model with a focus on type-flexibility and 16-bit number formats. ShallowWaters allows for Float64/32/16, 
+[Posit32/16/8](https://github.com/milankl/SoftPosit.jl), [BFloat16](https://github.com/JuliaComputing/BFloat16s.jl), 
+[LogFixPoint16](https://github.com/milankl/LogFixPoint16s.jl), [Sonum16](https://github.com/milankl/Sonums.jl), 
+[Float32/16 & BFloat16 with stochastic rounding](https://github.com/milankl/StochasticRounding.jl) and in 
+general every number format with arithmetics and conversions implemented. ShallowWaters also allows for
+mixed-precision and reduced precision communication.
+
+ShallowWaters uses an energy and enstrophy conserving advection scheme and a Smagorinsky-like biharmonic diffusion operator. 
+Tracer advection is implemented with a semi-Lagrangian advection scheme. Strong stability-preserving Runge-Kutta schemes of
+various orders and stages are used with a semi-implicit treatment of the continuity equation. Boundary conditions are either 
+periodic (only in x direction) or non-periodic super-slip, free-slip, partial-slip, or no-slip.
+Output via [NetCDF](https://github.com/JuliaGeo/NetCDF.jl).
+
+Please feel free to raise an [issue](https://github.com/milankl/ShallowWaters.jl/issues) if you discover bugs or have an idea how to improve ShallowWaters.
+
+Requires: Julia 1.2 or higher
+
+## How to use
+
+`RunModel` initialises the model, preallocates memory and starts the time integration. You find the options and default parameters in `src/DefaultParameters.jl` (or by typing `?Parameter`).
+```julia
+help?> Parameter
+search: Parameter
+
+  Creates a Parameter struct with following options and default values
+
+  T::DataType=Float32                 # number format
+
+  Tprog::DataType=T                   # number format for prognostic variables
+  Tcomm::DataType=Tprog               # number format for ghost-point copies
+
+  # DOMAIN RESOLUTION AND RATIO
+  nx::Int=100                         # number of grid cells in x-direction
+  Lx::Real=2000e3                     # length of the domain in x-direction [m]
+  L_ratio::Real=2                     # Domain aspect ratio of Lx/Ly
+  ...
+```
+They can be changed with keyword arguments. The number format `T` is defined as the first (but optional) argument of `RunModel(T,...)`
+```julia
+julia> Prog = run_model(Float32,Ndays=10,g=10,H=500,Fx0=0.12);
+Starting ShallowWaters on Sun, 20 Oct 2019 19:58:25 without output.
+100% Integration done in 4.65s.
+```
+or by creating a Parameter struct
+```julia
+julia> P = Parameter(bc="nonperiodic",wind_forcing_x="double_gyre",L_ratio=1,nx=128);
+julia> Prog = run_model(P);
+```
+The number formats can be different (aka mixed-precision) for different parts of the model. `Tprog` is the number type for the prognostic variables, `Tcomm` is used for communication of boundary values.
+
+## Double-gyre example
+
+You can for example run a double gyre simulation like this
+```julia
+julia> using ShallowWaters
+julia> P = run_model(Ndays=100,nx=100,L_ratio=1,bc="nonperiodic",wind_forcing_x="double_gyre",topography="seamount");
+Starting ShallowWaters on Sat, 15 Aug 2020 11:59:21 without output.
+100% Integration done in 13.7s.
+```
+Sea surface height can be visualised via
+```julia
+julia> using PyPlot
+julia> pcolormesh(P.η')
+```
+![Figure_1](https://user-images.githubusercontent.com/25530332/90311163-1ee40a00-def0-11ea-8911-810d7762cd3f.png)
+
+Or let's calculate the speed of the currents
+```julia
+julia> speed = sqrt.(Ix(P.u.^2)[:,2:end-1] + Iy(P.v.^2)[2:end-1,:])
+```
+`P.u` and `P.v` are the u,v velocity components on the Arakawa C-grid. To add them, we need to interpolate them with `Ix,Iy` (which are exported by `ShallowWaters.jl` too), then chopping off the edges to get two arrays of the same size.
+```julia
+julia> pcolormesh(speed')
+```
+![Figure_2](https://user-images.githubusercontent.com/25530332/90311211-88fcaf00-def0-11ea-8308-b4f438495152.png)
+
+Such that the currents are strongest around the two eddies, as expected in this quasi-geostrophic setup.
+
+## (Some) Features
+
+- Interpolation of initial conditions from low resolution / high resolution runs.
+- Output of relative vorticity, potential vorticity and tendencies du,dv,deta
+- (Pretty accurate) duration estimate
+- Can be run in ensemble mode with ordered non-conflicting output files
+- Runs at CFL=1 (RK4), and more with the strong stability-preserving Runge-Kutta methods
+- Solving the tracer advection comes at basically no cost, thanks to semi-Lagrangian advection scheme
+- Also outputs the gradient operators ∂/∂x,∂/∂y and interpolations Ix, Iy for easier post-processing.
+
+## Installation
+
+ShallowWaters.jl is a registered package, so simply do
+
+```julia
+julia> ] add ShallowWaters
+```
+
+## References
+
+ShallowWaters.jl was used and is described in more detail in  
+
+Klöwer M, Düben PD, Palmer TN. Number formats, error mitigation and scope for 16-bit arithmetics in weather and climate modelling analysed with a shallow water model. Journal of Advances in Modeling Earth Systems. doi: [10.1029/2020MS002246](https://dx.doi.org/10.1029/2020MS002246)
+
+Klöwer M, Düben PD, Palmer TN. Posits as an alternative to floats for weather and climate models. In: Proceedings of the Conference for Next Generation Arithmetic 2019. doi: [10.1145/3316279.3316281](https://dx.doi.org/10.1145/3316279.3316281)
+
+## The equations
+
+The non-linear shallow water model plus tracer equation is
+
+          ∂u/∂t + (u⃗⋅∇)u - f*v = -g*∂η/∂x - c_D*|u⃗|*u + ∇⋅ν*∇(∇²u) + Fx(x,y)     (1)
+          ∂v/∂t + (u⃗⋅∇)v + f*u = -g*∂η/∂y - c_D*|u⃗|*v + ∇⋅ν*∇(∇²v) + Fy(x,y)     (2)
+          ∂η/∂t = -∇⋅(u⃗h) + γ*(η_ref - η) + Fηt(t)*Fη(x,y)                       (3)
+          ∂ϕ/∂t = -u⃗⋅∇ϕ                                                          (4)
+
+with the prognostic variables velocity u⃗ = (u,v) and sea surface heigth η. The layer thickness is h = η + H(x,y). The Coriolis parameter is f = f₀ + βy with beta-plane approximation. The graviational acceleration is g. Bottom friction is either quadratic with drag coefficient c_D or linear with inverse time scale r. Diffusion is realized with a biharmonic diffusion operator, with either a constant viscosity coefficient ν, or a Smagorinsky-like coefficient that scales as ν = c_Smag*|D|, with deformation rate |D| = √((∂u/∂x - ∂v/∂y)² + (∂u/∂y + ∂v/∂x)²). Wind forcing Fx is constant in time, but may vary in space.
+
+The linear shallow water model equivalent is
+
+          ∂u/∂t - f*v = -g*∂η/∂x - r*u + ∇⋅ν*∇(∇²u) + Fx(x,y)     (1)
+          ∂v/∂t + f*u = -g*∂η/∂y - r*v + ∇⋅ν*∇(∇²v) + Fy(x,y)     (2)
+          ∂η/∂t = -H*∇⋅u⃗ + γ*(η_ref - η) + Fηt(t)*Fη(x,y)         (3)
+          ∂ϕ/∂t = -u⃗⋅∇ϕ                                           (4)
+
+ShallowWaters.jl discretises the equation on an equi-distant Arakawa C-grid, with 2nd order finite-difference operators. Boundary conditions are implemented via a ghost-point copy and each variable has a halo of variable size to account for different stencil sizes of various operators.

docs/src/reference.md

@@ -0,0 +1,17 @@
+# Reference
+​
+## Contents
+​
+```@contents
+Pages = ["reference.md"]
+```
+​
+## Index
+​
+```@index
+Pages = ["reference.md"]
+```
+​
+```@autodocs
+Modules = [ShallowWaters]
+```