Add 2D-SWE and move wet/dry functionality from Trixi.jl (#18)

* add equation SWEWetDry1D

* add tests and elixirs

* add OrdinaryDiffEq as a test dependency

* comment functions related to HR_chen_noelle

* fix tests

* fix tests #2

* fix tests #3

* fix tests 4

* fix tests 5

* fix tests 6

* add unit tests and namespace conflict test

* update unit tests, comment wave speed estimates

* update unit test

* add compat bounds

* add some comments

* clean some comments, remove wave speed estimates

* comment unused lines

* include examples folder for coverage test

* add upstream tests

* fix ci

* Adjust comment

* remove comments

* do using Trixi instead of using Trixi : Trixi

* update swe-1d

* add swe-2d

* add wet/dry functionality, tests and examples

* add unit tests

* Update comments in test files

* add upstream tests

* update comments

* fix typos

* add Downloads pkg, apply formatter

* move downloads dependency to test project.toml

* Update Trixi version to 0.7

* fix formatter to older version

* add Printf to Project.toml

* Change wavespeed estimate to fix tests with hll fluxes

* add unit tests, remove unused function

* apply formatter

* Change comments according to code review, add compat bound for LinearAlgebra

* reorganize single test files for each mesh type

* Add own flux hll and wave speed estimates; add unit tests

* remove unnecessary qualifier

* apply formatter

* Apply changes from code review and introduce new field equations.swe_trixi

* apply formatter

* Update name and type assignment of swe_trixi -> basic_swe

* add comment about duplication of g and H0

---------

Co-authored-by: Andrew Winters <[email protected]>

.github/workflows/FormatCheck.yml

@@ -25,8 +25,10 @@ jobs:
         # This will use the latest version by default but you can set the version like so:
         #
         # julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version = "0.13.0"))'
+        #
+        # TODO: remove version restriction once formatting has sorted itself out
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
           julia  -e 'using JuliaFormatter; format(".")'
       - name: Format check
         run: |

Project.toml

@@ -4,14 +4,18 @@ authors = ["Andrew R. Winters <[email protected]>", "Michael Schlottke-
 version = "0.1.0-pre"
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
+LinearAlgebra = "1"
 MuladdMacro = "0.2.2"
 Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
 StaticArrays = "1"
-Trixi = "0.5.17, 0.6"
+Trixi = "0.7"
 julia = "1.8"
+Printf = "1"

examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl

@@ -0,0 +1,113 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 1.4)
+
+"""
+    initial_condition_conical_island(x, t, equations::ShallowWaterEquationsWetDry2D)
+
+Initial condition for the [`ShallowWaterEquationsWetDry2D`](@ref) to test the [`hydrostatic_reconstruction_chen_noelle`](@ref)
+and its handling of discontinuous water heights at the start in combination with wetting and
+drying. The bottom topography is given by a conical island in the middle of the domain. Around that
+island, there is a cylindrical water column at t=0 and the rest of the domain is dry. This
+discontinuous water height is smoothed by a logistic function. This simulation uses periodic
+boundary conditions.
+"""
+function initial_condition_conical_island(x, t, equations::ShallowWaterEquationsWetDry2D)
+    # Set the background values
+
+    v1 = 0.0
+    v2 = 0.0
+
+    x1, x2 = x
+    b = max(0.1, 1.0 - 4.0 * sqrt(x1^2 + x2^2))
+
+    # use a logistic function to transfer water height value smoothly
+    L = equations.H0    # maximum of function
+    x0 = 0.3   # center point of function
+    k = -25.0 # sharpness of transfer
+
+    H = max(b, L / (1.0 + exp(-k * (sqrt(x1^2 + x2^2) - x0))))
+
+    # It is mandatory to shift the water level at dry areas to make sure the water height h
+    # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+    # the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold
+    # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above
+    # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0.
+    # This default value can be changed within the constructor call depending on the simulation setup.
+    H = max(H, b + equations.threshold_limiter)
+    return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_conical_island
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle,
+                                              hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(4)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max = 0.5,
+                                                     alpha_min = 0.001,
+                                                     alpha_smooth = true,
+                                                     variable = waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the StructuredMesh and setup a periodic mesh
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+cells_per_dimension = (16, 16)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback = callbacks);
+
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl

@@ -0,0 +1,118 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81)
+
+"""
+    initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquationsWetDry2D)
+
+Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its
+wet-dry mechanics. This test has an analytical solution. The initial condition is defined by the
+analytical solution at time t=0. The bottom topography defines a bowl and the water level is given
+by an oscillating lake.
+
+The original test and its analytical solution were first presented in
+- William C. Thacker (1981)
+  Some exact solutions to the nonlinear shallow-water wave equations
+  [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882).
+
+The particular setup below is taken from Section 6.2 of
+- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018)
+  An entropy stable discontinuous Galerkin method for the shallow water equations on
+  curvilinear meshes with wet/dry fronts accelerated by GPUs
+  [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038).
+"""
+function initial_condition_parabolic_bowl(x, t, equations::ShallowWaterEquationsWetDry2D)
+    a = 1.0
+    h_0 = 0.1
+    sigma = 0.5
+    ω = sqrt(2 * equations.gravity * h_0) / a
+
+    v1 = -sigma * ω * sin(ω * t)
+    v2 = sigma * ω * cos(ω * t)
+
+    b = h_0 * ((x[1])^2 + (x[2])^2) / a^2
+
+    H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) + 2 * x[2] * sin(ω * t) - sigma) + h_0
+
+    # It is mandatory to shift the water level at dry areas to make sure the water height h
+    # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+    # the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold
+    # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above
+    # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0.
+    # This default value can be changed within the constructor call depending on the simulation setup.
+    H = max(H, b + equations.threshold_limiter)
+    return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_parabolic_bowl
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle,
+                                              hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(4)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max = 0.6,
+                                                     alpha_min = 0.001,
+                                                     alpha_smooth = true,
+                                                     variable = waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+
+coordinates_min = (-2.0, -2.0)
+coordinates_max = (2.0, 2.0)
+
+cells_per_dimension = (150, 150)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (energy_kinetic,
+                                                                 energy_internal))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback = callbacks);
+
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl

@@ -0,0 +1,59 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (sqrt(2.0), sqrt(2.0))
+
+cells_per_dimension = (8, 8)
+
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 2.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary