Project Status: Active – The project has reached a stable, usable state and is being actively developed.TravisAppVeyorCirrus CI

DOI

ShallowWaters.jl - A type-flexible 16bit shallow water model

sst

A shallow water model with a focus on type-flexibility and 16bit number formats. ShallowWaters allows for Float64/32/16, BigFloat/ArbFloat/DoubleFloats, Posit32/16/8, BFloat16, Sonum16 and in general every number format with arithmetics and conversions implemented. ShallowWaters also allows for mixed-precision and reduced precision communication.

ShallowWaters is fully-explicit with an energy and enstrophy conserving advection scheme and a Smagorinsky-like biharmonic diffusion operator. Tracer advection is implemented with a semi-Lagrangian advection scheme. Runge-Kutta 4th-order is used for pressure, advective and Coriolis terms and the continuity equation. Semi-implicit time stepping for diffusion and bottom friction. Boundary conditions are either periodic (only in x direction) or non-periodic super-slip, free-slip, partial-slip, or no-slip. Output via NetCDF.

Please feel free to raise an issue if you discover bugs or have an idea how to improve ShallowWaters.

Requires: Julia 1.2

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).

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> Prog = RunModel(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> P = Parameter(bc="nonperiodic",wind_forcing_x="double_gyre",L_ratio=1,nx=128);
julia> Prog = RunModel(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.

(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
  • 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> ] add ShallowWaters

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.

ShallowWaters.jl splits the time steps for various terms: Runge Kutta 4th order scheme for the fast varying terms. The diffusive terms (bottom friction and diffusion) are solved semi-implicitly every n-th time step. The tracer equation is solved with a semi-Lagrangian scheme that uses much larger time steps.