diff --git a/src/solvers/dgmulti/sbp.jl b/src/solvers/dgmulti/sbp.jl
index 232555e18b..76444f6489 100644
--- a/src/solvers/dgmulti/sbp.jl
+++ b/src/solvers/dgmulti/sbp.jl
@@ -221,7 +221,10 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh,
                 for id in nzrange(A_base, i)
                     j = rows[id]
                     u_j = u[j]
-                    A_ij = vals[id]
+
+                    # we use the negative of A_ij since A is skew-symmetric, 
+                    # and we are accessing the transpose of A. 
+                    A_ij = -vals[id]
                     AF_ij = 2 * A_ij *
                             volume_flux(u_i, u_j, normal_direction, equations)
                     du_i = du_i + AF_ij
diff --git a/test/test_threaded.jl b/test/test_threaded.jl
index a3d52c1923..83544aa007 100644
--- a/test/test_threaded.jl
+++ b/test/test_threaded.jl
@@ -442,7 +442,7 @@ end
         end
     end
 
-    @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin
+    @trixi_testset "elixir_euler_fdsbp_periodic.jl (2D)" begin
         @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d",
                                      "elixir_euler_fdsbp_periodic.jl"),
                             l2=[
@@ -467,6 +467,33 @@ end
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 5000
         end
     end
+
+    @trixi_testset "elixir_euler_fdsbp_periodic.jl (3D)" begin
+        @test_trixi_include(joinpath(examples_dir(),
+                                     "dgmulti_3d/elixir_euler_fdsbp_periodic.jl"),
+                            l2=[
+                                7.561896970325353e-5,
+                                6.884047859361093e-5,
+                                6.884047859363204e-5,
+                                6.884047859361148e-5,
+                                0.000201107274617457
+                            ],
+                            linf=[
+                                0.0001337520020225913,
+                                0.00011571467799287305,
+                                0.0001157146779990903,
+                                0.00011571467799376123,
+                                0.0003446082308800058
+                            ])
+        # 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