[WIP] Implementing variable density for unsteady incompressible flow

SU2_CFD/include/output/COutput.hpp

@@ -300,7 +300,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"};
 
   /*----------------------------- Convergence monitoring ----------------------------*/
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1069,6 +1069,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_Solution_time_n1();
+    solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Density_time_n();
     if (rans) {
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n1();
@@ -1099,6 +1100,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

@@ -633,7 +633,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
           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_n = flowNodes->GetDensity_time_n(iPoint);
           Density_nP1 = flowNodes->GetDensity(iPoint);
         } else {
           Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0];
@@ -694,7 +694,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);  // Updated for transient density
         else
           Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
       }
@@ -750,7 +750,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);  // Updated for transient density
             else
               Density_n = flowNodes->GetSolution_time_n(iPoint, 0);
           }
@@ -796,7 +796,7 @@ void CScalarSolver<VariableType>::SetResidual_DualTime(CGeometry* geometry, CSol
           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_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

@@ -184,6 +184,14 @@ class CEulerVariable : public CFlowVariable {
     return Primitive(iPoint, indices.Density()) <= 0.0;
   }
 
+  inline void Set_Density_time_n(unsigned long iPoint, su2double val) {
+    Density_time_n[iPoint] = val;
+  }
+
+  inline void Set_Density_unsteady(unsigned long iPoint, su2double val) {
+    Density_unsteady[iPoint] = val;
+  }
+
   /*!
    * \brief Set the value of the temperature.
    * \param[in] temperature - how agitated the particles are :)

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -82,7 +82,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 temperature,
+  CIncEulerVariable(su2double density, su2double pressure, const su2double *velocity, su2double temperature,
                     unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config);
 
   /*!
@@ -99,7 +99,19 @@ class CIncEulerVariable : public CFlowVariable {
     Primitive(iPoint, indices.Density()) = val_density;
     return val_density <= 0.0;
   }
+
+  inline void Set_Density_time_n(unsigned long iPoint, su2double val) {
+    Density_time_n[iPoint] = val;
+  }
 
+  inline void Set_Density_unsteady(unsigned long iPoint, su2double val) {
+    Density_unsteady[iPoint] = val;
+  }
+
+  inline su2double GetDensity_time_n(unsigned long iPoint) const {
+    return Density_time_n[iPoint];
+  }
+
   /*!
    * \brief Set the value of the density for the incompressible flows.
    * \param[in] iPoint - Point index.

SU2_CFD/include/variables/CIncNSVariable.hpp

@@ -52,7 +52,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

@@ -68,6 +68,8 @@ 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;
+  VectorType Density_time_n;         /*!< \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. */
@@ -275,6 +277,11 @@ class CVariable {
    */
   void Set_Solution_time_n1();
 
+    /*!
+   * \brief Set the density at time n.
+   */
+  void Set_Density_time_n();
+
   /*!
    * \brief Set the variable solution at time n.
    * \param[in] iPoint - Point index.
@@ -308,6 +315,14 @@ class CVariable {
     Solution_time_n1(iPoint,iVar) = val_sol;
   }
 
+  inline void Set_Density_time_n(unsigned long iPoint, su2double val) {
+    Density_time_n[iPoint] = val;
+  }
+
+  inline void Set_Density_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.
@@ -502,6 +517,22 @@ 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 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.
@@ -2154,6 +2185,11 @@ class CVariable {
    */
   void RegisterSolution_time_n1();
 
+   /*!
+   * \brief Register the variables in the solution_time_n1 array as input/output variable.
+   */
+  void RegisterDensity_time_n();
+
   /*!
    * \brief Set the adjoint values of the solution.
    * \param[in] adj_sol - The adjoint values of the solution.

SU2_CFD/src/integration/CIntegration.cpp

@@ -248,6 +248,7 @@ 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()->Set_Density_time_n();
 
   SU2_OMP_SAFE_GLOBAL_ACCESS(solver->ResetCFLAdapt();)
 

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -375,6 +375,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
 
   AddCommonFVMOutputs(config);
 
+  AddVolumeOutput("DENSITY_TIME_N", "Density_time_n", "SOLUTION", "Density at previous time step");
   if (config->GetTime_Domain()) {
     SetTimeAveragedFields();
   }
@@ -418,6 +419,7 @@ 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));
+  SetVolumeOutputValue("DENSITY_TIME_N", iPoint, Node_Flow->GetDensity_time_n(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));

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -208,9 +208,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, Temperature_Inf, nPoint, nDim, nVar, config);
+    nodes = new CIncNSVariable(Density_Inf, Pressure_Inf, Velocity_Inf, Temperature_Inf, nPoint, nDim, nVar, config);
   } else {
-    nodes = new CIncEulerVariable(Pressure_Inf, Velocity_Inf, Temperature_Inf, nPoint, nDim, nVar, config);
+    nodes = new CIncEulerVariable(Density_Inf, Pressure_Inf, Velocity_Inf, Temperature_Inf, nPoint, nDim, nVar, config);
   }
   SetBaseClassPointerToNodes();
 
@@ -2670,7 +2670,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, Cp;
+  su2double Density, Cp, Density_time_n, Density_unsteady;
 
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST);
@@ -2710,13 +2710,15 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
 
       /*--- Access the density and Cp at this node (constant for now). ---*/
 
+      Density_time_n = nodes->GetDensity_time_n(iPoint);
+      Density_unsteady = nodes->GetDensity_unsteady(iPoint);
       Density = nodes->GetDensity(iPoint);
       Cp = nodes->GetSpecificHeatCp(iPoint);
-
+      
       /*--- Compute the conservative variable vector for all time levels. ---*/
 
       V2U(Density, Cp, V_time_nM1, U_time_nM1);
-      V2U(Density, Cp, V_time_n, U_time_n);
+      V2U(Density_time_n, Cp, V_time_n, U_time_n);
       V2U(Density, Cp, V_time_nP1, U_time_nP1);
 
       /*--- CV volume at time n+1. As we are on a static mesh, the volume

SU2_CFD/src/variables/CIncEulerVariable.cpp

@@ -28,7 +28,7 @@
 #include "../../include/variables/CIncEulerVariable.hpp"
 #include "../../include/fluid/CFluidModel.hpp"
 
-CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *velocity, su2double temperature,
+CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, const su2double *velocity, su2double temperature,
                                      unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config)
   : CFlowVariable(npoint, ndim, nvar, ndim + 9,
                   ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED ? 2 : 4), config),
@@ -46,7 +46,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci
   for(unsigned long iPoint=0; iPoint<nPoint; ++iPoint)
     for (unsigned long iVar = 0; iVar < nVar; iVar++)
       Solution(iPoint,iVar) = val_solution[iVar];
-
+  Density_unsteady = density;
   Solution_Old = Solution;
 
   if (classical_rk4) Solution_New = Solution;
@@ -56,6 +56,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci
   if (dual_time) {
     Solution_time_n = Solution;
     Solution_time_n1 = Solution;
+    Density_unsteady = density;
   }
 
   if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) {
@@ -65,8 +66,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci
   }
 }
 
-bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) {
-
+bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel)  {
   bool physical = true;
 
   /*--- Set the value of the pressure ---*/
@@ -89,6 +89,7 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
   /*--- Set the value of the density ---*/
 
   const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity());
+  Density_unsteady[iPoint] = FluidModel->GetDensity();
 
   /*--- Non-physical solution found. Revert to old values. ---*/
 

SU2_CFD/src/variables/CIncNSVariable.cpp

@@ -28,9 +28,9 @@
 #include "../../include/variables/CIncNSVariable.hpp"
 #include "../../include/fluid/CFluidModel.hpp"
 
-CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su2double temperature,
+CIncNSVariable::CIncNSVariable(su2double density, su2double pressure, const su2double *velocity, su2double temperature,
                                unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) :
-                               CIncEulerVariable(pressure, velocity, temperature, npoint, ndim, nvar, config) {
+                               CIncEulerVariable(density, pressure, velocity, temperature, npoint, ndim, nvar, config) {
 
   Vorticity.resize(nPoint,3);
   StrainMag.resize(nPoint);
@@ -78,6 +78,7 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do
   /*--- Set the value of the density ---*/
 
   const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity());
+  Density_unsteady[iPoint] = FluidModel->GetDensity();
 
   /*--- Non-physical solution found. Revert to old values. ---*/
 

SU2_CFD/src/variables/CVariable.cpp

@@ -39,6 +39,9 @@ CVariable::CVariable(unsigned long npoint, unsigned long nvar, const CConfig *co
   /*--- Allocate the solution array. ---*/
   Solution.resize(nPoint,nVar) = su2double(0.0);
 
+  Density_unsteady.resize(nPoint);
+  Density_time_n.resize(nPoint);
+
   if (config->GetMultizone_Problem())
     Solution_BGS_k.resize(nPoint,nVar) = su2double(0.0);
 
@@ -55,15 +58,19 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva
   /*--- Allocate fields common to all problems. Do not allocate fields
    that are specific to one solver, i.e. not common, in this class. ---*/
   Solution.resize(nPoint,nVar) = su2double(0.0);
-
   Solution_Old.resize(nPoint,nVar) = su2double(0.0);
 
   if (config->GetTime_Domain())
     Solution_time_n.resize(nPoint,nVar) = su2double(0.0);
 
+  if (config->GetTime_Domain()) {
+    Density_unsteady.resize(nPoint);
+    Density_time_n.resize(nPoint);
+  }
+
   if (config->GetTime_Marching() != TIME_MARCHING::STEADY)
     Solution_time_n1.resize(nPoint,nVar) = su2double(0.0);
-
+    
   if (config->GetDiscrete_Adjoint()) {
     if (adjoint && config->GetMultizone_Problem())
       External.resize(nPoint,nVar) = su2double(0.0);
@@ -93,6 +100,11 @@ void CVariable::Set_Solution_time_n() {
   parallelCopy(Solution.size(), Solution.data(), Solution_time_n.data());
 }
 
+void CVariable::Set_Density_time_n(){
+  assert(Density_time_n.size() == Density_unsteady.size());
+  parallelCopy(Density_unsteady.size(),Density_unsteady.data(), Density_time_n.data());
+}
+
 void CVariable::Set_Solution_time_n1() {
   assert(Solution_time_n1.size() == Solution_time_n.size());
   parallelCopy(Solution_time_n.size(), Solution_time_n.data(), Solution_time_n1.data());