From 7de4ce1b68c4d5aabb278a4de62236b3d65562c5 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Sat, 20 May 2023 19:20:57 +0200 Subject: [PATCH 01/17] Add flux_nonconservative_ersing_etal --- ...lixir_shallowwater_twolayer_convergence.jl | 4 +- .../elixir_shallowwater_twolayer_dam_break.jl | 4 +- ...xir_shallowwater_twolayer_well_balanced.jl | 4 +- ...lixir_shallowwater_twolayer_convergence.jl | 4 +- ...xir_shallowwater_twolayer_well_balanced.jl | 4 +- ...lixir_shallowwater_twolayer_convergence.jl | 4 +- .../elixir_shallowwater_twolayer_dam_break.jl | 4 +- ...xir_shallowwater_twolayer_well_balanced.jl | 4 +- src/Trixi.jl | 3 +- src/equations/shallow_water_two_layer_1d.jl | 46 ++++++++++- src/equations/shallow_water_two_layer_2d.jl | 81 ++++++++++++++++++- 11 files changed, 140 insertions(+), 22 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl index e9f444aed27..e95eba6a094 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -13,8 +13,8 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal), volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) 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 67c1098b178..ca2a68fd3ed 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -21,8 +21,8 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal), volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) 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 bec0b8ab69c..5398a6b47a2 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -33,8 +33,8 @@ initial_condition = initial_condition_fjordholm_well_balanced ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal), +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg=3, surface_flux=(flux__es_ersing_etal, flux_nonconservative_ersing_etal), volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index 402361c8823..d9fe9675de0 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -12,8 +12,8 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal), volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) 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 ba4bcd25774..e75a8b724b3 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -28,8 +28,8 @@ initial_condition = initial_condition_well_balanced ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg=3, surface_flux=surface_flux, volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index d9e71e4aa44..6fc0ed521a6 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -14,8 +14,8 @@ initial_condition = initial_condition_convergence_test ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg=6, surface_flux=surface_flux, volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) 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 4295f93b342..74114ba6b8e 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -40,8 +40,8 @@ boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_b ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux= (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux= (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg=6, surface_flux=surface_flux, volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) 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 40bc9f2ab42..40388a7c3cd 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -30,8 +30,8 @@ initial_condition = initial_condition_well_balanced ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (flux_es_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal) solver = DGSEM(polydeg=6, surface_flux=surface_flux, volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/src/Trixi.jl b/src/Trixi.jl index c0cecf86bd4..7ed9ef12610 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -147,8 +147,9 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, flux_kennedy_gruber, flux_shima_etal, flux_ec, - flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal, + flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, + flux_nonconservative_ersing_etal, flux_es_ersing_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index fd4fbc017ec..435f065c207 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -234,6 +234,46 @@ Further details are available in the paper: end +""" + 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 modified version of flux_nonconservative_wintermeyer_etal that gives entropy conservation +and well-balancedness in both the volume and surface when combined with flux_wintermeyer_etal. +""" +@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] + +h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) +h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) +b_avg = 0.5 * (b_ll + b_rr) + +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 * 2.0 * (b_avg + h_lower_avg), +z, +equations.gravity * h_lower_ll * 2.0 * (b_avg + equations.r * h_upper_avg), +z) +return f +end + + """ flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) @@ -376,7 +416,7 @@ end """ - flux_es_fjordholm_etal(u_ll, u_rr, orientation, + flux_es_ersing_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy @@ -391,11 +431,11 @@ 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_es_fjordholm_etal(u_ll, u_rr, +@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, + f_ec = flux_wintermeyer_etal(u_ll, u_rr, orientation, equations) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 60c389d8c4a..bab251c695b 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -343,6 +343,83 @@ end end + """ + flux_nonconservative_ersing_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 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 modified version of flux_nonconservative_wintermeyer_etal that gives entropy conservation + and well-balancedness in both the volume and surface when combined with flux_wintermeyer_etal. +""" +@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] + +h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) +h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) +b_avg = 0.5 * (b_ll + b_rr) + +z = zero(eltype(u_ll)) + +# 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)_y, 0) +if orientation == 1 + f = SVector(z, + equations.gravity * h_upper_ll * 2.0 * (b_avg + h_lower_avg), + z,z, + equations.gravity * h_lower_ll * 2.0 * (b_avg + equations.r * h_upper_avg), + z,z) +else # orientation == 2 + f = SVector(z, z, + equations.gravity * h_upper_ll * 2.0 * (b_avg + h_lower_avg), + z,z, + equations.gravity * h_lower_ll * 2.0 * (b_avg + equations.r * h_upper_avg), + z) +end + +return f +end + +@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] + +h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) +h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) +b_avg = 0.5 * (b_ll + b_rr) + +# 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 * 2.0 * (b_avg + h_lower_avg), + normal_direction_average[2] * equations.gravity * h_upper_ll * 2.0 * (b_avg + h_lower_avg), + zero(eltype(u_ll)), + normal_direction_average[1] * equations.gravity * h_lower_ll * 2.0 * (b_avg + + equations.r * h_upper_avg), + normal_direction_average[2] * equations.gravity * h_lower_ll * 2.0 * (b_avg + + equations.r * h_upper_avg), + zero(eltype(u_ll))) +end + + """ flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterTwoLayerEquations2D) @@ -649,11 +726,11 @@ 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_es_fjordholm_etal(u_ll, u_rr, +@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, + f_ec = flux_wintermeyer_etal(u_ll, u_rr, orientation_or_normal_direction, equations) From b2394080f27f9b350cec68b092e0e8742de69b77 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 5 Jun 2023 09:33:05 +0200 Subject: [PATCH 02/17] Change formulation of flux_nonconservative_ersing --- src/equations/shallow_water_two_layer_1d.jl | 11 ++++--- src/equations/shallow_water_two_layer_2d.jl | 34 +++++++++++---------- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 435f065c207..916b6b37a53 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -257,18 +257,19 @@ h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[5] b_ll = u_ll[5] -h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) -h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) -b_avg = 0.5 * (b_ll + b_rr) +# 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, -equations.gravity * h_upper_ll * 2.0 * (b_avg + h_lower_avg), +equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, -equations.gravity * h_lower_ll * 2.0 * (b_avg + equations.r * h_upper_avg), +equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), z) return f end diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index bab251c695b..996c321eba4 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -366,9 +366,10 @@ h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[7] b_ll = u_ll[7] -h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) -h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) -b_avg = 0.5 * (b_ll + b_rr) +# 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)) @@ -377,15 +378,15 @@ z = zero(eltype(u_ll)) # g*h_lower*(b + r*h_upper)_y, 0) if orientation == 1 f = SVector(z, - equations.gravity * h_upper_ll * 2.0 * (b_avg + h_lower_avg), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z,z, - equations.gravity * h_lower_ll * 2.0 * (b_avg + equations.r * h_upper_avg), + 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 * 2.0 * (b_avg + h_lower_avg), + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z,z, - equations.gravity * h_lower_ll * 2.0 * (b_avg + equations.r * h_upper_avg), + equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), z) end @@ -402,20 +403,21 @@ h_upper_rr, h_lower_rr = waterheight(u_rr, equations) b_rr = u_rr[7] b_ll = u_ll[7] -h_upper_avg = 0.5 * (h_upper_ll + h_upper_rr) -h_lower_avg = 0.5 * (h_lower_ll + h_lower_rr) -b_avg = 0.5 * (b_ll + b_rr) +# 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 * 2.0 * (b_avg + h_lower_avg), - normal_direction_average[2] * equations.gravity * h_upper_ll * 2.0 * (b_avg + h_lower_avg), + normal_direction_average[1] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), + normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), zero(eltype(u_ll)), - normal_direction_average[1] * equations.gravity * h_lower_ll * 2.0 * (b_avg + - equations.r * h_upper_avg), - normal_direction_average[2] * equations.gravity * h_lower_ll * 2.0 * (b_avg + - equations.r * h_upper_avg), + normal_direction_average[1] * equations.gravity * h_lower_ll * (b_jump + + equations.r * h_upper_jump), + normal_direction_average[2] * equations.gravity * h_lower_ll* (b_jump + + equations.r * h_upper_jump), zero(eltype(u_ll))) end From 298918ab74e9a32fc8d6c13bab2824b7ed7b30e6 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 5 Jun 2023 09:43:05 +0200 Subject: [PATCH 03/17] remove old fluxes (fjordholm, wintermeyer) --- ...xir_shallowwater_twolayer_well_balanced.jl | 2 +- ...xir_shallowwater_twolayer_well_balanced.jl | 2 +- src/equations/shallow_water_two_layer_1d.jl | 167 +--------- src/equations/shallow_water_two_layer_2d.jl | 311 ++---------------- test/test_tree_1d_shallowwater_twolayer.jl | 12 +- 5 files changed, 33 insertions(+), 461 deletions(-) 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 5398a6b47a2..10ecc739792 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -34,7 +34,7 @@ initial_condition = initial_condition_fjordholm_well_balanced # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux__es_ersing_etal, flux_nonconservative_ersing_etal), +solver = DGSEM(polydeg=3, surface_flux=(flux_es_ersing_etal, flux_nonconservative_ersing_etal), 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 40388a7c3cd..f089e01f48c 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -31,7 +31,7 @@ initial_condition = initial_condition_well_balanced # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_es_ersing_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/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 916b6b37a53..c74d423410c 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -195,45 +195,6 @@ end end -""" - flux_nonconservative_wintermeyer_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. - -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) -""" -@inline function flux_nonconservative_wintermeyer_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] - - 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_ersing_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterTwoLayerEquations1D) @@ -241,12 +202,13 @@ end !!! 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 +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. -This is a modified version of flux_nonconservative_wintermeyer_etal that gives entropy conservation -and well-balancedness in both the volume and surface when combined with flux_wintermeyer_etal. +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). """ @inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, @@ -275,111 +237,13 @@ 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 - - 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), - z, - g * h_lower_ll * (b_ll + equations.r * h_upper_ll) + g * h_lower_average * (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 +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. @@ -417,20 +281,11 @@ end """ - flux_es_ersing_etal(u_ll, u_rr, orientation, + 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_fjordholm_etal 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. +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +flux_wintermeyer_etal and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy +variables. """ @inline function flux_es_ersing_etal(u_ll, u_rr, orientation::Integer, diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 996c321eba4..995f4b3d7e2 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -273,89 +273,19 @@ end """ - flux_nonconservative_wintermeyer_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 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. - -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) -""" -@inline function flux_nonconservative_wintermeyer_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] - - z = zero(eltype(u_ll)) - - # 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)_y, 0) - if orientation == 1 - f = SVector(z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), - z,z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), - z,z) - else # orientation == 2 - f = SVector(z, z, - equations.gravity * h_upper_ll * (b_rr + h_lower_rr), - z,z, - equations.gravity * h_lower_ll * (b_rr + equations.r * h_upper_rr), - 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) - # 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] - - # 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), - normal_direction_average[2] * equations.gravity * h_upper_ll * (b_rr + h_lower_rr), - zero(eltype(u_ll)), - normal_direction_average[1] * equations.gravity * h_lower_ll * (b_rr + - equations.r * h_upper_rr), - normal_direction_average[2] * equations.gravity * h_lower_ll * (b_rr + - equations.r * h_upper_rr), - zero(eltype(u_ll))) - end - - - """ - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + flux_nonconservative_ersing_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 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 modified version of flux_nonconservative_wintermeyer_etal that gives entropy conservation - and well-balancedness in both the volume and surface when combined with flux_wintermeyer_etal. +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 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). """ @inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, @@ -403,6 +333,10 @@ 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) # Calculate jumps h_upper_jump = (h_upper_rr - h_upper_ll) h_lower_jump = (h_lower_rr - h_lower_ll) @@ -411,6 +345,8 @@ 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_jump + h_lower_jump), + normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), normal_direction_average[1] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), zero(eltype(u_ll)), @@ -422,214 +358,13 @@ return SVector(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) 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 +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. @@ -714,19 +449,11 @@ end """ - flux_es_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations1D) + 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_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy +flux_wintermeyer_etal 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. """ @inline function flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index 6c0ad2941cc..97d6152d49c 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -17,16 +17,6 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") tspan = (0.0, 0.25)) 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)) - 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, @@ -41,7 +31,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") l2 = [0.35490827242437256, 1.6715402155795918, 0.6960264969949427, 0.9351481433409805, 0.7938172946965545], linf = [0.6417127471419837, 1.9742107034120873, 1.135774587483082, 1.236125279347084, 1.1], - surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal), tspan = (0.0, 0.25)) end From e09b59cfbeda625b1d3b130d4c6d9ba23dfe41ed Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 5 Jun 2023 10:08:13 +0200 Subject: [PATCH 04/17] Update tests for new flux_ersing_etal --- src/equations/shallow_water_two_layer_2d.jl | 6 ---- test/test_tree_1d_shallowwater_twolayer.jl | 19 +++++------ test/test_tree_2d_shallowwater_twolayer.jl | 36 ++++++++------------- test/test_unstructured_2d.jl | 21 ++++++------ 4 files changed, 34 insertions(+), 48 deletions(-) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 995f4b3d7e2..4ef6bd5f672 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -333,10 +333,6 @@ 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) # Calculate jumps h_upper_jump = (h_upper_rr - h_upper_ll) h_lower_jump = (h_lower_rr - h_lower_ll) @@ -345,8 +341,6 @@ 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_jump + h_lower_jump), - normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), normal_direction_average[1] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), zero(eltype(u_ll)), diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl index 97d6152d49c..fe2e20ff0d1 100644 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ b/test/test_tree_1d_shallowwater_twolayer.jl @@ -10,27 +10,28 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @testset "Shallow Water Two layer" begin @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], + 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)) 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, + l2 = [8.949288784402005e-16, 4.0636427176237915e-17, 0.001002881985401548, 2.133351105037203e-16, 0.0010028819854016578], - linf = [2.6229018956769323e-15, 1.878451903240623e-16, 0.005119880996670156, + linf = [2.6229018956769323e-15, 1.878451903240623e-16, 0.005119880996670156, 8.003199803957679e-16, 0.005119880996670666], tspan = (0.0, 0.25)) 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.35490827242437256, 1.6715402155795918, 0.6960264969949427, - 0.9351481433409805, 0.7938172946965545], - linf = [0.6417127471419837, 1.9742107034120873, 1.135774587483082, 1.236125279347084, 1.1], + l2 = [0.35490300089633237, 1.6715035180645277, 0.6960151891515254, + 0.9351487046099882, 0.7938172946965545], + linf = [0.6417505887480972, 1.9743522856822249, 1.1357745874830805, + 1.2361800390346178, 1.1], surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal), tspan = (0.0, 0.25)) end diff --git a/test/test_tree_2d_shallowwater_twolayer.jl b/test/test_tree_2d_shallowwater_twolayer.jl index 4bb45064714..0d15a709871 100644 --- a/test/test_tree_2d_shallowwater_twolayer.jl +++ b/test/test_tree_2d_shallowwater_twolayer.jl @@ -10,33 +10,23 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") @testset "Two-Layer Shallow Water" begin @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)) - 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), + 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)) 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], + 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)) end @@ -48,7 +38,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") 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), + surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal), tspan = (0.0, 0.25)) end end diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index d4b0d150ca1..5894395626d 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -157,11 +157,12 @@ isdir(outdir) && rm(outdir, recursive=true) @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)) end @@ -178,12 +179,12 @@ isdir(outdir) && rm(outdir, recursive=true) @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], - surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + 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_ersing_etal), tspan = (0.0, 0.25)) end end From c6a441a8d6845fa34bc5a326c8aa13b6c128ee85 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 5 Jun 2023 10:51:24 +0200 Subject: [PATCH 05/17] Add flux_nonconservative_ersing for standard SWE --- src/equations/shallow_water_1d.jl | 33 ++++++++++++++++ src/equations/shallow_water_2d.jl | 62 +++++++++++++++++++++++++++++++ test/test_tree_1d_shallowwater.jl | 18 +++++++++ test/test_tree_2d_shallowwater.jl | 18 +++++++++ test/test_unstructured_2d.jl | 18 +++++++++ 5 files changed, 149 insertions(+) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 949c6576006..96ae6cbbedf 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -307,6 +307,39 @@ Further details on the hydrostatic reconstruction and its motivation can be foun 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). +""" +@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 b07fbfc739e..1e54ad5af2a 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -520,6 +520,68 @@ end 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). +""" +@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, diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index f8901a3dcb6..edbeb679498 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -30,6 +30,15 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") tspan = (0.0, 0.25)) 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.2427984842961743, 1.6265844732207998e-14, 1.2427984842961741], + linf = [1.6190414782447629, 3.30154831542628e-14, 1.6190414782447629], + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + tspan = (0.0, 0.25)) + end + @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2 = [0.0022363707373868713, 0.01576799981934617, 4.436491725585346e-5], @@ -44,6 +53,15 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") tspan = (0.0, 0.025), surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal)) 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)) + end + @trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2 = [0.0022851099219788917, 0.01560453773635554, 4.43649172558535e-5], diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index f465a177a67..d75ea720ab7 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -37,6 +37,15 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") tspan = (0.0, 0.25)) 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)) + end + @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2 = [0.001868474306068482, 0.01731687445878443, 0.017649083171490863, 6.274146767717023e-5], @@ -57,6 +66,15 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") linf = [0.015156105797771602, 0.07964811135780492, 0.0839787097210376, 0.0001819675955490041], tspan = (0.0, 0.025), surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal)) 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)) + end end end # module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 5894395626d..b67ac1dcb43 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -119,6 +119,15 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.2)) 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.4034438316613565e-12, 1.067841401808785e-12, 1.2164292510839079], + linf = [1.5138512282315757, 4.24227234830862e-11, 1.9765812927086113e-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)) + end + @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2 = [0.0011197623982310795, 0.04456344888447023, 0.014317376629669337, 5.089218476758975e-6], @@ -134,6 +143,15 @@ isdir(outdir) && rm(outdir, recursive=true) tspan = (0.0, 0.025)) 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)) + end + @trixi_testset "elixir_shallowwater_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), l2 = [1.1577518608940115e-5, 4.867189932537344e-13, 4.647273240470541e-13, 1.1577518608933468e-5], From 43d6608b5a68a652dc7b3cc74e5162d5eecf8a33 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 5 Jun 2023 11:01:36 +0200 Subject: [PATCH 06/17] Edit comments, change flux in elixir --- .../elixir_shallowwater_twolayer_well_balanced.jl | 2 +- src/equations/shallow_water_1d.jl | 2 +- src/equations/shallow_water_two_layer_1d.jl | 4 ++-- src/equations/shallow_water_two_layer_2d.jl | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) 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 f089e01f48c..40388a7c3cd 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -31,7 +31,7 @@ initial_condition = initial_condition_well_balanced # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_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/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 96ae6cbbedf..e9fb1a47775 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -334,7 +334,7 @@ conservation and well-balancedness in both the volume and surface when combined 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) + f = SVector(z, equations.gravity * h_ll * b_jump, z) return f end diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index c74d423410c..60c877568e5 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -284,8 +284,8 @@ end 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 and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy -variables. +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. """ @inline function flux_es_ersing_etal(u_ll, u_rr, orientation::Integer, diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 4ef6bd5f672..d1bfb0f1b4a 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -446,8 +446,8 @@ end 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 and adds a Lax-Friedrichs type dissipation dependent on the jump of entropy -variables. +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. """ @inline function flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, From adc68064f2be6e3bdf0894e65edd7c4569f25caf Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 5 Jun 2023 15:40:37 +0200 Subject: [PATCH 07/17] fix well-balanced test --- test/test_unstructured_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index b67ac1dcb43..b7b60e5c54d 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -121,8 +121,8 @@ isdir(outdir) && rm(outdir, recursive=true) @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.4034438316613565e-12, 1.067841401808785e-12, 1.2164292510839079], - linf = [1.5138512282315757, 4.24227234830862e-11, 1.9765812927086113e-11, 1.513851228231574], + 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)) From 28668fe2b354685832ef4a5537f2379aa7802484 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Thu, 15 Jun 2023 10:50:35 +0200 Subject: [PATCH 08/17] Update src/Trixi.jl Change order of exporting statements Co-authored-by: Andrew Winters --- src/Trixi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 92aea659b65..14e56900ad1 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -152,7 +152,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, - flux_nonconservative_ersing_etal, flux_es_ersing_etal, + flux_es_ersing_etal, flux_nonconservative_ersing_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, From 7d75332d1eb1396dc1d0070ec2f53c8d2501853c Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 15 Jun 2023 11:08:27 +0200 Subject: [PATCH 09/17] fix indentation --- src/equations/shallow_water_two_layer_1d.jl | 58 ++++----- src/equations/shallow_water_two_layer_2d.jl | 128 ++++++++++---------- 2 files changed, 94 insertions(+), 92 deletions(-) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 60c877568e5..17a85858298 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -197,7 +197,7 @@ end """ flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) + equations::ShallowWaterTwoLayerEquations1D) !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. @@ -211,29 +211,29 @@ conservation and well-balancedness in both the volume and surface when combined [`flux_wintermeyer_etal`](@ref). """ @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] - -# 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, -equations.gravity * h_upper_ll * (b_jump + h_lower_jump), -z, -equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), -z) -return f + 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] + + # 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, + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), + z, + equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), + z) + return f end @@ -282,18 +282,18 @@ end """ flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations1D) + 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. """ @inline function flux_es_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) # Compute entropy conservative flux but without the bottom topography f_ec = flux_wintermeyer_etal(u_ll, u_rr, - orientation, - equations) + orientation, + equations) # Get maximum signal velocity λ = max_abs_speed_naive(u_ll, u_rr, orientation, equations) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index d1bfb0f1b4a..c746bf51324 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -274,7 +274,7 @@ end """ flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) + equations::ShallowWaterTwoLayerEquations2D) !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. @@ -288,67 +288,69 @@ conservation and well-balancedness in both the volume and surface when combined [`flux_wintermeyer_etal`](@ref). """ @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)) - -# 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)_y, 0) -if orientation == 1 - f = SVector(z, - equations.gravity * h_upper_ll * (b_jump + h_lower_jump), - z,z, - 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_jump + h_lower_jump), - z,z, - equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), - z) -end + 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)) -return f + # 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)_y, 0) + if orientation == 1 + f = SVector(z, + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), + z,z, + 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_jump + h_lower_jump), + z,z, + equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_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::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_jump + h_lower_jump), - normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), - zero(eltype(u_ll)), - normal_direction_average[1] * equations.gravity * h_lower_ll * (b_jump + - equations.r * h_upper_jump), - normal_direction_average[2] * equations.gravity * h_lower_ll* (b_jump + - equations.r * h_upper_jump), - zero(eltype(u_ll))) + 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_jump + + h_lower_jump), + normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + + h_lower_jump), + zero(eltype(u_ll)), + normal_direction_average[1] * equations.gravity * h_lower_ll * (b_jump + + equations.r * h_upper_jump), + normal_direction_average[2] * equations.gravity * h_lower_ll * (b_jump + + equations.r * h_upper_jump), + zero(eltype(u_ll))) end @@ -444,18 +446,18 @@ end """ flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) + 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. """ @inline function flux_es_ersing_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) + orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) # Compute entropy conservative flux but without the bottom topography f_ec = flux_wintermeyer_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations) + orientation_or_normal_direction, + equations) # Get maximum signal velocity λ = max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations) From fb489349d0bef73756d4caa1666bb43fe5ab944a Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 15 Jun 2023 11:12:17 +0200 Subject: [PATCH 10/17] fix indentation II --- src/equations/shallow_water_1d.jl | 4 ++-- src/equations/shallow_water_2d.jl | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index e9fb1a47775..874d9c5ef0e 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -309,7 +309,7 @@ end """ flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations1D) + equations::ShallowWaterEquations1D) !!! warning "Experimental code" This numerical flux is experimental and may change in any future release. @@ -322,7 +322,7 @@ conservation and well-balancedness in both the volume and surface when combined [`flux_wintermeyer_etal`](@ref). """ @inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations1D) + equations::ShallowWaterEquations1D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) b_rr = u_rr[3] diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 1e54ad5af2a..7758b1eb1ee 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -544,7 +544,7 @@ conservation and well-balancedness in both the volume and surface when combined [`flux_wintermeyer_etal`](@ref). """ @inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations2D) + equations::ShallowWaterEquations2D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) b_rr = u_rr[4] @@ -564,9 +564,9 @@ conservation and well-balancedness in both the volume and surface when combined end @inline function flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterEquations2D) + 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] From 9274c3676f39a14298277735ef1546e705a13ee8 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Thu, 15 Jun 2023 14:33:42 +0200 Subject: [PATCH 11/17] Update src/equations/shallow_water_2d.jl Co-authored-by: Andrew Winters --- src/equations/shallow_water_2d.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 7758b1eb1ee..a9a7b67e8e1 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -522,11 +522,11 @@ end """ flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations2D) + equations::ShallowWaterEquations2D) flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterEquations2D) + 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. From ba45d05531a209e51e1b8c75a66fe7e0f63ff083 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 19 Jun 2023 09:42:14 +0200 Subject: [PATCH 12/17] fix energy total function --- src/equations/shallow_water_two_layer_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 17a85858298..7af26e98a4b 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -454,7 +454,7 @@ end # Calculate total energy for a conservative state `cons` @inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_lower, h_v1_upper, h_v2_lower, b = cons + h_upper, h_v1_upper, h_lower, h_v2_lower, b = cons # Set new variables for better readability g = equations.gravity rho_upper = equations.rho_upper From 2381c2b7bbe87b09092126d5229e2509e40981c8 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 21 Jun 2023 16:24:31 +0200 Subject: [PATCH 13/17] Apply formatter --- src/equations/shallow_water_1d.jl | 21 ++-- src/equations/shallow_water_2d.jl | 7 +- src/equations/shallow_water_two_layer_1d.jl | 17 ++- src/equations/shallow_water_two_layer_2d.jl | 121 ++++++++++---------- 4 files changed, 82 insertions(+), 84 deletions(-) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 76007c344f7..e93aa7dca24 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -323,23 +323,22 @@ conservation and well-balancedness in both the volume and surface when combined """ @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] + # 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 + # Calculate jump + b_jump = b_rr - b_ll - z = zero(eltype(u_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) + # Bottom gradient nonconservative term: (0, g h b_x, 0) + f = SVector(z, equations.gravity * h_ll * b_jump, z) - return f + 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 cd0d35dec93..f85e11f1b6b 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -545,10 +545,10 @@ conservation and well-balancedness in both the volume and surface when combined 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 @@ -556,7 +556,7 @@ conservation and well-balancedness in both the volume and surface when combined else # orientation == 2 f = SVector(z, z, equations.gravity * h_ll * b_jump, z) end -return f + return f end @inline function flux_nonconservative_ersing_etal(u_ll, u_rr, @@ -578,7 +578,6 @@ end 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 43ea3c88d8e..8b7765be6e5 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -235,9 +235,9 @@ conservation and well-balancedness in both the volume and surface when combined b_ll = u_ll[5] # 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) + 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)) @@ -251,7 +251,6 @@ conservation and well-balancedness in both the volume and surface when combined return f end - """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations1D) @@ -434,11 +433,11 @@ end 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) + 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) + w3 = (equations.rho_lower * + (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v1_lower^2)) w4 = equations.rho_lower * v1_lower return SVector(w1, w2, w3, w4, b) end @@ -481,7 +480,7 @@ end @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) + 0.5 * equations.rho_lower * h_v2_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 4ce011eccec..970bba2188a 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -341,9 +341,9 @@ conservation and well-balancedness in both the volume and surface when combined 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) + 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)) @@ -354,13 +354,15 @@ conservation and well-balancedness in both the volume and surface when combined f = SVector(z, equations.gravity * h_upper_ll * (b_jump + h_lower_jump), z, z, - equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), + 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_jump + h_lower_jump), z, z, - equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), z) end @@ -376,28 +378,27 @@ end 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) + 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 * + normal_direction_average[1] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), - normal_direction_average[2] * equations.gravity * h_upper_ll * + normal_direction_average[2] * equations.gravity * h_upper_ll * (b_jump + h_lower_jump), zero(eltype(u_ll)), - normal_direction_average[1] * equations.gravity * h_lower_ll * + normal_direction_average[1] * equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), - normal_direction_average[2] * equations.gravity * h_lower_ll * + normal_direction_average[2] * equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), zero(eltype(u_ll))) end - """ flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterTwoLayerEquations2D) @@ -417,50 +418,51 @@ 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, +@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 -end + # 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 +end """ flux_wintermeyer_etal(u_ll, u_rr, orientation, @@ -713,7 +715,6 @@ end return max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr) * norm(normal_direction) end - # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, @@ -776,12 +777,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) From 641a5923ede98807f9c9dd30af75c5bc57261455 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 27 Jun 2023 19:13:57 +0200 Subject: [PATCH 14/17] add docstring references --- src/equations/shallow_water_1d.jl | 6 ++++++ src/equations/shallow_water_2d.jl | 6 ++++++ src/equations/shallow_water_two_layer_1d.jl | 12 ++++++++++++ src/equations/shallow_water_two_layer_2d.jl | 12 ++++++++++++ 4 files changed, 36 insertions(+) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index e93aa7dca24..3071d96943e 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -320,6 +320,12 @@ that contains the gradient of the bottom topography [`ShallowWaterEquations1D`]( 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) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index f85e11f1b6b..283cc3e2467 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -538,6 +538,12 @@ 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) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 8b7765be6e5..a3ef3f82322 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -224,6 +224,12 @@ additional term that couples the momentum of both layers. 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, @@ -299,6 +305,12 @@ end 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_ersing_etal(u_ll, u_rr, orientation::Integer, diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 970bba2188a..2883a6fd680 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -330,6 +330,12 @@ additional term that couples the momentum of both layers. 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, @@ -560,6 +566,12 @@ end 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_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, From 66e9c60d8984231228f617b4681977fa37c3f0b8 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 28 Jun 2023 10:41:57 +0200 Subject: [PATCH 15/17] remove flux_nonconservative_fjordholm_etal --- src/equations/shallow_water_two_layer_2d.jl | 64 --------------------- 1 file changed, 64 deletions(-) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 2883a6fd680..9458e79b6fe 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -405,70 +405,6 @@ end 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 -end """ flux_wintermeyer_etal(u_ll, u_rr, orientation, From fae20c6493c7a85a65dbc9af85bb9485bbb94940 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 28 Jun 2023 10:46:49 +0200 Subject: [PATCH 16/17] apply formatter --- src/equations/shallow_water_two_layer_2d.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 9458e79b6fe..365c3df6bc9 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -405,7 +405,6 @@ end zero(eltype(u_ll))) end - """ flux_wintermeyer_etal(u_ll, u_rr, orientation, equations::ShallowWaterTwoLayerEquations2D) From 840f94e63bda9be4f73edd0ba3b44c253f4ddde9 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 6 Nov 2023 12:41:57 +0100 Subject: [PATCH 17/17] Format examples --- .../elixir_shallowwater_twolayer_convergence.jl | 6 +++--- .../tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl | 5 +++-- .../elixir_shallowwater_twolayer_well_balanced.jl | 5 +++-- .../elixir_shallowwater_twolayer_convergence.jl | 6 +++--- .../elixir_shallowwater_twolayer_well_balanced.jl | 4 ++-- .../elixir_shallowwater_twolayer_convergence.jl | 4 ++-- .../elixir_shallowwater_twolayer_dam_break.jl | 6 +++--- .../elixir_shallowwater_twolayer_well_balanced.jl | 4 ++-- 8 files changed, 21 insertions(+), 19 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl index 2098a92ce51..e6a01849852 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -14,9 +14,9 @@ initial_condition = initial_condition_convergence_test # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal), - volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) - +solver = DGSEM(polydeg = 3, + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # Get the TreeMesh and setup a periodic mesh 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 2be63b564f7..03b93754d0f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -37,8 +37,9 @@ initial_condition = initial_condition_dam_break # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal), - volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # Get the TreeMesh and setup a non-periodic mesh 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 9027960b643..098e3aaf601 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -36,8 +36,9 @@ initial_condition = initial_condition_fjordholm_well_balanced # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_es_ersing_etal, flux_nonconservative_ersing_etal), - volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # Get the TreeMesh and setup a periodic mesh diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index e9ae931c688..790916e4467 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -14,9 +14,9 @@ initial_condition = initial_condition_convergence_test # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg=3, surface_flux=(flux_wintermeyer_etal, flux_nonconservative_ersing_etal), - volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) - +solver = DGSEM(polydeg = 3, + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # Get the TreeMesh and setup a periodic mesh 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 f2d46b292eb..264c26390fe 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -33,8 +33,8 @@ initial_condition = initial_condition_well_balanced 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)) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # Get the TreeMesh and setup a periodic mesh diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl index 43e1e46184c..0b86095663a 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -17,8 +17,8 @@ initial_condition = initial_condition_convergence_test 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)) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # This setup is for the curved, split form convergence test on a periodic domain 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 2226f5115d7..4ad5f7e3201 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -42,9 +42,9 @@ boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_b # Get the DG approximation space 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)) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # Get the unstructured quad mesh from a file (downloads the file if not available locally) 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 f4123ca58cd..6a727df2502 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -35,8 +35,8 @@ initial_condition = initial_condition_well_balanced 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)) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### # This setup is for the curved, split form well-balancedness testing