Skip to content

Commit ef6181e

Browse files
authored
Merge pull request #145 from milankl/ssprk
strong-stability preserving R2/3 methods
2 parents 68b7f9c + a9fd9f2 commit ef6181e

File tree

10 files changed

+562
-115
lines changed

10 files changed

+562
-115
lines changed

docs/swm_equations.ipynb

Lines changed: 257 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,257 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# ShallowWaters.jl - The equations"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"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"
15+
]
16+
},
17+
{
18+
"cell_type": "markdown",
19+
"metadata": {},
20+
"source": [
21+
"$$\n",
22+
"\\begin{align}\n",
23+
"\\partial_t u &+ u\\partial_xu + v\\partial_yu - fv = -g\\partial_x\\eta + D_x(u,v,\\eta) + F_x \\\\\n",
24+
"\\partial_t v &+ u\\partial_xv + v\\partial_yv + fu = -g\\partial_y\\eta + D_y(u,v,\\eta) + F_y \\\\\n",
25+
"\\partial_t \\eta &+ \\partial_x(uh) + \\partial_y(vh) = 0.\n",
26+
"\\end{align}\n",
27+
"$$"
28+
]
29+
},
30+
{
31+
"cell_type": "markdown",
32+
"metadata": {},
33+
"source": [
34+
"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"
35+
]
36+
},
37+
{
38+
"cell_type": "markdown",
39+
"metadata": {},
40+
"source": [
41+
"$$\n",
42+
"\\mathbf{D} = -\\frac{c_D}{h}\\vert \\mathbf{u} \\vert \\mathbf{u} - \\nu \\nabla^4\\mathbf{u}\n",
43+
"$$"
44+
]
45+
},
46+
{
47+
"cell_type": "markdown",
48+
"metadata": {},
49+
"source": [
50+
"which are a sum of a quadratic bottom drag with dimensionless coefficient $c_D$ and a biharmonic diffusion with viscosity coefficient $\\nu$. $\\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)$."
51+
]
52+
},
53+
{
54+
"cell_type": "markdown",
55+
"metadata": {},
56+
"source": [
57+
"## The vector invariant formulation"
58+
]
59+
},
60+
{
61+
"cell_type": "markdown",
62+
"metadata": {},
63+
"source": [
64+
"The Bernoulli potential $p$ is introduced as"
65+
]
66+
},
67+
{
68+
"cell_type": "markdown",
69+
"metadata": {},
70+
"source": [
71+
"$$\n",
72+
"p = \\frac{1}{2}(u^2 + v^2) + gh\n",
73+
"$$"
74+
]
75+
},
76+
{
77+
"cell_type": "markdown",
78+
"metadata": {},
79+
"source": [
80+
"The relative vorticity $\\zeta = \\partial_xv + \\partial_yu$ lets us define the potential vorticity $q$ as"
81+
]
82+
},
83+
{
84+
"cell_type": "markdown",
85+
"metadata": {},
86+
"source": [
87+
"$$\n",
88+
"q = \\frac{f + \\zeta}{h}\n",
89+
"$$"
90+
]
91+
},
92+
{
93+
"cell_type": "markdown",
94+
"metadata": {},
95+
"source": [
96+
"such that we can rewrite the shallow water equations as"
97+
]
98+
},
99+
{
100+
"cell_type": "markdown",
101+
"metadata": {},
102+
"source": [
103+
"$$\n",
104+
"\\begin{align}\n",
105+
"\\partial_t u &= qhv -g\\partial_xp + D_x(u,v,\\eta) + F_x \\\\\n",
106+
"\\partial_t v &= -qhu -g\\partial_yp + D_y(u,v,\\eta) + F_y \\\\\n",
107+
"\\partial_t \\eta &= -\\partial_x(uh) -\\partial_y(vh)\n",
108+
"\\end{align}\n",
109+
"$$"
110+
]
111+
},
112+
{
113+
"cell_type": "markdown",
114+
"metadata": {},
115+
"source": [
116+
"## Runge-Kutta time discretisation"
117+
]
118+
},
119+
{
120+
"cell_type": "markdown",
121+
"metadata": {},
122+
"source": [
123+
"Let\n",
124+
"\n",
125+
"$$\n",
126+
"R(u,v,\\eta) = \\begin{pmatrix} \n",
127+
" qhv-g\\partial_xp+F_x \\\\ \n",
128+
" -qhu-g\\partial_yp+F_y \\\\\n",
129+
" -\\partial_x(uh) -\\partial_y(vh)\n",
130+
" \\end{pmatrix}\n",
131+
"$$\n",
132+
"\n",
133+
"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\n",
134+
"\n",
135+
"$$\n",
136+
"\\begin{align}\n",
137+
"\\mathbf{d}_1 &= R(\\mathbf{k}_n) \\\\\n",
138+
"\\mathbf{d}_2 &= R(\\mathbf{k}_n + \\frac{1}{2}\\Delta t \\mathbf{d}_1) \\\\\n",
139+
"\\mathbf{d}_3 &= R(\\mathbf{k}_n + \\frac{1}{2}\\Delta t \\mathbf{d}_2) \\\\\n",
140+
"\\mathbf{d}_4 &= R(\\mathbf{k}_n + \\Delta t \\mathbf{d}_3) \\\\\n",
141+
"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)\n",
142+
"\\end{align}\n",
143+
"$$"
144+
]
145+
},
146+
{
147+
"cell_type": "markdown",
148+
"metadata": {},
149+
"source": [
150+
"and the dissipative terms are then added semi-implictly."
151+
]
152+
},
153+
{
154+
"cell_type": "markdown",
155+
"metadata": {},
156+
"source": [
157+
"$$\n",
158+
"u_{n+1} = u_{n+1}^* + \\Delta t D_x(u_{n+1}^*,v_{n+1}^*,\\eta_{n+1}) \\\\\n",
159+
"v_{n+1} = v_{n+1}^* + \\Delta t D_y(u_{n+1}^*,v_{n+1}^*,\\eta_{n+1})\n",
160+
"$$"
161+
]
162+
},
163+
{
164+
"cell_type": "markdown",
165+
"metadata": {},
166+
"source": [
167+
"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$."
168+
]
169+
},
170+
{
171+
"cell_type": "markdown",
172+
"metadata": {},
173+
"source": [
174+
"## Strong stability preserving Runge-Kutta with semi-implicit continuity equation"
175+
]
176+
},
177+
{
178+
"cell_type": "markdown",
179+
"metadata": {},
180+
"source": [
181+
"We split the right-hand side into the momentum equations and the continuity equation"
182+
]
183+
},
184+
{
185+
"cell_type": "markdown",
186+
"metadata": {},
187+
"source": [
188+
"$$\n",
189+
"R_m(u,v,\\eta) = \\begin{pmatrix} \n",
190+
" qhv-g\\partial_xp+F_x \\\\ \n",
191+
" -qhu-g\\partial_yp+F_y \\\\\n",
192+
" \\end{pmatrix}, \\quad R_\\eta(u,v,\\eta) = -\\partial_x(uh) -\\partial_y(vh)\n",
193+
"$$"
194+
]
195+
},
196+
{
197+
"cell_type": "markdown",
198+
"metadata": {},
199+
"source": [
200+
"The 4-stage strong stability preserving Runge-Kutta scheme, with a semi-implicit treatment of the continuity equation then reads as"
201+
]
202+
},
203+
{
204+
"cell_type": "markdown",
205+
"metadata": {},
206+
"source": [
207+
"$$\n",
208+
"\\begin{align}\n",
209+
"\\mathbf{u}_1 &= \\mathbf{u}_n + \\frac{1}{2} \\Delta t R_m(u_n,v_n,\\eta_n), \\quad \\text{then}\n",
210+
"\\quad \\eta_1 = \\eta_n + \\frac{1}{2} \\Delta t R_\\eta(u_1,v_1,\\eta_n) \\\\\n",
211+
"\\mathbf{u}_2 &= \\mathbf{u}_1 + \\frac{1}{2} \\Delta t R_m(u_1,v_1,\\eta_1), \\quad \\text{then}\n",
212+
"\\quad \\eta_2 = \\eta_1 + \\frac{1}{2} \\Delta t R_\\eta(u_2,v_2,\\eta_1) \\\\\n",
213+
"\\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}\n",
214+
"\\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) \\\\\n",
215+
"\\mathbf{u}_{n+1} &= \\mathbf{u}_3 + \\frac{1}{2} \\Delta t R_m(u_3,v_3,\\eta_3), \\quad \\text{then}\n",
216+
"\\quad \\eta_{n+1} = \\eta_3 + \\frac{1}{2} \\Delta t R_\\eta(u_{n+1},v_{n+1},\\eta_3) \\\\\n",
217+
"\\end{align}\n",
218+
"$$"
219+
]
220+
},
221+
{
222+
"cell_type": "markdown",
223+
"metadata": {},
224+
"source": [
225+
"## Splitting the continuity equation"
226+
]
227+
},
228+
{
229+
"cell_type": "code",
230+
"execution_count": null,
231+
"metadata": {},
232+
"outputs": [],
233+
"source": [
234+
"From\n",
235+
"$$\n",
236+
"\\partial_t = -\\partial_x(uh) - \\partial_y(vh)\n",
237+
"$$\n",
238+
"\n"
239+
]
240+
}
241+
],
242+
"metadata": {
243+
"kernelspec": {
244+
"display_name": "Julia 1.5.2",
245+
"language": "julia",
246+
"name": "julia-1.5"
247+
},
248+
"language_info": {
249+
"file_extension": ".jl",
250+
"mimetype": "application/julia",
251+
"name": "julia",
252+
"version": "1.5.2"
253+
}
254+
},
255+
"nbformat": 4,
256+
"nbformat_minor": 4
257+
}

src/Constants.jl

Lines changed: 13 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,11 @@
11
struct Constants{T<:AbstractFloat,Tprog<:AbstractFloat}
22

3-
# RUNGE-KUTTA COEFFICIENTS 3rd/4th order including timestep Δt
3+
# RUNGE-KUTTA COEFFICIENTS 2nd/3rd/4th order including timestep Δt
44
RKaΔt::Array{Tprog,1}
55
RKbΔt::Array{Tprog,1}
6+
Δt_Δs::Tprog # Δt/(s-1) wher s the number of stages
7+
Δt_Δ::Tprog # Δt/Δ - timestep divided by grid spacing
8+
Δt_Δ_half::Tprog # 1/2 * Δt/Δ
69

710
# BOUNDARY CONDITIONS
811
one_minus_α::Tprog # tangential boundary condition for the ghost-point copy
@@ -39,6 +42,13 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
3942
RKbΔt = Tprog.([.5,.5,1.]*G.dtint/G.Δ)
4043
end
4144

45+
# Δt/(s-1) for SSPRK2
46+
Δt_Δs = Tprog(G.dtint/G.Δ/(P.RKs-1))
47+
48+
# time step and half the time step including the grid spacing as this is not included in the RHS
49+
Δt_Δ = Tprog(G.dtint/G.Δ)
50+
Δt_Δ_half = Tprog(G.dtint/G.Δ/2)
51+
4252
one_minus_α = Tprog(1-P.α) # for the ghost point copy/tangential boundary conditions
4353
g = T(P.g) # gravity - for Bernoulli potential
4454

@@ -65,7 +75,8 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
6575
ωFx = 2π*P.ωFx/24/365.25/3600
6676
ωFy = 2π*P.ωFy/24/365.25/3600
6777

68-
return Constants{T,Tprog}( RKaΔt,RKbΔt,one_minus_α,
78+
return Constants{T,Tprog}( RKaΔt,RKbΔt,Δt_Δs,Δt_Δ,Δt_Δ_half,
79+
one_minus_α,
6980
g,cD,rD,γ,cSmag,νB,rSST,
7081
jSST,SSTmin,ωFη,ωFx,ωFy)
7182
end

src/Continuity.jl

Lines changed: 19 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -69,16 +69,31 @@ end
6969

7070

7171
"""Transit function to call the specified continuity function."""
72-
function continuity!( η::AbstractMatrix,
72+
function continuity!( u::AbstractMatrix,
73+
v::AbstractMatrix,
74+
η::AbstractMatrix,
7375
Diag::DiagnosticVars,
7476
S::ModelSetup,
7577
t::Int)
78+
79+
@unpack U,V,dUdx,dVdy = Diag.VolumeFluxes
80+
@unpack nstep_advcor = S.grid
81+
@unpack time_scheme,surface_relax,surface_forcing = S.parameters
7682

77-
if S.parameters.surface_relax
83+
if nstep_advcor > 0 || # then UV wasn't calculated inside rhs!
84+
time_scheme != "RK" # then u,v changed before evaluating semi-implciti continuity
85+
UVfluxes!(u,v,η,Diag,S)
86+
end
87+
88+
# divergence of mass flux
89+
∂x!(dUdx,U)
90+
∂y!(dVdy,V)
91+
92+
if surface_relax
7893
continuity_surf_relax!(η,Diag,S,t)
79-
elseif S.parameters.surface_forcing
94+
elseif surface_forcing
8095
continuity_forcing!(Diag,S,t)
8196
else
8297
continuity_itself!(Diag,S,t)
8398
end
84-
end
99+
end

src/DefaultParameters.jl

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,9 @@
4747
wk::Real=10e3 # width [m] in y of Gaussian used for surface forcing
4848

4949
# TIME STEPPING OPTIONS
50-
RKo::Int=4 # Order of the RK time stepping scheme (3 or 4)
50+
time_scheme::String="RK" # Runge-Kutta ("RK") or strong-stability preserving RK ("SSPRK2","SSPRK3")
51+
RKo::Int=4 # Order of the RK time stepping scheme (2, 3 or 4)
52+
RKs::Int=2 # Number of stages for SSPRK2
5153
cfl::Real=1.0 # CFL number (1.0 recommended for RK4, 0.6 for RK3)
5254
Ndays::Real=10.0 # number of days to integrate for
5355
nstep_diff::Int=1 # diffusive part every nstep_diff time steps.
@@ -118,19 +120,21 @@
118120
@assert topo_width > 0.0 "topo_width has to be >0, $topo_width given."
119121
@assert t_relax > 0.0 "t_relax has to be >0, $t_relax given."
120122
@assert η_refw > 0.0 "η_refw has to be >0, $η_refw given."
121-
@assert RKo in [2,3,4] "RKo has to be 2,3 or 4, $RKo given."
123+
@assert time_scheme in ["RK","SSPRK2","SSPRK3"] "Time scheme $time_scheme unsupported."
124+
@assert RKo in [2,3,4] "RKo has to be 2,3 or 4; $RKo given."
125+
@assert RKs > 1 "RKs has to be >= 2; $RKs given."
122126
@assert Ndays > 0.0 "Ndays has to be >0, $Ndays given."
123127
@assert nstep_diff > 0 "nstep_diff has to be >0, $nstep_diff given."
124128
@assert nstep_advcor >= 0 "nstep_advcor has to be >=0, $nstep_advcor given."
125129
@assert bc in ["periodic","nonperiodic"] "boundary condition '$bc' unsupported."
126130
@assert α >= 0.0 && α <= 2.0 "Tangential boundary condition α has to be in (0,2), given."
127131
@assert adv_scheme in ["Sadourny","ArakawaHsu"] "Advection scheme '$adv_scheme' unsupported"
128-
@assert dynamics in ["linear","nonlinear"] "Dynamics '$dynamics' unsupoorted."
132+
@assert dynamics in ["linear","nonlinear"] "Dynamics '$dynamics' unsupported."
129133
@assert bottom_drag in ["quadratic","linear","none"] "Bottom drag '$bottom_drag' unsupported."
130134
@assert cD >= 0.0 "Bottom drag coefficient cD has to be >=0, $cD given."
131135
@assert τD >= 0.0 "Bottom drag coefficient τD has to be >=0, $τD given."
132-
@assert diffusion in ["Smagorinsky", "constant"] "Diffusion '$diffusion' unsupoorted."
133-
@assert νB > 0.0 "Diffusion scaling constant νB has to be > 0, $νB given."
136+
@assert diffusion in ["Smagorinsky", "constant"] "Diffusion '$diffusion' unsupported."
137+
@assert νB > 0.0 "Diffusion scaling constant νB has to be >0, $νB given."
134138
@assert cSmag > 0.0 "Smagorinsky coefficient cSmag has to be >0, $cSmag given."
135139
@assert injection_region in ["west","south"] "Injection region '$injection_region' unsupported."
136140
@assert Uadv > 0.0 "Advection velocity scale Uadv has to be >0, $Uadv given."

src/GhostPoints.jl

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -176,6 +176,24 @@ function ghost_points!( u::AbstractMatrix,
176176
end
177177
end
178178

179+
"""Decide on boundary condition P.bc which ghost point function to execute."""
180+
function ghost_points!( u::AbstractMatrix,
181+
v::AbstractMatrix,
182+
S::ModelSetup)
183+
184+
@unpack bc = S.parameters
185+
C = S.constants
186+
187+
if bc == "periodic"
188+
@unpack Tcomm = S.parameters
189+
ghost_points_u_periodic!(Tcomm,u,C)
190+
ghost_points_v_periodic!(Tcomm,v)
191+
else
192+
ghost_points_u_nonperiodic!(C,u)
193+
ghost_points_v_nonperiodic!(C,v)
194+
end
195+
end
196+
179197
"""Decide on boundary condition P.bc which ghost point function to execute."""
180198
function ghost_points_sst!(sst::AbstractMatrix,S::ModelSetup)
181199

0 commit comments

Comments
 (0)