Fix/inc unsteady variable density

[Soumyadipta-Banerjee]

Proposed Changes

This fixes #2480 by removing the constant density assumption in the unsteady incompressible solvers.
Previously, the solvers were just calling nodes->GetDensity(iPoint) for all three dual-time levels (n+1, n, n-1). This breaks variable-density simulations (e.g., when running with the energy equation or species transport). There was a comment in the code noting this as a "temporary fix," so this PR finally addresses it.
I added proper historical density storage (Density_time_n and Density_time_n1) to CIncEulerVariable. CIncEulerSolver and CScalarSolver now grab the correct density for each time step when computing dual-time residuals. I also made sure the density pushback routines in CIntegration and CFVMFlowSolverBase use SU2_OMP_FOR_STAT and parallelCopy to stay thread-safe with OpenMP.

Related Work

Fixes #2480.
This takes the core idea from the closed draft PR #2481 but implements it fully for both 1st and 2nd-order time stepping.

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.

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I used the pre-commit hook to prevent dirty commits and used pre-commit run --all to format old commits.
• 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.

[bigfooted]

we already have a pull request for this, why open another one?

[Soumyadipta-Banerjee]

My apologies, I thought the other PR was closed so I opened this