Improve the Euler wall Jacobian

[pcarruscag]

#2485

[bigfooted]

I've updated the serial regression tests so we can see what is going on. It looks like for most cases, there was a small improvement in the residual values. Usually a good sign.
The only case that changed residuals a lot was the rotating naca case (line 744).

[YairMO]

By the way, does the SU2 code use weak or strong boundary conditions?

[pcarruscag]

It varies, I think Euler walls can be seen as strongly imposed now, velocity on no-slip walls is imposed strongly but temperature on isothermal boundaries is imposed weakly (via a flux), most other compressible boundaries are imposed weakly, for the incompressible solver it is possible to select the treatment for some boundaries.

If you have recommendations from experience with other vertex-centered codes, we would certainly like to hear them.
Cell-centered makes many things easier.

[YairMO]

Thank you for clarifying this. My experience is with cell-centered solvers, and BC's implicit treatment is easier.
Nevertheless, I can understand how it may become more complex with vertex-centred solvers. I may ask a colleague who has his own vertex-centred code how it treats the BCs implicitly. Moreover, if it is acceptable, I'm willing to prepare a short note on how to implicitly implement BC in a vertex-centred solver?
Moreover, using weak boundary conditions makes it much simpler to treat the BCs implicitly; it becomes similar to a cell-centred approach.

[pcarruscag]

I think that would be valuable, thank you. I think for engineering applications the most important boundary is no slip wall, for the longest time I've wanted to try a weak no slip condition to see if it is more robust than the current approach of strongly imposing 0 velocity. It would be great if your colleague can share any notes with either approach.

[YairMO]

Hi,
I've prepared a short note for possible implicit treatment at the boundaries, assuming the use of strong boundary conditions:
Implicit_BC.pdf
Your comments (corrections) are most welcome.
I admit I'm a "person" of cell-centered code, and I may be all wrong here.

[pcarruscag]

Hi Yair,
Thank you for your email and for taking the time to write the document.
I think the broad strokes make sense, the treatment you propose is similar to some of the boundaries we have where we construct the boundary state based on the characteristics dictated by the boundary condition and characteristics extrapolated from the interior of the domain. Our farfield boundary is perhaps where this is done more rigorously.
In your example pressure is a characteristic extrapolated from the interior, and velocities are dictated by the boundary condition.
Regarding the implicit treatment, in some boundaries we do not consider the dependence of this boundary state on the information extrapolated from the interior (as you suggest in the document).
This is certainly worth exploring.

Now, there are some important details.
If we want to set pressure at walls such that the normal gradient is strongly imposed as 0, we should also consider the effect of other boundary points in this pressure gradient, in your document you don't include these connections but in most cases they are required otherwise the gradient operator is ill conditioned, consider that with prism boundary layer each boundary point only has one interior point.
If we consider boundary points it becomes a system of equations for boundary pressures, instead of a simpler extrapolation from the interior.
If we solved for primitive variables we could modify the Jacobian matrix to perform this. However, we solve for conservative variables and thus I don't think it is easy to modify the Jacobian to reflect dP/dn = 0.

Our current treatment for no-slip walls involves setting the velocity equal to the grid velocity (or the momentum = rho * v_grid) strongly, and doing nothing special for the mass and energy fluxes (when there is a Neuman condition on temperature).
I'm interested in an alternative where we define a boundary state (v = v_grid, P and T extrapolated based on d/dn, for testing purposes the extrapolation can be very simple if the interior point is perpendicular to the boundary) and then evaluate the flux between this state and the boundary point. We then need to pick some components of this flux, for example, there will be an artificial mass flux that we know should be zero and we probably want to impose strongly to avoid unphysical results. However some of these fluxes may be useful to increase robustness (which is a challenge for SU2) maybe some additional dissipation in the momentum equations for example.
Then we can include the Jacobian of this flux between the boundary state and boundary point in the way you suggest.

This would be a treatment similar to cell-centered codes, where values are set on boundary faces. So in what I'm thinking we don't need to impose values (e.g. velocity) on the boundary points, if they converge to slightly different values than the boundary states it is not a problem, they would just behave like the first CVs in cell centered codes.
Computing fluxes at boundaries that become part of the global system of conservation equations is the only way (that I know) to preserve the conservation properties of the finite volume method.

[YairMO]

Hi,

Thank you for the comments.

I understand that in most cases, the grid is generated with one grid point connected to one interior point, as in the near-wall prism layers. Hence, the solution is extrapolated from a single interior point if required. For example, extrapolating the pressure at a stationary no-slip BC. Moreover, extrapolating the pressure to the BC point at a perfectly orthogonal grid automatically satisfies the zero pressure gradient.
I understand that the SU2 code solves the conservative variables, but the zero gradient condition is on the pressure. In such a case, one obtains this condition and calculates the correct conservative variables.

My schematic illustration of the possible grid was to give a general idea, and it is less common in near-wall BC. However, such a grid may happen in other BCs, such as symmetry plans and far fields. Some of the interior information is extrapolated to these BCs and should be reflected in the implicit treatment suggested in the note.

If I understood correctly, you suggested testing weak BCs in the SU2. In such a case, one doesn't impose the BC at the point but solves this point and imposes the BCs at the faces on the boundaries. This will be similar to the cell-centered concept. It also circumvents the situation of a given point that is located on different BC surfaces. However, (although I don't have much experience with this), it always bothered me how the source term is treated in such a case. For example, a singular term happens in the axisymmetric case.

In my opinion, it is worth making the effort to consider the fully implicit treatment at the current BSc implementation, especially for solid walls (no-slip and impermeable) and symmetry planes.

[YairMO]

I've slightly modified the note (again), and I hope it will clarify the main idea better.
Implicit_BC-Mod2.pdf

[pcarruscag]

I understand, I don't think it would be a weak BC because we would still enforce what is known about the flux, for example, 0 mass flux at walls. A weak treatment would not enforce the mass flow, letting it converge to 0 as the grid is refined and the solution in the boundary point approaches the imposed boundary state.
The implicit treatment is orthogonal to the boundary treatment, we can have weak BCs with accurate Jacobians, strong BCs with approximate Jacobians.
In my experience all these things need to be tested because there is a balancing act between having an accurate Jacobian that becomes too ill-conditioned for linear solvers and an approximate Jacobian that may not capture enough dependencies between variables to be effective. Both situations can force a low CFL and poor convergence.

[YairMO]

I'm not fully following up on why using strong BCs would result in a non-accurate Jacobian. Could you please elaborate more on this? Naively, I would consider the formulation given in the note to be an accurate Jacobian (first-order accurate).

[pcarruscag]

I'm not saying it would result in a non-accurate Jacobian, but rather that the consistency of the Jacobian with the flux or boundary treatment is a choice.
Nothing says the Jacobian must be consistent, and sometimes the best convergence is not obtained with a consistent (accurate) Jacobian.
Especially because it is not trivial to include second-order terms in the Jacobian and thus the full Jacobian matrix is always approximate in some sense.
Even if we had a fully accurate Jacobian (which implies infinite CFL) and an incredibly powerful linear solver that could handle such an ill-conditioned matrix, we would very likely still need some relaxation of the Newton updates.
In our case, we are marching the second-order system to steady state using a first-order Jacobian and finite CFL, therefore there is no way to determine a priori what Jacobian will give us the best overall performance (time to solution).

In my experience the best performing Jacobian depends on the application (flow regime, primal or adjoint) and the different options need to be tested.

[YairMO]

Agreed. All I'm saying is that in order to understand if the BCs are treated appropriately, the first-order Jacobian and the first-order flux should be able to obtain a Newton method. This is just a way to check the implementation, and not a way for actual computations. In actual computations, we use second-order flux with first-order Jacobian, which, of course, can not achieve Newton iterations and has its limitations, as you mentioned. I don't suggest going for the second-order Jacobian, either.

[rois1995]

I have tried the changes in this PR on the SWBLI case. With previous implementation it was diverging after a few iterations (NaN in the residual), whereas now I can obtain deep residual convergence. However, at some point I get an oscillatory behavior of the residuals. Plus, I still cannot achieve deep convergence of the turbulence equations, as the residuals stagnate at certain values. I have tried with the combination WLS-LS instead of GG only for the gradient computation and it gets worse.

Here the files to replicate SWBLI.zip .

[rois1995]

I have tested it one the problem in #2390 and I still don't get zero-gradients on the symmetry planes. I have zero-gradients for the Pressure but not for any other variable. And I still don't know if it is an expected behavior or not.

[YairMO]

I would check how forcing positivity (this is missing in the code) of the reconstructed (the 2nd-order reconstruction) relevant primitive values (density and pressure) about the convergence

[YairMO]

In most cases the grid that are generated for high-Reynolds viscous flows, characterize by prismatic elements near the wall surface. Hence a single surface point point connected to a single interior point. However, there are situations where this is not the case, for example at singular geometry, such as in sharp trailing edge, where a single surface grid point connected to several interior points. How the pressure on this singular surface grid point is specified?