Add new non-conservative flux for standard and two-layer shallow water equations

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl

@@ -13,8 +13,8 @@ initial_condition = initial_condition_convergence_test
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl

@@ -35,8 +35,8 @@ initial_condition = initial_condition_dam_break
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

@@ -33,8 +33,8 @@ initial_condition = initial_condition_fjordholm_well_balanced
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg=3, surface_flux=(flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal),
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_es_ersing_etal, flux_nonconservative_ersing_etal),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl

@@ -12,8 +12,8 @@ initial_condition = initial_condition_convergence_test
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 

examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

@@ -28,8 +28,8 @@ initial_condition = initial_condition_well_balanced
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal)
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal)
 solver = DGSEM(polydeg=3, surface_flux=surface_flux,
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 

examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl

@@ -14,8 +14,8 @@ initial_condition = initial_condition_convergence_test
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
 solver = DGSEM(polydeg=6, surface_flux=surface_flux,
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 

examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl

@@ -40,8 +40,8 @@ boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_b
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-surface_flux= (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+surface_flux= (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
 solver = DGSEM(polydeg=6, surface_flux=surface_flux,
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 

examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

@@ -30,8 +30,8 @@ initial_condition = initial_condition_well_balanced
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
-surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal)
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal)
 solver = DGSEM(polydeg=6, surface_flux=surface_flux,
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 

src/Trixi.jl

@@ -159,8 +159,9 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner,
        flux_nonconservative_powell,
        flux_kennedy_gruber, flux_shima_etal, flux_ec,
-       flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal,
+       flux_fjordholm_etal, flux_nonconservative_fjordholm_etal,
        flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,
+       flux_es_ersing_etal, flux_nonconservative_ersing_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
        FluxLaxFriedrichs, max_abs_speed_naive,

src/equations/shallow_water_1d.jl

@@ -307,6 +307,44 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
                    z)
 end
 
+"""
+    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
+                                     equations::ShallowWaterEquations1D)
+
+!!! warning "Experimental code"
+    This numerical flux is experimental and may change in any future release.
+
+Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
+that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).
+
+This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy 
+conservation and well-balancedness in both the volume and surface when combined with 
+[`flux_wintermeyer_etal`](@ref).
+
+For further details see:
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
+"""
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
+                                                  equations::ShallowWaterEquations1D)
+    # Pull the necessary left and right state information
+    h_ll = waterheight(u_ll, equations)
+    b_rr = u_rr[3]
+    b_ll = u_ll[3]
+
+    # Calculate jump
+    b_jump = b_rr - b_ll
+
+    z = zero(eltype(u_ll))
+
+    # Bottom gradient nonconservative term: (0, g h b_x, 0)
+    f = SVector(z, equations.gravity * h_ll * b_jump, z)
+
+    return f
+end
+
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation,
                         equations::ShallowWaterEquations1D)

src/equations/shallow_water_2d.jl

@@ -516,6 +516,74 @@ end
     return SVector(f1, f2, f3, f4)
 end
 
+"""
+    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
+                                     equations::ShallowWaterEquations2D)
+    flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                     normal_direction_ll::AbstractVector,
+                                     normal_direction_average::AbstractVector,
+                                     equations::ShallowWaterEquations2D)
+
+!!! warning "Experimental code"
+    This numerical flux is experimental and may change in any future release.
+
+Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
+that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).
+
+On curvilinear meshes, this nonconservative flux depends on both the
+contravariant vector (normal direction) at the current node and the averaged
+one. This is different from numerical fluxes used to discretize conservative
+terms.
+
+This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy 
+conservation and well-balancedness in both the volume and surface when combined with 
+[`flux_wintermeyer_etal`](@ref).
+
+For further details see:
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
+"""
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
+                                                  equations::ShallowWaterEquations2D)
+    # Pull the necessary left and right state information
+    h_ll = waterheight(u_ll, equations)
+    b_rr = u_rr[4]
+    b_ll = u_ll[4]
+
+    # Calculate jump
+    b_jump = b_rr - b_ll
+
+    z = zero(eltype(u_ll))
+    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
+    if orientation == 1
+        f = SVector(z, equations.gravity * h_ll * b_jump, z, z)
+    else # orientation == 2
+        f = SVector(z, z, equations.gravity * h_ll * b_jump, z)
+    end
+    return f
+end
+
+@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
+                                                  normal_direction_ll::AbstractVector,
+                                                  normal_direction_average::AbstractVector,
+                                                  equations::ShallowWaterEquations2D)
+    # Pull the necessary left and right state information
+    h_ll = waterheight(u_ll, equations)
+    b_rr = u_rr[4]
+    b_ll = u_ll[4]
+
+    # Calculate jump
+    b_jump = b_rr - b_ll
+    # Note this routine only uses the `normal_direction_average` and the average of the
+    # bottom topography to get a quadratic split form DG gradient on curved elements
+    return SVector(zero(eltype(u_ll)),
+                   normal_direction_average[1] * equations.gravity * h_ll * b_jump,
+                   normal_direction_average[2] * equations.gravity * h_ll * b_jump,
+                   zero(eltype(u_ll)))
+end
+
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction,
                         equations::ShallowWaterEquations2D)