Skip to content

Commit 70d38d1

Browse files
milanklmilankl
committed
equations and stability analysis
1 parent 68b7f9c commit 70d38d1

File tree

2 files changed

+257
-0
lines changed

2 files changed

+257
-0
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+
}

‎wip/ssprk2.png

68.9 KB
Loading

0 commit comments

Comments
 (0)