Refactoring of convective numerics classes

[pcarruscag]

Proposed Changes

Hi Folks,
I would like your views on something. Many of the convective schemes we have are small variations of each other, and so I was thinking we can maybe have some intermediate numerics classes to make things a bit easier to maintain, quicker to compile etc.
So far I experimented with having a base class for central schemes, the ComputeResidual method is moved to this class and the children classes (Lax, JST, and JST_KE) only implement the artificial dissipation part specific to them.
Similarly AUSM+up(1/2) and SLAU(1/2) only differ in how the mass and pressure fluxes are computed so again a common base for these 4 could be devised. Isolating the computation of mass and pressure fluxes could be an interesting way to compute the Jacobians of these schemes in a hybrid way (currently the Roe Jacobian is used).
There is a very small performance penalty (maybe 1-2% for explicit solution methods) but my reason to be looking at the numerics in the first place is that with some small tweaks to the Jacobians we can probably run most schemes at higher CFL (namely increase the dissipative part of the Jacobian to make the system matrix more diagonal dominant) so for implicit solution methods there would be a speedup.
So, what do you think?
Cheers,
Pedro

PR Checklist

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags).
• My contribution is commented and consistent with SU2 style.
• I have added a test case that demonstrates my contribution, if necessary.

[WallyMaier]

Hi Pedro,
Thanks for this! I think eliminating duplicated code is a good idea. That being said, would this make adding new numerical schemes more difficult? (Thinking from a new developer standpoint).

[pcarruscag]

Hi Wally,
I think that would depend, if the new scheme is part of a family that already exists it should be easier to implement, if not you do it from scratch... There is also the possibility of making life hard by trying to find common ground where there is none, e.g. the way CUSP is implemented in SU2 has something in common with the central schemes, but it did not feel right to attach an upwind scheme to the central family.

[aeroamit]

Hi Pedro,

It is indeed a good idea to keep similar schemes (with minor variations) under one umbrella, especially central scheme as only the dissipation term calculation differs.

• In case of AUSM+-up1/2 and SLAU1/2, each can be reduced by having a separate pressure flux definition.

• But if we try to combine AUSM+-up and SLAU, will it be a clear/helpful implementation as mass flux definitions/calculations (which is substantial portion) are different (unlike the central scheme, where only dissipation term varies)?

• Also could you please shed some light on Jacobian part for these schemes (as you mentioned - “...Isolating the computation of mass and pressure fluxes could be an interesting way to compute the Jacobians of these schemes in a hybrid way ...”).

Best
Amit

[pcarruscag]

Hi Amit,
The idea is basically this:

F is the general form of the flux per unit area for those 4 schemes (according to the AUSM+-up2 paper).
The Jacobians of PSI with respect to the conservatives (U) are known analytically, the derivatives of the mass and pressure fluxes are harder to derive, so if we encapsulated their calculation in a function ( [m,p] = f(U) ) we could approximate them by finite differences (or even use AD).
So the base class would deal with the general form and Jacobians, the derived classes would implement this function for m and p.
Cheers,
Pedro

[aeroamit]

Hi Pedro,

Thanks for sharing the nice idea. If it helps in ramping up the CFL than really good....

One query - generally how Jacobian computation by finite difference or AD tend to be in comparison to analytical Jacobians in terms of cost (time). Probably Analytical jacobians are more efficient but they are difficult to compute for a given scheme (is this mostly true?).

Thanks
Amit

[pcarruscag]

Hi Amit,
That is also my intuition, I guess the only way to know if it is worth the extra cost is by doing.
Do you have a good supersonic case I can use for testing (all my work is subsonic)? Preferably something that converges well without initialization.
Cheers,
Pedro

[aeroamit]

Hi Pedro,

SU2 Testcase repo has 3 supersonic TestCases under euler folder (wedge, biparabolic and bluntbody), they all go well without initialization (you may be already aware of these euler cases).

But I do not have any specific supersonic case for RANS or Laminar flow problems (some of them ideally be like nozzle flow or SHOCK-BL interaction).
May be need to have a look at NASA CFD site for validation and verification cases.

Best
Amit

[pcarruscag]

Ok, I will start with those when I have a minute.
What about you @WallyMaier, do you have any reference supersonic RANS case?
Pedro

[WallyMaier]

@pcarruscag @aeroamit I do not have a good test case. But am currently looking for one. Let me know if y'all find a good.

[economon]

All for improving the abstractions if possible. We did something similar for the scalar upwinding and the viscous schemes already.

My only comments, which you can probably guess already, are to make sure we don't take a large performance hit by adding more layers (for example, one can make an argument to remove the entire set of CNumerics classes and implement the methods more efficiently directly in the solver class with vectorization, etc) and that we maintain the flexibility for folks to quickly add new convective schemes which was the original motivation for the existing structure (and more layers could complicate this). Sounds like you're already considering these things, but it is important to find the right balance.

[pcarruscag]

The Bogeyman came for the Roe classes. Roe, L2Roe, and LMRoe now share a parent that implements the fc_{1/2} = kappa*(fc_i+fc_j)*Normal part, the rest is then done by the derived classes. And with this I concluded this type of refactoring.

[pcarruscag]

As discussed in #705, according to the literature the CUSP flux should be expressed as the average of the left and right fluxes plus the artificial diffusion part, and not as a flux of averages + diffusion. This way the scheme becomes 100% upwind above Mach 1, and as a bonus the implementation also becomes more efficient.

Whatever I could do to reduce duplication and improve convergence is done, this PR is ready for reviews.

[EduardoMolina]

Hi Pedro,

Do you know why this PR is failing during Travis? I looked at the test cases and both are unsteady using Roe?
BTW, I will be happy to review this PR.
Eduardo

[pcarruscag]

Hi Eduardo,
The channel 2D and supersonic vortex shedding cases are failing due to small changes in residuals, the operations in Roe were re-ordered a bit, and these cases run for many iterations, I am guessing it is just due to numerical precision. But I will have a look before updating the residuals.
The rotating cylinders case is the infamous one that appears to be sensitive to minute, and seemingly insignificant, changes as reported in #700, at the reference residuals the case is actually diverged.
Finally there is an AD case using JST that also fails because it is a file diff and we output 15 significant digits to that file. The change is due to the refactoring of central schemes that required changing the order in which variables are registered for pre-accumulation.
I do appreciate you taking the time to review this PR, especially since you are "first author" of the low Mach Roe schemes and there are no regressions for them, nevertheless I tested them before and after and it seems I did not break anything.
Cheers,
Pedro

[pcarruscag]

@EduardoMolina
I ran the channel_2D and the
supersonic_vortex_shedding testcases with a tighter tolerance on the linear solver to make sure the differences were not due to some change in number of iterations.
Some differences do accumulate over time but they are very small, for the channel_2D the final density field differs by 1e-6 at most (which is the output precision).
For the supersonic vortex shedding differences are greater, 2e-4 on average and up to 0.01 maximum, but the case is also larger and runs for longer. If I restart the case for a 6th time-step (from the same restart files) the initial residuals are the same.
I will update the other two cases when #700 gets merged as they changed there too.
Cheers,
Pedro

[EduardoMolina]

Hi @pcarruscag
I am also running a couple of test cases, mainly the RANS NASA Hump. Do you think it is a good idea to add since we don't have any test case covering the low Mach Roe implementation? If so, I am happy to add.
Best,
Eduardo

[pcarruscag]

Hi @EduardoMolina ,
Thanks! Let me know how it goes.
I think having coverage is important but maybe it would be better to re-purpose an existing case since in our regressions Roe/JST are over-sampled and the regressions already take quite a bit of time to run.
Cheers,
Pedro

[EduardoMolina]

Hi @pcarruscag ,

Sorry for the late response.

I went through all your implementation and performed some tests.

Although, I didn't see an improvement using the accurate jacobians in the subsonic test cases, as already mentioned here, and even in the transonic OneraM6 case from the repo. This implementation is very clean and in my opinion is a great improvement.

I hope that in other cases, the use of accurate jacobian will have a positive effect in the convergence.

The only thing that I would like to bring is that in the future, you open a PR from su2 repo instead from your personal account. It is easier for the reviewers to compile and test.

I will wait Travis pass to merge this PR.

Thanks again.

Eduardo

[EduardoMolina]

LGTM!
Eduardo

[pcarruscag]

Thank you @EduardoMolina.
In the future I will make sure the branch exists in the SU2 repo.
Yes unfortunately the accurate Jacobians for SLAU and AUSM+ are not silver bullets and their success seems to be very case dependent.
The rotating_cylinders case stopped diverging after merging the changes from #700 so I have updated all reference residuals.

[aeroamit]

Hi @pcarruscag, All,

It's good to see one of this important PR coming to completion. It was a fruitful interaction, and I got to know more stuff.
I am not aware of detailed review process for SU2 PRs, and would like to know more details (I certainly go along with the fact that any member should be eligible for specific case/topic to review at first place, but useful details and ideas can be shared at any point - correct me if I missed something).
I could see the code changes for this PR in bits and pieces, as I didn't have my lappy with me for nearly last two weeks, and it's difficult/unpleasant on smartphone. Anyway, I will go through it in detail and extract the juice from it :) .

Thanks and Regards
Amit

[EduardoMolina]

Hi @pcarruscag and @aeroamit ,

Thanks for the discussion and for share the ideas/results. I think this is the beauty of open-source, everyone is welcome to jump in and review all the pull requests. Certainly, as you said @aeroamit, we will all learn something fruitful when we review and merge each PR.

Best regards,

Eduardo