Merge branch 'sg/wet-dry-swe-1d' into sg/wet-dry-swe-2d

examples/tree_1d_dgsem/elixir_shallowwater_beach.jl

@@ -0,0 +1,118 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations1D(gravity_constant=9.812)
+
+"""
+    initial_condition_beach(x, t, equations:: ShallowWaterEquations1D)
+Initial condition to simulate a wave running towards a beach and crashing. Difficult test
+including both wetting and drying in the domain using slip wall boundary conditions.
+The water height and speed function used here, are an adaption of the initial condition
+found in section 5.2 of the paper:
+  - Andreas Bollermann, Sebastian Noelle, Maria Lukáčová-Medvid’ová (2011)
+    Finite volume evolution Galerkin methods for the shallow water equations with dry beds\n
+    [DOI: 10.4208/cicp.220210.020710a](https://dx.doi.org/10.4208/cicp.220210.020710a)
+The used bottom topography differs from the one out of the paper to be differentiable
+"""
+function initial_condition_beach(x, t, equations:: ShallowWaterEquations1D)
+  D = 1
+  delta = 0.02
+  gamma = sqrt((3 * delta) / (4 * D))
+  x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20))
+
+  f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2
+
+  # steep wall
+  b = 0.01 + 99/409600 * 4^x[1]
+
+  if x[1] >= 6
+    H = b
+    v = 0.0
+  else
+    H = f
+    v = sqrt(equations.gravity / D) * H
+  end
+
+  # 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 v) / 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 ShallowWaterEquations 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, v, b), equations)
+end
+
+initial_condition = initial_condition_beach
+boundary_condition = boundary_condition_slip_wall
+
+###############################################################################
+# 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(3)
+
+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)
+
+###############################################################################
+# Create the TreeMesh for the domain [0, 8]
+
+coordinates_min = 0.0
+coordinates_max = 8.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=7,
+                n_cells_max=10_000,
+                periodicity=false)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_condition)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.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,
+                                                               lake_at_rest_error))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(dt=0.5,
+                                     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/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl

@@ -0,0 +1,115 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations1D(gravity_constant=9.81)
+
+"""
+    initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D)
+
+Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its
+wet-dry mechanics. This test has analytical solutions. 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 in two dimensions 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:: ShallowWaterEquations1D)
+  a = 1
+  h_0 = 0.1
+  sigma = 0.5
+  ω = sqrt(2 * equations.gravity * h_0) / a
+
+  v = -sigma * ω * sin(ω * t)
+
+  b = h_0 * x[1]^2 / a^2
+
+  H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * 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 v) / 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 ShallowWaterEquations 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, v, 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(5)
+
+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)
+
+###############################################################################
+# Create the TreeMesh for the domain [-2, 2]
+
+coordinates_min = -2.0
+coordinates_max = 2.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=6,
+                n_cells_max=10_000)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.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,
+                                                               lake_at_rest_error))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=1000,
+                                     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