turbulence->validate(); in pisoFoam.C changes nut, thereby making the flow different from the case without save-and-restart.
Run a simulation for 10 steps, save it, and continue for another 10 steps. Compare the outcome of the second 10 steps with the outcome of running 20 steps from the beginning. The difference will be in the 5th digit. But for chaotic turbulent flows, this difference will amplify over time, making simulations with and without restart totally different.
Chaotic turbulent flows are indeed chaotic and only the averaged statistics are reproducible. If you average the results for a long time, restart and restart the averaging do the averaged statistics agree between the two runs?
That's correct. I found this bug, which does not occur for other LES models, when experimenting with numerical methods for sensitivity analysis. Being able to restart and continue on exactly the same trajectory is important for me, and likely for other numerical methods, e.g., for computing Lyapunov exponents of the flow dynamics.
All turbulence->validate() does is call correctNut which updates the nut field in the same manner as it does after the calculation of nuTilda. This call is necessary to ensure the nut field is physical and consistent with nuTilda at the beginning of the run, particularly when starting for user-input fields in which nut may be specified as 0. There is no difference in this procedure between SpalartAllmaras and any of the other models.
You are right. turbulence->validate() should not have changed the nut that was computed in the last time step. But somehow it does when using SpalartAllmarasDES.
Incidentally SpalartAllmarasDDES and SpalartAllmarasIDDES are derived from SpalartAllmarasDES and inherit the same correctNut code.
I suspect that the original diagnosis isn't fully correct. The re-calculation of the "nut" field is very likely only a symptom of the real problem.
Furthermore, the original report refers to "SpalartAllmaras" and not "SpalartAllmarasDES", which increases my suspicion.
@qiqi: Can you please provide a test case that replicates the problem that you're seeing?
I ask this because the most likely suspect is actually the boundary conditions at the inlet(s), given that if it's using a random inlet source for flow speed and/or pressure values, that would justify the problem.
If this is the case, either the inlet values need to be pre-planned for each time step or a re-programmable random number generator would have to be implemented. The random generators in C++11 do have this feature, namely to get a seed for at a particular time step, to be re-use in the next run.
Note that all the boundary conditions are "clean" -- fixedValue, zeroGradient.
The test passes only with a modified version of pisoFoam, the only difference between which and the shipped version is the validate(); line.
The "nut" field was considerably different, although still in similar ranges of values.
OK, I believe I've got it. I see what's wrong. In "SpalartAllmarasDES<BasicTurbulenceModel>::correct()", "fv1" is calculated with the "nuTilda" from before "nuTildaEqn" was solved. "correctNut(fv1)" at the end of the method uses the outdated fv1 value, resulting in the discrepancy. Removing "fv1" from the call should re-establish consistency.
Attached is the file "SpalartAllmarasDES.C" that replaces "src/TurbulenceModels/turbulenceModels/LES/SpalartAllmarasDES/SpalartAllmarasDES.C".
OK... I'm not use to yet to the possibility to attach files along with the comment... please ignore the file "SpalartAllmarasDES-2.C", which is a duplicate of the previous upload "SpalartAllmarasDES.C".
@qiqi: I forgot to mention that after replacing the file, please run the "Allwmake" script in the folder "src/TurbulenceModels", because it needs to update at least a couple of libraries that depend on this template class.
@Henry: I've confirmed that the provided test case works as intended when using the attached file "SpalartAllmarasDES.C", for replacing "src/TurbulenceModels/turbulenceModels/LES/SpalartAllmarasDES/SpalartAllmarasDES.C"
so that the final calculated "nut" field within the method "correct()" is based on the latest "nuTilda" and not based on the one from before solving the equation.
I've also done a quick look through the other sibling template classes within LES and the only ones that use "correctNut()" with additional arguments are "dynamicKEqn" and "dynamicLagrangian", but those two rely on calculations done based on "U", which doesn't change within the "correct()" method.
Thanks for fixing this issue. Great job Henry.
I should be thanking @wyldckat. Thanks wyldckat. Great job.