diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index e60b1702e4b..c66dddbb1e1 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1073,6 +1073,8 @@ void CFVMFlowSolverBase::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(); + solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Density_time_n1(); if (rans) { solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n(); solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n1(); @@ -1103,6 +1105,7 @@ void CFVMFlowSolverBase::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(); diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 50927ccbcc0..7a4e1f1fc0b 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -635,12 +635,9 @@ void CScalarSolver::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); + /*--- Use density at each time level for variable-density flows. ---*/ + 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 +698,7 @@ void CScalarSolver::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,7 +754,7 @@ void CScalarSolver::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); } @@ -798,12 +795,9 @@ void CScalarSolver::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); + /*--- Use density at each time level for variable-density flows. ---*/ + 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]; diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 9cc13f79986..4b3261a91c8 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -72,6 +72,10 @@ 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. */ su2double TemperatureLimits[2]; /*!< \brief Temperature limits [K]. */ + + VectorType Density_time_n; /*!< \brief Density at time level n for dual-time stepping (variable-density flows). */ + VectorType Density_time_n1; /*!< \brief Density at time level n-1 for dual-time stepping (variable-density flows). */ + public: /*! * \brief Constructor of the class. @@ -282,6 +286,34 @@ class CIncEulerVariable : public CFlowVariable { return Streamwise_Periodic_RecoveredTemperature(iPoint); } + /*! + * \brief Get the density at time level n. + * \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(iPoint); + } + + /*! + * \brief Get the density at time level n-1. + * \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(iPoint); + } + + /*! + * \brief Store current density into time level n: Density_time_n = Density. + */ + void Set_Density_time_n() final; + + /*! + * \brief Store time level n density into time level n-1: Density_time_n1 = Density_time_n. + */ + void Set_Density_time_n1() final; + /*! * \brief Specify a vector to set the velocity components of the solution. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index bc4bc55e57c..718b2805fc2 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -970,6 +970,30 @@ class CVariable { */ inline virtual su2double GetDensity(unsigned long iPoint, unsigned long val_iSpecies) const { return 0.0; } + /*! + * \brief Get the density at time level n (for incompressible unsteady flows with variable density). + * \param[in] iPoint - Point index. + * \return Density at time level n. Defaults to current density. + */ + inline virtual su2double GetDensity_time_n(unsigned long iPoint) const { return GetDensity(iPoint); } + + /*! + * \brief Get the density at time level n-1 (for incompressible unsteady flows with variable density). + * \param[in] iPoint - Point index. + * \return Density at time level n-1. Defaults to current density. + */ + inline virtual su2double GetDensity_time_n1(unsigned long iPoint) const { return GetDensity(iPoint); } + + /*! + * \brief Store current density into the time level n container: Density_time_n = Density. + */ + inline virtual void Set_Density_time_n() {} + + /*! + * \brief Store time level n density into the time level n-1 container: Density_time_n1 = Density_time_n. + */ + inline virtual void Set_Density_time_n1() {} + /*! * \brief A virtual member. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index edb4248d8c1..2ff0154c795 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -249,6 +249,10 @@ void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver solver->GetNodes()->Set_Solution_time_n1(); solver->GetNodes()->Set_Solution_time_n(); + /*--- Store density at previous time levels (for incompressible variable-density flows). ---*/ + solver->GetNodes()->Set_Density_time_n1(); + solver->GetNodes()->Set_Density_time_n(); + SU2_OMP_SAFE_GLOBAL_ACCESS(solver->ResetCFLAdapt();) SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index c99f87896c2..006481557c1 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -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_nP1, Density_n, Density_nM1; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST); @@ -2841,15 +2841,17 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = nodes->GetSolution_time_n(iPoint); V_time_nP1 = nodes->GetSolution(iPoint); - /*--- Access the density at this node (constant for now). ---*/ + /*--- Access the density at this node for each time level. ---*/ - Density = nodes->GetDensity(iPoint); + Density_nP1 = nodes->GetDensity(iPoint); + Density_n = nodes->GetDensity_time_n(iPoint); + Density_nM1 = nodes->GetDensity_time_n1(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_nM1, V_time_nM1, U_time_nM1); + V2U(Density_n, V_time_n, U_time_n); + V2U(Density_nP1, V_time_nP1, U_time_nP1); /*--- 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. ---*/ @@ -2870,7 +2872,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ if (implicit) { - su2double delta = (second_order ? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; + su2double delta = (second_order ? 1.5 : 1.0) * Volume_nP1 * Density_nP1 / TimeStep; for (iVar = 1; iVar < nVar; ++iVar) Jacobian.AddVal2Diag(iPoint, iVar, delta); } @@ -2896,8 +2898,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the conservative variables. ---*/ V_time_n = nodes->GetSolution_time_n(iPoint); - Density = nodes->GetDensity(iPoint); - V2U(Density, V_time_n, U_time_n); + Density_n = nodes->GetDensity_time_n(iPoint); + V2U(Density_n, V_time_n, U_time_n); GridVel_i = geometry->nodes->GetGridVel(iPoint); @@ -2951,8 +2953,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the GCL component of the source term for node i ---*/ V_time_n = nodes->GetSolution_time_n(iPoint); - Density = nodes->GetDensity(iPoint); - V2U(Density, V_time_n, U_time_n); + Density_n = nodes->GetDensity_time_n(iPoint); + V2U(Density_n, V_time_n, U_time_n); for (iVar = 0; iVar < nVar-!energy; iVar++) LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; @@ -2978,15 +2980,17 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V_time_n = nodes->GetSolution_time_n(iPoint); V_time_nP1 = nodes->GetSolution(iPoint); - /*--- Access the density at this node (constant for now). ---*/ + /*--- Access the density at this node for each time level. ---*/ - Density = nodes->GetDensity(iPoint); + Density_nP1 = nodes->GetDensity(iPoint); + Density_n = nodes->GetDensity_time_n(iPoint); + Density_nM1 = nodes->GetDensity_time_n1(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_nM1, V_time_nM1, U_time_nM1); + V2U(Density_n, V_time_n, U_time_n); + V2U(Density_nP1, 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 @@ -3010,7 +3014,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ if (implicit) { - su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; + su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density_nP1 / TimeStep; for (iVar = 1; iVar < nVar; ++iVar) Jacobian.AddVal2Diag(iPoint, iVar, delta); diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 6d77af4454d..0d2e715c347 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -58,6 +58,12 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; + + /*--- Allocate density storage for previous time levels. + * This is needed for variable-density incompressible unsteady flows + * (e.g. with energy equation or species transport). ---*/ + Density_time_n.resize(nPoint) = su2double(0.0); + Density_time_n1.resize(nPoint) = su2double(0.0); } if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) { @@ -122,8 +128,22 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel /*--- Set enthalpy ---*/ - SetEnthalpy(iPoint, FluidModel->GetEnthalpy()); + SetEnthalpy(iPoint, FluidModel->GetEnthalpy()); return physical; } + +void CIncEulerVariable::Set_Density_time_n() { + const auto nPoint = Density_time_n.size(); + SU2_OMP_FOR_STAT(roundUpDiv(nPoint, omp_get_num_threads())) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + Density_time_n(iPoint) = GetDensity(iPoint); + } + END_SU2_OMP_FOR +} + +void CIncEulerVariable::Set_Density_time_n1() { + assert(Density_time_n1.size() == Density_time_n.size()); + parallelCopy(Density_time_n.size(), Density_time_n.data(), Density_time_n1.data()); +}