[WIP] Implementing variable density for unsteady incompressible flow

SU2_CFD/include/output/COutput.hpp

@@ -302,7 +302,7 @@ class COutput {
   unsigned short nRequestedVolumeFields;
 
   /*! \brief Minimum required volume fields for restart file. */
-  const std::vector<string> restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY"};
+  const std::vector<string> restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY", "DENSITY_TIME_N", "DENSITY_TIME_N1"};
 
   /*----------------------------- Convergence monitoring ----------------------------*/
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1071,9 +1071,13 @@ 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();
+    // nijso asks: we are in first order here: can this be removed?
     solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n1();
+    solver_container[iMesh][FLOW_SOL]->GetNodes()->SetDensity_time_n();
+    solver_container[iMesh][FLOW_SOL]->GetNodes()->SetDensity_time_n1();
     if (rans) {
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
+      // nijso asks: we are in first order here: can this be removed?
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n1();
     }
 
@@ -1101,7 +1105,10 @@ void CFVMFlowSolverBase<V, R>::PushSolutionBackInTime(unsigned long TimeIter, bo
     /*--- Push back this new solution to time level N. ---*/
 
     for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
+      // nijso asks: now we are at second order, but we push solution to time_n instead of time_n1 ?
       solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n();
+      solver_container[iMesh][FLOW_SOL]->GetNodes()->SetDensity_time_n();
+      solver_container[iMesh][FLOW_SOL]->GetNodes()->SetDensity_time_n1();
       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,8 @@ 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_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,9 +697,11 @@ 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);  // Updated for transient density
+          // nijso: no second order for dynamic meshes?
         else
           Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
+          // nijso: no second order for dynamic meshes?
       }
 
       for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) {
@@ -757,11 +755,11 @@ 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);  // Updated for transient density
             else
               Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
           }
-
+          // nijso: what about second order?
           for (iVar = 0; iVar < nVar; iVar++) LinSysRes(iPoint, iVar) += Density_n * U_time_n[iVar] * Residual_GCL;
         }
         END_SU2_OMP_FOR
@@ -802,8 +800,8 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
           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_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/CEulerVariable.hpp

@@ -189,6 +189,18 @@ class CEulerVariable : public CFlowVariable {
     return Primitive(iPoint, indices.Density()) <= 0.0;
   }
 
+  inline void SetDensity_Time_n(unsigned long iPoint, su2double val_density_time_n) {
+    Density_time_n[iPoint] = val_density_time_n;
+  }
+
+  inline void SetDensity_time_n1(unsigned long iPoint, su2double val_density_time_n1) {
+    Density_time_n1[iPoint] = val_density_time_n1;
+  }
+
+  inline void SetDensity_Unsteady(unsigned long iPoint, su2double val_density_unsteady) {
+    Density_unsteady[iPoint] = val_density_unsteady;
+  }
+
   /*!
    * \brief Set the value of the temperature.
    * \param[in] temperature - how agitated the particles are :)

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -83,8 +83,7 @@ class CIncEulerVariable : public CFlowVariable {
    * \param[in] nvar - Number of variables of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-
-  CIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy,
+  CIncEulerVariable(su2double density, su2double pressure, const su2double *velocity, su2double enthalpy,
                     unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config);
 
   /*!
@@ -101,6 +100,26 @@ class CIncEulerVariable : public CFlowVariable {
     Primitive(iPoint, indices.Density()) = val_density;
     return val_density <= 0.0;
   }
+
+  inline void SetDensity_time_n(unsigned long iPoint, su2double val) {
+    Density_time_n[iPoint] = val;
+  }
+
+  inline void SetDensity_time_n1(unsigned long iPoint, su2double val) {
+    Density_time_n1[iPoint] = val;
+  }
+
+  inline void SetDensity_Unsteady(unsigned long iPoint, su2double val) {
+    Density_unsteady[iPoint] = val;
+  }
+
+  inline su2double GetDensity_time_n(unsigned long iPoint) const {
+    return Density_time_n[iPoint];
+  }
+
+  inline su2double GetDensity_time_n1(unsigned long iPoint) const {
+    return Density_time_n1[iPoint];
+  }
 
   /*!
    * \brief Set the value of the density for the incompressible flows.

SU2_CFD/include/variables/CIncNSVariable.hpp

@@ -41,6 +41,7 @@ class CIncNSVariable final : public CIncEulerVariable {
   VectorType Tau_Wall;        /*!< \brief Magnitude of the wall shear stress from a wall function. */
   VectorType DES_LengthScale;
   const bool Energy;          /*!< \brief Flag for Energy equation in incompressible flows. */
+  bool Unsteady = false; /*!< \brief Flag for unsteady incompressible flows. */
 
 public:
   /*!
@@ -53,7 +54,7 @@ class CIncNSVariable final : public CIncEulerVariable {
    * \param[in] nvar - Number of variables of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  CIncNSVariable(su2double pressure, const su2double *velocity, su2double temperature,
+  CIncNSVariable(su2double density, su2double pressure, const su2double *velocity, su2double temperature,
                  unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config);
 
   /*!

SU2_CFD/include/variables/CVariable.hpp

@@ -70,6 +70,9 @@ class CVariable {
   MatrixType Solution_time_n;    /*!< \brief Solution of the problem at time n for dual-time stepping technique. */
   MatrixType Solution_time_n1;   /*!< \brief Solution of the problem at time n-1 for dual-time stepping technique. */
   VectorType Delta_Time;         /*!< \brief Time step. */
+  VectorType Density_unsteady;       /*!< \brief density for unsteady flows. */
+  VectorType Density_time_n;         /*!< \brief density at time n for dual-time stepping technique. */
+  VectorType Density_time_n1;        /*!< \brief density at time n for dual-time stepping technique. */
 
   CVectorOfMatrix Gradient;  /*!< \brief Gradient of the solution of the problem. */
   C3DDoubleMatrix Rmatrix;   /*!< \brief Geometry-based matrix for weighted least squares gradient calculations. */
@@ -291,6 +294,16 @@ class CVariable {
    */
   void Set_Solution_time_n1();
 
+    /*!
+   * \brief Set the density at time n.
+   */
+  void SetDensity_time_n();
+
+    /*!
+   * \brief Set the density at time n-1.
+   */
+  void SetDensity_time_n1();
+
   /*!
    * \brief Set the variable solution at time n.
    * \param[in] iPoint - Point index.
@@ -324,6 +337,18 @@ class CVariable {
     Solution_time_n1(iPoint,iVar) = val_sol;
   }
 
+  inline void SetDensity_time_n(unsigned long iPoint, su2double val) {
+    Density_time_n[iPoint] = val;
+  }
+
+   inline void SetDensity_time_n1(unsigned long iPoint, su2double val) {
+    Density_time_n1[iPoint] = val;
+  }
+
+  inline void SetDensity_unsteady(unsigned long iPoint, su2double val) {
+    Density_unsteady[iPoint] = val;
+  } 
+
   /*!
    * \brief Virtual Member. Specify a vector to set the velocity components of the solution.
    *        Multiplied by density for compressible cases.
@@ -525,6 +550,30 @@ class CVariable {
   inline su2double *GetSolution_time_n1(unsigned long iPoint) { return Solution_time_n1[iPoint]; }
   inline MatrixType& GetSolution_time_n1() { return Solution_time_n1; }
 
+    /*!
+   * \brief Get the solution at time n.
+   * \param[in] iPoint - Point index.
+   * \return Pointer to the solution (at time n) vector.
+   */
+  inline su2double GetDensity_time_n(unsigned long iPoint) const { return Density_time_n[iPoint]; }
+  inline VectorType& GetDensity_time_n() { return Density_time_n; }
+
+  /*!
+   * \brief Get the solution at time n.
+   * \param[in] iPoint - Point index.
+   * \return Pointer to the solution (at time n) vector.
+   */
+  inline su2double GetDensity_time_n1(unsigned long iPoint) const { return Density_time_n1[iPoint]; }
+  inline VectorType& GetDensity_time_n1() { return Density_time_n1; }
+
+  /*!
+   * \brief Get the density.
+   * \param[in] iPoint - Point index.
+   * \return Pointer to the solution (at time n) vector.
+   */
+  inline su2double GetDensity_unsteady(unsigned long iPoint) const { return Density_unsteady[iPoint]; }
+  inline VectorType& GetDensity_unsteady() { return Density_unsteady; }
+
   /*!
    * \brief Set the value of the old residual.
    * \param[in] iPoint - Point index.
@@ -2177,6 +2226,15 @@ class CVariable {
    */
   void RegisterSolution_time_n1();
 
+   /*!
+   * \brief Register the variables in the density_time_n array as input/output variable.
+   */
+  void RegisterDensity_time_n();
+   /*!
+   * \brief Register the variables in the density_time_n1 array as input/output variable.
+   */
+  void RegisterDensity_time_n1();
+
   /*!
    * \brief Register the variables in the user defined source array as input/output variable.
    */

SU2_CFD/src/integration/CIntegration.cpp

@@ -248,6 +248,8 @@ void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver
   /*--- Store old solution ---*/
   solver->GetNodes()->Set_Solution_time_n1();
   solver->GetNodes()->Set_Solution_time_n();
+  solver->GetNodes()->SetDensity_time_n();
+  solver->GetNodes()->SetDensity_time_n1();
 
   SU2_OMP_SAFE_GLOBAL_ACCESS(solver->ResetCFLAdapt();)
 

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -390,7 +390,10 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
   }
 
   AddCommonFVMOutputs(config);
-
+  if (config->GetTime_Marching() != TIME_MARCHING::STEADY) {
+    AddVolumeOutput("DENSITY_TIME_N", "Density_time_n", "SOLUTION", "Density at previous time step n");
+    AddVolumeOutput("DENSITY_TIME_N1", "Density_time_n1", "SOLUTION", "Density at previous time step n-1");
+  }
   if (config->GetTime_Domain()) {
     SetTimeAveragedFields();
   }
@@ -436,7 +439,10 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
   const su2double factor = solver[FLOW_SOL]->GetReferenceDynamicPressure();
   SetVolumeOutputValue("PRESSURE_COEFF", iPoint, (Node_Flow->GetPressure(iPoint) - solver[FLOW_SOL]->GetPressure_Inf())/factor);
   SetVolumeOutputValue("DENSITY", iPoint, Node_Flow->GetDensity(iPoint));
-
+  if (config->GetTime_Marching() != TIME_MARCHING::STEADY) {
+    SetVolumeOutputValue("DENSITY_TIME_N", iPoint, Node_Flow->GetDensity_time_n(iPoint));
+    SetVolumeOutputValue("DENSITY_TIME_N1", iPoint, Node_Flow->GetDensity_time_n1(iPoint));
+  }
   if (config->GetKind_Solver() == MAIN_SOLVER::INC_RANS || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES){
     SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint));
     SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint));

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -200,9 +200,9 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
   /*--- Initialize the solution to the far-field state everywhere. ---*/
 
   if (navier_stokes) {
-    nodes = new CIncNSVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config);
+    nodes = new CIncNSVariable(Density_Inf, Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config);
   } else {
-    nodes = new CIncEulerVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config);
+    nodes = new CIncEulerVariable(Density_Inf, Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config);
   }
   SetBaseClassPointerToNodes();
 
@@ -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, Density_time_n, Density_time_nM1;
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
@@ -2843,12 +2843,14 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
 
       /*--- Access the density at this node (constant for now). ---*/
 
+      Density_time_nM1 = nodes->GetDensity_time_n1(iPoint);
+      Density_time_n = nodes->GetDensity_time_n(iPoint);
       Density = nodes->GetDensity(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_time_nM1, V_time_nM1, U_time_nM1);
+      V2U(Density_time_n, V_time_n, U_time_n);
       V2U(Density, V_time_nP1, U_time_nP1);
 
       /*--- CV volume at time n+1. As we are on a static mesh, the volume
@@ -2897,7 +2899,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
 
       V_time_n = nodes->GetSolution_time_n(iPoint);
       Density = nodes->GetDensity(iPoint);
-      V2U(Density, V_time_n, U_time_n);
+Density_time_n = nodes->GetDensity_time_n(iPoint);
+V2U(Density_time_n, V_time_n, U_time_n);
 
       GridVel_i = geometry->nodes->GetGridVel(iPoint);
 
@@ -2981,12 +2984,16 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       /*--- Access the density at this node (constant for now). ---*/
 
       Density = nodes->GetDensity(iPoint);
+      Density_time_nM1 = nodes->GetDensity_time_n1(iPoint);
+      Density_time_n = nodes->GetDensity_time_n(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_time_nM1, V_time_nM1, U_time_nM1);
+V2U(Density_time_n, V_time_n, U_time_n);
+V2U(Density, 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

SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp

@@ -88,6 +88,7 @@
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
 
   SU2_OMP_FOR_STAT(omp_chunk_size)
+  //for (auto i_point = 0u; i_point < nPointDomain; i_point++) {
   for (auto i_point = 0u; i_point < nPoint; i_point++) {
     CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel();
     su2double* scalars = nodes->GetSolution(i_point);