Inconsistent results with the Symmetric Interior Penalty method for the Navier-Stokes equation

Summary

The Interior Penalty (IP) method was recently fixed in the MR !1303 (merged). However, the Symmetric IP (SIP) method still presents inconsistent results for three-dimensional simulations (geometries extruded in the third direction). Although the SIP method does not fail in the CI/CD system for two-dimensional test cases, there is a bug in the three-dimensional implementation. After finding this bug, in a second (MR !1377 (merged)), we made the Incomplete IP (IIP) method the default method, so that we switched from the SIP to the IIP method.

The main motivation to fix this bug is that the SIP method is numerically more stable compared to the IIP method, since it includes the symmetric term in the numerical formulation which makes it symmetric and adjoint-consistent, thus improving convergence properties.

Steps to reproduce

The easiest way to reproduce this issue is to set up a flat plate case with period or symmetric boundary condition in the z-direction. The code presented the same bug with both boundary conditions. So, in order to create the geometry for the flat plate case, the user will need to build a two-dimensional mesh and then extrude it in the third direction. In addition, the Navier-Stokes equations must be used for the numerical simulation. Also, to set up the SIP method, the user has to add the parameter <P> IPSymmFluxCoeff = 1.0 </P> in the XML session file.

What is the current bug behavior?

We observed that high frequency oscillations are created in the third direction when the symmetric term of the interior penalty method is used (SIP method) (Figure below). The figure shows the velocity in the z-direction. These high frequency oscillations were observed in a two-dimensional geometry which was extruded in the z-direction. Note that, in the figure below, the numerical simulation was performed using P=3 and was used three planes in the z-direction for the spatial discretisation. It is easy to see that the same pattern repeats for each plane, which is the interface between elements in the z-direction.

In addition, this issue only appears when the Navier-Stokes equation is used for the numerical simulation. When the Euler equation is employed this issue does not appear, which makes me believe that this bug is in the diffusion operator, more precisely in the class DiffusionIP.

What is the expected correct behavior?

The correct behavior for this numerical simulation is presented in the figure below. In this figure, we see that there is no velocity in the z-direction, since the flow is symmetric with respect to the z-direction. It should be mentioned that this numerical simulation was performed with the IIP method.

Tests

At the moment, the CI/CD system is not able to check this issue, so a three-dimensional small test case should be added to the CI/CD system in order to verify that the SIP method works properly for the Navier-Stokes equation.