Implementing variable density for unsteady incompressible flow#2641
Implementing variable density for unsteady incompressible flow#2641tkiymaz wants to merge 102 commits into
Conversation
|
@Cristopher-Morales can you review this pr please? |
Removed commented-out code for setting density at local point.
Updated the comment for density retrieval to reflect changes for transient density handling.
Removed unnecessary blank lines in CIncEulerVariable.hpp
Removed unnecessary blank lines in CVariable.hpp
Removed unnecessary blank lines to improve code readability.
Hi @pcarruscag ! Yes, I can help reviewing this PR. Please let me know if you need something else from my side |
Removed commented-out code and updated density calculation.
|
@Cristopher-Morales I guess the changes to the preconditioner can be removed from this PR since we now have your implementation |
Removed unnecessary blank lines in SetPrimVar function.
Removed unnecessary blank lines in CVariable.cpp.
Removed debug print statement from SetPrimVar function.
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com>
…/su2 into fix_inc_unsteady_density
| inline su2double GetDensity_time_n(unsigned long iPoint) const override { | ||
| return Density_time_n.size() > 0 ? Density_time_n(iPoint) : GetDensity(iPoint); | ||
| } |
|
@tkiymaz @pcarruscag |
| it was 0.0 before SetPrimitive_Variables evaluated the exact field. | ||
| Checking GetDensity_time_n(0) == 0.0 ensures we only evaluate this | ||
| once when the history arrays are uninitialized. ---*/ | ||
| if (dual_time && nPoint > 0 && nodes->GetDensity_time_n(0) == 0.0 && iMesh == MESH_0) { |
There was a problem hiding this comment.
This is not OpenMP-safe: CommonPreprocessing runs inside the solver's parallel region, and this worksharing loop sits behind a branch whose condition is invalidated by the loop itself. A fast thread can enter the omp for and write point 0 while a slower thread has not evaluated the condition yet; the slow thread then sees a nonzero value and skips the construct. A worksharing construct encountered by only some threads of a team is undefined behavior (in practice a hang, and at minimum a TSan report in the CI). The decision to initialize must be made in a single-threaded context (BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS, or a flag set once outside the parallel region), not by probing the data from all threads.
| SU2_OMP_ATOMIC | ||
| ErrorCounter += SetPrimitive_Variables(solver_container, config); | ||
|
|
||
| /*--- Provide a valid initial Density_time_n for restart or step 0 since |
There was a problem hiding this comment.
Three related problems with this initialization strategy:
- Multigrid: the
iMesh == MESH_0guard means the coarse-grid variable objects (which also allocateDensity_time_n = 0) are never initialized, andSetResidual_DualTimeruns on all MG levels — soMGLEVEL > 0+ variable density + dual time uses rho_n = 0 on the coarse levels for the first time step. Either initialize all levels or error out for this combination. - Restart consistency: on restart, rho_n and rho_n-1 are set to the current density rather than densities consistent with
Solution_time_n/Solution_time_n1, so the first step after restart is not a true continuation (and the 2nd-order history is wrong — the PR says only 1st order is implemented, but the 2nd-order branches exist and nothing guardsDUAL_TIME_STEPPING-2ND_ORDER). - The
== 0.0sentinel probed every inner iteration is fragile even once the race is fixed.
All three disappear with one change: recompute Density_time_n{,1} from the stored solution histories via the fluid model once per time step (after the push-back), instead of lazily cloning the current density. That is correct after restart by construction, needs no sentinel, works on every MG level it runs on, and — if the recomputation sits in the recorded path — the dual-time residual's dependence on rho_n flows back to the already-registered Solution_time_n inputs on the tape, closing the adjoint chain rule that is currently left as a TODO, with no manual extraction needed.
| /*! | ||
| * \brief Virtual function returning whether this is the species solver. | ||
| */ | ||
| inline virtual bool IsSpeciesSolver() const { return false; } |
There was a problem hiding this comment.
These two virtuals don't belong on the CSolver base: their only call site is CScalarSolver::SetResidual_DualTime, called on this. A protected bounded_scalar member of CScalarSolver, set by the derived constructors from GetBounded_Species() / GetBounded_Turb() (same pattern as the existing Conservative flag), does the job with zero base-class footprint and no per-call config queries.
| * \param[in] iPoint - Point index. | ||
| * \return Density at time level n. | ||
| */ | ||
| inline virtual su2double GetDensity_time_n(unsigned long iPoint) const { |
There was a problem hiding this comment.
Two things here:
- These override the
CVariableversions, so per the style guide they should be markedoverride(and then don't needvirtual). - More fundamentally, the base declarations on
CVariableare not needed at all. The only call site that goes throughCVariable*isCScalarSolver::SetResidual_DualTime, and the idiom for exactly this situation is already used in the same file (CScalarSolver.inl:156):su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes()). With that cast, these getters can be declared here onCFlowVariable(compressible default) with theCIncEulerVariableoverride, and the dangerous silent-0.0default onCVariabledisappears — which is exactly what makes the uninitialized coarse-MG case fail silently instead of loudly.
| * \param[in] iPoint - Point index. | ||
| * \return Density at time n; 0 for compressible nodes. | ||
| */ | ||
| inline virtual su2double GetDensity_time_n(unsigned long /*iPoint*/) const { return 0.0; } |
There was a problem hiding this comment.
See the comment on CFlowVariable — these declarations can be removed from CVariable entirely by casting to CFlowVariable* at the scalar-solver call site. A base default that silently returns 0.0 for the density turns any missed override or missed initialization into a silently zeroed unsteady term. The doc comments here ("incompressible only; 0 for compressible nodes") also contradict the CFlowVariable behavior.
| * \param[in] iPoint - Point index. | ||
| * \return Adjoint of the density at time n. | ||
| */ | ||
| inline su2double GetAdjointDensity_time_n(unsigned long iPoint) const override { |
There was a problem hiding this comment.
GetAdjointDensity_time_n{,1} are never called anywhere — the only references are the TODO comments in CDiscAdjSolver.cpp. Please remove them (here and the CVariable defaults), along with the RegisterSolution_time_n{,1} overrides that register inputs whose adjoints are never extracted, and add them together with the actual chain-rule extraction in a follow-up (or better, avoid the manual chain rule entirely by recomputing the density histories on tape — see the comment in CIncEulerSolver.cpp). Until the adjoint is closed, an error for unsteady discrete adjoint + non-constant INC_DENSITY_MODEL would prevent silently wrong gradients.
|
|
||
| #include "../../../Common/include/geometry/CGeometry.hpp" | ||
| #include "../../include/solvers/CSolver.hpp" | ||
| #include "../../include/variables/CIncEulerVariable.hpp" |
There was a problem hiding this comment.
This include appears to be unused — nothing else in this file changed and no CIncEulerVariable symbol is referenced.
5 participants
Proposed Changes
Implements variable density treatment for unsteady incompressible flow simulations. Currently, SU2 uses constant density in transient simulations, which is inaccurate for combustion cases where density varies significantly due to heat release and species composition changes. This contribution enables proper density updates during time-stepping for flamelet-based combustion modeling. For now, only 1st order time marching is implemented
Related Work
Related to incompressible flow solver and flamelet combustion modeling. No specific issue linked yet.
PR Checklist
pre-commit run --allto format old commits.