su2code · tkiymaz · Jun 17, 2025 · Jun 17, 2025 · Sep 1, 2025 · Dec 10, 2025
diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml
@@ -233,7 +233,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}}
+          args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}}
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
@@ -282,7 +282,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests"
+          args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} -a "--tapetests"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054
         with:
@@ -330,7 +330,7 @@ jobs:
           PMIX_MCA_gds: hash
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan"
+          args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} -a "--tsan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054
         with:
@@ -375,7 +375,7 @@ jobs:
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054
         with:
           # -t <Tutorials-branch> -c <Testcases-branch>
-          args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan"
+          args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} -a "--asan"
       - name: Cleanup
         uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054
         with:

diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl
@@ -28,6 +28,7 @@
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../../include/solvers/CScalarSolver.hpp"
 #include "../../include/variables/CFlowVariable.hpp"
+#include "../../include/variables/CIncEulerVariable.hpp"
 
 template <class VariableType>
 CScalarSolver<VariableType>::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative)
@@ -624,9 +625,14 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
   const bool second_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND);
   const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE);
 
+  /*--- The bounded-scalar correction is currently only used by the incompressible solvers. ---*/
+  const bool bounded_scalar = incompressible && ((IsSpeciesSolver() && config->GetBounded_Species()) ||
+                                                 (IsTurbSolver()    && config->GetBounded_Turb()));
+
   /*--- Flow solution, needed to get density. ---*/
 
   CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
+  auto* incFlowNodes = incompressible ? su2staticcast_p<CIncEulerVariable*>(flowNodes) : nullptr;
 
   /*--- Store the physical time step ---*/
 
@@ -655,13 +661,9 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
     for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
       if (Conservative) {
         if (incompressible) {
-          /*--- This is temporary and only valid for constant-density problems:
-          density could also be temperature dependent, but as it is not a part
-          of the solution vector it's neither stored for previous time steps
-          nor updated with the solution at the end of each iteration. */
-          Density_nM1 = flowNodes->GetDensity(iPoint);
-          Density_n = flowNodes->GetDensity(iPoint);
-          Density_nP1 = flowNodes->GetDensity(iPoint);
+          Density_nM1 = incFlowNodes->GetDensity_time_n1(iPoint);
+          Density_n = incFlowNodes->GetDensity_time_n(iPoint);
+          Density_nP1 = incFlowNodes->GetDensity(iPoint);
         } else {
           Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];
           Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
@@ -686,13 +688,20 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
        time discretization scheme (1st- or 2nd-order).---*/
 
       for (iVar = 0; iVar < nVar; iVar++) {
+        su2double unsteady_term = 0.0;
         if (first_order)
-          LinSysRes(iPoint, iVar) +=
-              (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * Volume_nP1 / TimeStep;
+          unsteady_term = (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * Volume_nP1 / TimeStep;
         if (second_order)
-          LinSysRes(iPoint, iVar) += (3.0 * Density_nP1 * U_time_nP1[iVar] - 4.0 * Density_n * U_time_n[iVar] +
+          unsteady_term = (3.0 * Density_nP1 * U_time_nP1[iVar] - 4.0 * Density_n * U_time_n[iVar] +
                                       1.0 * Density_nM1 * U_time_nM1[iVar]) *
                                      Volume_nP1 / (2.0 * TimeStep);
+
+        if (bounded_scalar) {
+          if (first_order) unsteady_term -= U_time_nP1[iVar] * (Density_nP1 - Density_n) * Volume_nP1 / TimeStep;
+          if (second_order) unsteady_term -= U_time_nP1[iVar] * (3.0 * Density_nP1 - 4.0 * Density_n + 1.0 * Density_nM1) * Volume_nP1 / (2.0 * TimeStep);
+        }
+
+        LinSysRes(iPoint, iVar) += unsteady_term;
       }
 
       /*--- Compute the Jacobian contribution due to the dual time source term. ---*/
@@ -721,7 +730,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
 
       if (Conservative) {
         if (incompressible)
-          Density_n = flowNodes->GetDensity(iPoint);  // Temporary fix
+          Density_n = incFlowNodes->GetDensity_time_n(iPoint);
         else
           Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
       }
@@ -777,7 +786,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
 
           if (Conservative) {
             if (incompressible)
-              Density_n = flowNodes->GetDensity(iPoint);  // Temporary fix
+              Density_n = incFlowNodes->GetDensity_time_n(iPoint);
             else
               Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
           }
@@ -818,13 +827,9 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
         /*--- If this is the SST model, we need to multiply by the density
          in order to get the conservative variables ---*/
         if (incompressible) {
-          /*--- This is temporary and only valid for constant-density problems:
-          density could also be temperature dependent, but as it is not a part
-          of the solution vector it's neither stored for previous time steps
-          nor updated with the solution at the end of each iteration. */
-          Density_nM1 = flowNodes->GetDensity(iPoint);
-          Density_n = flowNodes->GetDensity(iPoint);
-          Density_nP1 = flowNodes->GetDensity(iPoint);
+          Density_nM1 = incFlowNodes->GetDensity_time_n1(iPoint);
+          Density_n = incFlowNodes->GetDensity_time_n(iPoint);
+          Density_nP1 = incFlowNodes->GetDensity(iPoint);
         } else {
           Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];
           Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
@@ -833,19 +838,35 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
       }
 
       for (iVar = 0; iVar < nVar; iVar++) {
+        su2double unsteady_term = 0.0;
         if (first_order)
-          LinSysRes(iPoint, iVar) +=
-              (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nP1 / TimeStep);
+          unsteady_term = (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nP1 / TimeStep);
         if (second_order)
-          LinSysRes(iPoint, iVar) +=
-              (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) +
+          unsteady_term = (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) +
               (Density_nM1 * U_time_nM1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nM1 / (2.0 * TimeStep));
+
+        if (bounded_scalar) {
+          if (first_order) unsteady_term -= U_time_nP1[iVar] * (Density_nP1 - Density_n) * (Volume_nP1 / TimeStep);
+          if (second_order) unsteady_term -= U_time_nP1[iVar] * ((Density_nP1 - Density_n) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) +
+                                                                 (Density_nM1 - Density_n) * (Volume_nM1 / (2.0 * TimeStep)));
+        }
+
+        LinSysRes(iPoint, iVar) += unsteady_term;
       }
 
       /*--- Compute the Jacobian contribution due to the dual time source term. ---*/
       if (implicit) {
-        if (first_order) Jacobian.AddVal2Diag(iPoint, Volume_nP1 / TimeStep);
-        if (second_order) Jacobian.AddVal2Diag(iPoint, (Volume_nP1 * 3.0) / (2.0 * TimeStep));
+        su2double diag_factor = 1.0;
+        if (Conservative) {
+          if (bounded_scalar) {
+            if (first_order)  diag_factor = Density_n;
+            if (second_order) diag_factor = (4.0 * Density_n - Density_nM1) / 3.0;
+          } else {
+            diag_factor = Density_nP1;
+          }
+        }
+        if (first_order)  Jacobian.AddVal2Diag(iPoint, diag_factor * Volume_nP1 / TimeStep);
+        if (second_order) Jacobian.AddVal2Diag(iPoint, diag_factor * 3.0 * Volume_nP1 / (2.0 * TimeStep));
       }
     }
     END_SU2_OMP_FOR
@@ -884,4 +905,4 @@ void CScalarSolver<VariableType>::PushSolutionBackInTime(unsigned long TimeIter,
       nodes->Set_Solution_time_n();
     }
   }
-}
+}
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
@@ -347,6 +347,16 @@ class CSolver {
    */
   inline virtual CFluidModel* GetFluidModel(void) const { return nullptr;}
 
+  /*!
+   * \brief Virtual function returning whether this is the species solver.
+   */
+  inline virtual bool IsSpeciesSolver() const { return false; }
+
+  /*!
+   * \brief Virtual function returning whether this is the turbulent solver.
+   */
+  inline virtual bool IsTurbSolver() const { return false; }
+
   /*!
    * \brief Get number of linear solver iterations.
    * \return Number of linear solver iterations.

diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp
@@ -69,6 +69,11 @@ class CSpeciesSolver : public CScalarSolver<CSpeciesVariable> {
 
   void Initialize(CGeometry* geometry, CConfig* config, unsigned short iMesh, unsigned short nVar);
 
+  /*!
+   * \brief Virtual function returning whether this is the species solver.
+   */
+  inline bool IsSpeciesSolver() const final { return true; }
+
   /*!
    * \brief Restart residual and compute gradients.
    * \param[in] geometry - Geometrical definition of the problem.

diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp
@@ -55,6 +55,11 @@ class CTurbSolver : public CScalarSolver<CTurbVariable> {
    */
   CTurbSolver(CGeometry* geometry, CConfig *config, bool conservative);
 
+  /*!
+   * \brief Virtual function returning whether this is the turbulent solver.
+   */
+  inline bool IsTurbSolver() const final { return true; }
+
   /*!
    * \brief Impose via the residual the Euler wall boundary condition.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -147,7 +152,7 @@ class CTurbSolver : public CScalarSolver<CTurbVariable> {
    * \returns The number of extra variables.
    */
   unsigned long RegisterSolutionExtra(bool input, const CConfig* config) final;
-  
+
   /*!
    * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over
    * a nonlinear iteration for stability.

diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp
@@ -71,6 +71,8 @@ class CIncEulerVariable : public CFlowVariable {
 
   VectorType Streamwise_Periodic_RecoveredPressure,    /*!< \brief Recovered/Physical pressure [Pa] for streamwise periodic flow. */
              Streamwise_Periodic_RecoveredTemperature; /*!< \brief Recovered/Physical temperature [K] for streamwise periodic flow. */
+  VectorType Density_time_n,                           /*!< \brief Density at time n for dual-time stepping. */
+             Density_time_n1;                          /*!< \brief Density at time n-1 for dual-time stepping. */
   su2double TemperatureLimits[2];                      /*!< \brief Temperature limits [K]. */
  public:
   /*!
@@ -87,6 +89,16 @@ class CIncEulerVariable : public CFlowVariable {
   CIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy,
                     unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config);
 
+  /*!
+   * \brief Set all the solution at time level n to the current solution value (including density).
+   */
+  void Set_Solution_time_n() override;
+
+  /*!
+   * \brief Set all the solution at time level n-1 to the solution at time level n (including density).
+   */
+  void Set_Solution_time_n1() override;
+
   /*!
    * \brief Set the value of the pressure.
    * \param[in] iPoint - Point index.
@@ -291,4 +303,68 @@ class CIncEulerVariable : public CFlowVariable {
     for (unsigned long iDim = 0; iDim < nDim; iDim++) Solution(iPoint, iDim+1) = val_vector[iDim];
   }
 
+  /*!
+   * \brief Get the density at time level n for dual-time stepping.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n.
+   */
+  inline su2double GetDensity_time_n(unsigned long iPoint) const { return Density_time_n(iPoint); }
+
+  /*!
+   * \brief Get the density at time level n-1 for dual-time stepping.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n-1.
+   */
+  inline su2double GetDensity_time_n1(unsigned long iPoint) const { return Density_time_n1(iPoint); }
+
+  /*!
+   * \brief Set the density at time level n for dual-time stepping.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_density - Density value.
+   */
+  inline void SetDensity_time_n(unsigned long iPoint, su2double val_density) { Density_time_n(iPoint) = val_density; }
+
+  /*!
+   * \brief Set the density at time level n-1 for dual-time stepping.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_density - Density value.
+   */
+  inline void SetDensity_time_n1(unsigned long iPoint, su2double val_density) { Density_time_n1(iPoint) = val_density; }
+
+  /*!
+   * \brief Register the density at time n as an AD input variable (for discrete adjoint unsteady).
+   */
+  void RegisterDensity_time_n() override;
+
+  /*!
+   * \brief Register the density at time n-1 as an AD input variable (for discrete adjoint unsteady).
+   */
+  void RegisterDensity_time_n1() override;
+
+  /*!
+   * \brief Get the adjoint of density at time n.
+   *        Simplified version: directly extracts the AD derivative of Density_time_n,
+   *        analogous to GetAdjointSolution_time_n for the solution vector.
+   *        The chain-rule contribution d(rho)/d(h) is currently not applied (TODO).
+   * \param[in] iPoint - Point index.
+   * \return Adjoint of the density at time n.
+   */
+  inline su2double GetAdjointDensity_time_n(unsigned long iPoint) const override {
+    AD::Identifier index = AD::GetPassiveIndex();
+    AD::SetIndex(index, Density_time_n(iPoint));
+    return AD::GetDerivative(index);
+  }
+
+  /*!
+   * \brief Get the adjoint of density at time n-1.
+   *        See GetAdjointDensity_time_n.
+   * \param[in] iPoint - Point index.
+   * \return Adjoint of the density at time n-1.
+   */
+  inline su2double GetAdjointDensity_time_n1(unsigned long iPoint) const override {
+    AD::Identifier index = AD::GetPassiveIndex();
+    AD::SetIndex(index, Density_time_n1(iPoint));
+    return AD::GetDerivative(index);
+  }
+
 };
diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp
@@ -282,12 +282,12 @@ class CVariable {
   /*!
    * \brief Set the variable solution at time n.
    */
-  void Set_Solution_time_n();
+  virtual void Set_Solution_time_n();
 
   /*!
    * \brief Set the variable solution at time n-1.
    */
-  void Set_Solution_time_n1();
+  virtual void Set_Solution_time_n1();
 
   /*!
    * \brief Set the variable solution at time n.
@@ -2227,6 +2227,48 @@ class CVariable {
    */
   void RegisterUserDefinedSource();
 
+  /*!
+   * \brief Get the stored density at time n (incompressible only).
+   * \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; }
+
+  /*!
+   * \brief Get the stored density at time n-1 (incompressible only).
+   * \param[in] iPoint - Point index.
+   * \return Density at time n-1; 0 for compressible nodes.
+   */
+  inline virtual su2double GetDensity_time_n1(unsigned long /*iPoint*/) const { return 0.0; }
+
+  /*!
+   * \brief Register the density at time n as an AD input (no-op for non-incompressible nodes).
+   */
+  inline virtual void RegisterDensity_time_n() {}
+
+  /*!
+   * \brief Register the density at time n-1 as an AD input (no-op for non-incompressible nodes).
+   */
+  inline virtual void RegisterDensity_time_n1() {}
+
+  /*!
+   * \brief Get the adjoint values of the density at time n.
+   *        Implemented analogously to GetAdjointSolution_time_n but for a scalar
+   *        (Density_time_n is a VectorType, not a MatrixType).
+   *        Currently ignores the chain-rule d(rho)/d(h) dependence (TODO).
+   * \param[in] iPoint - Point index.
+   * \return Adjoint value of the density at time n.
+   */
+  inline virtual su2double GetAdjointDensity_time_n(unsigned long /*iPoint*/) const { return 0.0; }
+
+  /*!
+   * \brief Get the adjoint values of the density at time n-1.
+   *        See GetAdjointDensity_time_n.
+   * \param[in] iPoint - Point index.
+   * \return Adjoint value of the density at time n-1.
+   */
+  inline virtual su2double GetAdjointDensity_time_n1(unsigned long /*iPoint*/) const { return 0.0; }
+
   /*!
    * \brief Set the adjoint values of the solution.
    * \param[in] adj_sol - The adjoint values of the solution.