More robust LDG Switch for P4estMesh

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)

src/solvers/solvers_parabolic.jl

@@ -121,7 +121,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
 
@@ -157,12 +159,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()

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"),

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)