Add images and description of drag coefficient to model

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
@@ -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/LinearShallowWater.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

@@ -43,14 +43,32 @@ $$
 
 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. 
 
-<p align="center">
-  <img height="440px" src="./img/kelvin-wave-initial-erruption.png" />
-  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.
-</p>
+<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/LinearShallowWater.md

@@ -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.

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

examples/linear_shallow_water2d_kelvinwaves.f90

@@ -34,9 +34,9 @@ program linear_shallow_water2d_kelvinwaves
   integer,parameter :: controlDegree = 7 ! Degree of control polynomial
   integer,parameter :: targetDegree = 16 ! Degree of target polynomial
   real(prec),parameter :: dt = 0.0025_prec ! Time-step size
-  real(prec),parameter :: endtime = 30.0_prec ! (s);
+  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 :: Cd = 0.5_prec ! Linear drag coefficient (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
@@ -78,20 +78,7 @@ program linear_shallow_water2d_kelvinwaves
   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)
 
-  ! 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)*(1.0_prec+1.0_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%SetCoriolis(f0)
-  ! call modelobj%DiagnoseGeostrophicVelocity()
 
   call modelobj%WriteModel()
   call modelobj%IncrementIOCounter()

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