Added non-conservative terms to Cartesian SWE 3D

examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl

@@ -8,10 +8,11 @@ using TrixiAtmo
 
 equations = ShallowWaterEquations3D(gravity_constant = 9.81)
 
-# Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux
+# Create DG solver with polynomial degree = 3 and Fjordholm et al.'s flux as surface flux
 polydeg = 3
-volume_flux = flux_fjordholm_etal #flux_wintermeyer_etal
-surface_flux = flux_fjordholm_etal #flux_wintermeyer_etal  #flux_lax_friedrichs
+volume_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
+surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal)
+
 solver = DGSEM(polydeg = polydeg,
                surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
@@ -43,7 +44,7 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     # Compute the velocity of a rotating solid body
     # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system)
     v0 = 1.0 # Velocity at the "equator"
-    alpha = 0.0 #pi / 4
+    alpha = pi / 4
     v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
     vlat = -v0 * sin(phi) * sin(alpha)
 
@@ -52,7 +53,10 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long
     v3 = sin(colat) * vlat
 
-    return prim2cons(SVector(rho, v1, v2, v3, 0), equations)
+    # Non-constant topography
+    b = 0.1 * cos(10 * lat)
+
+    return prim2cons(SVector(rho, v1, v2, v3, b), equations)
 end
 
 # Source term function to apply the entropy correction term and the Lagrange multiplier to the semi-discretization.

examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl

@@ -12,8 +12,9 @@ equations = ShallowWaterEquations3D(gravity_constant = 9.81)
 
 # Create DG solver with polynomial degree = 3 and Wintemeyer et al.'s flux as surface flux
 polydeg = 3
-volume_flux = flux_wintermeyer_etal # flux_fjordholm_etal 
-surface_flux = flux_wintermeyer_etal # flux_fjordholm_etal #flux_lax_friedrichs
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+
 solver = DGSEM(polydeg = polydeg,
                surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
@@ -45,7 +46,7 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     # Compute the velocity of a rotating solid body
     # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system)
     v0 = 1.0 # Velocity at the "equator"
-    alpha = 0.0 #pi / 4
+    alpha = pi / 4
     v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha))
     vlat = -v0 * sin(phi) * sin(alpha)
 
@@ -54,7 +55,10 @@ function initial_condition_advection_sphere(x, t, equations::ShallowWaterEquatio
     v2 = -cos(colat) * sin(phi) * vlat + cos(phi) * v_long
     v3 = sin(colat) * vlat
 
-    return prim2cons(SVector(rho, v1, v2, v3, 0), equations)
+    # Non-constant topography
+    b = 0.1 * cos(10 * lat)
+
+    return prim2cons(SVector(rho, v1, v2, v3, b), equations)
 end
 
 initial_condition = initial_condition_advection_sphere

examples/elixir_shallowwater_cubed_sphere_shell_advection.jl

@@ -8,13 +8,18 @@ using TrixiAtmo
 
 cells_per_dimension = 5
 initial_condition = initial_condition_gaussian
+polydeg = 3
 
 # We use the ShallowWaterEquations3D equations structure but modify the rhs! function to
 # convert it to a variable-coefficient advection equation
 equations = ShallowWaterEquations3D(gravity_constant = 0.0)
 
-# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
-solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+volume_flux = (flux_central, flux_nonconservative_fjordholm_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal)
+
+solver = DGSEM(polydeg = polydeg,
+               surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 # Source term function to transform the Euler equations into a linear advection equation with variable advection velocity
 function source_terms_convert_to_linear_advection(u, du, x, t,

examples/elixir_shallowwater_cubed_sphere_shell_standard.jl

@@ -10,7 +10,12 @@ equations = ShallowWaterEquations3D(gravity_constant = 9.81)
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs as surface flux
 polydeg = 3
-solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs)
+volume_flux = (flux_central, flux_nonconservative_wintermeyer_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal)
+
+solver = DGSEM(polydeg = polydeg,
+               surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 # Initial condition for a Gaussian density profile with constant pressure
 # and the velocity of a rotating solid body

src/equations/shallow_water_3d.jl

@@ -61,7 +61,7 @@ function ShallowWaterEquations3D(; gravity_constant, H0 = zero(gravity_constant)
     ShallowWaterEquations3D(gravity_constant, H0)
 end
 
-Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = False() # Deactivate non-conservative terms for the moment...
+Trixi.have_nonconservative_terms(::ShallowWaterEquations3D) = True() # Deactivate non-conservative terms for the moment...
 Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquations3D) = ("h", "h_v1", "h_v2",
                                                                   "h_v3", "b")
 # Note, we use the total water height, H = h + b, as the first primitive variable for easier
@@ -159,6 +159,88 @@ Further details are available in Theorem 1 of the paper:
     return SVector(f1, f2, f3, f4, zero(eltype(u_ll)))
 end
 
+"""
+    flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
+                                          equations::ShallowWaterEquations2D)
+    flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
+                                          normal_direction::AbstractVector,
+                                          equations::ShallowWaterEquations2D)
+
+Non-symmetric two-point volume flux discretizing the nonconservative (source) term
+that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).
+
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
+
+Further details are available in the papers:
+- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
+  An entropy stable nodal discontinuous Galerkin method for the two dimensional
+  shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
+  [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
+- 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 Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
+                                                             normal_direction::AbstractVector,
+                                                             equations::ShallowWaterEquations3D)
+    # Pull the necessary left and right state information
+    h_ll = waterheight(u_ll, equations)
+    b_jump = u_rr[5] - u_ll[5]
+
+    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
+    return SVector(0,
+                   normal_direction[1] * equations.gravity * h_ll * b_jump,
+                   normal_direction[2] * equations.gravity * h_ll * b_jump,
+                   normal_direction[3] * equations.gravity * h_ll * b_jump,
+                   0)
+end
+
+"""
+    flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
+                                        equations::ShallowWaterEquations2D)
+    flux_nonconservative_fjordholm_etal(u_ll, u_rr,
+                                        normal_direction::AbstractVector,
+                                        equations::ShallowWaterEquations2D)
+
+Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
+that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).
+
+This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
+conservation and well-balancedness.
+
+Further details for the original finite volume formulation are available in
+- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011)
+  Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
+  [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
+and for curvilinear 2D case in the paper:
+- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
+  An entropy stable nodal discontinuous Galerkin method for the two dimensional
+  shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
+  [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
+"""
+@inline function Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr,
+                                                           normal_direction::AbstractVector,
+                                                           equations::ShallowWaterEquations3D)
+    # Pull the necessary left and right state information
+    h_ll, _, _, _, b_ll = u_ll
+    h_rr, _, _, _, b_rr = u_rr
+
+    h_average = 0.5f0 * (h_ll + h_rr)
+    b_jump = b_rr - b_ll
+
+    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
+    f2 = normal_direction[1] * equations.gravity * h_average * b_jump
+    f3 = normal_direction[2] * equations.gravity * h_average * b_jump
+    f4 = normal_direction[3] * equations.gravity * h_average * b_jump
+
+    # First and last equations do not have a nonconservative flux
+    f1 = f5 = 0
+
+    return SVector(f1, f2, f3, f4, f5)
+end
+
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation,
                         equations::ShallowWaterEquations1D)