diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl
new file mode 100644
index 00000000000..b8dbad50680
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl
@@ -0,0 +1,82 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25)
+
+# An initial condition with a bottom topography and a perturbation in the waterheight to test 
+# boundary_condition_slip_wall
+function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D)
+    # Set the background values
+    H = equations.H0
+    v1 = 0.0
+    v2 = 0.0
+
+    # Bottom topography
+    b = 1.5 * exp(-0.5 * ((x[1])^2 + (x[2])^2))
+    # Waterheight perturbation
+    H = H + 0.5 * exp(-10.0 * ((x[1])^2 + (x[2])^2))
+
+    return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_perturbation
+
+boundary_condition = boundary_condition_slip_wall
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_lax_friedrichs, 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
+
+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_condition)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.25)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     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
diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl
index e75c92a27d0..6728d7d5553 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -236,8 +236,12 @@ Should be used together with [`TreeMesh`](@ref).
         u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4])
     end
 
-    # compute and return the flux using `boundary_condition_slip_wall` routine above
-    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+    # Calculate boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
+    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
+    end
 
     return flux
 end
diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl
index 58db7c5f35f..1f3dfbf5267 100644
--- a/test/test_tree_2d_shallowwater.jl
+++ b/test/test_tree_2d_shallowwater.jl
@@ -327,6 +327,31 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_shallowwater_wall.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"),
+                        l2=[
+                            0.13517233723296504,
+                            0.20010876311162215,
+                            0.20010876311162223,
+                            2.719538414346464e-7,
+                        ],
+                        linf=[
+                            0.5303607982988336,
+                            0.5080989745682338,
+                            0.5080989745682352,
+                            1.1301675764130437e-6,
+                        ],
+                        tspan=(0.0, 0.25))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 end # module