Feature/rossby wave demo

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

@@ -31,12 +31,19 @@ $$
     \end{pmatrix}
 $$
 
-where $g$ is acceleration due to gravity and $H$ is uniform resting fluid depth. The source term is set to zero.
+where $g$ is acceleration due to gravity and $H$ is uniform resting fluid depth. The source term includes a coriolis force
 
 $$
-    \vec{q} = \vec{0}
+    \vec{q} = 
+        \begin{pmatrix}
+        -fv \\ 
+        fu \\ 
+        0
+    \end{pmatrix}
 $$
 
+where $f$ is the coriolis parameter.
+
 To track stability of the Euler equation, the total entropy function is
 
 $$
@@ -46,6 +53,81 @@ $$
 ## Implementation
 The 2D Linear Shallow Water model is implemented as a type extension of the `DGModel2d` class. The `LinearShallowWater2D_t` class adds parameters for acceleration due to gravity and the uniform resting fluid depth. It also overrides `SetMetaData`, `entropy_func`, `flux2d`, and `riemannflux2d` type-bound procedures.
 
+### Defining the coriolis parameter and geostrophic velocities
+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
+```fortran
+type(LinearShallowWater2D) :: modelobj
+real(prec), parameter :: f0 = 10.0_prec*(-4)
+...
+
+  call modelobj%SetCoriolis(f0)
+
+```
+
+#### 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
+```fortran
+type(LinearShallowWater2D) :: modelobj
+real(prec), parameter :: f0 = 10.0_prec*(-4)
+real(prec), parameter :: beta = 10.0_prec*(-11) 
+...
+
+  call modelobj%SetCoriolis(f0,beta)
+
+```
+
+#### 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
+
+
+```fortran
+type(LinearShallowWater2D) :: modelobj
+integer :: iel
+integer :: i
+integer :: j
+
+do concurrent(i=1:modelobj%solution%N+1,j=1:modelobj%solution%N+1,iel=1:modelobj%mesh%nElem)
+    x = modelobj%geometry%x%interior(i,j,1,iel,1) ! Get the x-coordinate
+    y = modelobj%geometry%x%interior(i,j,1,iel,2) ! Get the y-coordinate
+    this%fCori%interior(i,j,1,iel) =  ! Define the coriolis parameter here as a function of x and y
+enddo
+call this%fCori%UpdateDevice()
+```
+
+#### Defining Geostrophic velocities
+With the `fCori` attribute defined, you can define geostrophic velocities from an initial condition for the free-surface height.
+
+!!! note
+    Setting geostrophic velocities is only valid when $f \neq 0$ everywhere in the domain.
+
+To define geostrophic velocities, you can simply use the `DiagnoseGeostrophicVelocity` procedure. This will 
+
+* Reset the velocity field to zero,
+* Calculate the free-surface height gradients using the `CalculateTendency` method
+* Compute the `u` and `v` variables using geostrophic balance
+
+As an example,
+
+```fortran
+type(LinearShallowWater2D) :: modelobj
+real(prec), parameter :: f0 = 10.0_prec*(-4)
+real(prec), parameter :: beta = 10.0_prec*(-11) 
+...
+
+  call modelobj%SetCoriolis(f0,beta)  
+
+  ! Set the free-surface height using an equation string
+  call modelobj%solution%SetEquation(3,'f = 0.01*exp( -( (x-500000.0)^2 + (y-500000.0)^2 )/(2.0*(10.0^10)) )')
+  call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
+
+  ! Calculate u and v from the free surface height using geostrophy
+  call modelobj%DiagnoseGeostrophicVelocity()
+
+```
+
 ### 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

docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md

@@ -0,0 +1,53 @@
+# Planetary Rossby Wave
+This experiment is designed to show how a geostrophic monopole evolves on a $\beta$-plane. 
+
+## Configuration
+
+### Equations
+
+The equations solved are the linear shallow water equations, given by
+$$
+    u_t - fv = -g \eta_x
+$$
+$$
+    v_t + fu = -g \eta_y
+$$
+$$
+    \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, $H$ is a uniform resting fluid depth, and $\eta$ is the deviation of the fluid free surface relative to the resting fluid. In this model, the $x$ direction is similar to longitude and $y$ is similar to latitude. 
+
+A $\beta$-plane, in geophysical fluid dynamics, is an approximation that accounts for first order variability in the (vertical component of the) coriolis parameter with latitude, 
+
+$$
+    f = f_0 + \beta y
+$$
+
+The background variation in the planetary vorticity supports Rossby waves, which propagate "westward" with higher potential vorticity to the right of phase propagation.
+
+### Domain Discretization
+In this problem, the domain is a square with $(x,y) \in [-500km, 500km]^2$. The model domain is divided into $10\times 10$ elements of uniform size. 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.
+
+### Initial Condition
+The initial condition is defined by setting the free surface height to a Gaussian, centered at the origin, with a half width of 10 km and a height of 1 cm.
+$$
+    \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
+Radiation boundary conditions are applied by setting the external state to a motionless fluid with no free surface height variation  ( $u=v=0, \eta = 0$). The model is integrated forward in time using Williamson's $3^{rd}$ order low storage Runge-Kutta, with a time step of $\Delta t = 0.5 s$. 
+
+### Physical Parameters
+The remaining parameters for the problem are as follows
+
+* $g = 10 m s^{-2}$
+* $f_0 = 10^{-4} s^{-1}$
+* $\beta = 10^{-11} m^{-1} s^{-1}$
+* $H = 1000 m$
+
+

examples/CMakeLists.txt

@@ -66,5 +66,7 @@ add_fortran_tests (
   "linear_euler2d_planewave_reflection.f90"
   "linear_euler2d_planewave_propagation.f90"
   "linear_shallow_water2d_nonormalflow.f90"
+  "linear_shallow_water2d_planetaryrossby_wave.f90"
+  "linear_shallow_water2d_kelvinwaves.f90"
   "linear_euler3d_spherical_soundwave_radiation.f90"
     )

examples/linear_shallow_water2d_kelvinwaves.f90

@@ -0,0 +1,122 @@
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+!
+! Maintainers : support@fluidnumerics.com
+! Official Repository : https://github.com/FluidNumerics/self/
+!
+! Copyright © 2024 Fluid Numerics LLC
+!
+! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+!
+! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+!
+! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in
+!    the documentation and/or other materials provided with the distribution.
+!
+! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from
+!    this software without specific prior written permission.
+!
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+!
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+
+program linear_shallow_water2d_kelvinwaves
+  use self_data
+  use self_LinearShallowWater2D
+  use self_mesh_2d
+
+  implicit none
+  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 :: endtime = 1.0_prec !30.0_prec ! (s);
+  real(prec),parameter :: f0 = -10.0_prec ! reference coriolis parameter (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
+  type(LinearShallowWater2D) :: modelobj ! Shallow water model
+  type(Lagrange),target :: interp ! Interpolant
+  type(Mesh2D),target :: mesh ! Mesh class
+  type(SEMQuad),target :: geometry ! Geometry class
+  integer :: i,j,iel
+  real(prec),parameter :: g = 1.0_prec ! Acceleration due to gravity
+  real(prec),parameter :: H = 1.0_prec ! Uniform resting depth
+  character(LEN=255) :: WORKSPACE
+
+  ! Create a uniform block mesh
+  call get_environment_variable("WORKSPACE",WORKSPACE)
+  call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Circle/Circle_mesh.h5")
+  call mesh%ResetBoundaryConditionType(SELF_BC_NONORMALFLOW)
+
+  ! Create an interpolant
+  call interp%Init(N=controlDegree, &
+                   controlNodeType=GAUSS, &
+                   M=targetDegree, &
+                   targetNodeType=UNIFORM)
+
+  ! Generate geometry (metric terms) from the mesh elements
+  call geometry%Init(interp,mesh%nElem)
+  call geometry%GenerateFromMesh(mesh)
+
+  ! Initialize the model
+  call modelobj%Init(mesh,geometry)
+  modelobj%prescribed_bcs_enabled = .false. ! Disables prescribed boundary condition block for gpu accelerated implementations
+  modelobj%tecplot_enabled = .false. ! Disables tecplot output
+
+  ! Set the resting surface height and gravity
+  modelobj%H = H
+  modelobj%g = g
+
+  ! ! 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%SetCoriolis(f0)
+  call modelobj%DiagnoseGeostrophicVelocity()
+
+  call modelobj%WriteModel()
+  call modelobj%IncrementIOCounter()
+
+  call modelobj%CalculateEntropy()
+  call modelobj%ReportEntropy()
+  call modelobj%ReportMetrics()
+  e0 = modelobj%entropy
+
+  ! Set the model's time integration method
+  call modelobj%SetTimeIntegrator(integrator)
+
+  ! forward step the model to `endtime` using a time step
+  ! of `dt` and outputing model data every `iointerval`
+  call modelobj%ForwardStep(endtime,dt,iointerval)
+
+  ef = modelobj%entropy
+
+  print*,e0,ef
+  if(abs(ef-e0) > epsilon(e0)) then
+    print*,"Warning: Final entropy greater than initial entropy! ",e0,ef
+  endif
+
+  ! Clean up
+  call modelobj%free()
+  call mesh%free()
+  call geometry%free()
+  call interp%free()
+
+endprogram linear_shallow_water2d_kelvinwaves