update density implementation for flow solver

Author: bigfooted

Common/include/option_structure.hpp

@@ -2565,6 +2565,7 @@ enum class MPI_QUANTITIES {
   SOLUTION_FEA         ,  /*!< \brief FEA solution communication. */
   MESH_DISPLACEMENTS   ,  /*!< \brief Mesh displacements at the interface. */
   SOLUTION_TIME_N      ,  /*!< \brief Solution at time n. */
+ // DENSITY_TIME_N      ,  /*!< \brief Density at time n. */
   SOLUTION_TIME_N1     ,  /*!< \brief Solution at time n-1. */
 };
 

SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp

@@ -991,6 +991,9 @@ class CFVMFlowSolverBase : public CSolver {
         for (unsigned short iVar = 0; iVar < nVar; iVar++) {
           nodes->AddSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint)*LinSysSol[iPoint*nVar+iVar]);
         }
+        // nijso: here we can update the density (for incompressible euler including energy equation)
+        // but: we do not need the intermediate density fields actually.
+        nodes->AddDensity(iPoint,0,nodes->GetDensity(iPoint));
       }
       END_SU2_OMP_FOR
     }

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -613,6 +613,8 @@ void CFVMFlowSolverBase<V, R>::ImplicitEuler_Iteration(CGeometry *geometry, CSol
   END_SU2_OMP_SAFE_GLOBAL_ACCESS
 
   CompleteImplicitIteration(geometry, nullptr, config);
+
+  // can we update density here?
 }
 
 template <class V, ENUM_REGIME R>
@@ -1065,10 +1067,14 @@ void CFVMFlowSolverBase<V, R>::PushSolutionBackInTime(unsigned long TimeIter, bo
                                                       CConfig* config) {
   /*--- Push back the initial condition to previous solution containers
    for a 1st-order restart or when simply initializing to freestream. ---*/
-
+  cout << "restart for unsteady solution" << endl;
   for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
     solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n();
+    // nijso:note only for second order!
     solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n1();
+    //if ((config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) && (config->GetEnergy_Equation()==true))
+      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();
@@ -1574,7 +1580,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::SetResidual_DualTime(CGeometry *geometry
   TimeStep = config->GetDelta_UnstTimeND();
 
   /*--- Compute the dual time-stepping source term for static meshes ---*/
-
+  cout << "setresidual_dualtime cfvmflowsolverbase" << endl;
   if (!dynamic_grid) {
 
     /*--- Loop over all nodes (excluding halos) ---*/

SU2_CFD/include/variables/CIncEulerVariable.hpp

@@ -71,8 +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 timestep n for time-dependent flow. */
-  VectorType Cp_time_n;    /*!< \brief Cp at timestep n for time-dependent flow. */
+  //VectorType Density_time_n;    /*!< \brief density at timestep n for time-dependent flow. */
+  //VectorType Cp_time_n;    /*!< \brief Cp at timestep n for time-dependent flow. */
 
  public:
   /*!
@@ -85,7 +85,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);
 
   /*!

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

@@ -53,6 +53,8 @@ class CVariable {
   using MatrixType = C2DContainer<unsigned long, su2double, StorageType::RowMajor,    64, DynamicSize, DynamicSize>;
 
   MatrixType Solution;       /*!< \brief Solution of the problem. */
+  // necessary?
+  VectorType Density;       /*!< \brief Solution of the problem. */
   MatrixType Solution_Old;   /*!< \brief Old solution of the problem R-K. */
 
   MatrixType External;       /*!< \brief External (outer) contribution in discrete adjoint multizone problems. */
@@ -67,7 +69,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 Density_time_n;     /*!< \brief Solution of the problem at time n-1 for dual-time stepping technique. */
   VectorType Delta_Time;         /*!< \brief Time step. */
+  VectorType Cp_time_n;     /*!< \brief Solution of the problem at time n-1 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. */
@@ -270,6 +274,11 @@ class CVariable {
    */
   void Set_Solution_time_n();
 
+  /*!
+   * \brief Set the variable solution at time n.
+   */
+  void Set_Density_time_n();
+
   /*!
    * \brief Set the variable solution at time n-1.
    */
@@ -280,6 +289,7 @@ class CVariable {
    * \param[in] iPoint - Point index.
    */
   inline void Set_Solution_time_n(unsigned long iPoint, const su2double* val_sol) {
+    cout <<"cvariable: set_solution_time_n" << endl;
     for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_time_n(iPoint,iVar) = val_sol[iVar];
   }
 
@@ -297,8 +307,27 @@ class CVariable {
    * \param[in] iPoint - Point index.
    */
   inline void Set_Solution_time_n(unsigned long iPoint, unsigned long iVar, su2double val_sol) {
+    cout <<"cvariable: set_solution_time_n 111" << endl;
     Solution_time_n(iPoint,iVar) = val_sol;
   }
+  /*!
+   * \brief Set the variable solution at time n.
+   * \param[in] iPoint - Point index.
+   */
+
+  inline void Set_Density_time_n(unsigned long iPoint, su2double val_sol) {
+    Density_time_n(iPoint) = val_sol;
+  }
+
+  /*!
+   * \brief Add density.
+   * \param[in] iPoint - Point index.
+   * \param[in] iVar - Number of the variable.
+   * \param[in] solution - Value that we want to add to the solution.
+   */
+  inline void AddDensity(unsigned long iPoint, unsigned long iVar, su2double solution) {
+    Density(iPoint) = solution;
+  }
 
   /*!
    * \brief Set the variable solution at time n-1.
@@ -462,14 +491,18 @@ class CVariable {
    * \return Reference to the solution matrix.
    */
   inline const MatrixType& GetSolution() const { return Solution; }
-  inline MatrixType& GetSolution() { return Solution; }
+  inline MatrixType& GetSolution() { 
+    return Solution; 
+    }
 
   /*!
    * \brief Get the solution of the problem.
    * \param[in] iPoint - Point index.
    * \return Pointer to the solution vector.
    */
-  inline su2double *GetSolution(unsigned long iPoint) { return Solution[iPoint]; }
+  inline su2double *GetSolution(unsigned long iPoint) { 
+    return Solution[iPoint]; 
+    }
 
   /*!
    * \brief Get the old solution of the problem (Runge-Kutta method)
@@ -491,8 +524,23 @@ class CVariable {
    * \param[in] iPoint - Point index.
    * \return Pointer to the solution (at time n) vector.
    */
-  inline su2double *GetSolution_time_n(unsigned long iPoint) { return Solution_time_n[iPoint]; }
-  inline MatrixType& GetSolution_time_n() { return Solution_time_n; }
+  inline su2double *GetSolution_time_n(unsigned long iPoint) { 
+    return Solution_time_n[iPoint]; }
+  inline MatrixType& GetSolution_time_n() { 
+    cout << "getsolution"<<endl;
+    return Solution_time_n; }
+
+  /*!
+   * \brief Get the solution at time n.
+   * \param[in] iPoint - Point index.
+   * \return Pointer to the solution (at time n) vector.
+   */
+  // we do not use this, we use the one in cinceulervariable
+  //inline su2double GetDensity_time_n(unsigned long iPoint) { 
+  //  return Density_time_n[iPoint]; }
+  //inline VectorType& GetDensity_time_n() { 
+  //  cout << "getsolution"<<endl;
+  //  return Density_time_n; }
 
   /*!
    * \brief Get the solution at time n-1.
@@ -2154,6 +2202,17 @@ class CVariable {
    */
   void RegisterSolution_time_n1();
 
+ /*!
+   * \brief Register the Density_time_n array as input/output variable.
+   */
+  //void RegisterDensity_time_n();
+
+
+  /*!
+   * \brief Register the variables in the solution_time_n1 array as input/output variable.
+   */
+  //void RegisterDensity_time_n1();
+
   /*!
    * \brief Set the adjoint values of the solution.
    * \param[in] adj_sol - The adjoint values of the solution.

SU2_CFD/src/integration/CIntegration.cpp

@@ -246,9 +246,13 @@ void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver
   SU2_OMP_PARALLEL
   {
   /*--- Store old solution ---*/
+  // nijso: note that we only need time_n1 for second order timestepping, we can save some memory here.
   solver->GetNodes()->Set_Solution_time_n1();
   solver->GetNodes()->Set_Solution_time_n();
 
+  solver->GetNodes()->Set_Density_time_n();
+  //solver->GetNodes()->Set_Density_time_n1();
+
   SU2_OMP_SAFE_GLOBAL_ACCESS(solver->ResetCFLAdapt();)
 
   SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads()))

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -167,6 +167,10 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) {
   if (time_n_needed)
     direct_solver->GetNodes()->RegisterSolution_time_n();
 
+// nijso: we probably need this as well for adjoint unsteady incompressible with varying density...
+//if (time_n_needed)
+//    direct_solver->GetNodes()->RegisterDensity_time_n();
+
   if (time_n1_needed)
     direct_solver->GetNodes()->RegisterSolution_time_n1();
 }

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();
 
@@ -2717,16 +2717,18 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
       //Cp_time_n = nodes->GetSpecificHeatCp(iPoint);
 
       Density_time_n = nodes->GetDensity_time_n(iPoint);
-      Cp_time_n = nodes->GetCp_time_n(iPoint);
+      //Cp_time_n = nodes->GetCp_time_n(iPoint);
 
 
 
       /*--- Compute the conservative variable vector for all time levels. ---*/
-
+      
       V2U(Density, Cp, V_time_nM1, U_time_nM1);
-      V2U(Density_time_n, Cp_time_n, V_time_n, U_time_n);
+      //V2U(Density_time_n, Cp_time_n, 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);
-
+      if (iPoint == 10065)
+        cout << iPoint << ", " << Density << ", " << Density_time_n << ", " << V_time_nP1[0] << ", " << V_time_n[0] << endl;
       /*--- 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. ---*/
 
@@ -2909,15 +2911,15 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver
 
 
   // update the density and cp
-  SU2_OMP_FOR_STAT(omp_chunk_size)
-  for (iPoint = 0; iPoint < nPointDomain; iPoint++) {      // set density of time n to density of time n+1
-    Density = nodes->GetDensity(iPoint);
-    Cp = nodes->GetSpecificHeatCp(iPoint);
-
-    nodes->SetDensity_time_n(iPoint, Density);
-    nodes->SetCp_time_n(iPoint, Cp);
-  }
-  END_SU2_OMP_FOR
+  //SU2_OMP_FOR_STAT(omp_chunk_size)
+  //for (iPoint = 0; iPoint < nPointDomain; iPoint++) {      // set density of time n to density of time n+1
+  //  Density = nodes->GetDensity(iPoint);
+  //  //Cp = nodes->GetSpecificHeatCp(iPoint);
+
+  //  nodes->SetDensity_time_n(iPoint, Density);
+    //nodes->SetCp_time_n(iPoint, Cp);
+  //}
+  //END_SU2_OMP_FOR
 
 
 }

SU2_CFD/src/solvers/CMeshSolver.cpp

@@ -881,6 +881,7 @@ void CMeshSolver::RestartOldGeometry(CGeometry *geometry, const CConfig *config)
   for(unsigned short iStep = 1; iStep <= nSteps; ++iStep) {
 
     MPI_QUANTITIES CommType = (iStep == 1) ? MPI_QUANTITIES::SOLUTION_TIME_N : MPI_QUANTITIES::SOLUTION_TIME_N1;
+    //MPI_QUANTITIES CommType = (iStep == 1) ? MPI_QUANTITIES::DENSITY_TIME_N : MPI_QUANTITIES::DENSITY_TIME_N1;
 
     /*--- Modify file name for an unsteady restart ---*/
     int Unst_RestartIter;