diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 44f5afe1f893..cfde0c3b1f1e 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -708,8 +708,8 @@ class CSysMatrix { * \param[in] val_block - Block to add to the diagonal of the matrix. * \param[in] alpha - Scale factor. */ - template - inline void SetBlock2Diag(unsigned long block_i, const OtherType* const* val_block, OtherType alpha = 1.0) { + template + inline void SetBlock2Diag(unsigned long block_i, const OtherType& val_block, T alpha = 1.0) { auto mat_ii = &matrix[dia_ptr[block_i]*nVar*nEqn]; @@ -723,8 +723,8 @@ class CSysMatrix { /*! * \brief Non overwrite version of SetBlock2Diag, also with scaling. */ - template - inline void AddBlock2Diag(unsigned long block_i, const OtherType* const* val_block, OtherType alpha = 1.0) { + template + inline void AddBlock2Diag(unsigned long block_i, const OtherType& val_block, T alpha = 1.0) { SetBlock2Diag(block_i, val_block, alpha); } @@ -732,8 +732,8 @@ class CSysMatrix { * \brief Short-hand to AddBlock2Diag with alpha = -1, i.e. subtracts from the current diagonal. */ template - inline void SubtractBlock2Diag(unsigned long block_i, const OtherType* const* val_block) { - AddBlock2Diag(block_i, val_block, OtherType(-1)); + inline void SubtractBlock2Diag(unsigned long block_i, const OtherType& val_block) { + AddBlock2Diag(block_i, val_block, -1.0); } /*! @@ -748,6 +748,18 @@ class CSysMatrix { matrix[dia_ptr[block_i]*nVar*nVar + iVar*(nVar+1)] += PassiveAssign(val_matrix); } + /*! + * \brief Adds the specified value to the diagonal of the (i, i) subblock + * of the matrix-by-blocks structure. + * \param[in] block_i - Diagonal index. + * \param[in] iVar - Variable index. + * \param[in] val - Value to add to the diagonal elements of A(i, i). + */ + template + inline void AddVal2Diag(unsigned long block_i, unsigned long iVar, OtherType val) { + matrix[dia_ptr[block_i]*nVar*nVar + iVar*(nVar+1)] += PassiveAssign(val); + } + /*! * \brief Sets the specified value to the diagonal of the (i, i) subblock * of the matrix-by-blocks structure. diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 2c8efec9691c..5eb2935751d4 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -114,11 +114,6 @@ class CEulerSolver : public CFVMFlowSolverBase { vector FluidModel; /*!< \brief fluid model used in the solver. */ - unsigned long ErrorCounter = 0; /*!< \brief Counter for number of un-physical states. */ - - su2double Global_Delta_Time = 0.0, /*!< \brief Time-step for TIME_STEPPING time marching strategy. */ - Global_Delta_UnstTimeND = 0.0; /*!< \brief Unsteady time step for the dual time strategy. */ - /*--- Turbomachinery Solver Variables ---*/ su2double ***AverageFlux = nullptr, @@ -160,12 +155,6 @@ class CEulerSolver : public CFVMFlowSolverBase { /*--- End of Turbomachinery Solver Variables ---*/ - /*! - * \brief Generic implementation of explicit iterations (RK, Classic RK and EULER). - */ - template - void Explicit_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iRKStep); - /*! * \brief Preprocessing actions common to the Euler and NS solvers. * \param[in] geometry - Geometrical definition of the problem. @@ -239,7 +228,7 @@ class CEulerSolver : public CFVMFlowSolverBase { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void SetMax_Eigenvalue(CGeometry *geometry, CConfig *config); + void SetMax_Eigenvalue(CGeometry *geometry, const CConfig *config); /*! * \brief Compute the undivided laplacian for the solution. @@ -266,11 +255,10 @@ class CEulerSolver : public CFVMFlowSolverBase { * \brief Compute the velocity^2, SoundSpeed, Pressure, Enthalpy, Viscosity. * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. - * \param[in] Output - boolean to determine whether to print output. * \return - The number of non-physical points. */ virtual unsigned long SetPrimitive_Variables(CSolver **solver_container, - CConfig *config, bool Output); + const CConfig *config); /*! * \brief Set gradients of coefficients for fixed CL mode @@ -285,6 +273,13 @@ class CEulerSolver : public CFVMFlowSolverBase { */ void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) final; + /*! + * \brief Set the solver nondimensionalization. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void SetNondimensionalization(CConfig *config, unsigned short iMesh); + public: /*! * \brief Constructor of the class. @@ -305,13 +300,6 @@ class CEulerSolver : public CFVMFlowSolverBase { */ ~CEulerSolver(void) override; - /*! - * \brief Set the solver nondimensionalization. - * \param[in] config - Definition of the particular problem. - * \param[in] iMesh - Index of the mesh in multigrid computations. - */ - void SetNondimensionalization(CConfig *config, unsigned short iMesh) final; - /*! * \brief Compute the pressure at the infinity. * \return Value of the pressure at the infinity. @@ -459,8 +447,7 @@ class CEulerSolver : public CFVMFlowSolverBase { * \param[in,out] preconditioner - The preconditioner matrix, must be allocated outside. */ void SetPreconditioner(const CConfig *config, unsigned long iPoint, - su2double delta, su2double** preconditioner) const; - using CSolver::SetPreconditioner; /*--- Silence warning. ---*/ + su2double delta, su2activematrix& preconditioner) const; /*! * \brief Parallelization of Undivided Laplacian. @@ -1084,22 +1071,6 @@ class CEulerSolver : public CFVMFlowSolverBase { */ void UpdateCustomBoundaryConditions(CGeometry **geometry_container, CConfig *config) final; - /*! - * \brief Set the total residual adding the term that comes from the Dual Time Strategy. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] config - Definition of the particular problem. - * \param[in] iRKStep - Current step of the Runge-Kutta iteration. - * \param[in] iMesh - Index of the mesh in multigrid computations. - * \param[in] RunTime_EqSystem - System of equations which is going to be solved. - */ - void SetResidual_DualTime(CGeometry *geometry, - CSolver **solver_container, - CConfig *config, - unsigned short iRKStep, - unsigned short iMesh, - unsigned short RunTime_EqSystem) final; - /*! * \brief Load a solution from a restart file. * \param[in] geometry - Geometrical definition of the problem. @@ -1130,7 +1101,7 @@ class CEulerSolver : public CFVMFlowSolverBase { * \brief Set the solution using the Freestream values. * \param[in] config - Definition of the particular problem. */ - void SetFreeStream_Solution(CConfig *config) final; + void SetFreeStream_Solution(const CConfig *config) final; /*! * \brief Initilize turbo containers. diff --git a/SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp b/SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp index aa41b64a18f9..bff89172fe9d 100644 --- a/SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp +++ b/SU2_CFD/include/solvers/CFEM_DG_EulerSolver.hpp @@ -306,7 +306,6 @@ class CFEM_DG_EulerSolver : public CSolver { void SetNondimensionalization(CConfig *config, unsigned short iMesh, const bool writeOutput); - using CSolver::SetNondimensionalization; /*! * \brief Get a pointer to the vector of the solution degrees of freedom. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index 75492c114085..836f01659688 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -56,6 +56,11 @@ class CFVMFlowSolverBase : public CSolver { su2double StrainMag_Max; /*!< \brief Maximum Strain Rate magnitude. */ su2double Omega_Max; /*!< \brief Maximum Omega. */ + su2double Global_Delta_Time = 0.0, /*!< \brief Time-step for TIME_STEPPING time marching strategy. */ + Global_Delta_UnstTimeND = 0.0; /*!< \brief Unsteady time step for the dual time strategy. */ + + unsigned long ErrorCounter = 0; /*!< \brief Counter for number of un-physical states. */ + /*! * \brief Auxilary types to store common aero coefficients (avoids repeating oneself so much). */ @@ -244,6 +249,748 @@ class CFVMFlowSolverBase : public CSolver { */ inline virtual void InstantiateEdgeNumerics(const CSolver* const* solvers, const CConfig* config) {} + /*! + * \brief Generic implementation to compute the time step based on CFL and conv/visc eigenvalues. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] Iteration - Value of the current iteration. + * \tparam SoundSpeedFunc - Function object to compute speed of sound. + * \tparam LambdaViscFunc - Function object to compute the viscous lambda. + * \note Both functors need to implement (nodes,iPoint,jPoint) for edges, and (nodes,iPoint) for vertices. + */ + template + FORCEINLINE void SetTime_Step_impl(const SoundSpeedFunc& soundSpeed, + const LambdaViscFunc& lambdaVisc, + CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh, + unsigned long Iteration) { + + const bool viscous = config->GetViscous(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool time_stepping = (config->GetTime_Marching() == TIME_STEPPING); + const bool dual_time = (config->GetTime_Marching() == DT_STEPPING_1ST) || + (config->GetTime_Marching() == DT_STEPPING_2ND); + const su2double K_v = 0.25; + + /*--- Init thread-shared variables to compute min/max values. + * Critical sections are used for this instead of reduction + * clauses for compatibility with OpenMP 2.0 (Windows...). ---*/ + + SU2_OMP_MASTER + { + Min_Delta_Time = 1e30; + Max_Delta_Time = 0.0; + Global_Delta_UnstTimeND = 1e30; + } + SU2_OMP_BARRIER + + /*--- Loop domain points. ---*/ + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) { + + /*--- Set maximum eigenvalues to zero. ---*/ + + nodes->SetMax_Lambda_Inv(iPoint,0.0); + + if (viscous) + nodes->SetMax_Lambda_Visc(iPoint,0.0); + + /*--- Loop over the neighbors of point i. ---*/ + + for (unsigned short iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) + { + auto jPoint = geometry->nodes->GetPoint(iPoint,iNeigh); + + auto iEdge = geometry->nodes->GetEdge(iPoint,iNeigh); + auto Normal = geometry->edges->GetNormal(iEdge); + auto Area2 = GeometryToolbox::SquaredNorm(nDim, Normal); + + /*--- Mean Values ---*/ + + su2double Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); + su2double Mean_SoundSpeed = soundSpeed(*nodes, iPoint, jPoint) * sqrt(Area2); + + /*--- Adjustment for grid movement ---*/ + + if (dynamic_grid) { + const su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); + const su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); + + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Mean_ProjVel -= 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * Normal[iDim]; + } + + /*--- Inviscid contribution ---*/ + + su2double Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; + nodes->AddMax_Lambda_Inv(iPoint,Lambda); + + /*--- Viscous contribution ---*/ + + if (!viscous) continue; + + Lambda = lambdaVisc(*nodes, iPoint, jPoint) * Area2; + nodes->AddMax_Lambda_Visc(iPoint, Lambda); + } + + } + + /*--- Loop boundary edges ---*/ + + for (unsigned short iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && + (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + + /*--- Point identification, Normal vector and area ---*/ + + auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (!geometry->nodes->GetDomain(iPoint)) continue; + + auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + auto Area2 = GeometryToolbox::SquaredNorm(nDim, Normal); + + /*--- Mean Values ---*/ + + su2double ProjVel = nodes->GetProjVel(iPoint,Normal); + su2double SoundSpeed = soundSpeed(*nodes, iPoint) * sqrt(Area2); + + /*--- Adjustment for grid movement ---*/ + + if (dynamic_grid) { + ProjVel -= GeometryToolbox::DotProduct(nDim, Normal, geometry->nodes->GetGridVel(iPoint)); + } + + /*--- Inviscid contribution ---*/ + + su2double Lambda = fabs(ProjVel) + SoundSpeed; + nodes->AddMax_Lambda_Inv(iPoint, Lambda); + + /*--- Viscous contribution ---*/ + + if (!viscous) continue; + + Lambda = lambdaVisc(*nodes,iPoint) * Area2; + nodes->AddMax_Lambda_Visc(iPoint, Lambda); + } + } + } + + /*--- Each element uses their own speed, steady state simulation. ---*/ + { + /*--- Thread-local variables for min/max reduction. ---*/ + su2double minDt = 1e30, maxDt = 0.0; + + SU2_OMP(for schedule(static,omp_chunk_size) nowait) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + + su2double Vol = geometry->nodes->GetVolume(iPoint); + + if (Vol != 0.0) { + su2double Local_Delta_Time = nodes->GetLocalCFL(iPoint)*Vol / nodes->GetMax_Lambda_Inv(iPoint); + + if(viscous) { + su2double dt_visc = nodes->GetLocalCFL(iPoint)*K_v*Vol*Vol / nodes->GetMax_Lambda_Visc(iPoint); + Local_Delta_Time = min(Local_Delta_Time, dt_visc); + } + + minDt = min(minDt, Local_Delta_Time); + maxDt = max(maxDt, Local_Delta_Time); + + nodes->SetDelta_Time(iPoint, min(Local_Delta_Time, config->GetMax_DeltaTime())); + } + else { + nodes->SetDelta_Time(iPoint,0.0); + } + } + /*--- Min/max over threads. ---*/ + SU2_OMP_CRITICAL + { + Min_Delta_Time = min(Min_Delta_Time, minDt); + Max_Delta_Time = max(Max_Delta_Time, maxDt); + Global_Delta_Time = Min_Delta_Time; + } + SU2_OMP_BARRIER + } + + /*--- Compute the min/max dt (in parallel, now over mpi ranks). ---*/ + + SU2_OMP_MASTER + if (config->GetComm_Level() == COMM_FULL) { + su2double rbuf_time; + SU2_MPI::Allreduce(&Min_Delta_Time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + Min_Delta_Time = rbuf_time; + + SU2_MPI::Allreduce(&Max_Delta_Time, &rbuf_time, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + Max_Delta_Time = rbuf_time; + } + SU2_OMP_BARRIER + + /*--- For exact time solution use the minimum delta time of the whole mesh. ---*/ + if (time_stepping) { + + /*--- If the unsteady CFL is set to zero, it uses the defined unsteady time step, + * otherwise it computes the time step based on the unsteady CFL. ---*/ + + SU2_OMP_MASTER + { + if (config->GetUnst_CFL() == 0.0) { + Global_Delta_Time = config->GetDelta_UnstTime(); + } + else { + Global_Delta_Time = Min_Delta_Time; + } + Max_Delta_Time = Global_Delta_Time; + + config->SetDelta_UnstTimeND(Global_Delta_Time); + } + SU2_OMP_BARRIER + + /*--- Sets the regular CFL equal to the unsteady CFL. ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + nodes->SetLocalCFL(iPoint, config->GetUnst_CFL()); + nodes->SetDelta_Time(iPoint, Global_Delta_Time); + } + + } + + /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is diferent from 0. ---*/ + + if ((dual_time) && (Iteration == 0) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) { + + /*--- Thread-local variable for reduction. ---*/ + su2double glbDtND = 1e30; + + SU2_OMP(for schedule(static,omp_chunk_size) nowait) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + glbDtND = min(glbDtND, config->GetUnst_CFL()*Global_Delta_Time / nodes->GetLocalCFL(iPoint)); + } + SU2_OMP_CRITICAL + Global_Delta_UnstTimeND = min(Global_Delta_UnstTimeND, glbDtND); + SU2_OMP_BARRIER + + SU2_OMP_MASTER + { + SU2_MPI::Allreduce(&Global_Delta_UnstTimeND, &glbDtND, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); + Global_Delta_UnstTimeND = glbDtND; + + config->SetDelta_UnstTimeND(Global_Delta_UnstTimeND); + } + SU2_OMP_BARRIER + } + + /*--- The pseudo local time (explicit integration) cannot be greater than the physical time ---*/ + + if (dual_time && !implicit) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + su2double dt = min((2.0/3.0)*config->GetDelta_UnstTimeND(), nodes->GetDelta_Time(iPoint)); + nodes->SetDelta_Time(iPoint, dt); + } + } + } + + /*! + * \brief Compute the max eigenvalue, gemeric implementation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \tparam SoundSpeedFunc - Function object to compute speed of sound. + * \note Functor needs to implement (nodes,iPoint,jPoint) for edges, and (nodes,iPoint) for vertices. + */ + template + FORCEINLINE void SetMax_Eigenvalue_impl(const SoundSpeedFunc& soundSpeed, CGeometry *geometry, const CConfig *config) { + + /*--- Loop domain points. ---*/ + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { + + /*--- Set eigenvalues to zero. ---*/ + nodes->SetLambda(iPoint,0.0); + + /*--- Loop over the neighbors of point i. ---*/ + for (unsigned short iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) + { + auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + + auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); + auto Normal = geometry->edges->GetNormal(iEdge); + su2double Area = GeometryToolbox::Norm(nDim, Normal); + + /*--- Mean Values ---*/ + + su2double Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); + su2double Mean_SoundSpeed = soundSpeed(*nodes, iPoint, jPoint) * Area; + + /*--- Adjustment for grid movement ---*/ + + if (dynamic_grid) { + const su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); + const su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); + + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Mean_ProjVel -= 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * Normal[iDim]; + } + + /*--- Inviscid contribution ---*/ + + nodes->AddLambda(iPoint, fabs(Mean_ProjVel) + Mean_SoundSpeed); + } + } + + /*--- Loop boundary edges ---*/ + + for (unsigned short iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && + (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (unsigned long iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + + /*--- Point identification, Normal vector and area ---*/ + + auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (!geometry->nodes->GetDomain(iPoint)) continue; + + auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + su2double Area = GeometryToolbox::Norm(nDim, Normal); + + /*--- Mean Values ---*/ + + su2double Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); + su2double Mean_SoundSpeed = soundSpeed(*nodes, iPoint) * Area; + + /*--- Adjustment for grid movement ---*/ + + if (dynamic_grid) { + Mean_ProjVel -= GeometryToolbox::DotProduct(nDim, Normal, geometry->nodes->GetGridVel(iPoint)); + } + + /*--- Inviscid contribution ---*/ + + nodes->AddLambda(iPoint, fabs(Mean_ProjVel) + Mean_SoundSpeed); + } + } + } + + /*--- Correct the eigenvalue values across any periodic boundaries. ---*/ + + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_MAX_EIG); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_MAX_EIG); + } + + /*--- MPI parallelization ---*/ + + InitiateComms(geometry, config, MAX_EIGENVALUE); + CompleteComms(geometry, config, MAX_EIGENVALUE); + } + + /*! + * \brief Compute the dissipation sensor for centered schemes. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \tparam SensVarFunc - Function object implementing (nodes, iPoint) to return the sensor variable, e.g. pressure. + */ + template + FORCEINLINE void SetCentered_Dissipation_Sensor_impl(const SensVarFunc& sensVar, + CGeometry *geometry, const CConfig *config) { + + /*--- We can access memory more efficiently if there are no periodic boundaries. ---*/ + + const bool isPeriodic = (config->GetnMarker_Periodic() > 0); + + /*--- Loop domain points. ---*/ + + SU2_OMP_FOR_DYN(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { + + const bool boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint); + const su2double sensVar_i = sensVar(*nodes, iPoint); + + /*--- Initialize. ---*/ + iPoint_UndLapl[iPoint] = 0.0; + jPoint_UndLapl[iPoint] = 0.0; + + /*--- Loop over the neighbors of point i. ---*/ + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) + { + bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); + + /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ + if (boundary_i && !boundary_j) continue; + + su2double sensVar_j = sensVar(*nodes, jPoint); + + /*--- Dissipation sensor, add variable difference and variable sum. ---*/ + iPoint_UndLapl[iPoint] += sensVar_j - sensVar_i; + jPoint_UndLapl[iPoint] += sensVar_j + sensVar_i; + } + + if (!isPeriodic) { + /*--- Every neighbor is accounted for, sensor can be computed. ---*/ + nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); + } + } + + if (isPeriodic) { + /*--- Correct the sensor values across any periodic boundaries. ---*/ + + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); + } + + /*--- Set final pressure switch for each point ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) + nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); + } + + /*--- MPI parallelization ---*/ + + InitiateComms(geometry, config, SENSOR); + CompleteComms(geometry, config, SENSOR); + + } + + /*! + * \brief Generic implementation of explicit iterations with a preconditioner. + * \note The preconditioner is a functor implementing the methods: + * - compute(config, iPoint): Should prepare the preconditioner for iPoint. + * - apply(iVar, residual[], resTruncError[]): Apply it to compute the iVar update. + * See Explicit_Iteration for the general form of the preconditioner. + */ + template + void Explicit_Iteration_impl(ResidualPrecond& preconditioner, CGeometry *geometry, + CSolver **solver_container, CConfig *config, unsigned short iRKStep) { + + static_assert(IntegrationType == CLASSICAL_RK4_EXPLICIT || + IntegrationType == RUNGE_KUTTA_EXPLICIT || + IntegrationType == EULER_EXPLICIT, ""); + + const bool adjoint = config->GetContinuous_Adjoint(); + + const su2double RK_AlphaCoeff = config->Get_Alpha_RKStep(iRKStep); + + /*--- Hard-coded classical RK4 coefficients. Will be added to config. ---*/ + const su2double RK_FuncCoeff[] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0}; + const su2double RK_TimeCoeff[] = {0.5, 0.5, 1.0, 1.0}; + + /*--- Set shared residual variables to 0 and declare + * local ones for current thread to work on. ---*/ + + SU2_OMP_MASTER + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + SetRes_RMS(iVar, 0.0); + SetRes_Max(iVar, 0.0, 0); + } + SU2_OMP_BARRIER + + su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; + const su2double* coordMax[MAXNVAR] = {nullptr}; + unsigned long idxMax[MAXNVAR] = {0}; + + /*--- Update the solution and residuals ---*/ + + if (!adjoint) { + SU2_OMP(for schedule(static,omp_chunk_size) nowait) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); + su2double Delta = nodes->GetDelta_Time(iPoint) / Vol; + + const su2double* Res_TruncError = nodes->GetResTruncError(iPoint); + const su2double* Residual = LinSysRes.GetBlock(iPoint); + + preconditioner.compute(config, iPoint); + + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + + su2double Res = preconditioner.apply(iVar, Residual, Res_TruncError); + + /*--- "Static" switch which should be optimized at compile time. ---*/ + switch(IntegrationType) { + + case EULER_EXPLICIT: + nodes->AddSolution(iPoint,iVar, -Res*Delta); + break; + + case RUNGE_KUTTA_EXPLICIT: + nodes->AddSolution(iPoint, iVar, -Res*Delta*RK_AlphaCoeff); + break; + + case CLASSICAL_RK4_EXPLICIT: + { + su2double tmp_time = -1.0*RK_TimeCoeff[iRKStep]*Delta; + su2double tmp_func = -1.0*RK_FuncCoeff[iRKStep]*Delta; + + if (iRKStep < 3) { + /* Base Solution Update */ + nodes->AddSolution(iPoint,iVar, tmp_time*Res); + + /* New Solution Update */ + nodes->AddSolution_New(iPoint,iVar, tmp_func*Res); + } else { + nodes->SetSolution(iPoint, iVar, nodes->GetSolution_New(iPoint, iVar) + tmp_func*Res); + } + } + break; + } + + /*--- Update residual information for current thread. ---*/ + resRMS[iVar] += Res*Res; + if (fabs(Res) > resMax[iVar]) { + resMax[iVar] = fabs(Res); + idxMax[iVar] = iPoint; + coordMax[iVar] = geometry->nodes->GetCoord(iPoint); + } + } + } + /*--- Reduce residual information over all threads in this rank. ---*/ + SU2_OMP_CRITICAL + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + AddRes_RMS(iVar, resRMS[iVar]); + AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); + } + SU2_OMP_BARRIER + } + + /*--- MPI solution ---*/ + + InitiateComms(geometry, config, SOLUTION); + CompleteComms(geometry, config, SOLUTION); + + if (!adjoint) { + SU2_OMP_MASTER { + /*--- Compute the root mean square residual ---*/ + + SetResidual_RMS(geometry, config); + + /*--- For verification cases, compute the global error metrics. ---*/ + + ComputeVerificationError(geometry, config); + } + SU2_OMP_BARRIER + } + + } + + /*! + * \brief Generic implementation of explicit iterations without preconditioner. + */ + template + FORCEINLINE void Explicit_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iRKStep) { + struct Identity { + FORCEINLINE void compute(const CConfig*, unsigned long) {} + FORCEINLINE su2double apply(unsigned short iVar, const su2double* res, const su2double* resTrunc) { + return res[iVar] + resTrunc[iVar]; + } + } precond; + + Explicit_Iteration_impl(precond, geometry, solver_container, config, iRKStep); + } + + /*! + * \brief Generic implementation of implicit Euler iteration with an optional preconditioner applied to the diagonal. + * \param[in] compute_ur - Whether to use automatic under-relaxation for the update. + * \tparam DiagonalPrecond - A function object implementing: + * - active: A boolean variable to determine if the preconditioner should be used. + * - (config, iPoint, delta): Compute and return a matrix type compatible with the Jacobian matrix, + * where "delta" is V/dt. + */ + template + void ImplicitEuler_Iteration_impl(DiagonalPrecond& preconditioner, CGeometry *geometry, + CSolver **solver_container, CConfig *config, bool compute_ur) { + + const bool adjoint = config->GetContinuous_Adjoint(); + + /*--- Set shared residual variables to 0 and declare + * local ones for current thread to work on. ---*/ + + SU2_OMP_MASTER + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + SetRes_RMS(iVar, 0.0); + SetRes_Max(iVar, 0.0, 0); + } + SU2_OMP_BARRIER + + su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; + const su2double* coordMax[MAXNVAR] = {nullptr}; + unsigned long idxMax[MAXNVAR] = {0}; + + /*--- Build implicit system ---*/ + + SU2_OMP(for schedule(static,omp_chunk_size) nowait) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Read the residual ---*/ + + su2double* local_Res_TruncError = nodes->GetResTruncError(iPoint); + + /*--- Read the volume ---*/ + + su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); + + /*--- Modify matrix diagonal to assure diagonal dominance ---*/ + + if (nodes->GetDelta_Time(iPoint) != 0.0) { + + su2double Delta = Vol / nodes->GetDelta_Time(iPoint); + + if (preconditioner.active) { + Jacobian.AddBlock2Diag(iPoint, preconditioner(config, iPoint, Delta)); + } + else { + Jacobian.AddVal2Diag(iPoint, Delta); + } + } + else { + Jacobian.SetVal2Diag(iPoint, 1.0); + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + LinSysRes(iPoint,iVar) = 0.0; + local_Res_TruncError[iVar] = 0.0; + } + } + + /*--- Right hand side of the system (-Residual) and initial guess (x = 0) ---*/ + + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + unsigned long total_index = iPoint*nVar + iVar; + LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]); + LinSysSol[total_index] = 0.0; + + su2double Res = fabs(LinSysRes[total_index]); + resRMS[iVar] += Res*Res; + if (Res > resMax[iVar]) { + resMax[iVar] = Res; + idxMax[iVar] = iPoint; + coordMax[iVar] = geometry->nodes->GetCoord(iPoint); + } + } + } + SU2_OMP_CRITICAL + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + AddRes_RMS(iVar, resRMS[iVar]); + AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); + } + + /*--- Initialize residual and solution at the ghost points ---*/ + + SU2_OMP(sections nowait) + { + SU2_OMP(section) + for (unsigned long iPoint = nPointDomain; iPoint < nPoint; iPoint++) + LinSysRes.SetBlock_Zero(iPoint); + + SU2_OMP(section) + for (unsigned long iPoint = nPointDomain; iPoint < nPoint; iPoint++) + LinSysSol.SetBlock_Zero(iPoint); + } + + /*--- Solve or smooth the linear system. ---*/ + + auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); + SU2_OMP_MASTER + { + SetIterLinSolver(iter); + SetResLinSolver(System.GetResidual()); + } + SU2_OMP_BARRIER + + if (compute_ur) ComputeUnderRelaxationFactor(solver_container, config); + + /*--- Update solution (system written in terms of increments) ---*/ + + if (!adjoint) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + nodes->AddSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint)*LinSysSol[iPoint*nVar+iVar]); + } + } + } + + for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { + InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); + } + + /*--- MPI solution ---*/ + + InitiateComms(geometry, config, SOLUTION); + CompleteComms(geometry, config, SOLUTION); + + SU2_OMP_MASTER + { + /*--- Compute the root mean square residual ---*/ + + SetResidual_RMS(geometry, config); + + /*--- For verification cases, compute the global error metrics. ---*/ + + ComputeVerificationError(geometry, config); + } + SU2_OMP_BARRIER + + } + + /*! + * \brief Evaluate the vorticity and strain rate magnitude. + */ + inline void ComputeVorticityAndStrainMag(const CConfig& config, unsigned short iMesh) { + + SU2_OMP_MASTER { + StrainMag_Max = 0.0; + Omega_Max = 0.0; + } + SU2_OMP_BARRIER + + nodes->SetVorticity_StrainMag(); + + /*--- Min and Max are not really differentiable ---*/ + const bool wasActive = AD::BeginPassive(); + + su2double strainMax = 0.0, omegaMax = 0.0; + + SU2_OMP(for schedule(static,omp_chunk_size) nowait) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + strainMax = max(strainMax, nodes->GetStrainMag(iPoint)); + omegaMax = max(omegaMax, GeometryToolbox::Norm(3, nodes->GetVorticity(iPoint))); + } + SU2_OMP_CRITICAL { + StrainMag_Max = max(StrainMag_Max, strainMax); + Omega_Max = max(Omega_Max, omegaMax); + } + + if ((iMesh == MESH_0) && (config.GetComm_Level() == COMM_FULL)) { + SU2_OMP_BARRIER + SU2_OMP_MASTER + { + su2double MyOmega_Max = Omega_Max; + su2double MyStrainMag_Max = StrainMag_Max; + + SU2_MPI::Allreduce(&MyStrainMag_Max, &StrainMag_Max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + SU2_MPI::Allreduce(&MyOmega_Max, &Omega_Max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + } + SU2_OMP_BARRIER + } + + AD::EndPassive(wasActive); + } + /*! * \brief Destructor. */ @@ -282,6 +1029,18 @@ class CFVMFlowSolverBase : public CSolver { */ void ComputeUnderRelaxationFactor(CSolver** solver, const CConfig* config) final; + /*! + * \brief Set the total residual adding the term that comes from the Dual Time Strategy. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + */ + void SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iRKStep, + unsigned short iMesh, unsigned short RunTime_EqSystem) override; + /*! * \brief Set a uniform inlet profile * diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 8cfcdcd299b3..d2ccebcb9676 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1341,6 +1341,183 @@ void CFVMFlowSolverBase::SumEdgeFluxes(const CGeometry* geometry) { } } +template +void CFVMFlowSolverBase::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iRKStep, unsigned short iMesh, + unsigned short RunTime_EqSystem) { + /*--- Local variables ---*/ + + unsigned short iVar, iMarker, iDim, iNeigh; + unsigned long iPoint, jPoint, iEdge, iVertex; + + const su2double *U_time_nM1 = nullptr, *U_time_n = nullptr, *U_time_nP1 = nullptr; + su2double Volume_nM1, Volume_nP1, TimeStep; + const su2double *Normal = nullptr, *GridVel_i = nullptr, *GridVel_j = nullptr; + su2double Residual_GCL; + + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool first_order = (config->GetTime_Marching() == DT_STEPPING_1ST); + const bool second_order = (config->GetTime_Marching() == DT_STEPPING_2ND); + + /*--- Store the physical time step ---*/ + + TimeStep = config->GetDelta_UnstTimeND(); + + /*--- Compute the dual time-stepping source term for static meshes ---*/ + + if (!dynamic_grid) { + + /*--- Loop over all nodes (excluding halos) ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that + we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, + previous solutions that are stored in memory. ---*/ + + U_time_nM1 = nodes->GetSolution_time_n1(iPoint); + U_time_n = nodes->GetSolution_time_n(iPoint); + U_time_nP1 = nodes->GetSolution(iPoint); + + /*--- CV volume at time n+1. As we are on a static mesh, the volume + of the CV will remained fixed for all time steps. ---*/ + + Volume_nP1 = geometry->nodes->GetVolume(iPoint); + + /*--- Compute the dual time-stepping source term based on the chosen + time discretization scheme (1st- or 2nd-order).---*/ + + for (iVar = 0; iVar < nVar; iVar++) { + if (first_order) + LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*Volume_nP1 / TimeStep; + if (second_order) + LinSysRes(iPoint,iVar) += ( 3.0*U_time_nP1[iVar] - 4.0*U_time_n[iVar] + +1.0*U_time_nM1[iVar])*Volume_nP1 / (2.0*TimeStep); + } + + /*--- 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)); + } + } + + } + + else { + + /*--- For unsteady flows on dynamic meshes (rigidly transforming or + dynamically deforming), the Geometric Conservation Law (GCL) should be + satisfied in conjunction with the ALE formulation of the governing + equations. The GCL prevents accuracy issues caused by grid motion, i.e. + a uniform free-stream should be preserved through a moving grid. First, + we will loop over the edges and boundaries to compute the GCL component + of the dual time source term that depends on grid velocities. ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; ++iPoint) { + + GridVel_i = geometry->nodes->GetGridVel(iPoint); + U_time_n = nodes->GetSolution_time_n(iPoint); + + for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { + + iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); + Normal = geometry->edges->GetNormal(iEdge); + + jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + GridVel_j = geometry->nodes->GetGridVel(jPoint); + + /*--- Determine whether to consider the normal outward or inward. ---*/ + su2double dir = (geometry->edges->GetNode(iEdge,0) == iPoint)? 0.5 : -0.5; + + Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL += dir*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; + + for (iVar = 0; iVar < nVar; iVar++) + LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; + } + } + + /*--- Loop over the boundary edges ---*/ + + for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && + (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + + /*--- Get the index for node i plus the boundary face normal ---*/ + + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + + /*--- Grid velocities stored at boundary node i ---*/ + + GridVel_i = geometry->nodes->GetGridVel(iPoint); + + /*--- Compute the GCL term by dotting the grid velocity with the face + normal. The normal is negated to match the boundary convention. ---*/ + + Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; + + /*--- Compute the GCL component of the source term for node i ---*/ + + U_time_n = nodes->GetSolution_time_n(iPoint); + for (iVar = 0; iVar < nVar; iVar++) + LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; + } + } + } + + /*--- Loop over all nodes (excluding halos) to compute the remainder + of the dual time-stepping source term. ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; iPoint++) { + + /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that + we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, + previous solutions that are stored in memory. ---*/ + + U_time_nM1 = nodes->GetSolution_time_n1(iPoint); + U_time_n = nodes->GetSolution_time_n(iPoint); + U_time_nP1 = nodes->GetSolution(iPoint); + + /*--- 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 + volumes will remain constant. ---*/ + + Volume_nM1 = geometry->nodes->GetVolume_nM1(iPoint); + Volume_nP1 = geometry->nodes->GetVolume(iPoint); + + /*--- Compute the dual time-stepping source residual. Due to the + introduction of the GCL term above, the remainder of the source residual + due to the time discretization has a new form.---*/ + + for (iVar = 0; iVar < nVar; iVar++) { + if (first_order) + LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); + if (second_order) + LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) + + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); + } + + /*--- 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)); + } + } + } + +} + template void CFVMFlowSolverBase::Pressure_Forces(const CGeometry* geometry, const CConfig* config) { unsigned long iVertex, iPoint; diff --git a/SU2_CFD/include/solvers/CIncEulerSolver.hpp b/SU2_CFD/include/solvers/CIncEulerSolver.hpp index 85cd5b7805dc..9c1ab37d16d1 100644 --- a/SU2_CFD/include/solvers/CIncEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CIncEulerSolver.hpp @@ -38,13 +38,92 @@ */ class CIncEulerSolver : public CFVMFlowSolverBase { protected: - su2double - *Primitive = nullptr, /*!< \brief Auxiliary nPrimVar vector. */ - *Primitive_i = nullptr, /*!< \brief Auxiliary nPrimVar vector for storing the primitive at point i. */ - *Primitive_j = nullptr; /*!< \brief Auxiliary nPrimVar vector for storing the primitive at point j. */ - CFluidModel *FluidModel = nullptr; /*!< \brief fluid model used in the solver */ - su2double **Preconditioner = nullptr; /*!< \brief Auxiliary matrix for storing the low speed preconditioner. */ + + /*! + * \brief Preprocessing actions common to the Euler and NS solvers. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] Output - boolean to determine whether to print output. + */ + void CommonPreprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output); + + /*! + * \brief Compute the preconditioner for low-Mach flows. + * \param[in] iPoint - Index of the grid point + * \param[in] config - Definition of the particular problem. + * \param[in] delta - Volume over delta time, does not matter for explicit. + * \param[out] preconditioner - The preconditioner matrix. + */ + void SetPreconditioner(const CConfig *config, unsigned long iPoint, + su2double delta, su2activematrix& preconditioner) const; + + /*! + * \brief Compute a pressure sensor switch. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + */ + void SetCentered_Dissipation_Sensor(CGeometry *geometry, const CConfig *config); + + /*! + * \brief Compute the undivided laplacian for the solution, except the energy equation. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config); + + /*! + * \brief Compute the max eigenvalue. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetMax_Eigenvalue(CGeometry *geometry, const CConfig *config); + + /*! + * \brief Compute the velocity^2, SoundSpeed, Pressure, Enthalpy, Viscosity. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \return - The number of non-physical points. + */ + virtual unsigned long SetPrimitive_Variables(CSolver **solver_container, const CConfig *config); + + /*! + * \brief Update the Beta parameter for the incompressible preconditioner. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - current mesh level for the multigrid. + */ + void SetBeta_Parameter(CGeometry *geometry, + CSolver **solver_container, + CConfig *config, + unsigned short iMesh); + + /*! + * \brief A virtual member. + */ + void GetOutlet_Properties(CGeometry *geometry, + CConfig *config, + unsigned short iMesh, + bool Output); + + /*! + * \brief Set the solver nondimensionalization. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + */ + void SetNondimensionalization(CConfig *config, unsigned short iMesh); + + /*! + * \brief Generic implementation of explicit iterations with preconditioner. + */ + template + void Explicit_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iRKStep); public: /*! @@ -66,13 +145,6 @@ class CIncEulerSolver : public CFVMFlowSolverBase a,std::vector b); - /*! * \brief Update the solution using a Runge-Kutta scheme. * \param[in] geometry - Geometrical definition of the problem. @@ -289,6 +309,18 @@ class CIncEulerSolver : public CFVMFlowSolverBaseTRUE means that it is an adjoint solver. @@ -1667,13 +1660,6 @@ class CSolver { const CConfig *config, bool reconstruction = false) { } - /*! - * \brief A virtual member. - * \param[in] iPoint - Index of the grid point. - * \param[in] config - Definition of the particular problem. - */ - inline virtual void SetPreconditioner(CConfig *config, unsigned long iPoint) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. @@ -2235,18 +2221,6 @@ class CSolver { */ inline virtual su2double GetInflow_MassFlow(unsigned short val_marker) const { return 0; } - /*! - * \brief A virtual member. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] iMesh - current mesh level for the multigrid. - * \param[in] Output - boolean to determine whether to print output. - */ - inline virtual void GetOutlet_Properties(CGeometry *geometry, - CConfig *config, - unsigned short iMesh, - bool Output) { } - /*! * \brief A virtual member. * \param[in] config - Definition of the particular problem. @@ -3925,7 +3899,7 @@ class CSolver { * \brief A virtual member. * \param[in] config - Definition of the particular problem. */ - inline virtual void SetFreeStream_Solution(CConfig *config) { } + inline virtual void SetFreeStream_Solution(const CConfig *config) { } /*! * \brief A virtual member. @@ -4314,18 +4288,6 @@ class CSolver { */ inline virtual void SetFreeStream_TurboSolution(CConfig *config) { } - /*! - * \brief A virtual member. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] config - Definition of the particular problem. - * \param[in] iMesh - current mesh level for the multigrid. - */ - inline virtual void SetBeta_Parameter(CGeometry *geometry, - CSolver **solver_container, - CConfig *config, - unsigned short iMesh) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition. diff --git a/SU2_CFD/include/solvers/CTurbSASolver.hpp b/SU2_CFD/include/solvers/CTurbSASolver.hpp index 23c819a5e62f..31dfc46904f1 100644 --- a/SU2_CFD/include/solvers/CTurbSASolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSASolver.hpp @@ -353,7 +353,7 @@ class CTurbSASolver final : public CTurbSolver { * \brief Set the solution using the Freestream values. * \param[in] config - Definition of the particular problem. */ - inline void SetFreeStream_Solution(CConfig *config) override { + inline void SetFreeStream_Solution(const CConfig *config) override { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) nodes->SetSolution(iPoint, 0, nu_tilde_Inf); } diff --git a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp index a1f7df67848d..205194444cfa 100644 --- a/SU2_CFD/include/solvers/CTurbSSTSolver.hpp +++ b/SU2_CFD/include/solvers/CTurbSSTSolver.hpp @@ -239,7 +239,7 @@ class CTurbSSTSolver final : public CTurbSolver { * \brief Set the solution using the Freestream values. * \param[in] config - Definition of the particular problem. */ - inline void SetFreeStream_Solution(CConfig *config) override { + inline void SetFreeStream_Solution(const CConfig *config) override { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ nodes->SetSolution(iPoint, 0, kine_Inf); nodes->SetSolution(iPoint, 1, omega_Inf); diff --git a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp index 9e29f47ed552..dea4fbbbb67f 100644 --- a/SU2_CFD/include/variables/CNEMOEulerVariable.hpp +++ b/SU2_CFD/include/variables/CNEMOEulerVariable.hpp @@ -141,7 +141,7 @@ class CNEMOEulerVariable : public CVariable { * \param[in] iVar - Index of the variable. * \return Set the value of the primitive variable for the index iVar. */ - inline void SetPrimitive(unsigned long iPoint, unsigned long iVar, su2double val_prim) override { Primitive(iPoint,iVar) = val_prim; } + inline void SetPrimitive(unsigned long iPoint, unsigned long iVar, su2double val_prim) final { Primitive(iPoint,iVar) = val_prim; } /*! * \brief Set the value of the primitive variables. @@ -283,7 +283,7 @@ class CNEMOEulerVariable : public CVariable { * \param[in] iDim - Index of the dimension. * \param[in] value - Value of the reconstruction gradient component. */ - inline void SetGradient_Reconstruction(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) override { + inline void SetGradient_Reconstruction(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) final { Gradient_Reconstruction(iPoint,iVar,iDim) = value; } @@ -364,7 +364,7 @@ class CNEMOEulerVariable : public CVariable { * \brief Set the norm 2 of the velocity. * \return Norm 2 of the velocity vector. */ - void SetVelocity2(unsigned long iPoint) override; + void SetVelocity2(unsigned long iPoint) final; /*! * \brief Get the norm 2 of the velocity. @@ -459,28 +459,28 @@ class CNEMOEulerVariable : public CVariable { * \brief A virtual member. * \return Value of the vibrational-electronic temperature. */ - inline su2double GetTemperature_ve(unsigned long iPoint) const override + inline su2double GetTemperature_ve(unsigned long iPoint) const final { return Primitive(iPoint,TVE_INDEX); } /*! * \brief Sets the vibrational electronic temperature of the flow. * \return Value of the temperature of the flow. */ - inline bool SetTemperature_ve(unsigned long iPoint, su2double val_Tve) override + inline bool SetTemperature_ve(unsigned long iPoint, su2double val_Tve) final { Primitive(iPoint,TVE_INDEX) = val_Tve; return false; } /*! * \brief Get the mixture specific heat at constant volume (trans.-rot.). * \return \f$\rho C^{t-r}_{v} \f$ */ - inline su2double GetRhoCv_tr(unsigned long iPoint) const override + inline su2double GetRhoCv_tr(unsigned long iPoint) const final { return Primitive(iPoint,RHOCVTR_INDEX); } /*! * \brief Get the mixture specific heat at constant volume (vib.-el.). * \return \f$\rho C^{v-e}_{v} \f$ */ - inline su2double GetRhoCv_ve(unsigned long iPoint) const override + inline su2double GetRhoCv_ve(unsigned long iPoint) const final { return Primitive(iPoint,RHOCVVE_INDEX); } /*! @@ -496,24 +496,24 @@ class CNEMOEulerVariable : public CVariable { /*! * \brief Set partial derivative of pressure w.r.t. density \f$\frac{\partial P}{\partial \rho_s}\f$ */ - inline su2double *GetdPdU(unsigned long iPoint) override { return dPdU[iPoint]; } + inline su2double *GetdPdU(unsigned long iPoint) final { return dPdU[iPoint]; } /*! * \brief Set partial derivative of temperature w.r.t. density \f$\frac{\partial T}{\partial \rho_s}\f$ */ - inline su2double *GetdTdU(unsigned long iPoint) override { return dTdU[iPoint]; } + inline su2double *GetdTdU(unsigned long iPoint) final { return dTdU[iPoint]; } /*! * \brief Set partial derivative of vib.-el. temperature w.r.t. density \f$\frac{\partial T^{V-E}}{\partial \rho_s}\f$ */ - inline su2double *GetdTvedU(unsigned long iPoint) override { return dTvedU[iPoint]; } + inline su2double *GetdTvedU(unsigned long iPoint) final { return dTvedU[iPoint]; } /*! * \brief Get the mass fraction \f$\rho_s / \rho \f$ of species s. * \param[in] val_Species - Index of species s. * \return Value of the mass fraction of species s. */ - inline su2double GetMassFraction(unsigned long iPoint, unsigned long val_Species) const override { + inline su2double GetMassFraction(unsigned long iPoint, unsigned long val_Species) const final { return Primitive(iPoint,RHOS_INDEX+val_Species) / Primitive(iPoint,RHO_INDEX); } diff --git a/SU2_CFD/include/variables/CNEMONSVariable.hpp b/SU2_CFD/include/variables/CNEMONSVariable.hpp index 5531409ce8e0..657e5aac85ae 100644 --- a/SU2_CFD/include/variables/CNEMONSVariable.hpp +++ b/SU2_CFD/include/variables/CNEMONSVariable.hpp @@ -104,17 +104,6 @@ class CNEMONSVariable final : public CNEMOEulerVariable { */ inline const MatrixType& GetPrimitive_Aux(void) const { return Primitive_Aux; } - /*! - * \brief Set the value of the reconstruction variables gradient at a node. - * \param[in] iPoint - Index of the current node. - * \param[in] iVar - Index of the variable. - * \param[in] iDim - Index of the dimension. - * \param[in] value - Value of the reconstruction gradient component. - */ - /* Works as a dummy function for consistency since no reconstruction is needed for primitive variables*/ - inline void SetGradient_Reconstruction(unsigned long iPoint, unsigned long iVar, unsigned long iDim, su2double value) override { } - - /*! * \brief Set all the primitive variables for compressible flows. */ diff --git a/SU2_CFD/src/SU2_CFD.cpp b/SU2_CFD/src/SU2_CFD.cpp index 8ad8edc08fc2..25246cc1d21d 100644 --- a/SU2_CFD/src/SU2_CFD.cpp +++ b/SU2_CFD/src/SU2_CFD.cpp @@ -60,7 +60,7 @@ int main(int argc, char *argv[]) { /*--- MPI initialization, and buffer setting ---*/ -#ifdef HAVE_OMP +#if defined(HAVE_OMP) && defined(HAVE_MPI) int required = use_thread_mult? MPI_THREAD_MULTIPLE : MPI_THREAD_FUNNELED; int provided; SU2_MPI::Init_thread(&argc, &argv, required, &provided); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 83221d895ccc..99342009f6c4 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2200,7 +2200,7 @@ void CEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_con SU2_OMP_BARRIER SU2_OMP_ATOMIC - ErrorCounter += SetPrimitive_Variables(solver_container, config, Output); + ErrorCounter += SetPrimitive_Variables(solver_container, config); if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) { SU2_OMP_BARRIER @@ -2305,7 +2305,7 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container } } -unsigned long CEulerSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { +unsigned long CEulerSolver::SetPrimitive_Variables(CSolver **solver_container, const CConfig *config) { /*--- Number of non-physical points, local to the thread, needs * further reduction if function is called in parallel ---*/ @@ -2330,260 +2330,50 @@ unsigned long CEulerSolver::SetPrimitive_Variables(CSolver **solver_container, C void CEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned long Iteration) { - const bool viscous = config->GetViscous(); - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool time_stepping = (config->GetTime_Marching() == TIME_STEPPING); - const bool dual_time = (config->GetTime_Marching() == DT_STEPPING_1ST) || - (config->GetTime_Marching() == DT_STEPPING_2ND); - const su2double K_v = 0.25; - - /*--- Init thread-shared variables to compute min/max values. - * Critical sections are used for this instead of reduction - * clauses for compatibility with OpenMP 2.0 (Windows...). ---*/ - - SU2_OMP_MASTER - { - Min_Delta_Time = 1e30; - Max_Delta_Time = 0.0; - Global_Delta_UnstTimeND = 1e30; - } - SU2_OMP_BARRIER - - const su2double *Normal = nullptr; - su2double Area, Vol, Mean_SoundSpeed, Mean_ProjVel, Lambda, Local_Delta_Time, Local_Delta_Time_Visc; - su2double Mean_LaminarVisc, Mean_EddyVisc, Mean_Density, Lambda_1, Lambda_2; - unsigned long iEdge, iVertex, iPoint, jPoint; - unsigned short iDim, iMarker; - - /*--- Loop domain points. ---*/ - - SU2_OMP_FOR_DYN(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - /*--- Set maximum eigenvalues to zero. ---*/ - - nodes->SetMax_Lambda_Inv(iPoint,0.0); - - if (viscous) - nodes->SetMax_Lambda_Visc(iPoint,0.0); - - /*--- Loop over the neighbors of point i. ---*/ - - for (unsigned short iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) - { - jPoint = geometry->nodes->GetPoint(iPoint,iNeigh); - - iEdge = geometry->nodes->GetEdge(iPoint,iNeigh); - Normal = geometry->edges->GetNormal(iEdge); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); - Mean_SoundSpeed = 0.5 * (nodes->GetSoundSpeed(iPoint) + nodes->GetSoundSpeed(jPoint)) * Area; - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - const su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); - const su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); - - for (iDim = 0; iDim < nDim; iDim++) - Mean_ProjVel -= 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * Normal[iDim]; - } - - /*--- Inviscid contribution ---*/ - - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed ; - nodes->AddMax_Lambda_Inv(iPoint,Lambda); - - /*--- Viscous contribution ---*/ - - if (!viscous) continue; - - Mean_LaminarVisc = 0.5*(nodes->GetLaminarViscosity(iPoint) + nodes->GetLaminarViscosity(jPoint)); - Mean_EddyVisc = 0.5*(nodes->GetEddyViscosity(iPoint) + nodes->GetEddyViscosity(jPoint)); - Mean_Density = 0.5*(nodes->GetDensity(iPoint) + nodes->GetDensity(jPoint)); - - Lambda_1 = (4.0/3.0)*(Mean_LaminarVisc + Mean_EddyVisc); - //TODO (REAL_GAS) removing Gamma it cannot work with FLUIDPROP - Lambda_2 = (1.0 + (Prandtl_Lam/Prandtl_Turb)*(Mean_EddyVisc/Mean_LaminarVisc))*(Gamma*Mean_LaminarVisc/Prandtl_Lam); - - Lambda = (Lambda_1 + Lambda_2)*Area*Area/Mean_Density; - nodes->AddMax_Lambda_Visc(iPoint, Lambda); + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return 0.5 * (nodes.GetSoundSpeed(iPoint) + nodes.GetSoundSpeed(jPoint)); } - } - - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - - if (!geometry->nodes->GetDomain(iPoint)) continue; - - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - Mean_SoundSpeed = nodes->GetSoundSpeed(iPoint) * Area; - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - const su2double *GridVel = geometry->nodes->GetGridVel(iPoint); - - for (iDim = 0; iDim < nDim; iDim++) - Mean_ProjVel -= GridVel[iDim]*Normal[iDim]; - } - - /*--- Inviscid contribution ---*/ - - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - nodes->AddMax_Lambda_Inv(iPoint,Lambda); - - /*--- Viscous contribution ---*/ - - if (!viscous) continue; - - Mean_LaminarVisc = nodes->GetLaminarViscosity(iPoint); - Mean_EddyVisc = nodes->GetEddyViscosity(iPoint); - Mean_Density = nodes->GetDensity(iPoint); - - Lambda_1 = (4.0/3.0)*(Mean_LaminarVisc + Mean_EddyVisc); - Lambda_2 = (1.0 + (Prandtl_Lam/Prandtl_Turb)*(Mean_EddyVisc/Mean_LaminarVisc))*(Gamma*Mean_LaminarVisc/Prandtl_Lam); - Lambda = (Lambda_1 + Lambda_2)*Area*Area/Mean_Density; - - nodes->AddMax_Lambda_Visc(iPoint, Lambda); - - } + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint) const { + return nodes.GetSoundSpeed(iPoint); } - } - - /*--- Each element uses their own speed, steady state simulation. ---*/ - { - /*--- Thread-local variables for min/max reduction. ---*/ - su2double minDt = 1e30, maxDt = 0.0; - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - Vol = geometry->nodes->GetVolume(iPoint); + } soundSpeed; - if (Vol != 0.0) { - Local_Delta_Time = nodes->GetLocalCFL(iPoint)*Vol / nodes->GetMax_Lambda_Inv(iPoint); + /*--- Define an object to compute the viscous eigenvalue. ---*/ + struct LambdaVisc { + const su2double gamma, prandtlLam, prandtlTurb; - if(viscous) { - Local_Delta_Time_Visc = nodes->GetLocalCFL(iPoint)*K_v*Vol*Vol/ nodes->GetMax_Lambda_Visc(iPoint); - Local_Delta_Time = min(Local_Delta_Time, Local_Delta_Time_Visc); - } - - minDt = min(minDt, Local_Delta_Time); - maxDt = max(maxDt, Local_Delta_Time); + LambdaVisc(su2double g, su2double pl, su2double pt) : gamma(g), prandtlLam(pl), prandtlTurb(pt) {} - nodes->SetDelta_Time(iPoint, min(Local_Delta_Time, config->GetMax_DeltaTime())); - } - else { - nodes->SetDelta_Time(iPoint,0.0); - } - } - /*--- Min/max over threads. ---*/ - SU2_OMP_CRITICAL - { - Min_Delta_Time = min(Min_Delta_Time, minDt); - Max_Delta_Time = max(Max_Delta_Time, maxDt); - Global_Delta_Time = Min_Delta_Time; + FORCEINLINE su2double lambda(su2double laminarVisc, su2double eddyVisc, su2double density) const { + su2double Lambda_1 = (4.0/3.0)*(laminarVisc + eddyVisc); + /// TODO: (REAL_GAS) removing gamma as it cannot work with FLUIDPROP + su2double Lambda_2 = (1.0 + (prandtlLam/prandtlTurb)*(eddyVisc/laminarVisc))*(gamma*laminarVisc/prandtlLam); + return (Lambda_1 + Lambda_2) / density; } - SU2_OMP_BARRIER - } - - /*--- Compute the min/max dt (in parallel, now over mpi ranks). ---*/ - - SU2_OMP_MASTER - if (config->GetComm_Level() == COMM_FULL) { - su2double rbuf_time; - SU2_MPI::Allreduce(&Min_Delta_Time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - Min_Delta_Time = rbuf_time; - - SU2_MPI::Allreduce(&Max_Delta_Time, &rbuf_time, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - Max_Delta_Time = rbuf_time; - } - SU2_OMP_BARRIER - - /*--- For exact time solution use the minimum delta time of the whole mesh. ---*/ - if (time_stepping) { - - /*--- If the unsteady CFL is set to zero, it uses the defined unsteady time step, - * otherwise it computes the time step based on the unsteady CFL. ---*/ - - SU2_OMP_MASTER - { - if (config->GetUnst_CFL() == 0.0) { - Global_Delta_Time = config->GetDelta_UnstTime(); - } - else { - Global_Delta_Time = Min_Delta_Time; - } - Max_Delta_Time = Global_Delta_Time; - config->SetDelta_UnstTimeND(Global_Delta_Time); + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + su2double laminarVisc = 0.5*(nodes.GetLaminarViscosity(iPoint) + nodes.GetLaminarViscosity(jPoint)); + su2double eddyVisc = 0.5*(nodes.GetEddyViscosity(iPoint) + nodes.GetEddyViscosity(jPoint)); + su2double density = 0.5*(nodes.GetDensity(iPoint) + nodes.GetDensity(jPoint)); + return lambda(laminarVisc, eddyVisc, density); } - SU2_OMP_BARRIER - - /*--- Sets the regular CFL equal to the unsteady CFL. ---*/ - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - nodes->SetLocalCFL(iPoint, config->GetUnst_CFL()); - nodes->SetDelta_Time(iPoint, Global_Delta_Time); - } - - } - - /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is diferent from 0. ---*/ - - if ((dual_time) && (Iteration == 0) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) { - - /*--- Thread-local variable for reduction. ---*/ - su2double glbDtND = 1e30; - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - glbDtND = min(glbDtND, config->GetUnst_CFL()*Global_Delta_Time / nodes->GetLocalCFL(iPoint)); + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint) const { + su2double laminarVisc = nodes.GetLaminarViscosity(iPoint); + su2double eddyVisc = nodes.GetEddyViscosity(iPoint); + su2double density = nodes.GetDensity(iPoint); + return lambda(laminarVisc, eddyVisc, density); } - SU2_OMP_CRITICAL - Global_Delta_UnstTimeND = min(Global_Delta_UnstTimeND, glbDtND); - SU2_OMP_BARRIER - - SU2_OMP_MASTER - { - SU2_MPI::Allreduce(&Global_Delta_UnstTimeND, &glbDtND, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - Global_Delta_UnstTimeND = glbDtND; - config->SetDelta_UnstTimeND(Global_Delta_UnstTimeND); - } - SU2_OMP_BARRIER - } + } lambdaVisc(Gamma, Prandtl_Lam, Prandtl_Turb); - /*--- The pseudo local time (explicit integration) cannot be greater than the physical time ---*/ + /*--- Now instantiate the generic implementation with the two functors above. ---*/ - if (dual_time && !implicit) { - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - Local_Delta_Time = min((2.0/3.0)*config->GetDelta_UnstTimeND(), nodes->GetDelta_Time(iPoint)); - nodes->SetDelta_Time(iPoint, Local_Delta_Time); - } - } + SetTime_Step_impl(soundSpeed, lambdaVisc, geometry, solver_container, config, iMesh, Iteration); } @@ -3166,97 +2956,23 @@ void CEulerSolver::Source_Template(CGeometry *geometry, CSolver **solver_contain } -void CEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, CConfig *config) { - - /*--- Loop domain points. ---*/ - - SU2_OMP_FOR_DYN(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - /*--- Set eigenvalues to zero. ---*/ - nodes->SetLambda(iPoint,0.0); - - /*--- Loop over the neighbors of point i. ---*/ - for (unsigned short iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) - { - auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); - - auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); - auto Normal = geometry->edges->GetNormal(iEdge); - su2double Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - su2double Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); - su2double Mean_SoundSpeed = 0.5 * (nodes->GetSoundSpeed(iPoint) + nodes->GetSoundSpeed(jPoint)) * Area; - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - const su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); - const su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); - - for (unsigned short iDim = 0; iDim < nDim; iDim++) - Mean_ProjVel -= 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * Normal[iDim]; - } - - /*--- Inviscid contribution ---*/ +void CEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, const CConfig *config) { - su2double Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed ; - nodes->AddLambda(iPoint, Lambda); + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return 0.5 * (nodes.GetSoundSpeed(iPoint) + nodes.GetSoundSpeed(jPoint)); } - } - - /*--- Loop boundary edges ---*/ - - for (unsigned short iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (unsigned long iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - - auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - su2double Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - su2double Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - su2double Mean_SoundSpeed = nodes->GetSoundSpeed(iPoint) * Area; - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - auto GridVel = geometry->nodes->GetGridVel(iPoint); - for (unsigned short iDim = 0; iDim < nDim; iDim++) - Mean_ProjVel -= GridVel[iDim]*Normal[iDim]; - } - - /*--- Inviscid contribution ---*/ - - su2double Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) { - nodes->AddLambda(iPoint,Lambda); - } - } + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint) const { + return nodes.GetSoundSpeed(iPoint); } - } - - /*--- Correct the eigenvalue values across any periodic boundaries. ---*/ - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_MAX_EIG); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_MAX_EIG); - } + } soundSpeed; - /*--- MPI parallelization ---*/ + /*--- Instantiate generic implementation. ---*/ - InitiateComms(geometry, config, MAX_EIGENVALUE); - CompleteComms(geometry, config, MAX_EIGENVALUE); + SetMax_Eigenvalue_impl(soundSpeed, geometry, config); } @@ -3308,62 +3024,15 @@ void CEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *co void CEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, const CConfig *config) { - /*--- We can access memory more efficiently if there are no periodic boundaries. ---*/ - - const bool isPeriodic = (config->GetnMarker_Periodic() > 0); - - /*--- Loop domain points. ---*/ - - SU2_OMP_FOR_DYN(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - const bool boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint); - const su2double Pressure_i = nodes->GetPressure(iPoint); - - /*--- Initialize. ---*/ - iPoint_UndLapl[iPoint] = 0.0; - jPoint_UndLapl[iPoint] = 0.0; - - /*--- Loop over the neighbors of point i. ---*/ - for (auto jPoint : geometry->nodes->GetPoints(iPoint)) - { - bool boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); - - /*--- If iPoint is boundary it only takes contributions from other boundary points. ---*/ - if (boundary_i && !boundary_j) continue; - - su2double Pressure_j = nodes->GetPressure(jPoint); - - /*--- Dissipation sensor, add pressure difference and pressure sum. ---*/ - iPoint_UndLapl[iPoint] += Pressure_j - Pressure_i; - jPoint_UndLapl[iPoint] += Pressure_j + Pressure_i; - } - - if (!isPeriodic) { - nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); + /*--- Define an object for the sensor variable, pressure. ---*/ + struct SensVar { + FORCEINLINE su2double operator() (const CEulerVariable& nodes, unsigned long iPoint) const { + return nodes.GetPressure(iPoint); } - } - - if (isPeriodic) { - /*--- Correct the sensor values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); - } - - /*--- Set final pressure switch for each point ---*/ - - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) - nodes->SetSensor(iPoint, fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); - } - - /*--- MPI parallelization ---*/ - - InitiateComms(geometry, config, SENSOR); - CompleteComms(geometry, config, SENSOR); + } sensVar; + /*--- Instantiate generic implementation. ---*/ + SetCentered_Dissipation_Sensor_impl(sensVar, geometry, config); } void CEulerSolver::SetUpwind_Ducros_Sensor(CGeometry *geometry, CConfig *config){ @@ -3416,120 +3085,6 @@ void CEulerSolver::SetUpwind_Ducros_Sensor(CGeometry *geometry, CConfig *config) } -template -void CEulerSolver::Explicit_Iteration(CGeometry *geometry, CSolver **solver_container, - CConfig *config, unsigned short iRKStep) { - - static_assert(IntegrationType == CLASSICAL_RK4_EXPLICIT || - IntegrationType == RUNGE_KUTTA_EXPLICIT || - IntegrationType == EULER_EXPLICIT, ""); - - const bool adjoint = config->GetContinuous_Adjoint(); - - const su2double RK_AlphaCoeff = config->Get_Alpha_RKStep(iRKStep); - - /*--- Hard-coded classical RK4 coefficients. Will be added to config. ---*/ - const su2double RK_FuncCoeff[] = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0}; - const su2double RK_TimeCoeff[] = {0.5, 0.5, 1.0, 1.0}; - - /*--- Set shared residual variables to 0 and declare - * local ones for current thread to work on. ---*/ - - SU2_OMP_MASTER - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - SU2_OMP_BARRIER - - su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; - unsigned long idxMax[MAXNVAR] = {0}; - - /*--- Update the solution and residuals ---*/ - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - - su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); - su2double Delta = nodes->GetDelta_Time(iPoint) / Vol; - - const su2double* Res_TruncError = nodes->GetResTruncError(iPoint); - const su2double* Residual = LinSysRes.GetBlock(iPoint); - - if (!adjoint) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - - su2double Res = Residual[iVar] + Res_TruncError[iVar]; - - /*--- "Static" switch which should be optimized at compile time. ---*/ - switch(IntegrationType) { - - case EULER_EXPLICIT: - nodes->AddSolution(iPoint,iVar, -Res*Delta); - break; - - case RUNGE_KUTTA_EXPLICIT: - nodes->AddSolution(iPoint, iVar, -Res*Delta*RK_AlphaCoeff); - break; - - case CLASSICAL_RK4_EXPLICIT: - { - su2double tmp_time = -1.0*RK_TimeCoeff[iRKStep]*Delta; - su2double tmp_func = -1.0*RK_FuncCoeff[iRKStep]*Delta; - - if (iRKStep < 3) { - /* Base Solution Update */ - nodes->AddSolution(iPoint,iVar, tmp_time*Res); - - /* New Solution Update */ - nodes->AddSolution_New(iPoint,iVar, tmp_func*Res); - } else { - nodes->SetSolution(iPoint, iVar, nodes->GetSolution_New(iPoint, iVar) + tmp_func*Res); - } - } - break; - } - - /*--- Update residual information for current thread. ---*/ - resRMS[iVar] += Res*Res; - if (fabs(Res) > resMax[iVar]) { - resMax[iVar] = fabs(Res); - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } - } - } - } - if (!adjoint) { - /*--- Reduce residual information over all threads in this rank. ---*/ - SU2_OMP_CRITICAL - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - AddRes_RMS(iVar, resRMS[iVar]); - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - } - SU2_OMP_BARRIER - - /*--- MPI solution ---*/ - - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - - SU2_OMP_MASTER - { - /*--- Compute the root mean square residual ---*/ - - SetResidual_RMS(geometry, config); - - /*--- For verification cases, compute the global error metrics. ---*/ - - ComputeVerificationError(geometry, config); - } - SU2_OMP_BARRIER - -} - void CEulerSolver::ExplicitRK_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iRKStep) { @@ -3549,159 +3104,28 @@ void CEulerSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **solver void CEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - const bool adjoint = config->GetContinuous_Adjoint(); - const bool roe_turkel = config->GetKind_Upwind_Flow() == TURKEL; - const bool low_mach_prec = config->Low_Mach_Preconditioning(); - - /*--- Local matrix for preconditioning. ---*/ - su2double** LowMachPrec = nullptr; - if (roe_turkel || low_mach_prec) { - LowMachPrec = new su2double* [nVar]; - for(unsigned short iVar = 0; iVar < nVar; ++iVar) - LowMachPrec[iVar] = new su2double [nVar]; - } - - /*--- Set shared residual variables to 0 and declare - * local ones for current thread to work on. ---*/ - - SU2_OMP_MASTER - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - SU2_OMP_BARRIER - - su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; - const su2double* coordMax[MAXNVAR] = {nullptr}; - unsigned long idxMax[MAXNVAR] = {0}; - - /*--- Build implicit system ---*/ - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + struct LowMachPrec { + const CEulerSolver* solver; + const bool active; + su2activematrix matrix; - /*--- Read the residual ---*/ - - su2double* local_Res_TruncError = nodes->GetResTruncError(iPoint); - - /*--- Read the volume ---*/ - - su2double Vol = geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint); - - /*--- Modify matrix diagonal to assure diagonal dominance ---*/ - - if (nodes->GetDelta_Time(iPoint) != 0.0) { - - su2double Delta = Vol / nodes->GetDelta_Time(iPoint); - - if (roe_turkel || low_mach_prec) { - SetPreconditioner(config, iPoint, Delta, LowMachPrec); - Jacobian.AddBlock2Diag(iPoint, LowMachPrec); - } - else { - Jacobian.AddVal2Diag(iPoint, Delta); - } - } - else { - Jacobian.SetVal2Diag(iPoint, 1.0); - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - LinSysRes(iPoint,iVar) = 0.0; - local_Res_TruncError[iVar] = 0.0; - } + LowMachPrec(const CEulerSolver* s, bool a, unsigned short nVar) : solver(s), active(a) { + if (active) matrix.resize(nVar,nVar); } - /*--- Right hand side of the system (-Residual) and initial guess (x = 0) ---*/ - - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - unsigned long total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]); - LinSysSol[total_index] = 0.0; - - su2double Res = fabs(LinSysRes[total_index]); - resRMS[iVar] += Res*Res; - if (Res > resMax[iVar]) { - resMax[iVar] = Res; - idxMax[iVar] = iPoint; - coordMax[iVar] = geometry->nodes->GetCoord(iPoint); - } + FORCEINLINE const su2activematrix& operator() (const CConfig* config, unsigned long iPoint, su2double delta) { + solver->SetPreconditioner(config, iPoint, delta, matrix); + return matrix; } - } - SU2_OMP_CRITICAL - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - AddRes_RMS(iVar, resRMS[iVar]); - AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); - } - - /*--- Initialize residual and solution at the ghost points ---*/ - - SU2_OMP(sections nowait) - { - SU2_OMP(section) - for (unsigned long iPoint = nPointDomain; iPoint < nPoint; iPoint++) - LinSysRes.SetBlock_Zero(iPoint); - - SU2_OMP(section) - for (unsigned long iPoint = nPointDomain; iPoint < nPoint; iPoint++) - LinSysSol.SetBlock_Zero(iPoint); - } - - /*--- Free local preconditioner. ---*/ - if (LowMachPrec) { - for(unsigned short iVar = 0; iVar < nVar; ++iVar) - delete [] LowMachPrec[iVar]; - delete [] LowMachPrec; - } - - /*--- Solve or smooth the linear system. ---*/ - - auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); - SU2_OMP_MASTER - { - SetIterLinSolver(iter); - SetResLinSolver(System.GetResidual()); - } - SU2_OMP_BARRIER - - ComputeUnderRelaxationFactor(solver_container, config); + } precond(this, config->Low_Mach_Preconditioning() || (config->GetKind_Upwind_Flow() == TURKEL), nVar); - /*--- Update solution (system written in terms of increments) ---*/ - - if (!adjoint) { - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - nodes->AddSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint)*LinSysSol[iPoint*nVar+iVar]); - } - } - } - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); - } - - /*--- MPI solution ---*/ - - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - - SU2_OMP_MASTER - { - /*--- Compute the root mean square residual ---*/ - - SetResidual_RMS(geometry, config); - - /*--- For verification cases, compute the global error metrics. ---*/ - - ComputeVerificationError(geometry, config); - } - SU2_OMP_BARRIER + ImplicitEuler_Iteration_impl(precond, geometry, solver_container, config, true); } void CEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPoint, - su2double delta, su2double** preconditioner) const { + su2double delta, su2activematrix& preconditioner) const { unsigned short iDim, jDim, iVar, jVar; su2double local_Mach, rho, enthalpy, soundspeed, sq_vel; @@ -9760,182 +9184,6 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver } } -void CEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config, - unsigned short iRKStep, unsigned short iMesh, unsigned short RunTime_EqSystem) { - - /*--- Local variables ---*/ - - unsigned short iVar, iMarker, iDim, iNeigh; - unsigned long iPoint, jPoint, iEdge, iVertex; - - const su2double *U_time_nM1 = nullptr, *U_time_n = nullptr, *U_time_nP1 = nullptr; - su2double Volume_nM1, Volume_nP1, TimeStep; - const su2double *Normal = nullptr, *GridVel_i = nullptr, *GridVel_j = nullptr; - su2double Residual_GCL; - - const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - const bool first_order = (config->GetTime_Marching() == DT_STEPPING_1ST); - const bool second_order = (config->GetTime_Marching() == DT_STEPPING_2ND); - - /*--- Store the physical time step ---*/ - - TimeStep = config->GetDelta_UnstTimeND(); - - /*--- Compute the dual time-stepping source term for static meshes ---*/ - - if (!dynamic_grid) { - - /*--- Loop over all nodes (excluding halos) ---*/ - - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that - we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, - previous solutions that are stored in memory. ---*/ - - U_time_nM1 = nodes->GetSolution_time_n1(iPoint); - U_time_n = nodes->GetSolution_time_n(iPoint); - U_time_nP1 = nodes->GetSolution(iPoint); - - /*--- CV volume at time n+1. As we are on a static mesh, the volume - of the CV will remained fixed for all time steps. ---*/ - - Volume_nP1 = geometry->nodes->GetVolume(iPoint); - - /*--- Compute the dual time-stepping source term based on the chosen - time discretization scheme (1st- or 2nd-order).---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - if (first_order) - LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*Volume_nP1 / TimeStep; - if (second_order) - LinSysRes(iPoint,iVar) += ( 3.0*U_time_nP1[iVar] - 4.0*U_time_n[iVar] - +1.0*U_time_nM1[iVar])*Volume_nP1 / (2.0*TimeStep); - } - - /*--- 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)); - } - } - - } - - else { - - /*--- For unsteady flows on dynamic meshes (rigidly transforming or - dynamically deforming), the Geometric Conservation Law (GCL) should be - satisfied in conjunction with the ALE formulation of the governing - equations. The GCL prevents accuracy issues caused by grid motion, i.e. - a uniform free-stream should be preserved through a moving grid. First, - we will loop over the edges and boundaries to compute the GCL component - of the dual time source term that depends on grid velocities. ---*/ - - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - GridVel_i = geometry->nodes->GetGridVel(iPoint); - U_time_n = nodes->GetSolution_time_n(iPoint); - - for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { - - iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); - Normal = geometry->edges->GetNormal(iEdge); - - jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); - GridVel_j = geometry->nodes->GetGridVel(jPoint); - - /*--- Determine whether to consider the normal outward or inward. ---*/ - su2double dir = (geometry->edges->GetNode(iEdge,0) == iPoint)? 0.5 : -0.5; - - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL += dir*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - - for (iVar = 0; iVar < nVar; iVar++) - LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; - } - } - - /*--- Loop over the boundary edges ---*/ - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Get the index for node i plus the boundary face normal ---*/ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - - /*--- Grid velocities stored at boundary node i ---*/ - - GridVel_i = geometry->nodes->GetGridVel(iPoint); - - /*--- Compute the GCL term by dotting the grid velocity with the face - normal. The normal is negated to match the boundary convention. ---*/ - - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; - - /*--- Compute the GCL component of the source term for node i ---*/ - - U_time_n = nodes->GetSolution_time_n(iPoint); - for (iVar = 0; iVar < nVar; iVar++) - LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; - } - } - } - - /*--- Loop over all nodes (excluding halos) to compute the remainder - of the dual time-stepping source term. ---*/ - - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that - we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, - previous solutions that are stored in memory. ---*/ - - U_time_nM1 = nodes->GetSolution_time_n1(iPoint); - U_time_n = nodes->GetSolution_time_n(iPoint); - U_time_nP1 = nodes->GetSolution(iPoint); - - /*--- 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 - volumes will remain constant. ---*/ - - Volume_nM1 = geometry->nodes->GetVolume_nM1(iPoint); - Volume_nP1 = geometry->nodes->GetVolume(iPoint); - - /*--- Compute the dual time-stepping source residual. Due to the - introduction of the GCL term above, the remainder of the source residual - due to the time discretization has a new form.---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - if (first_order) - LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); - if (second_order) - LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) - + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); - } - - /*--- 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)); - } - } - } - -} - void CEulerSolver::PrintVerificationError(const CConfig *config) const { if ((rank != MASTER_NODE) || (MGLevel != MESH_0)) return; @@ -10183,7 +9431,7 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig } -void CEulerSolver::SetFreeStream_Solution(CConfig *config) { +void CEulerSolver::SetFreeStream_Solution(const CConfig *config) { unsigned long iPoint; unsigned short iDim; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index e3df3da16efb..f223f820c5f6 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -41,7 +41,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned * being called by itself, or by its derived class CIncNSSolver. ---*/ const string description = navier_stokes? "Navier-Stokes" : "Euler"; - unsigned short iVar, iMarker, nLineLets; + unsigned short iMarker, nLineLets; ifstream restart_file; unsigned short nZone = geometry->GetnZone(); bool restart = (config->GetRestart() || config->GetRestart_Flow()); @@ -132,22 +132,6 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned SetVerificationSolution(nDim, nVar, config); - /// TODO: This type of variables will be replaced. - - AllocateTerribleLegacyTemporaryVariables(); - - /*--- Define some auxiliary vectors related to the primitive solution ---*/ - - Primitive = new su2double[nPrimVar] (); - Primitive_i = new su2double[nPrimVar] (); - Primitive_j = new su2double[nPrimVar] (); - - /*--- Allocate preconditioning matrix. ---*/ - - Preconditioner = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar ++) - Preconditioner[iVar] = new su2double[nVar]; - /*--- Allocate base class members. ---*/ Allocate(*config); @@ -226,18 +210,6 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned CIncEulerSolver::~CIncEulerSolver(void) { - unsigned short iVar; - - delete [] Primitive; - delete [] Primitive_i; - delete [] Primitive_j; - - if (Preconditioner != nullptr) { - for (iVar = 0; iVar < nVar; iVar ++) - delete [] Preconditioner[iVar]; - delete [] Preconditioner; - } - delete FluidModel; } @@ -919,50 +891,39 @@ void CIncEulerSolver::SetInitialCondition(CGeometry **geometry, CSolver ***solve } } -void CIncEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { +void CIncEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - unsigned long ErrorCounter = 0; - - unsigned long InnerIter = config->GetInnerIter(); - bool cont_adjoint = config->GetContinuous_Adjoint(); - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool muscl = (config->GetMUSCL_Flow() || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == ROE)); - bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); - bool center = ((config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED)); - bool center_jst = center && (config->GetKind_Centered_Flow() == JST); - bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; - bool outlet = ((config->GetnMarker_Outlet() != 0)); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); + const bool center_jst = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); + const bool outlet = (config->GetnMarker_Outlet() != 0); /*--- Set the primitive variables ---*/ - ErrorCounter = SetPrimitive_Variables(solver_container, config, Output); - - /*--- Upwind second order reconstruction ---*/ - - if ((muscl && !center) && (iMesh == MESH_0) && !Output) { - - /*--- Gradient computation for MUSCL reconstruction. ---*/ - - if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) - SetPrimitive_Gradient_GG(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); + SU2_OMP_MASTER + ErrorCounter = 0; + SU2_OMP_BARRIER - /*--- Limiter computation ---*/ + SU2_OMP_ATOMIC + ErrorCounter += SetPrimitive_Variables(solver_container, config); - if ((limiter) && (iMesh == MESH_0) && !Output && !van_albada) { - SetPrimitive_Limiter(geometry, config); + if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) { + SU2_OMP_BARRIER + SU2_OMP_MASTER + { + unsigned long tmp = ErrorCounter; + SU2_MPI::Allreduce(&tmp, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); + config->SetNonphysical_Points(ErrorCounter); } - + SU2_OMP_BARRIER } /*--- Artificial dissipation ---*/ if (center && !Output) { SetMax_Eigenvalue(geometry, config); - if ((center_jst) && (iMesh == MESH_0)) { + if (center_jst) { SetCentered_Dissipation_Sensor(geometry, config); SetUndivided_Laplacian(geometry, config); } @@ -976,251 +937,120 @@ void CIncEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contai if (outlet) GetOutlet_Properties(geometry, config, iMesh, Output); - /*--- Initialize the Jacobian matrices ---*/ - - if (implicit && !Output) Jacobian.SetValZero(); + /*--- Initialize the Jacobian matrix and residual, not needed for the reducer strategy + * as we set blocks (including diagonal ones) and completely overwrite. ---*/ - /*--- Error message ---*/ - - if (config->GetComm_Level() == COMM_FULL) { -#ifdef HAVE_MPI - unsigned long MyErrorCounter = ErrorCounter; ErrorCounter = 0; - SU2_MPI::Allreduce(&MyErrorCounter, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); -#endif - if (iMesh == MESH_0) config->SetNonphysical_Points(ErrorCounter); + if(!ReducerStrategy && !Output) { + LinSysRes.SetValZero(); + if (implicit) Jacobian.SetValZero(); + else {SU2_OMP_BARRIER} // because of "nowait" in LinSysRes } - } -void CIncEulerSolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, - unsigned short iMesh) { } +void CIncEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { -unsigned long CIncEulerSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { + const auto InnerIter = config->GetInnerIter(); + const bool muscl = config->GetMUSCL_Flow() && (iMesh == MESH_0); + const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); + const bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); + const bool van_albada = (config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE); - unsigned long iPoint, nonPhysicalPoints = 0; - bool physical = true; - - for (iPoint = 0; iPoint < nPoint; iPoint ++) { + /*--- Common preprocessing steps. ---*/ - /*--- Incompressible flow, primitive variables ---*/ + CommonPreprocessing(geometry, solver_container, config, iMesh, iRKStep, RunTime_EqSystem, Output); - physical = nodes->SetPrimVar(iPoint,FluidModel); + /*--- Upwind second order reconstruction ---*/ - /* Check for non-realizable states for reporting. */ + if (!Output && muscl && !center) { - if (!physical) nonPhysicalPoints++; + /*--- Gradient computation for MUSCL reconstruction. ---*/ - /*--- Initialize the convective, source and viscous residual vector ---*/ + switch (config->GetKind_Gradient_Method_Recon()) { + case GREEN_GAUSS: + SetPrimitive_Gradient_GG(geometry, config, true); break; + case LEAST_SQUARES: + case WEIGHTED_LEAST_SQUARES: + SetPrimitive_Gradient_LS(geometry, config, true); break; + default: break; + } - if (!Output) LinSysRes.SetBlock_Zero(iPoint); + /*--- Limiter computation ---*/ + if (limiter && !van_albada) SetPrimitive_Limiter(geometry, config); } - - return nonPhysicalPoints; } -void CIncEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config, - unsigned short iMesh, unsigned long Iteration) { - - su2double Area, Vol, Mean_SoundSpeed = 0.0, Mean_ProjVel = 0.0, - Mean_BetaInc2, Lambda, Local_Delta_Time, - Global_Delta_Time = 1E6, Global_Delta_UnstTimeND, ProjVel, ProjVel_i, ProjVel_j; - const su2double* Normal; - - unsigned long iEdge, iVertex, iPoint, jPoint; - unsigned short iDim, iMarker; - - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool time_stepping = config->GetTime_Marching() == TIME_STEPPING; - bool dual_time = ((config->GetTime_Marching() == DT_STEPPING_1ST) || - (config->GetTime_Marching() == DT_STEPPING_2ND)); - - Min_Delta_Time = 1.E30; Max_Delta_Time = 0.0; - - /*--- Set maximum inviscid eigenvalue to zero, and compute sound speed ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) - nodes->SetMax_Lambda_Inv(iPoint,0.0); - - /*--- Loop interior edges ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->edges->GetNode(iEdge,0); - jPoint = geometry->edges->GetNode(iEdge,1); - - Normal = geometry->edges->GetNormal(iEdge); - - Area = GeometryToolbox::Norm(nDim, Normal); +unsigned long CIncEulerSolver::SetPrimitive_Variables(CSolver **solver_container, const CConfig *config) { - /*--- Mean Values ---*/ - - Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); - Mean_BetaInc2 = 0.5 * (nodes->GetBetaInc2(iPoint) + nodes->GetBetaInc2(jPoint)); - Mean_SoundSpeed = sqrt(Mean_BetaInc2*Area*Area); + unsigned long iPoint, nonPhysicalPoints = 0; - /*--- Adjustment for grid movement ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPoint; iPoint ++) { - if (dynamic_grid) { - su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); - su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); - ProjVel_i = 0.0; ProjVel_j = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjVel_i += GridVel_i[iDim]*Normal[iDim]; - ProjVel_j += GridVel_j[iDim]*Normal[iDim]; - } - Mean_ProjVel -= 0.5 * (ProjVel_i + ProjVel_j); - } + /*--- Incompressible flow, primitive variables ---*/ - /*--- Inviscid contribution ---*/ + auto physical = nodes->SetPrimVar(iPoint,FluidModel); - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) nodes->AddMax_Lambda_Inv(iPoint,Lambda); - if (geometry->nodes->GetDomain(jPoint)) nodes->AddMax_Lambda_Inv(jPoint,Lambda); + /* Check for non-realizable states for reporting. */ + if (!physical) nonPhysicalPoints++; } - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - Mean_BetaInc2 = nodes->GetBetaInc2(iPoint); - Mean_SoundSpeed = sqrt(Mean_BetaInc2*Area*Area); - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - su2double *GridVel = geometry->nodes->GetGridVel(iPoint); - ProjVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjVel += GridVel[iDim]*Normal[iDim]; - Mean_ProjVel -= ProjVel; - } - - /*--- Inviscid contribution ---*/ + return nonPhysicalPoints; +} - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) { - nodes->AddMax_Lambda_Inv(iPoint,Lambda); - } +void CIncEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config, + unsigned short iMesh, unsigned long Iteration) { + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return sqrt(0.5 * (nodes.GetBetaInc2(iPoint) + nodes.GetBetaInc2(jPoint))); } - } - } - - /*--- Local time-stepping: each element uses their own speed for steady state - simulations or for pseudo time steps in a dual time simulation. ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - Vol = geometry->nodes->GetVolume(iPoint); - if (Vol != 0.0) { - Local_Delta_Time = nodes->GetLocalCFL(iPoint)*Vol / nodes->GetMax_Lambda_Inv(iPoint); - Global_Delta_Time = min(Global_Delta_Time, Local_Delta_Time); - Min_Delta_Time = min(Min_Delta_Time, Local_Delta_Time); - Max_Delta_Time = max(Max_Delta_Time, Local_Delta_Time); - if (Local_Delta_Time > config->GetMax_DeltaTime()) - Local_Delta_Time = config->GetMax_DeltaTime(); - nodes->SetDelta_Time(iPoint,Local_Delta_Time); - } - else { - nodes->SetDelta_Time(iPoint,0.0); + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint) const { + return sqrt(nodes.GetBetaInc2(iPoint)); } - } + } soundSpeed; - /*--- Compute the max and the min dt (in parallel) ---*/ + /*--- Define an object to compute the viscous eigenvalue. ---*/ + struct LambdaVisc { + const bool energy; - if (config->GetComm_Level() == COMM_FULL) { -#ifdef HAVE_MPI - su2double rbuf_time, sbuf_time; - sbuf_time = Min_Delta_Time; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Min_Delta_Time = rbuf_time; - - sbuf_time = Max_Delta_Time; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MAX, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Max_Delta_Time = rbuf_time; -#endif - } + LambdaVisc(bool e) : energy(e) {} - /*--- For time-accurate simulations use the minimum delta time of the whole mesh (global) ---*/ - - if (time_stepping) { -#ifdef HAVE_MPI - su2double rbuf_time, sbuf_time; - sbuf_time = Global_Delta_Time; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Global_Delta_Time = rbuf_time; -#endif - /*--- If the unsteady CFL is set to zero, it uses the defined - unsteady time step, otherwise it computes the time step based - on the unsteady CFL ---*/ - - if (config->GetUnst_CFL() == 0.0) { - Global_Delta_Time = config->GetDelta_UnstTime(); + FORCEINLINE su2double lambda(su2double lamVisc, su2double eddyVisc, su2double rho, su2double k, su2double cv) const { + su2double Lambda_1 = (4.0/3.0)*(lamVisc + eddyVisc); + su2double Lambda_2 = 0.0; + if (energy) Lambda_2 = k / cv; + return (Lambda_1 + Lambda_2) / rho; } - config->SetDelta_UnstTimeND(Global_Delta_Time); - for (iPoint = 0; iPoint < nPointDomain; iPoint++){ - - /*--- Sets the regular CFL equal to the unsteady CFL ---*/ - - nodes->SetLocalCFL(iPoint, config->GetUnst_CFL()); - nodes->SetDelta_Time(iPoint, Global_Delta_Time); - Min_Delta_Time = Global_Delta_Time; - Max_Delta_Time = Global_Delta_Time; + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + su2double thermalCond = 0.5*(nodes.GetThermalConductivity(iPoint) + nodes.GetThermalConductivity(jPoint)); + su2double laminarVisc = 0.5*(nodes.GetLaminarViscosity(iPoint) + nodes.GetLaminarViscosity(jPoint)); + su2double eddyVisc = 0.5*(nodes.GetEddyViscosity(iPoint) + nodes.GetEddyViscosity(jPoint)); + su2double density = 0.5*(nodes.GetDensity(iPoint) + nodes.GetDensity(jPoint)); + su2double cv = 0.5*(nodes.GetSpecificHeatCv(iPoint) + nodes.GetSpecificHeatCv(jPoint)); + return lambda(laminarVisc, eddyVisc, density, thermalCond, cv); } - } - - /*--- Recompute the unsteady time step for the dual time strategy - if the unsteady CFL is diferent from 0 ---*/ - - if ((dual_time) && (Iteration == 0) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) { - Global_Delta_UnstTimeND = 1e30; - for (iPoint = 0; iPoint < nPointDomain; iPoint++){ - Global_Delta_UnstTimeND = min(Global_Delta_UnstTimeND,config->GetUnst_CFL()*Global_Delta_Time/nodes->GetLocalCFL(iPoint)); + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint) const { + su2double thermalCond = nodes.GetThermalConductivity(iPoint); + su2double laminarVisc = nodes.GetLaminarViscosity(iPoint); + su2double eddyVisc = nodes.GetEddyViscosity(iPoint); + su2double density = nodes.GetDensity(iPoint); + su2double cv = nodes.GetSpecificHeatCv(iPoint); + return lambda(laminarVisc, eddyVisc, density, thermalCond, cv); } -#ifdef HAVE_MPI - su2double rbuf_time, sbuf_time; - sbuf_time = Global_Delta_UnstTimeND; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Global_Delta_UnstTimeND = rbuf_time; -#endif - config->SetDelta_UnstTimeND(Global_Delta_UnstTimeND); - } + } lambdaVisc(config->GetEnergy_Equation()); - /*--- The pseudo local time (explicit integration) cannot be greater than the physical time ---*/ + /*--- Now instantiate the generic implementation with the two functors above. ---*/ - if (dual_time) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - if (!implicit) { - Local_Delta_Time = min((2.0/3.0)*config->GetDelta_UnstTimeND(), nodes->GetDelta_Time(iPoint)); - nodes->SetDelta_Time(iPoint,Local_Delta_Time); - } - } + SetTime_Step_impl(soundSpeed, lambdaVisc, geometry, solver_container, config, iMesh, Iteration); } @@ -1286,8 +1116,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont CNumerics* numerics = numerics_container[CONV_TERM]; - su2double **Gradient_i, **Gradient_j, Project_Grad_i, Project_Grad_j, - *V_i, *V_j, *S_i, *S_j, *Limiter_i = nullptr, *Limiter_j = nullptr; + /*--- Static arrays of MUSCL-reconstructed primitives and secondaries (thread safety). ---*/ + su2double Primitive_i[MAXNVAR] = {0.0}, Primitive_j[MAXNVAR] = {0.0}; unsigned long iEdge, iPoint, jPoint, counter_local = 0, counter_global = 0; unsigned short iDim, iVar; @@ -1314,20 +1144,25 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont /*--- Get primitive variables ---*/ - V_i = nodes->GetPrimitive(iPoint); V_j = nodes->GetPrimitive(jPoint); - S_i = nodes->GetSecondary(iPoint); S_j = nodes->GetSecondary(jPoint); + auto V_i = nodes->GetPrimitive(iPoint); + auto V_j = nodes->GetPrimitive(jPoint); /*--- High order reconstruction using MUSCL strategy ---*/ if (muscl) { + auto Coord_i = geometry->nodes->GetCoord(iPoint); + auto Coord_j = geometry->nodes->GetCoord(jPoint); + + su2double Vector_ij[MAXNDIM] = {0.0}; for (iDim = 0; iDim < nDim; iDim++) { - Vector_i[iDim] = 0.5*(geometry->nodes->GetCoord(jPoint, iDim) - geometry->nodes->GetCoord(iPoint, iDim)); - Vector_j[iDim] = 0.5*(geometry->nodes->GetCoord(iPoint, iDim) - geometry->nodes->GetCoord(jPoint, iDim)); + Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]); } - Gradient_i = nodes->GetGradient_Reconstruction(iPoint); - Gradient_j = nodes->GetGradient_Reconstruction(jPoint); + auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint); + auto Gradient_j = nodes->GetGradient_Reconstruction(jPoint); + + su2double *Limiter_i = nullptr, *Limiter_j = nullptr; if (limiter) { Limiter_i = nodes->GetLimiter_Primitive(iPoint); @@ -1335,10 +1170,10 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont } for (iVar = 0; iVar < nPrimVarGrad; iVar++) { - Project_Grad_i = 0.0; Project_Grad_j = 0.0; + su2double Project_Grad_i = 0.0, Project_Grad_j = 0.0; for (iDim = 0; iDim < nDim; iDim++) { - Project_Grad_i += Vector_i[iDim]*Gradient_i[iVar][iDim]; - Project_Grad_j += Vector_j[iDim]*Gradient_j[iVar][iDim]; + Project_Grad_i += Vector_ij[iDim]*Gradient_i[iVar][iDim]; + Project_Grad_j -= Vector_ij[iDim]*Gradient_j[iVar][iDim]; } if (limiter) { if (van_albada){ @@ -1373,17 +1208,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont bool neg_density_i = (Primitive_i[nDim+2] < 0.0); bool neg_density_j = (Primitive_j[nDim+2] < 0.0); - if (neg_density_i || neg_temperature_i) { - nodes->SetNon_Physical(iPoint, true); - } else { - nodes->SetNon_Physical(iPoint, false); - } - - if (neg_density_j || neg_temperature_j) { - nodes->SetNon_Physical(jPoint, true); - } else { - nodes->SetNon_Physical(jPoint, false); - } + nodes->SetNon_Physical(iPoint, neg_density_i || neg_temperature_i); + nodes->SetNon_Physical(jPoint, neg_density_j || neg_temperature_j); /* Lastly, check for existing first-order points still active from previous iterations. */ @@ -1404,10 +1230,9 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont } else { - /*--- Set conservative variables without reconstruction ---*/ + /*--- Set primitive variables without reconstruction ---*/ numerics->SetPrimitive(V_i, V_j); - numerics->SetSecondary(S_i, S_j); } @@ -1724,115 +1549,27 @@ void CIncEulerSolver::Source_Template(CGeometry *geometry, CSolver **solver_cont } -void CIncEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, CConfig *config) { - - su2double Area, Mean_SoundSpeed = 0.0, Mean_ProjVel = 0.0, - Mean_BetaInc2, Lambda, ProjVel, ProjVel_i, ProjVel_j, *GridVel, *GridVel_i, *GridVel_j; - const su2double* Normal; - - unsigned long iEdge, iVertex, iPoint, jPoint; - unsigned short iDim, iMarker; - - /*--- Set maximum inviscid eigenvalue to zero, and compute sound speed ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - nodes->SetLambda(iPoint,0.0); - } - - /*--- Loop interior edges ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->edges->GetNode(iEdge,0); - jPoint = geometry->edges->GetNode(iEdge,1); - - Normal = geometry->edges->GetNormal(iEdge); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); - Mean_BetaInc2 = 0.5 * (nodes->GetBetaInc2(iPoint) + nodes->GetBetaInc2(jPoint)); - Mean_SoundSpeed = sqrt(Mean_BetaInc2*Area*Area); +void CIncEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, const CConfig *config) { - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - GridVel_i = geometry->nodes->GetGridVel(iPoint); - GridVel_j = geometry->nodes->GetGridVel(jPoint); - ProjVel_i = 0.0; ProjVel_j =0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjVel_i += GridVel_i[iDim]*Normal[iDim]; - ProjVel_j += GridVel_j[iDim]*Normal[iDim]; - } - Mean_ProjVel -= 0.5 * (ProjVel_i + ProjVel_j); + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return sqrt(0.5 * (nodes.GetBetaInc2(iPoint) + nodes.GetBetaInc2(jPoint))); } - /*--- Inviscid contribution ---*/ - - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) nodes->AddLambda(iPoint,Lambda); - if (geometry->nodes->GetDomain(jPoint)) nodes->AddLambda(jPoint,Lambda); - - } - - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - Mean_BetaInc2 = nodes->GetBetaInc2(iPoint); - Mean_SoundSpeed = sqrt(Mean_BetaInc2*Area*Area); - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - ProjVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjVel += GridVel[iDim]*Normal[iDim]; - Mean_ProjVel -= ProjVel; - } - - /*--- Inviscid contribution ---*/ - - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) { - nodes->AddLambda(iPoint,Lambda); - } - - } + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint) const { + return sqrt(nodes.GetBetaInc2(iPoint)); } - } - /*--- Correct the eigenvalue values across any periodic boundaries. ---*/ + } soundSpeed; - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_MAX_EIG); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_MAX_EIG); - } + /*--- Instantiate generic implementation. ---*/ - /*--- MPI parallelization ---*/ - - InitiateComms(geometry, config, MAX_EIGENVALUE); - CompleteComms(geometry, config, MAX_EIGENVALUE); + SetMax_Eigenvalue_impl(soundSpeed, geometry, config); } -void CIncEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *config) { +void CIncEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, const CConfig *config) { unsigned long iPoint, jPoint, iEdge; su2double *Diff; @@ -1891,306 +1628,86 @@ void CIncEulerSolver::SetUndivided_Laplacian(CGeometry *geometry, CConfig *confi } -void CIncEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, CConfig *config) { - - unsigned long iEdge, iPoint, jPoint; - su2double Pressure_i = 0.0, Pressure_j = 0.0; - bool boundary_i, boundary_j; - - /*--- Reset variables to store the undivided pressure ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - iPoint_UndLapl[iPoint] = 0.0; - jPoint_UndLapl[iPoint] = 0.0; - } - - /*--- Evaluate the pressure sensor ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - iPoint = geometry->edges->GetNode(iEdge,0); - jPoint = geometry->edges->GetNode(iEdge,1); - - /*--- Get the pressure, or density for incompressible solvers ---*/ - - Pressure_i = nodes->GetDensity(iPoint); - Pressure_j = nodes->GetDensity(jPoint); - - boundary_i = geometry->nodes->GetPhysicalBoundary(iPoint); - boundary_j = geometry->nodes->GetPhysicalBoundary(jPoint); - - /*--- Both points inside the domain, or both on the boundary ---*/ - - if ((!boundary_i && !boundary_j) || (boundary_i && boundary_j)) { - - if (geometry->nodes->GetDomain(iPoint)) { - iPoint_UndLapl[iPoint] += (Pressure_j - Pressure_i); - jPoint_UndLapl[iPoint] += (Pressure_i + Pressure_j); - } - - if (geometry->nodes->GetDomain(jPoint)) { - iPoint_UndLapl[jPoint] += (Pressure_i - Pressure_j); - jPoint_UndLapl[jPoint] += (Pressure_i + Pressure_j); - } +void CIncEulerSolver::SetCentered_Dissipation_Sensor(CGeometry *geometry, const CConfig *config) { + /*--- Define an object for the sensor variable, density. ---*/ + struct SensVar { + FORCEINLINE su2double operator() (const CIncEulerVariable& nodes, unsigned long iPoint) const { + return nodes.GetDensity(iPoint); } + } sensVar; - /*--- iPoint inside the domain, jPoint on the boundary ---*/ - - if (!boundary_i && boundary_j) - if (geometry->nodes->GetDomain(iPoint)) { - iPoint_UndLapl[iPoint] += (Pressure_j - Pressure_i); - jPoint_UndLapl[iPoint] += (Pressure_i + Pressure_j); - } - - /*--- jPoint inside the domain, iPoint on the boundary ---*/ - - if (boundary_i && !boundary_j) - if (geometry->nodes->GetDomain(jPoint)) { - iPoint_UndLapl[jPoint] += (Pressure_i - Pressure_j); - jPoint_UndLapl[jPoint] += (Pressure_i + Pressure_j); - } - - } - - /*--- Correct the sensor values across any periodic boundaries. ---*/ - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_SENSOR); - } - - /*--- Set pressure switch for each point ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) - nodes->SetSensor(iPoint,fabs(iPoint_UndLapl[iPoint]) / jPoint_UndLapl[iPoint]); - - /*--- MPI parallelization ---*/ - - InitiateComms(geometry, config, SENSOR); - CompleteComms(geometry, config, SENSOR); + /*--- Instantiate generic implementation. ---*/ + SetCentered_Dissipation_Sensor_impl(sensVar, geometry, config); } -void CIncEulerSolver::ExplicitRK_Iteration(CGeometry *geometry, CSolver **solver_container, - CConfig *config, unsigned short iRKStep) { - - su2double *Residual, *Res_TruncError, Vol, Delta, Res; - unsigned short iVar, jVar; - unsigned long iPoint; - - su2double RK_AlphaCoeff = config->Get_Alpha_RKStep(iRKStep); - bool adjoint = config->GetContinuous_Adjoint(); +template +FORCEINLINE void CIncEulerSolver::Explicit_Iteration(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iRKStep) { + struct Precond { + const CIncEulerSolver* solver; + su2activematrix matrix; + unsigned short nVar; - for (iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - - /*--- Update the solution ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - Vol = (geometry->nodes->GetVolume(iPoint) + - geometry->nodes->GetPeriodicVolume(iPoint)); - Delta = nodes->GetDelta_Time(iPoint) / Vol; - - Res_TruncError = nodes->GetResTruncError(iPoint); - Residual = LinSysRes.GetBlock(iPoint); - - if (!adjoint) { - SetPreconditioner(config, iPoint); - for (iVar = 0; iVar < nVar; iVar ++ ) { - Res = 0.0; - for (jVar = 0; jVar < nVar; jVar ++ ) - Res += Preconditioner[iVar][jVar]*(Residual[jVar] + Res_TruncError[jVar]); - nodes->AddSolution(iPoint,iVar, -Res*Delta*RK_AlphaCoeff); - AddRes_RMS(iVar, Res*Res); - AddRes_Max(iVar, fabs(Res), geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetCoord(iPoint)); - } + Precond(const CIncEulerSolver* s, unsigned short n) : solver(s), nVar(n) { + matrix.resize(nVar,nVar); } - } - /*--- MPI solution ---*/ + FORCEINLINE void compute(const CConfig* config, unsigned long iPoint) { + solver->SetPreconditioner(config, iPoint, 1.0, matrix); + } - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + FORCEINLINE su2double apply(unsigned short iVar, const su2double* res, const su2double* resTrunc) const { + su2double resPrec = 0.0; + for (unsigned short jVar = 0; jVar < nVar; ++jVar) + resPrec += matrix(iVar,jVar) * (res[jVar] + resTrunc[jVar]); + return resPrec; + } + } precond(this, nVar); - /*--- Compute the root mean square residual ---*/ + Explicit_Iteration_impl(precond, geometry, solver_container, config, iRKStep); +} - SetResidual_RMS(geometry, config); +void CIncEulerSolver::ExplicitRK_Iteration(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iRKStep) { - /*--- For verification cases, compute the global error metrics. ---*/ + Explicit_Iteration(geometry, solver_container, config, iRKStep); +} - ComputeVerificationError(geometry, config); +void CIncEulerSolver::ClassicalRK4_Iteration(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iRKStep) { + Explicit_Iteration(geometry, solver_container, config, iRKStep); } void CIncEulerSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - su2double *local_Residual, *local_Res_TruncError, Vol, Delta, Res; - unsigned short iVar, jVar; - unsigned long iPoint; - - bool adjoint = config->GetContinuous_Adjoint(); - - for (iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - - /*--- Update the solution ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - Vol = (geometry->nodes->GetVolume(iPoint) + - geometry->nodes->GetPeriodicVolume(iPoint)); - Delta = nodes->GetDelta_Time(iPoint) / Vol; - - local_Res_TruncError = nodes->GetResTruncError(iPoint); - local_Residual = LinSysRes.GetBlock(iPoint); - - - if (!adjoint) { - SetPreconditioner(config, iPoint); - for (iVar = 0; iVar < nVar; iVar ++ ) { - Res = 0.0; - for (jVar = 0; jVar < nVar; jVar ++ ) - Res += Preconditioner[iVar][jVar]*(local_Residual[jVar] + local_Res_TruncError[jVar]); - nodes->AddSolution(iPoint,iVar, -Res*Delta); - AddRes_RMS(iVar, Res*Res); - AddRes_Max(iVar, fabs(Res), geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetCoord(iPoint)); - } - } - } - - /*--- MPI solution ---*/ - - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - - /*--- Compute the root mean square residual ---*/ - - SetResidual_RMS(geometry, config); - - /*--- For verification cases, compute the global error metrics. ---*/ - - ComputeVerificationError(geometry, config); - + Explicit_Iteration(geometry, solver_container, config, 0); } void CIncEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - unsigned short iVar, jVar; - unsigned long iPoint, total_index, IterLinSol = 0; - su2double Delta, *local_Res_TruncError, Vol; - - bool adjoint = config->GetContinuous_Adjoint(); - - /*--- Set maximum residual to zero ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - - /*--- Build implicit system ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Read the residual ---*/ - - local_Res_TruncError = nodes->GetResTruncError(iPoint); - - /*--- Read the volume ---*/ - - Vol = (geometry->nodes->GetVolume(iPoint) + - geometry->nodes->GetPeriodicVolume(iPoint)); - - /*--- Apply the preconditioner and add to the diagonal. ---*/ - - if (nodes->GetDelta_Time(iPoint) != 0.0) { - Delta = Vol / nodes->GetDelta_Time(iPoint); - SetPreconditioner(config, iPoint); - for (iVar = 0; iVar < nVar; iVar ++ ) { - for (jVar = 0; jVar < nVar; jVar ++ ) { - Preconditioner[iVar][jVar] = Delta*Preconditioner[iVar][jVar]; - } - } - Jacobian.AddBlock2Diag(iPoint, Preconditioner); - } else { - Jacobian.SetVal2Diag(iPoint, 1.0); - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = 0.0; - local_Res_TruncError[iVar] = 0.0; - } - } - - /*--- Right hand side of the system (-Residual) and initial guess (x = 0) ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]); - LinSysSol[total_index] = 0.0; - AddRes_RMS(iVar, LinSysRes[total_index]*LinSysRes[total_index]); - AddRes_Max(iVar, fabs(LinSysRes[total_index]), geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetCoord(iPoint)); - } - - } - - /*--- Initialize residual and solution at the ghost points ---*/ - - for (iPoint = nPointDomain; iPoint < nPoint; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - } - } - - /*--- Solve or smooth the linear system ---*/ - - IterLinSol = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); - - /*--- Store the value of the residual. ---*/ - - SetResLinSolver(System.GetResidual()); + struct IncPrec { + const CIncEulerSolver* solver; + const bool active = true; + su2activematrix matrix; - /*--- The the number of iterations of the linear solver ---*/ + IncPrec(const CIncEulerSolver* s, unsigned short nVar) : solver(s) { matrix.resize(nVar,nVar); } - SetIterLinSolver(IterLinSol); - - /*--- Update solution (system written in terms of increments) ---*/ - - if (!adjoint) { - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - nodes->AddSolution(iPoint, iVar, nodes->GetUnderRelaxation(iPoint)*LinSysSol[iPoint*nVar+iVar]); - } + FORCEINLINE const su2activematrix& operator() (const CConfig* config, unsigned long iPoint, su2double delta) { + solver->SetPreconditioner(config, iPoint, delta, matrix); + return matrix; } - } - - for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic()/2; iPeriodic++) { - InitiatePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); - CompletePeriodicComms(geometry, config, iPeriodic, PERIODIC_IMPLICIT); - } - - /*--- MPI solution ---*/ - - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); + } precond(this, nVar); - /*--- For verification cases, compute the global error metrics. ---*/ - - ComputeVerificationError(geometry, config); + ImplicitEuler_Iteration_impl(precond, geometry, solver_container, config, false); } void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_container, - CConfig *config, unsigned short iMesh) { + CConfig *config, unsigned short iMesh) { su2double epsilon2 = config->GetBeta_Factor(); su2double epsilon2_default = 4.1; @@ -2214,10 +1731,8 @@ void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_co /*--- Communicate the max globally to give a conservative estimate. ---*/ -#ifdef HAVE_MPI su2double myMaxVel2 = maxVel2; maxVel2 = 0.0; SU2_MPI::Allreduce(&myMaxVel2, &maxVel2, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); -#endif Beta = max(1e-10,maxVel2); config->SetMax_Vel2(Beta); @@ -2233,12 +1748,13 @@ void CIncEulerSolver::SetBeta_Parameter(CGeometry *geometry, CSolver **solver_co } -void CIncEulerSolver::SetPreconditioner(CConfig *config, unsigned long iPoint) { +void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPoint, + su2double delta, su2activematrix& Preconditioner) const { - unsigned short iDim, jDim; + unsigned short iDim, jDim, iVar, jVar; su2double BetaInc2, Density, dRhodT, Temperature, oneOverCp, Cp; - su2double Velocity[3] = {0.0,0.0,0.0}; + su2double Velocity[MAXNDIM] = {0.0}; bool variable_density = (config->GetKind_DensityModel() == VARIABLE); bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -2296,6 +1812,10 @@ void CIncEulerSolver::SetPreconditioner(CConfig *config, unsigned long iPoint) { if (energy) Preconditioner[nDim+1][nDim+1] = Cp*(dRhodT*Temperature + Density); else Preconditioner[nDim+1][nDim+1] = 1.0; + for (iVar = 0; iVar < nVar; iVar ++ ) + for (jVar = 0; jVar < nVar; jVar ++ ) + Preconditioner[iVar][jVar] = delta*Preconditioner[iVar][jVar]; + } else { /*--- For explicit calculations, we move the residual to the @@ -2337,132 +1857,123 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain unsigned short iDim; unsigned long iVertex, iPoint, Point_Normal; - su2double *V_infty, *V_domain; + const bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; + const bool viscous = config->GetViscous(); - bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; - bool viscous = config->GetViscous(); - - su2double *Normal = new su2double[nDim]; + su2double Normal[MAXNDIM] = {0.0}; /*--- Loop over all the vertices on this boundary marker ---*/ for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - /*--- Allocate the value at the infinity ---*/ - - V_infty = GetCharacPrimVar(val_marker, iVertex); - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Index of the closest interior node ---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + /*--- Allocate the value at the infinity ---*/ - /*--- Normal vector for this vertex (negate for outward convention) ---*/ + auto V_infty = GetCharacPrimVar(val_marker, iVertex); - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); + /*--- Index of the closest interior node ---*/ - /*--- Retrieve solution at the farfield boundary node ---*/ + Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - V_domain = nodes->GetPrimitive(iPoint); + /*--- Normal vector for this vertex (negate for outward convention) ---*/ - /*--- Recompute and store the velocity in the primitive variable vector. ---*/ + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) - V_infty[iDim+1] = GetVelocity_Inf(iDim); + /*--- Retrieve solution at the farfield boundary node ---*/ - /*--- Far-field pressure set to static pressure (0.0). ---*/ + auto V_domain = nodes->GetPrimitive(iPoint); - V_infty[0] = GetPressure_Inf(); + /*--- Recompute and store the velocity in the primitive variable vector. ---*/ - /*--- Dirichlet condition for temperature at far-field (if energy is active). ---*/ + for (iDim = 0; iDim < nDim; iDim++) + V_infty[iDim+1] = GetVelocity_Inf(iDim); - V_infty[nDim+1] = GetTemperature_Inf(); + /*--- Far-field pressure set to static pressure (0.0). ---*/ - /*--- Store the density. ---*/ + V_infty[0] = GetPressure_Inf(); - V_infty[nDim+2] = GetDensity_Inf(); + /*--- Dirichlet condition for temperature at far-field (if energy is active). ---*/ - /*--- Beta coefficient stored at the node ---*/ + V_infty[nDim+1] = GetTemperature_Inf(); - V_infty[nDim+3] = nodes->GetBetaInc2(iPoint); + /*--- Store the density. ---*/ - /*--- Cp is needed for Temperature equation. ---*/ + V_infty[nDim+2] = GetDensity_Inf(); - V_infty[nDim+7] = nodes->GetSpecificHeatCp(iPoint); + /*--- Beta coefficient stored at the node ---*/ - /*--- Set various quantities in the numerics class ---*/ + V_infty[nDim+3] = nodes->GetBetaInc2(iPoint); - conv_numerics->SetPrimitive(V_domain, V_infty); + /*--- Cp is needed for Temperature equation. ---*/ - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + V_infty[nDim+7] = nodes->GetSpecificHeatCp(iPoint); - /*--- Compute the convective residual using an upwind scheme ---*/ + /*--- Set various quantities in the numerics class ---*/ - auto residual = conv_numerics->ComputeResidual(config); + conv_numerics->SetPrimitive(V_domain, V_infty); - /*--- Update residual value ---*/ + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); - LinSysRes.AddBlock(iPoint, residual); + /*--- Compute the convective residual using an upwind scheme ---*/ - /*--- Convective Jacobian contribution for implicit integration ---*/ + auto residual = conv_numerics->ComputeResidual(config); - if (implicit) - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + /*--- Update residual value ---*/ - /*--- Viscous residual contribution ---*/ + LinSysRes.AddBlock(iPoint, residual); - if (viscous) { + /*--- Convective Jacobian contribution for implicit integration ---*/ - /*--- Set transport properties at infinity. ---*/ + if (implicit) + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - V_infty[nDim+4] = nodes->GetLaminarViscosity(iPoint); - V_infty[nDim+5] = nodes->GetEddyViscosity(iPoint); - V_infty[nDim+6] = nodes->GetThermalConductivity(iPoint); + /*--- Viscous residual contribution ---*/ - /*--- Set the normal vector and the coordinates ---*/ + if (!viscous) continue; - visc_numerics->SetNormal(Normal); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - geometry->nodes->GetCoord(Point_Normal)); + /*--- Set transport properties at infinity. ---*/ - /*--- Primitive variables, and gradient ---*/ + V_infty[nDim+4] = nodes->GetLaminarViscosity(iPoint); + V_infty[nDim+5] = nodes->GetEddyViscosity(iPoint); + V_infty[nDim+6] = nodes->GetThermalConductivity(iPoint); - visc_numerics->SetPrimitive(V_domain, V_infty); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - nodes->GetGradient_Primitive(iPoint)); + /*--- Set the normal vector and the coordinates ---*/ - /*--- Turbulent kinetic energy ---*/ + visc_numerics->SetNormal(Normal); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(Point_Normal)); - if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) - visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), - solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); + /*--- Primitive variables, and gradient ---*/ - /*--- Compute and update viscous residual ---*/ + visc_numerics->SetPrimitive(V_domain, V_infty); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), + nodes->GetGradient_Primitive(iPoint)); - auto residual = visc_numerics->ComputeResidual(config); - LinSysRes.SubtractBlock(iPoint, residual); + /*--- Turbulent kinetic energy ---*/ - /*--- Viscous Jacobian contribution for implicit integration ---*/ + if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) + visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), + solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); - if (implicit) - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + /*--- Compute and update viscous residual ---*/ - } + auto residual_v = visc_numerics->ComputeResidual(config); + LinSysRes.SubtractBlock(iPoint, residual_v); - } - } + /*--- Viscous Jacobian contribution for implicit integration ---*/ - /*--- Free locally allocated memory ---*/ + if (implicit) + Jacobian.SubtractBlock2Diag(iPoint, residual_v.jacobian_i); - delete [] Normal; + } } @@ -2473,248 +1984,234 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, unsigned long Point_Normal; su2double *Flow_Dir, Flow_Dir_Mag, Vel_Mag, Area, P_total, P_domain, Vn; su2double *V_inlet, *V_domain; - su2double UnitFlowDir[3] = {0.0,0.0,0.0}; - su2double dV[3] = {0.0,0.0,0.0}; + su2double UnitFlowDir[MAXNDIM] = {0.0}, dV[MAXNDIM] = {0.0}; su2double Damping = config->GetInc_Inlet_Damping(); - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool viscous = config->GetViscous(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool viscous = config->GetViscous(); - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + string Marker_Tag = config->GetMarker_All_TagBound(val_marker); unsigned short Kind_Inlet = config->GetKind_Inc_Inlet(Marker_Tag); - su2double *Normal = new su2double[nDim]; + su2double Normal[MAXNDIM] = {0.0}; /*--- Loop over all the vertices on this boundary marker ---*/ for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - /*--- Allocate the value at the inlet ---*/ - - V_inlet = GetCharacPrimVar(val_marker, iVertex); - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { + if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Index of the closest interior node ---*/ + /*--- Allocate the value at the inlet ---*/ - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + V_inlet = GetCharacPrimVar(val_marker, iVertex); - /*--- Normal vector for this vertex (negate for outward convention) ---*/ + /*--- Index of the closest interior node ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); + Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - Area = GeometryToolbox::Norm(nDim, Normal); + /*--- Normal vector for this vertex (negate for outward convention) ---*/ - /*--- Both types of inlets may use the prescribed flow direction. - Ensure that the flow direction is a unit vector. ---*/ + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); - Flow_Dir = Inlet_FlowDir[val_marker][iVertex]; - Flow_Dir_Mag = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Flow_Dir_Mag += Flow_Dir[iDim]*Flow_Dir[iDim]; - Flow_Dir_Mag = sqrt(Flow_Dir_Mag); + Area = GeometryToolbox::Norm(nDim, Normal); - /*--- Store the unit flow direction vector. ---*/ + /*--- Both types of inlets may use the prescribed flow direction. + Ensure that the flow direction is a unit vector. ---*/ - for (iDim = 0; iDim < nDim; iDim++) - UnitFlowDir[iDim] = Flow_Dir[iDim]/Flow_Dir_Mag; + Flow_Dir = Inlet_FlowDir[val_marker][iVertex]; + Flow_Dir_Mag = GeometryToolbox::Norm(nDim, Flow_Dir); - /*--- Retrieve solution at this boundary node. ---*/ + /*--- Store the unit flow direction vector. ---*/ - V_domain = nodes->GetPrimitive(iPoint); + for (iDim = 0; iDim < nDim; iDim++) + UnitFlowDir[iDim] = Flow_Dir[iDim]/Flow_Dir_Mag; - /*--- Neumann condition for dynamic pressure ---*/ + /*--- Retrieve solution at this boundary node. ---*/ - V_inlet[0] = nodes->GetPressure(iPoint); + V_domain = nodes->GetPrimitive(iPoint); - /*--- The velocity is either prescribed or computed from total pressure. ---*/ + /*--- Neumann condition for dynamic pressure ---*/ - switch (Kind_Inlet) { + V_inlet[0] = nodes->GetPressure(iPoint); - /*--- Velocity and temperature (if required) been specified at the inlet. ---*/ + /*--- The velocity is either prescribed or computed from total pressure. ---*/ - case VELOCITY_INLET: + switch (Kind_Inlet) { - /*--- Retrieve the specified velocity and temperature for the inlet. ---*/ + /*--- Velocity and temperature (if required) been specified at the inlet. ---*/ - Vel_Mag = Inlet_Ptotal[val_marker][iVertex]/config->GetVelocity_Ref(); + case VELOCITY_INLET: - /*--- Store the velocity in the primitive variable vector. ---*/ + /*--- Retrieve the specified velocity and temperature for the inlet. ---*/ - for (iDim = 0; iDim < nDim; iDim++) - V_inlet[iDim+1] = Vel_Mag*UnitFlowDir[iDim]; + Vel_Mag = Inlet_Ptotal[val_marker][iVertex]/config->GetVelocity_Ref(); - /*--- Dirichlet condition for temperature (if energy is active) ---*/ + /*--- Store the velocity in the primitive variable vector. ---*/ - V_inlet[nDim+1] = Inlet_Ttotal[val_marker][iVertex]/config->GetTemperature_Ref(); + for (iDim = 0; iDim < nDim; iDim++) + V_inlet[iDim+1] = Vel_Mag*UnitFlowDir[iDim]; - break; + /*--- Dirichlet condition for temperature (if energy is active) ---*/ - /*--- Stagnation pressure has been specified at the inlet. ---*/ + V_inlet[nDim+1] = Inlet_Ttotal[val_marker][iVertex]/config->GetTemperature_Ref(); - case PRESSURE_INLET: + break; - /*--- Retrieve the specified total pressure for the inlet. ---*/ + /*--- Stagnation pressure has been specified at the inlet. ---*/ - P_total = Inlet_Ptotal[val_marker][iVertex]/config->GetPressure_Ref(); + case PRESSURE_INLET: - /*--- Store the current static pressure for clarity. ---*/ + /*--- Retrieve the specified total pressure for the inlet. ---*/ - P_domain = nodes->GetPressure(iPoint); + P_total = Inlet_Ptotal[val_marker][iVertex]/config->GetPressure_Ref(); - /*--- Check for back flow through the inlet. ---*/ + /*--- Store the current static pressure for clarity. ---*/ - Vn = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - Vn += V_domain[iDim+1]*(-1.0*Normal[iDim]/Area); - } + P_domain = nodes->GetPressure(iPoint); - /*--- If the local static pressure is larger than the specified - total pressure or the velocity is directed upstream, we have a - back flow situation. The specified total pressure should be used - as a static pressure condition and the velocity from the domain - is used for the BC. ---*/ + /*--- Check for back flow through the inlet. ---*/ - if ((P_domain > P_total) || (Vn < 0.0)) { + Vn = 0.0; + for (iDim = 0; iDim < nDim; iDim++) { + Vn += V_domain[iDim+1]*(-1.0*Normal[iDim]/Area); + } - /*--- Back flow: use the prescribed P_total as static pressure. ---*/ + /*--- If the local static pressure is larger than the specified + total pressure or the velocity is directed upstream, we have a + back flow situation. The specified total pressure should be used + as a static pressure condition and the velocity from the domain + is used for the BC. ---*/ - V_inlet[0] = Inlet_Ptotal[val_marker][iVertex]/config->GetPressure_Ref(); + if ((P_domain > P_total) || (Vn < 0.0)) { - /*--- Neumann condition for velocity. ---*/ + /*--- Back flow: use the prescribed P_total as static pressure. ---*/ - for (iDim = 0; iDim < nDim; iDim++) - V_inlet[iDim+1] = V_domain[iDim+1]; + V_inlet[0] = Inlet_Ptotal[val_marker][iVertex]/config->GetPressure_Ref(); - /*--- Neumann condition for the temperature. ---*/ + /*--- Neumann condition for velocity. ---*/ - V_inlet[nDim+1] = nodes->GetTemperature(iPoint); + for (iDim = 0; iDim < nDim; iDim++) + V_inlet[iDim+1] = V_domain[iDim+1]; - } else { + /*--- Neumann condition for the temperature. ---*/ - /*--- Update the velocity magnitude using the total pressure. ---*/ + V_inlet[nDim+1] = nodes->GetTemperature(iPoint); - Vel_Mag = sqrt((P_total - P_domain)/(0.5*nodes->GetDensity(iPoint))); + } else { - /*--- If requested, use the local boundary normal (negative), - instead of the prescribed flow direction in the config. ---*/ + /*--- Update the velocity magnitude using the total pressure. ---*/ - if (config->GetInc_Inlet_UseNormal()) { - for (iDim = 0; iDim < nDim; iDim++) - UnitFlowDir[iDim] = -Normal[iDim]/Area; - } + Vel_Mag = sqrt((P_total - P_domain)/(0.5*nodes->GetDensity(iPoint))); - /*--- Compute the delta change in velocity in each direction. ---*/ + /*--- If requested, use the local boundary normal (negative), + instead of the prescribed flow direction in the config. ---*/ + if (config->GetInc_Inlet_UseNormal()) { for (iDim = 0; iDim < nDim; iDim++) - dV[iDim] = Vel_Mag*UnitFlowDir[iDim] - V_domain[iDim+1]; + UnitFlowDir[iDim] = -Normal[iDim]/Area; + } - /*--- Update the velocity in the primitive variable vector. - Note we use damping here to improve stability/convergence. ---*/ + /*--- Compute the delta change in velocity in each direction. ---*/ - for (iDim = 0; iDim < nDim; iDim++) - V_inlet[iDim+1] = V_domain[iDim+1] + Damping*dV[iDim]; + for (iDim = 0; iDim < nDim; iDim++) + dV[iDim] = Vel_Mag*UnitFlowDir[iDim] - V_domain[iDim+1]; - /*--- Dirichlet condition for temperature (if energy is active) ---*/ + /*--- Update the velocity in the primitive variable vector. + Note we use damping here to improve stability/convergence. ---*/ - V_inlet[nDim+1] = Inlet_Ttotal[val_marker][iVertex]/config->GetTemperature_Ref(); + for (iDim = 0; iDim < nDim; iDim++) + V_inlet[iDim+1] = V_domain[iDim+1] + Damping*dV[iDim]; - } + /*--- Dirichlet condition for temperature (if energy is active) ---*/ - break; + V_inlet[nDim+1] = Inlet_Ttotal[val_marker][iVertex]/config->GetTemperature_Ref(); - } + } - /*--- Access density at the node. This is either constant by - construction, or will be set fixed implicitly by the temperature - and equation of state. ---*/ + break; - V_inlet[nDim+2] = nodes->GetDensity(iPoint); + } - /*--- Beta coefficient from the config file ---*/ + /*--- Access density at the node. This is either constant by + construction, or will be set fixed implicitly by the temperature + and equation of state. ---*/ - V_inlet[nDim+3] = nodes->GetBetaInc2(iPoint); + V_inlet[nDim+2] = nodes->GetDensity(iPoint); - /*--- Cp is needed for Temperature equation. ---*/ + /*--- Beta coefficient from the config file ---*/ - V_inlet[nDim+7] = nodes->GetSpecificHeatCp(iPoint); + V_inlet[nDim+3] = nodes->GetBetaInc2(iPoint); - /*--- Set various quantities in the solver class ---*/ + /*--- Cp is needed for Temperature equation. ---*/ - conv_numerics->SetPrimitive(V_domain, V_inlet); + V_inlet[nDim+7] = nodes->GetSpecificHeatCp(iPoint); - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + /*--- Set various quantities in the solver class ---*/ - /*--- Compute the residual using an upwind scheme ---*/ + conv_numerics->SetPrimitive(V_domain, V_inlet); - auto residual = conv_numerics->ComputeResidual(config); + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); - /*--- Update residual value ---*/ + /*--- Compute the residual using an upwind scheme ---*/ - LinSysRes.AddBlock(iPoint, residual); + auto residual = conv_numerics->ComputeResidual(config); - /*--- Jacobian contribution for implicit integration ---*/ + /*--- Update residual value ---*/ - if (implicit) - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + LinSysRes.AddBlock(iPoint, residual); - /*--- Viscous contribution, commented out because serious convergence problems ---*/ + /*--- Jacobian contribution for implicit integration ---*/ - if (viscous) { + if (implicit) + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - /*--- Set transport properties at the inlet ---*/ + /*--- Viscous contribution, commented out because serious convergence problems ---*/ - V_inlet[nDim+4] = nodes->GetLaminarViscosity(iPoint); - V_inlet[nDim+5] = nodes->GetEddyViscosity(iPoint); - V_inlet[nDim+6] = nodes->GetThermalConductivity(iPoint); + if (!viscous) continue; - /*--- Set the normal vector and the coordinates ---*/ + /*--- Set transport properties at the inlet ---*/ - visc_numerics->SetNormal(Normal); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - geometry->nodes->GetCoord(Point_Normal)); + V_inlet[nDim+4] = nodes->GetLaminarViscosity(iPoint); + V_inlet[nDim+5] = nodes->GetEddyViscosity(iPoint); + V_inlet[nDim+6] = nodes->GetThermalConductivity(iPoint); - /*--- Primitive variables, and gradient ---*/ + /*--- Set the normal vector and the coordinates ---*/ - visc_numerics->SetPrimitive(V_domain, V_inlet); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - nodes->GetGradient_Primitive(iPoint)); + visc_numerics->SetNormal(Normal); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(Point_Normal)); - /*--- Turbulent kinetic energy ---*/ + /*--- Primitive variables, and gradient ---*/ - if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) - visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), - solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); + visc_numerics->SetPrimitive(V_domain, V_inlet); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), + nodes->GetGradient_Primitive(iPoint)); - /*--- Compute and update residual ---*/ + /*--- Turbulent kinetic energy ---*/ - auto residual = visc_numerics->ComputeResidual(config); + if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) + visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), + solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); - LinSysRes.SubtractBlock(iPoint, residual); + /*--- Compute and update residual ---*/ - /*--- Jacobian contribution for implicit integration ---*/ + auto residual_v = visc_numerics->ComputeResidual(config); - if (implicit) - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + LinSysRes.SubtractBlock(iPoint, residual_v); - } + /*--- Jacobian contribution for implicit integration ---*/ - } + if (implicit) + Jacobian.SubtractBlock2Diag(iPoint, residual_v.jacobian_i); } - - /*--- Free locally allocated memory ---*/ - - delete [] Normal; - } void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, @@ -2725,198 +2222,191 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, su2double mDot_Target, mDot_Old, dP, Density_Avg, Area_Outlet; su2double Damping = config->GetInc_Outlet_Damping(); - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool viscous = config->GetViscous(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool viscous = config->GetViscous(); string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - su2double *Normal = new su2double[nDim]; + su2double Normal[MAXNDIM] = {0.0}; unsigned short Kind_Outlet = config->GetKind_Inc_Outlet(Marker_Tag); /*--- Loop over all the vertices on this boundary marker ---*/ for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - /*--- Allocate the value at the outlet ---*/ - - V_outlet = GetCharacPrimVar(val_marker, iVertex); - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e., not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { + if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Index of the closest interior node ---*/ + /*--- Allocate the value at the outlet ---*/ - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + V_outlet = GetCharacPrimVar(val_marker, iVertex); - /*--- Normal vector for this vertex (negate for outward convention) ---*/ + /*--- Index of the closest interior node ---*/ - geometry->vertex[val_marker][iVertex]->GetNormal(Normal); - for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; - conv_numerics->SetNormal(Normal); + Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - /*--- Current solution at this boundary node ---*/ + /*--- Normal vector for this vertex (negate for outward convention) ---*/ - V_domain = nodes->GetPrimitive(iPoint); + geometry->vertex[val_marker][iVertex]->GetNormal(Normal); + for (iDim = 0; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim]; + conv_numerics->SetNormal(Normal); - /*--- Store the current static pressure for clarity. ---*/ + /*--- Current solution at this boundary node ---*/ - P_domain = nodes->GetPressure(iPoint); + V_domain = nodes->GetPrimitive(iPoint); - /*--- Compute a boundary value for the pressure depending on whether - we are prescribing a back pressure or a mass flow target. ---*/ + /*--- Store the current static pressure for clarity. ---*/ - switch (Kind_Outlet) { + P_domain = nodes->GetPressure(iPoint); - /*--- Velocity and temperature (if required) been specified at the inlet. ---*/ + /*--- Compute a boundary value for the pressure depending on whether + we are prescribing a back pressure or a mass flow target. ---*/ - case PRESSURE_OUTLET: + switch (Kind_Outlet) { - /*--- Retrieve the specified back pressure for this outlet. ---*/ + /*--- Velocity and temperature (if required) been specified at the inlet. ---*/ - P_Outlet = config->GetOutlet_Pressure(Marker_Tag)/config->GetPressure_Ref(); + case PRESSURE_OUTLET: - /*--- The pressure is prescribed at the outlet. ---*/ + /*--- Retrieve the specified back pressure for this outlet. ---*/ - V_outlet[0] = P_Outlet; + P_Outlet = config->GetOutlet_Pressure(Marker_Tag)/config->GetPressure_Ref(); - /*--- Neumann condition for the velocity. ---*/ + /*--- The pressure is prescribed at the outlet. ---*/ - for (iDim = 0; iDim < nDim; iDim++) { - V_outlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); - } + V_outlet[0] = P_Outlet; + + /*--- Neumann condition for the velocity. ---*/ - break; + for (iDim = 0; iDim < nDim; iDim++) { + V_outlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); + } - /*--- A mass flow target has been specified for the outlet. ---*/ + break; - case MASS_FLOW_OUTLET: + /*--- A mass flow target has been specified for the outlet. ---*/ - /*--- Retrieve the specified target mass flow at the outlet. ---*/ + case MASS_FLOW_OUTLET: - mDot_Target = config->GetOutlet_Pressure(Marker_Tag)/(config->GetDensity_Ref() * config->GetVelocity_Ref()); + /*--- Retrieve the specified target mass flow at the outlet. ---*/ - /*--- Retrieve the old mass flow, density, and area of the outlet, - which has been computed in a preprocessing step. These values - were stored in non-dim. form in the config container. ---*/ + mDot_Target = config->GetOutlet_Pressure(Marker_Tag)/(config->GetDensity_Ref() * config->GetVelocity_Ref()); - mDot_Old = config->GetOutlet_MassFlow(Marker_Tag); - Density_Avg = config->GetOutlet_Density(Marker_Tag); - Area_Outlet = config->GetOutlet_Area(Marker_Tag); + /*--- Retrieve the old mass flow, density, and area of the outlet, + which has been computed in a preprocessing step. These values + were stored in non-dim. form in the config container. ---*/ - /*--- Compute the pressure increment based on the difference - between the current and target mass flow. Note that increasing - pressure decreases flow speed. ---*/ + mDot_Old = config->GetOutlet_MassFlow(Marker_Tag); + Density_Avg = config->GetOutlet_Density(Marker_Tag); + Area_Outlet = config->GetOutlet_Area(Marker_Tag); - dP = 0.5*Density_Avg*(mDot_Old*mDot_Old - mDot_Target*mDot_Target)/((Density_Avg*Area_Outlet)*(Density_Avg*Area_Outlet)); + /*--- Compute the pressure increment based on the difference + between the current and target mass flow. Note that increasing + pressure decreases flow speed. ---*/ - /*--- Update the new outlet pressure. Note that we use damping - here to improve stability/convergence. ---*/ + dP = 0.5*Density_Avg*(mDot_Old*mDot_Old - mDot_Target*mDot_Target)/((Density_Avg*Area_Outlet)*(Density_Avg*Area_Outlet)); - P_Outlet = P_domain + Damping*dP; + /*--- Update the new outlet pressure. Note that we use damping + here to improve stability/convergence. ---*/ - /*--- The pressure is prescribed at the outlet. ---*/ + P_Outlet = P_domain + Damping*dP; - V_outlet[0] = P_Outlet; + /*--- The pressure is prescribed at the outlet. ---*/ - /*--- Neumann condition for the velocity ---*/ + V_outlet[0] = P_Outlet; - for (iDim = 0; iDim < nDim; iDim++) { - V_outlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); - } + /*--- Neumann condition for the velocity ---*/ - break; + for (iDim = 0; iDim < nDim; iDim++) { + V_outlet[iDim+1] = nodes->GetVelocity(iPoint,iDim); + } - } + break; - /*--- Neumann condition for the temperature. ---*/ + } - V_outlet[nDim+1] = nodes->GetTemperature(iPoint); + /*--- Neumann condition for the temperature. ---*/ - /*--- Access density at the interior node. This is either constant by - construction, or will be set fixed implicitly by the temperature - and equation of state. ---*/ + V_outlet[nDim+1] = nodes->GetTemperature(iPoint); - V_outlet[nDim+2] = nodes->GetDensity(iPoint); + /*--- Access density at the interior node. This is either constant by + construction, or will be set fixed implicitly by the temperature + and equation of state. ---*/ - /*--- Beta coefficient from the config file ---*/ + V_outlet[nDim+2] = nodes->GetDensity(iPoint); - V_outlet[nDim+3] = nodes->GetBetaInc2(iPoint); + /*--- Beta coefficient from the config file ---*/ - /*--- Cp is needed for Temperature equation. ---*/ + V_outlet[nDim+3] = nodes->GetBetaInc2(iPoint); - V_outlet[nDim+7] = nodes->GetSpecificHeatCp(iPoint); + /*--- Cp is needed for Temperature equation. ---*/ - /*--- Set various quantities in the solver class ---*/ + V_outlet[nDim+7] = nodes->GetSpecificHeatCp(iPoint); - conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- Set various quantities in the solver class ---*/ - if (dynamic_grid) - conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), - geometry->nodes->GetGridVel(iPoint)); + conv_numerics->SetPrimitive(V_domain, V_outlet); - /*--- Compute the residual using an upwind scheme ---*/ + if (dynamic_grid) + conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), + geometry->nodes->GetGridVel(iPoint)); - auto residual = conv_numerics->ComputeResidual(config); + /*--- Compute the residual using an upwind scheme ---*/ - /*--- Update residual value ---*/ + auto residual = conv_numerics->ComputeResidual(config); - LinSysRes.AddBlock(iPoint, residual); + /*--- Update residual value ---*/ - /*--- Jacobian contribution for implicit integration ---*/ + LinSysRes.AddBlock(iPoint, residual); - if (implicit) { - Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); - } + /*--- Jacobian contribution for implicit integration ---*/ - /*--- Viscous contribution, commented out because serious convergence problems ---*/ + if (implicit) { + Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i); + } - if (viscous) { + /*--- Viscous contribution, commented out because serious convergence problems ---*/ - /*--- Set transport properties at the outlet. ---*/ + if (!viscous) continue; - V_outlet[nDim+4] = nodes->GetLaminarViscosity(iPoint); - V_outlet[nDim+5] = nodes->GetEddyViscosity(iPoint); - V_outlet[nDim+6] = nodes->GetThermalConductivity(iPoint); + /*--- Set transport properties at the outlet. ---*/ - /*--- Set the normal vector and the coordinates ---*/ + V_outlet[nDim+4] = nodes->GetLaminarViscosity(iPoint); + V_outlet[nDim+5] = nodes->GetEddyViscosity(iPoint); + V_outlet[nDim+6] = nodes->GetThermalConductivity(iPoint); - visc_numerics->SetNormal(Normal); - visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), - geometry->nodes->GetCoord(Point_Normal)); + /*--- Set the normal vector and the coordinates ---*/ - /*--- Primitive variables, and gradient ---*/ + visc_numerics->SetNormal(Normal); + visc_numerics->SetCoord(geometry->nodes->GetCoord(iPoint), + geometry->nodes->GetCoord(Point_Normal)); - visc_numerics->SetPrimitive(V_domain, V_outlet); - visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), - nodes->GetGradient_Primitive(iPoint)); + /*--- Primitive variables, and gradient ---*/ - /*--- Turbulent kinetic energy ---*/ + visc_numerics->SetPrimitive(V_domain, V_outlet); + visc_numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), + nodes->GetGradient_Primitive(iPoint)); - if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) - visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), - solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); + /*--- Turbulent kinetic energy ---*/ - /*--- Compute and update residual ---*/ + if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) + visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), + solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); - auto residual = visc_numerics->ComputeResidual(config); + /*--- Compute and update residual ---*/ - LinSysRes.SubtractBlock(iPoint, residual); + auto residual_v = visc_numerics->ComputeResidual(config); - /*--- Jacobian contribution for implicit integration ---*/ - if (implicit) - Jacobian.SubtractBlock2Diag(iPoint, residual.jacobian_i); + LinSysRes.SubtractBlock(iPoint, residual_v); - } + /*--- Jacobian contribution for implicit integration ---*/ + if (implicit) + Jacobian.SubtractBlock2Diag(iPoint, residual_v.jacobian_i); - } } - /*--- Free locally allocated memory ---*/ - delete [] Normal; - } void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_container, CConfig *config, @@ -2924,18 +2414,26 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Local variables ---*/ - unsigned short iVar, jVar, iMarker, iDim; + unsigned short iVar, iMarker, iDim, iNeigh; unsigned long iPoint, jPoint, iEdge, iVertex; - su2double Density, Cp; - su2double *V_time_nM1, *V_time_n, *V_time_nP1; - su2double U_time_nM1[5], U_time_n[5], U_time_nP1[5]; + const su2double *V_time_nM1 = nullptr, *V_time_n = nullptr, *V_time_nP1 = nullptr; + su2double U_time_nM1[MAXNVAR], U_time_n[MAXNVAR], U_time_nP1[MAXNVAR]; su2double Volume_nM1, Volume_nP1, TimeStep; - su2double *GridVel_i = nullptr, *GridVel_j = nullptr, Residual_GCL; - const su2double* Normal; + const su2double *Normal = nullptr, *GridVel_i = nullptr, *GridVel_j = nullptr; + su2double Density, Cp; - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool energy = config->GetEnergy_Equation(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool first_order = (config->GetTime_Marching() == DT_STEPPING_1ST); + const bool second_order = (config->GetTime_Marching() == DT_STEPPING_2ND); + const bool energy = config->GetEnergy_Equation(); + + const int ndim = nDim; + auto V2U = [ndim](su2double Density, su2double Cp, const su2double* V, su2double* U) { + U[0] = Density; + for (int iDim = 0; iDim < ndim; iDim++) U[iDim+1] = Density*V[iDim+1]; + U[ndim+1] = Density*Cp*V[ndim+1]; + }; /*--- Store the physical time step ---*/ @@ -2947,18 +2445,9 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Loop over all nodes (excluding halos) ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - /*--- Initialize the Residual / Jacobian container to zero. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Residual[iVar] = 0.0; - if (implicit) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - } - /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, previous solutions that are stored in memory. These are actually @@ -2970,24 +2459,14 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Access the density and Cp at this node (constant for now). ---*/ - Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); + Density = nodes->GetDensity(iPoint); + Cp = nodes->GetSpecificHeatCp(iPoint); /*--- Compute the conservative variable vector for all time levels. ---*/ - U_time_nM1[0] = Density; - U_time_n[0] = Density; - U_time_nP1[0] = Density; - - for (iDim = 0; iDim < nDim; iDim++) { - U_time_nM1[iDim+1] = Density*V_time_nM1[iDim+1]; - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - U_time_nP1[iDim+1] = Density*V_time_nP1[iDim+1]; - } - - U_time_nM1[nDim+1] = Density*Cp*V_time_nM1[nDim+1]; - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - U_time_nP1[nDim+1] = Density*Cp*V_time_nP1[nDim+1]; + V2U(Density, Cp, V_time_nM1, U_time_nM1); + V2U(Density, Cp, V_time_n, U_time_n); + V2U(Density, Cp, V_time_nP1, U_time_nP1); /*--- CV volume at time n+1. As we are on a static mesh, the volume of the CV will remained fixed for all time steps. ---*/ @@ -2995,41 +2474,28 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver Volume_nP1 = geometry->nodes->GetVolume(iPoint); /*--- Compute the dual time-stepping source term based on the chosen - time discretization scheme (1st- or 2nd-order). Note that for an - incompressible problem, the pressure equation does not have a - contribution, as the time derivative should always be zero. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*Volume_nP1 / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Residual[iVar] = ( 3.0*U_time_nP1[iVar] - 4.0*U_time_n[iVar] - +1.0*U_time_nM1[iVar])*Volume_nP1 / (2.0*TimeStep); - } - - if (!energy) Residual[nDim+1] = 0.0; + time discretization scheme (1st- or 2nd-order).---*/ - /*--- Store the residual and compute the Jacobian contribution due - to the dual time source term. ---*/ + for (iVar = 0; iVar < nVar-!energy; iVar++) { + if (first_order) + LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*Volume_nP1 / TimeStep; + if (second_order) + LinSysRes(iPoint,iVar) += ( 3.0*U_time_nP1[iVar] - 4.0*U_time_n[iVar] + +1.0*U_time_nM1[iVar])*Volume_nP1 / (2.0*TimeStep); + } - LinSysRes.AddBlock(iPoint, Residual); + /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ if (implicit) { - for (iVar = 1; iVar < nVar; iVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep); - } + su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; + for (iDim = 0; iDim < nDim; iDim++) - Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1]; - if (energy) Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1]; + Jacobian.AddVal2Diag(iPoint, iDim+1, delta); - Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + if (energy) delta *= Cp; + Jacobian.AddVal2Diag(iPoint, nDim+1, delta); } - } - } else { @@ -3042,141 +2508,86 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver we will loop over the edges and boundaries to compute the GCL component of the dual time source term that depends on grid velocities. ---*/ - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Initialize the Residual / Jacobian container to zero. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - - /*--- Get indices for nodes i & j plus the face normal ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) + for (iPoint = 0; iPoint < nPointDomain; ++iPoint) { - iPoint = geometry->edges->GetNode(iEdge,0); - jPoint = geometry->edges->GetNode(iEdge,1); - Normal = geometry->edges->GetNormal(iEdge); - - /*--- Grid velocities stored at nodes i & j ---*/ - - GridVel_i = geometry->nodes->GetGridVel(iPoint); - GridVel_j = geometry->nodes->GetGridVel(jPoint); - - /*--- Compute the GCL term by averaging the grid velocities at the - edge mid-point and dotting with the face normal. ---*/ - - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL += 0.5*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - - /*--- Compute the GCL component of the source term for node i ---*/ + /*--- Compute the conservative variables. ---*/ V_time_n = nodes->GetSolution_time_n(iPoint); + Density = nodes->GetDensity(iPoint); + Cp = nodes->GetSpecificHeatCp(iPoint); + V2U(Density, Cp, V_time_n, U_time_n); - /*--- Access the density and Cp at this node (constant for now). ---*/ - - Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); + GridVel_i = geometry->nodes->GetGridVel(iPoint); - /*--- Compute the conservative variable vector for all time levels. ---*/ + for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); iNeigh++) { - U_time_n[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - } - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = U_time_n[iVar]*Residual_GCL; + iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); + Normal = geometry->edges->GetNormal(iEdge); - if (!energy) Residual[nDim+1] = 0.0; - LinSysRes.AddBlock(iPoint, Residual); + jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + GridVel_j = geometry->nodes->GetGridVel(jPoint); - /*--- Compute the GCL component of the source term for node j ---*/ + /*--- Determine whether to consider the normal outward or inward. ---*/ + su2double dir = (geometry->edges->GetNode(iEdge,0) == iPoint)? 0.5 : -0.5; - V_time_n = nodes->GetSolution_time_n(jPoint); + su2double Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL += dir*(GridVel_i[iDim]+GridVel_j[iDim])*Normal[iDim]; - U_time_n[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + for (iVar = 0; iVar < nVar-!energy; iVar++) + LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; } - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = U_time_n[iVar]*Residual_GCL; - - if (!energy) Residual[nDim+1] = 0.0; - LinSysRes.SubtractBlock(jPoint, Residual); - } - /*--- Loop over the boundary edges ---*/ + /*--- Loop over the boundary edges ---*/ for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && + if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Initialize the Residual / Jacobian container to zero. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - - /*--- Get the index for node i plus the boundary face normal ---*/ - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - /*--- Grid velocities stored at boundary node i ---*/ + /*--- Get the index for node i plus the boundary face normal ---*/ - GridVel_i = geometry->nodes->GetGridVel(iPoint); + iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - /*--- Compute the GCL term by dotting the grid velocity with the face - normal. The normal is negated to match the boundary convention. ---*/ + /*--- Grid velocities stored at boundary node i ---*/ - Residual_GCL = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; + GridVel_i = geometry->nodes->GetGridVel(iPoint); - /*--- Compute the GCL component of the source term for node i ---*/ + /*--- Compute the GCL term by dotting the grid velocity with the face + normal. The normal is negated to match the boundary convention. ---*/ - V_time_n = nodes->GetSolution_time_n(iPoint); + su2double Residual_GCL = 0.0; + for (iDim = 0; iDim < nDim; iDim++) + Residual_GCL -= 0.5*(GridVel_i[iDim]+GridVel_i[iDim])*Normal[iDim]; - /*--- Access the density and Cp at this node (constant for now). ---*/ + /*--- Compute the GCL component of the source term for node i ---*/ - Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); + V_time_n = nodes->GetSolution_time_n(iPoint); + Density = nodes->GetDensity(iPoint); + Cp = nodes->GetSpecificHeatCp(iPoint); + V2U(Density, Cp, V_time_n, U_time_n); - U_time_n[0] = Density; - for (iDim = 0; iDim < nDim; iDim++) { - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; + for (iVar = 0; iVar < nVar-!energy; iVar++) + LinSysRes(iPoint,iVar) += U_time_n[iVar]*Residual_GCL; } - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - - for (iVar = 0; iVar < nVar; iVar++) - Residual[iVar] = U_time_n[iVar]*Residual_GCL; - - if (!energy) Residual[nDim+1] = 0.0; - LinSysRes.AddBlock(iPoint, Residual); - - } } } /*--- Loop over all nodes (excluding halos) to compute the remainder of the dual time-stepping source term. ---*/ + SU2_OMP_FOR_STAT(omp_chunk_size) for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - /*--- Initialize the Residual / Jacobian container to zero. ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Residual[iVar] = 0.0; - if (implicit) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - } - /*--- Retrieve the solution at time levels n-1, n, and n+1. Note that we are currently iterating on U^n+1 and that U^n & U^n-1 are fixed, - previous solutions that are stored in memory. ---*/ + previous solutions that are stored in memory. These are actually + the primitive values, but we will convert to conservatives. ---*/ V_time_nM1 = nodes->GetSolution_time_n1(iPoint); V_time_n = nodes->GetSolution_time_n(iPoint); @@ -3184,24 +2595,14 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver /*--- Access the density and Cp at this node (constant for now). ---*/ - Density = nodes->GetDensity(iPoint); - Cp = nodes->GetSpecificHeatCp(iPoint); + Density = nodes->GetDensity(iPoint); + Cp = nodes->GetSpecificHeatCp(iPoint); /*--- Compute the conservative variable vector for all time levels. ---*/ - U_time_nM1[0] = Density; - U_time_n[0] = Density; - U_time_nP1[0] = Density; - - for (iDim = 0; iDim < nDim; iDim++) { - U_time_nM1[iDim+1] = Density*V_time_nM1[iDim+1]; - U_time_n[iDim+1] = Density*V_time_n[iDim+1]; - U_time_nP1[iDim+1] = Density*V_time_nP1[iDim+1]; - } - - U_time_nM1[nDim+1] = Density*Cp*V_time_nM1[nDim+1]; - U_time_n[nDim+1] = Density*Cp*V_time_n[nDim+1]; - U_time_nP1[nDim+1] = Density*Cp*V_time_nP1[nDim+1]; + V2U(Density, Cp, V_time_nM1, U_time_nM1); + V2U(Density, Cp, V_time_n, U_time_n); + V2U(Density, Cp, 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 @@ -3214,33 +2615,25 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver introduction of the GCL term above, the remainder of the source residual due to the time discretization has a new form.---*/ - for (iVar = 0; iVar < nVar; iVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Residual[iVar] = (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) - + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); + for (iVar = 0; iVar < nVar-!energy; iVar++) { + if (first_order) + LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*(Volume_nP1/TimeStep); + if (second_order) + LinSysRes(iPoint,iVar) += (U_time_nP1[iVar] - U_time_n[iVar])*(3.0*Volume_nP1/(2.0*TimeStep)) + + (U_time_nM1[iVar] - U_time_n[iVar])*(Volume_nM1/(2.0*TimeStep)); } - /*--- Store the residual and compute the Jacobian contribution due - to the dual time source term. ---*/ - if (!energy) Residual[nDim+1] = 0.0; - LinSysRes.AddBlock(iPoint, Residual); + /*--- Compute the Jacobian contribution due to the dual time source term. ---*/ if (implicit) { - for (iVar = 1; iVar < nVar; iVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep); - } + su2double delta = (second_order? 1.5 : 1.0) * Volume_nP1 * Density / TimeStep; + for (iDim = 0; iDim < nDim; iDim++) - Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1]; - if (energy) Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1]; + Jacobian.AddVal2Diag(iPoint, iDim+1, delta); - Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + if (energy) delta *= Cp; + Jacobian.AddVal2Diag(iPoint, nDim+1, delta); } - } } @@ -3254,6 +2647,7 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, Velocity2, Density, Area, AxiFactor; unsigned short iMarker_Outlet, nMarker_Outlet; string Inlet_TagBound, Outlet_TagBound; + su2double Vector[MAXNDIM] = {0.0}; bool axisymmetric = config->GetAxisymmetric(); @@ -3372,22 +2766,10 @@ void CIncEulerSolver::GetOutlet_Properties(CGeometry *geometry, CConfig *config, /*--- All the ranks to compute the total value ---*/ -#ifdef HAVE_MPI - SU2_MPI::Allreduce(Outlet_MassFlow_Local, Outlet_MassFlow_Total, nMarker_Outlet, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); SU2_MPI::Allreduce(Outlet_Density_Local, Outlet_Density_Total, nMarker_Outlet, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); SU2_MPI::Allreduce(Outlet_Area_Local, Outlet_Area_Total, nMarker_Outlet, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -#else - - for (iMarker_Outlet = 0; iMarker_Outlet < nMarker_Outlet; iMarker_Outlet++) { - Outlet_MassFlow_Total[iMarker_Outlet] = Outlet_MassFlow_Local[iMarker_Outlet]; - Outlet_Density_Total[iMarker_Outlet] = Outlet_Density_Local[iMarker_Outlet]; - Outlet_Area_Total[iMarker_Outlet] = Outlet_Area_Local[iMarker_Outlet]; - } - -#endif - for (iMarker_Outlet = 0; iMarker_Outlet < nMarker_Outlet; iMarker_Outlet++) { if (Outlet_Area_Total[iMarker_Outlet] != 0.0) { Outlet_Density_Total[iMarker_Outlet] /= Outlet_Area_Total[iMarker_Outlet]; @@ -3506,7 +2888,7 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf unsigned short iDim, iVar, iMesh, iMeshFine; unsigned long iPoint, index, iChildren, Point_Fine; unsigned short turb_model = config->GetKind_Turb_Model(); - su2double Area_Children, Area_Parent, Coord[3] = {0.0}, *Solution_Fine; + su2double Area_Children, Area_Parent, Coord[MAXNDIM] = {0.0}, *Solution_Fine; bool static_fsi = ((config->GetTime_Marching() == STEADY) && config->GetFSI_Simulation()); bool dual_time = ((config->GetTime_Marching() == DT_STEPPING_1ST) || (config->GetTime_Marching() == DT_STEPPING_2ND)); @@ -3541,6 +2923,7 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf unsigned short nVar_Restart = nVar; if ((!energy) && (!weakly_coupled_heat)) nVar_Restart--; + su2double Solution[MAXNVAR] = {0.0}; Solution[nVar-1] = GetTemperature_Inf(); /*--- Read the restart data from either an ASCII or binary SU2 file. ---*/ @@ -3584,7 +2967,7 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf index = counter*Restart_Vars[1]; for (iDim = 0; iDim < nDim; iDim++) { Coord[iDim] = Restart_Data[index+iDim]; } - su2double GridVel[3] = {0.0,0.0,0.0}; + su2double GridVel[MAXNDIM] = {0.0}; if (!steady_restart) { /*--- Move the index forward to get the grid velocities. ---*/ index = counter*Restart_Vars[1] + skipVars + nVar_Restart + turbVars; @@ -3710,7 +3093,7 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf } -void CIncEulerSolver::SetFreeStream_Solution(CConfig *config){ +void CIncEulerSolver::SetFreeStream_Solution(const CConfig *config){ unsigned long iPoint; unsigned short iDim; diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 976da459e9e3..5dd6b564bbaa 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -54,120 +54,60 @@ CIncNSSolver::CIncNSSolver(CGeometry *geometry, CConfig *config, unsigned short } } -void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - - unsigned long iPoint, ErrorCounter = 0; - su2double StrainMag = 0.0, Omega = 0.0, *Vorticity; - - unsigned long InnerIter = config->GetInnerIter(); - bool cont_adjoint = config->GetContinuous_Adjoint(); - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool center = ((config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED)); - bool center_jst = center && config->GetKind_Centered_Flow() == JST; - bool limiter_flow = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); - bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); - bool limiter_adjflow = (cont_adjoint && (config->GetKind_SlopeLimit_AdjFlow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter())); - bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; - bool outlet = ((config->GetnMarker_Outlet() != 0)); - - /*--- Set the primitive variables ---*/ - - ErrorCounter = SetPrimitive_Variables(solver_container, config, Output); - - /*--- Compute gradient for MUSCL reconstruction. ---*/ - - if (config->GetReconstructionGradientRequired() && (iMesh == MESH_0)) { - if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) - SetPrimitive_Gradient_GG(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); - if (config->GetKind_Gradient_Method_Recon() == WEIGHTED_LEAST_SQUARES) - SetPrimitive_Gradient_LS(geometry, config, true); - } +void CIncNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - /*--- Compute gradient of the primitive variables ---*/ + const auto InnerIter = config->GetInnerIter(); + const bool muscl = config->GetMUSCL_Flow() && (iMesh == MESH_0); + const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); + const bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); + const bool van_albada = (config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE); - if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { - SetPrimitive_Gradient_GG(geometry, config); - } - if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { - SetPrimitive_Gradient_LS(geometry, config); - } + /*--- Common preprocessing steps (implemented by CEulerSolver) ---*/ - /*--- Compute the limiter in case we need it in the turbulence model - or to limit the viscous terms (check this logic with JST and 2nd order turbulence model) ---*/ + CommonPreprocessing(geometry, solver_container, config, iMesh, iRKStep, RunTime_EqSystem, Output); - if ((iMesh == MESH_0) && (limiter_flow || limiter_turb || limiter_adjflow) - && !Output && !van_albada) { SetPrimitive_Limiter(geometry, config); } + /*--- Compute gradient for MUSCL reconstruction ---*/ - /*--- Artificial dissipation for centered schemes. ---*/ - - if (center && !Output) { - SetMax_Eigenvalue(geometry, config); - if ((center_jst) && (iMesh == MESH_0)) { - SetCentered_Dissipation_Sensor(geometry, config); - SetUndivided_Laplacian(geometry, config); + if (config->GetReconstructionGradientRequired() && muscl && !center) { + switch (config->GetKind_Gradient_Method_Recon()) { + case GREEN_GAUSS: + SetPrimitive_Gradient_GG(geometry, config, true); break; + case LEAST_SQUARES: + case WEIGHTED_LEAST_SQUARES: + SetPrimitive_Gradient_LS(geometry, config, true); break; + default: break; } } - /*--- Update the beta value based on the maximum velocity / viscosity. ---*/ - - SetBeta_Parameter(geometry, solver_container, config, iMesh); - - /*--- Compute properties needed for mass flow BCs. ---*/ - - if (outlet) GetOutlet_Properties(geometry, config, iMesh, Output); - - /*--- Evaluate the vorticity and strain rate magnitude ---*/ - - nodes->SetVorticity_StrainMag(); - - StrainMag_Max = 0.0; Omega_Max = 0.0; - for (iPoint = 0; iPoint < nPoint; iPoint++) { - - StrainMag = nodes->GetStrainMag(iPoint); - Vorticity = nodes->GetVorticity(iPoint); - Omega = sqrt(Vorticity[0]*Vorticity[0]+ Vorticity[1]*Vorticity[1]+ Vorticity[2]*Vorticity[2]); - - StrainMag_Max = max(StrainMag_Max, StrainMag); - Omega_Max = max(Omega_Max, Omega); + /*--- Compute gradient of the primitive variables ---*/ + if (config->GetKind_Gradient_Method() == GREEN_GAUSS) { + SetPrimitive_Gradient_GG(geometry, config); + } + else if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) { + SetPrimitive_Gradient_LS(geometry, config); } - /*--- Initialize the Jacobian matrices ---*/ - - if (implicit && !Output) Jacobian.SetValZero(); - - /*--- Error message ---*/ - - if (config->GetComm_Level() == COMM_FULL) { - -#ifdef HAVE_MPI - unsigned long MyErrorCounter = ErrorCounter; ErrorCounter = 0; - su2double MyOmega_Max = Omega_Max; Omega_Max = 0.0; - su2double MyStrainMag_Max = StrainMag_Max; StrainMag_Max = 0.0; - - SU2_MPI::Allreduce(&MyErrorCounter, &ErrorCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD); - SU2_MPI::Allreduce(&MyStrainMag_Max, &StrainMag_Max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - SU2_MPI::Allreduce(&MyOmega_Max, &Omega_Max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); -#endif - - if (iMesh == MESH_0) - config->SetNonphysical_Points(ErrorCounter); + /*--- Compute the limiters ---*/ + if (muscl && !center && limiter && !van_albada && !Output) { + SetPrimitive_Limiter(geometry, config); } + ComputeVorticityAndStrainMag(*config, iMesh); + } -unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { +unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, const CConfig *config) { unsigned long iPoint, nonPhysicalPoints = 0; su2double eddy_visc = 0.0, turb_ke = 0.0, DES_LengthScale = 0.0; unsigned short turb_model = config->GetKind_Turb_Model(); - bool physical = true; bool tkeNeeded = ((turb_model == SST) || (turb_model == SST_SUST)); + SU2_OMP_FOR_STAT(omp_chunk_size) for (iPoint = 0; iPoint < nPoint; iPoint++) { /*--- Retrieve the value of the kinetic energy (if needed) ---*/ @@ -183,7 +123,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, C /*--- Incompressible flow, primitive variables --- */ - physical = static_cast(nodes)->SetPrimVar(iPoint,eddy_visc, turb_ke, FluidModel); + bool physical = static_cast(nodes)->SetPrimVar(iPoint,eddy_visc, turb_ke, FluidModel); /* Check for non-realizable states for reporting. */ @@ -193,249 +133,12 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, C nodes->SetDES_LengthScale(iPoint,DES_LengthScale); - /*--- Initialize the convective, source and viscous residual vector ---*/ - - if (!Output) LinSysRes.SetBlock_Zero(iPoint); - } return nonPhysicalPoints; } -void CIncNSSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned long Iteration) { - - su2double Mean_BetaInc2, Area, Vol, Mean_SoundSpeed = 0.0, Mean_ProjVel = 0.0, Lambda, Local_Delta_Time, Local_Delta_Time_Visc, - Global_Delta_Time = 1E6, Mean_LaminarVisc = 0.0, Mean_EddyVisc = 0.0, Mean_Density = 0.0, Mean_Thermal_Conductivity = 0.0, Mean_Cv = 0.0, Lambda_1, Lambda_2, K_v = 0.25, Global_Delta_UnstTimeND; - unsigned long iEdge, iVertex, iPoint = 0, jPoint = 0; - unsigned short iDim, iMarker; - su2double ProjVel, ProjVel_i, ProjVel_j; - const su2double* Normal; - - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool dual_time = ((config->GetTime_Marching() == DT_STEPPING_1ST) || - (config->GetTime_Marching() == DT_STEPPING_2ND)); - bool energy = config->GetEnergy_Equation(); - - Min_Delta_Time = 1.E30; Max_Delta_Time = 0.0; - - /*--- Set maximum inviscid eigenvalue to zero, and compute sound speed and viscosity ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - nodes->SetMax_Lambda_Inv(iPoint,0.0); - nodes->SetMax_Lambda_Visc(iPoint,0.0); - } - - /*--- Loop interior edges ---*/ - - for (iEdge = 0; iEdge < geometry->GetnEdge(); iEdge++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->edges->GetNode(iEdge,0); - jPoint = geometry->edges->GetNode(iEdge,1); - - Normal = geometry->edges->GetNormal(iEdge); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); - Mean_BetaInc2 = 0.5 * (nodes->GetBetaInc2(iPoint) + nodes->GetBetaInc2(jPoint)); - Mean_Density = 0.5 * (nodes->GetDensity(iPoint) + nodes->GetDensity(jPoint)); - Mean_SoundSpeed = sqrt(Mean_BetaInc2*Area*Area); - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); - su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); - ProjVel_i = 0.0; ProjVel_j =0.0; - for (iDim = 0; iDim < nDim; iDim++) { - ProjVel_i += GridVel_i[iDim]*Normal[iDim]; - ProjVel_j += GridVel_j[iDim]*Normal[iDim]; - } - Mean_ProjVel -= 0.5 * (ProjVel_i + ProjVel_j); - } - - /*--- Inviscid contribution ---*/ - - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) nodes->AddMax_Lambda_Inv(iPoint,Lambda); - if (geometry->nodes->GetDomain(jPoint)) nodes->AddMax_Lambda_Inv(jPoint,Lambda); - - /*--- Viscous contribution ---*/ - - Mean_LaminarVisc = 0.5*(nodes->GetLaminarViscosity(iPoint) + nodes->GetLaminarViscosity(jPoint)); - Mean_EddyVisc = 0.5*(nodes->GetEddyViscosity(iPoint) + nodes->GetEddyViscosity(jPoint)); - Mean_Density = 0.5*(nodes->GetDensity(iPoint) + nodes->GetDensity(jPoint)); - Mean_Thermal_Conductivity = 0.5*(nodes->GetThermalConductivity(iPoint) + nodes->GetThermalConductivity(jPoint)); - Mean_Cv = 0.5*(nodes->GetSpecificHeatCv(iPoint) + nodes->GetSpecificHeatCv(jPoint)); - - Lambda_1 = (4.0/3.0)*(Mean_LaminarVisc + Mean_EddyVisc); - Lambda_2 = 0.0; - if (energy) Lambda_2 = (1.0/Mean_Cv)*Mean_Thermal_Conductivity; - Lambda = (Lambda_1 + Lambda_2)*Area*Area/Mean_Density; - - if (geometry->nodes->GetDomain(iPoint)) nodes->AddMax_Lambda_Visc(iPoint,Lambda); - if (geometry->nodes->GetDomain(jPoint)) nodes->AddMax_Lambda_Visc(jPoint,Lambda); - - } - - /*--- Loop boundary edges ---*/ - - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - - Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - Mean_BetaInc2 = nodes->GetBetaInc2(iPoint); - Mean_Density = nodes->GetDensity(iPoint); - Mean_SoundSpeed = sqrt(Mean_BetaInc2*Area*Area); - - /*--- Adjustment for grid movement ---*/ - - if (dynamic_grid) { - su2double *GridVel = geometry->nodes->GetGridVel(iPoint); - ProjVel = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - ProjVel += GridVel[iDim]*Normal[iDim]; - Mean_ProjVel -= ProjVel; - } - - /*--- Inviscid contribution ---*/ - - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) { - nodes->AddMax_Lambda_Inv(iPoint,Lambda); - } - - /*--- Viscous contribution ---*/ - - Mean_LaminarVisc = nodes->GetLaminarViscosity(iPoint); - Mean_EddyVisc = nodes->GetEddyViscosity(iPoint); - Mean_Density = nodes->GetDensity(iPoint); - Mean_Thermal_Conductivity = nodes->GetThermalConductivity(iPoint); - Mean_Cv = nodes->GetSpecificHeatCv(iPoint); - - Lambda_1 = (4.0/3.0)*(Mean_LaminarVisc + Mean_EddyVisc); - Lambda_2 = 0.0; - if (energy) Lambda_2 = (1.0/Mean_Cv)*Mean_Thermal_Conductivity; - Lambda = (Lambda_1 + Lambda_2)*Area*Area/Mean_Density; - - if (geometry->nodes->GetDomain(iPoint)) nodes->AddMax_Lambda_Visc(iPoint,Lambda); - - } - } - } - - /*--- Each element uses their own speed, steady state simulation ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - Vol = geometry->nodes->GetVolume(iPoint); - - if (Vol != 0.0) { - Local_Delta_Time = nodes->GetLocalCFL(iPoint)*Vol / nodes->GetMax_Lambda_Inv(iPoint); - Local_Delta_Time_Visc = nodes->GetLocalCFL(iPoint)*K_v*Vol*Vol/ nodes->GetMax_Lambda_Visc(iPoint); - Local_Delta_Time = min(Local_Delta_Time, Local_Delta_Time_Visc); - Global_Delta_Time = min(Global_Delta_Time, Local_Delta_Time); - Min_Delta_Time = min(Min_Delta_Time, Local_Delta_Time); - Max_Delta_Time = max(Max_Delta_Time, Local_Delta_Time); - if (Local_Delta_Time > config->GetMax_DeltaTime()) - Local_Delta_Time = config->GetMax_DeltaTime(); - nodes->SetDelta_Time(iPoint,Local_Delta_Time); - } - else { - nodes->SetDelta_Time(iPoint,0.0); - } - - } - - /*--- Compute the max and the min dt (in parallel) ---*/ - if (config->GetComm_Level() == COMM_FULL) { -#ifdef HAVE_MPI - su2double rbuf_time, sbuf_time; - sbuf_time = Min_Delta_Time; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Min_Delta_Time = rbuf_time; - - sbuf_time = Max_Delta_Time; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MAX, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Max_Delta_Time = rbuf_time; -#endif - } - - /*--- For exact time solution use the minimum delta time of the whole mesh ---*/ - if (config->GetTime_Marching() == TIME_STEPPING) { -#ifdef HAVE_MPI - su2double rbuf_time, sbuf_time; - sbuf_time = Global_Delta_Time; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Global_Delta_Time = rbuf_time; -#endif - /*--- If the unsteady CFL is set to zero, it uses the defined - unsteady time step, otherwise it computes the time step based - on the unsteady CFL ---*/ - - if (config->GetUnst_CFL() == 0.0) { - Global_Delta_Time = config->GetDelta_UnstTime(); - } - config->SetDelta_UnstTimeND(Global_Delta_Time); - for (iPoint = 0; iPoint < nPointDomain; iPoint++){ - - /*--- Sets the regular CFL equal to the unsteady CFL ---*/ - - nodes->SetLocalCFL(iPoint, config->GetUnst_CFL()); - nodes->SetDelta_Time(iPoint, Global_Delta_Time); - Min_Delta_Time = Global_Delta_Time; - Max_Delta_Time = Global_Delta_Time; - - } - } - - /*--- Recompute the unsteady time step for the dual time strategy - if the unsteady CFL is diferent from 0 ---*/ - if ((dual_time) && (Iteration == 0) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) { - - Global_Delta_UnstTimeND = 1e30; - for (iPoint = 0; iPoint < nPointDomain; iPoint++){ - Global_Delta_UnstTimeND = min(Global_Delta_UnstTimeND,config->GetUnst_CFL()*Global_Delta_Time/nodes->GetLocalCFL(iPoint)); - } - -#ifdef HAVE_MPI - su2double rbuf_time, sbuf_time; - sbuf_time = Global_Delta_UnstTimeND; - SU2_MPI::Reduce(&sbuf_time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MASTER_NODE, MPI_COMM_WORLD); - SU2_MPI::Bcast(&rbuf_time, 1, MPI_DOUBLE, MASTER_NODE, MPI_COMM_WORLD); - Global_Delta_UnstTimeND = rbuf_time; -#endif - config->SetDelta_UnstTimeND(Global_Delta_UnstTimeND); - } - - /*--- The pseudo local time (explicit integration) cannot be greater than the physical time ---*/ - if (dual_time) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - if (!implicit) { - Local_Delta_Time = min((2.0/3.0)*config->GetDelta_UnstTimeND(), nodes->GetDelta_Time(iPoint)); - nodes->SetDelta_Time(iPoint,Local_Delta_Time); - } - } - -} - void CIncNSSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep) { @@ -488,373 +191,222 @@ void CIncNSSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contai } -void CIncNSSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { - - unsigned short iDim, iVar, jVar;// Wall_Function; - unsigned long iVertex, iPoint, total_index; - - su2double *GridVel, *Normal, Area, Wall_HeatFlux; +void CIncNSSolver::BC_Wall_Generic(const CGeometry *geometry, const CConfig *config, + unsigned short val_marker, unsigned short kind_boundary) { - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool energy = config->GetEnergy_Equation(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool energy = config->GetEnergy_Equation(); /*--- Identify the boundary by string name ---*/ - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - - /*--- Get the specified wall heat flux from config ---*/ - - Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag)/config->GetHeat_Flux_Ref(); - -// /*--- Get wall function treatment from config. ---*/ -// -// Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); -// if (Wall_Function != NO_WALL_FUNCTION) { -// SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); -// } - - /*--- Loop over all of the vertices on this boundary marker ---*/ - - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - - /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Compute dual-grid area and boundary normal ---*/ - - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Initialize the convective & viscous residuals to zero ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Res_Conv[iVar] = 0.0; - Res_Visc[iVar] = 0.0; - if (implicit) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - } - - /*--- Store the corrected velocity at the wall which will - be zero (v = 0), unless there are moving walls (v = u_wall)---*/ - - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = GridVel[iDim]; - } else { - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = 0.0; - } - - /*--- Impose the value of the velocity as a strong boundary - condition (Dirichlet). Fix the velocity and remove any - contribution to the residual at this node. ---*/ - - nodes->SetVelocity_Old(iPoint,Vector); - - for (iDim = 0; iDim < nDim; iDim++) - LinSysRes(iPoint, iDim+1) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); - - if (energy) { - - /*--- Apply a weak boundary condition for the energy equation. - Compute the residual due to the prescribed heat flux. ---*/ - - Res_Visc[nDim+1] = Wall_HeatFlux*Area; - - /*--- Viscous contribution to the residual at the wall ---*/ - - LinSysRes.SubtractBlock(iPoint, Res_Visc); - - } - - /*--- Enforce the no-slip boundary condition in a strong way by - modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ - - if (implicit) { - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - } - - } - } -} - -void CIncNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CNumerics *visc_numerics, CConfig *config, unsigned short val_marker) { + const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); - unsigned short iDim, iVar, jVar, Wall_Function; - unsigned long iVertex, iPoint, Point_Normal, total_index; + /*--- Get the specified wall heat flux or temperature from config ---*/ - su2double *GridVel; - su2double *Normal, *Coord_i, *Coord_j, Area, dist_ij; - su2double Twall, dTdn; - su2double thermal_conductivity; + su2double Wall_HeatFlux = 0.0, Twall = 0.0; - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool energy = config->GetEnergy_Equation(); - - /*--- Identify the boundary by string name ---*/ - - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); - - /*--- Retrieve the specified wall temperature ---*/ - - Twall = config->GetIsothermal_Temperature(Marker_Tag)/config->GetTemperature_Ref(); + if (kind_boundary == HEAT_FLUX) + Wall_HeatFlux = config->GetWall_HeatFlux(Marker_Tag)/config->GetHeat_Flux_Ref(); + else if (kind_boundary == ISOTHERMAL) + Twall = config->GetIsothermal_Temperature(Marker_Tag)/config->GetTemperature_Ref(); + else + SU2_MPI::Error("Unknown type of boundary condition", CURRENT_FUNCTION); /*--- Get wall function treatment from config. ---*/ - Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); - if (Wall_Function != NO_WALL_FUNCTION) { - SU2_MPI::Error("Wall function treatment not implemented yet.", CURRENT_FUNCTION); - } + const auto Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); + if (Wall_Function != NO_WALL_FUNCTION) + SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); /*--- Loop over all of the vertices on this boundary marker ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { - - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Check if the node belongs to the domain (i.e, not a halo node) ---*/ - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Initialize the convective & viscous residuals to zero ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Res_Conv[iVar] = 0.0; - Res_Visc[iVar] = 0.0; - if (implicit) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - } - - /*--- Store the corrected velocity at the wall which will - be zero (v = 0), unless there are moving walls (v = u_wall)---*/ + if (!geometry->nodes->GetDomain(iPoint)) continue; - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = GridVel[iDim]; - } else { - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = 0.0; - } + /*--- Compute dual-grid area and boundary normal ---*/ - /*--- Impose the value of the velocity as a strong boundary - condition (Dirichlet). Fix the velocity and remove any - contribution to the residual at this node. ---*/ + const auto Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - nodes->SetVelocity_Old(iPoint,Vector); + const su2double Area = GeometryToolbox::Norm(nDim, Normal); - for (iDim = 0; iDim < nDim; iDim++) - LinSysRes(iPoint, iDim+1) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); + /*--- Impose the value of the velocity as a strong boundary + condition (Dirichlet). Fix the velocity and remove any + contribution to the residual at this node. ---*/ - if (energy) { - - /*--- Compute dual grid area and boundary normal ---*/ - - Normal = geometry->vertex[val_marker][iVertex]->GetNormal(); - - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Compute closest normal neighbor ---*/ + if (dynamic_grid) { + nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint)); + } else { + su2double zero[MAXNDIM] = {0.0}; + nodes->SetVelocity_Old(iPoint, zero); + } - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + LinSysRes(iPoint, iDim+1) = 0.0; + nodes->SetVel_ResTruncError_Zero(iPoint); - /*--- Get coordinates of i & nearest normal and compute distance ---*/ + /*--- Enforce the no-slip boundary condition in a strong way by + modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ - Coord_i = geometry->nodes->GetCoord(iPoint); - Coord_j = geometry->nodes->GetCoord(Point_Normal); - dist_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - dist_ij = sqrt(dist_ij); + if (implicit) { + for (unsigned short iVar = 1; iVar <= nDim; iVar++) + Jacobian.DeleteValsRowi(iPoint*nVar+iVar); + } - /*--- Compute the normal gradient in temperature using Twall ---*/ + if (!energy) continue; - dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij; + if (kind_boundary == HEAT_FLUX) { - /*--- Get thermal conductivity ---*/ + /*--- Apply a weak boundary condition for the energy equation. + Compute the residual due to the prescribed heat flux. ---*/ - thermal_conductivity = nodes->GetThermalConductivity(iPoint); + LinSysRes(iPoint, nDim+1) -= Wall_HeatFlux*Area; + } + else { // ISOTHERMAL - /*--- Apply a weak boundary condition for the energy equation. - Compute the residual due to the prescribed heat flux. ---*/ + auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - Res_Visc[nDim+1] = thermal_conductivity*dTdn*Area; + /*--- Get coordinates of i & nearest normal and compute distance ---*/ - /*--- Jacobian contribution for temperature equation. ---*/ + auto Coord_i = geometry->nodes->GetCoord(iPoint); + auto Coord_j = geometry->nodes->GetCoord(Point_Normal); + su2double Edge_Vector[MAXNDIM]; + GeometryToolbox::Distance(nDim, Coord_j, Coord_i, Edge_Vector); + su2double dist_ij_2 = GeometryToolbox::SquaredNorm(nDim, Edge_Vector); + su2double dist_ij = sqrt(dist_ij_2); - if (implicit) { - su2double Edge_Vector[3]; - su2double dist_ij_2 = 0, proj_vector_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) { - Edge_Vector[iDim] = Coord_j[iDim]-Coord_i[iDim]; - dist_ij_2 += Edge_Vector[iDim]*Edge_Vector[iDim]; - proj_vector_ij += Edge_Vector[iDim]*Normal[iDim]; - } - if (dist_ij_2 == 0.0) proj_vector_ij = 0.0; - else proj_vector_ij = proj_vector_ij/dist_ij_2; + /*--- Compute the normal gradient in temperature using Twall ---*/ - Jacobian_i[nDim+1][nDim+1] = -thermal_conductivity*proj_vector_ij; + su2double dTdn = -(nodes->GetTemperature(Point_Normal) - Twall)/dist_ij; - Jacobian.SubtractBlock2Diag(iPoint, Jacobian_i); - } + /*--- Get thermal conductivity ---*/ - /*--- Viscous contribution to the residual at the wall ---*/ + su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint); - LinSysRes.SubtractBlock(iPoint, Res_Visc); + /*--- Apply a weak boundary condition for the energy equation. + Compute the residual due to the prescribed heat flux. ---*/ - } + LinSysRes(iPoint, nDim+1) -= thermal_conductivity*dTdn*Area; - /*--- Enforce the no-slip boundary condition in a strong way by - modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ + /*--- Jacobian contribution for temperature equation. ---*/ if (implicit) { - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } + su2double proj_vector_ij = 0.0; + if (dist_ij_2 > 0.0) + proj_vector_ij = GeometryToolbox::DotProduct(nDim, Edge_Vector, Normal) / dist_ij_2; + Jacobian.AddVal2Diag(iPoint, nDim+1, thermal_conductivity*proj_vector_ij); } - } } } +void CIncNSSolver::BC_HeatFlux_Wall(CGeometry *geometry, CSolver**, CNumerics*, + CNumerics*, CConfig *config, unsigned short val_marker) { -void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, - CConfig *config, unsigned short val_marker) { + BC_Wall_Generic(geometry, config, val_marker, HEAT_FLUX); +} - unsigned short iVar, jVar, iDim, Wall_Function; - unsigned long iVertex, iPoint, total_index, Point_Normal; +void CIncNSSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver**, CNumerics*, + CNumerics*, CConfig *config, unsigned short val_marker) { - su2double *Coord_i, *Coord_j, dist_ij; - su2double *GridVel, There, Tconjugate, Twall= 0.0, Temperature_Ref, thermal_conductivity, HF_FactorHere, HF_FactorConjugate; + BC_Wall_Generic(geometry, config, val_marker, ISOTHERMAL); +} - Temperature_Ref = config->GetTemperature_Ref(); +void CIncNSSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solver_container, CNumerics *conv_numerics, + CConfig *config, unsigned short val_marker) { - bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool energy = config->GetEnergy_Equation(); + const su2double Temperature_Ref = config->GetTemperature_Ref(); + const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const bool energy = config->GetEnergy_Equation(); /*--- Identify the boundary ---*/ - string Marker_Tag = config->GetMarker_All_TagBound(val_marker); + const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker); /*--- Retrieve the specified wall function treatment.---*/ - Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); - if(Wall_Function != NO_WALL_FUNCTION) { - SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); + const auto Wall_Function = config->GetWallFunction_Treatment(Marker_Tag); + if (Wall_Function != NO_WALL_FUNCTION) { + SU2_MPI::Error("Wall function treament not implemented yet", CURRENT_FUNCTION); } /*--- Loop over boundary points ---*/ - for (iVertex = 0; iVertex < geometry->nVertex[val_marker]; iVertex++) { + for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { - iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); - if (geometry->nodes->GetDomain(iPoint)) { - - /*--- Initialize the convective & viscous residuals to zero ---*/ - - for (iVar = 0; iVar < nVar; iVar++) { - Res_Conv[iVar] = 0.0; - Res_Visc[iVar] = 0.0; - if (implicit) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - } - } - - /*--- Store the corrected velocity at the wall which will - be zero (v = 0), unless there are moving walls (v = u_wall)---*/ - - if (dynamic_grid) { - GridVel = geometry->nodes->GetGridVel(iPoint); - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = GridVel[iDim]; - } else { - for (iDim = 0; iDim < nDim; iDim++) Vector[iDim] = 0.0; - } + if (!geometry->nodes->GetDomain(iPoint)) continue; - /*--- Impose the value of the velocity as a strong boundary - condition (Dirichlet). Fix the velocity and remove any - contribution to the residual at this node. ---*/ + /*--- Impose the value of the velocity as a strong boundary + condition (Dirichlet). Fix the velocity and remove any + contribution to the residual at this node. ---*/ - nodes->SetVelocity_Old(iPoint,Vector); + if (dynamic_grid) { + nodes->SetVelocity_Old(iPoint, geometry->nodes->GetGridVel(iPoint)); + } else { + su2double zero[MAXNDIM] = {0.0}; + nodes->SetVelocity_Old(iPoint, zero); + } - for (iDim = 0; iDim < nDim; iDim++) - LinSysRes(iPoint, iDim+1) = 0.0; - nodes->SetVel_ResTruncError_Zero(iPoint); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + LinSysRes(iPoint, iDim+1) = 0.0; + nodes->SetVel_ResTruncError_Zero(iPoint); - if (energy) { + /*--- Enforce the no-slip boundary condition in a strong way by + modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ - Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0)/Temperature_Ref; + if (implicit) { + for (unsigned short iVar = 1; iVar <= nDim; iVar++) + Jacobian.DeleteValsRowi(iPoint*nVar+iVar); + if (energy) Jacobian.DeleteValsRowi(iPoint*nVar+nDim+1); + } - if ((config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX) || - (config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) { + if (!energy) continue; - /*--- Compute closest normal neighbor ---*/ + su2double Tconjugate = GetConjugateHeatVariable(val_marker, iVertex, 0) / Temperature_Ref; + su2double Twall = 0.0; - Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); + if ((config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_NEUMANN_HEATFLUX) || + (config->GetKind_CHT_Coupling() == AVERAGED_TEMPERATURE_ROBIN_HEATFLUX)) { - /*--- Get coordinates of i & nearest normal and compute distance ---*/ + /*--- Compute closest normal neighbor ---*/ - Coord_i = geometry->nodes->GetCoord(iPoint); - Coord_j = geometry->nodes->GetCoord(Point_Normal); - dist_ij = 0; - for (iDim = 0; iDim < nDim; iDim++) - dist_ij += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - dist_ij = sqrt(dist_ij); + auto Point_Normal = geometry->vertex[val_marker][iVertex]->GetNormal_Neighbor(); - /*--- Compute wall temperature from both temperatures ---*/ + /*--- Get coordinates of i & nearest normal and compute distance ---*/ - thermal_conductivity = nodes->GetThermalConductivity(iPoint); - There = nodes->GetTemperature(Point_Normal); - HF_FactorHere = thermal_conductivity*config->GetViscosity_Ref()/dist_ij; - HF_FactorConjugate = GetConjugateHeatVariable(val_marker, iVertex, 2); + auto Coord_i = geometry->nodes->GetCoord(iPoint); + auto Coord_j = geometry->nodes->GetCoord(Point_Normal); + su2double dist_ij = GeometryToolbox::Distance(nDim, Coord_j, Coord_i); - Twall = (There*HF_FactorHere + Tconjugate*HF_FactorConjugate)/(HF_FactorHere + HF_FactorConjugate); - } - else if ((config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_NEUMANN_HEATFLUX) || - (config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX)) { + /*--- Compute wall temperature from both temperatures ---*/ - /*--- (Directly) Set wall temperature to conjugate temperature. ---*/ + su2double thermal_conductivity = nodes->GetThermalConductivity(iPoint); + su2double There = nodes->GetTemperature(Point_Normal); + su2double HF_FactorHere = thermal_conductivity*config->GetViscosity_Ref()/dist_ij; + su2double HF_FactorConjugate = GetConjugateHeatVariable(val_marker, iVertex, 2); - Twall = Tconjugate; - } - else { - Twall = 0.0; - SU2_MPI::Error("Unknown CHT coupling method.", CURRENT_FUNCTION); - } + Twall = (There*HF_FactorHere + Tconjugate*HF_FactorConjugate)/(HF_FactorHere + HF_FactorConjugate); + } + else if ((config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_NEUMANN_HEATFLUX) || + (config->GetKind_CHT_Coupling() == DIRECT_TEMPERATURE_ROBIN_HEATFLUX)) { - /*--- Strong imposition of the temperature on the fluid zone. ---*/ + /*--- (Directly) Set wall temperature to conjugate temperature. ---*/ - LinSysRes(iPoint, nDim+1) = 0.0; - nodes->SetSolution_Old(iPoint, nDim+1, Twall); - nodes->SetEnergy_ResTruncError_Zero(iPoint); - } + Twall = Tconjugate; + } + else { + SU2_MPI::Error("Unknown CHT coupling method.", CURRENT_FUNCTION); + } - /*--- Enforce the no-slip boundary condition in a strong way by - modifying the velocity-rows of the Jacobian (1 on the diagonal). ---*/ + /*--- Strong imposition of the temperature on the fluid zone. ---*/ - if (implicit) { - for (iVar = 1; iVar <= nDim; iVar++) { - total_index = iPoint*nVar+iVar; - Jacobian.DeleteValsRowi(total_index); - } - if(energy) { - total_index = iPoint*nVar+nDim+1; - Jacobian.DeleteValsRowi(total_index); - } - } - } + LinSysRes(iPoint, nDim+1) = 0.0; + nodes->SetSolution_Old(iPoint, nDim+1, Twall); + nodes->SetEnergy_ResTruncError_Zero(iPoint); } } diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 86aa426bef54..8e18fd1f5bd8 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -432,314 +432,71 @@ unsigned long CNEMOEulerSolver::SetPrimitive_Variables(CSolver **solver_containe void CNEMOEulerSolver::SetTime_Step(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned long Iteration) { - const bool viscous = config->GetViscous(); - const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - const bool time_stepping = (config->GetTime_Marching() == TIME_STEPPING); - const bool dual_time = (config->GetTime_Marching() == DT_STEPPING_1ST) || - (config->GetTime_Marching() == DT_STEPPING_2ND); - const su2double K_v = 0.25; - - /*--- Init thread-shared variables to compute min/max values. - * Critical sections are used for this instead of reduction - * clauses for compatibility with OpenMP 2.0 (Windows...). ---*/ - SU2_OMP_MASTER - { - Min_Delta_Time = 1e30; - Max_Delta_Time = 0.0; - Global_Delta_UnstTimeND = 1e30; - } - SU2_OMP_BARRIER - - const su2double *Normal = nullptr; - su2double Area, Vol, Mean_SoundSpeed, Mean_ProjVel, Lambda, Local_Delta_Time, Local_Delta_Time_Visc; - su2double Mean_LaminarVisc, Mean_EddyVisc, Mean_Density, Lambda_1, Lambda_2; - su2double Mean_ThermalCond, Mean_ThermalCond_ve, cv; - unsigned long iEdge, iVertex, iPoint, jPoint; - unsigned short iDim, iMarker; - - /*--- Loop domain points. ---*/ - SU2_OMP_FOR_DYN(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - /*--- Set maximum eigenvalue to zero. ---*/ - nodes->SetMax_Lambda_Inv(iPoint, 0.0); - - if (viscous) - nodes->SetMax_Lambda_Visc(iPoint,0.0); - - /*--- Loop over the neighbors of point i. ---*/ - for (unsigned short iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) - { - jPoint = geometry->nodes->GetPoint(iPoint,iNeigh); - - iEdge = geometry->nodes->GetEdge(iPoint,iNeigh); - Normal = geometry->edges->GetNormal(iEdge); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint, Normal) + nodes->GetProjVel(jPoint,Normal)); - Mean_SoundSpeed = 0.5 * (nodes->GetSoundSpeed(iPoint) + nodes->GetSoundSpeed(jPoint)) * Area; - - /*--- Adjustment for grid movement ---*/ - if (dynamic_grid) { - const su2double *GridVel_i = geometry->nodes->GetGridVel(iPoint); - const su2double *GridVel_j = geometry->nodes->GetGridVel(jPoint); - - for (iDim = 0; iDim < nDim; iDim++) - Mean_ProjVel -= 0.5 * (GridVel_i[iDim] + GridVel_j[iDim]) * Normal[iDim]; - } - - /*--- Inviscid contribution ---*/ - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - nodes->AddMax_Lambda_Inv(iPoint,Lambda); - - /*--- Viscous contribution ---*/ - if (!viscous) continue; - - /*--- Calculate mean viscous quantities ---*/ - Mean_LaminarVisc = 0.5*(nodes->GetLaminarViscosity(iPoint) + - nodes->GetLaminarViscosity(jPoint)); - Mean_EddyVisc = 0.5*(nodes->GetEddyViscosity(iPoint) + - nodes->GetEddyViscosity(jPoint)); - Mean_ThermalCond = 0.5*(nodes->GetThermalConductivity(iPoint) + - nodes->GetThermalConductivity(jPoint)); - Mean_ThermalCond_ve = 0.5*(nodes->GetThermalConductivity_ve(iPoint) + - nodes->GetThermalConductivity_ve(jPoint)); - Mean_Density = 0.5*(nodes->GetDensity(iPoint) + - nodes->GetDensity(jPoint)); - cv = 0.5*(nodes->GetRhoCv_tr(iPoint) + nodes->GetRhoCv_ve(iPoint) + - nodes->GetRhoCv_tr(jPoint) + nodes->GetRhoCv_ve(jPoint) )/ Mean_Density; - - /*--- Determine the viscous spectral radius and apply it to the control volume ---*/ - Lambda_1 = (4.0/3.0)*(Mean_LaminarVisc + Mean_EddyVisc); - Lambda_2 = (Mean_ThermalCond+Mean_ThermalCond_ve)/cv; - - Lambda = (Lambda_1 + Lambda_2)*Area*Area/Mean_Density; - nodes->AddMax_Lambda_Visc(iPoint, Lambda); + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CNEMOEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return 0.5 * (nodes.GetSoundSpeed(iPoint) + nodes.GetSoundSpeed(jPoint)); } - } - - /*--- Loop boundary edges ---*/ - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - - if (!geometry->nodes->GetDomain(iPoint)) continue; - - Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - Mean_SoundSpeed = nodes->GetSoundSpeed(iPoint) * Area; - - /*--- Adjustment for grid movement ---*/ - if (dynamic_grid) { - const su2double *GridVel = geometry->nodes->GetGridVel(iPoint); - - for (iDim = 0; iDim < nDim; iDim++) - Mean_ProjVel -= GridVel[iDim]*Normal[iDim]; - } - - /*--- Inviscid contribution ---*/ - Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - nodes->AddMax_Lambda_Inv(iPoint,Lambda); - - /*--- Viscous contribution ---*/ - if (!viscous) continue; - - /*--- Calculate viscous mean quantities ---*/ - Mean_LaminarVisc = nodes->GetLaminarViscosity(iPoint); - Mean_EddyVisc = nodes->GetEddyViscosity(iPoint); - Mean_ThermalCond = nodes->GetThermalConductivity(iPoint); - Mean_ThermalCond_ve = nodes->GetThermalConductivity_ve(iPoint); - Mean_Density = nodes->GetDensity(iPoint); - cv = (nodes->GetRhoCv_tr(iPoint) + - nodes->GetRhoCv_ve(iPoint)) / Mean_Density; - - Lambda_1 = (4.0/3.0)*(Mean_LaminarVisc+Mean_EddyVisc); - Lambda_2 = (Mean_ThermalCond+Mean_ThermalCond_ve)/cv; - Lambda = (Lambda_1 + Lambda_2)*Area*Area/Mean_Density; - nodes->AddMax_Lambda_Visc(iPoint,Lambda); - - } + FORCEINLINE su2double operator() (const CNEMOEulerVariable& nodes, unsigned long iPoint) const { + return nodes.GetSoundSpeed(iPoint); } - } - - /*--- Each element uses their own speed, steady state simulation. ---*/ - { - /*--- Thread-local variables for min/max reduction. ---*/ - su2double minDt = 1e30, maxDt = 0.0; - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - Vol = geometry->nodes->GetVolume(iPoint); - - if (Vol != 0.0) { - Local_Delta_Time = nodes->GetLocalCFL(iPoint)*Vol / nodes->GetMax_Lambda_Inv(iPoint); - - if(viscous) { - Local_Delta_Time_Visc = nodes->GetLocalCFL(iPoint)*K_v*Vol*Vol/ nodes->GetMax_Lambda_Visc(iPoint); - Local_Delta_Time = min(Local_Delta_Time, Local_Delta_Time_Visc); - } - minDt = min(minDt, Local_Delta_Time); - maxDt = max(maxDt, Local_Delta_Time); + } soundSpeed; - nodes->SetDelta_Time(iPoint, min(Local_Delta_Time, config->GetMax_DeltaTime())); - } - else { - nodes->SetDelta_Time(iPoint,0.0); - } - } - /*--- Min/max over threads. ---*/ - SU2_OMP_CRITICAL - { - Min_Delta_Time = min(Min_Delta_Time, minDt); - Max_Delta_Time = max(Max_Delta_Time, maxDt); - Global_Delta_Time = Min_Delta_Time; + /*--- Define an object to compute the viscous eigenvalue. ---*/ + struct LambdaVisc { + FORCEINLINE su2double lambda(su2double lamVisc, su2double eddyVisc, su2double rho, su2double k, su2double cv) const { + /*--- Determine the viscous spectral radius and apply it to the control volume ---*/ + su2double Lambda_1 = (4.0/3.0)*(lamVisc + eddyVisc); + return (Lambda_1 + k/cv)/rho; } - SU2_OMP_BARRIER - } - - /*--- Compute the min/max dt (in parallel, now over mpi ranks). ---*/ - SU2_OMP_MASTER - if (config->GetComm_Level() == COMM_FULL) { - su2double rbuf_time; - SU2_MPI::Allreduce(&Min_Delta_Time, &rbuf_time, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - Min_Delta_Time = rbuf_time; - - SU2_MPI::Allreduce(&Max_Delta_Time, &rbuf_time, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - Max_Delta_Time = rbuf_time; - } - SU2_OMP_BARRIER - - /*--- For exact time solution use the minimum delta time of the whole mesh. ---*/ - if (time_stepping) { - /*--- If the unsteady CFL is set to zero, it uses the defined unsteady time step, - * otherwise it computes the time step based on the unsteady CFL. ---*/ - SU2_OMP_MASTER - { - if (config->GetUnst_CFL() == 0.0) { - Global_Delta_Time = config->GetDelta_UnstTime(); - } - else { - Global_Delta_Time = Min_Delta_Time; - } - Max_Delta_Time = Global_Delta_Time; - - config->SetDelta_UnstTimeND(Global_Delta_Time); + FORCEINLINE su2double operator() (const CNEMOEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + su2double lamVisc = 0.5*(nodes.GetLaminarViscosity(iPoint) + nodes.GetLaminarViscosity(jPoint)); + su2double eddyVisc = 0.5*(nodes.GetEddyViscosity(iPoint) + nodes.GetEddyViscosity(jPoint)); + su2double thermalCond = 0.5*(nodes.GetThermalConductivity(iPoint) + nodes.GetThermalConductivity(jPoint) + + nodes.GetThermalConductivity_ve(iPoint) + nodes.GetThermalConductivity_ve(jPoint)); + su2double density = 0.5*(nodes.GetDensity(iPoint) + nodes.GetDensity(jPoint)); + su2double cv = 0.5*(nodes.GetRhoCv_tr(iPoint) + nodes.GetRhoCv_ve(iPoint) + + nodes.GetRhoCv_tr(jPoint) + nodes.GetRhoCv_ve(jPoint))/ density; + return lambda(lamVisc, eddyVisc, density, thermalCond, cv); } - SU2_OMP_BARRIER - /*--- Sets the regular CFL equal to the unsteady CFL. ---*/ - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - nodes->SetLocalCFL(iPoint, config->GetUnst_CFL()); - nodes->SetDelta_Time(iPoint, Global_Delta_Time); + FORCEINLINE su2double operator() (const CNEMOEulerVariable& nodes, unsigned long iPoint) const { + su2double lamVisc = nodes.GetLaminarViscosity(iPoint); + su2double eddyVisc = nodes.GetEddyViscosity(iPoint); + su2double thermalCond = nodes.GetThermalConductivity(iPoint) + nodes.GetThermalConductivity_ve(iPoint); + su2double density = nodes.GetDensity(iPoint); + su2double cv = (nodes.GetRhoCv_tr(iPoint) + nodes.GetRhoCv_ve(iPoint))/ density; + return lambda(lamVisc, eddyVisc, density, thermalCond, cv); } - } + } lambdaVisc; - /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is diferent from 0. ---*/ - if ((dual_time) && (Iteration == 0) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) { + /*--- Now instantiate the generic implementation with the two functors above. ---*/ - /*--- Thread-local variable for reduction. ---*/ - su2double glbDtND = 1e30; - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - glbDtND = min(glbDtND, config->GetUnst_CFL()*Global_Delta_Time / nodes->GetLocalCFL(iPoint)); - } - SU2_OMP_CRITICAL - Global_Delta_UnstTimeND = min(Global_Delta_UnstTimeND, glbDtND); - SU2_OMP_BARRIER - - SU2_OMP_MASTER - { - SU2_MPI::Allreduce(&Global_Delta_UnstTimeND, &glbDtND, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD); - Global_Delta_UnstTimeND = glbDtND; - - config->SetDelta_UnstTimeND(Global_Delta_UnstTimeND); - } - SU2_OMP_BARRIER - } - - /*--- The pseudo local time (explicit integration) cannot be greater than the physical time ---*/ - if (dual_time && !implicit) { - SU2_OMP_FOR_STAT(omp_chunk_size) - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - Local_Delta_Time = min((2.0/3.0)*config->GetDelta_UnstTimeND(), nodes->GetDelta_Time(iPoint)); - nodes->SetDelta_Time(iPoint, Local_Delta_Time); - } - } + SetTime_Step_impl(soundSpeed, lambdaVisc, geometry, solver_container, config, iMesh, Iteration); } void CNEMOEulerSolver::SetMax_Eigenvalue(CGeometry *geometry, CConfig *config) { - /*--- Loop domain points. ---*/ - for ( unsigned long iPoint = 0; iPoint < nPointDomain; ++iPoint) { - - /*--- Set inviscid eigenvalues to zero. ---*/ - nodes->SetLambda(iPoint, 0.0); - - /*--- Loop over the neighbors of point i. ---*/ - for (unsigned short iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) - { - auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); - - auto iEdge = geometry->nodes->GetEdge(iPoint, iNeigh); - auto Normal = geometry->edges->GetNormal(iEdge); - su2double Area = GeometryToolbox::Norm(nDim, Normal); - - /*--- Mean Values ---*/ - su2double Mean_ProjVel = 0.5 * (nodes->GetProjVel(iPoint,Normal) + nodes->GetProjVel(jPoint,Normal)); - su2double Mean_SoundSpeed = 0.5 * (nodes->GetSoundSpeed(iPoint) + nodes->GetSoundSpeed(jPoint)) * Area; - - /*--- Inviscid contribution ---*/ - su2double Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - nodes->AddLambda(iPoint,Lambda); + /*--- Define an object to compute the speed of sound. ---*/ + struct SoundSpeed { + FORCEINLINE su2double operator() (const CNEMOEulerVariable& nodes, unsigned long iPoint, unsigned long jPoint) const { + return 0.5 * (nodes.GetSoundSpeed(iPoint) + nodes.GetSoundSpeed(jPoint)); } - } - - /*--- Loop boundary edges ---*/ - for (unsigned short iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { - if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && - (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { - for (unsigned long iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - - /*--- Point identification, Normal vector and area ---*/ - auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); - auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); - su2double Area = GeometryToolbox::Norm(nDim, Normal); + FORCEINLINE su2double operator() (const CNEMOEulerVariable& nodes, unsigned long iPoint) const { + return nodes.GetSoundSpeed(iPoint); + } - /*--- Mean Values ---*/ - su2double Mean_ProjVel = nodes->GetProjVel(iPoint,Normal); - su2double Mean_SoundSpeed = nodes->GetSoundSpeed(iPoint) * Area; + } soundSpeed; - /*--- Inviscid contribution ---*/ - su2double Lambda = fabs(Mean_ProjVel) + Mean_SoundSpeed; - if (geometry->nodes->GetDomain(iPoint)) { - nodes->AddLambda(iPoint,Lambda); - } - } - } - } + /*--- Instantiate generic implementation. ---*/ - /*--- Call the MPI routine ---*/ - InitiateComms(geometry, config, MAX_EIGENVALUE); - CompleteComms(geometry, config, MAX_EIGENVALUE); + SetMax_Eigenvalue_impl(soundSpeed, geometry, config); } @@ -1300,162 +1057,31 @@ void CNEMOEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_con } } -void CNEMOEulerSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - - su2double *local_Residual, *local_Res_TruncError, Vol, Delta, Res; - unsigned short iVar; - unsigned long iPoint; - - bool adjoint = config->GetContinuous_Adjoint(); - - for (iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - - /*--- Update the solution ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - Vol = (geometry->nodes->GetVolume(iPoint) + - geometry->nodes->GetPeriodicVolume(iPoint)); - - Delta = nodes->GetDelta_Time(iPoint) / Vol; - - local_Res_TruncError = nodes->GetResTruncError(iPoint); - local_Residual = LinSysRes.GetBlock(iPoint); - - if (!adjoint) { - for (iVar = 0; iVar < nVar; iVar++) { - - Res = local_Residual[iVar] + local_Res_TruncError[iVar]; - nodes->AddSolution(iPoint, iVar, -Res*Delta); - AddRes_RMS(iVar, Res*Res); - AddRes_Max(iVar, fabs(Res), geometry->nodes->GetGlobalIndex(iPoint), geometry->nodes->GetCoord(iPoint)); - - } - } - } - - /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); +void CNEMOEulerSolver::ExplicitRK_Iteration(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iRKStep) { + Explicit_Iteration(geometry, solver_container, config, iRKStep); } -void CNEMOEulerSolver::ExplicitRK_Iteration(CGeometry *geometry,CSolver **solver_container, CConfig *config, unsigned short iRKStep) { - - su2double *Residual, *Res_TruncError, Vol, Delta, Res; - unsigned short iVar; - unsigned long iPoint; - - su2double RK_AlphaCoeff = config->Get_Alpha_RKStep(iRKStep); - - for (iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } - - /*--- Update the solution ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - Vol = geometry-> nodes->GetVolume(iPoint); - Delta = nodes->GetDelta_Time(iPoint) / Vol; - - Res_TruncError = nodes->GetResTruncError(iPoint); - Residual = LinSysRes.GetBlock(iPoint); - - for (iVar = 0; iVar < nVar; iVar++) { - Res = Residual[iVar] + Res_TruncError[iVar]; - nodes->AddSolution(iPoint,iVar, -Res*Delta*RK_AlphaCoeff); - AddRes_RMS(iVar, Res*Res); - AddRes_Max(iVar, fabs(Res), geometry-> nodes->GetGlobalIndex(iPoint),geometry->nodes->GetCoord(iPoint)); - } - } +void CNEMOEulerSolver::ClassicalRK4_Iteration(CGeometry *geometry, CSolver **solver_container, + CConfig *config, unsigned short iRKStep) { - /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); + Explicit_Iteration(geometry, solver_container, config, iRKStep); +} - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); +void CNEMOEulerSolver::ExplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { + Explicit_Iteration(geometry, solver_container, config, 0); } void CNEMOEulerSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { - unsigned short iVar; - unsigned long iPoint, total_index, IterLinSol = 0; - su2double Delta, *local_Res_TruncError, Vol; - - /*--- Set maximum residual to zero ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - SetRes_RMS(iVar, 0.0); - SetRes_Max(iVar, 0.0, 0); - } + struct DummyPrec { + const bool active = false; + FORCEINLINE su2double** operator() (const CConfig*, unsigned long, su2double) const { return nullptr; } + } precond; - /*--- Build implicit system ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Read the residual ---*/ - local_Res_TruncError = nodes->GetResTruncError(iPoint); - - /*--- Read the volume ---*/ - Vol = geometry-> nodes->GetVolume(iPoint); - - /*--- Modify matrix diagonal to assure diagonal dominance ---*/ - if (nodes->GetDelta_Time(iPoint) != 0.0) { - Delta = Vol / nodes->GetDelta_Time(iPoint); - Jacobian.AddVal2Diag(iPoint, Delta); - } - else { - Jacobian.SetVal2Diag(iPoint, 1.0); - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = 0.0; - local_Res_TruncError[iVar] = 0.0; - } - } - - /*--- Right hand side of the system (-Residual) and initial guess (x = 0) ---*/ - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = - (LinSysRes[total_index] + local_Res_TruncError[iVar]); - LinSysSol[total_index] = 0.0; - AddRes_RMS(iVar, LinSysRes[total_index]*LinSysRes[total_index]); - AddRes_Max(iVar, fabs(LinSysRes[total_index]), geometry-> nodes->GetGlobalIndex(iPoint), geometry->nodes->GetCoord(iPoint)); - } - } - - /*--- Initialize residual and solution at the ghost points ---*/ - for (iPoint = nPointDomain; iPoint < nPoint; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar + iVar; - LinSysRes[total_index] = 0.0; - LinSysSol[total_index] = 0.0; - } - } - - /*--- Solve or smooth the linear system ---*/ - IterLinSol = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config); - - /*--- The the number of iterations of the linear solver ---*/ - SetIterLinSolver(IterLinSol); - - /*--- Update solution (system written in terms of increments) ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - nodes->AddSolution(iPoint,iVar, nodes->GetUnderRelaxation(iPoint)*LinSysSol[iPoint*nVar+iVar]); - } - } - - /*--- MPI solution ---*/ - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - - /*--- Compute the root mean square residual ---*/ - SetResidual_RMS(geometry, config); + ImplicitEuler_Iteration_impl(precond, geometry, solver_container, config, false); } void CNEMOEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMesh) { @@ -2914,69 +2540,6 @@ void CNEMOEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solut // //} -void CNEMOEulerSolver::SetResidual_DualTime(CGeometry *geometry, - CSolver **solution_container, - CConfig *config, - unsigned short iRKStep, - unsigned short iMesh, - unsigned short RunTime_EqSystem) { - unsigned short iVar, jVar; - unsigned long iPoint; - su2double *U_time_nM1, *U_time_n, *U_time_nP1, Volume_nM1, Volume_n, Volume_nP1, TimeStep; - - bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - bool dynamic_grid = config->GetGrid_Movement(); - - /*--- loop over points ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - - /*--- Solution at time n-1, n and n+1 ---*/ - U_time_nM1 = nodes->GetSolution_time_n1(iPoint); - U_time_n = nodes->GetSolution_time_n(iPoint); - U_time_nP1 = nodes->GetSolution(iPoint); - - /*--- Volume at time n-1 and n ---*/ - if (dynamic_grid) { - Volume_nM1 = geometry->nodes->GetVolume_nM1(iPoint); - Volume_n = geometry->nodes->GetVolume_n(iPoint); - Volume_nP1 = geometry->nodes->GetVolume(iPoint); - } - else { - Volume_nM1 = geometry->nodes->GetVolume(iPoint); - Volume_n = geometry->nodes->GetVolume(iPoint); - Volume_nP1 = geometry->nodes->GetVolume(iPoint); - } - - /*--- Time Step ---*/ - TimeStep = config->GetDelta_UnstTimeND(); - - /*--- Compute Residual ---*/ - for(iVar = 0; iVar < nVar; iVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Residual[iVar] = ( U_time_nP1[iVar]*Volume_nP1 - U_time_n[iVar]*Volume_n ) / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Residual[iVar] = ( 3.0*U_time_nP1[iVar]*Volume_nP1 - 4.0*U_time_n[iVar]*Volume_n - + 1.0*U_time_nM1[iVar]*Volume_nM1 ) / (2.0*TimeStep); - } - - /*--- Add Residual ---*/ - LinSysRes.AddBlock(iPoint, Residual); - - if (implicit) { - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) - Jacobian_i[iVar][jVar] = 0.0; - - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep); - } - Jacobian.AddBlock(iPoint, iPoint, Jacobian_i); - } - } -} - void CNEMOEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *config, int val_iter, bool val_update_geo) { /*--- Restart the solution from file information ---*/ diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index 9ea15c36e10a..8b36400de1ff 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -148,46 +148,7 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C SetPrimitive_Limiter(geometry, config); } - /*--- Evaluate the vorticity and strain rate magnitude ---*/ - - SU2_OMP_MASTER - { - StrainMag_Max = 0.0; - Omega_Max = 0.0; - } - SU2_OMP_BARRIER - - nodes->SetVorticity_StrainMag(); - - /*--- Min and Max are not really differentiable ---*/ - const bool wasActive = AD::BeginPassive(); - - su2double strainMax = 0.0, omegaMax = 0.0; - - SU2_OMP(for schedule(static,omp_chunk_size) nowait) - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - strainMax = max(strainMax, nodes->GetStrainMag(iPoint)); - omegaMax = max(omegaMax, GeometryToolbox::Norm(3, nodes->GetVorticity(iPoint))); - } - SU2_OMP_CRITICAL { - StrainMag_Max = max(StrainMag_Max, strainMax); - Omega_Max = max(Omega_Max, omegaMax); - } - - if ((iMesh == MESH_0) && (config->GetComm_Level() == COMM_FULL)) { - SU2_OMP_BARRIER - SU2_OMP_MASTER - { - su2double MyOmega_Max = Omega_Max; - su2double MyStrainMag_Max = StrainMag_Max; - - SU2_MPI::Allreduce(&MyStrainMag_Max, &StrainMag_Max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - SU2_MPI::Allreduce(&MyOmega_Max, &Omega_Max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - } - SU2_OMP_BARRIER - } - - AD::EndPassive(wasActive); + ComputeVorticityAndStrainMag(*config, iMesh); /*--- Compute the TauWall from the wall functions ---*/ @@ -197,7 +158,7 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C } -unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { +unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, const CConfig *config) { /*--- Number of non-physical points, local to the thread, needs * further reduction if function is called in parallel ---*/ diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 286a4cb1105e..5a01f8177fd6 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -987,7 +987,7 @@ def main(): channel_2D.cfg_dir = "sliding_interface/channel_2D" channel_2D.cfg_file = "channel_2D_WA.cfg" channel_2D.test_iter = 2 - channel_2D.test_vals = [2.000000, 0.000000, 0.398053, 0.352779, 0.405462] + channel_2D.test_vals = [2.000000, 0.000000, 0.397970, 0.352779, 0.405462] channel_2D.su2_exec = "parallel_computation.py -f" channel_2D.timeout = 100 channel_2D.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index ec181e1b07b7..36e1985f4a3a 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1144,7 +1144,7 @@ def main(): channel_2D.cfg_dir = "sliding_interface/channel_2D" channel_2D.cfg_file = "channel_2D_WA.cfg" channel_2D.test_iter = 2 - channel_2D.test_vals = [2.000000, 0.000000, 0.397985, 0.352786, 0.405475] #last 4 columns + channel_2D.test_vals = [2.000000, 0.000000, 0.398017, 0.352786, 0.405475] #last 4 columns channel_2D.su2_exec = "SU2_CFD" channel_2D.timeout = 100 channel_2D.tol = 0.00001