From 509a608293550637c2b038343865153b7a571314 Mon Sep 17 00:00:00 2001 From: Tahsin Berk Kiymaz Date: Tue, 17 Jun 2025 10:28:16 +0200 Subject: [PATCH 01/46] my changes --- SU2_CFD/include/output/COutput.hpp | 3 +- .../include/solvers/CFVMFlowSolverBase.inl | 7 ++++ SU2_CFD/include/solvers/CScalarSolver.inl | 8 ++-- SU2_CFD/include/variables/CEulerVariable.hpp | 9 ++++ .../include/variables/CIncEulerVariable.hpp | 18 +++++++- SU2_CFD/include/variables/CIncNSVariable.hpp | 2 +- .../variables/CSpeciesFlameletVariable.hpp | 2 +- .../include/variables/CSpeciesVariable.hpp | 2 +- SU2_CFD/include/variables/CVariable.hpp | 41 +++++++++++++++++++ SU2_CFD/src/integration/CIntegration.cpp | 1 + SU2_CFD/src/output/CFlowIncOutput.cpp | 6 +++ SU2_CFD/src/solvers/CIncEulerSolver.cpp | 25 +++++++++-- SU2_CFD/src/variables/CIncEulerVariable.cpp | 17 ++++++-- SU2_CFD/src/variables/CIncNSVariable.cpp | 7 +++- .../variables/CSpeciesFlameletVariable.cpp | 6 ++- SU2_CFD/src/variables/CSpeciesVariable.cpp | 4 +- SU2_CFD/src/variables/CVariable.cpp | 19 +++++++++ Tutorials | 1 + build_su2.slurm | 24 +++++++++++ build_su2_develop.slurm | 24 +++++++++++ 20 files changed, 205 insertions(+), 21 deletions(-) create mode 160000 Tutorials create mode 100644 build_su2.slurm create mode 100644 build_su2_develop.slurm diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index 14b0bbec9b6..90a3e4422fc 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -300,7 +300,8 @@ class COutput { unsigned short nRequestedVolumeFields; /*! \brief Minimum required volume fields for restart file. */ - const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY"}; + const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY","DENSITY_TIME_N"}; + /*----------------------------- Convergence monitoring ----------------------------*/ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 8a99712dc8e..121987ff038 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -882,6 +882,9 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** nodes->SetSolution(iPoint_Local, SolutionRestart); } + + // nodes->Set_Density_time_n(iPoint_Local, Restart_Data[index + nVar_Restart]); + /*--- For dynamic meshes, read in and store the grid coordinates and grid velocities for each node. ---*/ @@ -893,6 +896,7 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** /*--- Rewind the index to retrieve the Coords. ---*/ index = counter * Restart_Vars[1]; + const auto* Coord = &Restart_Data[index]; su2double GridVel[MAXNDIM] = {0.0}; @@ -914,6 +918,7 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** if (static_fsi && update_geo) { /*--- Rewind the index to retrieve the Coords. ---*/ index = counter*Restart_Vars[1]; + const auto* Coord = &Restart_Data[index]; for (auto iDim = 0u; iDim < nDim; iDim++) { @@ -1069,6 +1074,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_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 +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 1411de2ca78..33812ec035e 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -633,7 +633,7 @@ void CScalarSolver::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::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity(iPoint); // Temporary fix + Density_n = flowNodes->GetDensity_time_n(iPoint); // Temporary fix else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } @@ -750,7 +750,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); // Temporary fix else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } @@ -796,7 +796,7 @@ void CScalarSolver::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]; diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index c1dd2c45239..137f912fabc 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -184,6 +184,15 @@ 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 :) diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 7f398351f04..8f846ca0d8e 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/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,6 +99,20 @@ 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. @@ -147,6 +161,8 @@ class CIncEulerVariable : public CFlowVariable { */ inline su2double GetDensity(unsigned long iPoint) const final { return Primitive(iPoint, indices.Density()); } + + /*! * \brief Get the temperature of the flow. * \return Value of the temperature of the flow. diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp index 138026140b5..475b00e5412 100644 --- a/SU2_CFD/include/variables/CIncNSVariable.hpp +++ b/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); /*! diff --git a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp index a1eac80ae48..7d9cf63e850 100644 --- a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp @@ -48,7 +48,7 @@ class CSpeciesFlameletVariable final : public CSpeciesVariable { * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CSpeciesFlameletVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, + CSpeciesFlameletVariable(const su2double density, const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config); /*! diff --git a/SU2_CFD/include/variables/CSpeciesVariable.hpp b/SU2_CFD/include/variables/CSpeciesVariable.hpp index faeca8fa719..c6bc9d97612 100644 --- a/SU2_CFD/include/variables/CSpeciesVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesVariable.hpp @@ -48,7 +48,7 @@ class CSpeciesVariable : public CScalarVariable { * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CSpeciesVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, + CSpeciesVariable(su2double density,const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config); /*! diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 98e8c253819..1931dee89c1 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/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,12 @@ 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 +316,15 @@ 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 +519,23 @@ 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. @@ -1869,6 +1903,8 @@ class CVariable { */ inline su2double GetSolution_time_n1(unsigned long iPoint, unsigned long iVar) const { return Solution_time_n1(iPoint,iVar); } + + /*! * \brief Get the velocity (Structural Analysis). * \param[in] iVar - Index of the variable. @@ -2154,6 +2190,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. diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index a3af58d572b..8673ea4d223 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/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();) diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 7913d84f130..fe22d238660 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -293,6 +293,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("PRESSURE", "Pressure", "SOLUTION", "Pressure"); AddVolumeOutput("VELOCITY-X", "Velocity_x", "SOLUTION", "x-component of the velocity vector"); AddVolumeOutput("VELOCITY-Y", "Velocity_y", "SOLUTION", "y-component of the velocity vector"); + if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "SOLUTION", "z-component of the velocity vector"); if (heat || weakly_coupled_heat || flamelet) @@ -375,6 +376,9 @@ 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 +422,8 @@ 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)); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 6aaf09a5183..16a48ed770a 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/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(); @@ -1927,8 +1927,16 @@ void CIncEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, C void CIncEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) { CompleteImplicitIteration_impl(geometry, config); + + //if (config->GetTime_Domain()) { + // for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + // nodes->SetDensity_time_n(iPoint, nodes->GetDensity_unsteady(iPoint)); + // } + // } + } + void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { static su2double MaxVel2; @@ -2670,7 +2678,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,15 +2718,24 @@ 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); + if (iPoint==100) + cout <<"density="<GetKind_ConvNumScheme_Flow() == SPACE_CENTERED ? 2 : 4), config), @@ -46,7 +46,7 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci for(unsigned long iPoint=0; iPointGetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) { @@ -65,7 +66,14 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci } } -bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { +bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { + + static bool printed = false; + if (!printed) { + std::cout << "[DEBUG] SetPrimVar was called!" << std::endl; + printed = true; + } + bool physical = true; @@ -89,6 +97,8 @@ 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. ---*/ @@ -106,6 +116,7 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel FluidModel->SetTDState_T(Temperature); SetDensity(iPoint, FluidModel->GetDensity()); + /*--- Flag this point as non-physical. ---*/ physical = false; diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index cbcfa3b0aad..73412a341b8 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/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); @@ -52,6 +52,7 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do bool physical = true; + /*--- Set the value of the pressure ---*/ SetPressure(iPoint); @@ -78,6 +79,8 @@ 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. ---*/ diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp index 9b8d9b56212..960a3dbffee 100644 --- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -27,13 +27,14 @@ #include "../../include/variables/CSpeciesFlameletVariable.hpp" -CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, unsigned long npoint, +CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double density, const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config) - : CSpeciesVariable(species_inf, npoint, ndim, nvar, config) { + : CSpeciesVariable(density,species_inf, npoint, ndim, nvar, config) { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution(iPoint, iVar) = species_inf[iVar]; Solution_Old = Solution; + Density_unsteady = density; /*--- Allocate and initialize solution for the dual time strategy ---*/ bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || @@ -42,6 +43,7 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; + Density_unsteady = density; } /*--- Allocate residual structures ---*/ diff --git a/SU2_CFD/src/variables/CSpeciesVariable.cpp b/SU2_CFD/src/variables/CSpeciesVariable.cpp index 5715ae784ca..964595ba25f 100644 --- a/SU2_CFD/src/variables/CSpeciesVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesVariable.cpp @@ -27,7 +27,7 @@ #include "../../include/variables/CSpeciesVariable.hpp" -CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, +CSpeciesVariable::CSpeciesVariable(su2double density, const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config) : CScalarVariable(npoint, ndim, nvar, config) { /*--- Allocate space for the mass diffusivity. ---*/ @@ -38,6 +38,7 @@ CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long n Solution(iPoint, iVar) = species_inf[iVar]; Solution_Old = Solution; + Density_unsteady = density; /*--- Allocate and initialize solution for the dual time strategy ---*/ bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || @@ -46,5 +47,6 @@ CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long n if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; + Density_unsteady = density; } } diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 98ab867b2e9..3a36316106d 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -39,6 +39,11 @@ 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); @@ -56,13 +61,22 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva 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()) @@ -93,6 +107,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()); diff --git a/Tutorials b/Tutorials new file mode 160000 index 00000000000..e8d5cdbb960 --- /dev/null +++ b/Tutorials @@ -0,0 +1 @@ +Subproject commit e8d5cdbb96093e659f268b15cfe4a32572335256 diff --git a/build_su2.slurm b/build_su2.slurm new file mode 100644 index 00000000000..5185d731306 --- /dev/null +++ b/build_su2.slurm @@ -0,0 +1,24 @@ +#!/bin/bash + + +#SBATCH --job-name=build_su2 # Job name +#SBATCH -p genoa +#SBATCH --time=00:10:00 # Maximum runtime (adjust as needed) +#SBATCH --nodes=1 # Number of nodes +#SBATCH --ntasks=1 # Total tasks +#SBATCH --cpus-per-task=16 # CPU cores per task +#SBATCH --output=su2_build_%j.log # Output log file + +# Load required modules +module load 2022 +module load foss/2022a + +# Clone SU2 repository if not done already + + +/home/tberk/SU2/meson.py build -Dwarning_level=2 -Denable-autodiff=true -Denable-directdiff=false -Dwith-mpi=enabled -Denable-cgns=true -Denable-tecio=false -Denable-mlpcpp=true --prefix=/gpfs/home5/tberk/SU2 + +/home/tberk/SU2/ninja -C build install -j16 + + + diff --git a/build_su2_develop.slurm b/build_su2_develop.slurm new file mode 100644 index 00000000000..5185d731306 --- /dev/null +++ b/build_su2_develop.slurm @@ -0,0 +1,24 @@ +#!/bin/bash + + +#SBATCH --job-name=build_su2 # Job name +#SBATCH -p genoa +#SBATCH --time=00:10:00 # Maximum runtime (adjust as needed) +#SBATCH --nodes=1 # Number of nodes +#SBATCH --ntasks=1 # Total tasks +#SBATCH --cpus-per-task=16 # CPU cores per task +#SBATCH --output=su2_build_%j.log # Output log file + +# Load required modules +module load 2022 +module load foss/2022a + +# Clone SU2 repository if not done already + + +/home/tberk/SU2/meson.py build -Dwarning_level=2 -Denable-autodiff=true -Denable-directdiff=false -Dwith-mpi=enabled -Denable-cgns=true -Denable-tecio=false -Denable-mlpcpp=true --prefix=/gpfs/home5/tberk/SU2 + +/home/tberk/SU2/ninja -C build install -j16 + + + From 6aeb0a3598742447ed8bf49617a3358173c1fa9b Mon Sep 17 00:00:00 2001 From: Tahsin Berk Kiymaz Date: Tue, 17 Jun 2025 11:56:55 +0200 Subject: [PATCH 02/46] my changes --- Tutorials | 1 - 1 file changed, 1 deletion(-) delete mode 160000 Tutorials diff --git a/Tutorials b/Tutorials deleted file mode 160000 index e8d5cdbb960..00000000000 --- a/Tutorials +++ /dev/null @@ -1 +0,0 @@ -Subproject commit e8d5cdbb96093e659f268b15cfe4a32572335256 From 79e238ccf237f2e288387780279f90ddfc3e9be7 Mon Sep 17 00:00:00 2001 From: Tahsin Berk Kiymaz Date: Tue, 2 Sep 2025 00:12:36 +0200 Subject: [PATCH 03/46] specieschanges --- SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp | 2 +- SU2_CFD/include/variables/CSpeciesVariable.hpp | 2 +- SU2_CFD/src/numerics/flow/convection/fds.cpp | 6 +++--- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 6 +++--- SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp | 8 ++++---- SU2_CFD/src/variables/CSpeciesVariable.cpp | 6 +++--- 6 files changed, 15 insertions(+), 15 deletions(-) diff --git a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp index 7d9cf63e850..a1eac80ae48 100644 --- a/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesFlameletVariable.hpp @@ -48,7 +48,7 @@ class CSpeciesFlameletVariable final : public CSpeciesVariable { * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CSpeciesFlameletVariable(const su2double density, const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, + CSpeciesFlameletVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config); /*! diff --git a/SU2_CFD/include/variables/CSpeciesVariable.hpp b/SU2_CFD/include/variables/CSpeciesVariable.hpp index c6bc9d97612..faeca8fa719 100644 --- a/SU2_CFD/include/variables/CSpeciesVariable.hpp +++ b/SU2_CFD/include/variables/CSpeciesVariable.hpp @@ -48,7 +48,7 @@ class CSpeciesVariable : public CScalarVariable { * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CSpeciesVariable(su2double density,const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, + CSpeciesVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config); /*! diff --git a/SU2_CFD/src/numerics/flow/convection/fds.cpp b/SU2_CFD/src/numerics/flow/convection/fds.cpp index ae48e59f336..6a588ed9670 100644 --- a/SU2_CFD/src/numerics/flow/convection/fds.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fds.cpp @@ -157,9 +157,9 @@ CNumerics::ResidualType<> CUpwFDSInc_Flow::ComputeResidual(const CConfig *config MeandRhodT = 0.0; dRhodT_i = 0.0; dRhodT_j = 0.0; if (variable_density) { - MeandRhodT = -MeanDensity/MeanTemperature; - dRhodT_i = -DensityInc_i/Temperature_i; - dRhodT_j = -DensityInc_j/Temperature_j; + MeandRhodT = 0; + dRhodT_i = 0; + dRhodT_j = 0; } /*--- Compute ProjFlux_i ---*/ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 16a48ed770a..d32c475c522 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2008,7 +2008,7 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo law, but in the future, dRhodT should be in the fluid model. ---*/ if (variable_density) { - dRhodT = -Density/Temperature; + dRhodT = 0.0; } else { dRhodT = 0.0; } @@ -2732,8 +2732,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver V2U(Density_time_n, Cp, V_time_n, U_time_n); V2U(Density, Cp, V_time_nP1, U_time_nP1); - if (iPoint==100) - cout <<"density="<GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || @@ -43,7 +43,7 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double density, cons if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; - Density_unsteady = density; + } /*--- Allocate residual structures ---*/ diff --git a/SU2_CFD/src/variables/CSpeciesVariable.cpp b/SU2_CFD/src/variables/CSpeciesVariable.cpp index 964595ba25f..6f35e5bf9c3 100644 --- a/SU2_CFD/src/variables/CSpeciesVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesVariable.cpp @@ -27,7 +27,7 @@ #include "../../include/variables/CSpeciesVariable.hpp" -CSpeciesVariable::CSpeciesVariable(su2double density, const su2double* species_inf, unsigned long npoint, unsigned long ndim, +CSpeciesVariable::CSpeciesVariable( const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config) : CScalarVariable(npoint, ndim, nvar, config) { /*--- Allocate space for the mass diffusivity. ---*/ @@ -38,7 +38,7 @@ CSpeciesVariable::CSpeciesVariable(su2double density, const su2double* species_i Solution(iPoint, iVar) = species_inf[iVar]; Solution_Old = Solution; - Density_unsteady = density; + /*--- Allocate and initialize solution for the dual time strategy ---*/ bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || @@ -47,6 +47,6 @@ CSpeciesVariable::CSpeciesVariable(su2double density, const su2double* species_i if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; - Density_unsteady = density; + } } From 526b9a142892a8b4996651f328daa66685dd086b Mon Sep 17 00:00:00 2001 From: Tahsin Berk Kiymaz Date: Wed, 10 Dec 2025 14:56:00 +0100 Subject: [PATCH 04/46] Remove slurm files from repository --- build_su2.slurm | 24 ------------------------ build_su2_develop.slurm | 24 ------------------------ 2 files changed, 48 deletions(-) delete mode 100644 build_su2.slurm delete mode 100644 build_su2_develop.slurm diff --git a/build_su2.slurm b/build_su2.slurm deleted file mode 100644 index 5185d731306..00000000000 --- a/build_su2.slurm +++ /dev/null @@ -1,24 +0,0 @@ -#!/bin/bash - - -#SBATCH --job-name=build_su2 # Job name -#SBATCH -p genoa -#SBATCH --time=00:10:00 # Maximum runtime (adjust as needed) -#SBATCH --nodes=1 # Number of nodes -#SBATCH --ntasks=1 # Total tasks -#SBATCH --cpus-per-task=16 # CPU cores per task -#SBATCH --output=su2_build_%j.log # Output log file - -# Load required modules -module load 2022 -module load foss/2022a - -# Clone SU2 repository if not done already - - -/home/tberk/SU2/meson.py build -Dwarning_level=2 -Denable-autodiff=true -Denable-directdiff=false -Dwith-mpi=enabled -Denable-cgns=true -Denable-tecio=false -Denable-mlpcpp=true --prefix=/gpfs/home5/tberk/SU2 - -/home/tberk/SU2/ninja -C build install -j16 - - - diff --git a/build_su2_develop.slurm b/build_su2_develop.slurm deleted file mode 100644 index 5185d731306..00000000000 --- a/build_su2_develop.slurm +++ /dev/null @@ -1,24 +0,0 @@ -#!/bin/bash - - -#SBATCH --job-name=build_su2 # Job name -#SBATCH -p genoa -#SBATCH --time=00:10:00 # Maximum runtime (adjust as needed) -#SBATCH --nodes=1 # Number of nodes -#SBATCH --ntasks=1 # Total tasks -#SBATCH --cpus-per-task=16 # CPU cores per task -#SBATCH --output=su2_build_%j.log # Output log file - -# Load required modules -module load 2022 -module load foss/2022a - -# Clone SU2 repository if not done already - - -/home/tberk/SU2/meson.py build -Dwarning_level=2 -Denable-autodiff=true -Denable-directdiff=false -Dwith-mpi=enabled -Denable-cgns=true -Denable-tecio=false -Denable-mlpcpp=true --prefix=/gpfs/home5/tberk/SU2 - -/home/tberk/SU2/ninja -C build install -j16 - - - From 59979d0897198a54c96974e0a549b5e066250fbc Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:08:32 +0100 Subject: [PATCH 05/46] Remove unnecessary blank line in COutput.hpp --- SU2_CFD/include/output/COutput.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index 90a3e4422fc..67b30394eba 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -302,7 +302,6 @@ class COutput { /*! \brief Minimum required volume fields for restart file. */ const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY","DENSITY_TIME_N"}; - /*----------------------------- Convergence monitoring ----------------------------*/ su2double cauchyValue, /*!< \brief Summed value of the convergence indicator. */ From 99abfa429cf37603fd81b3a003b98657a1e44f18 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:13:22 +0100 Subject: [PATCH 06/46] Clean up commented-out density setting code Removed commented-out code for setting density at local point. --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 121987ff038..8ef8e422d18 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -882,9 +882,6 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** nodes->SetSolution(iPoint_Local, SolutionRestart); } - - // nodes->Set_Density_time_n(iPoint_Local, Restart_Data[index + nVar_Restart]); - /*--- For dynamic meshes, read in and store the grid coordinates and grid velocities for each node. ---*/ @@ -896,7 +893,6 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** /*--- Rewind the index to retrieve the Coords. ---*/ index = counter * Restart_Vars[1]; - const auto* Coord = &Restart_Data[index]; su2double GridVel[MAXNDIM] = {0.0}; @@ -918,7 +914,6 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** if (static_fsi && update_geo) { /*--- Rewind the index to retrieve the Coords. ---*/ index = counter*Restart_Vars[1]; - const auto* Coord = &Restart_Data[index]; for (auto iDim = 0u; iDim < nDim; iDim++) { From bee5a270010f2071d9d28db4437db65304cf7394 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:16:13 +0100 Subject: [PATCH 07/46] Update density retrieval comment for transient handling Updated the comment for density retrieval to reflect changes for transient density handling. --- SU2_CFD/include/solvers/CScalarSolver.inl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 33812ec035e..a023162463a 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -694,7 +694,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity_time_n(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::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity_time_n(iPoint); // Temporary fix + Density_n = flowNodes->GetDensity_time_n(iPoint); // Updated for transient density else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } From cea69cc22f6753e5029176a675cfe4a968467afd Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:17:10 +0100 Subject: [PATCH 08/46] Remove empty line before temperature setter comment --- SU2_CFD/include/variables/CEulerVariable.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 137f912fabc..44006d9d4fc 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -192,7 +192,6 @@ class CEulerVariable : public CFlowVariable { Density_unsteady[iPoint] = val; } - /*! * \brief Set the value of the temperature. * \param[in] temperature - how agitated the particles are :) From a721a2ec90b8df48778a067f3d87df0a3eb929f0 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:20:10 +0100 Subject: [PATCH 09/46] Clean up whitespace in CIncEulerVariable.hpp Removed unnecessary blank lines in CIncEulerVariable.hpp --- SU2_CFD/include/variables/CIncEulerVariable.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 8f846ca0d8e..7fe647ce44b 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -100,7 +100,6 @@ class CIncEulerVariable : public CFlowVariable { return val_density <= 0.0; } - inline void Set_Density_time_n(unsigned long iPoint, su2double val) { Density_time_n[iPoint] = val; } @@ -113,7 +112,6 @@ class CIncEulerVariable : public CFlowVariable { return Density_time_n[iPoint]; } - /*! * \brief Set the value of the density for the incompressible flows. * \param[in] iPoint - Point index. @@ -161,8 +159,6 @@ class CIncEulerVariable : public CFlowVariable { */ inline su2double GetDensity(unsigned long iPoint) const final { return Primitive(iPoint, indices.Density()); } - - /*! * \brief Get the temperature of the flow. * \return Value of the temperature of the flow. From 02c03bbfde11d7d6a0e4891246de92f456e9fedf Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:22:14 +0100 Subject: [PATCH 10/46] Clean up blank lines in CVariable.hpp Removed unnecessary blank lines in CVariable.hpp --- SU2_CFD/include/variables/CVariable.hpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 1931dee89c1..120d01f4cd2 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -282,7 +282,6 @@ class CVariable { */ void Set_Density_time_n(); - /*! * \brief Set the variable solution at time n. * \param[in] iPoint - Point index. @@ -324,7 +323,6 @@ class CVariable { Density_unsteady[iPoint] = val; } - /*! * \brief Virtual Member. Specify a vector to set the velocity components of the solution. * Multiplied by density for compressible cases. @@ -519,7 +517,6 @@ 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. @@ -1903,8 +1900,6 @@ class CVariable { */ inline su2double GetSolution_time_n1(unsigned long iPoint, unsigned long iVar) const { return Solution_time_n1(iPoint,iVar); } - - /*! * \brief Get the velocity (Structural Analysis). * \param[in] iVar - Index of the variable. From 0449bb47a3b730b398cb0179f41708f735f457ad Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:23:56 +0100 Subject: [PATCH 11/46] Clean up CFlowIncOutput.cpp by removing blank lines Removed unnecessary blank lines to improve code readability. --- SU2_CFD/src/output/CFlowIncOutput.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index fe22d238660..0476a1352fe 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -293,7 +293,6 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("PRESSURE", "Pressure", "SOLUTION", "Pressure"); AddVolumeOutput("VELOCITY-X", "Velocity_x", "SOLUTION", "x-component of the velocity vector"); AddVolumeOutput("VELOCITY-Y", "Velocity_y", "SOLUTION", "y-component of the velocity vector"); - if (nDim == 3) AddVolumeOutput("VELOCITY-Z", "Velocity_z", "SOLUTION", "z-component of the velocity vector"); if (heat || weakly_coupled_heat || flamelet) @@ -377,8 +376,6 @@ 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(); } @@ -424,7 +421,6 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve 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)); SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint)); From 572291f96e62e5669a927a699dd3ceef8068dc8c Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:28:50 +0100 Subject: [PATCH 12/46] Clean up CompleteImplicitIteration and update dRhodT Removed commented-out code and updated density calculation. --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index d32c475c522..26d93e5e942 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1927,16 +1927,8 @@ void CIncEulerSolver::PrepareImplicitIteration(CGeometry *geometry, CSolver**, C void CIncEulerSolver::CompleteImplicitIteration(CGeometry *geometry, CSolver**, CConfig *config) { CompleteImplicitIteration_impl(geometry, config); - - //if (config->GetTime_Domain()) { - // for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { - // nodes->SetDensity_time_n(iPoint, nodes->GetDensity_unsteady(iPoint)); - // } - // } - } - void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) { static su2double MaxVel2; @@ -2008,7 +2000,7 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo law, but in the future, dRhodT should be in the fluid model. ---*/ if (variable_density) { - dRhodT = 0.0; + dRhodT = -Density/Temperature; } else { dRhodT = 0.0; } @@ -2720,22 +2712,15 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver 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_time_n, Cp, V_time_n, U_time_n); V2U(Density, Cp, V_time_nP1, U_time_nP1); - /*if (iPoint==100) - cout <<"density="<GetDensity()); Density_unsteady[iPoint] = FluidModel->GetDensity(); - /*--- Non-physical solution found. Revert to old values. ---*/ if (check_dens || check_temp) { @@ -116,7 +113,6 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel FluidModel->SetTDState_T(Temperature); SetDensity(iPoint, FluidModel->GetDensity()); - /*--- Flag this point as non-physical. ---*/ physical = false; From 644a543f37bddb5885be5da80e08b2fe557879fb Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:37:50 +0100 Subject: [PATCH 15/46] Clean up blank lines in CIncNSVariable.cpp Removed unnecessary blank lines in SetPrimVar function. --- SU2_CFD/src/variables/CIncNSVariable.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 73412a341b8..904d4813335 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -52,7 +52,6 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do bool physical = true; - /*--- Set the value of the pressure ---*/ SetPressure(iPoint); @@ -81,7 +80,6 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity()); Density_unsteady[iPoint] = FluidModel->GetDensity(); - /*--- Non-physical solution found. Revert to old values. ---*/ if (check_dens || check_temp) { From fed4892e122b52e2f13c8c57ed50e4e0955ea099 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:40:28 +0100 Subject: [PATCH 16/46] Clean up whitespace in CVariable.cpp Removed unnecessary blank lines in CVariable.cpp. --- SU2_CFD/src/variables/CVariable.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 3a36316106d..04aa14bb34a 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -42,8 +42,6 @@ CVariable::CVariable(unsigned long npoint, unsigned long nvar, const CConfig *co Density_unsteady.resize(nPoint); Density_time_n.resize(nPoint); - - if (config->GetMultizone_Problem()) Solution_BGS_k.resize(nPoint,nVar) = su2double(0.0); @@ -60,10 +58,6 @@ 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()) @@ -77,7 +71,6 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva 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); From 8e09f7f8327ebf79450ac6ef1a34fb4a04d6da2e Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:43:10 +0100 Subject: [PATCH 17/46] Fix constructor signature for CSpeciesVariable --- SU2_CFD/src/variables/CSpeciesVariable.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_CFD/src/variables/CSpeciesVariable.cpp b/SU2_CFD/src/variables/CSpeciesVariable.cpp index 6f35e5bf9c3..a670d10f562 100644 --- a/SU2_CFD/src/variables/CSpeciesVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesVariable.cpp @@ -27,7 +27,7 @@ #include "../../include/variables/CSpeciesVariable.hpp" -CSpeciesVariable::CSpeciesVariable( const su2double* species_inf, unsigned long npoint, unsigned long ndim, +CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config) : CScalarVariable(npoint, ndim, nvar, config) { /*--- Allocate space for the mass diffusivity. ---*/ @@ -38,7 +38,6 @@ CSpeciesVariable::CSpeciesVariable( const su2double* species_inf, unsigned long Solution(iPoint, iVar) = species_inf[iVar]; Solution_Old = Solution; - /*--- Allocate and initialize solution for the dual time strategy ---*/ bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || From 1e5008a4ff3fdc2eaec8895c0b78ea60072961cd Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:43:53 +0100 Subject: [PATCH 18/46] Remove unnecessary blank line in CSpeciesVariable.cpp --- SU2_CFD/src/variables/CSpeciesVariable.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/src/variables/CSpeciesVariable.cpp b/SU2_CFD/src/variables/CSpeciesVariable.cpp index a670d10f562..5715ae784ca 100644 --- a/SU2_CFD/src/variables/CSpeciesVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesVariable.cpp @@ -46,6 +46,5 @@ CSpeciesVariable::CSpeciesVariable(const su2double* species_inf, unsigned long n if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; - } } From 2aebc5407f870cfda9c7c8ab6ce2201d45b8022c Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:45:11 +0100 Subject: [PATCH 19/46] Fix constructor signature for CSpeciesFlameletVariable --- SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp index 424cfc561f2..9b8d9b56212 100644 --- a/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp +++ b/SU2_CFD/src/variables/CSpeciesFlameletVariable.cpp @@ -27,14 +27,13 @@ #include "../../include/variables/CSpeciesFlameletVariable.hpp" -CSpeciesFlameletVariable::CSpeciesFlameletVariable( const su2double* species_inf, unsigned long npoint, +CSpeciesFlameletVariable::CSpeciesFlameletVariable(const su2double* species_inf, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig* config) : CSpeciesVariable(species_inf, npoint, ndim, nvar, config) { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution(iPoint, iVar) = species_inf[iVar]; Solution_Old = Solution; - /*--- Allocate and initialize solution for the dual time strategy ---*/ bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || @@ -43,7 +42,6 @@ CSpeciesFlameletVariable::CSpeciesFlameletVariable( const su2double* species_inf if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; - } /*--- Allocate residual structures ---*/ From 3f19c55fcdac688bd71329ebb7de1751b28b2cc5 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Wed, 10 Dec 2025 18:51:32 +0100 Subject: [PATCH 20/46] Remove debug output from SetPrimVar method Removed debug print statement from SetPrimVar function. --- SU2_CFD/src/variables/CIncEulerVariable.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 2d8fe1098b8..43c04053665 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -67,12 +67,6 @@ CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, cons } bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { - static bool printed = false; - if (!printed) { - std::cout << "[DEBUG] SetPrimVar was called!" << std::endl; - printed = true; - } - bool physical = true; /*--- Set the value of the pressure ---*/ From 735c27bd8f9126aba4576bb5e9eb6862b6211993 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:13:10 +0100 Subject: [PATCH 21/46] Update SU2_CFD/include/variables/CVariable.hpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/include/variables/CVariable.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 436125aba78..ce3d221db76 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -70,7 +70,7 @@ 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_unsteady; /*!< \brief density for unsteady flows. */ VectorType Density_time_n; /*!< \brief density at time n for dual-time stepping technique. */ CVectorOfMatrix Gradient; /*!< \brief Gradient of the solution of the problem. */ From f9f98a08223916f2f2cc033b9f6ddfca89277bde Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:13:23 +0100 Subject: [PATCH 22/46] Update SU2_CFD/include/output/COutput.hpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/include/output/COutput.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index be78782c0af..7a3e350d430 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/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 restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY","DENSITY_TIME_N"}; + const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY", "DENSITY_TIME_N"}; /*----------------------------- Convergence monitoring ----------------------------*/ From adf7581664c4a94f42fe7d01ce4cf79d8e049fc0 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:13:44 +0100 Subject: [PATCH 23/46] Update SU2_CFD/include/variables/CIncEulerVariable.hpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/include/variables/CIncEulerVariable.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index db07e055f5a..9c68044ebf7 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -105,7 +105,7 @@ class CIncEulerVariable : public CFlowVariable { Density_time_n[iPoint] = val; } - inline void Set_Density_unsteady(unsigned long iPoint, su2double val) { + inline void SetDensity_Unsteady(unsigned long iPoint, su2double val) { Density_unsteady[iPoint] = val; } From cc06d7e4037194e4a9bad61c415dd5466a46b41d Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:14:16 +0100 Subject: [PATCH 24/46] Update SU2_CFD/include/variables/CEulerVariable.hpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/include/variables/CEulerVariable.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 3053f1ff6c4..99a2237aa45 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -189,12 +189,12 @@ 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 SetDensity_Time_n(unsigned long iPoint, su2double val_density_time_n) { + Density_time_n[iPoint] = val_density_time_n; } - inline void Set_Density_unsteady(unsigned long iPoint, su2double val) { - Density_unsteady[iPoint] = val; + inline void SetDensity_Unsteady(unsigned long iPoint, su2double val_density_unsteady) { + Density_unsteady[iPoint] = val_density_unsteady; } /*! From 6f2a8eb32b04e4578e6aed21eac2a614af5bf8bb Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:16:13 +0100 Subject: [PATCH 25/46] Update SU2_CFD/src/solvers/CIncEulerSolver.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 737f1d737f9..10c9df3ac2e 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2754,7 +2754,6 @@ 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); From fa318d2d2d71352b9d92112cba0e71423870444e Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:16:29 +0100 Subject: [PATCH 26/46] Update SU2_CFD/src/variables/CVariable.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/variables/CVariable.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 63251cce304..64ae679a2ba 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -69,7 +69,7 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva } if (config->GetTime_Marching() != TIME_MARCHING::STEADY) - Solution_time_n1.resize(nPoint,nVar) = su2double(0.0); + Solution_time_n1.resize(nPoint,nVar) = su2double(0.0); /*--- User defined source terms ---*/ if (config->GetPyCustomSource()) UserDefinedSource.resize(nPoint,nVar) = su2double(0.0); From 3e1b898c328847a1dbfbb06018af04751652984d Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:17:07 +0100 Subject: [PATCH 27/46] Update SU2_CFD/src/variables/CVariable.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/variables/CVariable.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 64ae679a2ba..5ef080971da 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -60,10 +60,8 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva 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()) { + Solution_time_n.resize(nPoint,nVar) = su2double(0.0); Density_unsteady.resize(nPoint); Density_time_n.resize(nPoint); } From 87a33811e872b67684ea757c07906ff71dd46dfe Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:18:10 +0100 Subject: [PATCH 28/46] Update SU2_CFD/src/variables/CIncNSVariable.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/variables/CIncNSVariable.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index d2ecbfad083..660e9ed56ab 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -78,7 +78,6 @@ 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. ---*/ From b278ad68e5c948484eca06750da1374a339b1217 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:18:25 +0100 Subject: [PATCH 29/46] Update SU2_CFD/src/variables/CIncNSVariable.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/variables/CIncNSVariable.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 660e9ed56ab..58b7d8a6fc7 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -101,6 +101,9 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do } +/*--- Set density for unsteady problems ---*/ +Density_unsteady[iPoint] = FluidModel->GetDensity(); + /*--- Set the value of the velocity and velocity^2 (requires density) ---*/ SetVelocity(iPoint); From dda9271dba92887807150d2bf716db40465ca198 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:18:50 +0100 Subject: [PATCH 30/46] Update SU2_CFD/src/variables/CIncEulerVariable.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/variables/CIncEulerVariable.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 7bd87a051ad..ac31a8eaaa4 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -69,6 +69,7 @@ CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, cons } bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { + bool physical = true; /*--- Set the value of the pressure ---*/ From ff74a7ff35da1d27f5a86595cbeb0b98df3bf253 Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:19:06 +0100 Subject: [PATCH 31/46] Update SU2_CFD/include/solvers/CFVMFlowSolverBase.inl --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 1 + 1 file changed, 1 insertion(+) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 25f2a46d10f..65c6dd8b8a3 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1073,6 +1073,7 @@ void CFVMFlowSolverBase::PushSolutionBackInTime(unsigned long TimeIter, bo 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(); From c194c4f7d173fc9977f07110a8521facea9247af Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:19:36 +0100 Subject: [PATCH 32/46] Update SU2_CFD/src/variables/CIncEulerVariable.cpp --- SU2_CFD/src/variables/CIncEulerVariable.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index ac31a8eaaa4..cdbea02734e 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -58,7 +58,8 @@ CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, cons if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; - Density_unsteady = density; + Density_time_n = density; + Density_time_n1 = density; } if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) { From ac17146cd14775ee922c3865a5ff0fb93142aa0d Mon Sep 17 00:00:00 2001 From: Nijso Date: Tue, 16 Dec 2025 17:19:52 +0100 Subject: [PATCH 33/46] Update SU2_CFD/include/solvers/CFVMFlowSolverBase.inl --- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 1 + 1 file changed, 1 insertion(+) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 65c6dd8b8a3..2507188def3 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1105,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(); + solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Density_time_n1(); if (rans) solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n(); geometry[iMesh]->nodes->SetVolume_n(); From 35b018f7ef5c3594f63de73942a0236485972a70 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Tue, 16 Dec 2025 18:29:13 +0100 Subject: [PATCH 34/46] most of the second order stuff --- SU2_CFD/include/output/COutput.hpp | 2 +- .../include/solvers/CFVMFlowSolverBase.inl | 5 ++++ SU2_CFD/include/solvers/CScalarSolver.inl | 8 ++---- SU2_CFD/include/variables/CEulerVariable.hpp | 4 +++ .../include/variables/CIncEulerVariable.hpp | 10 ++++++- SU2_CFD/include/variables/CVariable.hpp | 26 +++++++++++++++++-- SU2_CFD/src/integration/CIntegration.cpp | 1 + SU2_CFD/src/output/CFlowIncOutput.cpp | 4 ++- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 6 +++-- SU2_CFD/src/variables/CVariable.cpp | 6 +++++ 10 files changed, 59 insertions(+), 13 deletions(-) diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index be78782c0af..ece3f3eb905 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/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 restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY","DENSITY_TIME_N"}; + const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY", "DENSITY_TIME_N", "DENSITY_TIME_N1"}; /*----------------------------- Convergence monitoring ----------------------------*/ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 25f2a46d10f..df66822c3c6 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1071,10 +1071,13 @@ 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(); + // 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()->Set_Density_time_n(); + 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(); } @@ -1102,8 +1105,10 @@ void CFVMFlowSolverBase::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()->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(); geometry[iMesh]->nodes->SetVolume_n(); diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 189024dfce8..49fc4ea927d 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -634,11 +634,7 @@ 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_nM1 = flowNodes->GetDensity_time_n1(iPoint); Density_n = flowNodes->GetDensity_time_n(iPoint); Density_nP1 = flowNodes->GetDensity(iPoint); } else { @@ -801,7 +797,7 @@ void CScalarSolver::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_nM1 = flowNodes->GetDensity_time_n1(iPoint); Density_n = flowNodes->GetDensity_time_n(iPoint); Density_nP1 = flowNodes->GetDensity(iPoint); } else { diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index 3053f1ff6c4..996134bb694 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -193,6 +193,10 @@ class CEulerVariable : public CFlowVariable { Density_time_n[iPoint] = val; } + inline void Set_Density_time_n1(unsigned long iPoint, su2double val) { + Density_time_n1[iPoint] = val; + } + inline void Set_Density_unsteady(unsigned long iPoint, su2double val) { Density_unsteady[iPoint] = val; } diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index db07e055f5a..6b425faeae0 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -105,6 +105,10 @@ class CIncEulerVariable : public CFlowVariable { Density_time_n[iPoint] = val; } + inline void Set_Density_time_n1(unsigned long iPoint, su2double val) { + Density_time_n1[iPoint] = val; + } + inline void Set_Density_unsteady(unsigned long iPoint, su2double val) { Density_unsteady[iPoint] = val; } @@ -112,7 +116,11 @@ class CIncEulerVariable : public CFlowVariable { 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. * \param[in] iPoint - Point index. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 436125aba78..d6485a98ce6 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -72,6 +72,7 @@ class CVariable { VectorType Delta_Time; /*!< \brief Time step. */ VectorType Density_unsteady; 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. */ @@ -298,6 +299,11 @@ class CVariable { */ void Set_Density_time_n(); + /*! + * \brief Set the density at time n-1. + */ + void Set_Density_time_n1(); + /*! * \brief Set the variable solution at time n. * \param[in] iPoint - Point index. @@ -334,7 +340,11 @@ class CVariable { inline void Set_Density_time_n(unsigned long iPoint, su2double val) { Density_time_n[iPoint] = val; } - + + inline void Set_Density_time_n1(unsigned long iPoint, su2double val) { + Density_time_n1[iPoint] = val; + } + inline void Set_Density_unsteady(unsigned long iPoint, su2double val) { Density_unsteady[iPoint] = val; } @@ -548,6 +558,14 @@ class CVariable { 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. @@ -2209,9 +2227,13 @@ class CVariable { void RegisterSolution_time_n1(); /*! - * \brief Register the variables in the solution_time_n1 array as input/output variable. + * \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. diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index fc1e269c726..199c377f7a0 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -249,6 +249,7 @@ void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver 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();) diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index d1d4257803a..53fb62b1478 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -375,7 +375,8 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddCommonFVMOutputs(config); - AddVolumeOutput("DENSITY_TIME_N", "Density_time_n", "SOLUTION", "Density at previous time step"); + 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(); } @@ -420,6 +421,7 @@ 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)); 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)); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 737f1d737f9..4e875f39e62 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2713,7 +2713,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, Density_time_n, Density_unsteady; + su2double Density, Cp, Density_time_n, Density_time_nM1, Density_unsteady; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST); @@ -2753,6 +2753,7 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Access the density and Cp at this node (constant for now). ---*/ + Density_time_nM1 = nodes->GetDensity_time_n1(iPoint); Density_time_n = nodes->GetDensity_time_n(iPoint); Density_unsteady = nodes->GetDensity_unsteady(iPoint); Density = nodes->GetDensity(iPoint); @@ -2760,8 +2761,9 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the conservative variable vector for all time levels. ---*/ - V2U(Density, Cp, V_time_nM1, U_time_nM1); + V2U(Density_time_nM1, Cp, V_time_nM1, U_time_nM1); V2U(Density_time_n, Cp, V_time_n, U_time_n); + // nijso asks: what is the difference between density_unsteady and density? V2U(Density, Cp, V_time_nP1, U_time_nP1); /*--- CV volume at time n+1. As we are on a static mesh, the volume diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 63251cce304..55ee4d007e7 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -41,6 +41,7 @@ CVariable::CVariable(unsigned long npoint, unsigned long nvar, const CConfig *co Density_unsteady.resize(nPoint); Density_time_n.resize(nPoint); + Density_time_n1.resize(nPoint); if (config->GetMultizone_Problem()) Solution_BGS_k.resize(nPoint,nVar) = su2double(0.0); @@ -66,6 +67,7 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva if (config->GetTime_Domain()) { Density_unsteady.resize(nPoint); Density_time_n.resize(nPoint); + Density_time_n1.resize(nPoint); } if (config->GetTime_Marching() != TIME_MARCHING::STEADY) @@ -108,6 +110,10 @@ void CVariable::Set_Density_time_n(){ parallelCopy(Density_unsteady.size(),Density_unsteady.data(), Density_time_n.data()); } +void CVariable::Set_Density_time_n1(){ + assert(Density_time_n.size() == Density_unsteady.size()); + parallelCopy(Density_time_n.size(),Density_time_n.data(), Density_time_n1.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()); From 630856b0477416e3c62c072c13fadbe2f8de36ac Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Thu, 18 Dec 2025 12:56:59 +0100 Subject: [PATCH 35/46] Update SU2_CFD/src/solvers/CIncEulerSolver.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index d62c91ae2f8..3d8ea5dec35 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2978,6 +2978,8 @@ 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. ---*/ From 35ad9763f190a4e0971bdcee321ae54c995f3467 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Thu, 18 Dec 2025 12:57:17 +0100 Subject: [PATCH 36/46] Update SU2_CFD/src/solvers/CIncEulerSolver.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 3d8ea5dec35..a5e69a3c39d 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2985,7 +2985,9 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver 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 From 31f8af3b107e3fd65b2b9c85ef80ccbe74ce786f Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Thu, 18 Dec 2025 12:57:57 +0100 Subject: [PATCH 37/46] Update SU2_CFD/src/solvers/CIncEulerSolver.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index a5e69a3c39d..5bf19323549 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2894,7 +2894,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); From 00ae736b85e0d3265727fbe74d1bc440ebf03a7a Mon Sep 17 00:00:00 2001 From: Nijso Date: Thu, 18 Dec 2025 16:23:23 +0100 Subject: [PATCH 38/46] Potential fix for code scanning alert no. 5888: Unused local variable Co-authored-by: Copilot Autofix powered by AI <62310815+github-advanced-security[bot]@users.noreply.github.com> --- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 5bf19323549..6c062b1397b 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2799,7 +2799,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, Density_time_n, Density_time_nM1, Density_unsteady; + 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); From df27a844a9f8a7f7f58be57388ae63e0e086b7af Mon Sep 17 00:00:00 2001 From: Nijso Date: Wed, 24 Dec 2025 18:04:41 +0100 Subject: [PATCH 39/46] Update SU2_CFD/src/variables/CIncNSVariable.cpp Co-authored-by: Cristopher Morales <98025159+Cristopher-Morales@users.noreply.github.com> --- SU2_CFD/src/variables/CIncNSVariable.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 0704e65d835..39c8ad55fa9 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -103,8 +103,9 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do } /*--- Set density for unsteady problems ---*/ -Density_unsteady[iPoint] = FluidModel->GetDensity(); - + if (Unsteady) { + SetDensity_unsteady(iPoint, FluidModel->GetDensity();) + } /*--- Set the value of the velocity and velocity^2 (requires density) ---*/ SetVelocity(iPoint); From adedf6bff7721a877b41790a2c006f451108bc00 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 16 Jan 2026 11:08:44 +0100 Subject: [PATCH 40/46] fix mpi and flame init --- SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 4247cea840a..79bc44949a1 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -88,7 +88,7 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver 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 < nPoint; i_point++) { + for (auto i_point = 0u; i_point < nPointDomain; i_point++) { CFluidModel* fluid_model_local = solver_container[FLOW_SOL]->GetFluidModel(); su2double* scalars = nodes->GetSolution(i_point); for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar]; @@ -154,14 +154,16 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** unsigned long ExtIter) { const bool restart = (config->GetRestart() || config->GetRestart_Flow()); - if ((!restart) && ExtIter == 0) { + bool flame_front_ignition = (flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::FLAME_FRONT); + + /*--- Also allow flame ignition when restarting. ---*/ + if (((!restart) && ExtIter == 0) || (restart && (flamelet_config_options.ignition_method != FLAMELET_INIT_TYPE::NONE))) { if (rank == MASTER_NODE) { cout << "Initializing progress variable and total enthalpy (using temperature)" << endl; } su2double flame_offset[3] = {0, 0, 0}, flame_normal[3] = {0, 0, 0}, flame_thickness = 0, flame_burnt_thickness = 0, flamenorm = 0; - bool flame_front_ignition = (flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::FLAME_FRONT); if (flame_front_ignition) { /*--- Collect flame front ignition parameters. ---*/ @@ -212,7 +214,6 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** for (unsigned long i_mesh = 0; i_mesh <= config->GetnMGLevels(); i_mesh++) { fluid_model_local = solver_container[i_mesh][FLOW_SOL]->GetFluidModel(); prog_burnt = GetBurntProgressVariable(fluid_model_local, scalar_init); - for (auto iVar = 0u; iVar < nVar; iVar++) scalar_init[iVar] = config->GetSpecies_Init()[iVar]; /*--- Set enthalpy based on initial temperature and scalars. ---*/ @@ -221,7 +222,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** prog_unburnt = config->GetSpecies_Init()[I_PROGVAR]; SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long i_point = 0; i_point < nPoint; i_point++) { + for (unsigned long i_point = 0; i_point < geometry[i_mesh]->GetnPointDomain(); i_point++) { auto coords = geometry[i_mesh]->nodes->GetCoord(i_point); if (flame_front_ignition) { @@ -270,6 +271,7 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** solver_container[i_mesh][SPECIES_SOL]->GetNodes()->SetSolution(i_point, scalar_init); } + END_SU2_OMP_FOR solver_container[i_mesh][SPECIES_SOL]->InitiateComms(geometry[i_mesh], config, MPI_QUANTITIES::SOLUTION); solver_container[i_mesh][SPECIES_SOL]->CompleteComms(geometry[i_mesh], config, MPI_QUANTITIES::SOLUTION); @@ -279,7 +281,6 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** solver_container[i_mesh][FLOW_SOL]->Preprocessing(geometry[i_mesh], solver_container[i_mesh], config, i_mesh, NO_RK_ITER, RUNTIME_FLOW_SYS, false); - END_SU2_OMP_FOR } /* --- Sum up some global counters over processes. --- */ From df281010aca8341258be3bf6ad18d9d1a7eafc04 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 16 Jan 2026 12:18:31 +0100 Subject: [PATCH 41/46] guards --- SU2_CFD/include/variables/CIncNSVariable.hpp | 1 + SU2_CFD/src/output/CFlowIncOutput.cpp | 14 ++++++++------ SU2_CFD/src/variables/CIncNSVariable.cpp | 15 ++++++++++++--- 3 files changed, 21 insertions(+), 9 deletions(-) diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp index a5bc25e5f6e..a76c768bfb6 100644 --- a/SU2_CFD/include/variables/CIncNSVariable.hpp +++ b/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: /*! diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index ee2e61c510e..96aef62b28f 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -390,9 +390,10 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ } AddCommonFVMOutputs(config); - - 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_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(); } @@ -438,9 +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)); - SetVolumeOutputValue("DENSITY_TIME_N", iPoint, Node_Flow->GetDensity_time_n(iPoint)); - SetVolumeOutputValue("DENSITY_TIME_N1", iPoint, Node_Flow->GetDensity_time_n1(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)); diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index ead0d59d8c5..53e83d43269 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -47,6 +47,15 @@ CIncNSVariable::CIncNSVariable(su2double density, su2double pressure, const su2d AuxVar.resize(nPoint,nAuxVar) = su2double(0.0); Grad_AuxVar.resize(nPoint,nAuxVar,nDim); } + + if (config->GetTime_Marching() != TIME_MARCHING::STEADY) { + Unsteady = true; + //Density_unsteady.resize(nPoint); + //Density_time_n.resize(nPoint); + //Density_time_n1.resize(nPoint); + } + + } bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel, const su2double *scalar) { @@ -94,9 +103,9 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do } /*--- Set density for unsteady problems ---*/ - if (Unsteady) { - SetDensity_unsteady(iPoint, FluidModel->GetDensity();) - } + if (Unsteady) + SetDensity_unsteady(iPoint, FluidModel->GetDensity()); + /*--- Set the value of the velocity and velocity^2 (requires density) ---*/ SetVelocity(iPoint); From 2308622e883a66d1ac3028559301e36283dc812b Mon Sep 17 00:00:00 2001 From: bigfooted Date: Fri, 16 Jan 2026 17:30:19 +0100 Subject: [PATCH 42/46] revert --- SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 79bc44949a1..d6e30dc4c24 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -88,7 +88,8 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver 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 < 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); for (auto iVar = 0u; iVar < nVar; iVar++) scalars_vector[iVar] = scalars[iVar]; @@ -222,7 +223,8 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** prog_unburnt = config->GetSpecies_Init()[I_PROGVAR]; SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long i_point = 0; i_point < geometry[i_mesh]->GetnPointDomain(); i_point++) { + //for (unsigned long i_point = 0; i_point < geometry[i_mesh]->GetnPointDomain(); i_point++) { + for (unsigned long i_point = 0; i_point < nPoint; i_point++) { auto coords = geometry[i_mesh]->nodes->GetCoord(i_point); if (flame_front_ignition) { From d529abeca74a7b5b955e77ffd5528069c3ca64c1 Mon Sep 17 00:00:00 2001 From: Nijso Date: Fri, 16 Jan 2026 19:03:17 +0100 Subject: [PATCH 43/46] Apply suggestion from @bigfooted --- SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index d6e30dc4c24..0a33623da76 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -223,7 +223,6 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** prog_unburnt = config->GetSpecies_Init()[I_PROGVAR]; SU2_OMP_FOR_STAT(omp_chunk_size) - //for (unsigned long i_point = 0; i_point < geometry[i_mesh]->GetnPointDomain(); i_point++) { for (unsigned long i_point = 0; i_point < nPoint; i_point++) { auto coords = geometry[i_mesh]->nodes->GetCoord(i_point); From 02c3311e84c24ab1e9587e3341b188cde7d69b57 Mon Sep 17 00:00:00 2001 From: Nijso Date: Fri, 16 Jan 2026 19:04:17 +0100 Subject: [PATCH 44/46] Apply suggestion from @bigfooted --- SU2_CFD/src/variables/CIncNSVariable.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 53e83d43269..0a74f95b6b9 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -50,9 +50,6 @@ CIncNSVariable::CIncNSVariable(su2double density, su2double pressure, const su2d if (config->GetTime_Marching() != TIME_MARCHING::STEADY) { Unsteady = true; - //Density_unsteady.resize(nPoint); - //Density_time_n.resize(nPoint); - //Density_time_n1.resize(nPoint); } From b355867cbfa3037f632fc5b8c6452239143113f0 Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 15 Mar 2026 21:59:03 +0100 Subject: [PATCH 45/46] update code --- Common/src/CConfig.cpp | 16 +++++ SU2_CFD/include/output/COutput.hpp | 2 +- .../include/solvers/CFVMFlowSolverBase.inl | 7 -- SU2_CFD/include/solvers/CScalarSolver.inl | 12 +--- SU2_CFD/include/variables/CEulerVariable.hpp | 12 ---- .../include/variables/CIncEulerVariable.hpp | 63 ++++++++++++------ SU2_CFD/include/variables/CIncNSVariable.hpp | 3 +- SU2_CFD/include/variables/CVariable.hpp | 62 +----------------- SU2_CFD/src/integration/CIntegration.cpp | 2 - SU2_CFD/src/output/CFlowIncOutput.cpp | 26 ++++++-- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 65 ++++++++++++------- .../src/solvers/CSpeciesFlameletSolver.cpp | 21 ++++-- SU2_CFD/src/variables/CIncEulerVariable.cpp | 41 ++++++++++-- SU2_CFD/src/variables/CIncNSVariable.cpp | 14 +--- SU2_CFD/src/variables/CVariable.cpp | 19 +----- 15 files changed, 177 insertions(+), 188 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 18884a52a35..d0edcc3de78 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5790,6 +5790,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()) { diff --git a/SU2_CFD/include/output/COutput.hpp b/SU2_CFD/include/output/COutput.hpp index 1c391a195e8..ce9e8737936 100644 --- a/SU2_CFD/include/output/COutput.hpp +++ b/SU2_CFD/include/output/COutput.hpp @@ -305,7 +305,7 @@ class COutput { unsigned short nRequestedVolumeFields; /*! \brief Minimum required volume fields for restart file. */ - const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY", "DENSITY_TIME_N", "DENSITY_TIME_N1"}; + const std::vector restartVolumeFields = {"COORDINATES", "SOLUTION", "SENSITIVITY", "GRID_VELOCITY"}; /*----------------------------- Convergence monitoring ----------------------------*/ diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 09d83fe5ffa..e60b1702e4b 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1072,13 +1072,9 @@ 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(); - // 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(); } @@ -1106,10 +1102,7 @@ void CFVMFlowSolverBase::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(); diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index e68e17c099f..c77aa8c3f8e 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -697,11 +697,9 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity_time_n(iPoint); // Updated for transient density - // nijso: no second order for dynamic meshes? + Density_n = flowNodes->GetDensity_time_n(iPoint); 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++) { @@ -755,7 +753,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity_time_n(iPoint); // Updated for transient density + Density_n = flowNodes->GetDensity_time_n(iPoint); else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } @@ -796,10 +794,6 @@ 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_time_n1(iPoint); Density_n = flowNodes->GetDensity_time_n(iPoint); Density_nP1 = flowNodes->GetDensity(iPoint); @@ -861,4 +855,4 @@ void CScalarSolver::PushSolutionBackInTime(unsigned long TimeIter, nodes->Set_Solution_time_n(); } } -} \ No newline at end of file +} diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index f117ca30774..e0fbcc5ce7b 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -189,18 +189,6 @@ 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 :) diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 15c859b5369..5c30e7b209a 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/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: /*! @@ -83,9 +85,20 @@ class CIncEulerVariable : public CFlowVariable { * \param[in] nvar - Number of variables of the problem. * \param[in] config - Definition of the particular problem. */ - CIncEulerVariable(su2double density, su2double pressure, const su2double *velocity, su2double enthalpy, + + 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. @@ -100,26 +113,6 @@ 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. @@ -310,4 +303,32 @@ 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; } + }; diff --git a/SU2_CFD/include/variables/CIncNSVariable.hpp b/SU2_CFD/include/variables/CIncNSVariable.hpp index a76c768bfb6..7dd082cdca2 100644 --- a/SU2_CFD/include/variables/CIncNSVariable.hpp +++ b/SU2_CFD/include/variables/CIncNSVariable.hpp @@ -41,7 +41,6 @@ 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: /*! @@ -54,7 +53,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 density, su2double pressure, const su2double *velocity, su2double temperature, + CIncNSVariable(su2double pressure, const su2double *velocity, su2double temperature, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config); /*! diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 4599dcc3461..63a05fe4116 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -70,9 +70,6 @@ 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. */ @@ -287,22 +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(); - - /*! - * \brief Set the density at time n. - */ - void SetDensity_time_n(); - - /*! - * \brief Set the density at time n-1. - */ - void SetDensity_time_n1(); + virtual void Set_Solution_time_n1(); /*! * \brief Set the variable solution at time n. @@ -337,18 +324,6 @@ 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. @@ -550,30 +525,6 @@ 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. @@ -2226,15 +2177,6 @@ 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. */ diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index e58f4868331..edb4248d8c1 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -248,8 +248,6 @@ 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();) diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 96aef62b28f..de894f70fe4 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/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"); @@ -390,10 +399,7 @@ 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(); } @@ -402,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(Node_Flow); const CVariable* Node_Heat = nullptr; const CVariable* Node_Rad = nullptr; auto* Node_Geo = geometry->nodes; @@ -439,10 +446,15 @@ 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)); + + // 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)); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 7bbb838563f..2b9b2898ad7 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/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(Density_Inf, Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); + nodes = new CIncNSVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); } else { - nodes = new CIncEulerVariable(Density_Inf, Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); + nodes = new CIncEulerVariable(Pressure_Inf, Velocity_Inf, Enthalpy_Inf, nPoint, nDim, nVar, config); } SetBaseClassPointerToNodes(); @@ -961,6 +961,10 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); const bool center_jst = (config->GetKind_Centered_Flow() == CENTERED::JST) && (iMesh == MESH_0); const bool outlet = (config->GetnMarker_Outlet() != 0); + const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); + const bool restart = config->GetRestart() || config->GetRestart_Flow(); + const auto TimeIter = config->GetTimeIter(); /*--- Set the primitive variables ---*/ @@ -969,6 +973,18 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ SU2_OMP_ATOMIC ErrorCounter += SetPrimitive_Variables(solver_container, config); + /*--- For the first time iteration without restart, update density time-levels after computing + actual density from SetPrimitive_Variables. ---*/ + if (dual_time && TimeIter == 0 && !restart && iRKStep == 0 && iMesh == MESH_0) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + su2double density = nodes->GetDensity(iPoint); + nodes->SetDensity_time_n(iPoint, density); + nodes->SetDensity_time_n1(iPoint, density); + } + END_SU2_OMP_FOR + } + if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { @@ -2804,7 +2820,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, Density_time_n, Density_time_nM1; + su2double Density; const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); const bool first_order = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST); @@ -2841,16 +2857,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 different time levels for non-constant density. ---*/ - Density_time_nM1 = nodes->GetDensity_time_n1(iPoint); - Density_time_n = nodes->GetDensity_time_n(iPoint); - Density = nodes->GetDensity(iPoint); + su2double Density_nM1 = nodes->GetDensity_time_n1(iPoint); + su2double Density_n = nodes->GetDensity_time_n(iPoint); + Density = nodes->GetDensity(iPoint); // Density at n+1 - /*--- Compute the conservative variable vector for all time levels. ---*/ + /*--- Compute the conservative variable vector for all time levels. + Use the density from the corresponding time level. ---*/ - V2U(Density_time_nM1, V_time_nM1, U_time_nM1); - V2U(Density_time_n, V_time_n, U_time_n); + V2U(Density_nM1, V_time_nM1, U_time_nM1); + V2U(Density_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 @@ -2898,9 +2915,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); -Density_time_n = nodes->GetDensity_time_n(iPoint); -V2U(Density_time_n, V_time_n, U_time_n); + su2double Density_n = nodes->GetDensity_time_n(iPoint); + V2U(Density_n, V_time_n, U_time_n); GridVel_i = geometry->nodes->GetGridVel(iPoint); @@ -2954,8 +2970,8 @@ V2U(Density_time_n, V_time_n, U_time_n); /*--- 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); + su2double 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; @@ -2981,19 +2997,18 @@ V2U(Density_time_n, V_time_n, U_time_n); 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 different time levels for non-constant density. ---*/ - Density = nodes->GetDensity(iPoint); - Density_time_nM1 = nodes->GetDensity_time_n1(iPoint); - Density_time_n = nodes->GetDensity_time_n(iPoint); + su2double Density_nM1 = nodes->GetDensity_time_n1(iPoint); + su2double Density_n = nodes->GetDensity_time_n(iPoint); + Density = nodes->GetDensity(iPoint); // Density at n+1 - /*--- Compute the conservative variable vector for all time levels. ---*/ + /*--- Compute the conservative variable vector for all time levels. + Use the density from the corresponding time level. ---*/ - 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); + V2U(Density_nM1, V_time_nM1, U_time_nM1); + V2U(Density_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 diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 990c0a9bc4a..2bb8b3f4a81 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -81,14 +81,18 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver auto spark_init = flamelet_config_options.spark_init; spark_iter_start = ceil(spark_init[4]); spark_duration = ceil(spark_init[5]); - unsigned long iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter(); + unsigned long iter; + if (config->GetTime_Domain()) { + iter = config->GetTimeIter(); // Use time step counter for unsteady problems + } else { + iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter(); + } ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration))); } 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); @@ -103,8 +107,12 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver spark_radius = flamelet_config_options.spark_init[3]; dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data()); if (dist_from_center < pow(spark_radius,2)) { - for (auto iVar = 0u; iVar < nVar; iVar++) - nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar]); + /*--- Add spark reaction rates to the sources that were just set by SetScalarSources ---*/ + const su2double* current_sources = nodes->GetScalarSources(i_point); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + nodes->SetScalarSource(i_point, iVar, current_sources[iVar] + flamelet_config_options.spark_reaction_rates[iVar]); + } } } @@ -195,10 +203,11 @@ void CSpeciesFlameletSolver::SetInitialCondition(CGeometry** geometry, CSolver** cout << "Ignition with a straight flame front" << endl; break; case FLAMELET_INIT_TYPE::SPARK: - cout << "Ignition with an artificial spark" << endl; + cout << "Ignition with an artificial spark at iteration "<< flamelet_config_options.spark_init[4] + << " for a duration of " << flamelet_config_options.spark_init[5] << " iterations." << endl; break; case FLAMELET_INIT_TYPE::NONE: - cout << "No solution ignition (cold flow)" << endl; + cout << "No solution ignition (cold flow or restart)" << endl; break; default: break; diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index cbe5e409b09..a6f004fd54b 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -27,8 +27,9 @@ #include "../../include/variables/CIncEulerVariable.hpp" #include "../../include/fluid/CFluidModel.hpp" +#include "../../../Common/include/parallelization/omp_structure.hpp" -CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, const su2double *velocity, su2double enthalpy, +CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) : CFlowVariable(npoint, ndim, nvar, ndim + 10, ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED ? 2 : 4), config), @@ -48,7 +49,7 @@ CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, cons for(unsigned long iPoint=0; iPointGetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) { @@ -69,7 +72,7 @@ CIncEulerVariable::CIncEulerVariable(su2double density, su2double pressure, cons } } -bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { +bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel) { bool physical = true; @@ -90,7 +93,6 @@ 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. ---*/ @@ -130,3 +132,30 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel return physical; } +void CIncEulerVariable::Set_Solution_time_n() { + /*--- Set the solution at time n ---*/ + CVariable::Set_Solution_time_n(); + + /*--- Store the current density at time n ---*/ + if (Density_time_n.size() > 0) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + Density_time_n(iPoint) = GetDensity(iPoint); + } + END_SU2_OMP_FOR + } +} + +void CIncEulerVariable::Set_Solution_time_n1() { + /*--- Set the solution at time n-1 ---*/ + CVariable::Set_Solution_time_n1(); + + /*--- Store the density at time n-1 by copying from time n ---*/ + if (Density_time_n1.size() > 0 && Density_time_n.size() > 0) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + Density_time_n1(iPoint) = Density_time_n(iPoint); + } + END_SU2_OMP_FOR + } +} diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp index 0a74f95b6b9..f5201c92cb5 100644 --- a/SU2_CFD/src/variables/CIncNSVariable.cpp +++ b/SU2_CFD/src/variables/CIncNSVariable.cpp @@ -28,9 +28,9 @@ #include "../../include/variables/CIncNSVariable.hpp" #include "../../include/fluid/CFluidModel.hpp" -CIncNSVariable::CIncNSVariable(su2double density, su2double pressure, const su2double *velocity, su2double enthalpy, +CIncNSVariable::CIncNSVariable(su2double pressure, const su2double *velocity, su2double enthalpy, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) : - CIncEulerVariable(density, pressure, velocity, enthalpy, npoint, ndim, nvar, config), + CIncEulerVariable(pressure, velocity, enthalpy, npoint, ndim, nvar, config), Energy(config->GetEnergy_Equation()) { Vorticity.resize(nPoint,3); @@ -47,12 +47,6 @@ CIncNSVariable::CIncNSVariable(su2double density, su2double pressure, const su2d AuxVar.resize(nPoint,nAuxVar) = su2double(0.0); Grad_AuxVar.resize(nPoint,nAuxVar,nDim); } - - if (config->GetTime_Marching() != TIME_MARCHING::STEADY) { - Unsteady = true; - } - - } bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2double turb_ke, CFluidModel *FluidModel, const su2double *scalar) { @@ -99,10 +93,6 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do } -/*--- Set density for unsteady problems ---*/ - if (Unsteady) - SetDensity_unsteady(iPoint, FluidModel->GetDensity()); - /*--- Set the value of the velocity and velocity^2 (requires density) ---*/ SetVelocity(iPoint); diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 24b6dc7325d..5cde8161b5e 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -39,10 +39,6 @@ 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); - Density_time_n1.resize(nPoint); - if (config->GetMultizone_Problem()) Solution_BGS_k.resize(nPoint,nVar) = su2double(0.0); @@ -61,12 +57,8 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva Solution.resize(nPoint,nVar) = su2double(0.0); Solution_Old.resize(nPoint,nVar) = su2double(0.0); - if (config->GetTime_Domain()) { + if (config->GetTime_Domain()) Solution_time_n.resize(nPoint,nVar) = su2double(0.0); - Density_unsteady.resize(nPoint); - Density_time_n.resize(nPoint); - Density_time_n1.resize(nPoint); - } if (config->GetTime_Marching() != TIME_MARCHING::STEADY) Solution_time_n1.resize(nPoint,nVar) = su2double(0.0); @@ -103,15 +95,6 @@ void CVariable::Set_Solution_time_n() { parallelCopy(Solution.size(), Solution.data(), Solution_time_n.data()); } -void CVariable::SetDensity_time_n(){ - assert(Density_time_n.size() == Density_unsteady.size()); - parallelCopy(Density_unsteady.size(),Density_unsteady.data(), Density_time_n.data()); -} - -void CVariable::SetDensity_time_n1(){ - assert(Density_time_n1.size() == Density_time_n.size()); - parallelCopy(Density_time_n.size(),Density_time_n.data(), Density_time_n1.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()); From 76d652ce73d2e5ecce6013ecd50d7aa8e57e8e9b Mon Sep 17 00:00:00 2001 From: bigfooted Date: Sun, 15 Mar 2026 22:38:32 +0100 Subject: [PATCH 46/46] fix density registration for AD --- .../include/variables/CIncEulerVariable.hpp | 36 ++++++++++++++++ SU2_CFD/include/variables/CVariable.hpp | 42 +++++++++++++++++++ SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 15 ++++++- SU2_CFD/src/variables/CIncEulerVariable.cpp | 10 ++++- 4 files changed, 100 insertions(+), 3 deletions(-) diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 5c30e7b209a..7d43605bb36 100644 --- a/SU2_CFD/include/variables/CIncEulerVariable.hpp +++ b/SU2_CFD/include/variables/CIncEulerVariable.hpp @@ -331,4 +331,40 @@ class CIncEulerVariable : public CFlowVariable { */ 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); + } + }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 63a05fe4116..9c931d430c8 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -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. diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 5f6d8beecac..d380d8e6b79 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/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 diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index a6f004fd54b..a8809504afe 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -150,7 +150,7 @@ void CIncEulerVariable::Set_Solution_time_n1() { /*--- Set the solution at time n-1 ---*/ CVariable::Set_Solution_time_n1(); - /*--- Store the density at time n-1 by copying from time n ---*/ + /*--- Store the density at n-1 by copying from n ---*/ if (Density_time_n1.size() > 0 && Density_time_n.size() > 0) { SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { @@ -159,3 +159,11 @@ void CIncEulerVariable::Set_Solution_time_n1() { END_SU2_OMP_FOR } } + +void CIncEulerVariable::RegisterDensity_time_n() { + RegisterContainer(true, Density_time_n); +} + +void CIncEulerVariable::RegisterDensity_time_n1() { + RegisterContainer(true, Density_time_n1); +}