Implementing variable density for unsteady incompressible flow

Common/src/CConfig.cpp

@@ -5764,6 +5764,14 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
         Kind_Solver != MAIN_SOLVER::MULTIPHYSICS)
       SU2_MPI::Error("Species transport currently only available for compressible and incompressible flow.", CURRENT_FUNCTION);
 
+    /*--- The dual-time density history is recomputed via the fluid model, which needs the species
+          solution; the species solver only exists on the finest grid. ---*/
+    if ((Kind_Regime == ENUM_REGIME::INCOMPRESSIBLE) && (Kind_DensityModel != INC_DENSITYMODEL::CONSTANT) &&
+        (TimeMarching == TIME_MARCHING::DT_STEPPING_1ST || TimeMarching == TIME_MARCHING::DT_STEPPING_2ND) &&
+        (nMGLevels > 0))
+      SU2_MPI::Error("Dual-time stepping with species-dependent variable density does not support MGLEVEL > 0.",
+                     CURRENT_FUNCTION);
+
     /*--- Species specific OF currently can only handle one entry in Marker_Analyze. ---*/
     for (unsigned short iObj = 0; iObj < nObj; iObj++) {
       if ((Kind_ObjFunc[iObj] == SURFACE_SPECIES_0 ||

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -86,6 +86,18 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
    */
   virtual unsigned long SetPrimitive_Variables(CSolver **solver_container, const CConfig *config);
 
+  /*!
+   * \brief Recompute the dual-time density history (rho_n, rho_n-1) from the stored
+   *        solution histories (Solution_time_n[,1], and the species histories for
+   *        flamelet/variable-density) via the fluid model.
+   * \details Called once per physical time step (InnerIter == 0), after the solution
+   *          push-back. This keeps the density history consistent with the primitive
+   *          state at time n (also after restart).
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void RecomputeDensity_time_n(CSolver **solver_container, const CConfig *config);
+
   /*!
    * \brief Update the Beta parameter for the incompressible preconditioner.
    * \param[in] geometry - Geometrical definition of the problem.

SU2_CFD/include/solvers/CScalarSolver.hpp

@@ -58,6 +58,7 @@ class CScalarSolver : public CSolver {
   su2double Solution_Inf[MAXNVAR]; /*!< \brief Far-field solution. */
 
   const bool Conservative; /*!< \brief Transported Variable is conservative. Solution has to be multiplied with rho. */
+  const bool BoundedScalar; /*!< \brief Whether the derived solver uses the bounded-scalar convective scheme. */
 
   const CPrimitiveIndices<unsigned short> prim_idx; /*!< \brief Indices of the primitive flow variables. */
 
@@ -440,7 +441,7 @@ class CScalarSolver : public CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative);
+  CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative, bool bounded_scalar);
 
   /*!
    * \brief Compute the spatial integration using a upwind scheme.

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -30,8 +30,8 @@
 #include "../../include/variables/CFlowVariable.hpp"
 
 template <class VariableType>
-CScalarSolver<VariableType>::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative)
-    : CSolver(), Conservative(conservative),
+CScalarSolver<VariableType>::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative, bool bounded_scalar)
+    : CSolver(), Conservative(conservative), BoundedScalar(bounded_scalar),
       prim_idx(config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE,
                config->GetNEMOProblem(), geometry->GetnDim(), config->GetnSpecies()) {
   SU2_ZONE_SCOPED
@@ -639,11 +639,11 @@ 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);
 
-  /*--- Flow solution, needed to get density. ---*/
+  const bool bounded_scalar = BoundedScalar;
 
-  CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes();
+  /*--- Flow solution, needed to get density. ---*/
+  auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
   /*--- Store the physical time step ---*/
 
@@ -671,19 +671,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
@@ -703,13 +693,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. ---*/
@@ -737,10 +734,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++) {
@@ -793,10 +787,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;
@@ -832,37 +823,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

SU2_CFD/include/solvers/CTurbSolver.hpp

@@ -133,7 +133,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/CEulerVariable.hpp

@@ -97,6 +97,20 @@ class CEulerVariable : public CFlowVariable {
   CEulerVariable(su2double density, const su2double *velocity, su2double energy,
                  unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config);
 
+  /*!
+   * \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 final { return GetSolution_time_n(iPoint, 0); }
+
+  /*!
+   * \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 final { return GetSolution_time_n1(iPoint, 0); }
+
   /*!
    * \brief A virtual member.
    */

SU2_CFD/include/variables/CFlowVariable.hpp

@@ -270,4 +270,18 @@ 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.
+   * \param[in] iPoint - Point index.
+   * \return Density at time level n.
+   */
+  virtual su2double GetDensity_time_n(unsigned long iPoint) const = 0;
+
+  /*!
+   * \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.
+   */
+  virtual su2double GetDensity_time_n1(unsigned long iPoint) const = 0;
 };

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:
   /*!
@@ -291,4 +293,36 @@ 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 final {
+    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 final {
+    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; }
+
 };

SU2_CFD/include/variables/CNEMOEulerVariable.hpp

@@ -122,6 +122,20 @@ class CNEMOEulerVariable : public CFlowVariable {
                      unsigned long nvar, unsigned long nvalprim,
                      unsigned long nvarprimgrad, const CConfig *config, CNEMOGas *fluidmodel);
 
+  /*!
+   * \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 final { return GetSolution_time_n(iPoint, 0); }
+
+  /*!
+   * \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 final { return GetSolution_time_n1(iPoint, 0); }
+
   /*---------------------------------------*/
   /*---          U,V,S Routines         ---*/
   /*---------------------------------------*/

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -168,11 +168,13 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) {
   /*--- Register quantities that are no solver variables but further inputs/outputs of the (outer) iteration. ---*/
   direct_solver->RegisterSolutionExtra(true, config);
 
-  if (time_n_needed)
+  if (time_n_needed) {
     direct_solver->GetNodes()->RegisterSolution_time_n();
+  }
 
-  if (time_n1_needed)
+  if (time_n1_needed) {
     direct_solver->GetNodes()->RegisterSolution_time_n1();
+  }
 }
 
 void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, bool reset) {

SU2_CFD/src/solvers/CHeatSolver.cpp

@@ -34,7 +34,7 @@
 template class CScalarSolver<CHeatVariable>;
 
 CHeatSolver::CHeatSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh)
-  : CScalarSolver<CHeatVariable>(geometry, config, false),
+  : CScalarSolver<CHeatVariable>(geometry, config, false, false),
     flow(config->GetFluidProblem()) {
   SU2_ZONE_SCOPED