Merge pull request #107 from FluidNumerics/demo/kelvin-wave

Demo/kelvin wave

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

@@ -36,15 +36,15 @@ where $g$ is acceleration due to gravity and $H$ is uniform resting fluid depth.
 $$
     \vec{q} = 
         \begin{pmatrix}
-        -fv \\ 
-        fu \\ 
+        -fv - C_d u\\ 
+        fu - C_d v\\ 
         0
     \end{pmatrix}
 $$
 
-where $f$ is the coriolis parameter.
+where $f$ is the coriolis parameter and $C_d$ is the linear drag coefficient.
 
-To track stability of the Euler equation, the total entropy function is
+To track stability of the shallow water equations, the total entropy function is taken to be the total (kinetic plus potential) energy
 
 $$
     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
 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.
 
 #### Setting up an f-plane
-Assuming you've created interpolant ,mesh, geometry objects, and model objects you can define a constant value for the coriolis parameter using the following
+Assuming you've created interpolant, mesh, geometry objects, and model objects you can define a constant value for the coriolis parameter using the following
 ```fortran
 type(LinearShallowWater2D) :: modelobj
 real(prec), parameter :: f0 = 10.0_prec*(-4)
@@ -68,7 +68,7 @@ real(prec), parameter :: f0 = 10.0_prec*(-4)
 ```
 
 #### Setting up a beta-plane
-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
+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
 ```fortran
 type(LinearShallowWater2D) :: modelobj
 real(prec), parameter :: f0 = 10.0_prec*(-4)
@@ -80,7 +80,7 @@ real(prec), parameter :: beta = 10.0_prec*(-11)
 ```
 
 #### Setting arbitrary spatially varying coriolis parameter
-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
+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
 
 
 ```fortran
@@ -128,6 +128,17 @@ real(prec), parameter :: beta = 10.0_prec*(-11)
 
 ```
 
+### Setting the Drag coefficient
+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. 
+
+```fortran
+type(LinearShallowWater2D) :: modelobj
+real(prec), parameter :: fCd = 0.25
+...
+
+  modelobj % Cd = Cd ! Set the drag coefficient
+
+```
 ### Riemann Solver
 The `LinearShallowWater2D` class is defined using the advective form.
 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,&
 
 For examples, see any of the following
 
-* [`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
+* [Gravity waves in closed square domain](../Tutorials/LinearShallowWater/ReflectingWave.md)
+* [Kelvin waves in a closed circular rotating domain (f-plane)](../Tutorials/LinearShallowWater/KelvinWaves.md)
+* [Planetary Rossby waves in an open square domain (beta-plane)](../Tutorials/LinearShallowWater/PlanetaryRossbyWave.md)

docs/Tutorials/LinearShallowWater/KelvinWaves.md

@@ -0,0 +1,74 @@
+# Kelvin waves
+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.
+
+## Configuration
+
+### Equations
+
+The equations solved are the linear shallow water equations, given by
+$$
+    u_t - fv = -g \eta_x - C_d u
+$$
+$$
+    v_t + fu = -g \eta_y - C_d v
+$$
+$$
+    \eta_t + (Hu)_x + (Hv)_y = 0
+$$
+
+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.
+
+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.
+
+For this simulation, we use the following parameters
+
+* $g = 1 m s^{-2}$
+* $f_0 = 10 s^{-1}$
+* $\beta = 0 m^{-1} s^{-1}$
+* $H = 1 m$
+* $C_d = 0.25 s^{-1}$
+
+### Domain Discretization
+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.
+
+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$ . 
+
+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$.
+
+### Initial Condition and evolution of the model
+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.
+$$
+    \eta(t=0) = 0.001e^{ -( ( (x-1.0)^2 + y^2 )/(0.02) )}
+$$
+
+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. 
+
+<figure markdown="span">
+  ![Geostrophic adjustment releasing unbalanced flows](./img/kelvin-wave-initial-erruption.png){ width="500" }
+  <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.
+  </figcaption>
+</figure>
+
+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.
+
+
+
+
+## Running this example
+
+To run this example, simply execute
+
+```shell
+# Set WORKSPACE to the path to the SELF source code
+export WORKSPACE=/path/to/SELF
+${SELF_ROOT}/examples/linear_shallow_water2d_kelvinwaves
+```
+
+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
+From the SELF source code directory
+
+```shell
+# Assuming you are in the SELF source code directory
+pip install . --upgrade
+```
+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

@@ -35,8 +35,6 @@ $$
     \eta(t=0) = 0.01e^{ -( (x^2 + y^2 )/(2.0*10.0^{10}) )}
 $$
 
-![Rossby Wave Initial Condition](./rossbywave_initialcondition.png){ align=center }
-
 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.
 
 ### Boundary Conditions

docs/Tutorials/LinearShallowWater/LinearShallowWater.md → docs/Tutorials/LinearShallowWater/ReflectingWave.md

@@ -1,4 +1,4 @@
-# Linear Shallow Water No Normal Flow Tutorial
+# Reflecting Wave
 
 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.
 
@@ -82,15 +82,18 @@ $$
 
 The model is integrated forward in time using $3^{rd}$ order Runge-Kutta with a time step of $\Delta t = 0.5 s$. 
 
-<p align="center">
-  <img height="440px" src="img/shallow-water.0000.png" />
-  Free surface height (<code>eta</code>) at time <code>t=0</code>.
-</p>
+<figure markdown="span">
+  ![Initial condition](./img/shallow-water.0000.png){ width="500" }
+  <figcaption>  Free surface height (<code>eta</code>) at time <code>t=0</code>.
+  </figcaption>
+</figure>
+
+<figure markdown="span">
+  ![Gravity wave reflection](./img/shallow-water.0019.png){ width="500" }
+  <figcaption>  Free surface height (<code>eta</code>) at time <code>t=1</code>.
+  </figcaption>
+</figure>
 
-<p align="center">
-  <img height="440px" src="img/shallow-water.0019.png" />
-  Free surface height (<code>eta</code>) at time <code>t=1</code>.
-</p>
 
 ## How we implement this
 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
 
 ## Running this example
 
-<p align="center">
-  <div align="center">
-    <img height="360px" src="img/shallow-water.gif" />
-  </div>
-  <div align="center">
-    Free surface height (<code>eta</code>) for the full duration (1 second) of the problem.
-  </div>
-</p>
+<figure markdown="span">
+ ![Geostrophic adjustment releasing unbalanced flows](./img/shallow-water.gif){ width="500" }
+</figure>
 
 !!! note
     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.

examples/linear_shallow_water2d_kelvinwaves.f90

@@ -33,9 +33,10 @@ program linear_shallow_water2d_kelvinwaves
   character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' ! Which integrator method
   integer,parameter :: controlDegree = 7 ! Degree of control polynomial
   integer,parameter :: targetDegree = 16 ! Degree of target polynomial
-  real(prec),parameter :: dt = 0.001_prec ! Time-step size
+  real(prec),parameter :: dt = 0.0025_prec ! Time-step size
   real(prec),parameter :: endtime = 1.0_prec !30.0_prec ! (s);
-  real(prec),parameter :: f0 = -10.0_prec ! reference coriolis parameter (1/s)
+  real(prec),parameter :: f0 = 10.0_prec ! reference coriolis parameter (1/s)
+  real(prec),parameter :: Cd = 0.25_prec ! Linear drag coefficient (1/s)
   real(prec),parameter :: iointerval = 0.05 ! Write files 20 times per characteristic time scale
   real(prec) :: r
   real(prec) :: e0,ef ! Initial and final entropy
@@ -71,25 +72,13 @@ program linear_shallow_water2d_kelvinwaves
   ! Set the resting surface height and gravity
   modelobj%H = H
   modelobj%g = g
+  modelobj%Cd = Cd
 
   ! ! Set the initial conditions
-  !call modelobj%solution%SetEquation(3,'f = 0.001*exp( -( x^2 + y^2 )/0.02 ) ')
-  !call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
-
-  do iel = 1,modelobj%mesh%nElem
-    do j = 1,modelobj%solution%N+1
-      do i = 1,modelobj%solution%N+1
-        call random_number(r)
-        !  modelobj%solution%interior(i,j,iel,3) = modelobj%solution%interior(i,j,iel,3) + 0.0001_prec*(r-0.5)
-        modelobj%solution%interior(i,j,iel,3) = 0.0001_prec*(r-0.5)
-
-      enddo
-    enddo
-  enddo
-  call modelobj%solution%UpdateDevice()
+  call modelobj%solution%SetEquation(3,'f = 0.001*exp( -( (x-1.0)^2 + y^2 )/0.02 ) ')
+  call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
 
   call modelobj%SetCoriolis(f0)
-  call modelobj%DiagnoseGeostrophicVelocity()
 
   call modelobj%WriteModel()
   call modelobj%IncrementIOCounter()
@@ -108,9 +97,8 @@ program linear_shallow_water2d_kelvinwaves
 
   ef = modelobj%entropy
 
-  print*,e0,ef
-  if(abs(ef-e0) > epsilon(e0)) then
-    print*,"Warning: Final entropy greater than initial entropy! ",e0,ef
+  if(ef > e0) then
+    print*,"Warning: Final entropy greater than initial entropy! ",ef,e0
   endif
 
   ! Clean up

mkdocs.yml

@@ -31,6 +31,7 @@ nav:
       - Spherical sound wave: Tutorials/LinearEuler3D/SphericalSoundWave.md
     - Linear Shallow Water (2D):
       - Reflecting wave: Tutorials/LinearShallowWater/LinearShallowWater.md
+      - Kelvin waves: Tutorials/LinearShallowWater/KelvinWaves.md
       - Planetary Rossby wave: Tutorials/LinearShallowWater/PlanetaryRossbyWave.md
   - Models:
     - Viscous Burgers Equation: Models/burgers-equation-model.md

src/SELF_LinearShallowWater2D_t.f90

@@ -34,6 +34,7 @@ module self_LinearShallowWater2D_t
   type,extends(dgmodel2d) :: LinearShallowWater2D_t
     real(prec) :: H = 0.0_prec ! uniform resting depth
     real(prec) :: g = 0.0_prec ! acceleration due to gravity
+    real(prec) :: Cd = 0.0_prec ! Linear drag coefficient (1/s)
     type(MappedScalar2D) :: fCori ! The coriolis parameter
 
   contains
@@ -246,8 +247,8 @@ subroutine sourcemethod_LinearShallowWater2D_t(this)
 
       s = this%solution%interior(i,j,iel,1:this%nvar)
 
-      this%source%interior(i,j,iel,1) = this%fCori%interior(i,j,iel,1)*s(2) ! du/dt = f*v
-      this%source%interior(i,j,iel,2) = -this%fCori%interior(i,j,iel,1)*s(1) ! dv/dt = -f*u
+      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
+      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
 
     enddo
 

src/gpu/SELF_LinearShallowWater2D.cpp

@@ -144,25 +144,25 @@ extern "C"
   }
 }
 
-__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, int ndof){
+__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, real Cd, int ndof){
   uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x;
 
   if( idof < ndof ){
     real u = solution[idof];
     real v = solution[idof + ndof];
 
-    source[idof] = fCori[idof]*v; // du/dt = fv
-    source[idof+ndof] = -fCori[idof]*u;   // dv/dt  = -fu
+    source[idof] = fCori[idof]*v - Cd*u; // du/dt = fv - Cd*u
+    source[idof+ndof] = -fCori[idof]*u - Cd*v;   // dv/dt  = -fu - Cd*v
 
   }
 
 }
 extern "C"
 {
-  void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, int N, int nel, int nvar){
+  void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, real Cd, int N, int nel, int nvar){
     int ndof = (N+1)*(N+1)*nel;
     int threads_per_block = 256;
     int nblocks_x = ndof/threads_per_block +1;
-    sourcemethod_LinearShallowWater2D_gpukernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(solution,source,fCori,ndof);
+    sourcemethod_LinearShallowWater2D_gpukernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(solution,source,fCori,Cd,ndof);
   }
 }

src/gpu/SELF_LinearShallowWater2D.f90

@@ -71,11 +71,12 @@ subroutine fluxmethod_LinearShallowWater2D_gpu(solution,flux,g,H,N,nel,nvar) &
   endinterface
 
   interface
-    subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,N,nel,nvar) &
+    subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,Cd,N,nel,nvar) &
       bind(c,name="sourcemethod_LinearShallowWater2D_gpu")
       use iso_c_binding
       use SELF_Constants
       type(c_ptr),value :: solution,source,fCori
+      real(c_prec),value :: Cd
       integer(c_int),value :: N,nel,nvar
     endsubroutine sourcemethod_LinearShallowWater2D_gpu
   endinterface
@@ -165,6 +166,7 @@ subroutine sourcemethod_LinearShallowWater2D(this)
     call sourcemethod_LinearShallowWater2D_gpu(this%solution%interior_gpu, &
                                                this%source%interior_gpu, &
                                                this%fCori%interior_gpu, &
+                                               this%Cd, &
                                                this%solution%interp%N, &
                                                this%solution%nelem, &
                                                this%solution%nvar)