Introduce a function to compute stress tensor

[maxaehle]

Proposed Changes

Provide a function CNumerics::ComputeStressTensor to compute the laminar and turbulent stress tensor. Replace the code blocks currently computing it, which exist at many locations inside the code, by function calls.

Related Work

?

To be discussed

Is CNumerics the best location for this function?

4 of the 5 rans_uq regression tests fail (in addition to those that failed before the change already). This seems to be due to floating-point errors in CSourcePieceWise_TurbSST::SetReynoldsStressMatrix, CAvgGrad_Base::SetReynoldsStressMatrix. These little deviations (relative difference around 1e-16) are probably amplified over the 10 iterations, so that the solver takes a notably different path to the steady state solution and the intermediate logarithmic residuals differ. I checked (for the 1c and 2c test cases) that the converged solutions are nearly equal though. (I can supply vtu files and plots if requested.)
Should I update the stored residual values in serial_regression.py then?

Also I am not sure whether I have found all code blocks that can be replaced by the new function.

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

• [X ] I am submitting my contribution to the develop branch.
• [X ] My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags, or simply --warnlevel=2 when using meson).
• [X ] My contribution is commented and consistent with SU2 style.
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp) , if necessary.

[pcarruscag]

Absolutely awesome, this is my favorite kind of pull request!

If you feel like it, I think the QCR modification is also repeated in some places, and there is a bug in how it is computed (see #992).

See if the template suggestion makes sense to simplify things even a bit more. It should be fine to update the failing testcases (they are very sensitive ones).

[pcarruscag]

That is a good question, you could try adding it to see how much it breaks the regressions. Maybe someone knows if it is for stability reasons.

[clarkpede]

Just noticed this. Where to include the turbulent kinetic energy is a messy issue. If I'm not mistaken, @economon removed the turb_ke from this calculation to fix the pressure problems described in: #797

[pcarruscag]

Regarding finding all code blocks, there is numerics_simd/flow/diffusion/common.hpp :: stressTensor.
But it is not very easy to conciliate pointers with the data types used in those vectorized functions.

[clarkpede]

Thanks for the contribution! This is a much needed change. I wanted to perform a refactoring like this myself. The viscous (or turbulent) stress calculation is repeated many times in the code. Changing that calculation (e.g. adding QCR or UQ) currently requires changing many files.

In my humble opinion, the stress tensor calculation should have its own class at the lowest level of the code, like the gas model (e.g. ideal gas) or the viscosity calculation (e.g. Sutherland's law). CNumerics classes need access to it, and CNumerics is meant to have minimal dependencies. But the stress isn't really "numerics." Its part of the underlying equations we're trying to solve. Changing the stress tensor isn't the same flavor as changing from JST to Roe. Because the stress is part of the underlying equations, it is sometimes needed in the CSolver classes. Previously, I've seen the stress tensor used in CNSSolver boundary conditions and axisymmetric source terms. That may still be true in the current version of the code.

The best way to make sure that both CSolver and CNumerics have access to it is to make it an independent class. Unless, of course, we want to be manually setting the six tensor elements in CNumerics, like we do with primitive variables. I don't like that alternative.

[pcarruscag]

Maybe CNumerics is not the perfect place, but it is good enough for government work (there are much more misplaced things in there).

The "everything is class" OO approach applied to the lowest level of abstraction is... well I think it is terrible -- and it has taken me a mighty long time to get rid of it in CPoint and CVariable and to design alternative Numerics -- because:

• Boilerplate: Set this, get that, constructor, destructor;
• Thread safety: Those classes always end up having some mutable state that renders them thread-unsafe;
• Correctness: Many of the classes we have follow this paradigm of "pass by member variable" - I like to know what are the inputs and outputs of something just by looking at the signature;
• Slowww: Too much virtual;
• Unnecessary complexity: A case of using a canon to kill the mosquito, good code should be as simple as possible, if a function does the job then that is the level of abstraction we should use.

My introduction to C++ was also the "everything is class", then one day I read "From Mathematics to Generic Programming" and well, I started liking C++ a whole lot more.
The standard template library is incredibly successful, and it "just" provides some containers and generic algorithms which are functions. That is what we need in SU2, some decent containers and generic algorithms to operate on them.

OO and its patterns are very good high level tools to achieve encapsulation and to isolate code, which are very important for projects with millions of lines of code, but for low level things they are overkill.

Sorry for the rant, I guess I have strong feelings about tiny classes.

[pcarruscag]

Thank you for the changes @maxaehle , LGTM

[clarkpede]

@pcarruscag I agree about OO being overkill for the stress tensor. My main concerns are:

• We find a way to share the stress tensor calc between solver and numerics classes.
• We make switching between different turbulent stress models (like normal vs QCR vs UQ) painless, without if/then branches being hardcoded for every time the stress is computed.

That can be accomplished in a number of ways, including functional or data-oriented programming. And, since I'm not the one submitting the PR, I understand that my opinions are just opinions. Pay as much or as little attention to them as you like.

[pcarruscag]

Ah! uff glad we agree (I do pay attention to what you think).
I think the first concern is more or less covered, it can always be cut and pasted somewhere else, we have sgs and wall models floating around, maybe we should have a physical_models directory.
Regarding your second point I would just have the stressTensor function deciding what modifications to perform, QCR and UQ, based on input parameters.

[clarkpede]

I just poked around in flow_sources.cpp and CNSSolver. All the places where I remember there being a hard-coded stress tensor calculation have been updated in this PR. I couldn't find any missing calculations. Great work.

In CNSSolver::AddDynamicGridResidualContribution, there is still a hardcoded Jacobian calculation. That will break if the stress tensor calculation changes independently, without also changing the Jacobian calculation. I don't think that affects this PR. I just wanted to document that here, in case later changes are made.

[maxaehle]

I searched for TWO3 in the project and this revealed two other locations. I adapted the one in CSolver.cpp in commit 64bd077. There is another one in python_wrapper_structure.cpp, lines 400-415. There I hesitate a bit because in the incompressible case, div_vel is manually set to zero. As long as the solver has not converged, this has an effect, right?

[maxaehle]

I'll wait for the regression tests to fail and then update serial_regression.py (and hybrid, parallel) for the rans_uq tests.
EDIT: Also stat_fsi in the hybrid and poiseuille_profile in the serial regression tests were affected.
EDIT: And poiseuille_profile in parallel

[maxaehle]

Now the tests are fine. Can I merge?

[pcarruscag]

You can replace the computation in python_wrapper_structure.cpp, the benefit of having a central definition outweighs the velocity divergence thing.
Also the way that method is used expects a reasonably converged solution, so it should not make a big difference.
Moreover that assumption for incompressible flow is not correct for variable density flows, so yes please replace it if you can.
I guess the TWO3 had a purpose in the end xD

[maxaehle]

I observed that sometimes only the product of the stress tensor with a vector is of interest, e.g.:

CNumerics::ComputeStressTensor(nDim, Tau, ...)
for (iDim = 0; iDim < nDim; iDim++) {
  for (jDim = 0 ; jDim < nDim; jDim++) {
    someforce[iDim] += Tau[iDim][jDim]*Normal[jDim];
  }
}
// never use Tau again

Do you think it is reasonable to have another function, which takes someforce and Normal instead of Tau as an argument and performs the addition? This would save from allocating Tau and going through two pairs of nested for loops.

[pcarruscag]

Maybe if that was in the hot path of the code it would be worth trying.
But so long as tau is of static size the compiler seems to do a good job at unrolling those loops https://gcc.godbolt.org/z/9rcrsM

[maxaehle]

Ok, so I will merge the pull request now?

[pcarruscag]

I think you can merge since this is a much needed change (but usually we give it some time so everyone has a chance to review).
But in this case I would say it's fine to merge. Thanks again!

[TobiKattmann]

But so long as tau is of static size the compiler seems to do a good job at unrolling those loops https://gcc.godbolt.org/z/9rcrsM

oh man, I totally saw the unrolled loops in assembler as well 🤥 ... maybe one day 😢

[TobiKattmann]

Hi @maxaehle , I just looked over this PR (after merging 🐌 ) and thanks a lot, it looks really tidy and I appreciate the detailed function documentation.
Unfortunately I cant help with the TODO why ignore turb_ke? question