add multilayer SWE in 2D

examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl

@@ -0,0 +1,61 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with three layers
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0,
+                                              rhos = (0.9, 1.0, 1.1))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_ersing_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 = 4,
+                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_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl

@@ -0,0 +1,92 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break test with a 
+# discontinuous bottom topography function to test entropy conservation
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0,
+                                              rhos = (0.9, 0.95, 1.0))
+
+# This academic testcase sets up a discontinuous bottom topography 
+# function and initial condition to test entropy conservation. 
+
+function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D)
+    # Bottom topography
+    b = 0.3 * exp(-0.5 * ((x[1])^2 + (x[2])^2))
+
+    if x[1] < 0.0
+        H = SVector(1.0, 0.8, 0.6)
+    else
+        H = SVector(0.9, 0.7, 0.5)
+        b += 0.1
+    end
+
+    v1 = zero(H)
+    v2 = zero(H)
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break)
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_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 = (-1.0, -1.0)
+coordinates_max = (1.0, 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 solver
+
+tspan = (0.0, 2.0)
+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))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.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

examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl

@@ -0,0 +1,79 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations with a bottom topography function 
+# to test well-balancedness
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6,
+                                              rhos = (0.7, 0.8, 0.9, 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::ShallowWaterMultiLayerEquations2D)
+    H = SVector(0.6, 0.55, 0.5, 0.45)
+    v1 = zero(H)
+    v2 = zero(H)
+    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..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_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

examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl

@@ -0,0 +1,62 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the two-layer shallow water equations with a periodic
+# bottom topography function (set in the initial conditions)
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 10.0,
+                                              rhos = (0.9, 1.0, 1.1))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 8, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# This setup is for the curved, split form convergence test on a periodic domain
+
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
+                           joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = true)
+
+# Create the semidiscretization 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/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl

@@ -0,0 +1,97 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break test with a 
+# discontinuous bottom topography function to test entropy conservation
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0,
+                                              rhos = (0.9, 0.95, 1.0))
+
+# This academic testcase sets up a discontinuous bottom topography 
+# function and initial condition to test entropy conservation. 
+
+function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D)
+    # Bottom topography
+    b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2))
+
+    if x[1] < sqrt(2) / 2
+        H = SVector(1.0, 0.8, 0.6)
+    else
+        H = SVector(0.9, 0.7, 0.5)
+        b += 0.1
+    end
+
+    v1 = zero(H)
+    v2 = zero(H)
+    return prim2cons(SVector(H..., v1..., v2..., b),
+                     equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break)
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 6, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
+                           joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh"))
+
+mesh = UnstructuredMesh2D(mesh_file, periodicity = false)
+
+# Boundary conditions
+boundary_condition = Dict(:Top => boundary_condition_slip_wall,
+                          :Left => boundary_condition_slip_wall,
+                          :Right => boundary_condition_slip_wall,
+                          :Bottom => boundary_condition_slip_wall)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
+                                    solver, boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 0.5)
+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))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.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