Fix/inc unsteady variable density

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1073,6 +1073,8 @@ void CFVMFlowSolverBase<V, R>::PushSolutionBackInTime(unsigned long TimeIter, bo
   for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
     solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n();
     solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n1();
+    solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Density_time_n();
+    solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Density_time_n1();
     if (rans) {
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n1();
@@ -1103,6 +1105,7 @@ void CFVMFlowSolverBase<V, R>::PushSolutionBackInTime(unsigned long TimeIter, bo
 
     for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
       solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n();
+      solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Density_time_n();
       if (rans) solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
 
       geometry[iMesh]->nodes->SetVolume_n();

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -635,12 +635,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);
+          /*--- Use density at each time level for variable-density flows. ---*/
+          Density_nM1 = flowNodes->GetDensity_time_n1(iPoint);
+          Density_n = flowNodes->GetDensity_time_n(iPoint);
           Density_nP1 = flowNodes->GetDensity(iPoint);
         } else {
           Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];
@@ -701,7 +698,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
 
       if (Conservative) {
         if (incompressible)
-          Density_n = flowNodes->GetDensity(iPoint);  // Temporary fix
+          Density_n = flowNodes->GetDensity_time_n(iPoint);
         else
           Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
       }
@@ -757,7 +754,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
 
           if (Conservative) {
             if (incompressible)
-              Density_n = flowNodes->GetDensity(iPoint);  // Temporary fix
+              Density_n = flowNodes->GetDensity_time_n(iPoint);
             else
               Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
           }
@@ -798,12 +795,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);
+          /*--- Use density at each time level for variable-density flows. ---*/
+          Density_nM1 = flowNodes->GetDensity_time_n1(iPoint);
+          Density_n = flowNodes->GetDensity_time_n(iPoint);
           Density_nP1 = flowNodes->GetDensity(iPoint);
         } else {
           Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -72,6 +72,10 @@ 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. */
   su2double TemperatureLimits[2];                      /*!< \brief Temperature limits [K]. */
+
+  VectorType Density_time_n;   /*!< \brief Density at time level n for dual-time stepping (variable-density flows). */
+  VectorType Density_time_n1;  /*!< \brief Density at time level n-1 for dual-time stepping (variable-density flows). */
+
  public:
   /*!
    * \brief Constructor of the class.
@@ -282,6 +286,34 @@ class CIncEulerVariable : public CFlowVariable {
     return Streamwise_Periodic_RecoveredTemperature(iPoint);
   }
 
+  /*!
+   * \brief Get the density at time level n.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n.
+   */
+  inline su2double GetDensity_time_n(unsigned long iPoint) const final {
+    return Density_time_n(iPoint);
+  }
+
+  /*!
+   * \brief Get the density at time level n-1.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n-1.
+   */
+  inline su2double GetDensity_time_n1(unsigned long iPoint) const final {
+    return Density_time_n1(iPoint);
+  }
+
+  /*!
+   * \brief Store current density into time level n: Density_time_n = Density.
+   */
+  void Set_Density_time_n() final;
+
+  /*!
+   * \brief Store time level n density into time level n-1: Density_time_n1 = Density_time_n.
+   */
+  void Set_Density_time_n1() final;
+
   /*!
    * \brief Specify a vector to set the velocity components of the solution.
    * \param[in] iPoint - Point index.

SU2_CFD/include/variables/CVariable.hpp

@@ -970,6 +970,30 @@ class CVariable {
    */
   inline virtual su2double GetDensity(unsigned long iPoint, unsigned long val_iSpecies) const { return 0.0; }
 
+  /*!
+   * \brief Get the density at time level n (for incompressible unsteady flows with variable density).
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n. Defaults to current density.
+   */
+  inline virtual su2double GetDensity_time_n(unsigned long iPoint) const { return GetDensity(iPoint); }
+
+  /*!
+   * \brief Get the density at time level n-1 (for incompressible unsteady flows with variable density).
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n-1. Defaults to current density.
+   */
+  inline virtual su2double GetDensity_time_n1(unsigned long iPoint) const { return GetDensity(iPoint); }
+
+  /*!
+   * \brief Store current density into the time level n container: Density_time_n = Density.
+   */
+  inline virtual void Set_Density_time_n() {}
+
+  /*!
+   * \brief Store time level n density into the time level n-1 container: Density_time_n1 = Density_time_n.
+   */
+  inline virtual void Set_Density_time_n1() {}
+
   /*!
    * \brief A virtual member.
    * \param[in] iPoint - Point index.

SU2_CFD/src/integration/CIntegration.cpp

@@ -249,6 +249,10 @@ void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver
   solver->GetNodes()->Set_Solution_time_n1();
   solver->GetNodes()->Set_Solution_time_n();
 
+  /*--- Store density at previous time levels (for incompressible variable-density flows). ---*/
+  solver->GetNodes()->Set_Density_time_n1();
+  solver->GetNodes()->Set_Density_time_n();
+
   SU2_OMP_SAFE_GLOBAL_ACCESS(solver->ResetCFLAdapt();)
 
   SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads()))

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -2804,7 +2804,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
   su2double U_time_nM1[MAXNVAR], U_time_n[MAXNVAR], U_time_nP1[MAXNVAR];
   su2double Volume_nM1, Volume_nP1, TimeStep;
   const su2double *Normal = nullptr, *GridVel_i = nullptr, *GridVel_j = nullptr;
-  su2double Density;
+  su2double Density_nP1, Density_n, Density_nM1;
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
@@ -2841,15 +2841,17 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       V_time_n   = nodes->GetSolution_time_n(iPoint);
       V_time_nP1 = nodes->GetSolution(iPoint);
 
-      /*--- Access the density at this node (constant for now). ---*/
+      /*--- Access the density at this node for each time level. ---*/
 
-      Density = nodes->GetDensity(iPoint);
+      Density_nP1 = nodes->GetDensity(iPoint);
+      Density_n   = nodes->GetDensity_time_n(iPoint);
+      Density_nM1 = nodes->GetDensity_time_n1(iPoint);
 
       /*--- Compute the conservative variable vector for all time levels. ---*/
 
-      V2U(Density, V_time_nM1, U_time_nM1);
-      V2U(Density, V_time_n, U_time_n);
-      V2U(Density, V_time_nP1, U_time_nP1);
+      V2U(Density_nM1, V_time_nM1, U_time_nM1);
+      V2U(Density_n, V_time_n, U_time_n);
+      V2U(Density_nP1, V_time_nP1, U_time_nP1);
 
       /*--- CV volume at time n+1. As we are on a static mesh, the volume
        of the CV will remained fixed for all time steps. ---*/
@@ -2870,7 +2872,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       /*--- Compute the Jacobian contribution due to the dual time source term. ---*/
 
       if (implicit) {
-        su2double delta = (second_order ? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep;
+        su2double delta = (second_order ? 1.5 : 1.0) * Volume_nP1 * Density_nP1 / TimeStep;
 
         for (iVar = 1; iVar < nVar; ++iVar) Jacobian.AddVal2Diag(iPoint, iVar, delta);
       }
@@ -2896,8 +2898,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       /*--- Compute the conservative variables. ---*/
 
       V_time_n = nodes->GetSolution_time_n(iPoint);
-      Density = nodes->GetDensity(iPoint);
-      V2U(Density, V_time_n, U_time_n);
+      Density_n = nodes->GetDensity_time_n(iPoint);
+      V2U(Density_n, V_time_n, U_time_n);
 
       GridVel_i = geometry->nodes->GetGridVel(iPoint);
 
@@ -2951,8 +2953,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
           /*--- Compute the GCL component of the source term for node i ---*/
 
           V_time_n = nodes->GetSolution_time_n(iPoint);
-          Density = nodes->GetDensity(iPoint);
-          V2U(Density, V_time_n, U_time_n);
+          Density_n = nodes->GetDensity_time_n(iPoint);
+          V2U(Density_n, V_time_n, U_time_n);
 
           for (iVar = 0; iVar < nVar-!energy; iVar++)
             LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL;
@@ -2978,15 +2980,17 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       V_time_n   = nodes->GetSolution_time_n(iPoint);
       V_time_nP1 = nodes->GetSolution(iPoint);
 
-      /*--- Access the density at this node (constant for now). ---*/
+      /*--- Access the density at this node for each time level. ---*/
 
-      Density = nodes->GetDensity(iPoint);
+      Density_nP1 = nodes->GetDensity(iPoint);
+      Density_n   = nodes->GetDensity_time_n(iPoint);
+      Density_nM1 = nodes->GetDensity_time_n1(iPoint);
 
       /*--- Compute the conservative variable vector for all time levels. ---*/
 
-      V2U(Density, V_time_nM1, U_time_nM1);
-      V2U(Density, V_time_n, U_time_n);
-      V2U(Density, V_time_nP1, U_time_nP1);
+      V2U(Density_nM1, V_time_nM1, U_time_nM1);
+      V2U(Density_n, V_time_n, U_time_n);
+      V2U(Density_nP1, V_time_nP1, U_time_nP1);
 
       /*--- CV volume at time n-1 and n+1. In the case of dynamically deforming
        grids, the volumes will change. On rigidly transforming grids, the
@@ -3010,7 +3014,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       /*--- Compute the Jacobian contribution due to the dual time source term. ---*/
 
       if (implicit) {
-        su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep;
+        su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density_nP1 / TimeStep;
 
         for (iVar = 1; iVar < nVar; ++iVar)
           Jacobian.AddVal2Diag(iPoint, iVar, delta);

SU2_CFD/src/variables/CIncEulerVariable.cpp

@@ -58,6 +58,12 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci
   if (dual_time) {
     Solution_time_n = Solution;
     Solution_time_n1 = Solution;
+
+    /*--- Allocate density storage for previous time levels.
+     *    This is needed for variable-density incompressible unsteady flows
+     *    (e.g. with energy equation or species transport). ---*/
+    Density_time_n.resize(nPoint) = su2double(0.0);
+    Density_time_n1.resize(nPoint) = su2double(0.0);
   }
 
   if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) {
@@ -122,8 +128,22 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
   /*--- Set enthalpy ---*/
 
-  SetEnthalpy(iPoint, FluidModel->GetEnthalpy());
+    SetEnthalpy(iPoint, FluidModel->GetEnthalpy());
 
   return physical;
 
 }
+
+void CIncEulerVariable::Set_Density_time_n() {
+  const auto nPoint = Density_time_n.size();
+  SU2_OMP_FOR_STAT(roundUpDiv(nPoint, omp_get_num_threads()))
+  for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
+    Density_time_n(iPoint) = GetDensity(iPoint);
+  }
+  END_SU2_OMP_FOR
+}
+
+void CIncEulerVariable::Set_Density_time_n1() {
+  assert(Density_time_n1.size() == Density_time_n.size());
+  parallelCopy(Density_time_n.size(), Density_time_n.data(), Density_time_n1.data());
+}