[WIP] Implementing variable density for unsteady incompressible flow

Common/src/CConfig.cpp

@@ -5792,6 +5792,22 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Number of initial species incompatible with number of controlling variables and user scalars.", CURRENT_FUNCTION);
     /*--- We can have additional user defined transported scalars ---*/
     flamelet_ParsedOptions.n_scalars = flamelet_ParsedOptions.n_control_vars + flamelet_ParsedOptions.n_user_scalars;
+
+    /*--- Check that spark ignition has required parameters defined ---*/
+    if (flamelet_ParsedOptions.ignition_method == FLAMELET_INIT_TYPE::SPARK) {
+      /*--- Check if SPARK_INIT was explicitly set in config file ---*/
+      if (all_options.find("SPARK_INIT") != all_options.end()) {
+        SU2_MPI::Error("FLAME_INIT_METHOD=SPARK requires SPARK_INIT to be defined in the config file.", CURRENT_FUNCTION);
+      }
+      /*--- Check if SPARK_REACTION_RATES was explicitly set in config file ---*/
+      if (all_options.find("SPARK_REACTION_RATES") != all_options.end()) {
+        SU2_MPI::Error("FLAME_INIT_METHOD=SPARK requires SPARK_REACTION_RATES to be defined in the config file.", CURRENT_FUNCTION);
+      }
+      if (flamelet_ParsedOptions.nspark < flamelet_ParsedOptions.n_scalars) {
+        SU2_MPI::Error("SPARK_REACTION_RATES must have at least " + to_string(flamelet_ParsedOptions.n_scalars) +
+                       " values (one for each scalar variable), but only " + to_string(flamelet_ParsedOptions.nspark) + " were provided.", CURRENT_FUNCTION);
+      }
+    }
   }
 
   if (Kind_Regime == ENUM_REGIME::COMPRESSIBLE && GetBounded_Scalar()) {

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,7 +697,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,11 +753,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);
             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
@@ -798,12 +794,8 @@ 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);
+          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];
@@ -863,4 +855,4 @@ void CScalarSolver<VariableType>::PushSolutionBackInTime(unsigned long TimeIter,
       nodes->Set_Solution_time_n();
     }
   }
-}
+}

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,68 @@ 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 { return Density_time_n(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 { return Density_time_n1(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 the density at time n as an AD input variable (for discrete adjoint unsteady).
+   */
+  void RegisterDensity_time_n() override;
+
+  /*!
+   * \brief Register the density at time n-1 as an AD input variable (for discrete adjoint unsteady).
+   */
+  void RegisterDensity_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 {
+    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 {
+    AD::Identifier index = AD::GetPassiveIndex();
+    AD::SetIndex(index, Density_time_n1(iPoint));
+    return AD::GetDerivative(index);
+  }
+
 };

SU2_CFD/include/variables/CVariable.hpp

@@ -284,12 +284,12 @@ class CVariable {
   /*!
    * \brief Set the variable solution at time n.
    */
-  void Set_Solution_time_n();
+  virtual void Set_Solution_time_n();
 
   /*!
    * \brief Set the variable solution at time n-1.
    */
-  void Set_Solution_time_n1();
+  virtual void Set_Solution_time_n1();
 
   /*!
    * \brief Set the variable solution at time n.
@@ -2182,6 +2182,48 @@ class CVariable {
    */
   void RegisterUserDefinedSource();
 
+  /*!
+   * \brief Get the stored density at time n (incompressible only).
+   * \param[in] iPoint - Point index.
+   * \return Density at time n; 0 for compressible nodes.
+   */
+  inline virtual su2double GetDensity_time_n(unsigned long /*iPoint*/) const { return 0.0; }
+
+  /*!
+   * \brief Get the stored density at time n-1 (incompressible only).
+   * \param[in] iPoint - Point index.
+   * \return Density at time n-1; 0 for compressible nodes.
+   */
+  inline virtual su2double GetDensity_time_n1(unsigned long /*iPoint*/) const { return 0.0; }
+
+  /*!
+   * \brief Register the density at time n as an AD input (no-op for non-incompressible nodes).
+   */
+  inline virtual void RegisterDensity_time_n() {}
+
+  /*!
+   * \brief Register the density at time n-1 as an AD input (no-op for non-incompressible nodes).
+   */
+  inline virtual void RegisterDensity_time_n1() {}
+
+  /*!
+   * \brief Get the adjoint values of the density at time n.
+   *        Implemented analogously to GetAdjointSolution_time_n but for a scalar
+   *        (Density_time_n is a VectorType, not a MatrixType).
+   *        Currently ignores the chain-rule d(rho)/d(h) dependence (TODO).
+   * \param[in] iPoint - Point index.
+   * \return Adjoint value of the density at time n.
+   */
+  inline virtual su2double GetAdjointDensity_time_n(unsigned long /*iPoint*/) const { return 0.0; }
+
+  /*!
+   * \brief Get the adjoint values of the density at time n-1.
+   *        See GetAdjointDensity_time_n.
+   * \param[in] iPoint - Point index.
+   * \return Adjoint value of the density at time n-1.
+   */
+  inline virtual su2double GetAdjointDensity_time_n1(unsigned long /*iPoint*/) const { return 0.0; }
+
   /*!
    * \brief Set the adjoint values of the solution.
    * \param[in] adj_sol - The adjoint values of the solution.

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -30,6 +30,7 @@
 
 #include "../../../Common/include/geometry/CGeometry.hpp"
 #include "../../include/solvers/CSolver.hpp"
+#include "../../include/variables/CIncEulerVariable.hpp"
 
 CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutput(config, nDim, false) {
 
@@ -326,6 +327,14 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("PRESSURE_COEFF", "Pressure_Coefficient", "PRIMITIVE", "Pressure coefficient");
   AddVolumeOutput("DENSITY",        "Density",              "PRIMITIVE", "Density");
 
+  // Density at previous timesteps for restart (unsteady with non-constant density)
+  if (config->GetTime_Domain() && config->GetKind_FluidModel() != CONSTANT_DENSITY) {
+    AddVolumeOutput("DENSITY_TIME_N", "Density_time_n", "SOLUTION", "Density at previous time step n");
+    if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+      AddVolumeOutput("DENSITY_TIME_N1", "Density_time_n1", "SOLUTION", "Density at previous time step n-1");
+    }
+  }
+
   if (config->GetKind_Solver() == MAIN_SOLVER::INC_RANS || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES){
     AddVolumeOutput("LAMINAR_VISCOSITY", "Laminar_Viscosity", "PRIMITIVE", "Laminar viscosity");
     AddVolumeOutput("HEAT_CAPACITY", "Heat_Capacity", "PRIMITIVE", "Heat capacity");
@@ -399,6 +408,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
 void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){
 
   const auto* Node_Flow = solver[FLOW_SOL]->GetNodes();
+  const auto* Node_Flow_Inc = static_cast<const CIncEulerVariable*>(Node_Flow);
   const CVariable* Node_Heat = nullptr;
   const CVariable* Node_Rad = nullptr;
   auto* Node_Geo = geometry->nodes;
@@ -437,6 +447,14 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
   SetVolumeOutputValue("PRESSURE_COEFF", iPoint, (Node_Flow->GetPressure(iPoint) - solver[FLOW_SOL]->GetPressure_Inf())/factor);
   SetVolumeOutputValue("DENSITY", iPoint, Node_Flow->GetDensity(iPoint));
 
+  // Load density at previous timesteps for restart (unsteady with non-constant density)
+  if (config->GetTime_Domain() && config->GetKind_FluidModel() != CONSTANT_DENSITY) {
+    SetVolumeOutputValue("DENSITY_TIME_N", iPoint, Node_Flow_Inc->GetDensity_time_n(iPoint));
+    if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) {
+      SetVolumeOutputValue("DENSITY_TIME_N1", iPoint, Node_Flow_Inc->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/CDiscAdjSolver.cpp

@@ -165,11 +165,18 @@ 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();
+    /*--- Density is stored separately from the solution vector for incompressible flows. ---*/
+    if (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE)
+      direct_solver->GetNodes()->RegisterDensity_time_n();
+  }
 
-  if (time_n1_needed)
+  if (time_n1_needed) {
     direct_solver->GetNodes()->RegisterSolution_time_n1();
+    if (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE)
+      direct_solver->GetNodes()->RegisterDensity_time_n1();
+  }
 }
 
 void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, bool reset) {
@@ -367,6 +374,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
     for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
       su2double Solution[MAXNVAR] = {0.0};
       direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution);
+      /*--- TODO: For incompressible flow, accumulate GetAdjointDensity_time_n
+       *    into the enthalpy adjoint via d(rho)/d(h) chain rule. ---*/
       nodes->Set_Solution_time_n(iPoint,Solution);
     }
     END_SU2_OMP_FOR
@@ -382,6 +391,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
     for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
       su2double Solution[MAXNVAR] = {0.0};
       direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution);
+      /*--- TODO: For incompressible flow, accumulate GetAdjointDensity_time_n1
+       *    into the enthalpy adjoint via d(rho)/d(h) chain rule. ---*/
       nodes->Set_Solution_time_n1(iPoint,Solution);
     }
     END_SU2_OMP_FOR