diff --git a/docs/Models/linear-shallow-water-model.md b/docs/Models/linear-shallow-water-model.md index 025f6587f..11e29272a 100644 --- a/docs/Models/linear-shallow-water-model.md +++ b/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 \ No newline at end of file +* [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) \ No newline at end of file diff --git a/docs/Tutorials/LinearShallowWater/KelvinWaves.md b/docs/Tutorials/LinearShallowWater/KelvinWaves.md new file mode 100644 index 000000000..3504fc9ad --- /dev/null +++ b/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. + +
+ ![Geostrophic adjustment releasing unbalanced flows](./img/kelvin-wave-initial-erruption.png){ width="500" } +
Free surface height (eta) 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. +
+
+ +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. diff --git a/docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md b/docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md index eb0b4de39..5f96b7dab 100644 --- a/docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md +++ b/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 diff --git a/docs/Tutorials/LinearShallowWater/LinearShallowWater.md b/docs/Tutorials/LinearShallowWater/ReflectingWave.md similarity index 94% rename from docs/Tutorials/LinearShallowWater/LinearShallowWater.md rename to docs/Tutorials/LinearShallowWater/ReflectingWave.md index 5b31b0bec..5e1e91386 100644 --- a/docs/Tutorials/LinearShallowWater/LinearShallowWater.md +++ b/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$. -

- - Free surface height (eta) at time t=0. -

+
+ ![Initial condition](./img/shallow-water.0000.png){ width="500" } +
Free surface height (eta) at time t=0. +
+
+ +
+ ![Gravity wave reflection](./img/shallow-water.0019.png){ width="500" } +
Free surface height (eta) at time t=1. +
+
-

- - Free surface height (eta) at time t=1. -

## 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 -

-

- -
-
- Free surface height (eta) for the full duration (1 second) of the problem. -
-

+
+ ![Geostrophic adjustment releasing unbalanced flows](./img/shallow-water.gif){ width="500" } +
!!! 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. diff --git a/docs/Tutorials/LinearShallowWater/img/kelvin-wave-initial-erruption.png b/docs/Tutorials/LinearShallowWater/img/kelvin-wave-initial-erruption.png new file mode 100644 index 000000000..4cd45c239 Binary files /dev/null and b/docs/Tutorials/LinearShallowWater/img/kelvin-wave-initial-erruption.png differ diff --git a/examples/linear_shallow_water2d_kelvinwaves.f90 b/examples/linear_shallow_water2d_kelvinwaves.f90 index f7965fb5b..c4985ec6c 100644 --- a/examples/linear_shallow_water2d_kelvinwaves.f90 +++ b/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 diff --git a/mkdocs.yml b/mkdocs.yml index 2c783fbb8..b3838f01c 100644 --- a/mkdocs.yml +++ b/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 diff --git a/src/SELF_LinearShallowWater2D_t.f90 b/src/SELF_LinearShallowWater2D_t.f90 index 1dfd75bfd..69121d622 100644 --- a/src/SELF_LinearShallowWater2D_t.f90 +++ b/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 diff --git a/src/gpu/SELF_LinearShallowWater2D.cpp b/src/gpu/SELF_LinearShallowWater2D.cpp index 0cc80fabf..7d68dbeab 100644 --- a/src/gpu/SELF_LinearShallowWater2D.cpp +++ b/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<<>>(solution,source,fCori,ndof); + sourcemethod_LinearShallowWater2D_gpukernel<<>>(solution,source,fCori,Cd,ndof); } } \ No newline at end of file diff --git a/src/gpu/SELF_LinearShallowWater2D.f90 b/src/gpu/SELF_LinearShallowWater2D.f90 index 97b31e8b7..9ae072d4f 100644 --- a/src/gpu/SELF_LinearShallowWater2D.f90 +++ b/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)