Different results on the same testcase of unsteady ONERA (v6.0.1 vs v6.2.0)

[RoccoBombardieri]

Description of the issue
Hello everybody,

I'm currently running the same testcase (Unsteady Euler ONERAM6 AoA=3deg) on two different versions of SU2 ( version 6.0.1 and 6.2.0). What I noticed is a different evolution in time of coefficients, Cl and Cd as given in the figure. Both the codes look at convergence and the input files (despite the version of the code being different) are the same. I attach the two screen outputs for comparison. I'm running SU2_CFD through the parallel_computation.py with only one core (using different number of cores doesn't affect the result anyway).

Summary_singlecore_6_2_0.txt
Summary_singlecore_6_0_1.txt

To Reproduce the problem
To reproduce the problem I attach the config files. The mesh is available in the repository: mesh_ONERAM6_inv.su2

Note: both versions are/were on the develop branches. Latest one is up to date.

config_CFD_6_0_1.txt
config_CFD_6_2_0.txt

Im running the two versions on different machines with same specs:
Details for the new version

• OS: [Ubuntu 16.04.6 LTS]
• C++ compiler and version: [g++ 5.4.0]
• MPI implementation and version: [Open MPI 1.10.2]
• SU2 Version: [v6.2.0]

Details for the old version

• OS: [Ubuntu 16.04.6 LTS]
• C++ compiler and version: [g++ 5.4.0]
• MPI implementation and version: [Open MPI 1.10.2]
• SU2 Version: [v6.0.1]

Thank you and best,
Rocco

[TobiKattmann]

Hi Rocco,
I used this mesh <Testcases>/control_surface/mesh_ONERAM6_inv.su2 and the boundary marker names are a bit different in the mesh, compared to your provided config (WING vs LOWER_SIDE, UPPER_SIDE, TIP + SYMMETRY vs SYMMETRY_FACE). Can you shed some light on this?

Will the airfoil emit some vortex-shedding (or other unsteadiness) or will this converge in time to a steady-state solution? I am asking because I would like to see some longer running results until e.g. steady state, cylic shedding or maybe we'll see some chaotic behaviour in lift and drag. The initial transient through the domain based on the initial values is unphysical to some extend and if e.g. 'steady state' results are the same for both code versions .... the problem is not that big. Although it differs quite a bit in the temporal evolution tbh

In the results (also in the mail) it looks a bit like it approaches a steady state for Cl and CD. So maybe 1000 timesteps will tell us a bit more 🤔

Another thing: There is no tag 6.0.2 😕 ... 6.0.1 and then 6.1.0..
And both versions are/were on the develop branches -> what exactly do you refer to here: config-files, code, mesh

Cheers 🐌
Tobi

[RoccoBombardieri]

Hello Tobi,

Thanks for your quick reply! I'll address you doubts as follows.

I used this mesh /control_surface/mesh_ONERAM6_inv.su2 and the boundary marker names are a bit different in the mesh, compared to your provided config (WING vs LOWER_SIDE, UPPER_SIDE, TIP + SYMMETRY vs SYMMETRY_FACE).

Yeah, the mesh you are mentioning is the same. I simply united the TIP, LOWER_SIDE and UPPER_SIDE under the marker tag WING.

The initial transient through the domain based on the initial values is unphysical to some extend and if e.g. 'steady state' results are the same for both code versions .... the problem is not that big. Although it differs quite a bit in the temporal evolution tbh
In the results (also in the mail) it looks a bit like it approaches a steady state for Cl and CD. So maybe 1000 timesteps will tell us a bit more 🤔

I guess they converge eventually to the same values. I tried to compare the same case but steady (I attach the summary). Results in this case are fairly similar and (in my opinion) raise no warnings.
Summary_steady_6_0_1.txt
Summary_steady_6_2_0.txt

Regarding the unsteady case, I'm aware that this test may have little physics behind (especially if comparing the first timesteps) but I wanted to investigate the reason :). Results are too different and I think these might be one of the reasons behind some discrepancies I'm encountering in my FSI framework.
My framework for FSI features a restart from a steady solution at t=0 and an unsteady simulation with imposed boundary of the wing marker. The discrepancies in that case are huge and are clearly wrong. I went back to the root and found this weird behaviour comparing the two versions so I thought this might be one cause.

Another thing: There is no tag 6.0.2 😕 ... 6.0.1 and then 6.1.0..

Yep you are right, the version is definitely 6.0.1 as it can be seen from the Summaries: I'm updating the issue.

And both versions are/were on the develop branches -> what exactly do you refer to here: config-files, code, mesh

Didn't quite understand you here.

Let me know if you need any other info regarding the topic. Looking forward to hear from you!

Best,
Rocco

[TobiKattmann]

Hey Rocco,

And both versions are/were on the develop branches -> what exactly do you refer to here: config-files, code, mesh

The code tags 6.0.1 and 6.2.0 refer to specific master-commits of the code. So I want to know what is exactly of the develop branch, can't be the code if you understand the tags as I do. Maybe you refer to the develop code at the time where these tags were introduced. But maybe sth different.

I guess they converge eventually to the same values. I tried to compare the same case but steady (I attach the summary). Results in this case are fairly similar and (in my opinion) raise no warnings.

That is good to hear.
Luckily @cvencro (Thanks for that 👍 ) ran some tests and found that #740 (which rewrites Euler Wall boundaries for compressible and incompressible flows ... which was me btw) introduces the changes you see. The Euler wall for compressible is changed a bit in the way the boundary condition is enforced. Upon convergence both version should recover (close to) the same value. If the results in a steady state are the same/very-similar and if an unsteady simulation at low enough Reynolds-number converges to the same steady result (or oscillate in the same bounds) -> I would say it is simply some purely numerical effect based on the unphysical initialization and different treatment of the euler wall boundary.

1 thing one could test additionally is to force a very high convergence in each physical timestep, i.e. usually one converges ~3 orders of magnitude in each ts -> now try converging until residual stalling in each ts. If the code versions get the same/closer results in the temporal evolution my euler wall boundary explanation would be quite reasonable.
Together with the long running unsteady simulation proposed above, I would feel more comfortable to judge.

@cvencro (bringing the conversation to this place 😉 ) concerning the pitching case: the if (dynamic_grid) statements in the 'old' version are not necessary as that is handled by the numerics class. In the new version an appropriate 'reflected state' is constructed and the numerics container is called to compute the residual. Before, the code of one numerics ->ComputeResidual Routine was simply copied and slightly modified.

But of course there is always room for errors 🐛

Cheers, Tobi

[cvencro]

Hi Tobi,

I ran some tests using the latest code in develop, as of Friday 6th December, but swapping between the symmetry plane boundary condition and the previous Euler wall boundary condition.

Results are below, where firstly the steady results were very similar as Rocco mentioned. In time-domain, they look to converge to similar values but the difference in transient behaviour can be seen. The solutions were converged by about 8 orders of magnitude at each time step for both cases. I reran the current develop code with symmetry boundary condition forcing all the way to residual stalling as suggested for a few hundred time steps but I think the case was pretty well converged already and there is no difference with increased convergence.

Although it is good that the converged behaviour is similar for both boundary conditions, the difference in transient behaviour will have a large effect for coupled unsteady problems.

About the dynamic grid, the reason I spotted this was because the gradient verification against rigid motion with symmetry plane boundary condition was good but not as good as the gradient verification with Euler wall boundary condition (as well as having slightly different inherent sensitivity magnitudes between the two boundary conditions caused by the difference in transient behaviour). Therefore, this is perhaps missing something as I would have expected to see the same level of agreement in the sensitivities for the symmetry plane boundary condition as is seen for the Euler wall. I’ve added the results using both boundary conditions for a pitching airfoil case using efficiency as the objective function and FFD control points as the design variables. The differences are seen at the trailing edge (points 10, 11, 21, 22).

Should the Euler wall boundary condition be replaced by the symmetry plane condition or could we have both as options to choose from?

Best wishes,
Charanya

[RoccoBombardieri]

Charanya,

thanks for the detailed answer. Let me ask you some more info. Apparently, you where able to reproduce the results I had comparing the two code versions simply using different boundary conditions on the latest code version. Can you specify me how?

Tobi,

In the meanwhile, to narrow down the problem, I ran also a couple tests in 2D to seek confirmation. For the (Euler) pitching_NACA64A010.cfg test case in the repository I rerun the same test at AoA of 1 deg and removing the pitching (no mesh deforming). I attach config files and summaries relative to the test cases.

config_CFD_6_2_0.txt
Summary_6_2_0.txt
config_CFD_6_0_1.txt
Summary_6_0_1.txt

In this case the situation is definitely better as the two solvers give the same results. It looks like it is an issue related to the 3D case.

I also quote Charanya in saying that it is definitely good that convergence to same values is reached but the transient response is fundamental and holds physical/mathematical values (I think of Wagner et similia) :). I would also suggest to add a regression test in this sense.

Best,
Rocco

[cvencro]

Hi Rocco,

My approach to test with both boundary conditions was rather crude. I just pasted back the Euler wall BC that was removed in #740 . You can see what I did if you look at the change to solver_direct_mean.cpp as part of the draft PR #833 . While debugging, it is easy to rebuild the code with true/false if statement as needed to swap between the two boundary conditions.

Thanks for the 2D case with no grid motion comparison you’ve run. The two boundary conditions do seem to give the same result in this case and helps narrow down the search. However, I don’t think the differences are restricted only to 3D cases because with either non-zero rigid motion grid movement or deforming mesh, there are still differences in the 2D solution.

Best wishes,
Charanya

[TobiKattmann]

Hey,

@cvencro and I were discussing this issue quite a bit this morning and here an attempt of a summary:

We are looking at the following cases:

1. 3D Onera m6, compressible euler, including euler_wall and sym_plane
   a. steady state
   b. unsteady (no pitching, deforming)
2. 2D NACA64A010, compressible euler, including euler_wall
   a. unsteady (no pitching/ deformation)
   b. pitching (with rigid and with deforming mesh -> used for the gradient validation of @cvencro 's post )

We are rather certain that the differences between the code-versions are due to the new euler_wall boundary which was introduced in #740 (by me :) ). @cvencro did a test where the old euler_wall was simply pasted into the newer function body (which currently directly calls the sym_plane boundary) -> that recovered the results obtained with the 'older' version.

The steady state results (onera m6) show that both codes deliver the same results (of course with some minor differences). So far so good.

Unsteady cases are a bit trickier: the latest 2D results of @RoccoBombardieri show pretty much no difference including the first timesteps. The onera m6 case in contrast converges to the same (steady-state) results for both code version eventually but their initial transient phase differs quite significantly.

A possible explanation for the different initial transient phase: In both implementations a "reflected state" is constructed where the normal component of the velocity is subtracted [once in the older code version & twice in the newer version]. I.e. subtracting the normal component of velocity twice means mirroring the velocity along the symmetry plane. The newer code version is a bit "stronger" in enforcing the boundary conditions which could lead to the reduced oscillation in CD and CL.

If we agree that both code versions produce physically correct results for steady cases and unsteady cases that were integrated long enough in time to get rid of the influence of the initial condition ... then we are left with this initial transient phase that differs. Of course if there is more evidence that one or the other version produces physically "better" results I would love to see that. Until then, both initial oscillations are equally correct/uncorrect in my opinion.

[I'll add a few words about the adjoint here later]

[I'll add the idea of an FSI case of @cvencro here later where the initial phase has a major impact]
Consider an FSI computation of an airfoil where the trailing edge can exhibit flutter (periodic up-and-down-movement of the trailing edge) in certain flow regimes. Now if during the initial transient the forces on the airfoil are higher than in the converged state then the fluttering can be excited where the initial transient can be seen as an activation energy. If you were to e.g. ramp up flow speed/conditions slowly up to the same magnitude as before one might end up without fluttering. So based on the flow history there might exist multiple "stable" states.

Cheers, Tobi

[TobiKattmann]

@cvencro @RoccoBombardieri

I guess I found the root of the difference for the different CL & CD (case: 3D compressible Euler ONERAM6):

1. As already said, the difference was introduced in Allow BC_Sym_Plane to be used as Slip wall for viscous flow #740 where I kicked out the "old" Euler implementation in favour of the symmetry plane
   --> old Euler BC: This is basically a JST scheme without dissipation term on the BC. A "reflected state" is created where scalars are copied and velocity is mirrored along the plane. Then the quantities (domain & reflected) are averaged (on the edge) and CNumerics::GetInviscidProjFlux is called with these values. The returned Proj_Flux is directly taken as Residual. For JST dissipation would be added in a step afterwards.
   -> new Euler BC: Here for compressible flow always (!) CUpwRoe_Flow is taken as numerical convective scheme. Boundaries have their own numerics container in CONV_BOUND_TERM. Look in CDriver.cpp for allocation of SPACE_CENTERED compressible schemes (as JST was used in all examples in this issue, but the problem is still there if upwind schemes were taken). And independent of the users choice of CONV_NUM_METHOD_FLOW, as said above, CUpwRoe_Flow is hard coded as numerical convective scheme for boundary nodes. As in the old Euler BC a "reflected state" is created and handed to the ComputeResidual Routine of the numerical scheme, together with the boundary vertex state.
   -> So one major difference between both formulations is that in the new Euler wall left & right state are created and the respective convective scheme class does the rest via ComputeResidual. In the old implementation everything done in the ComputeResidual routine is "copied" into the BC_Euler_Wall and modified to its needs... and it looks like central convective scheme (JST) without dissipation term

2. When changing the CUpwRoe_Flow scheme, which is used only for convective schemes on boundaries CONV_BOUND_TERM, such that it behaves like JST without dissipation term -> then old Euler and new Euler wall BC return the same results in CD and CL.

3. So the Convective Residual is somewhat different on the surface where the Objectives CL and CD are evaluated. JST without dissipation term vs CUpwRoeBase_Flow::ComputeResidual. Which (and if) one of both approaches is the "corect" is not something I am able to say with great confidence. I tend towards: "Well, just as you get (slightly) different results if you choose JST or ROE in the domain due to the different residual formulation, you get different results if different convective schemes on the boundary vertices are chosen. In the domain the choice is ~ Upwind if discontinuities like shocks are present but JST is a computationally cheaper. How this translates to boundary nodes ... I dont know."

I would like to get some opinion of someone who has more knowledge of Convective schemes. Maybe @economon :) ? And of course anyone else.

Cheers Tobi

[cvencro]

Hi Tobi,

Great to see the same dynamic behaviour from both boundary conditions for the ONERA M6 case after the changes you've tested for the convective flux residuals. The differences seen in the dynamic results are important and since the change in boundary condition wasn't intended to cause these changes in the dynamic solution, I think we should be able to reproduce what we had originally. I think the same is also true for the problems with grid motion, as a minimum we should be able to reproduce the same result with the new boundary condition as we had with the old Euler wall.

If it then comes down to choices of different convective schemes that we can employ with the new boundary condition then perhaps as further work the different available choices could be investigated to explore which option gives the most accurate result compared against references in literature?

I ran some tests for 2D cases with grid motion using both boundary conditions for us to use here. I started with the NACA64A010 airfoil @RoccoBombardieri tested, now with the pitching motion re-included, and looked at the periodic behaviour.

pitching_NACA64A010.zip

Using the Euler wall boundary condition, the results still match exactly to @economon 's thesis (no real surprise as I suspect the code hasn't changed much for this comparison) and this agrees pretty closely with other inviscid reference results like Jamesons's unsteady flow past airfoils. With the symmetry plane, the results differ from the other inviscid solutions but the agreement to the experiment is better. I couldn't get hold of reference drag or moment results to compare against, therefore the comparison is only against lift.

I also tested the pitching NACA 0012 airfoil using Batina's "Unsteady Euler Airfoil Solutions Using Unstructured Dynamic Meshes" as reference. In this case with the NACA 0012 airfoil, the reduced frequency (k=0.0814) is lower than the previous case and the amplitude is slightly higher at 2.51deg.

pitching_NACA0012.zip

Once again, using reference from computational and experimental results, the lift results are very similar for both boundary conditions. For this case, moment results were also available.

The sign convention was inverted in the paper and neither BC matched exactly to the reference which had a coarser mesh than the one used here but the Euler wall looks to be in better agreement with the reference and captures the overall shape. I've attached the case files for this as well and it would be interesting to see the results when run with the change you've made to reproduce the same convective flux strategy at the boundary.

As a final test I ran the rotating frame case, again from @economon 's thesis which is a steady problem with simulated grid motion. I was able to run and replicate the flow characteristics using the Euler wall boundary condition but with symmetry plane the residuals were straightaway at machine precision and the flow never develops around the airfoil.

rotating_frame.zip

Overall from these 2D airfoil cases, there are some unexpected variations between the two boundary conditions. A closer look at grid motion with the symmetry plane is needed. It would also be interesting to run the cases with symmetry plane using the convective flux residuals calculated in the same way as the Euler wall to see if that takes them closer to the original results, similarly as seen in the fixed grid problem.

Best wishes,
Charanya

[RoccoBombardieri]

Hello Charanya, Tobi,

I've finally performed some validation runs trying to tackle the Boundary condition issue which was recognized by @TobiKattmann as the source of discrepancy in unsteady Euler simulations between the old version and the new version.

My validation test case relies on the ONERA M6 model I originally used to point out the issue. The test consists in imposing a deformation of the wing boundary which resembles a bending mode of the wing (starting from an already converged steady state). Such input deformation is provided in time in the shape of a blended-step. For the considered mode, forces are recorded in the fashion of the popular (within an aeroelastic context) Generalized Aerodynamic Force (GAF) coefficient in time. A Fourier Transform allows me to recast the GAF in the frequency domain where I can compare the same value with the ones extracted by the Doublet Lattice Method on a geometrically equivalent panel mesh. Simulation is performed at mach 0.8(~ V_inf = 272.3 m/s ) and AoA = 0

In the figure here you can see the time domain GAF vs time as calculated by old SU2 (v 6.0.1), new SU2 (v6.2.0) enforcing the old boundary condition as suggested previously by @cvencro and new SU2 with the new boundary condition.

You can notice how at approx 0.0023 sec the blended step is fully developed. Now, despite some small differences between the old SU2 version and the new with enforced old boundary conditions (that are due, among other things, to small changes in the mesh deformation algorithm that I have noticed) you can notice the big difference in the application of the two boundary conditions.

I also attach a couple zooms of the signal: first at time 0 which gives the GAF for the steady initial simulation and second for the converged values. It can be seen how (even if less important) some discrepancies are found here as well.

To check which one of the two time evolution is physically representative you can observe the GAF in the frequency domain. The following figure plots module of the Fourier transform of the preceding signals together with the DLM reference as a function of the reduced frequency (proportional to omega/V_inf).

Old SU2 has been validated vs DLM in another test which I don't show here for sake of brevity and differences with respect to DLM are justified by the thickness of the ONERA M6 and high speed. It can be seen how new SU2 with old boundary conditions gives back an almost identical result to old SU2 (my reference). Anyway for the new boundary conditions values are too distant from the reference.
Little note: the Module of the GAF in frequency domain at K=0 is proportional to the difference of the value of the gaf in time domain between t=0 and t=infinite: namely GAF(t=infinite) - GAG(t=0). The same value sums other coefficients due to the unsteadiness of the signal at K different than 0. We can see how for both the boundary conditions the value approximates the original one (as we already pointed out at the beginning, noticing that both steady and converged values were very very similar for the two boundary conditions) Anyway, by increasing the frequency, we notice that differences increase rapidly pointing out the problem must be somehow related to the rate of unsteadiness of the flow.

In conclusion, now I can say with good confidence that old boundary conditions are the correct ones. I can't, anyway highlight precisely the source of the problem.

[TobiKattmann]

Hi @RoccoBombardieri ,

Can you retry your case with the branch testing_symmetry_plane and add the following settings:

% Surf; YES=oldEuler; NO=NewEuler
VISUALIZE_SURFACE_DEF= NO 
% Vol; YES=UpwRoe; NO=JST_BC
VISUALIZE_VOLUME_DEF= NO
% YES=JST; NO=see VOLUME_DEF
MUSCL_ADJFLOW= NO

In that branch I added the possibility to run the old or new Euler_Wall BC on MARKER_EULER walls (MARKER_SYM always uses the the new implementation) using VISUALIZE_SURFACE_DEF as a hooked boolean.
With VISUALIZE_VOLUME_DEF one can change which convective scheme is used on the boundary, i.e. what is in the numerics container CONV_BOUND_TERM. The Default is CUpwRoe_Flow (ROE) and I added CCentJST_Flow_BC which is the normal JST without dissipation.
The third option MUSCL_ADJFLOW has to be left on NO, this is work in progress :)

With the options given above you will use the new Euler wall implementation but using JST without dissipation term as a convective scheme for boundary nodes. This (as I say in point 2. in my above post) basically recovers the old euler wall implementation.

For the OneraM6 case the combination of new Euler + JST_BC recovers the CD & CL results of the old code.

Can you also try running CONV_NUM_METHOD_FLOW= ROE with the code you labelled as 'SU2 FSI new (new bound condition)' and 'SU2 FSI new (oldbound condition)'. Because by doing so you switch the numerical method in the volume and to the same ROE scheme on the boundary -> just as the current default in the develop. So my current assumption is that it is a matter of convective scheme (on the boundary).

I'll try to attach my cases as a zip and two images here.
Here the green line represents the current develop if one uses MARKER_EULER or MARKER_SYM. Another info here is: by using the oldEuler wall (i.e. line 1&2) the choice of convective scheme on the boundary is irrelevant for the Euler_wall as the old implementation "skips" that step. The red line 4 is the setup I want you to try because this recovers the old results in this case.

Here the testcase folder if someone wants to use the exact setup
oneram6test.tar.gz