diff --git a/NEWS.md b/NEWS.md
index 711bad308be..e30364b1dd3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,6 +5,16 @@ Trixi.jl follows the interpretation of
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes when updating to v0.16 from v0.15.x
+
+- The implementation of the local DG (`ViscousFormulationLocalDG`) `solver_parabolic` has been changed for the `P4estMesh`.
+In particular, instead of computing the `ldg_switch` as the dot product of the normal direction with ones,
+i.e., summing up the normal components, the `ldg_switch` is now selected as 
+the sign of the maximum (in absolute value sense) normal direction component,
+which corresponds to the dominant direction of the interface normal.
+This might change results slightly for some meshes where the sum of the normal might be close to zero,
+thus introducing some spurious switch assignments ([#2871]).
+
 ## Changes in the v0.15 lifecycle
 
 #### Added
diff --git a/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl b/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl
new file mode 100644
index 00000000000..5bf67a09000
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_advection_diffusion_rotated.jl
@@ -0,0 +1,60 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection-diffusion equation
+
+diffusivity() = 1.0e-2
+advection_velocity = (-1.0, 1.0)
+equations = LinearScalarAdvectionEquation2D(advection_velocity)
+equations_parabolic = LaplaceDiffusion2D(diffusivity(), equations)
+
+function initial_condition_gauss_damped(x, t, equations)
+    damping_factor = 1 + 4 * diffusivity() * t
+    return SVector(exp(-(x[1]^2 + x[2]^2) / damping_factor) / damping_factor)
+end
+initial_condition = initial_condition_gauss_damped
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+# This maps the domain [-1, 1]^2 to a 45-degree rotated increased square
+square_size() = 5.0
+function mapping(xi, eta)
+    x = square_size() * xi
+    y = square_size() * eta
+    return SVector((x - y) / sqrt(2), (x + y) / sqrt(2))
+end
+
+trees_per_dimension = (23, 23)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 mapping = mapping, periodicity = true)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             solver_parabolic = ViscousFormulationLocalDG(),
+                                             boundary_conditions = (boundary_condition_periodic,
+                                                                    boundary_condition_periodic))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+n_passes = 2
+tspan = (0.0, n_passes * square_size() * sqrt(2))
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1.0e-6
+sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            ode_default_options()..., callback = callbacks)
diff --git a/src/solvers/solvers_parabolic.jl b/src/solvers/solvers_parabolic.jl
index 77ee2c64c9b..9f502547d36 100644
--- a/src/solvers/solvers_parabolic.jl
+++ b/src/solvers/solvers_parabolic.jl
@@ -100,11 +100,12 @@ f_{\text{gradient}} = u_{L}
 ```
 on the Cartesian [`TreeMesh`](@ref).
 
-For the [`P4estMesh`](@ref), the `normal_direction` is used to compute the LDG "switch" ``\sigma`` for the upwinding/downwinding.
-This is realized by taking the sign of the dot product of the normal and positive-coordinate direction vector:
+For the [`P4estMesh`](@ref), the `normal_direction` is used to compute the LDG "switch" ``\sigma`` for the upwinding.
+This is realized by selecting the sign of the maximum (in absolute value sense) normal direction component,
+which corresponds to the "dominant" direction of the interface normal.
 ```math
-\sigma = \text{sign}(\vec{n} \cdot \vec{1})
-f = \frac{1}{2}\Big(f(u_{L}) + f(u_{R}) - \sigma \big[f(u_{R}) - f(u_{L})\big]\Big)
+i = \text{argmax} \{ \begin{pmatrix} \vert n_1 \vert \\ \vert n_2 \vert \\ \dots \end{pmatrix} \}
+\sigma = \text{sign} (n_i)
 ```
 """
 function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh`
@@ -121,7 +122,9 @@ end
 function flux_parabolic(u_ll, u_rr, normal_direction,
                         ::Gradient, equations_parabolic,
                         parabolic_scheme::ViscousFormulationLocalDG)
-    ldg_switch = sign(sum(normal_direction)) # equivalent to sign(dot(normal_direction, ones))
+    # Use "Upwind in dominant direction" for LDG switch
+    abs_max_dir = argmax(abs.(normal_direction))
+    ldg_switch = sign(normal_direction[abs_max_dir])
     return 0.5f0 * (u_ll + u_rr - ldg_switch * (u_rr - u_ll))
 end
 
@@ -145,11 +148,12 @@ f_{\text{divergence}} = u_{R}
 ```
 on the Cartesian [`TreeMesh`](@ref).
 
-For the [`P4estMesh`](@ref), the `normal_direction` is used to compute the LDG "switch" ``\sigma`` for the upwinding/downwinding.
-This is realized by taking the sign of the dot product of the normal and positive-coordinate direction vector:
+For the [`P4estMesh`](@ref), the `normal_direction` is used to compute the LDG "switch" ``\sigma`` for the downwinding.
+This is realized by selecting the sign of the maximum (in absolute value sense) normal direction component,
+which corresponds to the "dominant" direction of the interface normal.
 ```math
-\sigma = \text{sign}(\vec{n} \cdot \vec{1})
-f = \frac{1}{2}\Big(f(u_{L}) + f(u_{R}) + \sigma \big[f(u_{R}) - f(u_{L})\big]\Big)
+i = \text{argmax} \{ \begin{pmatrix} \vert n_1 \vert \\ \vert n_2 \vert \\ \dots \end{pmatrix} \}
+\sigma = -\text{sign} (n_i)
 ```
 """
 function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh`
@@ -157,12 +161,14 @@ function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh`
                         parabolic_scheme::ViscousFormulationLocalDG)
     return u_rr # Use the downwind value for the divergence interface flux
 end
-# Version or `P4estMesh`
+# Version for `P4estMesh`
 function flux_parabolic(u_ll, u_rr, normal_direction,
                         ::Divergence, equations_parabolic,
                         parabolic_scheme::ViscousFormulationLocalDG)
-    ldg_switch = sign(sum(normal_direction)) # equivalent to sign(dot(normal_direction, ones))
-    return 0.5f0 * (u_ll + u_rr + ldg_switch * (u_rr - u_ll))
+    # Use "Downwind in dominant direction" for LDG switch
+    abs_max_dir = argmax(abs.(normal_direction))
+    ldg_switch = -sign(normal_direction[abs_max_dir])
+    return 0.5f0 * (u_ll + u_rr - ldg_switch * (u_rr - u_ll))
 end
 
 default_parabolic_solver() = ViscousFormulationBassiRebay1()
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 19953579506..cafb1812161 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -750,6 +750,16 @@ end
     @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
 end
 
+@trixi_testset "P4estMesh2D: elixir_advection_diffusion_rotated.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
+                                 "elixir_advection_diffusion_rotated.jl"),
+                        l2=[4.8533724384822306e-5], linf=[0.0006284491001110615])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test_allocations(Trixi.rhs_parabolic!, semi, sol, 1000)
+end
+
 @trixi_testset "P4estMesh2D: elixir_advection_diffusion_periodic_curved.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
                                  "elixir_advection_diffusion_periodic_curved.jl"),
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index f8aa7e85afa..4e86de3cfa7 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -527,7 +527,7 @@ end
                                  "elixir_advection_diffusion_nonperiodic.jl"),
                         solver_parabolic=ViscousFormulationLocalDG(),
                         cfl_diffusive=0.07,
-                        l2=[0.0041854757843498725], linf=[0.05166356737492643])
+                        l2=[0.004185076476662267], linf=[0.05166349548111486])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     @test_allocations(Trixi.rhs!, semi, sol, 1000)