From 8b2bab919c9e389113b529880552294cc0dcb8ba Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 4 Jan 2024 15:52:16 +0100 Subject: [PATCH 1/5] fix_wall_bc --- .../tree_2d_dgsem/elixir_shallowwater_wall.jl | 85 +++++++++++++++++++ src/equations/shallow_water_2d.jl | 10 ++- 2 files changed, 92 insertions(+), 3 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_wall.jl diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl new file mode 100644 index 00000000000..28a27b60577 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -0,0 +1,85 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) + +# An initial condition with a perturbation in the waterheight to test boundary_condition_slip_wall +function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + + x1, x2 = x + b = 0.0 + + # Waterheight perturbation + r = (x[1])^2 + (x[2])^2 + H = H + 0.5 * exp(-10 * r) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.25) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index e75c92a27d0..850543e6f82 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -236,9 +236,13 @@ Should be used together with [`TreeMesh`](@ref). u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4]) end - # compute and return the flux using `boundary_condition_slip_wall` routine above - flux = surface_flux_function(u_inner, u_boundary, orientation, equations) - + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + return flux end From d6ce8f8e46e196ce4cd6f4114e4006c8d5aa8cfe Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 4 Jan 2024 16:17:36 +0100 Subject: [PATCH 2/5] add test --- test/test_tree_2d_shallowwater.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 58db7c5f35f..835d77909d1 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -327,6 +327,21 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_shallowwater_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), + l2 = [0.10982800604156373, 0.2696758586410491, 0.26967585864104954, 0.0], + linf = [0.4962808878041378, 0.6577212830580099, 0.6577212830580086, 0.0], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end end # module From 5412981a58f9d536cd83fe4773f75957fa68d5b5 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 4 Jan 2024 16:18:59 +0100 Subject: [PATCH 3/5] apply formatter --- src/equations/shallow_water_2d.jl | 2 +- test/test_tree_2d_shallowwater.jl | 16 +++++++++++++--- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 850543e6f82..6728d7d5553 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -242,7 +242,7 @@ Should be used together with [`TreeMesh`](@ref). else # u_boundary is "left" of boundary, u_inner is "right" of boundary flux = surface_flux_function(u_boundary, u_inner, orientation, equations) end - + return flux end diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 835d77909d1..8e5bdb69334 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -330,9 +330,19 @@ end @trixi_testset "elixir_shallowwater_wall.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), - l2 = [0.10982800604156373, 0.2696758586410491, 0.26967585864104954, 0.0], - linf = [0.4962808878041378, 0.6577212830580099, 0.6577212830580086, 0.0], - tspan=(0.0, 0.25)) + l2=[ + 0.10982800604156373, + 0.2696758586410491, + 0.26967585864104954, + 0.0, + ], + linf=[ + 0.4962808878041378, + 0.6577212830580099, + 0.6577212830580086, + 0.0, + ], + tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 5133ca09e232942d5447428bc7a797ec0cc36a42 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 5 Jan 2024 08:36:21 +0100 Subject: [PATCH 4/5] adjust some comments --- examples/tree_2d_dgsem/elixir_shallowwater_wall.jl | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl index 28a27b60577..4f4be8083b6 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -2,8 +2,7 @@ using OrdinaryDiffEq using Trixi ############################################################################### -# semidiscretization of the shallow water equations with a discontinuous -# bottom topography function +# Semidiscretization of the shallow water equations equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) @@ -37,7 +36,7 @@ solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) ############################################################################### -# Get the TreeMesh and setup a periodic mesh +# Get the TreeMesh and setup a non-periodic mesh coordinates_min = (-1.0, -1.0) coordinates_max = (1.0, 1.0) @@ -62,8 +61,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - extra_analysis_integrals = (lake_at_rest_error,)) +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) From 7d03b4a65dfa508d5367916abd9a7497d1fa6407 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 8 Jan 2024 08:52:55 +0100 Subject: [PATCH 5/5] update elixir and test; add topography --- .../tree_2d_dgsem/elixir_shallowwater_wall.jl | 15 +++++++-------- test/test_tree_2d_shallowwater.jl | 16 ++++++++-------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl index 4f4be8083b6..b8dbad50680 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -6,24 +6,23 @@ using Trixi equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) -# An initial condition with a perturbation in the waterheight to test boundary_condition_slip_wall -function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D) +# An initial condition with a bottom topography and a perturbation in the waterheight to test +# boundary_condition_slip_wall +function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D) # Set the background values H = equations.H0 v1 = 0.0 v2 = 0.0 - x1, x2 = x - b = 0.0 - + # Bottom topography + b = 1.5 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) # Waterheight perturbation - r = (x[1])^2 + (x[2])^2 - H = H + 0.5 * exp(-10 * r) + H = H + 0.5 * exp(-10.0 * ((x[1])^2 + (x[2])^2)) return prim2cons(SVector(H, v1, v2, b), equations) end -initial_condition = initial_condition_well_balancedness +initial_condition = initial_condition_perturbation boundary_condition = boundary_condition_slip_wall diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 8e5bdb69334..1f3dfbf5267 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -331,16 +331,16 @@ end @trixi_testset "elixir_shallowwater_wall.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), l2=[ - 0.10982800604156373, - 0.2696758586410491, - 0.26967585864104954, - 0.0, + 0.13517233723296504, + 0.20010876311162215, + 0.20010876311162223, + 2.719538414346464e-7, ], linf=[ - 0.4962808878041378, - 0.6577212830580099, - 0.6577212830580086, - 0.0, + 0.5303607982988336, + 0.5080989745682338, + 0.5080989745682352, + 1.1301675764130437e-6, ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations