Implementing variable density for unsteady incompressible flow

[tkiymaz]

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

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

[pcarruscag]

@Cristopher-Morales can you review this pr please?

[Cristopher-Morales]

@Cristopher-Morales can you review this pr please?

Hi @pcarruscag !

Yes, I can help reviewing this PR.

Please let me know if you need something else from my side

[bigfooted]

@Cristopher-Morales I guess the changes to the preconditioner can be removed from this PR since we now have your implementation

[bigfooted]

@pcarruscag like this?

[bigfooted]

@tkiymaz @pcarruscag
Done from my side, please have a look...

[pcarruscag]

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.

[tkiymaz]

Folded into the fix in abe9e19. The initialization is now guarded by InnerIter == 0, which is a uniform predicate, so the worksharing loop inside is entered by every thread of the team.

[pcarruscag]

Three related problems with this initialization strategy:

1. Multigrid: the iMesh == MESH_0 guard means the coarse-grid variable objects (which also allocate Density_time_n = 0) are never initialized, and SetResidual_DualTime runs on all MG levels — so MGLEVEL > 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.
2. 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 guards DUAL_TIME_STEPPING-2ND_ORDER).
3. The == 0.0 sentinel 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.

[tkiymaz]

Done in abe9e19. Density_time_n[,1] are now recomputed from Solution_time_n[,1] (plus the species histories on the flamelet path) with the fluid model. It happens once per physical time step in CommonPreprocessing.
I also updated the flamelet regression test_vals because the pre-fix values were tuned with the previous runs.

[pcarruscag]

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.

[tkiymaz]

Done in de8eb53. Both virtuals pulled off CSolver. CScalarSolver now stores a const bool BoundedScalar set from the derived constructors, same pattern as Conservative — CSpeciesSolver passes GetBounded_Species(), CTurbSolver passes GetBounded_Turb(), CHeatSolver passes false.

[pcarruscag]

Two things here:

1. These override the CVariable versions, so per the style guide they should be marked override (and then don't need virtual).
2. More fundamentally, the base declarations on CVariable are not needed at all. The only call site that goes through CVariable* is CScalarSolver::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 on CFlowVariable (compressible default) with the CIncEulerVariable override, and the dangerous silent-0.0 default on CVariable disappears — which is exactly what makes the uninitialized coarse-MG case fail silently instead of loudly.

[tkiymaz]

Done in 8120bda. GetDensity_time_n[,1] now live on CFlowVariable with the CIncEulerVariable override, and CScalarSolver::SetResidual_DualTime casts to CFlowVariable* using the same su2staticcast_p<> idiom. Thanks for pointing at the existing pattern.

[pcarruscag]

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.

[tkiymaz]

Fixed in same commit (8120bda). Silent-0 defaults are gone — a missed override would now fail to compile instead of quietly zeroing the unsteady term.

[pcarruscag]

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.

[tkiymaz]

Done in 324d94d — GetAdjointDensity_time_n[,1], the RegisterSolution_time_n[,1] density overrides, and the CVariable defaults are all gone, along with the two TODOs in CDiscAdjSolver::ExtractAdjoint_Solution. Added an SU2_MPI::Error in the discrete-adjoint constructor for unsteady + non-CONSTANT INC_DENSITY_MODEL.

[pcarruscag]

This include appears to be unused — nothing else in this file changed and no CIncEulerVariable symbol is referenced.

[tkiymaz]

Done in 12fd1f6.

[pcarruscag]

Suggested change

	* state at time n (also after restart), on every MG level where the inputs
	* exist, without the fragile per-iteration sentinel that used to seed it.
	* state at time n (also after restart).

[pcarruscag]

This type of AI comment is utterly useless (the strong words are for Claude to pay attention).

Suggested change

/*--- The bounded-scalar correction (set from the derived solver constructor). ---*/

[pcarruscag]

This is not the place to document another class hierarchy.

Suggested change

	/*--- Flow solution, needed to get density. The static cast to CFlowVariable*
	exposes GetDensity_time_n[,1] (declared on CFlowVariable, overridden by
	CIncEulerVariable) without a silent 0.0-returning base on CVariable.
	Same idiom as Upwind_Residual above. ---*/
	/*--- Flow solution, needed to get density. ---*/

[pcarruscag]

Suggested change

	inline su2double GetDensity_time_n(unsigned long iPoint) const override {
	inline su2double GetDensity_time_n(unsigned long iPoint) const final {

[pcarruscag]

Suggested change

	inline su2double GetDensity_time_n1(unsigned long iPoint) const override {
	inline su2double GetDensity_time_n1(unsigned long iPoint) const final {

[pcarruscag]

changes in this file are no longer needed, revert

[pcarruscag]

Probably not needed anymore since you are recomputing densities like I suggested in the first place.
Either way settings validation is done in CConfig not in the solvers so that dry run mode can catch issues like these.

[pcarruscag]

What is the purpose of repeating the function documentation?
What is the purpose of making a reference to an iteration of the code that will not be kept after merging into develop?

Suggested change

	/*--- Recompute the dual-time density history from the stored solution history
	via the fluid model, once per physical time step after the push-back.
	This is consistent with the primitive state at time n by construction
	(also after restart, since Solution_time_n[,1] are the loaded histories)
	and needs no sentinel. The InnerIter == 0 guard is a uniform, non
	data-dependent condition, so the worksharing loop inside is entered by
	every thread of the team -- OpenMP-safe, unlike the old probe of
	GetDensity_time_n(0) == 0.0. ---*/

[pcarruscag]

no it's not

[pcarruscag]

Suggested change

	/*--- For constant-density flows, save memory by using GetDensity() on demand. ---*/
	const bool variable_density = (config->GetKind_DensityModel() != INC_DENSITYMODEL::CONSTANT);
	if (variable_density) {
	/*--- Allocate density storage for time levels.
	Note: Actual density values will be set after SetPrimVar is called in Preprocessing. ---*/
	Density_time_n.resize(nPoint) = su2double(0.0);
	Density_time_n1.resize(nPoint) = su2double(0.0);
	}
	if (config->GetKind_DensityModel() != INC_DENSITYMODEL::CONSTANT) {
	Density_time_n.resize(nPoint) = su2double(0.0);
	Density_time_n1.resize(nPoint) = su2double(0.0);
	}