|
| 1 | +# ShallowWaters.jl - A type-flexible 16-bit shallow water model |
| 2 | +[](https://github.com/milankl/ShallowWaters.jl/actions/workflows/CI.yml) |
| 3 | +[](https://zenodo.org/badge/latestdoi/132787050) |
| 4 | + |
| 5 | + |
| 6 | +A shallow water model with a focus on type-flexibility and 16-bit number formats. ShallowWaters allows for Float64/32/16, |
| 7 | +[Posit32/16/8](https://github.com/milankl/SoftPosit.jl), [BFloat16](https://github.com/JuliaComputing/BFloat16s.jl), |
| 8 | +[LogFixPoint16](https://github.com/milankl/LogFixPoint16s.jl), [Sonum16](https://github.com/milankl/Sonums.jl), |
| 9 | +[Float32/16 & BFloat16 with stochastic rounding](https://github.com/milankl/StochasticRounding.jl) and in |
| 10 | +general every number format with arithmetics and conversions implemented. ShallowWaters also allows for |
| 11 | +mixed-precision and reduced precision communication. |
| 12 | + |
| 13 | +ShallowWaters uses an energy and enstrophy conserving advection scheme and a Smagorinsky-like biharmonic diffusion operator. |
| 14 | +Tracer advection is implemented with a semi-Lagrangian advection scheme. Strong stability-preserving Runge-Kutta schemes of |
| 15 | +various orders and stages are used with a semi-implicit treatment of the continuity equation. Boundary conditions are either |
| 16 | +periodic (only in x direction) or non-periodic super-slip, free-slip, partial-slip, or no-slip. |
| 17 | +Output via [NetCDF](https://github.com/JuliaGeo/NetCDF.jl). |
| 18 | + |
| 19 | +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. |
| 20 | + |
| 21 | +Requires: Julia 1.2 or higher |
| 22 | + |
| 23 | +## How to use |
| 24 | + |
| 25 | +`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`). |
| 26 | +```julia |
| 27 | +help?> Parameter |
| 28 | +search: Parameter |
| 29 | + |
| 30 | + Creates a Parameter struct with following options and default values |
| 31 | + |
| 32 | + T::DataType=Float32 # number format |
| 33 | + |
| 34 | + Tprog::DataType=T # number format for prognostic variables |
| 35 | + Tcomm::DataType=Tprog # number format for ghost-point copies |
| 36 | + |
| 37 | + # DOMAIN RESOLUTION AND RATIO |
| 38 | + nx::Int=100 # number of grid cells in x-direction |
| 39 | + Lx::Real=2000e3 # length of the domain in x-direction [m] |
| 40 | + L_ratio::Real=2 # Domain aspect ratio of Lx/Ly |
| 41 | + ... |
| 42 | +``` |
| 43 | +They can be changed with keyword arguments. The number format `T` is defined as the first (but optional) argument of `RunModel(T,...)` |
| 44 | +```julia |
| 45 | +julia> Prog = run_model(Float32,Ndays=10,g=10,H=500,Fx0=0.12); |
| 46 | +Starting ShallowWaters on Sun, 20 Oct 2019 19:58:25 without output. |
| 47 | +100% Integration done in 4.65s. |
| 48 | +``` |
| 49 | +or by creating a Parameter struct |
| 50 | +```julia |
| 51 | +julia> P = Parameter(bc="nonperiodic",wind_forcing_x="double_gyre",L_ratio=1,nx=128); |
| 52 | +julia> Prog = run_model(P); |
| 53 | +``` |
| 54 | +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. |
| 55 | + |
| 56 | +## Double-gyre example |
| 57 | + |
| 58 | +You can for example run a double gyre simulation like this |
| 59 | +```julia |
| 60 | +julia> using ShallowWaters |
| 61 | +julia> P = run_model(Ndays=100,nx=100,L_ratio=1,bc="nonperiodic",wind_forcing_x="double_gyre",topography="seamount"); |
| 62 | +Starting ShallowWaters on Sat, 15 Aug 2020 11:59:21 without output. |
| 63 | +100% Integration done in 13.7s. |
| 64 | +``` |
| 65 | +Sea surface height can be visualised via |
| 66 | +```julia |
| 67 | +julia> using PyPlot |
| 68 | +julia> pcolormesh(P.η') |
| 69 | +``` |
| 70 | + |
| 71 | + |
| 72 | +Or let's calculate the speed of the currents |
| 73 | +```julia |
| 74 | +julia> speed = sqrt.(Ix(P.u.^2)[:,2:end-1] + Iy(P.v.^2)[2:end-1,:]) |
| 75 | +``` |
| 76 | +`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. |
| 77 | +```julia |
| 78 | +julia> pcolormesh(speed') |
| 79 | +``` |
| 80 | + |
| 81 | + |
| 82 | +Such that the currents are strongest around the two eddies, as expected in this quasi-geostrophic setup. |
| 83 | + |
| 84 | +## (Some) Features |
| 85 | + |
| 86 | +- Interpolation of initial conditions from low resolution / high resolution runs. |
| 87 | +- Output of relative vorticity, potential vorticity and tendencies du,dv,deta |
| 88 | +- (Pretty accurate) duration estimate |
| 89 | +- Can be run in ensemble mode with ordered non-conflicting output files |
| 90 | +- Runs at CFL=1 (RK4), and more with the strong stability-preserving Runge-Kutta methods |
| 91 | +- Solving the tracer advection comes at basically no cost, thanks to semi-Lagrangian advection scheme |
| 92 | +- Also outputs the gradient operators ∂/∂x,∂/∂y and interpolations Ix, Iy for easier post-processing. |
| 93 | + |
| 94 | +## Installation |
| 95 | + |
| 96 | +ShallowWaters.jl is a registered package, so simply do |
| 97 | + |
| 98 | +```julia |
| 99 | +julia> ] add ShallowWaters |
| 100 | +``` |
| 101 | + |
| 102 | +## References |
| 103 | + |
| 104 | +ShallowWaters.jl was used and is described in more detail in |
| 105 | + |
| 106 | +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) |
| 107 | + |
| 108 | +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) |
| 109 | + |
| 110 | +## The equations |
| 111 | + |
| 112 | +The non-linear shallow water model plus tracer equation is |
| 113 | + |
| 114 | + ∂u/∂t + (u⃗⋅∇)u - f*v = -g*∂η/∂x - c_D*|u⃗|*u + ∇⋅ν*∇(∇²u) + Fx(x,y) (1) |
| 115 | + ∂v/∂t + (u⃗⋅∇)v + f*u = -g*∂η/∂y - c_D*|u⃗|*v + ∇⋅ν*∇(∇²v) + Fy(x,y) (2) |
| 116 | + ∂η/∂t = -∇⋅(u⃗h) + γ*(η_ref - η) + Fηt(t)*Fη(x,y) (3) |
| 117 | + ∂ϕ/∂t = -u⃗⋅∇ϕ (4) |
| 118 | + |
| 119 | +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. |
| 120 | + |
| 121 | +The linear shallow water model equivalent is |
| 122 | + |
| 123 | + ∂u/∂t - f*v = -g*∂η/∂x - r*u + ∇⋅ν*∇(∇²u) + Fx(x,y) (1) |
| 124 | + ∂v/∂t + f*u = -g*∂η/∂y - r*v + ∇⋅ν*∇(∇²v) + Fy(x,y) (2) |
| 125 | + ∂η/∂t = -H*∇⋅u⃗ + γ*(η_ref - η) + Fηt(t)*Fη(x,y) (3) |
| 126 | + ∂ϕ/∂t = -u⃗⋅∇ϕ (4) |
| 127 | + |
| 128 | +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. |
0 commit comments