Bugfix BC sym plane (compr. and incompr. solver)

[TobiKattmann]

Proposed Changes

Hi all,

In neither solver (incompressible and compressible) viscous fluxes are considered in the symmetry plane boundary condition. In both codes the BC_sym_plane links directly to Euler_Wall.

Compressible: A "reflected state" approach is done. Create a reflected halo node which enforces symmetry and compute the occurring flux just as in the interior. Problem: the primitive gradients are not adapted nor a viscous flux based on that is computed.

Incompressible: Here a "direct flux" approach is chosen in Euler Wall. Insert the symmetry conditions directly in the FV derivation on the boundary part and thereby eliminate most terms. For the convective terms only a pressure contribution in the momentum equation remains. No Viscous terms are currently considered which is not correct for the momentum equation.

For both solvers I implemented a "reflected state" approach in C*Solver::BC_sym_plane as I didn't want to touch the Euler_wall implementations especially from the compressible solver.

I (informally) wrote a bit down as explanation in Latex how I did the derivation.
Symmetry_BC.pdf
Note that in the DG-FEM solver of @vdweide the same "reflected state" approach is implemented although it is not clear to me how the derivation was done there. See the bottom of the pdf for more info. I would be happy to get some insight here :)

I also did a bit of testing with some cases and will showcase the differences between current develop and my implementation with pictures below. I attached two easy 2D cases (compr. & incompr. half cylinder) which you could quickly compare the current develop and this PR with.
Sym_plane_testcases.gz
tar xzf Sym_plane_testcases.gz

1. Compressible half cylinder, navier-stokes, no Turb-model, Pin has a solid wall, farfield around. You see pressure contour lines:
   

2. Incompressible half cylinder in channel, navier-stokes, no turb-model, velocity inlet pressure outlet. You see zebra shaded pressure contour lines:
   

3. 3D incompressible pin setup, SST turb model. You see Pressure-contour lines but more important are the velocity vectors. In the current code they have a normal component and therefore a (non-irrelevant) flux over the symmetry boundary:
   

Additional to that I tested in rotated domains such that the symmetry plane does not coincide with one of the axis. I tested the Navier-Stokes ONERA-M6 with SA turb-model and although this PR achieves a considerable lower velocity in symmetry-normal direction than current develop, it is still quite large. I suspect it is due to general "bad" convergence (2-3 orders of magnitude) of that case.
So, I encourage everyone that has a well converging compressibe 3D case with symmetry plane to test that or to send it to me.

Every regression test with symmetry plane in it will fail, so I didn't change the values yet. From what I have seen, only these are failing. I will correct these (fml) and commit if everyone is alright with the implementation otherwise.

Let me know if you have any comments/suggestions etc., or if you tested with your cases and can confirm (or not) that this PR is correct (that would be awesome).

Cheers, Tobi

Related Work

None known to me

PR Checklist

I will check the warnings based on travis output and report. No new testcase is necessary in my opinion.

Edit 25.02.: I checked the warnings for a serial and parallel build using the log-files of travis, comparing mine with the latest merged develop. No new compiler warnings are generated.

• 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.

[economon]

This is an important improvement.. thanks for the effort and extra analysis @TobiKattmann

[TobiKattmann]

Concerning the Regression tests:
All failing test have a symmetry plane and are based on the inc-/compressible solver (so far so good).
They are all based on the comparing the SU2 output, not by a file-diff like some AD cases do with the of_grad.dat file. That makes updating the regression test quite easy :) .
The difference between computed and stored values mostly quite small, apart from cases where a restart is done. I could change the solution_flow.dat where necessary if wanted.

[economon]

Nice job! I can confirm that everything is working ok for me with some incompressible cases. For example, even extruding the 2D laminar cylinder mesh from the test cases by one cell in the z-direction and running a quasi-2D calculation with symmetry planes in the z-direction gives almost indistinguishable results from the true 2D calculation.

Only thing not addressed yet is the symmetry BC for the turbulence models (we prob need to add something similar). @TobiKattmann : do you have time to look at that as part of this PR? Or shall we do that separately?

Just a few comments in the review.

[TobiKattmann]

Only thing not addressed yet is the symmetry BC for the turbulence models (we prob need to add something similar). TobiKattmann : do you have time to look at that as part of this PR? Or shall we do that separately?

The implementation for the turbulence models is correct as is. No flux is considered for the convective and viscous fluxes by inserting the symmetry conditions directly in the FV derivation. Convective fluxes vanish as velocity*normal=0, viscous fluxes vanish by grad(transported scalar)*n=0 on the symmetry faces of the boundary. Above I denoted that as "direct flux" approach.
(I will attach an updated pdf here later)

Although the "reflected state" approach is a more proactive approach and might be beneficial for convergence. Could be worth a shot if somebody really has an issue here, for now I leave it as is.

NEXT STEP I will update the Regression test values at the failing cases. For cases with a restart (either primal or adjoint computations) I will update the restart_file and upload these in a seperate Testcases branch.

[TobiKattmann]

All regression test values are corrected and for restarted/adjoint cases (euler/oneram6, disc_adj_heat, disc_adj_euler/arina2k) the solution_flow.dat was replaced. The respective Testcases branch is here https://github.com/su2code/Testcases/tree/bugfix_BC-sym-plane . From my side this PR could be merged (after merging the Testcases branch and reverting .travis.yml).

[economon]

LGTM! Thanks for the clean up.

Last chance to others for questions/comments. Otherwise, we can merge this in within a day or two.

[pcarruscag]

Hi @TobiKattmann,
I cannot help but notice that the implementation is identical for the compressible and incompressible case.
Is it not possible to avoid duplicating the code?
Cheers,
Pedro

[talbring]

@pcarruscag Indeed both implementations are very similar. But this is the case for a lot of routines in CEulerSolver and CIncEulerSolver at the moment. We can think about introducing a common base class for both in the future to remove this code duplication. But for now there is no other way except for moving the implementation to CSolver which we should avoid in my opinion.

[talbring]

Thanks for this. Very nice coding style with a lot of useful comments. I am happy and we can merge this is in.

Just one quick note below :)

[pcarruscag]

@talbring, if you are concerned with implementing BC_Sym_Plane in CSolver and every derived class inheriting a default implementation, it can be named BC_Sym_Plane_Flow and then the compressible and incompressible solvers can call it from their implementation of BC_Sym_Plane. Anything is better than copy pasting code, as that will only make a future refactoring more difficult.

By the way, I too think the implementation is very good and I really appreciate @TobiKattmann 's effort to document it from the theoretical and implementation point of view.

[talbring]

@pcarruscag I just had a chat with @TobiKattmann. Essentially there are two points that would, in our opinion, speak against moving the implementation to CSolver.

1. We dont know yet whether they might be some future differences in the implementations.
2. This defeats somehow the purpose of the class structure, as the base class should be free of specific implementations for a certain solver ... Although the intention would be to have it there only temporarily, we never know how long it actually stays there in the end. I don't mind having a little bit of code copy, if the structure is clear and easy to understand.

[pcarruscag]

Yes, well maybe there should be something between CSolver and the C(Inc)Euler solvers. If you run a file diff between the two cpp's you'll see what all the small copies add up to. Like you said no one knows how long these things stay in the code. If in the future you need different implementations you specialise them, but until then, you have only one piece of code to maintain / improve / compile. Anyway...

bigger problems in the world

[TobiKattmann]

Hi @pcarruscag ,

Yes, well maybe there should be something between CSolver and the C(Inc)Euler solvers. If you run a file diff between the two cpp's you'll see what all the small copies add up to.

I agree, introducing an (e.g.) CFlowSolver between CSolver and C(Inc)EulerSolver could be a solution to the 1:1 code duplications of multiple functions including this BC_sym_plane. It is worth to keep that idea in mind, although it is probalby not practical to implement it in this PR.

Thanks for the review!
Tobi

[economon]

I think we are all in agreement then.. as @pcarruscag points out, we should always keep code duplication in mind and try to abstract when possible. But, let's do it gradually.. sometimes we need to balance this against keeping duplicate code if it helps reduce conflicts and allows for different parts of the code to evolve separately and in parallel by different developers.

I think a good first step is what @talbring has in mind for the flow output.. exactly this type of mid-level flow output class will show up there I believe. We can see how well that goes and then port to the solver classes too.

Will merge this in today so we can keep moving forward. Thanks all