From 514ef46db664bcadc676353fb92e94f9bd838213 Mon Sep 17 00:00:00 2001 From: maxb <101979498+maxbertrand1996@users.noreply.github.com> Date: Wed, 30 Nov 2022 17:13:12 +0100 Subject: [PATCH] Tree Mesh 2D Dirichlet (#1251) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * First Version Shallow Water 1D * Testdatei * Testnachricht gelöscht * Testdatei geloescht * added examples tree_1d well balanced & blast wave * Add Elixier Tree_1D SWE well ball. & blast wave * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl Co-authored-by: Andrew Winters * Update examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl Co-authored-by: Andrew Winters * Update src/Trixi.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters * WIP: create test files for SWE 1D * Working test file 1d shallowwater * SWE dirichlet WIP: no errors run, EOC h_v bad * SWE Dirichlet Source Term + Well Balanced * WIP: SWE Dirichlet EC + Slip Wall * Well balanced property for SWE slip wall * SWE non-periodic remaining files * Change requests from PR except second last * Removed normal_direction fluxes + update docs * Combine well-bal. nonperiodic, code style update * Boundary Tree 2d fix * 2d elixirs SWE sourceterms dirichlet added * Well balanced wall example * Minor docstring adjustment slip wall * wall elixir comment update Co-authored-by: Andrew Winters * Updated comment wall elixir Co-authored-by: Andrew Winters * Suggested changes changes slip wall * Possible solution TreeMesh? * Well balanced test Tree 2d * Fix convergence_test * Exclude strucctured mesh part Co-authored-by: Andrew Winters * Updated TreeMesh Wall Co-authored-by: svengoldberg Co-authored-by: svengoldberg <102215246+svengoldberg@users.noreply.github.com> Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- ...xir_shallowwater_source_terms_dirichlet.jl | 61 +++++++++ .../elixir_shallowwater_well_balanced_wall.jl | 119 ++++++++++++++++++ src/equations/shallow_water_2d.jl | 27 +++- src/solvers/dgsem_tree/dg_2d.jl | 43 ++++++- test/test_tree_2d_shallowwater.jl | 14 +++ 5 files changed, 258 insertions(+), 6 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl new file mode 100644 index 00000000000..8506334c146 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl @@ -0,0 +1,61 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquations2D(gravity_constant=9.81) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg=3, surface_flux=(flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + volume_integral=VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=3, + n_cells_max=10_000, + periodicity=false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition, + source_terms=source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=200, + save_initial_solution=true, + save_final_solution=true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-8, reltol=1.0e-8, + save_everystep=false, callback=callbacks); +summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl new file mode 100644 index 00000000000..66cd27f6864 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl @@ -0,0 +1,119 @@ +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 constant total water height and zero velocities to test well-balancedness. +# Note, this routine is used to compute errors in the analysis callback but the initialization is +# overwritten by `initial_condition_discontinuous_well_balancedness` below. +function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = ( 1.5 / exp( 0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2) ) + + 0.75 / exp( 0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2) ) ) + 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_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg=4, 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=2, + 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, 100.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase of well-balancedness. +# The errors from the analysis callback are not important but the error for this lake at rest test case +# `∑|H0-(h+b)|` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the TreeMesh2D with initial_refinement_level=2. +function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + b = 0.0 + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[2]) + end + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for j in eachnode(semi.solver), i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, j, element) + u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# 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=3.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 \ No newline at end of file diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 6d17c00d51a..d98afd8151e 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -166,11 +166,9 @@ end """ boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, equations::ShallowWaterEquations2D) - Create a boundary state by reflecting the normal velocity component and keep the tangential velocity component unchanged. The boundary water height is taken from the internal value. - For details see Section 9.2.5 of the book: - Eleuterio F. Toro (2001) Shock-Capturing Methods for Free-Surface Shallow Flows @@ -200,6 +198,29 @@ For details see Section 9.2.5 of the book: end +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::ShallowWaterEquations2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::ShallowWaterEquations2D) + ## get the appropriate normal vector from the orientation + if orientation == 1 + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4]) + else # orientation == 2 + 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) + + return flux +end + # Calculate 1D flux for a single point # Note, the bottom topography has no flux @inline function flux(u, orientation::Integer, equations::ShallowWaterEquations2D) @@ -555,7 +576,7 @@ end v1_avg = 0.5 * (v1_ll + v1_rr ) v2_avg = 0.5 * (v2_ll + v2_rr ) h2_avg = 0.5 * (h_ll^2 + h_rr^2) - p_avg = 0.5 * equations.gravity * h2_avg + p_avg = 0.5 * equations.gravity * h2_avg v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr) # Calculate fluxes depending on normal_direction diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index d725a8d828e..9a3b31ad9f0 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -490,7 +490,6 @@ end end -# We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache, u, mesh::TreeMesh{2}, equations, surface_integral, dg::DG) @unpack interfaces = cache @@ -651,21 +650,25 @@ function calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple, # Calc boundary fluxes in each direction calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[1], + have_nonconservative_terms(equations), equations, surface_integral, dg, cache, 1, firsts[1], lasts[1]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[2], + have_nonconservative_terms(equations), equations, surface_integral, dg, cache, 2, firsts[2], lasts[2]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[3], + have_nonconservative_terms(equations), equations, surface_integral, dg, cache, 3, firsts[3], lasts[3]) calc_boundary_flux_by_direction!(surface_flux_values, t, boundary_conditions[4], + have_nonconservative_terms(equations), equations, surface_integral, dg, cache, 4, firsts[4], lasts[4]) end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, - boundary_condition, equations, + boundary_condition, nonconservative_terms::Val{false}, equations, surface_integral ,dg::DG, cache, direction, first_boundary, last_boundary) @unpack surface_flux = surface_integral @@ -697,6 +700,41 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A return nothing end +function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, + boundary_condition, nonconservative_terms::Val{true}, equations, + surface_integral ,dg::DG, cache, + direction, first_boundary, last_boundary) + surface_flux, nonconservative_flux = surface_integral.surface_flux + @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries + + @threaded for boundary in first_boundary:last_boundary + # Get neighboring element + neighbor = neighbor_ids[boundary] + + for i in eachnode(dg) + # Get boundary flux + u_ll, u_rr = get_surface_node_vars(u, equations, dg, i, boundary) + if neighbor_sides[boundary] == 1 # Element is on the left, boundary on the right + u_inner = u_ll + else # Element is on the right, boundary on the left + u_inner = u_rr + end + x = get_node_coords(node_coordinates, equations, dg, i, boundary) + flux = boundary_condition(u_inner, orientations[boundary], direction, x, t, surface_flux, + equations) + noncons_flux = boundary_condition(u_inner, orientations[boundary], direction, x, t, nonconservative_flux, + equations) + + # Copy flux to left and right element storage + for v in eachvariable(equations) + surface_flux_values[v, i, direction, neighbor] = flux[v] + 0.5 * noncons_flux[v] + end + end + end + + return nothing +end + function prolong2mortars!(cache, u, mesh::TreeMesh{2}, equations, @@ -962,7 +1000,6 @@ end end -# we pass in the hyperbolic `dg.surface_integral` as a dummy argument for dispatch function calc_surface_integral!(du, u, mesh::Union{TreeMesh{2}, StructuredMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DG, cache) diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 15578e77569..8b8af739e8b 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -22,6 +22,13 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") tspan = (0.0, 0.25)) end + @trixi_testset "elixir_shallowwater_well_balanced_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wall.jl"), + l2 = [0.9130579602987144, 1.0602847041965408e-14, 1.082225645390032e-14, 0.9130579602987147], + linf = [2.113062037615659, 4.6613606802974e-14, 5.4225772771633196e-14, 2.1130620376156584], + tspan = (0.0, 0.25)) + end + @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2 = [0.9130579602987147, 9.68729463970494e-15, 9.694538537436981e-15, 0.9130579602987147], @@ -37,6 +44,13 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") 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.0018746929418489125, 0.017332321628469628, 0.01634953679145536, 6.274146767717023e-5], + linf = [0.016262353691956388, 0.08726160620859424, 0.09043621801418844, 0.0001819675955490041], + tspan = (0.0, 0.025)) + end + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2 = [0.0018957692481057034, 0.016943229710439864, 0.01755623297390675, 6.274146767717414e-5],