Tree Mesh 2D Dirichlet (#1251)

* 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 <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_ec.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl Co-authored-by: Andrew Winters <[email protected]> * Update examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/Trixi.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * Update src/equations/shallow_water_1d.jl Co-authored-by: Andrew Winters <[email protected]> * 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 <[email protected]> * Updated comment wall elixir Co-authored-by: Andrew Winters <[email protected]> * 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 <[email protected]> * Updated TreeMesh Wall Co-authored-by: svengoldberg <[email protected]> Co-authored-by: svengoldberg <[email protected]> Co-authored-by: Andrew Winters <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Nov 30, 2022 · 514ef46 · 514ef46
1 parent 06c826e
commit 514ef46
Show file tree

Hide file tree

Showing 5 changed files with 258 additions and 6 deletions.
diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl 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
diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl 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
diff --git 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
@@ -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
@@ -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],