New slip wall boundary condition for compressible Euler 2D/3D and APE

NEWS.md

@@ -16,6 +16,8 @@ for human readability.
 
 - Implementation of acoustic perturbation equations now uses the conservative form, i.e. the
   perturbed pressure `p_prime` has been replaced with `p_prime_scaled = p_prime / c_mean^2`.
+- Removed the experimental `BoundaryConditionWall` and instead directly compute slip wall boundary
+  condition flux term using the function `boundary_condition_slip_wall`.
 
 #### Deprecated
 

docs/src/hohqmesh_tutorial.md

@@ -300,7 +300,6 @@ initial_condition = uniform_flow_state
 
 # boundary condition types
 boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
-boundary_condition_slip_wall = BoundaryConditionWall(boundary_state_slip_wall)
 
 # boundary condition dictionary
 boundary_conditions = Dict( :Bottom     => boundary_condition_uniform_flow,

examples/dgmulti_2d/elixir_euler_rayleigh_taylor_instability.jl

@@ -70,7 +70,7 @@ vy = map(x-> 0.5 * (1+x), vy) # map [-1, 1] to [0, 1] for single mode RTI
 vertex_coordinates = (vx, vy)
 mesh = VertexMappedMesh(vertex_coordinates, EToV, dg; is_periodic=(true,false))
 
-boundary_conditions = (; :entire_boundary => BoundaryConditionWall(boundary_state_slip_wall))
+boundary_conditions = (; :entire_boundary => boundary_condition_slip_wall)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_rayleigh_taylor_instability, dg;
                                     source_terms = source_terms_rayleigh_taylor_instability,

examples/p4est_3d_dgsem/elixir_euler_ec.jl

@@ -10,7 +10,6 @@ equations = CompressibleEulerEquations3D(5/3)
 
 initial_condition = initial_condition_weak_blast_wave
 
-boundary_condition_slip_wall = BoundaryConditionWall(boundary_state_slip_wall)
 boundary_conditions = Dict(
   :all => boundary_condition_slip_wall
 )

examples/structured_2d_dgsem/elixir_euler_rayleigh_taylor_instability.jl

@@ -64,29 +64,29 @@ solver = DGSEM(polydeg=3, surface_flux=flux_hll,
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
 # The domain is [0, 0.25] x [0, 1]
-mapping(xi, eta) = SVector(0.25 * 0.5 * (1 + xi), 0.5 * (1 + eta))
+mapping(xi, eta) = SVector(0.25 * 0.5 * (1.0 + xi), 0.5 * (1.0 + eta))
 
 num_elements_per_dimension = 32
 cells_per_dimension = (num_elements_per_dimension, num_elements_per_dimension * 4)
 mesh = StructuredMesh(cells_per_dimension, mapping)
 
-boundary_condition_wall = BoundaryConditionWall(boundary_state_slip_wall)
+initial_condition = initial_condition_rayleigh_taylor_instability
+
 boundary_conditions = (
-                       x_neg=boundary_condition_wall,
-                       x_pos=boundary_condition_wall,
-                       y_neg=boundary_condition_wall,
-                       y_pos=boundary_condition_wall,
+                       x_neg=boundary_condition_slip_wall,
+                       x_pos=boundary_condition_slip_wall,
+                       y_neg=boundary_condition_slip_wall,
+                       y_pos=boundary_condition_slip_wall,
                       )
 
-## Alternative setup: left/right periodic BCs and Dirichlet BCs on the top/bottom.
+# # Alternative setup: left/right periodic BCs and Dirichlet BCs on the top/bottom.
 # boundary_conditions = (
 #                        x_neg=boundary_condition_periodic,
 #                        x_pos=boundary_condition_periodic,
 #                        y_neg=BoundaryConditionDirichlet(initial_condition),
 #                        y_pos=BoundaryConditionDirichlet(initial_condition),
 #                       )
 
-initial_condition = initial_condition_rayleigh_taylor_instability
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                     source_terms=source_terms_rayleigh_taylor_instability,
                                     boundary_conditions=boundary_conditions)

examples/tree_2d_dgsem/elixir_euler_ec.jl

@@ -12,14 +12,16 @@ volume_flux = flux_ranocha
 solver = DGSEM(polydeg=3, surface_flux=flux_ranocha,
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
-coordinates_min = (-2, -2)
-coordinates_max = ( 2,  2)
+coordinates_min = (-2.0, -2.0)
+coordinates_max = ( 2.0,  2.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level=5,
-                n_cells_max=10_000)
+                n_cells_max=10_000,
+                periodicity=true)
 
 
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_condition_periodic)
 
 
 ###############################################################################

examples/tree_2d_dgsem/elixir_lbm_couette.jl

@@ -11,7 +11,7 @@ initial_condition = initial_condition_couette_unsteady
 boundary_conditions = (
                        x_neg=boundary_condition_periodic,
                        x_pos=boundary_condition_periodic,
-                       y_neg=boundary_condition_wall_noslip,
+                       y_neg=boundary_condition_noslip_wall,
                        y_pos=boundary_condition_couette,
                       )
 

examples/tree_2d_dgsem/elixir_lbm_lid_driven_cavity.jl

@@ -9,9 +9,9 @@ equations = LatticeBoltzmannEquations2D(Ma=0.1, Re=1000)
 
 initial_condition = initial_condition_lid_driven_cavity
 boundary_conditions = (
-                       x_neg=boundary_condition_wall_noslip,
-                       x_pos=boundary_condition_wall_noslip,
-                       y_neg=boundary_condition_wall_noslip,
+                       x_neg=boundary_condition_noslip_wall,
+                       x_pos=boundary_condition_noslip_wall,
+                       y_neg=boundary_condition_noslip_wall,
                        y_pos=boundary_condition_lid_driven_cavity,
                       )
 

examples/unstructured_2d_dgsem/elixir_ape_gauss_wall.jl

@@ -22,7 +22,6 @@ mesh = UnstructuredMesh2D(mesh_file)
 
 initial_condition = initial_condition_gauss_wall
 
-boundary_condition_slip_wall = BoundaryConditionWall(boundary_state_slip_wall)
 boundary_conditions = Dict( :OuterCircle  => boundary_condition_slip_wall,
                             :InnerCircle1 => boundary_condition_slip_wall,
                             :InnerCircle2 => boundary_condition_slip_wall,

examples/unstructured_2d_dgsem/elixir_euler_wall_bc.jl

@@ -27,7 +27,6 @@ end
 initial_condition = uniform_flow_state
 
 boundary_condition_uniform_flow = BoundaryConditionDirichlet(uniform_flow_state)
-boundary_condition_slip_wall = BoundaryConditionWall(boundary_state_slip_wall)
 boundary_conditions = Dict( :Bottom => boundary_condition_uniform_flow,
                             :Top    => boundary_condition_uniform_flow,
                             :Right  => boundary_condition_uniform_flow,

src/Trixi.jl

@@ -142,11 +142,10 @@ export initial_condition_constant,
 
 export boundary_condition_periodic,
        BoundaryConditionDirichlet,
-       boundary_condition_wall_noslip,
+       boundary_condition_noslip_wall,
+       boundary_condition_slip_wall,
        boundary_condition_wall,
-       boundary_condition_zero,
-       BoundaryConditionWall,
-       boundary_state_slip_wall
+       boundary_condition_zero
 
 export initial_condition_convergence_test, source_terms_convergence_test
 export initial_condition_harmonic_nonperiodic, source_terms_harmonic

src/equations/acoustic_perturbation_2d.jl

@@ -194,6 +194,7 @@ function initial_condition_gauss_wall(x, t, equations::AcousticPerturbationEquat
   return prim2cons(prim, equations)
 end
 
+
 """
     boundary_condition_wall(u_inner, orientation, direction, x, t, surface_flux_function,
                             equations::AcousticPerturbationEquations2D)
@@ -246,6 +247,7 @@ function initial_condition_monopole(x, t, equations::AcousticPerturbationEquatio
   return prim2cons(prim, equations)
 end
 
+
 """
   boundary_condition_monopole(u_inner, orientation, direction, x, t, surface_flux_function,
                               equations::AcousticPerturbationEquations2D)
@@ -282,6 +284,7 @@ function boundary_condition_monopole(u_inner, orientation, direction, x, t, surf
   return flux
 end
 
+
 """
     boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                             equations::AcousticPerturbationEquations2D)
@@ -304,6 +307,38 @@ function boundary_condition_zero(u_inner, orientation, direction, x, t, surface_
   return flux
 end
 
+
+"""
+    boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
+                                 equations::AcousticPerturbationEquations2D)
+
+Use an orthogonal projection of the perturbed velocities to zero out the normal velocity
+while retaining the possibility of a tangential velocity in the boundary state.
+Further details are available in the paper:
+- Marcus Bauer, Jürgen Dierke and Roland Ewert (2011)
+  Application of a discontinuous Galerkin method to discretize acoustic perturbation equations
+  [DOI: 10.2514/1.J050333](https://doi.org/10.2514/1.J050333)
+"""
+function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, x, t,
+                                      surface_flux_function, equations::AcousticPerturbationEquations2D)
+  # normalize the outward pointing direction
+  normal = normal_direction / norm(normal_direction)
+
+  # compute the normal perturbed velocity
+  u_normal = normal[1] * u_inner[1] + normal[2] * u_inner[2]
+
+  # create the "external" boundary solution state
+  u_boundary = SVector(u_inner[1] - 2.0 * u_normal * normal[1],
+                       u_inner[2] - 2.0 * u_normal * normal[2],
+                       u_inner[3], cons2mean(u_inner, equations)...)
+
+  # calculate the boundary flux
+  flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+
+  return flux
+end
+
+
 # Calculate 1D flux for a single point
 @inline function flux(u, orientation::Integer, equations::AcousticPerturbationEquations2D)
   v1_prime, v2_prime, p_prime_scaled = cons2state(u, equations)
@@ -385,34 +420,6 @@ end
 end
 
 
-"""
-    boundary_state_slip_wall(u_inner, normal_direction::AbstractVector,
-                             equations::AcousticPertubationEquations2D)
-
-Idea behind this boundary condition is to use an orthogonal projection of the perturbed velocities
-to zero out the normal velocity while retaining the possibility of a tangential velocity
-in the boundary state. Further details are available in the paper:
-- Marcus Bauer, Jürgen Dierke and Roland Ewert (2011)
-  Application of a discontinuous Galerkin method to discretize acoustic perturbation equations
-  [DOI: 10.2514/1.J050333](https://doi.org/10.2514/1.J050333)
-"""
-function boundary_state_slip_wall(u_inner, normal_direction::AbstractVector,
-                                  equations::AcousticPerturbationEquations2D)
-  # normalize the outward pointing direction
-  normal = normal_direction / norm(normal_direction)
-
-  # compute the normal and tangential components of the velocity
-  u_normal  = normal[1] * u_inner[1] + normal[2] * u_inner[2]
-  u_tangent = (u_inner[1] - u_normal * normal[1], u_inner[2] - u_normal * normal[2])
-
-  return SVector(u_tangent[1] - u_normal * normal[1],
-                 u_tangent[2] - u_normal * normal[2],
-                 u_inner[3],
-                 cons2mean(u_inner, equations)...)
-end
-
-
-
 @inline have_constant_speed(::AcousticPerturbationEquations2D) = Val(false)
 
 @inline function max_abs_speeds(u, equations::AcousticPerturbationEquations2D)

src/equations/compressible_euler_2d.jl

@@ -553,29 +553,106 @@ end
 
 
 """
-    boundary_state_slip_wall(u_internal, normal_direction::AbstractVector,
-                             equations::CompressibleEulerEquations2D)
+    boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
+                                 equations::CompressibleEulerEquations2D)
+
+Determine the boundary numerical surface flux for a slip wall condition.
+Imposes a zero normal velocity at the wall.
+Density is taken from the internal solution state and pressure is computed as an
+exact solution of a 1D Riemann problem. Further details about this boundary state
+are available in the paper:
+- J. J. W. van der Vegt and H. van der Ven (2002)
+  Slip flow boundary conditions in discontinuous Galerkin discretizations of
+  the Euler equations of gas dynamics
+  [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1)
+
+Details about the 1D pressure Riemann solution can be found in Section 6.3.3 of the book
+- Eleuterio F. Toro (2009)
+  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Pratical Introduction
+  3rd edition
+  [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
+
+Should be used together with [`UnstructuredMesh2D`](@ref).
 
-Determine the external solution value for a slip wall condition. Sets the normal
-velocity of the the exterior fictitious element to the negative of the internal value.
+!!! warning "Experimental code"
+    This wall function can change any time.
+"""
+function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, x, t,
+                                      surface_flux_function, equations::CompressibleEulerEquations2D)
+
+  norm_ = norm(normal_direction)
+  # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
+  normal = normal_direction / norm_
+
+  # rotate the internal solution state
+  u_local = rotate_to_x(u_inner, normal, equations)
+
+  # compute the primitive variables
+  rho_local, v_normal, v_tangent, p_local = cons2prim(u_local, equations)
+
+  # Get the solution of the pressure Riemann problem
+  # See Section 6.3.3 of
+  # Eleuterio F. Toro (2009)
+  # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Pratical Introduction
+  # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
+  if v_normal <= 0.0
+    sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed
+    p_star = p_local * (1.0 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2.0 * equations.gamma * equations.inv_gamma_minus_one)
+  else # v_normal > 0.0
+    A = 2.0 / ((equations.gamma + 1) * rho_local)
+    B = p_local * (equations.gamma - 1) / (equations.gamma + 1)
+    p_star = p_local + 0.5 * v_normal / A * (v_normal + sqrt(v_normal^2 + 4.0 * A * (p_local + B)))
+  end
+
+  # For the slip wall we directly set the flux as the normal velocity is zero
+  return SVector(zero(eltype(u_inner)),
+                 p_star * normal[1],
+                 p_star * normal[2],
+                 zero(eltype(u_inner))) * norm_
+end
+
+"""
+    boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
+                                 surface_flux_function, equations::CompressibleEulerEquations2D)
+Should be used together with [`TreeMesh`](@ref).
 
 !!! warning "Experimental code"
     This wall function can change any time.
 """
-@inline function boundary_state_slip_wall(u_internal, normal_direction::AbstractVector,
-                                          equations::CompressibleEulerEquations2D)
+function boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
+                                      surface_flux_function, equations::CompressibleEulerEquations2D)
+  # get the appropriate normal vector from the orientation
+  if orientation == 1
+    normal = SVector(1, 0)
+  else # orientation == 2
+    normal = SVector(0, 1)
+  end
 
-  # normalize the outward pointing direction
-  normal = normal_direction / norm(normal_direction)
+  # compute and return the flux using `boundary_condition_slip_wall` routine above
+  return boundary_condition_slip_wall(u_inner, normal, x, t, surface_flux_function, equations)
+end
+
+"""
+    boundary_condition_slip_wall(u_inner, normal_direction, direction, x, t,
+                                 surface_flux_function, equations::CompressibleEulerEquations2D)
+Should be used together with [`StructuredMesh`](@ref).
 
-  # compute the normal and tangential components of the velocity
-  u_normal  = normal[1] * u_internal[2] + normal[2] * u_internal[3]
-  u_tangent = (u_internal[2] - u_normal * normal[1], u_internal[3] - u_normal * normal[2])
+!!! warning "Experimental code"
+    This wall function can change any time.
+"""
+function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, direction, x, t,
+                                      surface_flux_function, equations::CompressibleEulerEquations2D)
+  # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back
+  # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh
+  if isodd(direction)
+    boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction,
+                                                  x, t, surface_flux_function, equations)
+  else
+    boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction,
+                                                 x, t, surface_flux_function, equations)
+  end
 
-  return SVector(u_internal[1],
-                 u_tangent[1] - u_normal * normal[1],
-                 u_tangent[2] - u_normal * normal[2],
-                 u_internal[4])
+  return boundary_flux
 end