Two-layer shallow water equations (trixi-framework#1297)

* Add new equation system for two-layer swe in 1d * Added slip_wall_bc, testing and elixirs. Updated MMS to depend on time. Equation file restructured to match 1d-SWE. * Changed the system to match literature: - Switched equation order to (h1,h1v1,h2,h2v2) - upper layer = 1; lower layer = 2 * Add new equation system for two-layer swe in 2d * Add flux functions for unstructured & slip_wall bc * Add entropy stable flux and tests Edit code format and comments * update test reference for 2l-swe * Change assignment of nonconservative terms Correct some code formatting and comments * add new field "r" to the equation struct * Fix bug in max_abs_speeds formula, add description * Add error msg for density ratios > 1 * Change equation name to match file names * Update the elixirs and tests * fix tests, add flux description, optimize H matrix * add experimental warning, indent and update latex * additional test to improve coverage * Fix some typos and spacing * fix spacing, replace non-ascii chars * change name of flux_es to flux_es_fjordholm_etal * Change variables names to upper / lower * Change array type to MArray * change position of flux_es_fjordholm_etal * remove allocations in flux_es_fjordholm * fix typos in docstring * add name to contributor list * fix consistency issues and comment elixir --------- Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Andrew Winters <[email protected]>
Xinyu-Li-1997 · Apr 19, 2023 · 5e28148 · 5e28148
1 parent cc2cda4
commit 5e28148
Show file tree

Hide file tree

Showing 18 changed files with 2,355 additions and 2 deletions.
diff --git a/AUTHORS.md b/AUTHORS.md
@@ -27,6 +27,7 @@ are listed in alphabetical order:
 * Jesse Chan
 * Lars Christmann
 * Christof Czernik
+* Patrick Ersing
 * Erik Faulhaber
 * Gregor Gassner
 * Lucas Gemein

diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl
@@ -0,0 +1,60 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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_wintermeyer_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_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
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl
@@ -0,0 +1,108 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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)
+
+###############################################################################
+# Workaround to set a discontinuous bottom topography and initial condition.
+
+# This test case uses a special work around to setup a truly discontinuous bottom topography 
+# function and initial condition for this academic testcase of entropy conservation. First, a 
+# dummy initial condition is introduced to create the semidiscretization. Then the initial condition 
+# is reset with the true discontinuous values from initial_condition_dam_break. Note, that this
+# initial condition only works for TreeMesh1D with `initial_refinement_level=5`.
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_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, callbacks etc.
+
+tspan = (0.0,0.4)
+ode = semidiscretize(semi, tspan)
+
+# Initial conditions dam break test case
+function initial_condition_dam_break(x, t, element_id, equations::ShallowWaterTwoLayerEquations1D)
+  v1_upper = 0.0
+  v1_lower = 0.0
+
+  # Set the discontinuity
+  if element_id <= 16
+    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
+
+
+# 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)
+    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates,
+      equations, semi.solver, i, element)
+    u_node = initial_condition_dam_break(x_node, first(tspan), element, equations)
+    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
+  end
+end
+
+
+
+
+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
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl
@@ -0,0 +1,82 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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_wintermeyer_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_es_fjordholm_etal, flux_nonconservative_fjordholm_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
diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl
@@ -0,0 +1,59 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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_wintermeyer_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_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