Add rossby wave demo

*Currently disabled the coriolis source terms to debug an issue on mi210

Author: fluidnumerics-joe

docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md

@@ -0,0 +1,63 @@
+# 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$
+
+## Runtimes
+
+To benchmark this example, we run the simulation for 10 days of simulation time, using a time-step of a half-second.
+The solution after 10 days looks like the image shown below.
+
+![Rossby Wave Initial Condition](./rossbywave_day10.png){ align=center }
+
+```
+$SELF_PREFIX/examples/LinearShallowWater_Rossbywave -dt 0.5 -tn 172800 -int "rk3"
+```
+

examples/CMakeLists.txt

@@ -66,5 +66,6 @@ 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_euler3d_spherical_soundwave_radiation.f90"
     )

examples/linear_shallow_water2d_planetaryrossby_wave.f90

@@ -0,0 +1,124 @@
+! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// !
+!
+! Maintainers : [email protected]
+! 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_rossbywave
+  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.5_prec ! Time-step size
+  real(prec),parameter :: endtime = 1000.0_prec ! (s); 1-day
+  real(prec),parameter :: f0 = 0.0_prec !10.0_prec**(-4) ! reference coriolis parameter (1/s)
+  real(prec),parameter :: beta = 0.0_prec ! 10.0_prec**(-11) ! beta parameter (1/ms)
+  real(prec),parameter :: iointerval = 10.0 ! Write files 10 times per day
+
+  real(prec) :: e0,ef ! Initial and final entropy
+  type(LinearShallowWater2D) :: modelobj ! Shallow water model
+  type(Lagrange),target :: interp ! Interpolant
+  integer :: bcids(1:4) ! Boundary conditions for structured mesh
+  type(Mesh2D),target :: mesh ! Mesh class
+  type(SEMQuad),target :: geometry ! Geometry class
+
+  real(prec),parameter :: g = 10.0_prec ! Acceleration due to gravity
+  real(prec),parameter :: H = 1000.0_prec ! Uniform resting depth
+  real(prec),parameter :: dx = 10.0_prec**(5) ! grid spacing in x-direction (m)
+  real(prec),parameter :: dy = 10.0_prec**(5) ! grid spacing in y-direction (m)
+  integer,parameter :: ntx = 2 ! number of tiles in x-direction
+  integer,parameter :: nty = 2 ! number of tiles in y-direction
+  integer,parameter :: nxp = 5 ! number of element per tile in x-direction
+  integer,parameter :: nyp = 5 ! number of element per tile in y-direction
+
+  ! Set no normal flow boundary conditions
+  bcids(1:4) = [SELF_BC_NONORMALFLOW, & ! South
+                SELF_BC_NONORMALFLOW, & ! East
+                SELF_BC_NONORMALFLOW, & ! North
+                SELF_BC_NONORMALFLOW] ! West
+
+  ! Create a uniform block mesh
+  call mesh%StructuredMesh(nxp,nyp,ntx,nty,dx,dy,bcids)
+
+  ! 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(1,'f = 0')
+  call modelobj%solution%SetEquation(2,'f = 0')
+  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)
+
+  ! call modelobj%SetCoriolis(f0,beta)
+  !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*,"Error: Final entropy greater than initial entropy! ",e0,ef
+    stop 1
+  endif
+
+  ! Clean up
+  call modelobj%free()
+  call mesh%free()
+  call geometry%free()
+  call interp%free()
+
+endprogram linear_shallow_water2d_rossbywave

examples/shallow_water_plot.py

@@ -0,0 +1,99 @@
+# //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #
+#
+# Maintainers : [email protected]
+# 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.
+#
+# //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// #
+#
+# NOTE #
+#
+# If you encounter the following error on debian/ubuntu:
+#
+# ```
+# libGL error: MESA-LOADER: failed to open swrast: /usr/lib/dri/swrast_dri.so: cannot open shared object file: No such file or directory (search paths /usr/lib/x86_64-linux-gnu/dri:\$${ORIGIN}/dri:/usr/lib/dri, suffix _dri)
+# libGL error: failed to load driver: swrast
+# 2024-10-18 01:47:34.339 (   3.885s) [    754DF00AD740]vtkXOpenGLRenderWindow.:651    ERR| vtkXOpenGLRenderWindow (0x3b03f80): Cannot create GLX context.  Aborting.
+# ERROR:root:Cannot create GLX context.  Aborting.
+# ```
+# See https://askubuntu.com/questions/1352158/libgl-error-failed-to-load-drivers-iris-and-swrast-in-ubuntu-20-04 for resolution
+
+import numpy as np
+import matplotlib
+import pyself.model2d as model2d
+import pyvista as pv
+import glob
+import os
+import sys
+import subprocess
+
+# Specify the directory you want to search in
+directory_path = "/home/joe/SELF/build/examples/" 
+# output video name
+video_name = "rossbywave.mp4"
+
+etamin = -1.0e-3
+etamax = 1.0e-3
+
+colormap = matplotlib.colormaps['RdBu']
+
+def get_sorted_files(pattern):
+    files = glob.glob(pattern)
+    files.sort(key=os.path.getmtime)  # Sort by modification time
+    return files
+
+# Filter the list to include only .dat files
+pickup_files = get_sorted_files(f"{directory_path}/solution.*.h5")
+
+# If you want to get the full path of each file
+pickup_files = [os.path.join(directory_path, f) for f in pickup_files]
+
+model = model2d.model()
+pv.start_xvfb()
+
+k = 0
+# Example usage of reading each file
+for pickup_file in pickup_files:
+    print(pickup_file)
+    if(k==0):
+        model.load(pickup_file)
+        pl = pv.Plotter(off_screen=True)
+        pl.add_mesh(model.pvdata,
+                    scalars="eta",
+                    cmap=colormap,
+                    clim=[etamin,etamax])
+        # pl.add_mesh(model.pvdata,
+        #             scalars="u",
+        #             cmap=colormap,
+        #             clim=[-1e-3, 0])
+        pl.show(auto_close=False)
+        pl.camera_position = 'xy'
+        pl.open_movie(video_name)
+    else:
+        model.update_from_file(pickup_file)
+
+    pl.write_frame()
+    k+=1
+
+pl.close()
+
+
+

pyself/geometry.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 #
 
-import self.lagrange as lagrange
+import pyself.lagrange as lagrange
 
 
 # class semline:

pyself/model2d.py

@@ -136,3 +136,58 @@ def generate_pyvista(self):
                     k+=1
 
         print(self.pvdata)
+
+    def update_from_file(self, hdf5File):
+        """Loads in 2-D model from SELF model output"""
+        import h5py
+        import dask.array as da
+
+        f = h5py.File(hdf5File, 'r')
+            
+        if 'controlgrid' in list(f.keys()):
+
+            controlgrid = f['controlgrid']
+            for group_name in controlgrid.keys():
+
+                if( group_name == 'geometry' ):
+                    continue
+
+                group = controlgrid[group_name]
+                # Create a list to hold data for this group
+                group_data = getattr(self, group_name)
+                print(f"Loading {group_name} group")
+                
+                for k in group.keys():
+                    k_decoded = k.encode('utf-8').decode('utf-8')
+                    if k == 'metadata':
+                        continue
+                    else:
+                        print(f"Loading {k_decoded} field")
+                        # Load the actual data
+                        d = group[k]
+                        N = d.shape[2]
+
+                        # Find index for this field
+                        i = 0
+                        for sol in group_data:
+                            if sol['name'] == k_decoded:
+                                break
+                            else:
+                                i += 1
+
+                        group_data[i]['data'] = da.from_array(d, chunks=(self.geom.daskChunkSize, N, N))
+
+            # # Load fields into pvdata
+            k = 0
+            for attr in self.__dict__:
+                if not attr in ['pvdata','varnames','varunits','geom']:
+                    controlgroup = getattr(self, attr)
+                    for var in controlgroup:
+                        self.pvdata.point_data.set_array(var['data'].flatten(),var['name'])
+                        k+=1
+
+        else:
+            print(f"Error: /controlgrid group not found in {hdf5File}.")
+            return 1
+
+        return 0