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

[andrewwinters5000]

This addresses #765 . Implements a treatment of a compressible Euler slip wall boundary condition that zeros out the normal velocity component and computes the pressure as the exact solution of a 1D Riemann problem.

NOTE: Due to the current implementation of the generic BoundaryConditionWall this new boundary treatment returns an external boundary state that will then cancel the normal velocity contribution and use the Riemann solution pressure p_star in a call to the surface_flux_function. This is slightly different compared to how FLUXO or the paper from van der Vegt does this. They avoid an additional call to the surface flux by directly setting the boundary flux to be something like

flux = SVector(0.0, p_star * normal_vector, 0.0)

UPDATE: The new BoundaryConditionFlux formulation that this PR uses mimics the FLUXO strategy

[sloede]

Thanks for adding this feature! Does it make sense to change the signature for BCWall type and instead do the "special" treatment for the no-slip wall? Would that be conceptually simpler/faster?

Also, is a similar change useful/necessary for the no-slip condition?

[andrewwinters5000]

Does it make sense to change the signature for BCWall type and instead do the "special" treatment for the no-slip wall?

I am not sure TBH. I struggled with this a bit because the acoustic perturbation equations uses this BCWall type with the current signature. I thought about the possibility of having several versions of function (boundary_condition::BoundaryConditionWall) that dispatches on the particular equation system, but I decided against it as it seemed too messy.

Would that be conceptually simpler/faster?

It would probably be faster because it would just set the flux value at the boundary instead of spending the effort to do another, e.g., flux_hll. Conceptually it is pretty easy to derive on paper why the way it is currently implemented in Trixi is equivalent to directly setting the flux a la van der Vegt (it takes about a page of calculations)

[ranocha]

I am not sure TBH. I struggled with this a bit because the acoustic perturbation equations uses this BCWall type with the current signature. I thought about the possibility of having several versions of function (boundary_condition::BoundaryConditionWall) that dispatches on the particular equation system, but I decided against it as it seemed too messy.

Then it's good that we didn't declare this one as stable. From my point of view, we should

• change the name of BoundaryConditionWall since it doesn't work as we would like
• implement the new behavior of a corrected wall BC as done here
• re-create the old behavior for the acoustic perturbation equations (can we get it using the same signature as done here, i.e., by setting the boundary flux directly?)

[andrewwinters5000]

Also, is a similar change useful/necessary for the no-slip condition?

I missed this earlier. I think it would be useful for no-slip. I am only speculating with how FLUXO does it, but basically we could get the advective Euler part using this new way and then you add in the diffusive parts (at least for an isothermal wall)

[andrewwinters5000]

• change the name of BoundaryConditionWall since it doesn't work as we would like

Do you mean to be something more specific like BoundaryConditionSlipWall? Then the function boundary_state_slip_wall would be renamed to something like pressure_state_slip_wall because it is the only "important" quantity in this case.

• implement the new behavior of a corrected wall BC as done here

So the altered and newly named BoundaryConditionSlipWall would call this pressure state function and then set the flux directly?

• re-create the old behavior for the acoustic perturbation equations (can we get it using the same signature as done here, i.e., by setting the boundary flux directly?)

I believe this will be straightforward. I am re-reading the reference I found for this and it looks like a similar trick can be used to avoid the additional call to a surface flux for the APE slip wall

[sloede]

• change the name of BoundaryConditionWall since it doesn't work as we would like

Do you mean to be something more specific like BoundaryConditionSlipWall? Then the function boundary_state_slip_wall would be renamed to something like pressure_state_slip_wall because it is the only "important" quantity in this case.

Out of curiosity: Is it in this case still useful to have BoundaryConditionSlipWall at all, or could one just implement a boundary_condition_slip_wall that does Everything Right?

[andrewwinters5000]

Out of curiosity: Is it in this case still useful to have BoundaryConditionSlipWall at all, or could one just implement a boundary_condition_slip_wall that does Everything Right?

Would this mean then that boundary_condition_slip_wall would fully live, e.g., in compressible_euler_2d.jl and then the BoundaryCondition functionality from equations.jl is only there to separate between orientations versus normals or the different mesh types?

[sloede]

Out of curiosity: Is it in this case still useful to have BoundaryConditionSlipWall at all, or could one just implement a boundary_condition_slip_wall that does Everything Right?

Would this mean then that boundary_condition_slip_wall would fully live, e.g., in compressible_euler_2d.jl and then the BoundaryCondition functionality from equations.jl is only there to separate between orientations versus normals or the different mesh types?

I don't know :-) I got confused by the question regarding which functionality now lives in which part of the implementation.

Let me try to rephrase my point: Does the BoundaryConditionX type still bring enough capabilities/flexibility to the table to justify its existence (and the added complexity in understanding the code)? Or could we - with only small redundancies - just define the wall BCs as we would define other BCs (just a single function), which might not be as elegant but easier to grasp for new users?

[codecov]

Codecov Report

Merging #830 (6921de2) into main (b239e09) will increase coverage by 0.02%.
The diff coverage is 95.74%.

@@            Coverage Diff             @@
##             main     #830      +/-   ##
==========================================
+ Coverage   92.84%   92.87%   +0.02%     
==========================================
  Files         189      189              
  Lines       17757    17775      +18     
==========================================
+ Hits        16486    16507      +21     
+ Misses       1271     1268       -3     

Flag	Coverage Δ
unittests	92.87% <95.74%> (+0.02%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
...lti_2d/elixir_euler_rayleigh_taylor_instability.jl	88.24% <ø> (ø)
examples/p4est_3d_dgsem/elixir_euler_ec.jl	100.00% <ø> (ø)
examples/tree_2d_dgsem/elixir_lbm_couette.jl	80.00% <ø> (ø)
...ples/unstructured_2d_dgsem/elixir_euler_wall_bc.jl	90.00% <ø> (ø)
src/Trixi.jl	66.67% <ø> (ø)
src/equations/equations.jl	92.68% <ø> (+2.30%)	⬆️
src/equations/compressible_euler_2d.jl	94.17% <90.91%> (-<0.01%)	⬇️
..._dgsem/elixir_euler_rayleigh_taylor_instability.jl	87.50% <100.00%> (ø)
src/equations/acoustic_perturbation_2d.jl	89.94% <100.00%> (+0.06%)	⬆️
src/equations/compressible_euler_3d.jl	95.20% <100.00%> (+0.25%)	⬆️
... and 3 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b239e09...6921de2. Read the comment docs.

[andrewwinters5000]

Does the BoundaryConditionX type still bring enough capabilities/flexibility to the table to justify its existence

I will keep this type around for the time being as I prototype a new way of doing this slip wall BC. We can always streamline this once it works

[sloede]

I will keep this type around for the time being as I prototype a new way of doing this slip wall BC. We can always streamline this once it works

Sounds good!

[ranocha]

That's a good plan - let's prototype what we need ("make it work"). Then, we can see how general this is and decide how to adapt it if necessary ("make it nice"). Finally, we should check whether we broke performance ("make it fast").

[andrewwinters5000]

I put in an initial proposition of how to set the boundary surface flux directly using the pressure values and the normal vector for the 2d UnstructuredMesh2D testcase elixir_euler_wall_bc.jl. I am not confident that it is easily extendable to have this work for compressible Euler and APE equations due to the need of setting p_star * normal in the velocity flux components and zeros everywhere else. It might be easier to just remove the BoundaryConditionWall type and implement the necessary routines equation by equation.

[jlchan]

The assumption on the existence of a scalar pressure state as well as the positions of the velocity definitely makes the current approach with pressure_state_slip_wall harder to generalize to APE or MHD.

I don't understand the issues with the old BoundaryConditionWall implementation where the flux is evaluated at the boundary state. It's more extensible - is the problem that it's inefficient? If so, we could treat it as a fallback and specialize either BoundaryConditionWall or boundary_condition_slip_wall for specific equations. We wouldn't have to do all of them in this PR - just specializing 2D Euler would provide a proof of concept, and we could leave the rest for future PRs.

The changes to boundary_state_slip_wall for the APE equations looks good (computing the normal projection without needing to compute tangentials).

[ranocha]

I don't understand the issues with the old BoundaryConditionWall implementation where the flux is evaluated at the boundary state.

The issue is that it's not doing that. Instead, it computes

flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)

In general, this is different from what we would like to have. That's basically the gist of what @andrewwinters5000 said above, isn't it?

NOTE: Due to the current implementation of the generic BoundaryConditionWall this new boundary treatment returns an external boundary state that will then cancel the normal velocity contribution and use the Riemann solution pressure p_star in a call to the surface_flux_function. This is slightly different compared to how FLUXO or the paper from van der Vegt does this. They avoid an additional call to the surface flux by directly setting the boundary flux to be something like

flux = SVector(0.0, p_star * normal_vector, 0.0)

but something like this (that is so specialized to this particular boundary case) is not really tractable in the current BC ecosystem of Trixi.

[ranocha]

It might be easier to just remove the BoundaryConditionWall type and implement the necessary routines equation by equation.

That's completely correct, the current setup doesn't scale. What I had in mind was more a high-level interface such as

# Flux-based boundary condition for use with unstructured meshes such as
# `UnstructuredMesh2D` and `P4estMesh2D`.
# Note: For unstructured meshes, we lose the concept of an "absolute direction"
#       that is used for Cartesian and structured meshes.
@inline function (boundary_condition::BoundaryConditionFlux)(u_inner,
                                                             normal_direction::AbstractVector,
                                                             x, t,
                                                             surface_flux_function, equations)
  return boundary_condition.boundary_flux_function(u_boundary, normal_direction, equations)
end

with equation-specific implementation of boundary_flux_function.

[sloede]

I've left some notes but nothing major. Overall, this LGTM. While I do not really like the "smart" reuse of wall BC functionality for the different mesh types, I also do not have a better suggestion except to maybe align the normal orientations between the structured and unstructured meshes (but that's a bigger issue by itself).

[sloede]

LGTM! Once tests pass & coverage is OK, this can be merged.

[andrewwinters5000]

LGTM! Once tests pass & coverage is OK, this can be merged.

Okay, but I am expecting the Rayleigh-Taylor structured mesh test to fail. I still cannot work out why the CI is giving slightly different numbers than Roci or my local machine.

[jlchan]

Okay, but I am expecting the Rayleigh-Taylor structured mesh test to fail. I still cannot work out why the CI is giving slightly different numbers than Roci or my local machine.

Does this only happen for structured mesh RTI but not for DGMulti?

Edit: just looked at CI - apparently so.

[andrewwinters5000]

Does this only happen for structured mesh RTI but not for DGMulti?

Yes the CI only fails for the structured mesh RTI, the DGMulti version passes. Everything passes for structured/unstructured/DGmulti on my local machine.

[jlchan]

@andrewwinters5000 that is really odd. Would you mind trying CI (for the RTI elixir on structured meshes) using the commented out BC setup? This avoids a slip wall BC altogether.

If this modified test passes, we could create an issue on the weird CI behavior for RTI under slip walls to move it out of this PR.

[andrewwinters5000]

@andrewwinters5000 that is really odd. Would you mind trying CI (for the RTI elixir on structured meshes) using the commented out BC setup? This avoids a slip wall BC altogether.

If this modified test passes, we could create an issue on the weird CI behavior for RTI under slip walls to move it out of this PR.

I tried this swap and the new test passes locally for me. We will see what the CI says. (for now I disabled auto-merge such that we can discuss the plan going forward)

[andrewwinters5000]

If this modified test passes, we could create an issue on the weird CI behavior for RTI under slip walls to move it out of this PR.

Okay, so this modification passes the CI. Any ideas what to do next @sloede , @ranocha ?

[ranocha]

I reduced the resolution of the StructuredMesh test case to match the one of the DGMulti test since it was very expensive (ca. 500 seconds in CI). Let's see whether this looks better.

[andrewwinters5000]

I reduced the resolution of the StructuredMesh test case to match the one of the DGMulti test since it was very expensive (ca. 500 seconds in CI). Let's see whether this looks better.

Sounds good. Do you think that using a coarser resolution might help with the weird test value mismatch between local and CI?

[ranocha]

Do you think that using a coarser resolution might help with the weird test value mismatch between local and CI?

Yes. We will use far fewer time steps and floating point operations in general, which will reduce chances for floating point differences to add up.

[ranocha]

CI passes 🎉