diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 02c94dd1198..5fd02fbd669 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -233,7 +233,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -282,7 +282,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} -a "--tapetests" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:260405-0054 with: @@ -330,7 +330,7 @@ jobs: PMIX_MCA_gds: hash with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:260405-0054 with: @@ -375,7 +375,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c fgm_unsteady -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:260405-0054 with: diff --git a/SU2_CFD/include/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index 7d6871f8ca4..dc005809634 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -28,6 +28,7 @@ #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../include/solvers/CScalarSolver.hpp" #include "../../include/variables/CFlowVariable.hpp" +#include "../../include/variables/CIncEulerVariable.hpp" template CScalarSolver::CScalarSolver(CGeometry* geometry, CConfig* config, bool conservative) @@ -627,6 +628,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol /*--- Flow solution, needed to get density. ---*/ CVariable* flowNodes = solver_container[FLOW_SOL]->GetNodes(); + auto* incFlowNodes = incompressible ? su2staticcast_p(flowNodes) : nullptr; /*--- Store the physical time step ---*/ @@ -655,13 +657,9 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol for (iPoint = 0; iPoint < nPointDomain; iPoint++) { if (Conservative) { if (incompressible) { - /*--- This is temporary and only valid for constant-density problems: - density could also be temperature dependent, but as it is not a part - of the solution vector it's neither stored for previous time steps - nor updated with the solution at the end of each iteration. */ - Density_nM1 = flowNodes->GetDensity(iPoint); - Density_n = flowNodes->GetDensity(iPoint); - Density_nP1 = flowNodes->GetDensity(iPoint); + Density_nM1 = incFlowNodes->GetDensity_time_n1(iPoint); + Density_n = incFlowNodes->GetDensity_time_n(iPoint); + Density_nP1 = incFlowNodes->GetDensity(iPoint); } else { Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0]; Density_n = flowNodes->GetSolution_time_n(iPoint, 0); @@ -686,13 +684,26 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol time discretization scheme (1st- or 2nd-order).---*/ for (iVar = 0; iVar < nVar; iVar++) { + su2double unsteady_term = 0.0; if (first_order) - LinSysRes(iPoint, iVar) += - (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * Volume_nP1 / TimeStep; + unsteady_term = (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * Volume_nP1 / TimeStep; if (second_order) - LinSysRes(iPoint, iVar) += (3.0 * Density_nP1 * U_time_nP1[iVar] - 4.0 * Density_n * U_time_n[iVar] + + unsteady_term = (3.0 * Density_nP1 * U_time_nP1[iVar] - 4.0 * Density_n * U_time_n[iVar] + 1.0 * Density_nM1 * U_time_nM1[iVar]) * Volume_nP1 / (2.0 * TimeStep); + + bool bounded_scalar = false; + if (config->GetKind_Solver() == MAIN_SOLVER::INC_EULER || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES || + config->GetKind_Solver() == MAIN_SOLVER::INC_RANS) { + if (IsSpeciesSolver()) bounded_scalar = config->GetBounded_Species(); + else if (IsTurbSolver()) bounded_scalar = config->GetBounded_Turb(); + } + if (bounded_scalar) { + if (first_order) unsteady_term -= U_time_nP1[iVar] * (Density_nP1 - Density_n) * Volume_nP1 / TimeStep; + if (second_order) unsteady_term -= U_time_nP1[iVar] * (3.0 * Density_nP1 - 4.0 * Density_n + 1.0 * Density_nM1) * Volume_nP1 / (2.0 * TimeStep); + } + + LinSysRes(iPoint, iVar) += unsteady_term; } /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ @@ -721,7 +732,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity(iPoint); // Temporary fix + Density_n = incFlowNodes->GetDensity_time_n(iPoint); else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } @@ -777,7 +788,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity(iPoint); // Temporary fix + Density_n = incFlowNodes->GetDensity_time_n(iPoint); else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } @@ -818,13 +829,9 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol /*--- If this is the SST model, we need to multiply by the density in order to get the conservative variables ---*/ if (incompressible) { - /*--- This is temporary and only valid for constant-density problems: - density could also be temperature dependent, but as it is not a part - of the solution vector it's neither stored for previous time steps - nor updated with the solution at the end of each iteration. */ - Density_nM1 = flowNodes->GetDensity(iPoint); - Density_n = flowNodes->GetDensity(iPoint); - Density_nP1 = flowNodes->GetDensity(iPoint); + Density_nM1 = incFlowNodes->GetDensity_time_n1(iPoint); + Density_n = incFlowNodes->GetDensity_time_n(iPoint); + Density_nP1 = incFlowNodes->GetDensity(iPoint); } else { Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0]; Density_n = flowNodes->GetSolution_time_n(iPoint, 0); @@ -833,20 +840,49 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol } for (iVar = 0; iVar < nVar; iVar++) { + su2double unsteady_term = 0.0; if (first_order) - LinSysRes(iPoint, iVar) += - (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nP1 / TimeStep); + unsteady_term = (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nP1 / TimeStep); if (second_order) - LinSysRes(iPoint, iVar) += - (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) + + unsteady_term = (Density_nP1 * U_time_nP1[iVar] - Density_n * U_time_n[iVar]) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) + (Density_nM1 * U_time_nM1[iVar] - Density_n * U_time_n[iVar]) * (Volume_nM1 / (2.0 * TimeStep)); + + bool bounded_scalar = false; + if (config->GetKind_Solver() == MAIN_SOLVER::INC_EULER || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES || + config->GetKind_Solver() == MAIN_SOLVER::INC_RANS) { + if (IsSpeciesSolver()) bounded_scalar = config->GetBounded_Species(); + else if (IsTurbSolver()) bounded_scalar = config->GetBounded_Turb(); + } + + if (bounded_scalar) { + if (first_order) unsteady_term -= U_time_nP1[iVar] * (Density_nP1 - Density_n) * (Volume_nP1 / TimeStep); + if (second_order) unsteady_term -= U_time_nP1[iVar] * ((Density_nP1 - Density_n) * (3.0 * Volume_nP1 / (2.0 * TimeStep)) + + (Density_nM1 - Density_n) * (Volume_nM1 / (2.0 * TimeStep))); + } + + LinSysRes(iPoint, iVar) += unsteady_term; } /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ if (implicit) { - if (first_order) Jacobian.AddVal2Diag(iPoint, Volume_nP1 / TimeStep); - if (second_order) Jacobian.AddVal2Diag(iPoint, (Volume_nP1 * 3.0) / (2.0 * TimeStep)); - } + bool bounded_scalar = false; + if (config->GetKind_Solver() == MAIN_SOLVER::INC_EULER || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES || + config->GetKind_Solver() == MAIN_SOLVER::INC_RANS) { + if (IsSpeciesSolver()) bounded_scalar = config->GetBounded_Species(); + else if (IsTurbSolver()) bounded_scalar = config->GetBounded_Turb(); + } + su2double diag_factor = 1.0; + if (Conservative) { + if (bounded_scalar) { + if (first_order) diag_factor = Density_n; + if (second_order) diag_factor = (4.0 * Density_n - Density_nM1) / 3.0; + } else { + diag_factor = Density_nP1; + } + } + if (first_order) Jacobian.AddVal2Diag(iPoint, diag_factor * Volume_nP1 / TimeStep); + if (second_order) Jacobian.AddVal2Diag(iPoint, diag_factor * 3.0 * Volume_nP1 / (2.0 * TimeStep)); + } } END_SU2_OMP_FOR @@ -884,4 +920,4 @@ void CScalarSolver::PushSolutionBackInTime(unsigned long TimeIter, nodes->Set_Solution_time_n(); } } -} \ No newline at end of file +} diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 6161e2f0fc3..f75de5b543e 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -347,6 +347,16 @@ class CSolver { */ inline virtual CFluidModel* GetFluidModel(void) const { return nullptr;} + /*! + * \brief Virtual function returning whether this is the species solver. + */ + inline virtual bool IsSpeciesSolver() const { return false; } + + /*! + * \brief Virtual function returning whether this is the turbulent solver. + */ + inline virtual bool IsTurbSolver() const { return false; } + /*! * \brief Get number of linear solver iterations. * \return Number of linear solver iterations. diff --git a/SU2_CFD/include/solvers/CSpeciesSolver.hpp b/SU2_CFD/include/solvers/CSpeciesSolver.hpp index 3a657d97c01..5b0cc08c10a 100644 --- a/SU2_CFD/include/solvers/CSpeciesSolver.hpp +++ b/SU2_CFD/include/solvers/CSpeciesSolver.hpp @@ -69,6 +69,11 @@ class CSpeciesSolver : public CScalarSolver { void Initialize(CGeometry* geometry, CConfig* config, unsigned short iMesh, unsigned short nVar); + /*! + * \brief Virtual function returning whether this is the species solver. + */ + inline bool IsSpeciesSolver() const override { return true; } + /*! * \brief Restart residual and compute gradients. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CTurbSolver.hpp b/SU2_CFD/include/solvers/CTurbSolver.hpp index f2d1c789f7d..22734d15514 100644 --- a/SU2_CFD/include/solvers/CTurbSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSolver.hpp @@ -55,6 +55,11 @@ class CTurbSolver : public CScalarSolver { */ CTurbSolver(CGeometry* geometry, CConfig *config, bool conservative); + /*! + * \brief Virtual function returning whether this is the turbulent solver. + */ + inline bool IsTurbSolver() const override { return true; } + /*! * \brief Impose via the residual the Euler wall boundary condition. * \param[in] geometry - Geometrical definition of the problem. @@ -147,7 +152,7 @@ class CTurbSolver : public CScalarSolver { * \returns The number of extra variables. */ unsigned long RegisterSolutionExtra(bool input, const CConfig* config) final; - + /*! * \brief Compute a suitable under-relaxation parameter to limit the change in the solution variables over * a nonlinear iteration for stability. diff --git a/SU2_CFD/include/variables/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 46be6d3dc2b..bedc421bbeb 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: /*! @@ -87,6 +89,16 @@ class CIncEulerVariable : public CFlowVariable { CIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config); + /*! + * \brief Set all the solution at time level n to the current solution value (including density). + */ + void Set_Solution_time_n() override; + + /*! + * \brief Set all the solution at time level n-1 to the solution at time level n (including density). + */ + void Set_Solution_time_n1() override; + /*! * \brief Set the value of the pressure. * \param[in] iPoint - Point index. @@ -291,4 +303,68 @@ class CIncEulerVariable : public CFlowVariable { for (unsigned long iDim = 0; iDim < nDim; iDim++) Solution(iPoint, iDim+1) = val_vector[iDim]; } + /*! + * \brief Get the density at time level n for dual-time stepping. + * \param[in] iPoint - Point index. + * \return Density at time level n. + */ + inline su2double GetDensity_time_n(unsigned long iPoint) const { return Density_time_n(iPoint); } + + /*! + * \brief Get the density at time level n-1 for dual-time stepping. + * \param[in] iPoint - Point index. + * \return Density at time level n-1. + */ + inline su2double GetDensity_time_n1(unsigned long iPoint) const { return Density_time_n1(iPoint); } + + /*! + * \brief Set the density at time level n for dual-time stepping. + * \param[in] iPoint - Point index. + * \param[in] val_density - Density value. + */ + inline void SetDensity_time_n(unsigned long iPoint, su2double val_density) { Density_time_n(iPoint) = val_density; } + + /*! + * \brief Set the density at time level n-1 for dual-time stepping. + * \param[in] iPoint - Point index. + * \param[in] val_density - Density value. + */ + inline void SetDensity_time_n1(unsigned long iPoint, su2double val_density) { Density_time_n1(iPoint) = val_density; } + + /*! + * \brief Register the density at time n as an AD input variable (for discrete adjoint unsteady). + */ + void RegisterDensity_time_n() override; + + /*! + * \brief Register the density at time n-1 as an AD input variable (for discrete adjoint unsteady). + */ + void RegisterDensity_time_n1() override; + + /*! + * \brief Get the adjoint of density at time n. + * Simplified version: directly extracts the AD derivative of Density_time_n, + * analogous to GetAdjointSolution_time_n for the solution vector. + * The chain-rule contribution d(rho)/d(h) is currently not applied (TODO). + * \param[in] iPoint - Point index. + * \return Adjoint of the density at time n. + */ + inline su2double GetAdjointDensity_time_n(unsigned long iPoint) const override { + AD::Identifier index = AD::GetPassiveIndex(); + AD::SetIndex(index, Density_time_n(iPoint)); + return AD::GetDerivative(index); + } + + /*! + * \brief Get the adjoint of density at time n-1. + * See GetAdjointDensity_time_n. + * \param[in] iPoint - Point index. + * \return Adjoint of the density at time n-1. + */ + inline su2double GetAdjointDensity_time_n1(unsigned long iPoint) const override { + AD::Identifier index = AD::GetPassiveIndex(); + AD::SetIndex(index, Density_time_n1(iPoint)); + return AD::GetDerivative(index); + } + }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 528adc139db..9a60d16c0fd 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -282,12 +282,12 @@ class CVariable { /*! * \brief Set the variable solution at time n. */ - void Set_Solution_time_n(); + virtual void Set_Solution_time_n(); /*! * \brief Set the variable solution at time n-1. */ - void Set_Solution_time_n1(); + virtual void Set_Solution_time_n1(); /*! * \brief Set the variable solution at time n. @@ -2227,6 +2227,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/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 23e5f066f15..dd8c195f8af 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"); @@ -399,6 +408,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){ const auto* Node_Flow = solver[FLOW_SOL]->GetNodes(); + const auto* Node_Flow_Inc = static_cast(Node_Flow); const CVariable* Node_Turb = nullptr; if (config->GetKind_Turb_Model() != TURB_MODEL::NONE) Node_Turb = solver[TURB_SOL]->GetNodes(); const CVariable* Node_Heat = nullptr; @@ -439,6 +449,14 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("PRESSURE_COEFF", iPoint, (Node_Flow->GetPressure(iPoint) - solver[FLOW_SOL]->GetPressure_Inf())/factor); SetVolumeOutputValue("DENSITY", iPoint, Node_Flow->GetDensity(iPoint)); + // Load density at previous timesteps for restart (unsteady with non-constant density) + if (config->GetTime_Domain() && config->GetKind_FluidModel() != CONSTANT_DENSITY) { + SetVolumeOutputValue("DENSITY_TIME_N", iPoint, Node_Flow_Inc->GetDensity_time_n(iPoint)); + if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) { + SetVolumeOutputValue("DENSITY_TIME_N1", iPoint, Node_Flow_Inc->GetDensity_time_n1(iPoint)); + } + } + if (config->GetKind_Solver() == MAIN_SOLVER::INC_RANS || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES){ SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint)); SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint)); diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 3194e28b6d3..22d1ffb9bd3 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -168,11 +168,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) { @@ -373,6 +380,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 @@ -388,6 +397,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/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 1cb4a1951bd..3f822d42070 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -966,6 +966,8 @@ 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 || config->GetKind_Centered_Flow() == CENTERED::LD2) && (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); /*--- Set the primitive variables ---*/ @@ -974,6 +976,20 @@ void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_ SU2_OMP_ATOMIC ErrorCounter += SetPrimitive_Variables(solver_container, config); + /*--- Provide a valid initial Density_time_n for restart or step 0 since + it was 0.0 before SetPrimitive_Variables evaluated the exact field. + Checking GetDensity_time_n(0) == 0.0 ensures we only evaluate this + once when the history arrays are uninitialized. ---*/ + if (dual_time && nPoint > 0 && nodes->GetDensity_time_n(0) == 0.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 { @@ -1183,7 +1199,7 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co if (LD2_Scheme) { numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(jPoint)); - if (!geometry->nodes->GetPeriodicBoundary(iPoint) || (geometry->nodes->GetPeriodicBoundary(iPoint) + if (!geometry->nodes->GetPeriodicBoundary(iPoint) || (geometry->nodes->GetPeriodicBoundary(iPoint) && !geometry->nodes->GetPeriodicBoundary(jPoint))) { numerics->SetCoord(geometry->nodes->GetCoord(iPoint), geometry->nodes->GetCoord(jPoint)); } else { @@ -2874,14 +2890,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 = 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, V_time_nM1, U_time_nM1); - V2U(Density, 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 @@ -2929,8 +2948,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the conservative variables. ---*/ V_time_n = nodes->GetSolution_time_n(iPoint); - Density = nodes->GetDensity(iPoint); - V2U(Density, V_time_n, U_time_n); + su2double Density_n = nodes->GetDensity_time_n(iPoint); + V2U(Density_n, V_time_n, U_time_n); GridVel_i = geometry->nodes->GetGridVel(iPoint); @@ -2984,8 +3003,8 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Compute the GCL component of the source term for node i ---*/ V_time_n = nodes->GetSolution_time_n(iPoint); - Density = nodes->GetDensity(iPoint); - V2U(Density, V_time_n, U_time_n); + 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; @@ -3011,14 +3030,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 = 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, V_time_nM1, U_time_nM1); - V2U(Density, 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 and n+1. In the case of dynamically deforming diff --git a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp index 10fb3c49ddb..e25e2fc0880 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -83,7 +83,12 @@ 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))); } diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp index 355a1587488..bab2b45d907 100644 --- a/SU2_CFD/src/variables/CIncEulerVariable.cpp +++ b/SU2_CFD/src/variables/CIncEulerVariable.cpp @@ -27,6 +27,7 @@ #include "../../include/variables/CIncEulerVariable.hpp" #include "../../include/fluid/CFluidModel.hpp" +#include "../../../Common/include/parallelization/omp_structure.hpp" CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *velocity, su2double enthalpy, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) @@ -58,6 +59,10 @@ CIncEulerVariable::CIncEulerVariable(su2double pressure, const su2double *veloci if (dual_time) { Solution_time_n = Solution; Solution_time_n1 = Solution; + /*--- Allocate density storage for time levels. + Note: Actual density values will be set after SetPrimVar is called in Preprocessing. ---*/ + Density_time_n.resize(nPoint) = su2double(0.0); + Density_time_n1.resize(nPoint) = su2double(0.0); } if (config->GetKind_Streamwise_Periodic() != ENUM_STREAMWISE_PERIODIC::NONE) { @@ -127,3 +132,38 @@ 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 + } +} + +void CIncEulerVariable::RegisterDensity_time_n() { + RegisterContainer(true, Density_time_n); +} + +void CIncEulerVariable::RegisterDensity_time_n1() { + RegisterContainer(true, Density_time_n1); +} diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index ac396113b75..6339a25518b 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -55,7 +55,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()) diff --git a/TestCases/flamelet/09_laminar_premixed_ch4_unsteady/lam_prem_ch4_unsteady.cfg b/TestCases/flamelet/09_laminar_premixed_ch4_unsteady/lam_prem_ch4_unsteady.cfg new file mode 100644 index 00000000000..aef2a4179da --- /dev/null +++ b/TestCases/flamelet/09_laminar_premixed_ch4_unsteady/lam_prem_ch4_unsteady.cfg @@ -0,0 +1,157 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Laminar premixed transient flame propagating in a channel % +% Author: Nijso Beishuizen % +% Institution: TU Eindhoven % +% Date: 30/04/2026 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= INC_NAVIER_STOKES +KIND_TURB_MODEL= NONE +MATH_PROBLEM= DIRECT +RESTART_SOL= YES +% +% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------% +% +%INC_DENSITY_MODEL= FLAMELET +INC_DENSITY_MODEL= VARIABLE +INC_DENSITY_INIT= 1.00 +INC_VELOCITY_INIT= (0.5, 0.0, 0.0 ) +INC_TEMPERATURE_INIT= 300.0 +INC_NONDIM= DIMENSIONAL +% +% -------------------- FLUID MODEL --------------------------------------- % +% +FLUID_MODEL= FLUID_FLAMELET +PREFERENTIAL_DIFFUSION= NO +FILENAMES_INTERPOLATOR= (LUT_ch4.drg) +%FILENAMES_INTERPOLATOR= (LUT_methane_3D.drg) +CONTROLLING_VARIABLE_NAMES= (ProgressVariable, EnthalpyTot, MixtureFraction) +CONTROLLING_VARIABLE_SOURCE_NAMES= (ProdRateTot_PV, NULL, NULL) +% +% -------------------- SCALAR TRANSPORT ---------------------------------------% +% +KIND_SCALAR_MODEL= FLAMELET +DIFFUSIVITY_MODEL= FLAMELET +VISCOSITY_MODEL= FLAMELET +CONDUCTIVITY_MODEL= FLAMELET +%FLAME_INIT_METHOD= FLAME_FRONT +FLAME_INIT_METHOD= NONE + +FLAME_INIT= (0.009, 0.00, 0.00, 1.0, 0.0, 0.0, 5.0e-4, 0.1) +% # progvar, enthalpy, mixture fraction phi=0.75 +SPECIES_INIT= (0.0, -193150, 0.04197) +CONV_NUM_METHOD_SPECIES= BOUNDED_SCALAR +%CONV_NUM_METHOD_SPECIES= SCALAR_UPWIND +MUSCL_SPECIES= YES +SLOPE_LIMITER_SPECIES= NONE +TIME_DISCRE_SPECIES= EULER_IMPLICIT +% SCALAR CLIPPING +SPECIES_CLIPPING= YES +SPECIES_CLIPPING_MIN= 0.000 -1e6 0.03 +SPECIES_CLIPPING_MAX= 0.300 +1e5 0.05 +% +MARKER_INLET_SPECIES= (inlet, 0, -193154.0, 0.04197) +CFL_REDUCTION_SPECIES= 1.0 +MARKER_SPECIES_STRONG_BC= (inlet, outlet) +LOOKUP_NAMES= (MolarWeightMix, Conductivity, Cp) +%USER_SCALAR_NAMES= (CO, NOx) +%USER_SOURCE_NAMES= ( ProdRatePos_Y-CO, ProdRateNegScaled_Y-CO , ProdRateTot_Y-NOx, NULL) +% +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X= 0.25 +REF_ORIGIN_MOMENT_Y= 0.00 +REF_ORIGIN_MOMENT_Z= 0.00 +REF_LENGTH= 1.0 +REF_AREA= 1.0 +% +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_SYM= (symmetry_bottom) +MARKER_EULER= wall +INC_INLET_TYPE= VELOCITY_INLET + +% switch from 'steady' flame to 'propagating with 0.1 m/s' +%MARKER_INLET= (inlet, 300.0, 0.227, 1.0, 0.0, 0.0) +MARKER_INLET= (inlet, 300.0, 0.127, 1.0, 0.0, 0.0) + +INC_OUTLET_TYPE= PRESSURE_OUTLET +INC_INLET_DAMPING= 0.1 +INC_OUTLET_DAMPING= 0.1 +MARKER_OUTLET= (outlet, 0.0) +MARKER_PLOTTING= ( inlet ) +MARKER_MONITORING= ( inlet ) +MARKER_ANALYZE= ( inlet,outlet ) +MARKER_ANALYZE_AVERAGE= AREA +% +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CFL_NUMBER= 10000 +CFL_ADAPT= NO +OUTPUT_WRT_FREQ= 10,10 +% +% ------------------------- UNSTEADY SIMULATION -------------------------------% +% +TIME_DOMAIN= YES +% +TIME_MARCHING= DUAL_TIME_STEPPING-1ST_ORDER +%TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER +% +RESTART_ITER=1 +TIME_ITER= 50 +TIME_STEP= 2e-4 +% time= 0.1 s +MAX_TIME= 0.1 +% +UNST_CFL_NUMBER= 0.0 +INNER_ITER= 100 +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +% no impact on flame speed +LINEAR_SOLVER_ERROR= 1e-3 +LINEAR_SOLVER_ITER= 10 +% +% -------------------------- MULTIGRID PARAMETERS -----------------------------% +% +MGLEVEL= 0 +% +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= FDS +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +TIME_DISCRE_FLOW= EULER_IMPLICIT +% +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_RESIDUAL_MINVAL= -10.0 +CONV_FIELD= RMS_ProgressVariable, RMS_VELOCITY-X, RMS_MixtureFraction +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-6 +SCREEN_OUTPUT= TIME_ITER INNER_ITER RMS_VELOCITY-X RMS_PRESSURE RMS_ProgressVariable RMS_EnthalpyTot RMS_MixtureFraction +HISTORY_OUTPUT= RMS_RES AERO_COEFF FLOW_COEFF FLOW_COEFF_SURF +VOLUME_OUTPUT= SOLUTION PRIMITIVE SOURCE RESIDUAL SENSITIVITY LOOKUP TIMESTEP +% +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= mesh_finer.su2 +MESH_OUT_FILENAME= mesh_out +SOLUTION_FILENAME= solution +RESTART_FILENAME= restart +OUTPUT_FILES= (RESTART,PARAVIEW) +%OUTPUT_FILES= (RESTART) +TABULAR_FORMAT= CSV +CONV_FILENAME= history +VOLUME_FILENAME= ch4_flame_cfd +SURFACE_FILENAME= surface_flow +WRT_PERFORMANCE= YES +SCREEN_WRT_FREQ_INNER= 1 +SCREEN_WRT_FREQ_OUTER= 1 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 7f10753db6e..25f1368a0a4 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -86,6 +86,14 @@ def main(): flame_init_methods.new_output = True test_list.append(flame_init_methods) + # 2D laminar premixed ch4-air flame, transient flame propagation + cfd_flamelet_ch4_unsteady = TestCase('cfd_flamelet_ch4_unsteady') + cfd_flamelet_ch4_unsteady.cfg_dir = "flamelet/09_laminar_premixed_ch4_unsteady" + cfd_flamelet_ch4_unsteady.cfg_file = "lam_prem_ch4_unsteady.cfg" + cfd_flamelet_ch4_unsteady.test_iter = 5 + cfd_flamelet_ch4_unsteady.test_vals = [-8.036958, -8.372668, -1.842800, -9.388446] + test_list.append(cfd_flamelet_ch4_unsteady) + ######################### ## NEMO solver ### #########################