diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 65a84ed437d..849007f62b3 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5707,6 +5707,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/solvers/CScalarSolver.inl b/SU2_CFD/include/solvers/CScalarSolver.inl index b02cfacf064..8365a720194 100644 --- a/SU2_CFD/include/solvers/CScalarSolver.inl +++ b/SU2_CFD/include/solvers/CScalarSolver.inl @@ -634,12 +634,8 @@ 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_nM1 = flowNodes->GetDensity_time_n1(iPoint); + Density_n = flowNodes->GetDensity_time_n(iPoint); Density_nP1 = flowNodes->GetDensity(iPoint); } else { Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0]; @@ -700,7 +696,7 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity(iPoint); // Temporary fix + Density_n = flowNodes->GetDensity_time_n(iPoint); else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } @@ -756,11 +752,11 @@ void CScalarSolver::SetResidual_DualTime(CGeometry* geometry, CSol if (Conservative) { if (incompressible) - Density_n = flowNodes->GetDensity(iPoint); // Temporary fix + Density_n = flowNodes->GetDensity_time_n(iPoint); else Density_n = flowNodes->GetSolution_time_n(iPoint, 0); } - + // nijso: what about second order? for (iVar = 0; iVar < nVar; iVar++) LinSysRes(iPoint, iVar) += Density_n * U_time_n[iVar] * Residual_GCL; } END_SU2_OMP_FOR @@ -797,12 +793,8 @@ 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_nM1 = flowNodes->GetDensity_time_n1(iPoint); + Density_n = flowNodes->GetDensity_time_n(iPoint); Density_nP1 = flowNodes->GetDensity(iPoint); } else { Density_nM1 = flowNodes->GetSolution_time_n1(iPoint)[0]; @@ -862,4 +854,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/CIncEulerVariable.hpp b/SU2_CFD/include/variables/CIncEulerVariable.hpp index 9cc13f79986..7d43605bb36 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 bc4bc55e57c..9c931d430c8 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -284,12 +284,12 @@ class CVariable { /*! * \brief Set the variable solution at time n. */ - void Set_Solution_time_n(); + virtual void Set_Solution_time_n(); /*! * \brief Set the variable solution at time n-1. */ - void Set_Solution_time_n1(); + virtual void Set_Solution_time_n1(); /*! * \brief Set the variable solution at time n. @@ -2182,6 +2182,48 @@ class CVariable { */ void RegisterUserDefinedSource(); + /*! + * \brief Get the stored density at time n (incompressible only). + * \param[in] iPoint - Point index. + * \return Density at time n; 0 for compressible nodes. + */ + inline virtual su2double GetDensity_time_n(unsigned long /*iPoint*/) const { return 0.0; } + + /*! + * \brief Get the stored density at time n-1 (incompressible only). + * \param[in] iPoint - Point index. + * \return Density at time n-1; 0 for compressible nodes. + */ + inline virtual su2double GetDensity_time_n1(unsigned long /*iPoint*/) const { return 0.0; } + + /*! + * \brief Register the density at time n as an AD input (no-op for non-incompressible nodes). + */ + inline virtual void RegisterDensity_time_n() {} + + /*! + * \brief Register the density at time n-1 as an AD input (no-op for non-incompressible nodes). + */ + inline virtual void RegisterDensity_time_n1() {} + + /*! + * \brief Get the adjoint values of the density at time n. + * Implemented analogously to GetAdjointSolution_time_n but for a scalar + * (Density_time_n is a VectorType, not a MatrixType). + * Currently ignores the chain-rule d(rho)/d(h) dependence (TODO). + * \param[in] iPoint - Point index. + * \return Adjoint value of the density at time n. + */ + inline virtual su2double GetAdjointDensity_time_n(unsigned long /*iPoint*/) const { return 0.0; } + + /*! + * \brief Get the adjoint values of the density at time n-1. + * See GetAdjointDensity_time_n. + * \param[in] iPoint - Point index. + * \return Adjoint value of the density at time n-1. + */ + inline virtual su2double GetAdjointDensity_time_n1(unsigned long /*iPoint*/) const { return 0.0; } + /*! * \brief Set the adjoint values of the solution. * \param[in] adj_sol - The adjoint values of the solution. diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index b13edb7732f..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"); @@ -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_Heat = nullptr; const CVariable* Node_Rad = nullptr; auto* Node_Geo = geometry->nodes; @@ -437,6 +447,14 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve SetVolumeOutputValue("PRESSURE_COEFF", iPoint, (Node_Flow->GetPressure(iPoint) - solver[FLOW_SOL]->GetPressure_Inf())/factor); SetVolumeOutputValue("DENSITY", iPoint, Node_Flow->GetDensity(iPoint)); + // Load density at previous timesteps for restart (unsteady with non-constant density) + if (config->GetTime_Domain() && config->GetKind_FluidModel() != CONSTANT_DENSITY) { + SetVolumeOutputValue("DENSITY_TIME_N", iPoint, Node_Flow_Inc->GetDensity_time_n(iPoint)); + if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) { + SetVolumeOutputValue("DENSITY_TIME_N1", iPoint, Node_Flow_Inc->GetDensity_time_n1(iPoint)); + } + } + if (config->GetKind_Solver() == MAIN_SOLVER::INC_RANS || config->GetKind_Solver() == MAIN_SOLVER::INC_NAVIER_STOKES){ SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint)); SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint)); 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/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 5ba375e7489..8715e0e2507 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -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 { @@ -2835,14 +2851,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 @@ -2890,8 +2909,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); @@ -2945,8 +2964,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; @@ -2972,14 +2991,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 d2143ba8d1d..c19c71c1291 100644 --- a/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp +++ b/SU2_CFD/src/solvers/CSpeciesFlameletSolver.cpp @@ -81,7 +81,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))); } @@ -102,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]); + } } } @@ -194,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 6d77af4454d..a8809504afe 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 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++) { + 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 7de950de622..5cde8161b5e 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())