Fix convergence problems for flux_nonconservative_ersing_etal_local_jump (#131)

patrickersing · amrueda · web-flow · commit cb3deba6243c · 2026-03-09T14:42:22.000+01:00
* use local normal direction for flux_ersing_local_jump

* add elixirs to test convergence and EC on curvilinear meshes

* apply formatter

* fix testing

* format examples

* make ec test case reproducible

* add note for flux_nonconservative_ersing_local_jump

* make use of local normal direction explicit

---------

Co-authored-by: Andrés Rueda-Ramírez &lt;aruedara@uni-koeln.de&gt;
diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl
@@ -0,0 +1,134 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+using TrixiShallowWater
+using Symbolics
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with a single layer to test
+# convergence. The initial condition and source term are created using the 
+# method of manufactured solutions (MMS) with the help of Symbolics.jl.
+
+equations = ShallowWaterMultiLayerEquations2D(gravity = 1.1, rhos = (1.0))
+
+### Create manufactured solution for method of manufactured solutions (MMS)
+
+# Symbolic Variables
+@variables x_sym[1:2], t_sym, g
+
+# Define Differentials
+Dt, Dx, Dy = Differential(t_sym), Differential(x_sym[1]), Differential(x_sym[2])
+
+## Initial condition
+###############################################################################
+# Primitive Variables
+ω = pi * 1.0
+H = 4.0 + 0.2 * cos(ω * x_sym[1] + t_sym) + 0.2 * cos(ω * x_sym[2] + t_sym)
+v1 = 0.5
+v2 = 0.5
+b = 1.0 + 0.2 * cos(ω * x_sym[1]) + 0.2 * cos(ω * x_sym[2])
+h = H - b
+
+init = [H, v1, v2, b]
+
+## PDE
+###############################################################################
+eqs = [
+    Dt(h) + Dx(h * v1) + Dy(h * v2),
+    Dt(h * v1) + Dx(h * v1^2 + 0.5 * g * h^2) + Dy(h * v1 * v2) + g * h * Dx(b),
+    Dt(h * v2) + Dx(h * v1 * v2) + Dy(h * v2^2 + 0.5 * g * h^2) + g * h * Dy(b),
+    0
+]
+
+## Create the functions for the manufactured solution
+########################################################################################
+# Expand derivatives
+du_exprs = expand_derivatives.(eqs)
+
+# Build functions
+du_funcs = build_function.(du_exprs, Ref(x_sym), t_sym, g, expression = Val(false))
+
+init_funcs = build_function.(init, Ref(x_sym), t_sym, expression = Val(false))
+
+function initial_condition_convergence_mms(x,
+                                           t,
+                                           equations::ShallowWaterMultiLayerEquations2D)
+    prim = SVector{4, Float64}([f(x, t) for f in init_funcs]...)
+    return prim2cons(prim, equations)
+end
+
+function source_terms_convergence_mms(u,
+                                      x,
+                                      t,
+                                      equations::ShallowWaterMultiLayerEquations2D)
+    g = equations.gravity
+    return SVector{4, Float64}([f(x, t, g) for f in du_funcs]...)
+end
+
+initial_condition = initial_condition_convergence_mms
+
+###############################################################################
+# Get the DG approximation space
+
+polydeg = 3
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
+
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;)
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the P4estMesh and setup a periodic mesh on the domain [-1,1]^2 with an affine type mapping to
+# obtain a warped curvilinear mesh.
+function mapping_twist(xi, eta)
+    y = eta + 0.1 * sin(pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.1 * sin(pi * eta) * cos(0.5 * pi * xi)
+    return SVector(x, y)
+end
+
+# Create warped P4estMesh with 8 x 8 elements
+trees_per_dimension = (2, 2)
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 initial_refinement_level = 2,
+                 periodicity = true,
+                 mapping = mapping_twist)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_mms,
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+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)
+
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(),)
+
+sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  ode_default_options()...,
+                  callback = callbacks, adaptive = false);
diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_ec_sc_subcell.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_ec_sc_subcell.jl
@@ -0,0 +1,110 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with a single layer and a bottom
+# topography function for a blast wave test with discontinuous initial conditions to test entropy
+# conservation with the subcell limiting volume integral on a curvilinear mesh.
+
+equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, rhos = (1.0))
+
+function initial_condition_weak_blast_wave(x, t,
+                                           equations::ShallowWaterMultiLayerEquations2D)
+    # Set up polar coordinates
+    RealT = eltype(x)
+    inicenter = SVector(convert(RealT, 0.0), convert(RealT, 0.0))
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+    phi = atan(y_norm, x_norm)
+    sin_phi, cos_phi = sincos(phi)
+
+    # Calculate primitive variables
+    H = r > 0.2f0 ? 3.8f0 : 4.0f0
+    v1 = r > 0.2f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi
+    v2 = r > 0.2f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi
+    v1 = 0.0
+    v2 = 0.0
+
+    # Setup a continuous bottom topography
+    r = 0.4
+    b = (((x[1])^2 + (x[2])^2) < r^2 ?
+         0.5 * (cos(1 / r * pi * sqrt((x[1])^2 + (x[2])^2)) + 1) : 0.0)
+
+    return prim2cons(SVector(H, v1, v2, b), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+###############################################################################
+# Get the DG approximation space
+
+polydeg = 3
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal_local_jump)
+
+basis = LobattoLegendreBasis(polydeg)
+# For reproducibility set the limiting coefficients to pure DG,
+# see https://github.com/trixi-framework/Trixi.jl/pull/2007
+limiter_idp = SubcellLimiterIDP(equations, basis;)
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the P4estMesh and setup a periodic mesh on the domain [-1,1]^2 with an affine type mapping to
+# obtain a warped curvilinear mesh.
+function mapping_twist(xi, eta)
+    y = eta + 0.1 * sin(pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.1 * sin(pi * eta) * cos(0.5 * pi * xi)
+    return SVector(x, y)
+end
+
+# Create P4estMesh with 8 x 8 elements
+trees_per_dimension = (2, 2)
+mesh = P4estMesh(trees_per_dimension, polydeg = 2,
+                 initial_refinement_level = 2,
+                 periodicity = true,
+                 mapping = mapping_twist)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 0.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_polydeg = polydeg,
+                                     extra_analysis_errors = (:conservation_error,))
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.1,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     extra_node_variables = (:limiting_coefficient,))
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+
+# Setup the stage callbacks
+stage_callbacks = (SubcellLimiterIDPCorrection(),)
+
+# run the simulation
+sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  ode_default_options()...,
+                  callback = callbacks, adaptive = false);
diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl
@@ -44,7 +44,7 @@ ode = semidiscretize(semi, tspan)
 ###############################################################################
 # Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.
 
-# alternative version of the initial conditinon used to setup a truly discontinuous
+# alternative version of the initial condition used to setup a truly discontinuous
 # bottom topography function and initial condition for this academic testcase of entropy conservation.
 # The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff
 # In contrast to the usual signature of initial conditions, this one get passed the
diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl
@@ -624,6 +624,12 @@ or jump portion of the non-conservative flux based on the type of the
 nonconservative_type argument, employing multiple dispatch. They are used to
 compute the subcell fluxes in [`Trixi.VolumeIntegralSubcellLimiting`](@extref).
 
+!!! note
+    While `flux_nonconservative_ersing_etal_local_jump` and [`flux_nonconservative_ersing_etal`](@ref)
+    are equivalent at interfaces, the volume formulation for [`Trixi.VolumeIntegralSubcellLimiting`](@extref)
+    differs as `flux_nonconservative_ersing_etal_local_jump` applies the local normal direction 
+    instead of the averaged one.
+
 ## References
 - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
   Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.
@@ -653,7 +659,7 @@ end
                                            nonconservative_type::Trixi.NonConservativeLocal,
                                            nonconservative_term::Integer)
     flux_nonconservative_ersing_local_jump(u_ll, u_rr,
-                                           normal_direction::AbstractVector,
+                                           normal_direction_ll::AbstractVector,
                                            equations::ShallowWaterMultiLayerEquations2D,
                                            nonconservative_type::Trixi.NonConservativeLocal,
                                            nonconservative_term::Integer)
@@ -695,7 +701,7 @@ Local part of the nonconservative term needed for the calculation of the non-con
 end
 
 @inline function flux_nonconservative_ersing_etal_local_jump(u_ll,
-                                                             normal_direction::AbstractVector,
+                                                             normal_direction_ll::AbstractVector,
                                                              equations::ShallowWaterMultiLayerEquations2D,
                                                              nonconservative_type::Trixi.NonConservativeLocal,
                                                              nonconservative_term::Integer)
@@ -714,7 +720,9 @@ end
         f_h = zero(real(equations))
         f_hv = g * h_ll[i]
 
-        setlayer!(f, f_h, f_hv, f_hv, i, equations)
+        setlayer!(f, f_h, f_hv * normal_direction_ll[1], f_hv * normal_direction_ll[2],
+                  i,
+                  equations)
     end
 
     return SVector(f)
@@ -812,8 +820,8 @@ end
                 f_hv += h_jump[j]
             end
         end
-        setlayer!(f, f_h, f_hv * normal_direction[1],
-                  f_hv * normal_direction[2], i, equations)
+        setlayer!(f, f_h, f_hv,
+                  f_hv, i, equations)
     end
 
     return SVector(f)
diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
@@ -165,6 +165,25 @@ end # SWE
         @test_allocations(Trixi.rhs!, semi, sol, 15000)
     end
 
+    @trixi_testset "elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_convergence_sc_subcell_curved.jl"),
+                            l2=[
+                                0.00013275624153302434,
+                                0.0001586600356395913,
+                                0.000158660035639501,
+                                2.9583315272612922e-5
+                            ],
+                            linf=[
+                                0.0007544272991792944,
+                                0.0007250877164874936,
+                                0.00072508771648927,
+                                9.31948456051046e-5
+                            ])
+        # Allocation testing is disabled as the symbolic source term computation is known to cause
+        # allocations.
+    end
+
     @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl"),
@@ -319,6 +338,30 @@ end # SWE
         # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
         @test_allocations(Trixi.rhs!, semi, sol, 15000)
     end
+
+    @trixi_testset "elixir_shallowwater_multilayer_ec_sc_subcell.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_multilayer_ec_sc_subcell.jl"),
+                            l2=[
+                                0.04338581073333642,
+                                0.12068863355590975,
+                                0.12068863355590971,
+                                7.906244739074657e-18
+                            ],
+                            linf=[
+                                0.2550462084344076,
+                                0.42004261172048024,
+                                0.4200426117204863,
+                                1.1102230246251565e-16
+                            ])
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        # Larger values for allowed allocations due to usage of custom
+        # integrator which are not *recorded* for the methods from
+        # OrdinaryDiffEq.jl
+        # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
+        @test_allocations(Trixi.rhs!, semi, sol, 15000)
+    end
 end # MLSWE
 end # P4estMesh2D
 
