Move two-layer shallow water equations from Trixi.jl

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl

@@ -0,0 +1,61 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations
+
+equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 10.0, rho_upper = 0.9,
+                                            rho_lower = 1.0)
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = sqrt(2.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# create the semi discretization object
+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 = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl

@@ -0,0 +1,95 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations for a dam break
+# test with a discontinuous bottom topography function to test entropy conservation
+
+equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 9.81, H0 = 2.0,
+                                            rho_upper = 0.9, rho_lower = 1.0)
+
+# Initial condition of a dam break with a discontinuous water heights and bottom topography.
+# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations1D)
+    v1_upper = 0.0
+    v1_lower = 0.0
+
+    # Set the discontinuity
+    if x[1] <= 10.0
+        H_lower = 2.0
+        H_upper = 4.0
+        b = 0.0
+    else
+        H_lower = 1.5
+        H_upper = 3.0
+        b = 0.5
+    end
+
+    return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a non-periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 20.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 5,
+                n_cells_max = 10000,
+                periodicity = false)
+
+boundary_condition = boundary_condition_slip_wall
+
+# create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solvers
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (energy_total,
+                                                                 energy_kinetic,
+                                                                 energy_internal))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

@@ -0,0 +1,87 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations to test well-balancedness
+
+equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 1.0, H0 = 0.6,
+                                            rho_upper = 0.9, rho_lower = 1.0)
+
+"""
+    initial_condition_fjordholm_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations1D)
+
+Initial condition to test well balanced with a bottom topography from Fjordholm
+- Ulrik Skre Fjordholm (2012)
+  Energy conservative and stable schemes for the two-layer shallow water equations.
+  [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039)
+"""
+function initial_condition_fjordholm_well_balanced(x, t,
+                                                   equations::ShallowWaterTwoLayerEquations1D)
+    inicenter = 0.5
+    x_norm = x[1] - inicenter
+    r = abs(x_norm)
+
+    H_lower = 0.5
+    H_upper = 0.6
+    v1_upper = 0.0
+    v1_lower = 0.0
+    b = r <= 0.1 ? 0.2 * (cos(10 * pi * (x[1] - 0.5)) + 1) : 0.0
+    return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations)
+end
+
+initial_condition = initial_condition_fjordholm_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# 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 = (lake_at_rest_error,))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+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,
+                        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

examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl

@@ -0,0 +1,61 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations
+
+equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 10.0, rho_upper = 0.9,
+                                            rho_lower = 1.0)
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (sqrt(2.0), sqrt(2.0))
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 20_000,
+                periodicity = true)
+
+# Create the semi discretization object
+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 = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

@@ -0,0 +1,82 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations with a bottom topography function 
+# to test well-balancedness
+
+equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 9.81, H0 = 0.6,
+                                            rho_upper = 0.9, rho_lower = 1.0)
+
+# An initial condition with constant total water height, zero velocities and a bottom topography to 
+# test well-balancedness
+function initial_condition_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations2D)
+    H_lower = 0.5
+    H_upper = 0.6
+    v1_upper = 0.0
+    v2_upper = 0.0
+    v1_lower = 0.0
+    v2_lower = 0.0
+    b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ?
+         0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] +
+                                                    -0.5)^2)) + 1) : 0.0)
+
+    return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b),
+                     equations)
+end
+
+initial_condition = initial_condition_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (lake_at_rest_error,))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+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,
+                        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