Add new non-conservative flux for standard and two-layer shallow wate…

…r equations (#1508)

* Add flux_nonconservative_ersing_etal

* Change formulation of flux_nonconservative_ersing

* remove old fluxes (fjordholm, wintermeyer)

* Update tests for new flux_ersing_etal

* Add flux_nonconservative_ersing for standard SWE

* Edit comments, change flux in elixir

* fix well-balanced test

* Update src/Trixi.jl

Change order of exporting statements

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>

* fix indentation

* fix indentation II

* Update src/equations/shallow_water_2d.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>

* fix energy total function

* Apply formatter

* add docstring references

* remove flux_nonconservative_fjordholm_etal

* apply formatter

* Format examples

---------

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl

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

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl

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

examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

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

examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl

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

examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl

@@ -31,8 +31,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

@@ -15,8 +15,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

@@ -41,8 +41,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

@@ -33,8 +33,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

@@ -164,8 +164,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_nonconservative_powell_local_symmetric,
        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,
        flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
 # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file

src/equations/shallow_water_1d.jl

@@ -380,6 +380,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

@@ -702,6 +702,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)