Skip to content

Commit d593a2d

Browse files
Merge pull request #107 from FluidNumerics/demo/kelvin-wave
Demo/kelvin wave
2 parents a995422 + 1c75bdc commit d593a2d

File tree

10 files changed

+130
-55
lines changed

10 files changed

+130
-55
lines changed

docs/Models/linear-shallow-water-model.md

Lines changed: 21 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -36,15 +36,15 @@ where $g$ is acceleration due to gravity and $H$ is uniform resting fluid depth.
3636
$$
3737
\vec{q} =
3838
\begin{pmatrix}
39-
-fv \\
40-
fu \\
39+
-fv - C_d u\\
40+
fu - C_d v\\
4141
0
4242
\end{pmatrix}
4343
$$
4444

45-
where $f$ is the coriolis parameter.
45+
where $f$ is the coriolis parameter and $C_d$ is the linear drag coefficient.
4646

47-
To track stability of the Euler equation, the total entropy function is
47+
To track stability of the shallow water equations, the total entropy function is taken to be the total (kinetic plus potential) energy
4848

4949
$$
5050
e = \frac{1}{2} \int_V H u^2 + H v^2 + g \eta^2 \hspace{1mm} dV
@@ -57,7 +57,7 @@ The 2D Linear Shallow Water model is implemented as a type extension of the `DGM
5757
The `LinearShallowWater2D` class has a generic method (`SetCoriolis`) that can be used for defining the coriolis parameter at each location in the model domain. The `SetCoriolis` method can be used for either setting an $f$ or $beta$ plane.
5858

5959
#### Setting up an f-plane
60-
Assuming you've created interpolant ,mesh, geometry objects, and model objects you can define a constant value for the coriolis parameter using the following
60+
Assuming you've created interpolant, mesh, geometry objects, and model objects you can define a constant value for the coriolis parameter using the following
6161
```fortran
6262
type(LinearShallowWater2D) :: modelobj
6363
real(prec), parameter :: f0 = 10.0_prec*(-4)
@@ -68,7 +68,7 @@ real(prec), parameter :: f0 = 10.0_prec*(-4)
6868
```
6969

7070
#### Setting up a beta-plane
71-
Assuming you've created interpolant ,mesh, geometry objects, and model objects you can define the coriolis so that it varies with the `y` coordinate in the geometry using
71+
Assuming you've created interpolant, mesh, geometry objects, and model objects you can define the coriolis so that it varies with the `y` coordinate in the geometry using
7272
```fortran
7373
type(LinearShallowWater2D) :: modelobj
7474
real(prec), parameter :: f0 = 10.0_prec*(-4)
@@ -80,7 +80,7 @@ real(prec), parameter :: beta = 10.0_prec*(-11)
8080
```
8181

8282
#### Setting arbitrary spatially varying coriolis parameter
83-
Perhaps you find that f-plane and beta-plane scenarios are just too boring, or their not an appropriate model for what you're considering. In this case, you can easily set the `fCori%interior` attribute of the `LinearShallowWater2D` class directly
83+
Perhaps you find that f-plane and beta-plane scenarios are just too boring, or they're not an appropriate model for what you're considering. In this case, you can easily set the `fCori%interior` attribute of the `LinearShallowWater2D` class directly
8484

8585

8686
```fortran
@@ -128,6 +128,17 @@ real(prec), parameter :: beta = 10.0_prec*(-11)
128128
129129
```
130130

131+
### Setting the Drag coefficient
132+
Assuming you've created interpolant, mesh, geometry objects, and model objects you can define a constant value for the linear drag coefficient by setting the constant parameter `Cd`, e.g.
133+
134+
```fortran
135+
type(LinearShallowWater2D) :: modelobj
136+
real(prec), parameter :: fCd = 0.25
137+
...
138+
139+
modelobj % Cd = Cd ! Set the drag coefficient
140+
141+
```
131142
### Riemann Solver
132143
The `LinearShallowWater2D` class is defined using the advective form.
133144
The Riemann solver for the hyperbolic part of the shallow water equations is the local Lax-Friedrichs upwind Riemann solver
@@ -212,4 +223,6 @@ call mesh%StructuredMesh(nxPerTile=5,nyPerTile=5,&
212223

213224
For examples, see any of the following
214225

215-
* [`examples/LinearShallowWater2D.f90`](https://github.com/FluidNumerics/SELF/blob/main/examples/LinearShallowWater2D.f90) - Implements the 2D shallow water equations with no normal flow boundary conditions
226+
* [Gravity waves in closed square domain](../Tutorials/LinearShallowWater/ReflectingWave.md)
227+
* [Kelvin waves in a closed circular rotating domain (f-plane)](../Tutorials/LinearShallowWater/KelvinWaves.md)
228+
* [Planetary Rossby waves in an open square domain (beta-plane)](../Tutorials/LinearShallowWater/PlanetaryRossbyWave.md)
Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
# Kelvin waves
2+
This experiment is designed to demonstrate the preferred direction of phase propagation exhibited by Kelvin Waves. We use a circular domain in a rotating reference frame with a no-normal-flow boundary condition and an initial disturbance in the free-surface height placed on the domain boundary. The free surface height disturbance adjusts into geostrophic balance and in the process radiates gravity waves and a Kelvin wave. This demonstration uses a constant negative value for the coriolis parameter which results in a Kelvin wave that propagates in a clockwise direction, with the domain boundary to the left of the propagation direction.
3+
4+
## Configuration
5+
6+
### Equations
7+
8+
The equations solved are the linear shallow water equations, given by
9+
$$
10+
u_t - fv = -g \eta_x - C_d u
11+
$$
12+
$$
13+
v_t + fu = -g \eta_y - C_d v
14+
$$
15+
$$
16+
\eta_t + (Hu)_x + (Hv)_y = 0
17+
$$
18+
19+
where $\vec{u} = u \hat{x} + v \hat{y}$ is the barotropic velocity, $g$ is the acceleration of gravity, $C_d$ is the linear drag coefficient, $H$ is a uniform resting fluid depth, and $\eta$ is the deviation of the fluid free surface relative to the resting fluid.
20+
21+
An $f$-plane, in geophysical fluid dynamics, is an approximation that represents the vertical component of the coriolis force using a fixed coriolis frequency. The presence of a constant, non-zero coriolis frequency permits inertial oscillations and Kelvin waves.
22+
23+
For this simulation, we use the following parameters
24+
25+
* $g = 1 m s^{-2}$
26+
* $f_0 = 10 s^{-1}$
27+
* $\beta = 0 m^{-1} s^{-1}$
28+
* $H = 1 m$
29+
* $C_d = 0.25 s^{-1}$
30+
31+
### Domain Discretization
32+
In this problem, the domain is a circle of radius 1m. The model domain meshed using [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) and processed with [HOPr](https://github.com/hopr-framework/hopr). Within each element, the solution is approximated as a Lagrange interpolating polynomial of degree 7, using the Legendre-Gauss quadrature points as interpolating knots. To exchange momentum and mass fluxes between neighboring elements, we use a local upwind (Lax-Friedrich's) Riemann solver.
33+
34+
The physical parameters result in a gravity wave speed of $c= 1 m s^{-1}$. For this mesh, the elements are roughly isotropic with a length scale of about $\Delta x_e \approx 0.2 m$; with a polynomial degree of 7, the smallest resolvable length scale is roughly $\Delta x = \frac{0.2 m}{7^2} \approx 0.004 m$ .
35+
36+
For time integration, we use Williamson's low storage third order Runge Kutta and we choose a time step size of $\Delta t = 0.0025 s$ so that the CFL number associated with the gravity wave speed is $C_g = \frac{c \Delta t}{\Delta x} \approx 0.61 < 1$.
37+
38+
### Initial Condition and evolution of the model
39+
The initial condition is defined by setting the free surface height to a Gaussian, centered at $(x_c,y_c) = (1,0)$, with a half width of 10 cm and a height of 1 mm.
40+
$$
41+
\eta(t=0) = 0.001e^{ -( ( (x-1.0)^2 + y^2 )/(0.02) )}
42+
$$
43+
44+
This initial condition is initially out of balance, which causes an erruption of unbalanced flows, including gravity waves, inertia gravity waves, and Kelvin waves. The Kelvin waves are the result of the unbalanced flow up against the no-normal flow wall. Since the coriolis parameter is positive in this demonstration, the Kelvin waves propagate with the boundary (the "coast") on its right. For this circular domain, the Kelvin waves propagate in a counter-clockwise directtion.
45+
46+
<figure markdown="span">
47+
![Geostrophic adjustment releasing unbalanced flows](./img/kelvin-wave-initial-erruption.png){ width="500" }
48+
<figcaption> Free surface height (<code>eta</code>) shortly after the initial condition. Here, we see a train of gravity waves propagating into the domain and a single peak Kelvin wave traveling along the boundary in a counter-clockwise direction. The initial disturbance is now adjusted into geostrophic balance.
49+
</figcaption>
50+
</figure>
51+
52+
The release of energy into the unbalanced flows allows the initial disturbance to come into geostrophic balance. As a result, in the vicinity of the initial disturbance, we see a stationary high pressure signal that remains in geostrophic balance.
53+
54+
55+
56+
57+
## Running this example
58+
59+
To run this example, simply execute
60+
61+
```shell
62+
# Set WORKSPACE to the path to the SELF source code
63+
export WORKSPACE=/path/to/SELF
64+
${SELF_ROOT}/examples/linear_shallow_water2d_kelvinwaves
65+
```
66+
67+
This will run the simulation from $t=0$ to $t=1.0$ and write model output at intervals of $Δ t_{io} = 0.05$. Model output can be visualized using `pyself` in python
68+
From the SELF source code directory
69+
70+
```shell
71+
# Assuming you are in the SELF source code directory
72+
pip install . --upgrade
73+
```
74+
You can use the `examples/shallow_water_plot.py` script to plot the model output and generate movie of the free surface height.

docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -35,8 +35,6 @@ $$
3535
\eta(t=0) = 0.01e^{ -( (x^2 + y^2 )/(2.0*10.0^{10}) )}
3636
$$
3737

38-
![Rossby Wave Initial Condition](./rossbywave_initialcondition.png){ align=center }
39-
4038
The initial velocity field is calculated by using the pressure gradient force and using geostrophic balance; in SELF, this is handled by the `LinearShallowWater % DiagnoseGeostrophicVelocity` type bound procedure after setting the initial free surface height.
4139

4240
### Boundary Conditions

docs/Tutorials/LinearShallowWater/LinearShallowWater.md renamed to docs/Tutorials/LinearShallowWater/ReflectingWave.md

Lines changed: 15 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
# Linear Shallow Water No Normal Flow Tutorial
1+
# Reflecting Wave
22

33
This tutorial will walk you through using an example program that uses the `LinearShallowWater2D` class to run a simulation with the linear shallow water equations in 2-D. This example is configured using the built in structured mesh generator with no normal flow boundary conditions on all domain boundaries.
44

@@ -82,15 +82,18 @@ $$
8282

8383
The model is integrated forward in time using $3^{rd}$ order Runge-Kutta with a time step of $\Delta t = 0.5 s$.
8484

85-
<p align="center">
86-
<img height="440px" src="img/shallow-water.0000.png" />
87-
Free surface height (<code>eta</code>) at time <code>t=0</code>.
88-
</p>
85+
<figure markdown="span">
86+
![Initial condition](./img/shallow-water.0000.png){ width="500" }
87+
<figcaption> Free surface height (<code>eta</code>) at time <code>t=0</code>.
88+
</figcaption>
89+
</figure>
90+
91+
<figure markdown="span">
92+
![Gravity wave reflection](./img/shallow-water.0019.png){ width="500" }
93+
<figcaption> Free surface height (<code>eta</code>) at time <code>t=1</code>.
94+
</figcaption>
95+
</figure>
8996

90-
<p align="center">
91-
<img height="440px" src="img/shallow-water.0019.png" />
92-
Free surface height (<code>eta</code>) at time <code>t=1</code>.
93-
</p>
9497

9598
## How we implement this
9699
You can find the example file for this demo in the `examples/linear_shallow_water2d_nonormalflow.f90` file. This file uses the `LinearShallowWater2D` module from `src/SELF_LinearShallowWater2D_t.f90`.
@@ -243,14 +246,9 @@ Running this program should output twenty `shallow-water.00XX.tec` in the build
243246

244247
## Running this example
245248

246-
<p align="center">
247-
<div align="center">
248-
<img height="360px" src="img/shallow-water.gif" />
249-
</div>
250-
<div align="center">
251-
Free surface height (<code>eta</code>) for the full duration (1 second) of the problem.
252-
</div>
253-
</p>
249+
<figure markdown="span">
250+
![Geostrophic adjustment releasing unbalanced flows](./img/shallow-water.gif){ width="500" }
251+
</figure>
254252

255253
!!! note
256254
To run this example, you must first [install SELF](../../GettingStarted/install.md) . We assume that SELF is installed in path referenced by the `SELF_ROOT` environment variable.
261 KB
Loading

examples/linear_shallow_water2d_kelvinwaves.f90

Lines changed: 8 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -33,9 +33,10 @@ program linear_shallow_water2d_kelvinwaves
3333
character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' ! Which integrator method
3434
integer,parameter :: controlDegree = 7 ! Degree of control polynomial
3535
integer,parameter :: targetDegree = 16 ! Degree of target polynomial
36-
real(prec),parameter :: dt = 0.001_prec ! Time-step size
36+
real(prec),parameter :: dt = 0.0025_prec ! Time-step size
3737
real(prec),parameter :: endtime = 1.0_prec !30.0_prec ! (s);
38-
real(prec),parameter :: f0 = -10.0_prec ! reference coriolis parameter (1/s)
38+
real(prec),parameter :: f0 = 10.0_prec ! reference coriolis parameter (1/s)
39+
real(prec),parameter :: Cd = 0.25_prec ! Linear drag coefficient (1/s)
3940
real(prec),parameter :: iointerval = 0.05 ! Write files 20 times per characteristic time scale
4041
real(prec) :: r
4142
real(prec) :: e0,ef ! Initial and final entropy
@@ -71,25 +72,13 @@ program linear_shallow_water2d_kelvinwaves
7172
! Set the resting surface height and gravity
7273
modelobj%H = H
7374
modelobj%g = g
75+
modelobj%Cd = Cd
7476

7577
! ! Set the initial conditions
76-
!call modelobj%solution%SetEquation(3,'f = 0.001*exp( -( x^2 + y^2 )/0.02 ) ')
77-
!call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
78-
79-
do iel = 1,modelobj%mesh%nElem
80-
do j = 1,modelobj%solution%N+1
81-
do i = 1,modelobj%solution%N+1
82-
call random_number(r)
83-
! modelobj%solution%interior(i,j,iel,3) = modelobj%solution%interior(i,j,iel,3) + 0.0001_prec*(r-0.5)
84-
modelobj%solution%interior(i,j,iel,3) = 0.0001_prec*(r-0.5)
85-
86-
enddo
87-
enddo
88-
enddo
89-
call modelobj%solution%UpdateDevice()
78+
call modelobj%solution%SetEquation(3,'f = 0.001*exp( -( (x-1.0)^2 + y^2 )/0.02 ) ')
79+
call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
9080

9181
call modelobj%SetCoriolis(f0)
92-
call modelobj%DiagnoseGeostrophicVelocity()
9382

9483
call modelobj%WriteModel()
9584
call modelobj%IncrementIOCounter()
@@ -108,9 +97,8 @@ program linear_shallow_water2d_kelvinwaves
10897

10998
ef = modelobj%entropy
11099

111-
print*,e0,ef
112-
if(abs(ef-e0) > epsilon(e0)) then
113-
print*,"Warning: Final entropy greater than initial entropy! ",e0,ef
100+
if(ef > e0) then
101+
print*,"Warning: Final entropy greater than initial entropy! ",ef,e0
114102
endif
115103

116104
! Clean up

mkdocs.yml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ nav:
3131
- Spherical sound wave: Tutorials/LinearEuler3D/SphericalSoundWave.md
3232
- Linear Shallow Water (2D):
3333
- Reflecting wave: Tutorials/LinearShallowWater/LinearShallowWater.md
34+
- Kelvin waves: Tutorials/LinearShallowWater/KelvinWaves.md
3435
- Planetary Rossby wave: Tutorials/LinearShallowWater/PlanetaryRossbyWave.md
3536
- Models:
3637
- Viscous Burgers Equation: Models/burgers-equation-model.md

src/SELF_LinearShallowWater2D_t.f90

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ module self_LinearShallowWater2D_t
3434
type,extends(dgmodel2d) :: LinearShallowWater2D_t
3535
real(prec) :: H = 0.0_prec ! uniform resting depth
3636
real(prec) :: g = 0.0_prec ! acceleration due to gravity
37+
real(prec) :: Cd = 0.0_prec ! Linear drag coefficient (1/s)
3738
type(MappedScalar2D) :: fCori ! The coriolis parameter
3839

3940
contains
@@ -246,8 +247,8 @@ subroutine sourcemethod_LinearShallowWater2D_t(this)
246247

247248
s = this%solution%interior(i,j,iel,1:this%nvar)
248249

249-
this%source%interior(i,j,iel,1) = this%fCori%interior(i,j,iel,1)*s(2) ! du/dt = f*v
250-
this%source%interior(i,j,iel,2) = -this%fCori%interior(i,j,iel,1)*s(1) ! dv/dt = -f*u
250+
this%source%interior(i,j,iel,1) = this%fCori%interior(i,j,iel,1)*s(2)-this%Cd*s(1) ! du/dt = f*v - Cd*u
251+
this%source%interior(i,j,iel,2) = -this%fCori%interior(i,j,iel,1)*s(1)-this%Cd*s(2) ! dv/dt = -f*u - Cd*v
251252

252253
enddo
253254

src/gpu/SELF_LinearShallowWater2D.cpp

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -144,25 +144,25 @@ extern "C"
144144
}
145145
}
146146

147-
__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, int ndof){
147+
__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, real Cd, int ndof){
148148
uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x;
149149

150150
if( idof < ndof ){
151151
real u = solution[idof];
152152
real v = solution[idof + ndof];
153153

154-
source[idof] = fCori[idof]*v; // du/dt = fv
155-
source[idof+ndof] = -fCori[idof]*u; // dv/dt = -fu
154+
source[idof] = fCori[idof]*v - Cd*u; // du/dt = fv - Cd*u
155+
source[idof+ndof] = -fCori[idof]*u - Cd*v; // dv/dt = -fu - Cd*v
156156

157157
}
158158

159159
}
160160
extern "C"
161161
{
162-
void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, int N, int nel, int nvar){
162+
void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, real Cd, int N, int nel, int nvar){
163163
int ndof = (N+1)*(N+1)*nel;
164164
int threads_per_block = 256;
165165
int nblocks_x = ndof/threads_per_block +1;
166-
sourcemethod_LinearShallowWater2D_gpukernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(solution,source,fCori,ndof);
166+
sourcemethod_LinearShallowWater2D_gpukernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(solution,source,fCori,Cd,ndof);
167167
}
168168
}

src/gpu/SELF_LinearShallowWater2D.f90

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,11 +71,12 @@ subroutine fluxmethod_LinearShallowWater2D_gpu(solution,flux,g,H,N,nel,nvar) &
7171
endinterface
7272

7373
interface
74-
subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,N,nel,nvar) &
74+
subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,Cd,N,nel,nvar) &
7575
bind(c,name="sourcemethod_LinearShallowWater2D_gpu")
7676
use iso_c_binding
7777
use SELF_Constants
7878
type(c_ptr),value :: solution,source,fCori
79+
real(c_prec),value :: Cd
7980
integer(c_int),value :: N,nel,nvar
8081
endsubroutine sourcemethod_LinearShallowWater2D_gpu
8182
endinterface
@@ -165,6 +166,7 @@ subroutine sourcemethod_LinearShallowWater2D(this)
165166
call sourcemethod_LinearShallowWater2D_gpu(this%solution%interior_gpu, &
166167
this%source%interior_gpu, &
167168
this%fCori%interior_gpu, &
169+
this%Cd, &
168170
this%solution%interp%N, &
169171
this%solution%nelem, &
170172
this%solution%nvar)

0 commit comments

Comments
 (0)