Issue with slip wall BCs for CompressibleEulerEquations2D using TreeMesh

[jlchan]

I'm observing an issue with slip wall BCs for TreeMesh. Here's a MWE:

using OrdinaryDiffEq
using Trixi

equations = CompressibleEulerEquations2D(1.4)
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# Create a uniformly refined mesh with semi-periodic boundaries
coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                periodicity=(true, false),
                n_cells_max=30_000) # set maximum capacity of tree data structure

function initial_condition_gaussian(x, t, equations)
  rho = 1.0 + .1 * exp(-25 * (x[1]^2 + x[2]^2))
  v1 = 0.0
  v2 = 0.0
  p = rho^equations.gamma
  return prim2cons(SVector(rho, v1, v2, p), equations)
end

initial_condition = initial_condition_gaussian

# define inviscid boundary conditions
boundary_conditions = (; x_neg = boundary_condition_periodic,
                         x_pos = boundary_condition_periodic,
                         y_neg = boundary_condition_slip_wall,
                         y_pos = boundary_condition_slip_wall)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.5
tspan = (0.0, .50)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
callbacks = CallbackSet(summary_callback, alive_callback)

###############################################################################
# run the simulation

time_int_tol = 1e-8
sol = solve(ode, RDPK3SpFSAL49(), dt = 1e-5, abstol=time_int_tol, reltol=time_int_tol,
            save_everystep=false, callback=callbacks)
summary_callback() # print the timer summary

The code crashes shortly after time 3.8826e-01. If I run the code up to this point, I get the following result

If I replace the slip wall BCs with periodic, the elixir works fine.

[ranocha]

I can confirm that it crashes for me, too. Everything is fine if I switch the mesh to

mesh = StructuredMesh((16, 16), coordinates_min, coordinates_max; periodicity=(true, false))

[jlchan]

Aha, StructuredMesh flips the sign of the flux

Trixi.jl/src/equations/compressible_euler_2d.jl

Lines 353 to 368 in 5ae580b

	@inline 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 boundary_flux
	end

while the TreeMesh implementation does not

Trixi.jl/src/equations/compressible_euler_2d.jl

Lines 332 to 345 in 5ae580b

	@inline 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
	# 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

I'd guess this is the root of the issue (and the reason why TreeMesh isn't working for #1165)