diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl index fc76b4f034b..e6a01849852 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/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)) ############################################################################### diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl index b2e6a81401b..03b93754d0f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/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)) ############################################################################### diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl index 7236f1697d0..098e3aaf601 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/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)) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index f7c8ab3a249..790916e4467 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/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)) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl index 1495e6d8568..264c26390fe 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/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)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index 2cab68b1cb5..0b86095663a 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/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)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl index 9d70e9287cf..4ad5f7e3201 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/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)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl index 35b027c3a81..6a727df2502 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/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)) diff --git a/src/Trixi.jl b/src/Trixi.jl index 97d518d5b78..a224dec884e 100644 --- a/src/Trixi.jl +++ b/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 diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 32782d5478c..25ce0fa79fe 100644 --- a/src/equations/shallow_water_1d.jl +++ b/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) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index a81fddeed49..e75c92a27d0 100644 --- a/src/equations/shallow_water_2d.jl +++ b/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) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 4b64481cca3..42ff393593e 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -87,15 +87,15 @@ end have_nonconservative_terms(::ShallowWaterTwoLayerEquations1D) = True() function varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations1D) - ("h_upper", "h_v1_upper", - "h_lower", "h_v1_lower", "b") + ("h_upper", "h_v_upper", + "h_lower", "h_v_lower", "b") end # Note, we use the total water height, H_lower = h_upper + h_lower + b, and first layer total height # H_upper = h_upper + b as the first primitive variable for easier visualization and setting initial # conditions function varnames(::typeof(cons2prim), ::ShallowWaterTwoLayerEquations1D) - ("H_upper", "v1_upper", - "H_lower", "v1_lower", "b") + ("H_upper", "v_upper", + "H_lower", "v_lower", "b") end # Set initial conditions at physical location `x` for time `t` @@ -113,11 +113,11 @@ function initial_condition_convergence_test(x, t, H_lower = 2.0 + 0.1 * sin(ω * x[1] + t) H_upper = 4.0 + 0.1 * cos(ω * x[1] + t) - v1_lower = 1.0 - v1_upper = 0.9 + v_lower = 1.0 + v_upper = 0.9 b = 1.0 + 0.1 * cos(2.0 * ω * x[1]) - return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) + return prim2cons(SVector(H_upper, v_upper, H_lower, v_lower, b), equations) end """ @@ -196,165 +196,77 @@ end # Note, the bottom topography has no flux @inline function flux(u, orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u + h_upper, h_v_upper, h_lower, h_v_lower, _ = u # Calculate velocities - v1_upper, v1_lower = velocity(u, equations) + v_upper, v_lower = velocity(u, equations) # Calculate pressure - p1 = 0.5 * equations.gravity * h_upper^2 - p2 = 0.5 * equations.gravity * h_lower^2 + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 - f1 = h_v1_upper - f2 = h_v1_upper * v1_upper + p1 - f3 = h_v2_lower - f4 = h_v2_lower * v1_lower + p2 + f1 = h_v_upper + f2 = h_v_upper * v_upper + p_upper + f3 = h_v_lower + f4 = h_v_lower * v_lower + p_lower return SVector(f1, f2, f3, f4, zero(eltype(u))) end """ - flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. -Non-symmetric two-point volume flux discretizing the nonconservative (source) term -that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version -to account for the additional source term compared to the standard SWE described in the paper. +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations1D`](@ref) and an +additional term that couples the momentum of both layers. -Further details are available 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) +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_wintermeyer_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) # Pull the necessary left and right state information h_upper_ll, h_lower_ll = waterheight(u_ll, equations) h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[5] + b_ll = u_ll[5] - z = zero(eltype(u_ll)) - - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, - # 0, g*h_lower*(b+r*h_upper)_x, 0) - f = SVector(z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), - z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), - z) - return f -end - -""" - flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both -layers [`ShallowWaterTwoLayerEquations2D`](@ref). - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Pull the necessary left and right state information - h_upper_ll, _, h_lower_ll, _, b_ll = u_ll - h_upper_rr, _, h_lower_rr, _, b_rr = u_rr - - # Create average and jump values - h_upper_average = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_average = 0.5 * (h_lower_ll + h_lower_rr) - h_upper_jump = h_upper_rr - h_upper_ll - h_lower_jump = h_lower_rr - h_lower_ll - b_jump = b_rr - b_ll - - # Assign variables for constants for better readability - g = equations.gravity + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) z = zero(eltype(u_ll)) # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, # 0, g*h_lower*(b+r*h_upper)_x, 0) f = SVector(z, - g * h_upper_ll * (b_ll + h_lower_ll) + - g * h_upper_average * (b_jump + h_lower_jump), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + - g * h_lower_average * (b_jump + - equations.r * h_upper_jump), + equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), z) return f end -""" - flux_fjordholm_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations1D) - -Total energy conservative (mathematical entropy for shallow water equations). When the bottom -topography is nonzero this should only be used as a surface flux otherwise the scheme will not be -well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). - -Details are available in Eq. (4.1) in the paper: -- 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 the application to two layers is shown in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Unpack left and right state - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - v1_ll, v2_ll = velocity(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - v1_rr, v2_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p1_avg = 0.25 * equations.gravity * (h_upper_ll^2 + h_upper_rr^2) - p2_avg = 0.25 * equations.gravity * (h_lower_ll^2 + h_lower_rr^2) - - # Calculate fluxes - f1 = h_upper_avg * v1_avg - f2 = f1 * v1_avg + p1_avg - f3 = h_lower_avg * v2_avg - f4 = f3 * v2_avg + p2_avg - - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) -end - """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. -When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. -The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +When the bottom topography is nonzero this scheme will be well-balanced when used with the +nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -367,51 +279,48 @@ Further details are available in Theorem 1 of the paper: orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_lower_rr, h_v2_lower_rr, _ = u_rr + h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll + h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr # Get the velocities on either side - v1_ll, v2_ll = velocity(u_ll, equations) - v1_rr, v2_rr = velocity(u_rr, equations) + v_upper_ll, v_lower_ll = velocity(u_ll, equations) + v_upper_rr, v_lower_rr = velocity(u_rr, equations) # Average each factor of products in flux - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - p1_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p2_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + v_upper_avg = 0.5 * (v_upper_ll + v_upper_rr) + v_lower_avg = 0.5 * (v_lower_ll + v_lower_rr) + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr # Calculate fluxes - f1 = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) - f2 = f1 * v1_avg + p1_avg - f3 = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) - f4 = f3 * v2_avg + p2_avg + f1 = 0.5 * (h_v_upper_ll + h_v_upper_rr) + f2 = f1 * v_upper_avg + p_upper_avg + f3 = 0.5 * (h_v_lower_ll + h_v_lower_rr) + f4 = f3 * v_lower_avg + p_lower_avg return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) end """ - flux_es_fjordholm_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations1D) - -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy -conservative [`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump -of entropy variables. - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. + flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations1D) +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. + +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_es_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) +@inline function flux_es_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) # Compute entropy conservative flux but without the bottom topography - f_ec = flux_fjordholm_etal(u_ll, u_rr, - orientation, - equations) + f_ec = flux_wintermeyer_etal(u_ll, u_rr, + orientation, + equations) # Get maximum signal velocity λ = max_abs_speed_naive(u_ll, u_rr, orientation, equations) @@ -474,12 +383,12 @@ end orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_lower_rr, h_v2_lower_rr, _ = u_rr + h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll + h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr # Get the averaged velocity - v_m_ll = (h_v1_upper_ll + h_v2_lower_ll) / (h_upper_ll + h_lower_ll) - v_m_rr = (h_v1_upper_rr + h_v2_lower_rr) / (h_upper_rr + h_lower_rr) + v_m_ll = (h_v_upper_ll + h_v_lower_ll) / (h_upper_ll + h_lower_ll) + v_m_rr = (h_v_upper_rr + h_v_lower_rr) / (h_upper_rr + h_lower_rr) # Calculate the wave celerity on the left and right h_upper_ll, h_lower_ll = waterheight(u_ll, equations) @@ -503,10 +412,10 @@ end # Absolute speed of the barotropic mode @inline function max_abs_speeds(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u + h_upper, h_v_upper, h_lower, h_v_lower, _ = u # Calculate averaged velocity of both layers - v_m = (h_v1_upper + h_v2_lower) / (h_upper + h_lower) + v_m = (h_v_upper + h_v_lower) / (h_upper + h_lower) c = sqrt(equations.gravity * (h_upper + h_lower)) return (abs(v_m) + c) @@ -514,11 +423,11 @@ end # Helper function to extract the velocity vector from the conservative variables @inline function velocity(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u + h_upper, h_v_upper, h_lower, h_v_lower, _ = u - v1_upper = h_v1_upper / h_upper - v1_lower = h_v2_lower / h_lower - return SVector(v1_upper, v1_lower) + v_upper = h_v_upper / h_upper + v_lower = h_v_lower / h_lower + return SVector(v_upper, v_lower) end # Convert conservative variables to primitive @@ -527,8 +436,8 @@ end H_lower = h_lower + b H_upper = h_lower + h_upper + b - v1_upper, v1_lower = velocity(u, equations) - return SVector(H_upper, v1_upper, H_lower, v1_lower, b) + v_upper, v_lower = velocity(u, equations) + return SVector(H_upper, v_upper, H_lower, v_lower, b) end # Convert conservative variables to entropy variables @@ -536,26 +445,26 @@ end # bottom topography values for convenience @inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D) h_upper, _, h_lower, _, b = u - v1_upper, v1_lower = velocity(u, equations) - - w1 = equations.rho_upper * - (equations.gravity * (h_upper + h_lower + b) - 0.5 * v1_upper^2) - w2 = equations.rho_upper * v1_upper - w3 = equations.rho_lower * - (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v1_lower^2) - w4 = equations.rho_lower * v1_lower + v_upper, v_lower = velocity(u, equations) + + w1 = (equations.rho_upper * + (equations.gravity * (h_upper + h_lower + b) - 0.5 * v_upper^2)) + w2 = equations.rho_upper * v_upper + w3 = (equations.rho_lower * + (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v_lower^2)) + w4 = equations.rho_lower * v_lower return SVector(w1, w2, w3, w4, b) end # Convert primitive to conservative variables @inline function prim2cons(prim, equations::ShallowWaterTwoLayerEquations1D) - H_upper, v1_upper, H_lower, v1_lower, b = prim + H_upper, v_upper, H_lower, v_lower, b = prim h_lower = H_lower - b h_upper = H_upper - h_lower - b - h_v1_upper = h_upper * v1_upper - h_v2_lower = h_lower * v1_lower - return SVector(h_upper, h_v1_upper, h_lower, h_v2_lower, b) + h_v_upper = h_upper * v_upper + h_v_lower = h_lower * v_lower + return SVector(h_upper, h_v_upper, h_lower, h_v_lower, b) end @inline function waterheight(u, equations::ShallowWaterTwoLayerEquations1D) @@ -569,23 +478,23 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, b = cons + h_upper, h_v_upper, h_lower, h_v_lower, b = cons # Set new variables for better readability g = equations.gravity rho_upper = equations.rho_upper rho_lower = equations.rho_lower - e = (0.5 * rho_upper * (h_v1_upper^2 / h_upper + g * h_upper^2) + - 0.5 * rho_lower * (h_v2_lower^2 / h_lower + g * h_lower^2) + + e = (0.5 * rho_upper * (h_v_upper^2 / h_upper + g * h_upper^2) + + 0.5 * rho_lower * (h_v_lower^2 / h_lower + g * h_lower^2) + g * rho_lower * h_lower * b + g * rho_upper * h_upper * (h_lower + b)) return e end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v1_upper, h_lower, h_v2_lower, _ = u - return 0.5 * equations.rho_upper * h_v1_upper^2 / h_upper + - 0.5 * equations.rho_lower * h_v2_lower^2 / h_lower + h_upper, h_v_upper, h_lower, h_v_lower, _ = u + return (0.5 * equations.rho_upper * h_v_upper^2 / h_upper + + 0.5 * equations.rho_lower * h_v_lower^2 / h_lower) end # Calculate potential energy for a conservative state `cons` diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 87249e91948..a31d881f2ef 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -271,24 +271,24 @@ end v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) # Calculate pressure - p1 = 0.5 * equations.gravity * h_upper^2 - p2 = 0.5 * equations.gravity * h_lower^2 + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 # Calculate fluxes depending on orientation if orientation == 1 f1 = h_v1_upper - f2 = h_v1_upper * v1_upper + p1 + f2 = h_v1_upper * v1_upper + p_upper f3 = h_v1_upper * v2_upper f4 = h_v1_lower - f5 = h_v1_lower * v1_lower + p2 + f5 = h_v1_lower * v1_lower + p_lower f6 = h_v1_lower * v2_lower else f1 = h_v2_upper f2 = h_v2_upper * v1_upper - f3 = h_v2_upper * v2_upper + p1 + f3 = h_v2_upper * v2_upper + p_upper f4 = h_v2_lower f5 = h_v2_lower * v1_lower - f6 = h_v2_lower * v2_lower + p2 + f6 = h_v2_lower * v2_lower + p_lower end return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) end @@ -305,44 +305,57 @@ end h_v_upper_normal = h_upper * v_normal_upper h_v_lower_normal = h_lower * v_normal_lower - p1 = 0.5 * equations.gravity * h_upper^2 - p2 = 0.5 * equations.gravity * h_lower^2 + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 f1 = h_v_upper_normal - f2 = h_v_upper_normal * v1_upper + p1 * normal_direction[1] - f3 = h_v_upper_normal * v2_upper + p1 * normal_direction[2] + f2 = h_v_upper_normal * v1_upper + p_upper * normal_direction[1] + f3 = h_v_upper_normal * v2_upper + p_upper * normal_direction[2] f4 = h_v_lower_normal - f5 = h_v_lower_normal * v1_lower + p2 * normal_direction[1] - f6 = h_v_lower_normal * v2_lower + p2 * normal_direction[2] + f5 = h_v_lower_normal * v1_lower + p_lower * normal_direction[1] + f6 = h_v_lower_normal * v2_lower + p_lower * normal_direction[2] return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) end """ - flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) !!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. + This numerical flux is experimental and may change in any future release. -Non-symmetric two-point volume flux discretizing the nonconservative (source) term +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. This is a slightly modified version -to account for the additional source term compared to the standard SWE described in the paper. +additional term that couples the momentum of both layers. -Further details are available 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) +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_wintermeyer_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) # Pull the necessary left and right state information h_upper_ll, h_lower_ll = waterheight(u_ll, equations) h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[7] + b_ll = u_ll[7] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) z = zero(eltype(u_ll)) @@ -351,268 +364,64 @@ Further details are available in the paper: # g*h_lower*(b + r*h_upper)_y, 0) if orientation == 1 f = SVector(z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), z, z) else # orientation == 2 f = SVector(z, z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), z) end return f end -@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) # Pull the necessary left and right state information h_upper_ll, h_lower_ll = waterheight(u_ll, equations) h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[7] + b_ll = u_ll[7] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + 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_upper_ll * - (b_rr + h_lower_rr), + (b_jump + h_lower_jump), normal_direction_average[2] * equations.gravity * h_upper_ll * - (b_rr + h_lower_rr), + (b_jump + h_lower_jump), zero(eltype(u_ll)), normal_direction_average[1] * equations.gravity * h_lower_ll * - (b_rr + - equations.r * h_upper_rr), + (b_jump + equations.r * h_upper_jump), normal_direction_average[2] * equations.gravity * h_lower_ll * - (b_rr + - equations.r * h_upper_rr), + (b_jump + equations.r * h_upper_jump), zero(eltype(u_ll))) end -""" - flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains -the gradients of the bottom topography and an additional term that couples the momentum of both -layers [`ShallowWaterTwoLayerEquations2D`](@ref). - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Pull the necessary left and right state information - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, b_ll = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, b_rr = u_rr - - # Create average and jump values - h_upper_average = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_average = 0.5 * (h_lower_ll + h_lower_rr) - h_upper_jump = h_upper_rr - h_upper_ll - h_lower_jump = h_lower_rr - h_lower_ll - b_jump = b_rr - b_ll - - # Assign variables for constants for better readability - g = equations.gravity - - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, g*h_upper*(b+h_lower)_y, 0, - # g*h_lower*(b+r*h_upper)_x, g*h_lower*(b+r*h_upper)_x, 0) - - # Includes two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid - # cross-averaging across a discontinuous bottom topography - # (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry - z = zero(eltype(u_ll)) - if orientation == 1 - f = SVector(z, - g * h_upper_ll * (b_ll + h_lower_ll) + - g * h_upper_average * (b_jump + h_lower_jump), - z, z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + - g * h_lower_average * (b_jump + - equations.r * h_upper_jump), - z, z) - else # orientation == 2 - f = SVector(z, z, - g * h_upper_ll * (b_ll + h_lower_ll) + - g * h_upper_average * (b_jump + h_lower_jump), - z, z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + - g * h_lower_average * (b_jump + - equations.r * h_upper_jump), - z) - end - - return f -end - -@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Pull the necessary left and right state information - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, b_ll = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, b_rr = u_rr - - # Create average and jump values - h_upper_average = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_average = 0.5 * (h_lower_ll + h_lower_rr) - h_upper_jump = h_upper_rr - h_upper_ll - h_lower_jump = h_lower_rr - h_lower_ll - b_jump = b_rr - b_ll - - # Comes in two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average` - # but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography - f2 = normal_direction_average[1] * equations.gravity * h_upper_ll * - (b_ll + h_lower_ll) - f3 = normal_direction_average[2] * equations.gravity * h_upper_ll * - (b_ll + h_lower_ll) - f5 = normal_direction_average[1] * equations.gravity * h_lower_ll * - (b_ll + equations.r * h_upper_ll) - f6 = normal_direction_average[2] * equations.gravity * h_lower_ll * - (b_ll + equations.r * h_upper_ll) - # (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump` - # to handle discontinuous bathymetry - f2 += normal_direction_ll[1] * equations.gravity * h_upper_average * - (b_jump + h_lower_jump) - f3 += normal_direction_ll[2] * equations.gravity * h_upper_average * - (b_jump + h_lower_jump) - f5 += normal_direction_ll[1] * equations.gravity * h_lower_average * - (b_jump + - equations.r * h_upper_jump) - f6 += normal_direction_ll[2] * equations.gravity * h_lower_average * - (b_jump + - equations.r * h_upper_jump) - - # Continuity equations do not have a nonconservative flux - f1 = f4 = zero(eltype(u_ll)) - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - -""" - flux_fjordholm_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations2D) - -Total energy conservative (mathematical entropy for two-layer shallow water equations). When the -bottom topography is nonzero this should only be used as a surface flux otherwise the scheme will -not be well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). - -Details are available in Eq. (4.1) in the paper: -- 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 the application to two layers is shown in the paper: -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. -""" -@inline function flux_fjordholm_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) - v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) - v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) - v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) - v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.25 * equations.gravity * (h_upper_ll^2 + h_upper_rr^2) - p2_avg = 0.25 * equations.gravity * (h_lower_ll^2 + h_lower_rr^2) - - # Calculate fluxes depending on orientation - if orientation == 1 - f1 = h_upper_avg * v1_upper_avg - f2 = f1 * v1_upper_avg + p1_avg - f3 = f1 * v2_upper_avg - f4 = h_lower_avg * v1_lower_avg - f5 = f4 * v1_lower_avg + p2_avg - f6 = f4 * v2_lower_avg - else - f1 = h_upper_avg * v2_upper_avg - f2 = f1 * v1_upper_avg - f3 = f1 * v2_upper_avg + p1_avg - f4 = h_lower_avg * v2_lower_avg - f5 = f4 * v1_lower_avg - f6 = f4 * v2_lower_avg + p2_avg - end - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - -@inline function flux_fjordholm_etal(u_ll, u_rr, - normal_direction::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - # Compute velocity in normal direction - v_upper_dot_n_ll = v1_upper_ll * normal_direction[1] + - v2_upper_ll * normal_direction[2] - v_upper_dot_n_rr = v1_upper_rr * normal_direction[1] + - v2_upper_rr * normal_direction[2] - v_lower_dot_n_ll = v1_lower_ll * normal_direction[1] + - v2_lower_ll * normal_direction[2] - v_lower_dot_n_rr = v1_lower_rr * normal_direction[1] + - v2_lower_rr * normal_direction[2] - - # Average each factor of products in flux - h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) - h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) - v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) - v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) - v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) - v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.25 * equations.gravity * (h_upper_ll^2 + h_upper_rr^2) - p2_avg = 0.25 * equations.gravity * (h_lower_ll^2 + h_lower_rr^2) - v_upper_dot_n_avg = 0.5 * (v_upper_dot_n_ll + v_upper_dot_n_rr) - v_lower_dot_n_avg = 0.5 * (v_lower_dot_n_ll + v_lower_dot_n_rr) - - # Calculate fluxes depending on normal_direction - f1 = h_upper_avg * v_upper_dot_n_avg - f2 = f1 * v1_upper_avg + p1_avg * normal_direction[1] - f3 = f1 * v2_upper_avg + p1_avg * normal_direction[2] - f4 = h_lower_avg * v_lower_dot_n_avg - f5 = f4 * v1_lower_avg + p2_avg * normal_direction[1] - f6 = f4 * v2_lower_avg + p2_avg * normal_direction[2] - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations2D) + flux_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. -When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. -The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow +When the bottom topography is nonzero this scheme will be well-balanced when used with the +nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the +two-layer shallow water equations the flux that is described in the paper for the normal shallow water equations is used within each layer. Further details are available in Theorem 1 of the paper: @@ -637,24 +446,24 @@ Further details are available in Theorem 1 of the paper: v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p2_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr # Calculate fluxes depending on orientation if orientation == 1 f1 = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) - f2 = f1 * v1_upper_avg + p1_avg + f2 = f1 * v1_upper_avg + p_upper_avg f3 = f1 * v2_upper_avg f4 = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) - f5 = f4 * v1_lower_avg + p2_avg + f5 = f4 * v1_lower_avg + p_lower_avg f6 = f4 * v2_lower_avg else f1 = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) f2 = f1 * v1_upper_avg - f3 = f1 * v2_upper_avg + p1_avg + f3 = f1 * v2_upper_avg + p_upper_avg f4 = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) f5 = f4 * v1_lower_avg - f6 = f4 * v2_lower_avg + p2_avg + f6 = f4 * v2_lower_avg + p_lower_avg end return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) @@ -676,8 +485,8 @@ end v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p1_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p2_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr h_v1_upper_avg = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) h_v2_upper_avg = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) h_v1_lower_avg = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) @@ -685,38 +494,36 @@ end # Calculate fluxes depending on normal_direction f1 = h_v1_upper_avg * normal_direction[1] + h_v2_upper_avg * normal_direction[2] - f2 = f1 * v1_upper_avg + p1_avg * normal_direction[1] - f3 = f1 * v2_upper_avg + p1_avg * normal_direction[2] + f2 = f1 * v1_upper_avg + p_upper_avg * normal_direction[1] + f3 = f1 * v2_upper_avg + p_upper_avg * normal_direction[2] f4 = h_v1_lower_avg * normal_direction[1] + h_v2_lower_avg * normal_direction[2] - f5 = f4 * v1_lower_avg + p2_avg * normal_direction[1] - f6 = f4 * v2_lower_avg + p2_avg * normal_direction[2] + f5 = f4 * v1_lower_avg + p_lower_avg * normal_direction[1] + f6 = f4 * v2_lower_avg + p_lower_avg * normal_direction[2] return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) end """ - flux_es_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations1D) - -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative -[`flux_fjordholm_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy -variables. - -Further details are available in the paper: -- Ulrik Skre Fjordholm (2012) -Energy conservative and stable schemes for the two-layer shallow water equations. -[DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -It should be noted that the equations are ordered differently and the -designation of the upper and lower layer has been changed which leads to a slightly different -formulation. + flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) + +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. + +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_es_fjordholm_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) +@inline function flux_es_ersing_etal(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) # Compute entropy conservative flux but without the bottom topography - f_ec = flux_fjordholm_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations) + f_ec = flux_wintermeyer_etal(u_ll, u_rr, + orientation_or_normal_direction, + equations) # Get maximum signal velocity λ = max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations) @@ -926,12 +733,12 @@ end rho_lower = equations.rho_lower v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - w1 = rho_upper * (equations.gravity * (h_upper + h_lower + b) + - -0.5 * (v1_upper^2 + v2_upper^2)) + w1 = (rho_upper * (equations.gravity * (h_upper + h_lower + b) + + -0.5 * (v1_upper^2 + v2_upper^2))) w2 = rho_upper * v1_upper w3 = rho_upper * v2_upper - w4 = rho_lower * (equations.gravity * (equations.r * h_upper + h_lower + b) + - -0.5 * (v1_lower^2 + v2_lower^2)) + w4 = (rho_lower * (equations.gravity * (equations.r * h_upper + h_lower + b) + + -0.5 * (v1_lower^2 + v2_lower^2))) w5 = rho_lower * v1_lower w6 = rho_lower * v2_lower return SVector(w1, w2, w3, w4, w5, w6, b) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 658f178c941..3e2db7ca67f 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -96,6 +96,29 @@ end end end +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254838, + 1.6657566141935285e-14, + 0.10416666834254838, + ], + linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), @@ -167,6 +190,33 @@ end end end +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.005774284062933275, + 0.017408601639513584, + 4.43649172561843e-5, + ], + linf=[ + 0.01639116193303547, + 0.05102877460799604, + 9.098379777450205e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index a504f4f93a6..180fb3ec3b3 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -10,95 +10,65 @@ include("test_trixi.jl") EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Shallow Water Two layer" begin -#! format: noindent - -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0050681532925156945, 0.002089013899370176, - 0.005105544300292713, 0.002526442122643306, - 0.0004744186597732706], - linf=[0.022256679217306008, 0.005421833004652266, - 0.02233993939574197, 0.008765261497422516, - 0.0008992474511784199], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.005012009872109003, 0.002091035326731071, + 0.005049271397924551, + 0.0024633066562966574, 0.0004744186597732739], + linf=[0.0213772149343594, 0.005385752427290447, + 0.02175023787351349, + 0.008212004668840978, 0.0008992474511784199], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl with flux_es_fjordholm_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0027681377074701345, 0.0018007543202559165, - 0.0028036917433720576, - 0.0013980358596935737, 0.0004744186597732706], - linf=[0.005699303919826093, 0.006432952918256296, - 0.0058507082844360125, 0.002717615543961216, - 0.0008992474511784199], - surface_flux=(flux_es_fjordholm_etal, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[8.949288784402005e-16, 4.0636427176237915e-17, + 0.001002881985401548, + 2.133351105037203e-16, 0.0010028819854016578], + linf=[2.6229018956769323e-15, 1.878451903240623e-16, + 0.005119880996670156, + 8.003199803957679e-16, 0.005119880996670666], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[8.949288784402005e-16, 4.0636427176237915e-17, - 0.001002881985401548, - 2.133351105037203e-16, 0.0010028819854016578], - linf=[2.6229018956769323e-15, 1.878451903240623e-16, - 0.005119880996670156, - 8.003199803957679e-16, 0.005119880996670666], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_dam_break.jl"), + l2=[0.1000774903431289, 0.5670692949571057, 0.08764242501014498, + 0.45412307886094555, 0.013638618139749523], + linf=[0.586718937495144, 2.1215606128311584, 0.5185911311186155, + 1.820382495072612, 0.5], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end end -@trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_dam_break.jl"), - l2=[0.10010269243463918, 0.5668733957648654, - 0.08759617327649398, - 0.4538443183566172, 0.013638618139749523], - linf=[ - 0.5854202777756559, - 2.1278930820498934, - 0.5193686074348809, - 1.8071213168086229, - 0.5, - ], - surface_flux=(flux_lax_friedrichs, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end -end - end # module diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index d280e380192..58db7c5f35f 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -116,6 +116,35 @@ end end end +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987146, + 1.0323158914614244e-14, + 1.0276096319430528e-14, + 0.9130579602987147, + ], + linf=[ + 2.11306203761566, + 4.063916419044386e-14, + 3.694484044448245e-14, + 2.1130620376156584, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), @@ -219,6 +248,35 @@ end end end +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.002471853426064005, + 0.05619168608950033, + 0.11844727575152562, + 6.274146767730281e-5, + ], + linf=[ + 0.014332922987500218, + 0.2141204806174546, + 0.5392313755637872, + 0.0001819675955490041, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_conical_island.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_conical_island.jl"), l2=[ diff --git a/test/test_tree_2d_shallowwater_twolayer.jl b/test/test_tree_2d_shallowwater_twolayer.jl index 5959e7ed882..802bf4e021c 100644 --- a/test/test_tree_2d_shallowwater_twolayer.jl +++ b/test/test_tree_2d_shallowwater_twolayer.jl @@ -10,105 +10,79 @@ include("test_trixi.jl") EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @testset "Two-Layer Shallow Water" begin -#! format: noindent - -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0004040147445601598, 0.005466848793475609, - 0.006149138398472166, 0.0002908599437447256, - 0.003011817461911792, 0.0026806180089700674, - 8.873630921431545e-6], - linf=[0.002822006686981293, 0.014859895905040332, - 0.017590546190827894, 0.0016323702636176218, - 0.009361402900653015, 0.008411036357379165, - 3.361991620143279e-5], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.0004016779699408397, 0.005466339651545468, + 0.006148841330156112, + 0.0002882339012602492, 0.0030120142442780313, + 0.002680752838455618, + 8.873630921431545e-6], + linf=[0.002788654460984752, 0.01484602033450666, + 0.017572229756493973, + 0.0016010835493927011, 0.009369847995372549, + 0.008407961775489636, + 3.361991620143279e-5], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_convergence.jl with flux_es_fjordholm_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.00024709443131137236, 0.0019215286339769443, - 0.0023833298173254447, - 0.00021258247976270914, 0.0011299428031136195, - 0.0009191313765262401, - 8.873630921431545e-6], - linf=[0.0016099763244645793, 0.007659242165565017, - 0.009123320235427057, - 0.0013496983982568267, 0.0035573687287770994, - 0.00296823235874899, - 3.361991620143279e-5], - surface_flux=(flux_es_fjordholm_etal, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[3.2935164267930016e-16, 4.6800825611195103e-17, + 4.843057532147818e-17, + 0.0030769233188015013, 1.4809161150389857e-16, + 1.509071695038043e-16, + 0.0030769233188014935], + linf=[2.248201624865942e-15, 2.346382070278936e-16, + 2.208565017494899e-16, + 0.026474051138910493, 9.237568031609006e-16, + 7.520758026187046e-16, + 0.026474051138910267], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -end -@trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[3.2935164267930016e-16, 4.6800825611195103e-17, - 4.843057532147818e-17, - 0.0030769233188015013, 1.4809161150389857e-16, - 1.509071695038043e-16, - 0.0030769233188014935], - linf=[2.248201624865942e-15, 2.346382070278936e-16, - 2.208565017494899e-16, - 0.026474051138910493, 9.237568031609006e-16, - 7.520758026187046e-16, - 0.026474051138910267], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced with flux_lax_friedrichs.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[2.0525741072929735e-16, 6.000589392730905e-17, + 6.102759428478984e-17, + 0.0030769233188014905, 1.8421386173122792e-16, + 1.8473184927121752e-16, + 0.0030769233188014935], + linf=[7.355227538141662e-16, 2.960836949170518e-16, + 4.2726562436938764e-16, + 0.02647405113891016, 1.038795478061861e-15, + 1.0401789378532516e-15, + 0.026474051138910267], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end end -@trixi_testset "elixir_shallowwater_twolayer_well_balanced with flux_lax_friedrichs.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[2.0525741072929735e-16, 6.000589392730905e-17, - 6.102759428478984e-17, - 0.0030769233188014905, 1.8421386173122792e-16, - 1.8473184927121752e-16, - 0.0030769233188014935], - linf=[7.355227538141662e-16, 2.960836949170518e-16, - 4.2726562436938764e-16, - 0.02647405113891016, 1.038795478061861e-15, - 1.0401789378532516e-15, - 0.026474051138910267], - surface_flux=(flux_lax_friedrichs, - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end -end - end # module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 26483931cf3..d4416ac5b6a 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -349,6 +349,35 @@ end end end +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839083, + 2.590643638636187e-12, + 1.0945471514840143e-12, + 1.2164292510839079, + ], + linf=[ + 1.5138512282315792, + 5.0276441977281156e-11, + 1.9816934589292803e-11, + 1.513851228231574, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -402,6 +431,35 @@ end end end +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011196687776346434, + 0.044562672453443995, + 0.014306265289763618, + 5.089218476759981e-6, + ], + linf=[ + 0.007825021762002393, + 0.348550815397918, + 0.1115517935018282, + 2.6407324614341476e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -535,15 +593,15 @@ end @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0007953969898161991, 0.00882074628714633, - 0.0024322572528892934, - 0.0007597425017400447, 0.004501238950166439, - 0.0015784803573661104, + l2=[0.0007935561625451243, 0.008825315509943844, + 0.002429969315645897, + 0.0007580145888686304, 0.004495741879625235, + 0.0015758146898767814, 6.849532064729749e-6], - linf=[0.00592559068081977, 0.08072451118697077, - 0.0344854497419107, 0.005892196680485795, - 0.04262651217675306, 0.014006223513881366, - 2.5829318284764646e-5], + linf=[0.0059205195991136605, 0.08072126590166251, + 0.03463806075399023, + 0.005884818649227186, 0.042658506561995546, + 0.014125956138838602, 2.5829318284764646e-5], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -582,17 +640,17 @@ end @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_dam_break.jl"), - l2=[0.012471300561905669, 0.012363413819726868, - 0.0009541478004413331, - 0.09120260327331643, 0.015269590815749993, - 0.0012064657396853422, - 0.09991983966647647], - linf=[0.04497814714937959, 0.03286959000796511, - 0.010746094385294369, - 0.11138723974511211, 0.03640850605444494, - 0.014368386516056392, 0.10000000000000003], + l2=[0.012447632879122346, 0.012361250464676683, + 0.0009551519536340908, + 0.09119400061322577, 0.015276216721920347, + 0.0012126995108983853, 0.09991983966647647], + linf=[0.044305765721807444, 0.03279620980615845, + 0.010754320388190101, + 0.111309922939555, 0.03663360204931427, + 0.014332822306649284, + 0.10000000000000003], surface_flux=(flux_lax_friedrichs, - flux_nonconservative_fjordholm_etal), + flux_nonconservative_ersing_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)