Add entropy stable LLF dissipation for the SWE (#132)

* add entropy stable LLF dissipation

* fix comments

* fix unit test

Author: patrickersing

src/equations/shallow_water_2d.jl

@@ -1130,7 +1130,9 @@ end
     return SVector(f1, f2, f3, 0)
 end
 
-# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography
+# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography.
+# For discontinuous bottom topography [`Trixi.DissipationLaxFriedrichsEntropyVariables`](@extref)
+# should be used instead as this version is not well-balanced.
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,
                                                               equations::ShallowWaterEquations2D)
@@ -1140,6 +1142,35 @@ end
     return SVector(diss[1], diss[2], diss[3], 0)
 end
 
+# Provably entropy stable and well-balanced local Lax-Friedrichs dissipation that avoids
+# spurious dissipation in the bottom topography.
+function (dissipation::DissipationLaxFriedrichsEntropyVariables)(u_ll, u_rr,
+                                                                 orientation_or_normal_direction,
+                                                                 equations::ShallowWaterEquations2D)
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations)
+    g = equations.gravity
+
+    # Compute averages
+    h_avg = 0.5f0 * (u_ll[1] + u_rr[1])
+    v1_avg = 0.5f0 * (u_ll[2] / u_ll[1] + u_rr[2] / u_rr[1])
+    v2_avg = 0.5f0 * (u_ll[3] / u_ll[1] + u_rr[3] / u_rr[1])
+
+    # Compute the jump in entropy variables
+    w_ll = cons2entropy(u_ll, equations)
+    w_rr = cons2entropy(u_rr, equations)
+    w_jump = SVector{3}(w_rr[1] - w_ll[1], w_rr[2] - w_ll[2], w_rr[3] - w_ll[3])
+
+    # Compute the H := du/dw
+    H = 1 / g * @SMatrix [1 v1_avg v2_avg;
+                  v1_avg g * h_avg+v1_avg^2 v1_avg*v2_avg;
+                  v2_avg v1_avg*v2_avg g * h_avg+v2_avg^2]
+
+    diss = SVector(-0.5f0 * λ * H * w_jump)
+
+    return SVector(diss..., 0)
+end
+
 # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography
 @inline function (numflux::FluxHLL)(u_ll, u_rr, orientation_or_normal_direction,
                                     equations::ShallowWaterEquations2D)

src/equations/shallow_water_multilayer_2d.jl

@@ -993,7 +993,8 @@ Use in combination with the generic numerical flux routine [`FluxHydrostaticReco
 end
 
 # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom
-# topography
+# topography. For discontinuous bottom topography [`Trixi.DissipationLaxFriedrichsEntropyVariables`](@extref)
+# should be used instead as this version is not well-balanced.
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,
                                                               equations::ShallowWaterMultiLayerEquations2D)
@@ -1003,6 +1004,37 @@ end
     return SVector(@views diss[1:(end - 1)]..., zero(eltype(u_ll)))
 end
 
+# Provably entropy stable and well-balanced local Lax-Friedrichs dissipation that avoids
+# spurious dissipation in the bottom topography.
+function (dissipation::DissipationLaxFriedrichsEntropyVariables)(u_ll, u_rr,
+                                                                 orientation_or_normal_direction,
+                                                                 equations::ShallowWaterMultiLayerEquations2D{4,
+                                                                                                              1,
+                                                                                                              T}) where {T}
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations)
+    g = equations.gravity
+
+    # Compute averages
+    h_avg = 0.5f0 * (u_ll[1] + u_rr[1])
+    v1_avg = 0.5f0 * (u_ll[2] / u_ll[1] + u_rr[2] / u_rr[1])
+    v2_avg = 0.5f0 * (u_ll[3] / u_ll[1] + u_rr[3] / u_rr[1])
+
+    # Compute the jump in entropy variables
+    w_ll = cons2entropy(u_ll, equations)
+    w_rr = cons2entropy(u_rr, equations)
+    w_jump = SVector{3}(w_rr[1] - w_ll[1], w_rr[2] - w_ll[2], w_rr[3] - w_ll[3])
+
+    # Compute the H := du/dw
+    H = 1 / g * @SMatrix [1 v1_avg v2_avg;
+                  v1_avg g * h_avg+v1_avg^2 v1_avg*v2_avg;
+                  v2_avg v1_avg*v2_avg g * h_avg+v2_avg^2]
+
+    diss = SVector(-0.5f0 * λ * H * w_jump)
+
+    return SVector(diss..., 0)
+end
+
 @inline function Trixi.max_abs_speeds(u, equations::ShallowWaterMultiLayerEquations2D)
     h = waterheight(u, equations)
     h_v1, h_v2 = momentum(u, equations)

test/test_tree_2d.jl

@@ -288,6 +288,30 @@ isdir(outdir) && rm(outdir, recursive = true)
         @test_allocations(Trixi.rhs!, semi, sol, 1000)
     end
 
+    @trixi_testset "elixir_shallowwater_source_terms.jl with es-llf dissipation" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_shallowwater_source_terms.jl"),
+                            l2=[
+                                0.0003733798777192089,
+                                0.02043287183047546,
+                                0.052383482460735806,
+                                6.274146767724067e-5
+                            ],
+                            linf=[
+                                0.0022222461533480953,
+                                0.06863676418052789,
+                                0.13355361762228668,
+                                0.0001819675955490041
+                            ],
+                            surface_flux=(FluxPlusDissipation(flux_fjordholm_etal,
+                                                              DissipationLaxFriedrichsEntropyVariables()),
+                                          flux_nonconservative_fjordholm_etal),
+                            tspan=(0.0, 0.25))
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    end
+
     @trixi_testset "elixir_shallowwater_conical_island.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_shallowwater_conical_island.jl"),

test/test_unit.jl

@@ -634,22 +634,6 @@ end
     end
 end
 
-@timed_testset "Consistency check for SWME dissipation terms" begin
-    # Test that for b=constant, the dissipation terms return the same output.
-    equations = ShallowWaterMomentEquations1D(gravity = 9.81, n_moments = 2)
-    equations_lin = ShallowWaterLinearizedMomentEquations1D(gravity = 9.81,
-                                                            n_moments = 2)
-    u_ll = SVector(1.0, 0.3, 0.15, 0.15, 0.1)
-    u_rr = SVector(1.5, 0.1, 0.25, 0.35, 0.1)
-
-    diss_lf = DissipationLocalLaxFriedrichs()
-    diss_ec = DissipationLaxFriedrichsEntropyVariables()
-
-    @test diss_lf(u_ll, u_rr, 1, equations) ≈ diss_ec(u_ll, u_rr, 1, equations)
-    @test diss_lf(u_ll, u_rr, 1, equations_lin) ≈
-          diss_ec(u_ll, u_rr, 1, equations_lin)
-end
-
 # Check consistency for conservative two-point fluxes (f(u, u) = f(u))
 @timed_testset "Consistency check for SWME fluxes" begin
     equations = ShallowWaterMomentEquations1D(gravity = 9.81, n_moments = 2)
@@ -669,6 +653,47 @@ end
     @test poly ≈ 1
     @test deriv ≈ 0
 end
+
+# Test that for b=constant, the dissipation terms return the same output.
+@testset "Consistency check for DissipationLaxFriedrichsEntropyVariables" begin
+    @timed_testset "ShallowWaterEquations2D" begin
+        equations = ShallowWaterEquations2D(gravity = 9.81)
+        u_ll = SVector(1.0, 0.3, 0.2, 0.2)
+        u_rr = SVector(1.5, 0.1, -0.1, 0.2)
+        @test DissipationLaxFriedrichsEntropyVariables()(u_ll, u_rr, 1, equations) ≈
+              DissipationLocalLaxFriedrichs()(u_ll, u_rr, 1, equations)
+    end
+
+    @timed_testset "ShalloWaterMultiLayerEquations2D (single layer)" begin
+        equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, rhos = (1.0))
+        u_ll = SVector(1.0, 0.3, 0.2, 0.1)
+        u_rr = SVector(1.5, 0.1, -0.1, 0.1)
+        @test DissipationLaxFriedrichsEntropyVariables()(u_ll, u_rr, 1, equations) ≈
+              DissipationLocalLaxFriedrichs()(u_ll, u_rr, 1, equations)
+
+        # For a single layer the dissipation should be consistent with the 2D SWE
+        u_ll = SVector(1.0, 0.3, 0.2, 0.1)
+        u_rr = SVector(1.5, 0.1, -0.1, 0.2)
+        equations_swe = ShallowWaterEquations2D(gravity = 9.81)
+        @test DissipationLaxFriedrichsEntropyVariables()(u_ll, u_rr, 1, equations) ≈
+              DissipationLaxFriedrichsEntropyVariables()(u_ll, u_rr, 1, equations_swe)
+    end
+
+    @timed_testset "ShallowWaterMomentEquations1D" begin
+        equations = ShallowWaterMomentEquations1D(gravity = 9.81, n_moments = 2)
+        equations_lin = ShallowWaterLinearizedMomentEquations1D(gravity = 9.81,
+                                                                n_moments = 2)
+        u_ll = SVector(1.0, 0.3, 0.15, 0.15, 0.1)
+        u_rr = SVector(1.5, 0.1, 0.25, 0.35, 0.1)
+
+        diss_lf = DissipationLocalLaxFriedrichs()
+        diss_ec = DissipationLaxFriedrichsEntropyVariables()
+
+        @test diss_lf(u_ll, u_rr, 1, equations) ≈ diss_ec(u_ll, u_rr, 1, equations)
+        @test diss_lf(u_ll, u_rr, 1, equations_lin) ≈
+              diss_ec(u_ll, u_rr, 1, equations_lin)
+    end
+end
 end # Unit tests
 
 end # module