Correct gradients at boundaries

[GiuSirianni]

Is your feature request related to a problem?

At any boundary that is not a Neumann-zero boundary ($\nabla \vec{V} \cdot \hat{n} = 0$) the gradients computed are incorrect, as they are all computed assuming that at boundaries $\vec{q}_i = \vec{q}_j$.

For example, for a Riemann boundary condition we have a certain imposed external state $\vec{q}_e$ but this is not considered in the computation of gradients.

In my experience, this issue is exacerbated when using MUSCL reconstruction at corners shared by walls and Riemann boundaries, particularly in NICFD cases, as it can cause divergence of the solver and/or unphysical solutions.

Describe the solution you'd like

The solution is not easy as ideally it would require storing all "external" states for the computation of gradients. One may need to choose if swapping the order of the computation of boundary residuals and of gradients to avoid double computations, but it may be inefficient for MPI.

Describe alternatives you've considered

One could also disable MUSCL at corner/boundary nodes which would require saving which points are corners. This is a 1st order approximation locally, but my opinion is that the current assumption that $\vec{q}_i = \vec{q}_j$ is also lower order.

Another possibility is correcting boundaries after they have been computed by summation/subtraction, before the upwind residual computations.

Additional context

[bigfooted]

Thanks for creating an issue for this.
The correct computation of gradients for symmetry planes and Euler walls is in PR #2194

For no-slip walls, we know that velocity gradients tangential to the wall are zero, this could be taken into account, or in general as you said we should take into account the gradient information from the type of boundary condition.

[bigfooted]

Hi @GiuSirianni, considering your experience with this, it would be great if you can contribute to this improvement. We have weekly developer meetings, do you have time to talk about this perhaps next week Wednesday? I think it is an important issue that needs to be addressed.

[GiuSirianni]

Hi @GiuSirianni, considering your experience with this, it would be great if you can contribute to this improvement. We have weekly developer meetings, do you have time to talk about this perhaps next week Wednesday? I think it is an important issue that needs to be addressed.

Sure, I can join! I'm not sure if/when I will be able to actively contribute to a fix (I'm not even sure what it would be)

[GiuSirianni]

To add some data I'll put a test where you can see the issue:

This is an inviscid Euler simulation using the Span-Wagner EoS (CoolProp) of a non-ideal MDM nozzle, see https://su2code.github.io/tutorials/NICFD_nozzle using GREEN_GAUSS for the gradient computation. Convergence can still be achieved but the solution at the inlet corner is completely wrong, while the outlet is ok most likely because characteristics are outgoing as it is supersonic. The artifacts disappear if we disable MUSCL everywhere (1st order solution) or disable it only on boundaries (not ideal solution, but disabling only on corners would still be good enough for now I believe).

Using WEIGHTED_LEAST_SQUARES seems to not present the same issue, in this test case at least, as the stencil "does not care" about the boundary states.

The boundary conditions are:

• MARKER_SYMMETRY at centerline
• MARKER_EULER at wall
• MARKER_RIEMANN= (INLET, TOTAL_CONDITIONS_PT, 904388, 542.13, 1.0, 0.0, 0.0) at inlet
• MARKER_RIEMANN= (OUTLET, STATIC_PRESSURE, 200000.0, 0.0, 0.0, 0.0, 0.0) at outlet

I tried both with and without a slope limiter as there are no discontinuities, but it makes no difference on the artifacts:

SLOPE_LIMITER_FLOW= VENKATAKRISHNAN_WANG
VENKAT_LIMITER_COEFF= 0.1

Complete test case ZIP

mdm_coolprop_nozzle.zip

Inlet pressure zoom

Notation:

• "1st order" no MUSCL
• "2nd order" MUSCL as implemented in SU2
• "2nd order (BC 1st order)" I simply disabled MUSCL on ALL physical boundaries in the upwind residual computations, see code snippet below

Proof of concept code modification for "2nd order (BC 1st order)":

To show that the error lies in MUSCL/gradients at boundaries I added these two lines of code in the upwind gradient computation

Residuals:

Mesh:

[pcarruscag]

Have you tried reconstructing from the interior point but not the one on the boundary?
For example, by setting the limiters of the interior point to 0.

[GiuSirianni]

@pcarruscag I'm not sure if I misunderstood, you mean setting the limiter to zero on the boundary point right?
I tried as shown below, you can still see the first order solution at boundaries is visible, but less. As said before, using weighted-least-squares does not present this issue as the stencil doesn't take the boundary into consideration of course.

Disable MUSCL on all points involved with a boundary

Disable MUSCL only on boundary points

[pcarruscag]

Thanks for testing. Both gradient methods use all the neighboring points. Green Gauss has special treatment for boundary surfaces and it must be particularly less accurate than least squares at these points.
With that said the way we handle intersecting boundaries is very important. For this case https://su2code.github.io/vandv/swbli/ I had to add an Euler wall after the inlet to prevent the kind of issues you are seeing.
My concern is that making things locally first order is just masking the problem with numerical dissipation.

[GiuSirianni]

I agree that the "proposed" solution is no solution at all, it is just to prove the source of (a) problem. I do not agree that it is just masking the issue, it is removing it and just introducing another. Unless there is something I'm missing on the treatment of corners, which is very likely.

I think the real solution in a perfect world would be to have ghost nodes with the boundary state saved and updated appropriately, as it would also fix symmetry boundaries, etc. I understand this is borderline non-implementable.

[bigfooted]

PR #2194 Seems to fix this specific issue because a symmetry plane and an Euler wall are used here. We can think about similar improvements for viscous walls etc.

[GiuSirianni]

That's great! @bigfooted can I ask if you see any issue in allowing the velocity index iVel used in the gradient computation to be more than one scalar value? From scouring the code in the PR I think it should not be a problem. I am asking as I am using SU2 as a playground for a full-disequilbrium two-phase solver so I effectively have 2 velocity vectors in my solution.

[bigfooted]

The current SU2 implementation checks if the array of variables contains the velocity and gives the index to the gradient calculator. You can give the gradient calculator the starting index and ending index as well.. So suppose you have an array of M variables, you could split it into 2 arrays, each containing a velocity vector, so array [0,..,N] and array [N+1,..,M], and you call greengauss twice.

[GiuSirianni]

Good to know, thank you!

[bigfooted]

@GiuSirianni @pcarruscag I think the correction of the symmetry plane and Euler wall has (mostly?) fixed this issue. But is there anything left in this discussion that we can use to improve quality and/or convergence?
After re-reading this issue, I conclude that we do not want to switch from second order to first order on boundaries?

Can this issue be closed?

[GiuSirianni]

No, switching to first order is not a proper solution so I agree on avoiding doing this.

The only possible thing I would add is to take into account the inlet and outlet boundary conditions in the computation of gradients at boundary points, but this could be a separate issue altogether.

[pcarruscag]

There are other boundary-variable combinations where the gradients are known and can be enforced.