diff --git a/docs/Models/linear-shallow-water-model.md b/docs/Models/linear-shallow-water-model.md index a9629f9a6..025f6587f 100644 --- a/docs/Models/linear-shallow-water-model.md +++ b/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 diff --git a/docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md b/docs/Tutorials/LinearShallowWater/PlanetaryRossbyWave.md new file mode 100644 index 000000000..eb0b4de39 --- /dev/null +++ b/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$ + + diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 888824053..ef6ba7d42 100644 --- a/examples/CMakeLists.txt +++ b/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" ) diff --git a/examples/linear_shallow_water2d_kelvinwaves.f90 b/examples/linear_shallow_water2d_kelvinwaves.f90 new file mode 100644 index 000000000..f7965fb5b --- /dev/null +++ b/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 diff --git a/examples/linear_shallow_water2d_planetaryrossby_wave.f90 b/examples/linear_shallow_water2d_planetaryrossby_wave.f90 new file mode 100644 index 000000000..5ddc20660 --- /dev/null +++ b/examples/linear_shallow_water2d_planetaryrossby_wave.f90 @@ -0,0 +1,123 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! 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_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 = 500.0_prec !1000000.0_prec ! (s); 1-day + real(prec),parameter :: f0 = 10.0_prec**(-4) ! reference coriolis parameter (1/s) + real(prec),parameter :: beta = 10.0_prec**(-11) ! beta parameter (1/ms) + real(prec),parameter :: iointerval = 500.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*,"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_rossbywave diff --git a/examples/shallow_water_plot.py b/examples/shallow_water_plot.py new file mode 100644 index 000000000..7eabea18e --- /dev/null +++ b/examples/shallow_water_plot.py @@ -0,0 +1,99 @@ +# //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// # +# +# 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. +# +# //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// # +# +# 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() + + + diff --git a/mkdocs.yml b/mkdocs.yml index b707fd103..2c783fbb8 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 + - Planetary Rossby wave: Tutorials/LinearShallowWater/PlanetaryRossbyWave.md - Models: - Viscous Burgers Equation: Models/burgers-equation-model.md - Linear Euler (2D) : Models/linear-euler-2d-model.md diff --git a/pyself/geometry.py b/pyself/geometry.py index 334c147e7..d00be626b 100644 --- a/pyself/geometry.py +++ b/pyself/geometry.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -import self.lagrange as lagrange +import pyself.lagrange as lagrange # class semline: diff --git a/pyself/model2d.py b/pyself/model2d.py index 23a9262c1..03186ed89 100644 --- a/pyself/model2d.py +++ b/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 \ No newline at end of file diff --git a/src/SELF_DGModel1D_t.f90 b/src/SELF_DGModel1D_t.f90 index a0d6b5f14..a82f55cbe 100644 --- a/src/SELF_DGModel1D_t.f90 +++ b/src/SELF_DGModel1D_t.f90 @@ -106,6 +106,8 @@ subroutine Init_DGModel1D_t(this,mesh,geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) + call this%AdditionalInit() + call this%SetMetadata() endsubroutine Init_DGModel1D_t @@ -143,6 +145,7 @@ subroutine Free_DGModel1D_t(this) call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() + call this%AdditionalFree() endsubroutine Free_DGModel1D_t diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index dfb4a041a..c3e0f84aa 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/src/SELF_DGModel2D_t.f90 @@ -113,6 +113,8 @@ subroutine Init_DGModel2D_t(this,mesh,geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) + call this%AdditionalInit() + call this%SetMetadata() endsubroutine Init_DGModel2D_t @@ -145,6 +147,7 @@ subroutine Free_DGModel2D_t(this) call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() + call this%AdditionalFree() endsubroutine Free_DGModel2D_t @@ -551,7 +554,6 @@ subroutine CalculateTendency_DGModel2D_t(this) call this%SetBoundaryCondition() ! User-supplied if(this%gradient_enabled) then - call this%solution%AverageSides() call this%CalculateSolutionGradient() call this%SetGradientBoundaryCondition() ! User-supplied call this%solutionGradient%AverageSides() diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index 33e9d77c3..dd6a303ef 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -113,6 +113,8 @@ subroutine Init_DGModel3D_t(this,mesh,geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) + call this%AdditionalInit() + call this%SetMetadata() endsubroutine Init_DGModel3D_t @@ -145,6 +147,7 @@ subroutine Free_DGModel3D_t(this) call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() + call this%AdditionalFree() endsubroutine Free_DGModel3D_t diff --git a/src/SELF_LinearShallowWater2D_t.f90 b/src/SELF_LinearShallowWater2D_t.f90 index ec41cdcd0..1dfd75bfd 100644 --- a/src/SELF_LinearShallowWater2D_t.f90 +++ b/src/SELF_LinearShallowWater2D_t.f90 @@ -34,18 +34,47 @@ 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 + type(MappedScalar2D) :: fCori ! The coriolis parameter contains + procedure :: AdditionalInit => AdditionalInit_LinearShallowWater2D_t + procedure :: AdditionalFree => AdditionalFree_LinearShallowWater2D_t procedure :: SetNumberOfVariables => SetNumberOfVariables_LinearShallowWater2D_t procedure :: SetMetadata => SetMetadata_LinearShallowWater2D_t procedure :: entropy_func => entropy_func_LinearShallowWater2D_t procedure :: flux2d => flux2d_LinearShallowWater2D_t procedure :: riemannflux2d => riemannflux2d_LinearShallowWater2D_t procedure :: hbc2d_NoNormalFlow => hbc2d_NoNormalFlow_LinearShallowWater2D_t + procedure :: sourcemethod => sourcemethod_LinearShallowWater2D_t + ! Custom methods + generic,public :: SetCoriolis => SetCoriolis_fplane_LinearShallowWater2D_t, & + SetCoriolis_betaplane_LinearShallowWater2D_t + procedure,private :: SetCoriolis_fplane_LinearShallowWater2D_t + procedure,private :: SetCoriolis_betaplane_LinearShallowWater2D_t + + procedure,public :: DiagnoseGeostrophicVelocity => DiagnoseGeostrophicVelocity_LinearShallowWater2D_t + endtype LinearShallowWater2D_t contains + subroutine AdditionalInit_LinearShallowWater2D_t(this) + implicit none + class(LinearShallowWater2D_t),intent(inout) :: this + + call this%fCori%Init(this%geometry%x%interp, & + 1,this%mesh%nElem) + + endsubroutine AdditionalInit_LinearShallowWater2D_t + + subroutine AdditionalFree_LinearShallowWater2D_t(this) + implicit none + class(LinearShallowWater2D_t),intent(inout) :: this + + call this%fCori%Free() + + endsubroutine AdditionalFree_LinearShallowWater2D_t + subroutine SetNumberOfVariables_LinearShallowWater2D_t(this) implicit none class(LinearShallowWater2D_t),intent(inout) :: this @@ -54,16 +83,89 @@ subroutine SetNumberOfVariables_LinearShallowWater2D_t(this) endsubroutine SetNumberOfVariables_LinearShallowWater2D_t + subroutine SetCoriolis_fplane_LinearShallowWater2D_t(this,f0) + implicit none + class(LinearShallowWater2D_t),intent(inout) :: this + real(prec),intent(in) :: f0 + ! Local + integer :: iel + integer :: i + integer :: j + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) + this%fCori%interior(i,j,iel,1) = f0 + enddo + call this%fCori%UpdateDevice() + + endsubroutine SetCoriolis_fplane_LinearShallowWater2D_t + + subroutine SetCoriolis_betaplane_LinearShallowWater2D_t(this,f0,beta) + implicit none + class(LinearShallowWater2D_t),intent(inout) :: this + real(prec),intent(in) :: f0 + real(prec),intent(in) :: beta + ! Local + integer :: iel + integer :: i + integer :: j + real(prec) :: y + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) + y = this%geometry%x%interior(i,j,iel,1,2) + this%fCori%interior(i,j,iel,1) = f0+beta*y + enddo + call this%fCori%UpdateDevice() + + endsubroutine SetCoriolis_betaplane_LinearShallowWater2D_t + + subroutine DiagnoseGeostrophicVelocity_LinearShallowWater2D_t(this) + implicit none + class(LinearShallowWater2D_t),intent(inout) :: this + ! Local + integer :: iel + integer :: i + integer :: j + real(prec) :: dpdx,dpdy,f + + ! We assume here that the velocity field is identically zero + ! everywhere and the only field that is set is the free surface height + ! with a non-zero coriolis parameter. + ! In this case, we have that the tendency calculation will give + ! the gradient in the free surface, consistent with the DG approximation + this%solution%interior(:,:,:,1) = 0.0_prec ! Set u=0 + this%solution%interior(:,:,:,2) = 0.0_prec ! Set v=0 + call this%solution%UpdateDevice() + call this%CalculateTendency() + call this%dSdt%UpdateHost() + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) + + dpdx = -this%dSdt%interior(i,j,iel,1) + dpdy = -this%dSdt%interior(i,j,iel,2) + f = this%fCori%interior(i,j,iel,1) + this%solution%interior(i,j,iel,1) = -dpdy/f ! u + this%solution%interior(i,j,iel,2) = dpdx/f ! v + enddo + + call this%solution%UpdateDevice() + + endsubroutine DiagnoseGeostrophicVelocity_LinearShallowWater2D_t + subroutine SetMetadata_LinearShallowWater2D_t(this) implicit none class(LinearShallowWater2D_t),intent(inout) :: this call this%solution%SetName(1,"u") - call this%solution%SetUnits(1,"[null]") + call this%solution%SetUnits(1,"m/s") call this%solution%SetName(2,"v") - call this%solution%SetUnits(2,"[null]") + call this%solution%SetUnits(2,"m/s") call this%solution%SetName(3,"eta") - call this%solution%SetUnits(3,"[null]") + call this%solution%SetUnits(3,"m") + call this%fCori%SetName(1,"f") + call this%fCori%SetUnits(1,"1/s") endsubroutine SetMetadata_LinearShallowWater2D_t @@ -130,4 +232,25 @@ pure function hbc2d_NoNormalFlow_LinearShallowWater2D_t(this,s,nhat) result(exts endfunction hbc2d_NoNormalFlow_LinearShallowWater2D_t + subroutine sourcemethod_LinearShallowWater2D_t(this) + implicit none + class(LinearShallowWater2D_t),intent(inout) :: this + ! Local + integer :: iel + integer :: i + integer :: j + real(prec) :: s(1:this%nvar) + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iel=1:this%mesh%nElem) + + 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 + + enddo + + endsubroutine sourcemethod_LinearShallowWater2D_t + endmodule self_LinearShallowWater2D_t diff --git a/src/gpu/SELF_LinearShallowWater2D.cpp b/src/gpu/SELF_LinearShallowWater2D.cpp index 589ab46d3..0cc80fabf 100644 --- a/src/gpu/SELF_LinearShallowWater2D.cpp +++ b/src/gpu/SELF_LinearShallowWater2D.cpp @@ -142,4 +142,27 @@ extern "C" setboundarycondition_LinearShallowWater2D_gpukernel<<>>(extBoundary,boundary,sideInfo,nhat,N,nel,nvar); } +} + +__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, 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 + + } + +} +extern "C" +{ + void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, 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); + } } \ No newline at end of file diff --git a/src/gpu/SELF_LinearShallowWater2D.f90 b/src/gpu/SELF_LinearShallowWater2D.f90 index 55e06f6e7..97b31e8b7 100644 --- a/src/gpu/SELF_LinearShallowWater2D.f90 +++ b/src/gpu/SELF_LinearShallowWater2D.f90 @@ -35,7 +35,7 @@ module self_LinearShallowWater2D procedure :: setboundarycondition => setboundarycondition_LinearShallowWater2D procedure :: boundaryflux => boundaryflux_LinearShallowWater2D procedure :: fluxmethod => fluxmethod_LinearShallowWater2D - ! 'procedure :: sourcemethod => sourcemethod_LinearShallowWater2D + procedure :: sourcemethod => sourcemethod_LinearShallowWater2D endtype LinearShallowWater2D @@ -70,6 +70,16 @@ subroutine fluxmethod_LinearShallowWater2D_gpu(solution,flux,g,H,N,nel,nvar) & endsubroutine fluxmethod_LinearShallowWater2D_gpu endinterface + interface + subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,N,nel,nvar) & + bind(c,name="sourcemethod_LinearShallowWater2D_gpu") + use iso_c_binding + use SELF_Constants + type(c_ptr),value :: solution,source,fCori + integer(c_int),value :: N,nel,nvar + endsubroutine sourcemethod_LinearShallowWater2D_gpu + endinterface + contains subroutine boundaryflux_LinearShallowWater2D(this) @@ -148,12 +158,17 @@ subroutine fluxmethod_LinearShallowWater2D(this) endsubroutine fluxmethod_LinearShallowWater2D - ! subroutine sourcemethod_LinearShallowWater2D(this) - ! implicit none - ! class(LinearShallowWater2D),intent(inout) :: this + subroutine sourcemethod_LinearShallowWater2D(this) + implicit none + class(LinearShallowWater2D),intent(inout) :: this - ! return + call sourcemethod_LinearShallowWater2D_gpu(this%solution%interior_gpu, & + this%source%interior_gpu, & + this%fCori%interior_gpu, & + this%solution%interp%N, & + this%solution%nelem, & + this%solution%nvar) - ! endsubroutine sourcemethod_LinearShallowWater2D + endsubroutine sourcemethod_LinearShallowWater2D endmodule self_LinearShallowWater2D diff --git a/src/gpu/SELF_MappedScalar_2D.f90 b/src/gpu/SELF_MappedScalar_2D.f90 index c396e74a5..9e705a0db 100644 --- a/src/gpu/SELF_MappedScalar_2D.f90 +++ b/src/gpu/SELF_MappedScalar_2D.f90 @@ -118,6 +118,8 @@ subroutine Init_MappedScalar2D(this,interp,nVar,nElem) call hipblasCheck(hipblasCreate(this%blas_handle)) + call this%UpdateDevice() + endsubroutine Init_MappedScalar2D subroutine Free_MappedScalar2D(this) diff --git a/test/mappedscalarbrgradient_2d_linear_structuredmesh.f90 b/test/mappedscalarbrgradient_2d_linear_structuredmesh.f90 index 6994cf977..86ecf145b 100644 --- a/test/mappedscalarbrgradient_2d_linear_structuredmesh.f90 +++ b/test/mappedscalarbrgradient_2d_linear_structuredmesh.f90 @@ -131,28 +131,6 @@ integer function mappedscalarbrgradient_2d_linear() result(r) call f%SetName(1,"f") call f%SetUnits(1,"[null]") - call Open_HDF5('output.h5',H5F_ACC_TRUNC_F,fileId) - - ! Write the interpolant to the file - print*,"Writing interpolant data to file" - call f%interp%WriteHDF5(fileId) - - ! Write the model state to file - print*,"Writing control grid solution to file" - call CreateGroup_HDF5(fileId,'/controlgrid') - call f%WriteHDF5(fileId,'/controlgrid/solution') - - print*,"Writing control grid solution gradient to file" - call CreateGroup_HDF5(fileId,'/controlgrid') - call df%WriteHDF5(fileId,'/controlgrid/solution_gradient') - - ! Write the geometry to file - print*,"Writing control grid geometry to file" - call CreateGroup_HDF5(fileId,'/controlgrid/geometry') - call geometry%x%WriteHDF5(fileId,'/controlgrid/geometry/x') - - call Close_HDF5(fileId) - ! Calculate diff from exact do iel = 1,mesh%nelem do j = 1,controlDegree+1