diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl index 956c7ac..bd1b789 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_EC_correction.jl +++ b/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. diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl index 21934c5..bbb80e7 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_EC_projection.jl +++ b/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 diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl b/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl index 0749c83..7de23ce 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_advection.jl +++ b/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, diff --git a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl b/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl index ba1c3bf..71d597c 100644 --- a/examples/elixir_shallowwater_cubed_sphere_shell_standard.jl +++ b/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 diff --git a/src/equations/shallow_water_3d.jl b/src/equations/shallow_water_3d.jl index 8340868..897f467 100644 --- a/src/equations/shallow_water_3d.jl +++ b/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)