Add multilayer shallow water equations in 2D

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl

@@ -6,7 +6,7 @@ using TrixiShallowWater
 ###############################################################################
 # Semidiscretization of the multilayer shallow water equations with three layers
 
-equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 10.0,
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.1,
                                               rhos = (0.9, 1.0, 1.1))
 
 initial_condition = initial_condition_convergence_test
@@ -25,7 +25,7 @@ solver = DGSEM(polydeg = 3,
 coordinates_min = 0.0
 coordinates_max = sqrt(2.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level = 4,
+                initial_refinement_level = 2,
                 n_cells_max = 10_000,
                 periodicity = true)
 
@@ -50,12 +50,16 @@ save_solution = SaveSolutionCallback(interval = 500,
                                      save_initial_solution = true,
                                      save_final_solution = true)
 
-callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
 
 ###############################################################################
 # 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,
+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_convergence.jl

@@ -0,0 +1,64 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with three layers
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.1,
+                                              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 = 2,
+                n_cells_max = 100_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)
+
+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_dam_break.jl

@@ -0,0 +1,134 @@
+
+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_conditions = boundary_condition_slip_wall
+
+###############################################################################
+# 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 non-periodic mesh with wall boundary conditions
+
+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 = false)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+###############################################################################
+#= 
+Workaround for TreeMesh2D to set true discontinuities for debugging and testing. 
+Essentially, this is a slight augmentation of the `compute_coefficients` where the `x` node values
+passed here are slightly perturbed in order to set a true discontinuity that avoids the doubled 
+value of the LGL nodes at a particular element interface.
+=#
+
+# Point to the data we want to augment
+u = Trixi.wrap_array(ode.u0, semi)
+# Reset the initial condition
+for element in eachelement(semi.solver, semi.cache)
+    for i in eachnode(semi.solver), j in eachnode(semi.solver)
+        x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations,
+                                       semi.solver, i, j, element)
+        # Changing the node positions passed to the initial condition by the minimum
+        # amount possible with the current type of floating point numbers allows setting
+        # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump`
+        # works if the jump location `x_jump` is at the position of an interface.
+        if i == 1 && j == 1 # bottom left corner
+            x_node = SVector(nextfloat(x_node[1]), nextfloat(x_node[2]))
+        elseif i == 1 && j == nnodes(semi.solver) # top left corner
+            x_node = SVector(nextfloat(x_node[1]), prevfloat(x_node[2]))
+        elseif i == nnodes(semi.solver) && j == 1 # bottom right corner
+            x_node = SVector(prevfloat(x_node[1]), nextfloat(x_node[2]))
+        elseif i == nnodes(semi.solver) && j == nnodes(semi.solver) # top right corner
+            x_node = SVector(prevfloat(x_node[1]), prevfloat(x_node[2]))
+        elseif i == 1 # left boundary
+            x_node = SVector(nextfloat(x_node[1]), x_node[2])
+        elseif j == 1 # bottom boundary
+            x_node = SVector(x_node[1], nextfloat(x_node[2]))
+        elseif i == nnodes(semi.solver) # right boundary
+            x_node = SVector(prevfloat(x_node[1]), x_node[2])
+        elseif j == nnodes(semi.solver) # top boundary
+            x_node = SVector(x_node[1], prevfloat(x_node[2]))
+        end
+
+        u_node = initial_condition_dam_break(x_node, first(tspan), equations)
+        Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element)
+    end
+end
+
+###############################################################################
+# 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 multilayer 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,65 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with a periodic
+# bottom topography function (set in the initial conditions)
+
+equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.1,
+                                              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 = 6, 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)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+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