Implementing variable density for unsteady incompressible flow

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

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -622,11 +622,14 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
   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. ---*/
+  const bool bounded_scalar = (IsSpeciesSolver() && config->GetBounded_Species()) ||
+                              (IsTurbSolver()    && config->GetBounded_Turb());
 
   /*--- Flow solution, needed to get density. ---*/
 
-  CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
+  auto* flowNodes = solver_container[FLOW_SOL]->GetNodes();
 
   /*--- Store the physical time step ---*/
 
@@ -654,19 +657,9 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
     SU2_OMP_FOR_STAT(omp_chunk_size)
     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);
-        } else {
-          Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];
-          Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
-          Density_nP1 = flowNodes->GetSolution(iPoint, 0);
-        }
+        Density_nM1 = flowNodes->GetDensity_time_n1(iPoint);
+        Density_n = flowNodes->GetDensity_time_n(iPoint);
+        Density_nP1 = flowNodes->GetDensity(iPoint);
       }
 
       /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that
@@ -686,13 +679,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. ---*/
@@ -720,10 +720,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
       U_time_n = nodes->GetSolution_time_n(iPoint);
 
       if (Conservative) {
-        if (incompressible)
-          Density_n = flowNodes->GetDensity(iPoint);  // Temporary fix
-        else
-          Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
+        Density_n = flowNodes->GetDensity_time_n(iPoint);
       }
 
       for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) {
@@ -776,10 +773,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
           /*--- Multiply by density at node i for the SST model ---*/
 
           if (Conservative) {
-            if (incompressible)
-              Density_n = flowNodes->GetDensity(iPoint);  // Temporary fix
-            else
-              Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
+            Density_n = flowNodes->GetDensity_time_n(iPoint);
           }
 
           for (iVar = 0; iVar < nVar; iVar++) LinSysRes(iPoint, iVar) += Density_n * U_time_n[iVar] * Residual_GCL;
@@ -815,37 +809,42 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
        due to the time discretization has a new form.---*/
 
       if (Conservative) {
-        /*--- 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);
-        } else {
-          Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];
-          Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
-          Density_nP1 = flowNodes->GetSolution(iPoint, 0);
-        }
+        /*--- Get density at different time levels via virtual methods ---*/
+        Density_nM1 = flowNodes->GetDensity_time_n1(iPoint);
+        Density_n = flowNodes->GetDensity_time_n(iPoint);
+        Density_nP1 = flowNodes->GetDensity(iPoint);
       }
 
       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 +883,4 @@ void CScalarSolver<VariableType>::PushSolutionBackInTime(unsigned long TimeIter,
       nodes->Set_Solution_time_n();
     }
   }
-}
+}

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.

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.

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.

SU2_CFD/include/variables/CFlowVariable.hpp

@@ -270,4 +270,26 @@ class CFlowVariable : public CVariable {
    * \return Vector of magnitudes.
    */
   inline su2activevector& GetStrainMag() { return StrainMag; }
+
+  /*!
+   * \brief Get the density at time level n for dual-time stepping (compressible flows).
+   * For compressible flows, density is extracted from the solution container.
+   * For incompressible flows, this is overridden.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n.
+   */
+  inline virtual su2double GetDensity_time_n(unsigned long iPoint) const {
+    return GetSolution_time_n(iPoint, 0);  // Density is the first component of Solution
+  }
+
+  /*!
+   * \brief Get the density at time level n-1 for dual-time stepping (compressible flows).
+   * For compressible flows, density is extracted from the solution container.
+   * For incompressible flows, this is overridden.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n-1.
+   */
+  inline virtual su2double GetDensity_time_n1(unsigned long iPoint) const {
+    return GetSolution_time_n1(iPoint, 0);  // Density is the first component of Solution
+  }
 };

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,74 @@ 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 override {
+    return Density_time_n.size() > 0 ? Density_time_n(iPoint) : GetDensity(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 override {
+    return Density_time_n1.size() > 0 ? Density_time_n1(iPoint) : GetDensity(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 solution at time n, including the density container for incompressible flows.
+   */
+  void RegisterSolution_time_n() override;
+
+  /*!
+   * \brief Register solution at time n-1, including the density container for incompressible flows.
+   */
+  void RegisterSolution_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 {
+    if (Density_time_n.size() == 0) return 0.0;  // Constant density - no adjoint contribution
+    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 {
+    if (Density_time_n1.size() == 0) return 0.0;  // Constant density - no adjoint contribution
+    AD::Identifier index = AD::GetPassiveIndex();
+    AD::SetIndex(index, Density_time_n1(iPoint));
+    return AD::GetDerivative(index);
+  }
+
 };