[WIP] Wet/dry capabilities for 2D shallow water equations

examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl

@@ -0,0 +1,100 @@
+
+using OrdinaryDiffEq
+using Trixi
+using Plots
+
+ ###############################################################################
+ # Semidiscretization of the shallow water equations
+equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=1.4)
+cfl = 0.5                                    
+
+function initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D)
+  # Set the background values
+
+  v1 = 0.0
+  v2 = 0.0
+
+  x1, x2 = x
+  b = max(0.1, - sqrt(16*x1^2+16*x2^2) + 1)
+
+  # use a logistic function to tranfer 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))))
+
+  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_cn, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(4)
+
+indicator_sc = IndicatorHennemannGassner(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., -1.)
+coordinates_max = (1.,  1.)
+
+cells_per_dimension = (20, 20)
+
+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.)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_integrals=(lake_at_rest_error,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=100,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,),
+                                                     variables=(Trixi.waterheight,))
+
+sol = solve(ode, SSPRK43(stage_limiter!),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, adaptive=false);
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl

@@ -0,0 +1,109 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations2D(gravity_constant=9.81)
+cfl = 0.6
+
+# Implemented based on Wintermeyer (scaled on half of domain)
+function parabolic_bowl_analytic_2D_H(gravity, x,t)
+  a = 1
+  h_0 = 0.1
+  σ = 0.5
+  ω = sqrt(2*gravity*h_0)/a
+
+  H = 0.5 * (σ * h_0/a^2 * (2*2*x[1]*cos(ω*t) + 2*2*x[2]*sin(ω*t)- σ) + h_0)
+  return H
+end
+
+# Implemented based on Wintermeyer (scaled on half of domain)
+function initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D)
+  a = 1
+  h_0 = 0.1
+  σ = 0.5
+  ω = sqrt(2*equations.gravity*h_0)/a
+
+  v1 = -σ*ω*sin(ω*t)
+  v2 = -σ*ω*cos(ω*t)
+
+  b = 0.5 * (h_0 * ((2*x[1])^2 + (2*x[2])^2)/a^2)
+
+  H = max(b, parabolic_bowl_analytic_2D_H(equations.gravity, x, 0))
+
+  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_cn, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(6)
+
+indicator_sc = IndicatorHennemannGassner(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 = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+cells_per_dimension = (60,60)
+
+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.)
+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)
+
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,),
+                                                     variables=(Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0,
+            save_everystep=false, callback=callbacks, adaptive=false);
+
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl

@@ -0,0 +1,103 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+
+equations = ShallowWaterEquations2D(gravity_constant=9.812)
+cfl = 1
+
+function initial_condition_complex_bottom_wb_CN(x, t, equations:: ShallowWaterEquations2D)
+  v1 = 0
+  v2 = 0
+  b = sin(4*pi*x[1]) + 3
+
+  if x[1] >= 0.5
+    b = sin(4*pi*x[1]) + 1
+  end
+
+  H = max(b, 2.5)
+
+  if x[1] >= 0.5
+    H = max(b, 1.5)
+  end
+
+  H = max(H, b + equations.threshold_limiter)
+  return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_complex_bottom_wb_CN
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+
+surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle),
+                flux_nonconservative_chen_noelle)
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassner(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 StructuredMesh for the domain [0, 1]^2
+
+coordinates_min = (0., 0.)
+coordinates_max = (1., 1.)
+
+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 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)
+
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,),
+                                                     variables=(Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0,
+            save_everystep=false, callback=callbacks, adaptive=false);
+
+summary_callback() # print the timer summary