From f7f1b967c3e6cec4ca7635eda6862436d139ae60 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 9 May 2021 10:46:37 +0100 Subject: [PATCH 01/21] local indices to allow parallel registration --- Common/include/geometry/CGeometry.hpp | 3 +- Common/include/geometry/dual_grid/CPoint.hpp | 59 ++++------------ Common/src/geometry/CGeometry.cpp | 10 +-- Common/src/geometry/dual_grid/CPoint.cpp | 13 +--- SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp | 19 ++--- SU2_CFD/include/variables/CVariable.hpp | 45 ++---------- .../src/iteration/CDiscAdjFEAIteration.cpp | 2 +- .../src/iteration/CDiscAdjFluidIteration.cpp | 2 +- .../src/iteration/CDiscAdjHeatIteration.cpp | 2 +- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 69 +++++-------------- SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 53 +++----------- SU2_CFD/src/variables/CVariable.cpp | 26 +++---- 12 files changed, 71 insertions(+), 232 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index cc8aa2867b2..ad4f6d80e67 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1250,9 +1250,8 @@ class CGeometry { /*! * \brief Register the coordinates of the mesh nodes. - * \param[in] config */ - void RegisterCoordinates(const CConfig *config) const; + void RegisterCoordinates() const; /*! * \brief Update the multi-grid structure and the wall-distance. diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 1ef38f5cd24..83f1b89b7ad 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -811,50 +811,21 @@ class CPoint { } /*! - * \brief Set the adjoint values of the coordinates. - * \param[in] iPoint - Index of the point. - * \param[in] adj_sol - The adjoint values of the coordinates. - */ - inline void SetAdjointCoord(unsigned long iPoint, const su2double *adj_coor) { - for (unsigned long iDim = 0; iDim < nDim; iDim++) - SU2_TYPE::SetDerivative(Coord(iPoint,iDim), SU2_TYPE::GetValue(adj_coor[iDim])); + * \brief Register coordinates of a point. + * \param[in] iPoint - Index of the point. + * \param[in] input - Register as input or output. + */ + inline void RegisterCoordinates(unsigned long iPoint, bool input) { + for (unsigned long iDim = 0; iDim < nDim; iDim++) { + if(input) { + AD::RegisterInput(Coord(iPoint,iDim),false); + AD::SetIndex(AD_InputIndex(iPoint,iDim), Coord(iPoint,iDim)); + } + else { + AD::RegisterOutput(Coord(iPoint,iDim)); + AD::SetIndex(AD_OutputIndex(iPoint,iDim), Coord(iPoint,iDim)); + } + } } - /*! - * \brief Set the adjoint values of the coordinates. - * \param[in] iPoint - Index of the point. - * \param[in] adj_sol - The adjoint values of the coordinates. - */ - inline void SetAdjointCoord_LocalIndex(unsigned long iPoint, const su2double *adj_coor) { - for (unsigned long iDim = 0; iDim < nDim; iDim++) - AD::SetDerivative(AD_OutputIndex(iPoint,iDim), SU2_TYPE::GetValue(adj_coor[iDim])); - } - - /*! - * \brief Get the adjoint values of the coordinates. - * \param[in] iPoint - Index of the point. - * \param[in] adj_sol - The adjoint values of the coordinates. - */ - inline void GetAdjointCoord(unsigned long iPoint, su2double *adj_coor) const { - for (unsigned long iDim = 0; iDim < nDim; iDim++) - adj_coor[iDim] = SU2_TYPE::GetDerivative(Coord(iPoint,iDim)); - } - - /*! - * \brief Get the adjoint values of the coordinates. - * \param[in] iPoint - Index of the point. - * \param[in] adj_sol - The adjoint values of the coordinates. - */ - inline void GetAdjointCoord_LocalIndex(unsigned long iPoint, su2double *adj_coor) const { - for (unsigned long iDim = 0; iDim < nDim; iDim++) - adj_coor[iDim] = AD::GetDerivative(AD_InputIndex(iPoint,iDim)); - } - - /*! - * \brief Set the adjoint vector indices of Coord vector. - * \param[in] iPoint - Index of the point. - * \param[in] input - Save them to the input or output indices vector. - */ - void SetIndex(unsigned long iPoint, bool input); - }; diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index cf708b8a112..d4aab89d4c0 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -2475,18 +2475,12 @@ void CGeometry::ComputeAirfoil_Section(su2double *Plane_P0, su2double *Plane_Nor } -void CGeometry::RegisterCoordinates(const CConfig *config) const { +void CGeometry::RegisterCoordinates() const { const bool input = true; - const bool push_index = config->GetMultizone_Problem()? false : true; SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - for (auto iDim = 0u; iDim < nDim; iDim++) { - AD::RegisterInput(nodes->GetCoord(iPoint)[iDim], push_index); - } - if(!push_index) { - nodes->SetIndex(iPoint, input); - } + nodes->RegisterCoordinates(iPoint, input); } END_SU2_OMP_FOR } diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index bc568662729..0c5e1714252 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -73,7 +73,7 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig *config) { } } - if(config->GetAD_Mode() && config->GetMultizone_Problem()) { + if (config->GetDiscrete_Adjoint()) { AD_InputIndex.resize(npoint,nDim) = 0; AD_OutputIndex.resize(npoint,nDim) = 0; } @@ -147,17 +147,6 @@ void CPoint::SetPoints(const vector >& pointsMatrix) { Edge = CCompressedSparsePatternL(Point.outerPtr(), Point.outerPtr()+Point.getOuterSize()+1, long(-1)); } -void CPoint::SetIndex(unsigned long iPoint, bool input) { - for (unsigned long iDim = 0; iDim < nDim; iDim++) { - if(input) { - AD::SetIndex(AD_InputIndex(iPoint,iDim), Coord(iPoint,iDim)); - } - else { - AD::SetIndex(AD_OutputIndex(iPoint,iDim), Coord(iPoint,iDim)); - } - } -} - void CPoint::SetVolume_n() { assert(Volume_n.size() == Volume.size()); parallelCopy(Volume.size(), Volume.data(), Volume_n.data()); diff --git a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp index 613b5a7fb35..e81ac660eb2 100644 --- a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp @@ -49,7 +49,6 @@ class CDiscAdjFEASolver final : public CSolver { unsigned short size = 0; su2double* val = nullptr; /*!< \brief Value of the variable. */ int* AD_Idx = nullptr; /*!< \brief Derivative index in the AD tape. */ - bool localIdx = false; su2double* LocalSens = nullptr; /*!< \brief Local sensitivity (domain). */ su2double* GlobalSens = nullptr; /*!< \brief Global sensitivity (mpi). */ su2double* TotalSens = nullptr; /*!< \brief Total sensitivity (time domain). */ @@ -69,7 +68,6 @@ class CDiscAdjFEASolver final : public CSolver { void clear() { size = 0; - localIdx = false; delete [] val; delete [] AD_Idx; delete [] LocalSens; @@ -77,20 +75,15 @@ class CDiscAdjFEASolver final : public CSolver { delete [] TotalSens; } - void Register(bool push_index) { - for (auto i = 0u; i < size; ++i) AD::RegisterInput(val[i], push_index); - } - - void SetIndex() { - for (auto i = 0u; i < size; ++i) AD::SetIndex(AD_Idx[i], val[i]); - localIdx = true; + void Register() { + for (auto i = 0u; i < size; ++i) { + AD::RegisterInput(val[i], false); + AD::SetIndex(AD_Idx[i], val[i]); + } } void GetDerivative() { - if (localIdx) - for (auto i = 0u; i < size; ++i) LocalSens[i] = AD::GetDerivative(AD_Idx[i]); - else - for (auto i = 0u; i < size; ++i) LocalSens[i] = SU2_TYPE::GetDerivative(val[i]); + for (auto i = 0u; i < size; ++i) LocalSens[i] = AD::GetDerivative(AD_Idx[i]); SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); } diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index db37de80848..12826d23988 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2112,80 +2112,43 @@ class CVariable { /*! * \brief Register the variables in the solution array as input/output variable. * \param[in] input - input or output variables. - * \param[in] push_index - boolean whether we want to push the index or save it in a member variable. */ - void RegisterSolution(bool input, bool push_index = true); + void RegisterSolution(bool input); /*! * \brief Register the variables in the solution_time_n array as input/output variable. */ - void RegisterSolution_time_n(bool push_index); + void RegisterSolution_time_n(); /*! * \brief Register the variables in the solution_time_n1 array as input/output variable. */ - void RegisterSolution_time_n1(bool push_index); + void RegisterSolution_time_n1(); /*! * \brief Set the adjoint values of the solution. * \param[in] adj_sol - The adjoint values of the solution. */ inline void SetAdjointSolution(unsigned long iPoint, const su2double *adj_sol) { - for (unsigned long iVar = 0; iVar < Solution.cols(); iVar++) - SU2_TYPE::SetDerivative(Solution(iPoint,iVar), SU2_TYPE::GetValue(adj_sol[iVar])); - } - - /*! - * \brief Set the adjoint values of the solution. - * \param[in] adj_sol - The adjoint values of the solution. - */ - inline void SetAdjointSolution_LocalIndex(unsigned long iPoint, const su2double *adj_sol) { for (unsigned long iVar = 0; iVar < AD_OutputIndex.cols(); iVar++) AD::SetDerivative(AD_OutputIndex(iPoint,iVar), SU2_TYPE::GetValue(adj_sol[iVar])); } - /*! - * \brief Get the adjoint values of the solution. - * \param[out] adj_sol - The adjoint values of the solution. - */ - inline void GetAdjointSolution(unsigned long iPoint, su2double *adj_sol) const { - for (unsigned long iVar = 0; iVar < Solution.cols(); iVar++) - adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution(iPoint,iVar)); - } - /*! * \brief Get the adjoint values of the solution. * \param[in] adj_sol - The adjoint values of the solution. */ - inline void GetAdjointSolution_LocalIndex(unsigned long iPoint, su2double *adj_sol) const { + inline void GetAdjointSolution(unsigned long iPoint, su2double *adj_sol) const { for (unsigned long iVar = 0; iVar < AD_InputIndex.cols(); iVar++) adj_sol[iVar] = AD::GetDerivative(AD_InputIndex(iPoint,iVar)); } - /*! - * \brief Get the adjoint values of the solution at time n. - * \param[out] adj_sol - The adjoint values of the solution. - */ inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const { - for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++) - adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n(iPoint,iVar)); - } - - inline void GetAdjointSolution_time_n_LocalIndex(unsigned long iPoint, su2double *adj_sol) const { for (unsigned long iVar = 0; iVar < AD_Time_n_InputIndex.cols(); iVar++) adj_sol[iVar] = AD::GetDerivative(AD_Time_n_InputIndex(iPoint,iVar)); } - /*! - * \brief Get the adjoint values of the solution at time n-1. - * \param[out] adj_sol - The adjoint values of the solution. - */ inline void GetAdjointSolution_time_n1(unsigned long iPoint, su2double *adj_sol) const { - for (unsigned long iVar = 0; iVar < nVar; iVar++) - adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n1(iPoint,iVar)); - } - - inline void GetAdjointSolution_time_n1_LocalIndex(unsigned long iPoint, su2double *adj_sol) const { for (unsigned long iVar = 0; iVar < nVar; iVar++) adj_sol[iVar] = AD::GetDerivative(AD_Time_n1_InputIndex(iPoint,iVar)); } diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index de4b30d7ee7..fdd5a885e17 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -178,7 +178,7 @@ void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geom } /*--- Register mesh coordinates for geometric sensitivities ---*/ - geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]); + geometry[iZone][iInst][MESH_0]->RegisterCoordinates(); } void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 9016d1e03fe..1c3ac43204e 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -456,7 +456,7 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge if (kind_recording == RECORDING::MESH_COORDS) { /*--- Register node coordinates as input ---*/ - geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]); + geometry[iZone][iInst][MESH_0]->RegisterCoordinates(); } if (config[iZone]->GetDeform_Mesh()) { diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp index 987ca9a48a8..a60cf7bb23e 100644 --- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp @@ -187,7 +187,7 @@ void CDiscAdjHeatIteration::RegisterInput(CSolver***** solver, CGeometry**** geo else if (kind_recording == RECORDING::MESH_COORDS) { /*--- Register node coordinates as input ---*/ - geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]); + geometry[iZone][iInst][MESH_0]->RegisterCoordinates(); } else if (kind_recording == RECORDING::MESH_DEFORM) { /*--- Register the variables of the mesh deformation ---*/ diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 11148489c3b..ab8b9fb3b5c 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -164,17 +164,16 @@ void CDiscAdjFEASolver::RegisterSolution(CGeometry *geometry, CConfig *config){ const bool input = true; const bool dynamic = config->GetTime_Domain(); - const bool push_index = !config->GetMultizone_Problem(); /*--- Register solution at all necessary time instances and other variables on the tape ---*/ - direct_solver->GetNodes()->RegisterSolution(input, push_index); + direct_solver->GetNodes()->RegisterSolution(input); if (dynamic) { /*--- Register solution (u), acceleration (u'') and velocity (u') at time step n-1 ---*/ - direct_solver->GetNodes()->RegisterSolution_time_n(push_index); + direct_solver->GetNodes()->RegisterSolution_time_n(); } } @@ -208,25 +207,12 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config, } if (!reset) { - const bool local_index = config->GetMultizone_Problem(); - const bool push_index = !local_index; - - E.Register(push_index); - Nu.Register(push_index); - Rho.Register(push_index); - Rho_DL.Register(push_index); - if (de_effects) EField.Register(push_index); - if (fea_dv) DV.Register(push_index); - - /*--- Explicitly store the tape indices for when we extract the derivatives ---*/ - if (local_index) { - E.SetIndex(); - Nu.SetIndex(); - Rho.SetIndex(); - Rho_DL.SetIndex(); - if (de_effects) EField.SetIndex(); - if (fea_dv) DV.SetIndex(); - } + E.Register(); + Nu.Register(); + Rho.Register(); + Rho_DL.Register(); + if (de_effects) EField.Register(); + if (fea_dv) DV.Register(); /*--- Register the flow tractions ---*/ if (config->GetnMarker_Fluid_Load() > 0) @@ -243,11 +229,10 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config, void CDiscAdjFEASolver::RegisterOutput(CGeometry *geometry, CConfig *config){ const bool input = false; - const bool push_index = !config->GetMultizone_Problem(); /*--- Register variables as output of the solver iteration ---*/ - direct_solver->GetNodes()->RegisterSolution(input, push_index); + direct_solver->GetNodes()->RegisterSolution(input); } @@ -276,12 +261,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co /*--- Extract the adjoint solution ---*/ - if(multizone) { - direct_solver->GetNodes()->GetAdjointSolution_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); - } + direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); /*--- Store the adjoint solution ---*/ @@ -298,12 +278,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co /*--- Extract the adjoint solution at time n ---*/ - if(multizone) { - direct_solver->GetNodes()->GetAdjointSolution_time_n_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - } + direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); /*--- Store the adjoint solution at time n ---*/ @@ -387,10 +362,7 @@ void CDiscAdjFEASolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config){ } } - if (multizone) - direct_solver->GetNodes()->SetAdjointSolution_LocalIndex(iPoint,Solution); - else - direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution); + direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution); } } @@ -418,24 +390,15 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo /*--- Extract the geometric sensitivities ---*/ - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + const bool time_domain = config->GetTime_Domain(); - su2double *Coord = geometry->nodes->GetCoord(iPoint); + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { for (unsigned short iDim = 0; iDim < nDim; iDim++) { - su2double Sensitivity; - - if(config->GetMultizone_Problem()) { - Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); - } - else { - Sensitivity = SU2_TYPE::GetDerivative(Coord[iDim]); - /*--- Set the index manually to zero. ---*/ - AD::ResetInput(Coord[iDim]); - } + su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); - if (!config->GetTime_Domain()) { + if (!time_domain) { nodes->SetSensitivity(iPoint, iDim, Sensitivity); } else { nodes->SetSensitivity(iPoint, iDim, nodes->GetSensitivity(iPoint, iDim) + Sensitivity); diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index ce0dc2fff6f..9f590a2522f 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -154,18 +154,17 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) { const bool time_n1_needed = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); const bool time_n_needed = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || time_n1_needed; - const bool push_index = !config->GetMultizone_Problem(); /*--- Register solution at all necessary time instances and other variables on the tape ---*/ /*--- Boolean true indicates that an input is registered ---*/ - direct_solver->GetNodes()->RegisterSolution(true, push_index); + direct_solver->GetNodes()->RegisterSolution(true); if (time_n_needed) - direct_solver->GetNodes()->RegisterSolution_time_n(push_index); + direct_solver->GetNodes()->RegisterSolution_time_n(); if (time_n1_needed) - direct_solver->GetNodes()->RegisterSolution_time_n1(push_index); + direct_solver->GetNodes()->RegisterSolution_time_n1(); } void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, bool reset) { @@ -293,11 +292,9 @@ void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, boo void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) { - const bool push_index = !config->GetMultizone_Problem(); - /*--- Register variables as output of the solver iteration. Boolean false indicates that an output is registered ---*/ - direct_solver->GetNodes()->RegisterSolution(false, push_index); + direct_solver->GetNodes()->RegisterSolution(false); } void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){ @@ -326,12 +323,7 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Extract the adjoint solution ---*/ - if(config->GetMultizone_Problem()) { - direct_solver->GetNodes()->GetAdjointSolution_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); - } + direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); /*--- Relax and store the adjoint solution, compute the residuals. ---*/ @@ -374,12 +366,7 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - if(config->GetMultizone_Problem()) { - direct_solver->GetNodes()->GetAdjointSolution_time_n_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - } + direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); } @@ -391,12 +378,7 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - if(config->GetMultizone_Problem()) { - direct_solver->GetNodes()->GetAdjointSolution_time_n1_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); - } + direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); if (!CrossTerm) nodes->Set_Solution_time_n1(iPoint,Solution); } @@ -496,10 +478,8 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) { } /*--- Set the adjoint values of the primal solution. ---*/ - if (multizone) - direct_solver->GetNodes()->SetAdjointSolution_LocalIndex(iPoint,Solution); - else - direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution); + + direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution); } END_SU2_OMP_FOR } @@ -514,22 +494,9 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - auto Coord = geometry->nodes->GetCoord(iPoint); - for (auto iDim = 0u; iDim < nDim; iDim++) { - su2double Sensitivity = 0.0; - - if(config->GetMultizone_Problem()) { - Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); - } - else { - Sensitivity = SU2_TYPE::GetDerivative(Coord[iDim]); - } - - /*--- Set the index manually to zero. ---*/ - - AD::ResetInput(Coord[iDim]); + su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); /*--- If sharp edge, set the sensitivity to 0 on that region ---*/ diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 37622e12441..e66991aa0a1 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -65,11 +65,11 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva if (config->GetTime_Marching() != TIME_MARCHING::STEADY) Solution_time_n1.resize(nPoint,nVar); - if (config->GetMultizone_Problem() && config->GetDiscrete_Adjoint()) { - if (adjoint) { + if (config->GetDiscrete_Adjoint()) { + if (adjoint && config->GetMultizone_Problem()) External.resize(nPoint,nVar) = su2double(0.0); - } - else { + + if (!adjoint) { AD_InputIndex.resize(nPoint,nVar) = -1; AD_OutputIndex.resize(nPoint,nVar) = -1; @@ -118,30 +118,30 @@ void CVariable::Restore_BGSSolution_k() { void CVariable::SetExternalZero() { parallelSet(External.size(), 0.0, External.data()); } namespace { - void RegisterContainer(bool input, bool push_index, su2activematrix& variable, su2matrix& ad_index) { + void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { const auto nPoint = variable.rows(); SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { for(unsigned long iVar=0; iVar Date: Sun, 9 May 2021 12:50:17 +0100 Subject: [PATCH 02/21] reset inputs as they are extracted --- SU2_CFD/include/variables/CMeshVariable.hpp | 8 ++++--- SU2_CFD/include/variables/CVariable.hpp | 17 +++++++++++++- .../src/iteration/CDiscAdjFEAIteration.cpp | 6 ++--- .../src/iteration/CDiscAdjFluidIteration.cpp | 7 +++--- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 3 +++ SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp | 23 +++++++++---------- SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 3 +++ SU2_CFD/src/variables/CMeshVariable.cpp | 12 ++++------ SU2_CFD/src/variables/CVariable.cpp | 17 -------------- TestCases/disc_adj_fsi/Airfoil_2d/config.cfg | 2 +- 10 files changed, 50 insertions(+), 48 deletions(-) diff --git a/SU2_CFD/include/variables/CMeshVariable.hpp b/SU2_CFD/include/variables/CMeshVariable.hpp index 7b917f51ffc..c420e523206 100644 --- a/SU2_CFD/include/variables/CMeshVariable.hpp +++ b/SU2_CFD/include/variables/CMeshVariable.hpp @@ -93,9 +93,11 @@ class CMeshVariable : public CVariable { /*! * \brief Recover the value of the adjoint of the mesh coordinates. */ - inline void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) const final { - for (unsigned long iDim = 0; iDim < nDim; iDim++) - adj_mesh[iDim] = SU2_TYPE::GetDerivative(Mesh_Coord(iPoint,iDim)); + inline void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) final { + for (unsigned long iDim = 0; iDim < nDim; iDim++) { + adj_mesh[iDim] = AD::GetDerivative(AD_InputIndex(iPoint,iDim)); + AD::ResetInput(Mesh_Coord(iPoint,iDim)); + } } }; diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 12826d23988..567e453021f 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -114,6 +114,21 @@ class CVariable { assert(false && "A base method of CVariable was used, but it should have been overridden by the derived class."); } + void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { + const auto nPoint = variable.rows(); + SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + for(unsigned long iVar=0; iVarRegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); - } - /*--- Register mesh coordinates for geometric sensitivities ---*/ + /*--- Register mesh coordinates for geometric sensitivities ---*/ - geometry[iZone][iInst][MESH_0]->RegisterCoordinates(); + geometry[iZone][iInst][MESH_0]->RegisterCoordinates(); + } } void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 1c3ac43204e..0e6f8a5a37d 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -418,7 +418,7 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver***** solver, CGeometry*** solver[iZone][iInst][MESH_0][ADJRAD_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (config[iZone]->GetFluidProblem()) { + if (config[iZone]->GetFluidProblem() && config[iZone]->GetSinglezone_Driver()) { solver[iZone][iInst][MESH_0][FLOW_SOL]->SetVertexTractionsAdjoint(geometry[iZone][iInst][MESH_0], config[iZone]); } @@ -459,8 +459,7 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge geometry[iZone][iInst][MESH_0]->RegisterCoordinates(); } - if (config[iZone]->GetDeform_Mesh()) { - /*--- Register the variables of the mesh deformation ---*/ + if (kind_recording == RECORDING::MESH_DEFORM) { /*--- Undeformed mesh coordinates ---*/ solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); @@ -563,7 +562,7 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver***** solver, CGeometry**** g if (config[iZone]->AddRadiation()) { solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (config[iZone]->GetFluidProblem()) { + if (config[iZone]->GetFluidProblem() && config[iZone]->GetSinglezone_Driver()) { solver[iZone][iInst][MESH_0][FLOW_SOL]->RegisterVertexTractions(geometry[iZone][iInst][MESH_0], config[iZone]); } diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index ab8b9fb3b5c..8115e187d0e 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -394,9 +394,12 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + auto Coord = geometry->nodes->GetCoord(iPoint); + for (unsigned short iDim = 0; iDim < nDim; iDim++) { su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); + AD::ResetInput(Coord[iDim]); if (!time_domain) { nodes->SetSensitivity(iPoint, iDim, Sensitivity); diff --git a/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp index 0c81be0c379..98a93919dc0 100644 --- a/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp @@ -98,18 +98,19 @@ void CDiscAdjMeshSolver::SetRecording(CGeometry* geometry, CConfig *config){ void CDiscAdjMeshSolver::RegisterSolution(CGeometry *geometry, CConfig *config){ - SU2_OMP_MASTER { - /*--- Register reference mesh coordinates ---*/ - direct_solver->GetNodes()->Register_MeshCoord(); - } - END_SU2_OMP_MASTER - SU2_OMP_BARRIER + /*--- Register reference mesh coordinates ---*/ + direct_solver->GetNodes()->Register_MeshCoord(); + } void CDiscAdjMeshSolver::RegisterVariables(CGeometry *geometry, CConfig *config, bool reset){ + /*--- Register boundary displacements as input. + * Except for FSI, where they are determined by the FEA solver. ---*/ + + if (config->GetFSI_Simulation()) return; + SU2_OMP_MASTER { - /*--- Register boundary displacements as input. ---*/ direct_solver->GetNodes()->Register_BoundDisp(); } END_SU2_OMP_MASTER @@ -139,18 +140,16 @@ void CDiscAdjMeshSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *c void CDiscAdjMeshSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config){ - /*--- Extract the sensitivities of the boundary displacements ---*/ + /*--- Extract the sensitivities of the boundary displacements, except for FSI. ---*/ + + if (config->GetFSI_Simulation()) return; SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++){ - /*--- Extract the adjoint solution of the boundary displacements ---*/ - su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjoint_BoundDisp(iPoint,Solution); - /*--- Store the sensitivities of the boundary displacements ---*/ - nodes->SetBoundDisp_Sens(iPoint,Solution); } diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 9f590a2522f..31259d78856 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -494,9 +494,12 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + auto Coord = geometry->nodes->GetCoord(iPoint); + for (auto iDim = 0u; iDim < nDim; iDim++) { su2double Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); + AD::ResetInput(Coord[iDim]); /*--- If sharp edge, set the sensitivity to 0 on that region ---*/ diff --git a/SU2_CFD/src/variables/CMeshVariable.cpp b/SU2_CFD/src/variables/CMeshVariable.cpp index 95aa08ab1aa..aed8430db5e 100644 --- a/SU2_CFD/src/variables/CMeshVariable.cpp +++ b/SU2_CFD/src/variables/CMeshVariable.cpp @@ -31,9 +31,6 @@ CMeshVariable::CMeshVariable(unsigned long npoint, unsigned long ndim, CConfig *config) : CVariable(npoint, ndim, config) { - /*--- Booleans that determine the kind of problems ---*/ - bool time_domain = config->GetTime_Domain(); - /*--- Store the dimensionality of the problem ---*/ nDim = ndim; @@ -41,15 +38,16 @@ CMeshVariable::CMeshVariable(unsigned long npoint, unsigned long ndim, CConfig * Mesh_Coord.resize(nPoint,nDim) = su2double(0.0); WallDistance.resize(nPoint) = su2double(1e-9); + if (config->GetDiscrete_Adjoint()) + AD_InputIndex.resize(nPoint,nDim) = -1; + /*--- Initialize the variables necessary when the problem is time domain ---*/ - if (time_domain) { + if (config->GetTime_Domain()) { Solution_time_n.resize(nPoint,nDim) = su2double(0.0); Solution_time_n1.resize(nPoint,nDim) = su2double(0.0); } } void CMeshVariable::Register_MeshCoord() { - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) - for (unsigned long iDim = 0; iDim < nDim; iDim++) - AD::RegisterInput(Mesh_Coord(iPoint,iDim)); + RegisterContainer(true, Mesh_Coord, AD_InputIndex); } diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index e66991aa0a1..bd39065c7a4 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -117,23 +117,6 @@ void CVariable::Restore_BGSSolution_k() { void CVariable::SetExternalZero() { parallelSet(External.size(), 0.0, External.data()); } -namespace { - void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { - const auto nPoint = variable.rows(); - SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) - for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { - for(unsigned long iVar=0; iVar Date: Mon, 10 May 2021 11:07:34 +0100 Subject: [PATCH 03/21] fix #1285 --- SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index d18a33b2333..ae9f2dd8487 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -943,7 +943,7 @@ void CDiscAdjMultizoneDriver::HandleDataTransfer() { for (unsigned short jZone = 0; jZone < nZone; jZone++){ /*--- The target zone is iZone ---*/ if (jZone != iZone && interface_container[jZone][iZone] != nullptr) { - DeformMesh = DeformMesh || Transfer_Data(jZone, iZone); + DeformMesh |= Transfer_Data(jZone, iZone); } } From 7f9f6479040eb2de9e3b490c1da6c6045c3c742f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 10 May 2021 11:15:09 +0100 Subject: [PATCH 04/21] address #1273 --- Common/include/code_config.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 377432ee945..17315154c34 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -116,9 +116,13 @@ using su2mixedfloat = passivedouble; /*--- Detect if OpDiLib has to be used. ---*/ #if defined(HAVE_OMP) && defined(CODI_REVERSE_TYPE) +#ifndef __INTEL_COMPILER #define HAVE_OPDI +#else +#warning Hybrid parallel reverse mode AD cannot be used with Intel compilers. #endif #if (_OPENMP >= 201811 && !defined(FORCE_OPDI_MACRO_BACKEND)) || defined(FORCE_OPDI_OMPT_BACKEND) #define HAVE_OMPT #endif +#endif From 2144acde713241f87a70c5f4e4d3264517ea97b4 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Thu, 13 May 2021 22:19:45 +0100 Subject: [PATCH 05/21] try to fix restarts --- Common/include/CConfig.hpp | 2 +- Common/src/CConfig.cpp | 4 +- .../include/solvers/CFVMFlowSolverBase.inl | 5 +- SU2_CFD/src/drivers/CDriver.cpp | 185 +++--------------- SU2_CFD/src/drivers/CMultizoneDriver.cpp | 2 +- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 12 +- SU2_CFD/src/solvers/CMeshSolver.cpp | 12 +- 7 files changed, 50 insertions(+), 172 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index c8a9819435f..f403284c37c 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -5235,7 +5235,7 @@ class CConfig { * \param[in] ext - the extension to be added. * \return The new filename */ - string GetFilename(string filename, string ext, unsigned long Iter) const; + string GetFilename(string filename, string ext, int Iter) const; /*! * \brief Append the zone index to the restart or the solution files. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index abe3f235992..8a51da7a728 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -7820,7 +7820,7 @@ CConfig::~CConfig(void) { } -string CConfig::GetFilename(string filename, string ext, unsigned long Iter) const { +string CConfig::GetFilename(string filename, string ext, int Iter) const { /*--- Remove any extension --- */ @@ -7840,7 +7840,7 @@ string CConfig::GetFilename(string filename, string ext, unsigned long Iter) con filename = GetMultiInstance_FileName(filename, GetiInst(), ext); if (GetTime_Domain()){ - filename = GetUnsteady_FileName(filename, (int)Iter, ext); + filename = GetUnsteady_FileName(filename, Iter, ext); } return filename; diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index fa533937244..64d3934c539 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -999,10 +999,13 @@ void CFVMFlowSolverBase::PushSolutionBackInTime(unsigned long TimeIter, bo } if (restart && (TimeIter == config->GetRestart_Iter()) && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { + + const bool update_geo = !config->GetFSI_Simulation(); + /*--- Load an additional restart file for a 2nd-order restart ---*/ solver_container[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver_container, config, config->GetRestart_Iter() - 1, - true); + update_geo); /*--- Load an additional restart file for the turbulence model ---*/ if (rans) diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index e78890f1fd6..d88af7dbc3b 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -250,17 +250,6 @@ CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunica interface_types, interface_container, interpolator_container); } - if(fsi && (config_container[ZONE_0]->GetRestart() || config_container[ZONE_0]->GetDiscrete_Adjoint())){ - if (rank == MASTER_NODE)cout << endl <<"Restarting Fluid and Structural Solvers." << endl; - - for (iZone = 0; iZone < nZone; iZone++) { - for (iInst = 0; iInst < nInst[iZone]; iInst++){ - Solver_Restart(solver_container[iZone][iInst], geometry_container[iZone][iInst], - config_container[iZone], true); - } - } - } - if (config_container[ZONE_0]->GetBoolTurbomachinery()){ if (rank == MASTER_NODE)cout << endl <<"---------------------- Turbomachinery Preprocessing ---------------------" << endl; Turbomachinery_Preprocessing(config_container, geometry_container, solver_container, interface_container); @@ -1137,170 +1126,56 @@ void CDriver::Inlet_Preprocessing(CSolver ***solver, CGeometry **geometry, void CDriver::Solver_Restart(CSolver ***solver, CGeometry **geometry, CConfig *config, bool update_geo) { - bool euler, ns, turbulent, - NEMO_euler, NEMO_ns, - adj_euler, adj_ns, adj_turb, - heat, fem, fem_euler, fem_ns, fem_dg_flow, - template_solver, disc_adj, disc_adj_fem, disc_adj_turb, disc_adj_heat; - int val_iter = 0; - - /*--- Initialize some useful booleans ---*/ - - euler = false; ns = false; turbulent = false; - NEMO_euler = false; NEMO_ns = false; - adj_euler = false; adj_ns = false; adj_turb = false; - fem_euler = false; fem_ns = false; fem_dg_flow = false; - disc_adj = false; - fem = false; disc_adj_fem = false; - disc_adj_turb = false; - heat = false; disc_adj_heat = false; - template_solver = false; - /*--- Check for restarts and use the LoadRestart() routines. ---*/ - bool restart = config->GetRestart(); - bool restart_flow = config->GetRestart_Flow(); - bool no_restart = false; + const bool restart = config->GetRestart(); + const bool restart_flow = config->GetRestart_Flow(); /*--- Adjust iteration number for unsteady restarts. ---*/ - bool dual_time = ((config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || - (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)); - bool time_stepping = config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING; - bool adjoint = (config->GetDiscrete_Adjoint() || config->GetContinuous_Adjoint()); - bool time_domain = (config->GetTime_Domain()); // Dynamic simulation (FSI). + int val_iter = 0; - if (dual_time) { - if (adjoint) val_iter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; - else if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) - val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; - else val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-2; - } + const bool adjoint = (config->GetDiscrete_Adjoint() || config->GetContinuous_Adjoint()); + const bool time_domain = config->GetTime_Domain(); + const bool dt_step_2nd = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) && + !config->GetStructuralProblem() && !config->GetFEMSolver() && + !adjoint && time_domain; - if (time_stepping) { - if (adjoint) val_iter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; - else val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + if (time_domain) { + if (adjoint) val_iter = config->GetUnst_AdjointIter() - 1; + else val_iter = config->GetRestart_Iter() - 1 - dt_step_2nd; } - /*--- Assign booleans ---*/ - - switch (config->GetKind_Solver()) { - case TEMPLATE_SOLVER: template_solver = true; break; - case EULER : case INC_EULER: euler = true; break; - case NEMO_EULER : NEMO_euler = true; break; - case NAVIER_STOKES: case INC_NAVIER_STOKES: ns = true; heat = config->GetWeakly_Coupled_Heat(); break; - case NEMO_NAVIER_STOKES: NEMO_ns = true; break; - case RANS : case INC_RANS: ns = true; turbulent = true; heat = config->GetWeakly_Coupled_Heat(); break; - case FEM_EULER : fem_euler = true; break; - case FEM_NAVIER_STOKES: fem_ns = true; break; - case FEM_RANS : fem_ns = true; break; - case FEM_LES : fem_ns = true; break; - case HEAT_EQUATION: heat = true; break; - case FEM_ELASTICITY: fem = true; break; - case ADJ_EULER : euler = true; adj_euler = true; break; - case ADJ_NAVIER_STOKES : ns = true; turbulent = (config->GetKind_Turb_Model() != NONE); adj_ns = true; break; - case ADJ_RANS : ns = true; turbulent = true; adj_ns = true; adj_turb = (!config->GetFrozen_Visc_Cont()); break; - case DISC_ADJ_EULER: case DISC_ADJ_INC_EULER: euler = true; disc_adj = true; break; - case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_INC_NAVIER_STOKES: ns = true; disc_adj = true; heat = config->GetWeakly_Coupled_Heat(); break; - case DISC_ADJ_RANS: case DISC_ADJ_INC_RANS: ns = true; turbulent = true; disc_adj = true; disc_adj_turb = (!config->GetFrozen_Visc_Disc()); heat = config->GetWeakly_Coupled_Heat(); break; - case DISC_ADJ_FEM_EULER: fem_euler = true; disc_adj = true; break; - case DISC_ADJ_FEM_NS: fem_ns = true; disc_adj = true; break; - case DISC_ADJ_FEM_RANS: fem_ns = true; turbulent = true; disc_adj = true; disc_adj_turb = (!config->GetFrozen_Visc_Disc()); break; - case DISC_ADJ_FEM: fem = true; disc_adj_fem = true; break; - case DISC_ADJ_HEAT: heat = true; disc_adj_heat = true; break; - - } - - /*--- Determine the kind of FEM solver used for the flow. ---*/ - - switch( config->GetKind_FEM_Flow() ) { - case DG: fem_dg_flow = true; break; - } - - /*--- Load restarts for any of the active solver containers. Note that - these restart routines fill the fine grid and interpolate to all MG levels. ---*/ + /*--- Restart direct solvers. ---*/ if (restart || restart_flow) { - if (euler || ns) { - SU2_OMP_PARALLEL_(if(solver[MESH_0][FLOW_SOL]->GetHasHybridParallel())) - solver[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - END_SU2_OMP_PARALLEL - } - if (NEMO_euler || NEMO_ns) { - solver[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (turbulent) { - SU2_OMP_PARALLEL_(if(solver[MESH_0][TURB_SOL]->GetHasHybridParallel())) - solver[MESH_0][TURB_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - END_SU2_OMP_PARALLEL - } - if (config->AddRadiation()) { - solver[MESH_0][RAD_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (fem) { - if (time_domain) { - if (config->GetRestart()) val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; - else val_iter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; + for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { + auto sol = solver[MESH_0][iSol]; + if (sol && !sol->GetAdjoint()) { + /*--- Note that the mesh solver always loads the most recent file (and not -2). ---*/ + SU2_OMP_PARALLEL_(if(sol->GetHasHybridParallel())) + sol->LoadRestart(geometry, solver, config, val_iter + (iSol==MESH_SOL && dt_step_2nd), update_geo); + END_SU2_OMP_PARALLEL } - solver[MESH_0][FEA_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (fem_euler || fem_ns) { - if (fem_dg_flow) - solver[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (heat) { - solver[MESH_0][HEAT_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); } } + /*--- Restart adjoint solvers. ---*/ + if (restart) { - if (template_solver) { - no_restart = true; - } - if (heat) { - solver[MESH_0][HEAT_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (adj_euler || adj_ns) { - solver[MESH_0][ADJFLOW_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (adj_turb) { - no_restart = true; - } - if (disc_adj) { - solver[MESH_0][ADJFLOW_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - if (disc_adj_turb) - solver[MESH_0][ADJTURB_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - if (disc_adj_heat) - solver[MESH_0][ADJHEAT_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - if (config->AddRadiation()) - solver[MESH_0][ADJRAD_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - if (disc_adj_fem) { - if (time_domain) val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; - solver[MESH_0][ADJFEA_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); + if ((config->GetKind_Solver() == TEMPLATE_SOLVER) || + (config->GetKind_Solver() == ADJ_RANS && !config->GetFrozen_Visc_Cont())) { + SU2_MPI::Error("A restart capability has not been implemented yet for this solver.\n" + "Please set RESTART_SOL= NO and try again.", CURRENT_FUNCTION); } - if (disc_adj_heat) { - solver[MESH_0][ADJHEAT_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - } - - if ((restart || restart_flow) && config->GetDeform_Mesh() && update_geo){ - /*--- Always restart with the last state ---*/ - if (config->GetRestart()) val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; - else val_iter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; - solver[MESH_0][MESH_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); - } - /*--- Exit if a restart was requested for a solver that is not available. ---*/ - - if (no_restart) { - SU2_MPI::Error(string("A restart capability has not been implemented yet for this solver.\n") + - string("Please set RESTART_SOL= NO and try again."), CURRENT_FUNCTION); + for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { + auto sol = solver[MESH_0][iSol]; + if (sol && sol->GetAdjoint()) + sol->LoadRestart(geometry, solver, config, val_iter, update_geo); + } } - /*--- Think about calls to pre / post-processing here, plus realizability checks. ---*/ - - } void CDriver::Solver_Postprocessing(CSolver ****solver, CGeometry **geometry, diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 41e5d901f98..db7b9599c33 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -239,7 +239,7 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) { /*--- Set the initial condition for EULER/N-S/RANS ---------------------------------------------*/ /*--- For FSI, the initial conditions are set, after the mesh has been moved. --------------------------------------*/ - if (!fsi && !config_container[iZone]->GetDiscrete_Adjoint() && config_container[iZone]->GetFluidProblem()) { + if (!fsi && config_container[iZone]->GetFluidProblem()) { solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->SetInitialCondition(geometry_container[iZone][INST_0], solver_container[iZone][INST_0], config_container[iZone], TimeIter); diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index f3d3e4b46f3..965665148d9 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -158,14 +158,14 @@ void CSinglezoneDriver::Run() { void CSinglezoneDriver::Postprocess() { - iteration_container[ZONE_0][INST_0]->Postprocess(output_container[ZONE_0], integration_container, geometry_container, solver_container, - numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); + iteration_container[ZONE_0][INST_0]->Postprocess(output_container[ZONE_0], integration_container, geometry_container, solver_container, + numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); - /*--- A corrector step can help preventing numerical instabilities ---*/ + /*--- A corrector step can help preventing numerical instabilities ---*/ - if (config_container[ZONE_0]->GetRelaxation()) - iteration_container[ZONE_0][INST_0]->Relaxation(output_container[ZONE_0], integration_container, geometry_container, solver_container, - numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); + if (config_container[ZONE_0]->GetRelaxation()) + iteration_container[ZONE_0][INST_0]->Relaxation(output_container[ZONE_0], integration_container, geometry_container, solver_container, + numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); } diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 06af9b8408c..62d3b2c7cf8 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -830,8 +830,6 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * minus the coordinates of the reference mesh file ---*/ su2double displ = curr_coord - nodes->GetMesh_Coord(iPoint_Local, iDim); nodes->SetSolution(iPoint_Local, iDim, displ); - su2double vel = Restart_Data[index+iDim+6]; - if (time_domain && config->GetFSI_Simulation()) geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, vel); } /*--- Increment the overall counter for how many points have been loaded. ---*/ @@ -866,8 +864,8 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * domain and on the boundaries. ---*/ UpdateDualGrid(geometry[MESH_0], config); - /*--- For time-domain problems, we need to compute the grid velocities ---*/ - if (time_domain){ + /*--- For time-domain problems, we need to compute the grid velocities. ---*/ + if (time_domain && !config->GetFSI_Simulation()) { /*--- Update the old geometry (coordinates n and n-1) ---*/ Restart_OldGeometry(geometry[MESH_0], config); /*--- Once Displacement_n and Displacement_n1 are filled, @@ -928,8 +926,10 @@ void CMeshSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { /*--- Modify file name for an unsteady restart ---*/ int Unst_RestartIter; - if (config->GetRestart()) Unst_RestartIter = SU2_TYPE::Int(config->GetRestart_Iter()) - iStep; - else Unst_RestartIter = SU2_TYPE::Int(config->GetUnst_AdjointIter()) - SU2_TYPE::Int(config->GetTimeIter())-iStep-1; + if (!config->GetDiscrete_Adjoint()) + Unst_RestartIter = static_cast(config->GetRestart_Iter()) - iStep; + else + Unst_RestartIter = static_cast(config->GetUnst_AdjointIter()) - config->GetTimeIter() - iStep - 1; if (Unst_RestartIter < 0) { From 0d55909d8eaf29855b284aa067ac114c31be5cc7 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 18 May 2021 00:09:36 +0100 Subject: [PATCH 06/21] update regressions --- TestCases/hybrid_regression.py | 6 +++--- TestCases/parallel_regression.py | 6 +++--- TestCases/serial_regression.py | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index fb06b8d5a6a..48b223f8024 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -646,7 +646,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4.000000, 0.000000, -3.768521, -4.159940] + fsi2d.test_vals = [4, 0, -3.771897, -4.163473] fsi2d.multizone= True fsi2d.unsteady = True test_list.append(fsi2d) @@ -656,7 +656,7 @@ def main(): stat_fsi.cfg_dir = "fea_fsi/stat_fsi" stat_fsi.cfg_file = "config.cfg" stat_fsi.test_iter = 7 - stat_fsi.test_vals = [-3.242851, -4.866383, 0.000000, 11.000000] + stat_fsi.test_vals = [-3.242851, -4.866383, 0.000000, 11] stat_fsi.multizone = True test_list.append(stat_fsi) @@ -665,7 +665,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.398527, -4.086735, 0.000000, 121.000000] + dyn_fsi.test_vals = [-4.400852, -4.085149, 5.3839e-08, 121] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index b0c92bcb386..abee592c8eb 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1175,7 +1175,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4, 0, -3.768496, -4.159961] #last 4 columns + fsi2d.test_vals = [4, 0, -3.771873, -4.163491] #last 4 columns fsi2d.su2_exec = "parallel_computation_fsi.py -f" fsi2d.timeout = 1600 fsi2d.multizone= True @@ -1200,7 +1200,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.398550, -4.086739, 0.000000, 117.000000] + dyn_fsi.test_vals = [-4.400874, -4.085152, 0.000000, 117] dyn_fsi.multizone = True dyn_fsi.unsteady = True dyn_fsi.su2_exec = "mpirun -n 2 SU2_CFD" @@ -1356,7 +1356,7 @@ def main(): pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" pywrapper_fsi2d.cfg_file = "configFSI.cfg" pywrapper_fsi2d.test_iter = 4 - pywrapper_fsi2d.test_vals = [4, 0, -3.768496, -4.159961] #last 4 columns + pywrapper_fsi2d.test_vals = [4, 0, -3.771873, -4.163491] #last 4 columns pywrapper_fsi2d.su2_exec = "mpirun -np 2 SU2_CFD.py --nZone 2 --fsi True --parallel -f" pywrapper_fsi2d.timeout = 1600 pywrapper_fsi2d.unsteady = True diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 98fd0feb499..feb426f565e 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1319,7 +1319,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4, 0, -3.768501, -4.159959] #last 4 columns + fsi2d.test_vals = [4, 0, -3.771878, -4.163489] #last 4 columns fsi2d.su2_exec = "SU2_CFD" fsi2d.timeout = 1600 fsi2d.multizone = True @@ -1356,7 +1356,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.398530, -4.086741, 0.000000, 101.000000] #last 4 columns + dyn_fsi.test_vals = [-4.400854, -4.085153, 0.000000, 101] #last 4 columns dyn_fsi.multizone = True dyn_fsi.unsteady = True dyn_fsi.su2_exec = "SU2_CFD" @@ -1874,7 +1874,7 @@ def main(): pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" pywrapper_fsi2d.cfg_file = "configFSI.cfg" pywrapper_fsi2d.test_iter = 4 - pywrapper_fsi2d.test_vals = [4, 0, -3.768501, -4.159959] #last 4 columns + pywrapper_fsi2d.test_vals = [4, 0, -3.771878, -4.163489] #last 4 columns pywrapper_fsi2d.su2_exec = "SU2_CFD.py --nZone 2 --fsi True -f" pywrapper_fsi2d.new_output = True pywrapper_fsi2d.unsteady = True From 8c6233a3a8d180d093eb32f4295e2a184b9cd610 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Tue, 18 May 2021 12:37:07 +0200 Subject: [PATCH 07/21] Move \!Crossterm bool in front of loop for unsteady adjoint extraction. --- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 4 ++-- SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 8115e187d0e..59d2a277e9b 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -271,7 +271,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co /*--- Solution for acceleration (u'') and velocity (u') at time n ---*/ - if (dynamic){ + if (dynamic && !CrossTerm){ /*--- NOW: The solution at time n ---*/ for (iPoint = 0; iPoint < nPoint; iPoint++){ @@ -282,7 +282,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co /*--- Store the adjoint solution at time n ---*/ - if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); + nodes->Set_Solution_time_n(iPoint,Solution); } } diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 31259d78856..680e61670af 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -362,25 +362,25 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi END_SU2_OMP_MASTER /*--- Extract and store the adjoint of the primal solution at time n ---*/ - if (time_n_needed) { + if (time_n_needed && !CrossTerm) { SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); + nodes->Set_Solution_time_n(iPoint,Solution); } END_SU2_OMP_FOR } /*--- Extract and store the adjoint of the primal solution at time n-1 ---*/ - if (time_n1_needed) { + if (time_n1_needed && !CrossTerm) { SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); - if (!CrossTerm) nodes->Set_Solution_time_n1(iPoint,Solution); + nodes->Set_Solution_time_n1(iPoint,Solution); } END_SU2_OMP_FOR } From 7422c0003398393c419d335308686f2cf5224fac Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 18 May 2021 12:07:09 +0100 Subject: [PATCH 08/21] dont extract some terms during cross term evaluation --- SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp | 55 +++++++---------------- SU2_CFD/src/solvers/CDiscAdjSolver.cpp | 26 +++++------ 2 files changed, 29 insertions(+), 52 deletions(-) diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index e044ae9659c..3f0e634230e 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -236,64 +236,41 @@ void CDiscAdjFEASolver::RegisterOutput(CGeometry *geometry, CConfig *config){ } -void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){ - - const bool dynamic = config->GetTime_Domain(); - const bool multizone = config->GetMultizone_Problem(); - - unsigned short iVar; - unsigned long iPoint; - su2double residual; - - su2double Solution[MAXNVAR] = {0.0}; - - /*--- Set Residuals to zero ---*/ - - SetResToZero(); +void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) { /*--- Set the old solution, for multi-zone problems this is done after computing the * residuals, otherwise the per-zone-residuals do not make sense, as on entry Solution * contains contributions from other zones but on extraction it does not. ---*/ - if(!multizone) nodes->Set_OldSolution(); - - for (iPoint = 0; iPoint < nPoint; iPoint++){ + if (!config->GetMultizone_Problem()) nodes->Set_OldSolution(); - /*--- Extract the adjoint solution ---*/ + /*--- Extract and store the adjoint solution ---*/ + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); - - /*--- Store the adjoint solution ---*/ - nodes->SetSolution(iPoint,Solution); - } - /*--- Solution for acceleration (u'') and velocity (u') at time n ---*/ + if (CrossTerm) return; - if (dynamic){ - - /*--- NOW: The solution at time n ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++){ - - /*--- Extract the adjoint solution at time n ---*/ + /*--- Extract and store the adjoint solution at time n (including accel. and velocity) ---*/ + if (config->GetTime_Domain()) { + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - - /*--- Store the adjoint solution at time n ---*/ - - if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); + nodes->Set_Solution_time_n(iPoint,Solution); } - } - /*--- TODO: Need to set the MPI solution in the previous TS ---*/ - /*--- Set the residuals ---*/ - for (iPoint = 0; iPoint < nPointDomain; iPoint++){ - for (iVar = 0; iVar < nVar; iVar++){ - residual = nodes->GetSolution(iPoint, iVar) - nodes->GetSolution_Old(iPoint, iVar); + SetResToZero(); + + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + for (auto iVar = 0u; iVar < nVar; iVar++){ + su2double residual = nodes->GetSolution(iPoint, iVar) - nodes->GetSolution_Old(iPoint, iVar); Residual_RMS[iVar] += residual*residual; AddRes_Max(iVar,fabs(residual),geometry->nodes->GetGlobalIndex(iPoint),geometry->nodes->GetCoord(iPoint)); diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index 31259d78856..84ce1b3906c 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -297,18 +297,14 @@ void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) { direct_solver->GetNodes()->RegisterSolution(false); } -void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){ +void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) { const bool time_n1_needed = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND; const bool time_n_needed = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || time_n1_needed; const su2double relax = (config->GetInnerIter()==0) ? 1.0 : config->GetRelaxation_Factor_Adjoint(); - su2double Solution[MAXNVAR] = {0.0}; - - /*--- Set Residuals to zero ---*/ - - SetResToZero(); + /*--- Thread-local residual variables. ---*/ su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; const su2double* coordMax[MAXNVAR] = {nullptr}; @@ -316,13 +312,14 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Set the old solution and compute residuals. ---*/ - if(!config->GetMultizone_Problem()) nodes->Set_OldSolution(); + if (!config->GetMultizone_Problem()) nodes->Set_OldSolution(); SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { /*--- Extract the adjoint solution ---*/ + su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); /*--- Relax and store the adjoint solution, compute the residuals. ---*/ @@ -344,6 +341,11 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi } END_SU2_OMP_FOR + /*--- Residuals and time_n terms are not needed when evaluating multizone cross terms. ---*/ + if (CrossTerm) return; + + SetResToZero(); + /*--- Reduce residual information over all threads in this rank. ---*/ SU2_OMP_CRITICAL for (auto iVar = 0u; iVar < nVar; iVar++) { @@ -365,10 +367,9 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi if (time_n_needed) { SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - + su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - - if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); + nodes->Set_Solution_time_n(iPoint,Solution); } END_SU2_OMP_FOR } @@ -377,10 +378,9 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi if (time_n1_needed) { SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { - + su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); - - if (!CrossTerm) nodes->Set_Solution_time_n1(iPoint,Solution); + nodes->Set_Solution_time_n1(iPoint,Solution); } END_SU2_OMP_FOR } From d949cbed776ab1088b5ac6d22e708b0f5a891d19 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 18 May 2021 19:31:19 +0100 Subject: [PATCH 09/21] apply BGS relaxation also to velocities --- SU2_CFD/include/variables/CFEAVariable.hpp | 4 ++-- SU2_CFD/include/variables/CVariable.hpp | 2 +- SU2_CFD/src/solvers/CFEASolver.cpp | 17 ++++++++++++----- 3 files changed, 15 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/include/variables/CFEAVariable.hpp b/SU2_CFD/include/variables/CFEAVariable.hpp index 5dc26957717..975f4b24f56 100644 --- a/SU2_CFD/include/variables/CFEAVariable.hpp +++ b/SU2_CFD/include/variables/CFEAVariable.hpp @@ -266,8 +266,8 @@ class CFEAVariable : public CVariable { /*! * \brief Set the value of the solution velocity predictor. */ - inline void SetSolution_Vel_Pred(unsigned long iPoint) final { - for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_Vel_Pred(iPoint,iVar) = Solution_Vel(iPoint,iVar); + inline void SetSolution_Vel_Pred(unsigned long iPoint, const su2double *val_solution_pred) final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_Vel_Pred(iPoint,iVar) = val_solution_pred[iVar]; } /*! diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 567e453021f..8f049516587 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1941,7 +1941,7 @@ class CVariable { /*! * \brief A virtual member. Set the value of the velocity solution predictor. */ - inline virtual void SetSolution_Vel_Pred(unsigned long iPoint) {} + inline virtual void SetSolution_Vel_Pred(unsigned long iPoint, const su2double *val_solution_pred) { } /*! * \brief A virtual member. Set the value of the old solution. diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 30307741607..7763d01f1b1 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2628,7 +2628,7 @@ void CFEASolver::PredictStruct_Displacement(CGeometry *geometry, const CConfig * } break; } - if (dynamic) nodes->SetSolution_Vel_Pred(iPoint); + if (dynamic) nodes->SetSolution_Vel_Pred(iPoint, nodes->GetSolution_Vel(iPoint)); } END_SU2_OMP_PARALLEL @@ -2744,15 +2744,22 @@ void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, const CConfig *config /*--- Set calculated solution as the old solution (needed for dynamic Aitken relaxation) ---*/ nodes->SetSolution_Old(iPoint, dispCalc); - /*--- Set predicted velocity to update in multizone iterations ---*/ - if (dynamic) nodes->SetSolution_Vel_Pred(iPoint); - /*--- Apply the Aitken relaxation ---*/ su2double newDispPred[MAXNVAR] = {0.0}; for (unsigned short iDim=0; iDim < nDim; iDim++) newDispPred[iDim] = (1.0 - WAitken)*dispPred[iDim] + WAitken*dispCalc[iDim]; nodes->SetSolution_Pred(iPoint, newDispPred); + + /*--- Set predicted velocity to update in multizone iterations ---*/ + if (dynamic) { + su2double newVelPred[MAXNVAR] = {0.0}; + const su2double* velPred = nodes->GetSolution_Vel_Pred(iPoint); + const su2double* velCalc = nodes->GetSolution_Vel(iPoint); + for (unsigned short iDim=0; iDim < nDim; iDim++) + newVelPred[iDim] = (1.0 - WAitken)*velPred[iDim] + WAitken*velCalc[iDim]; + nodes->SetSolution_Vel_Pred(iPoint, newVelPred); + } } END_SU2_OMP_PARALLEL @@ -3199,7 +3206,7 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c if (dynamic) { for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) - nodes->SetSolution_Vel_Pred(iPoint); + nodes->SetSolution_Vel_Pred(iPoint, nodes->GetSolution_Vel(iPoint)); } if (discrete_adjoint) { From b7f5b05dc28a49d389abffd1b41c5a2173581d67 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 19 May 2021 14:16:46 +0100 Subject: [PATCH 10/21] cleanup overzealous SU2_OMP_MASTER --- SU2_CFD/src/drivers/CDriver.cpp | 158 +++++++--------------------- SU2_CFD/src/solvers/CMeshSolver.cpp | 32 +++--- 2 files changed, 60 insertions(+), 130 deletions(-) diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index b383c006778..25f56322426 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -140,28 +140,22 @@ CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunica nInst[iZone] = config_container[iZone]->GetnTimeInstances(); - geometry_container[iZone] = new CGeometry** [nInst[iZone]]; - iteration_container[iZone] = new CIteration* [nInst[iZone]]; - solver_container[iZone] = new CSolver*** [nInst[iZone]]; - integration_container[iZone] = new CIntegration** [nInst[iZone]]; - numerics_container[iZone] = new CNumerics**** [nInst[iZone]]; - grid_movement[iZone] = new CVolumetricMovement* [nInst[iZone]]; + geometry_container[iZone] = new CGeometry** [nInst[iZone]] (); + iteration_container[iZone] = new CIteration* [nInst[iZone]] (); + solver_container[iZone] = new CSolver*** [nInst[iZone]] (); + integration_container[iZone] = new CIntegration** [nInst[iZone]] (); + numerics_container[iZone] = new CNumerics**** [nInst[iZone]] (); + grid_movement[iZone] = new CVolumetricMovement* [nInst[iZone]] (); /*--- Allocate transfer and interpolation container --- */ interface_container[iZone] = new CInterface*[nZone] (); interpolator_container[iZone].resize(nZone); - for (iInst = 0; iInst < nInst[iZone]; iInst++){ + for (iInst = 0; iInst < nInst[iZone]; iInst++) { config_container[iZone]->SetiInst(iInst); - geometry_container[iZone][iInst] = nullptr; - iteration_container[iZone][iInst] = nullptr; - solver_container[iZone][iInst] = nullptr; - integration_container[iZone][iInst] = nullptr; - grid_movement[iZone][iInst] = nullptr; - /*--- Preprocessing of the geometry for all zones. In this routine, the edge- based data structure is constructed, i.e. node and cell neighbors are identified and linked, face areas and volumes of the dual mesh cells are @@ -319,40 +313,28 @@ void CDriver::SetContainers_Null(){ interface_types = nullptr; nInst = nullptr; - /*--- Definition and of the containers for all possible zones. ---*/ - iteration_container = new CIteration**[nZone]; - solver_container = new CSolver****[nZone]; - integration_container = new CIntegration***[nZone]; - numerics_container = new CNumerics*****[nZone]; - config_container = new CConfig*[nZone]; - geometry_container = new CGeometry***[nZone]; - surface_movement = new CSurfaceMovement*[nZone]; - grid_movement = new CVolumetricMovement**[nZone]; - FFDBox = new CFreeFormDefBox**[nZone]; + iteration_container = new CIteration**[nZone] (); + solver_container = new CSolver****[nZone] (); + integration_container = new CIntegration***[nZone] (); + numerics_container = new CNumerics*****[nZone] (); + config_container = new CConfig*[nZone] (); + geometry_container = new CGeometry***[nZone] (); + surface_movement = new CSurfaceMovement*[nZone] (); + grid_movement = new CVolumetricMovement**[nZone] (); + FFDBox = new CFreeFormDefBox**[nZone] (); interpolator_container.resize(nZone); - interface_container = new CInterface**[nZone]; - interface_types = new unsigned short*[nZone]; - output_container = new COutput*[nZone]; - nInst = new unsigned short[nZone]; + interface_container = new CInterface**[nZone] (); + interface_types = new unsigned short*[nZone] (); + output_container = new COutput*[nZone] (); + nInst = new unsigned short[nZone] (); driver_config = nullptr; driver_output = nullptr; - for (iZone = 0; iZone < nZone; iZone++) { - solver_container[iZone] = nullptr; - integration_container[iZone] = nullptr; - numerics_container[iZone] = nullptr; - config_container[iZone] = nullptr; - geometry_container[iZone] = nullptr; - surface_movement[iZone] = nullptr; - grid_movement[iZone] = nullptr; - FFDBox[iZone] = nullptr; - interface_container[iZone] = nullptr; - interface_types[iZone] = new unsigned short[nZone]; - output_container[iZone] = nullptr; - nInst[iZone] = 1; + interface_types[iZone] = new unsigned short[nZone]; + nInst[iZone] = 1; } strcpy(runtime_file_name, "runtime.dat"); @@ -425,8 +407,7 @@ void CDriver::Postprocessing() { for (iZone = 0; iZone < nZone; iZone++) { if (interface_container[iZone] != nullptr) { for (unsigned short jZone = 0; jZone < nZone; jZone++) - if (interface_container[iZone][jZone] != nullptr) - delete interface_container[iZone][jZone]; + delete interface_container[iZone][jZone]; delete [] interface_container[iZone]; } } @@ -435,20 +416,17 @@ void CDriver::Postprocessing() { } if (interface_types != nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (interface_types[iZone] != nullptr) + for (iZone = 0; iZone < nZone; iZone++) delete [] interface_types[iZone]; - } delete [] interface_types; } for (iZone = 0; iZone < nZone; iZone++) { if (geometry_container[iZone] != nullptr) { for (iInst = 0; iInst < nInst[iZone]; iInst++){ - for (unsigned short iMGlevel = 0; iMGlevel < config_container[iZone]->GetnMGLevels()+1; iMGlevel++) { - if (geometry_container[iZone][iInst][iMGlevel] != nullptr) delete geometry_container[iZone][iInst][iMGlevel]; - } - if (geometry_container[iZone][iInst] != nullptr) delete [] geometry_container[iZone][iInst]; + for (unsigned short iMGlevel = 0; iMGlevel < config_container[iZone]->GetnMGLevels()+1; iMGlevel++) + delete geometry_container[iZone][iInst][iMGlevel]; + delete [] geometry_container[iZone][iInst]; } delete [] geometry_container[iZone]; } @@ -469,10 +447,9 @@ void CDriver::Postprocessing() { if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl; for (iZone = 0; iZone < nZone; iZone++) { - for (iInst = 0; iInst < nInst[iZone]; iInst++){ - if (grid_movement[iZone][iInst] != nullptr) delete grid_movement[iZone][iInst]; - } - if (grid_movement[iZone] != nullptr) delete [] grid_movement[iZone]; + for (iInst = 0; iInst < nInst[iZone]; iInst++) + delete grid_movement[iZone][iInst]; + delete [] grid_movement[iZone]; } delete [] grid_movement; if (rank == MASTER_NODE) cout << "Deleted CVolumetricMovement class." << endl; @@ -486,11 +463,8 @@ void CDriver::Postprocessing() { /*--- Deallocate config container ---*/ if (config_container!= nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (config_container[iZone] != nullptr) { - delete config_container[iZone]; - } - } + for (iZone = 0; iZone < nZone; iZone++) + delete config_container[iZone]; delete [] config_container; } delete driver_config; @@ -502,18 +476,13 @@ void CDriver::Postprocessing() { /*--- Deallocate output container ---*/ if (output_container!= nullptr) { - for (iZone = 0; iZone < nZone; iZone++) { - if (output_container[iZone] != nullptr) { - delete output_container[iZone]; - } - } + for (iZone = 0; iZone < nZone; iZone++) + delete output_container[iZone]; delete [] output_container; } - delete driver_output; - if (rank == MASTER_NODE) cout << "Deleted COutput class." << endl; if (rank == MASTER_NODE) cout << "-------------------------------------------------------------------------" << endl; @@ -647,7 +616,7 @@ void CDriver::Geometrical_Preprocessing(CConfig* config, CGeometry **&geometry, unsigned short iMGlevel; - geometry = new CGeometry*[config->GetnMGLevels()+1]; + geometry = new CGeometry*[config->GetnMGLevels()+1] (); if (!fem_solver){ for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { @@ -745,7 +714,7 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet /*--- Allocate the memory of the current domain, and divide the grid between the ranks. ---*/ - geometry = new CGeometry *[config->GetnMGLevels()+1]; + geometry = new CGeometry *[config->GetnMGLevels()+1] (); /*--- Build the grid data structures using the ParMETIS coloring. ---*/ @@ -881,6 +850,7 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet if (config->GetnMGLevels() != requestedMGlevels) { delete geometry[iMGlevel]; + geometry[iMGlevel] = nullptr; break; } @@ -953,7 +923,7 @@ void CDriver::Geometrical_Preprocessing_DGFEM(CConfig* config, CGeometry **&geom /*--- Allocate the memory of the current domain, and divide the grid between the ranks. ---*/ - geometry = new CGeometry *[config->GetnMGLevels()+1]; + geometry = new CGeometry *[config->GetnMGLevels()+1] (); geometry[MESH_0] = new CMeshFEM_DG(geometry_aux, config); @@ -1029,7 +999,7 @@ void CDriver::Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolve if (rank == MASTER_NODE) cout << endl <<"-------------------- Solver Preprocessing ( Zone " << config->GetiZone() <<" ) --------------------" << endl; - solver = new CSolver**[config->GetnMGLevels()+1]; + solver = new CSolver**[config->GetnMGLevels()+1] (); for (iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++){ solver[iMesh] = CSolverFactory::CreateSolverContainer(kindSolver, config, geometry[iMesh], iMesh); @@ -1234,7 +1204,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol nVar_Rad = 0, nVar_Heat = 0; - numerics = new CNumerics***[config->GetnMGLevels()+1]; + numerics = new CNumerics***[config->GetnMGLevels()+1] (); const su2double *constants = nullptr; su2double kine_Inf = 0.0, omega_Inf = 0.0; @@ -1421,9 +1391,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol numerics[iMGlevel][TEMPLATE_SOL][conv_term] = new CConvective_Template(nDim, nVar_Template, config); break; default: - SU2_OMP_MASTER SU2_MPI::Error("Convective scheme not implemented (template_solver).", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1448,9 +1416,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Flow()) { case NO_CONVECTIVE : - SU2_OMP_MASTER SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_FLOW option.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; case SPACE_CENTERED : @@ -1468,9 +1434,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol case LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLaxInc_Flow(nDim, nVar_Flow, config); break; case JST : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentJSTInc_Flow(nDim, nVar_Flow, config); break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid centered scheme or not implemented.\n Currently, only JST and LAX-FRIEDRICH are available for incompressible flows.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } for (iMGlevel = 1; iMGlevel <= config->GetnMGLevels(); iMGlevel++) @@ -1587,9 +1551,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid upwind scheme or not implemented.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1604,18 +1566,14 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid upwind scheme or not implemented.\n Currently, only FDS is available for incompressible flows.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } } break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid convective scheme for the Euler / Navier-Stokes equations.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1714,9 +1672,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Flow()) { case NO_CONVECTIVE : - SU2_OMP_MASTER SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_FLOW option.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; case SPACE_CENTERED : @@ -1725,9 +1681,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol switch (config->GetKind_Centered_Flow()) { case LAX : numerics[MESH_0][FLOW_SOL][conv_term] = new CCentLax_NEMO(nDim, nVar_NEMO, nPrimVar_NEMO, nPrimVarGrad_NEMO, config); break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid centered scheme or not implemented.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1779,9 +1733,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid upwind scheme or not implemented.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1789,9 +1741,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid convective scheme for the NEMO Euler / Navier-Stokes equations.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1861,9 +1811,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol break; default: - SU2_OMP_MASTER SU2_MPI::Error("Riemann solver not implemented.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1877,9 +1825,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol switch (config->GetKind_ConvNumScheme_Turb()) { case NO_UPWIND: - SU2_OMP_MASTER SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; case SPACE_UPWIND : for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { @@ -1890,9 +1836,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid convective scheme for the turbulence equations.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1942,9 +1886,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_Turb()) { case NO_UPWIND: - SU2_OMP_MASTER SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_TURB option.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; case SPACE_UPWIND: for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { @@ -1952,9 +1894,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid convective scheme for the transition equations.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -1997,9 +1937,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid convective scheme for the heat transfer equations.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } } @@ -2023,17 +1961,13 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (adj_euler || adj_ns) { if (incompressible) - SU2_OMP_MASTER SU2_MPI::Error("Convective schemes not implemented for incompressible continuous adjoint.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_AdjFlow()) { case NO_CONVECTIVE: - SU2_OMP_MASTER SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_ADJFLOW option.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; case SPACE_CENTERED : @@ -2046,9 +1980,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol case LAX : numerics[MESH_0][ADJFLOW_SOL][conv_term] = new CCentLax_AdjFlow(nDim, nVar_Adj_Flow, config); break; case JST : numerics[MESH_0][ADJFLOW_SOL][conv_term] = new CCentJST_AdjFlow(nDim, nVar_Adj_Flow, config); break; default: - SU2_OMP_MASTER SU2_MPI::Error("Centered scheme not implemented.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -2075,18 +2007,14 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol } break; default: - SU2_OMP_MASTER SU2_MPI::Error("Upwind scheme not implemented.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } } break; default: - SU2_OMP_MASTER SU2_MPI::Error("Invalid convective scheme for the continuous adjoint Euler / Navier-Stokes equations.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } @@ -2148,25 +2076,19 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol if (adj_turb) { if (!spalart_allmaras) - SU2_OMP_MASTER SU2_MPI::Error("Only the SA turbulence model can be used with the continuous adjoint solver.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER /*--- Definition of the convective scheme for each equation and mesh level ---*/ switch (config->GetKind_ConvNumScheme_AdjTurb()) { case NO_CONVECTIVE: - SU2_OMP_MASTER SU2_MPI::Error("Config file is missing the CONV_NUM_METHOD_ADJTURB option.", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; case SPACE_UPWIND : for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) numerics[iMGlevel][ADJTURB_SOL][conv_term] = new CUpwSca_AdjTurb(nDim, nVar_Adj_Turb, config); break; default: - SU2_OMP_MASTER SU2_MPI::Error("Convective scheme not implemented (adjoint turbulence).", CURRENT_FUNCTION); - END_SU2_OMP_MASTER break; } diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 62d3b2c7cf8..8a8d5ad1461 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -502,7 +502,14 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig /*--- Clear residual (loses AD info), we do not want an incremental solution. ---*/ SU2_OMP_PARALLEL { LinSysRes.SetValZero(); - if (time_domain && config->GetFSI_Simulation()) LinSysSol.SetValZero(); + + if (time_domain && config->GetFSI_Simulation()) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) + for (unsigned short iDim = 0; iDim < nDim; ++iDim) + LinSysSol(iPoint, iDim) = nodes->GetSolution(iPoint, iDim); + END_SU2_OMP_FOR + } } END_SU2_OMP_PARALLEL @@ -593,10 +600,18 @@ void CMeshSolver::ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumeri AD::EndPassive(wasActive); + const su2double velRef = config->GetVelocity_Ref(); + const su2double invVelRef = 1.0 / velRef; + /*--- Clear residual (loses AD info), we do not want an incremental solution. ---*/ SU2_OMP_PARALLEL { LinSysRes.SetValZero(); - LinSysSol.SetValZero(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + for (unsigned short iDim = 0; iDim < nDim; iDim++) + LinSysSol(iPoint, iDim) = geometry[MESH_0]->nodes->GetGridVel(iPoint)[iDim] * velRef; + END_SU2_OMP_FOR } END_SU2_OMP_PARALLEL @@ -607,16 +622,9 @@ void CMeshSolver::ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumeri Solve_System(geometry[MESH_0], config); SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - for (unsigned short iDim = 0; iDim < nDim; iDim++) { - su2double val_vel = LinSysSol(iPoint, iDim); - - /*--- Non-dimensionalize velocity ---*/ - val_vel /= config->GetVelocity_Ref(); - - geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, val_vel); - } - } + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + for (unsigned short iDim = 0; iDim < nDim; iDim++) + geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, LinSysSol(iPoint,iDim)*invVelRef); END_SU2_OMP_PARALLEL } From 6364cbf235c950023b558b5bef942457c4ffb420 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 24 May 2021 14:31:04 +0100 Subject: [PATCH 11/21] less duplication, better unsteady FSI restarts --- Common/include/geometry/CGeometry.hpp | 2 +- Common/include/geometry/dual_grid/CPoint.hpp | 2 +- .../include/solvers/CFVMFlowSolverBase.inl | 69 +++++------ SU2_CFD/include/solvers/CMeshSolver.hpp | 30 ++--- SU2_CFD/include/solvers/CSolver.hpp | 2 +- SU2_CFD/src/drivers/CDriver.cpp | 47 ++++---- SU2_CFD/src/drivers/CMultizoneDriver.cpp | 2 +- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 2 +- .../src/iteration/CDiscAdjFluidIteration.cpp | 2 +- .../src/iteration/CDiscAdjHeatIteration.cpp | 2 +- SU2_CFD/src/output/CFlowCompOutput.cpp | 9 +- SU2_CFD/src/output/CFlowIncOutput.cpp | 11 +- SU2_CFD/src/output/CNEMOCompOutput.cpp | 14 ++- SU2_CFD/src/solvers/CMeshSolver.cpp | 114 ++++++------------ 14 files changed, 134 insertions(+), 174 deletions(-) diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index ad4f6d80e67..c0d1ac313ba 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1258,7 +1258,7 @@ class CGeometry { * \param geometry_container - Geometrical definition. * \param config - Config */ - void UpdateGeometry(CGeometry **geometry_container, CConfig *config); + static void UpdateGeometry(CGeometry **geometry_container, CConfig *config); /*! * \brief Update the multi-grid structure for the customized boundary conditions diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 83f1b89b7ad..56a21fc1ad9 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -543,7 +543,7 @@ class CPoint { void SetVolume_n(); /*! - * \brief Set the volume of the control volume at time n+1. + * \brief Set the volume of the control volume at time n-1. */ void SetVolume_nM1(); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 64d3934c539..374a6ec3917 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -707,7 +707,7 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** /*--- Restart the solution from file information ---*/ - unsigned short iDim, iVar, iMesh, iMeshFine; + unsigned short iDim, iVar, iMesh; unsigned long iPoint, index, iChildren, Point_Fine; unsigned short turb_model = config->GetKind_Turb_Model(); su2double Area_Children, Area_Parent; @@ -746,6 +746,15 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename); } + if (update_geo && dynamic_grid) { + auto notFound = fields.end(); + if (find(fields.begin(), notFound, string("\"Grid_Velocity_x\"")) == notFound) { + if (rank == MASTER_NODE) + cout << "\nWARNING: The restart file does not contain grid velocities, these will be set to zero.\n" << endl; + steady_restart = true; + } + } + /*--- Load data from the restart into correct containers. ---*/ unsigned long counter = 0, iPoint_Global = 0; @@ -829,39 +838,26 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** END_SU2_OMP_MASTER SU2_OMP_BARRIER - /*--- Update the geometry for flows on deforming meshes ---*/ + /*--- Update the geometry for flows on deforming meshes. ---*/ if ((dynamic_grid || static_fsi) && update_geo) { - /*--- Communicate the new coordinates and grid velocities at the halos ---*/ - - geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, COORDINATES); - geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, COORDINATES); + CGeometry::UpdateGeometry(geometry, config); if (dynamic_grid) { - geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY); - geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY); - } - - /*--- Recompute the edges and dual mesh control volumes in the - domain and on the boundaries. ---*/ + for (iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) { - geometry[MESH_0]->SetControlVolume(config, UPDATE); - geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); - geometry[MESH_0]->SetMaxLength(config); + /*--- Compute the grid velocities on the coarser levels. ---*/ + if (iMesh) geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMesh-1], config); - /*--- Update the multigrid structure after setting up the finest grid, - including computing the grid velocities on the coarser levels. ---*/ + geometry[iMesh]->nodes->SetVolume_n(); + geometry[iMesh]->nodes->SetVolume_nM1(); - for (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { - iMeshFine = iMesh-1; - geometry[iMesh]->SetControlVolume(config, geometry[iMeshFine], UPDATE); - geometry[iMesh]->SetBoundControlVolume(config, geometry[iMeshFine],UPDATE); - geometry[iMesh]->SetCoord(geometry[iMeshFine]); - if (dynamic_grid) { - geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMeshFine], config); + if (config->GetGrid_Movement()) { + geometry[iMesh]->nodes->SetCoord_n(); + geometry[iMesh]->nodes->SetCoord_n1(); + } } - geometry[iMesh]->SetMaxLength(config); } } @@ -973,7 +969,7 @@ void CFVMFlowSolverBase::SetInitialCondition(CGeometry **geometry, CSolver /*--- The value of the solution for the first iteration of the dual time ---*/ - if (dual_time && (TimeIter == 0 || (restart && TimeIter == config->GetRestart_Iter()))) { + if (dual_time && TimeIter == config->GetRestart_Iter()) { PushSolutionBackInTime(TimeIter, restart, rans, solver_container, geometry, config); } @@ -998,27 +994,24 @@ void CFVMFlowSolverBase::PushSolutionBackInTime(unsigned long TimeIter, bo } } - if (restart && (TimeIter == config->GetRestart_Iter()) && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { - - const bool update_geo = !config->GetFSI_Simulation(); + if (restart && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { - /*--- Load an additional restart file for a 2nd-order restart ---*/ + /*--- Load an additional restart file for a 2nd-order restart. ---*/ - solver_container[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver_container, config, config->GetRestart_Iter() - 1, - update_geo); + solver_container[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver_container, config, TimeIter-1, true); - /*--- Load an additional restart file for the turbulence model ---*/ + /*--- Load an additional restart file for the turbulence model. ---*/ if (rans) - solver_container[MESH_0][TURB_SOL]->LoadRestart(geometry, solver_container, config, config->GetRestart_Iter() - 1, - false); + solver_container[MESH_0][TURB_SOL]->LoadRestart(geometry, solver_container, config, TimeIter-1, false); /*--- Push back this new solution to time level N. ---*/ for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) { solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n(); - if (rans) { - solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n(); - } + if (rans) solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n(); + + geometry[iMesh]->nodes->SetVolume_n(); + if (config->GetGrid_Movement()) geometry[iMesh]->nodes->SetCoord_n(); } } } diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index 5a4a6d927f0..27fd06dd39c 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -72,21 +72,14 @@ class CMeshSolver final : public CFEASolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void UpdateGridCoord(CGeometry *geometry, CConfig *config); - - /*! - * \brief Update the dual grid after the grid movement (edges and control volumes). - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void UpdateDualGrid(CGeometry *geometry, CConfig *config); + void UpdateGridCoord(CGeometry *geometry, const CConfig *config); /*! * \brief Compute the grid velocity form the displacements of the mesh. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void ComputeGridVelocity(CGeometry *geometry, CConfig *config); + void ComputeGridVelocity(CGeometry **geometry, const CConfig *config) const; /*! * \brief Compute the grid velocity form the velocity at deformable boundary. @@ -96,13 +89,6 @@ class CMeshSolver final : public CFEASolver { */ void ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumerics **numerics, CConfig *config); - /*! - * \brief Update the coarse multigrid levels after the grid movement. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void UpdateMultiGrid(CGeometry **geometry, CConfig *config) const; - /*! * \brief Check the boundary vertex that are going to be moved. * \param[in] geometry - Geometrical definition of the problem. @@ -120,6 +106,12 @@ class CMeshSolver final : public CFEASolver { */ void BC_Deforming(CGeometry *geometry, const CConfig *config, unsigned short val_marker, bool velocity); + /*! + * \brief Load the geometries at the previous time states n and nM1. + * \param[in] geometry - Geometrical definition of the problem. + */ + void RestartOldGeometry(CGeometry *geometry, const CConfig *config); + public: /*! * \brief Constructor of the class. @@ -174,12 +166,6 @@ class CMeshSolver final : public CFEASolver { int val_iter, bool val_update_geo) override; - /*! - * \brief Load the geometries at the previous time states n and nM1. - * \param[in] geometry - Geometrical definition of the problem. - */ - void Restart_OldGeometry(CGeometry *geometry, CConfig *config) override; - /*! * \brief Get minimun volume in the mesh * \return diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 946c4e36502..ccf0bbf7408 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -630,7 +630,7 @@ class CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - virtual void Restart_OldGeometry(CGeometry *geometry, CConfig *config); + void Restart_OldGeometry(CGeometry *geometry, CConfig *config); /*! * \brief A virtual member. diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 25f56322426..ee98c0b25a3 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -229,14 +229,9 @@ CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunica CGeometry::ComputeWallDistance(config_container, geometry_container); } + /*--- Definition of the interface and transfer conditions between different zones. ---*/ - /*--- Definition of the interface and transfer conditions between different zones. - *--- The transfer container is defined for zones paired one to one. - *--- This only works for a multizone FSI problem (nZone > 1). - *--- Also, at the moment this capability is limited to two zones (nZone < 3). - *--- This will change in the future. ---*/ - - if ( nZone > 1 ) { + if (nZone > 1) { if (rank == MASTER_NODE) cout << endl <<"------------------- Multizone Interface Preprocessing -------------------" << endl; @@ -244,8 +239,19 @@ CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunica interface_types, interface_container, interpolator_container); } + if (fsi) { + for (iZone = 0; iZone < nZone; iZone++) { + for (iInst = 0; iInst < nInst[iZone]; iInst++){ + Solver_Restart(solver_container[iZone][iInst], geometry_container[iZone][iInst], + config_container[iZone], true); + } + } + } + if (config_container[ZONE_0]->GetBoolTurbomachinery()){ - if (rank == MASTER_NODE)cout << endl <<"---------------------- Turbomachinery Preprocessing ---------------------" << endl; + if (rank == MASTER_NODE) + cout << endl <<"---------------------- Turbomachinery Preprocessing ---------------------" << endl; + Turbomachinery_Preprocessing(config_container, geometry_container, solver_container, interface_container); } @@ -861,16 +867,18 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet /*--- For unsteady simulations, initialize the grid volumes and coordinates for previous solutions. Loop over all zones/grids ---*/ - if ((config->GetTime_Marching() != TIME_MARCHING::STEADY) && config->GetGrid_Movement()) { + if ((config->GetTime_Marching() != TIME_MARCHING::STEADY) && config->GetDynamic_Grid()) { for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { /*--- Update cell volume ---*/ geometry[iMGlevel]->nodes->SetVolume_n(); geometry[iMGlevel]->nodes->SetVolume_nM1(); - /*--- Update point coordinates ---*/ - geometry[iMGlevel]->nodes->SetCoord_n(); - geometry[iMGlevel]->nodes->SetCoord_n1(); + if (config->GetGrid_Movement()) { + /*--- Update point coordinates ---*/ + geometry[iMGlevel]->nodes->SetCoord_n(); + geometry[iMGlevel]->nodes->SetCoord_n1(); + } } } @@ -1009,16 +1017,13 @@ void CDriver::Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolve DOFsPerPoint = 0; - for (unsigned int iSol = 0; iSol < MAX_SOLS; iSol++){ - if (solver[MESH_0][iSol] != nullptr){ - DOFsPerPoint += solver[MESH_0][iSol]->GetnVar(); - } - } + for (unsigned int iSol = 0; iSol < MAX_SOLS; iSol++) + if (solver[MESH_0][iSol]) DOFsPerPoint += solver[MESH_0][iSol]->GetnVar(); - bool update_geo = true; - if (config->GetFSI_Simulation()) update_geo = false; + /*--- Restart solvers, for FSI the geometry cannot be updated because the interpolation classes + * should always use the undeformed mesh (otherwise the results would not be repeatable). ---*/ - Solver_Restart(solver, geometry, config, update_geo); + if (!fsi) Solver_Restart(solver, geometry, config, true); /*--- Set up any necessary inlet profiles ---*/ @@ -2251,7 +2256,7 @@ void CDriver::DynamicMesh_Preprocessing(CConfig *config, CGeometry **geometry, C /*--- Update the multi-grid structure to propagate the derivative information to the coarser levels ---*/ - geometry[MESH_0]->UpdateGeometry(geometry,config); + CGeometry::UpdateGeometry(geometry,config); } diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 707f5bf264d..314f21daace 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -244,7 +244,7 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) { solver_container[iZone][INST_0], config_container[iZone], TimeIter); } - else if (!fsi && !config_container[iZone]->GetDiscrete_Adjoint() && config_container[iZone]->GetHeatProblem()) { + else if (!fsi && config_container[iZone]->GetHeatProblem()) { /*--- Set the initial condition for HEAT equation ---------------------------------------------*/ solver_container[iZone][INST_0][MESH_0][HEAT_SOL]->SetInitialCondition(geometry_container[iZone][INST_0], solver_container[iZone][INST_0], diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index 05090c168c3..2b9f488d493 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -129,7 +129,7 @@ void CSinglezoneDriver::Preprocess(unsigned long TimeIter) { solver_container[ZONE_0][INST_0], config_container[ZONE_0], TimeIter); } - else if ( config_container[ZONE_0]->GetHeatProblem()) { + else if (config_container[ZONE_0]->GetHeatProblem()) { /*--- Set the initial condition for HEAT equation ---------------------------------------------*/ solver_container[ZONE_0][INST_0][MESH_0][HEAT_SOL]->SetInitialCondition(geometry_container[ZONE_0][INST_0], solver_container[ZONE_0][INST_0], diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 8cfd548399c..87449a46cc4 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -503,7 +503,7 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver***** solver, CGeometry**** (kind_recording == RECORDING::SOLUTION_AND_MESH)) { /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ - geometry[iZone][iInst][MESH_0]->UpdateGeometry(geometry[iZone][iInst], config[iZone]); + CGeometry::UpdateGeometry(geometry[iZone][iInst], config[iZone]); CGeometry::ComputeWallDistance(config, geometry); } diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp index aa6f9471a01..be710ab4b2e 100644 --- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp @@ -212,7 +212,7 @@ void CDiscAdjHeatIteration::SetDependencies(CSolver***** solver, CGeometry**** g (kind_recording == RECORDING::SOLUTION_AND_MESH)) { /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ - geometries[MESH_0]->UpdateGeometry(geometries, config[iZone]); + CGeometry::UpdateGeometry(geometries, config[iZone]); CGeometry::ComputeWallDistance(config, geometry); } diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index 077aa5b988f..5f053a80275 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -58,10 +58,17 @@ CFlowCompOutput::CFlowCompOutput(CConfig *config, unsigned short nDim) : CFlowOu requestedVolumeFields.emplace_back("COORDINATES"); requestedVolumeFields.emplace_back("SOLUTION"); requestedVolumeFields.emplace_back("PRIMITIVE"); - if (gridMovement) requestedVolumeFields.emplace_back("GRID_VELOCITY"); nRequestedVolumeFields = requestedVolumeFields.size(); } + if (gridMovement) { + auto notFound = requestedVolumeFields.end(); + if (find(requestedVolumeFields.begin(), notFound, string("GRID_VELOCITY")) == notFound) { + requestedVolumeFields.emplace_back("GRID_VELOCITY"); + nRequestedVolumeFields ++; + } + } + stringstream ss; ss << "Zone " << config->GetiZone() << " (Comp. Fluid)"; multiZoneHeaderString = ss.str(); diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index b1b23e74c42..229d61f866c 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -64,10 +64,17 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp requestedVolumeFields.emplace_back("COORDINATES"); requestedVolumeFields.emplace_back("SOLUTION"); requestedVolumeFields.emplace_back("PRIMITIVE"); - if (gridMovement) requestedVolumeFields.emplace_back("GRID_VELOCITY"); nRequestedVolumeFields = requestedVolumeFields.size(); } + if (gridMovement) { + auto notFound = requestedVolumeFields.end(); + if (find(requestedVolumeFields.begin(), notFound, string("GRID_VELOCITY")) == notFound) { + requestedVolumeFields.emplace_back("GRID_VELOCITY"); + nRequestedVolumeFields ++; + } + } + stringstream ss; ss << "Zone " << config->GetiZone() << " (Incomp. Fluid)"; multiZoneHeaderString = ss.str(); @@ -436,7 +443,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("RES_VELOCITY-Y", "Residual_Velocity_y", "RESIDUAL", "Residual of the y-velocity component"); if (nDim == 3) AddVolumeOutput("RES_VELOCITY-Z", "Residual_Velocity_z", "RESIDUAL", "Residual of the z-velocity component"); - if (config->GetEnergy_Equation()) + if (config->GetEnergy_Equation()) AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature"); switch(config->GetKind_Turb_Model()){ diff --git a/SU2_CFD/src/output/CNEMOCompOutput.cpp b/SU2_CFD/src/output/CNEMOCompOutput.cpp index c3de5fbf3a5..291127f4e08 100644 --- a/SU2_CFD/src/output/CNEMOCompOutput.cpp +++ b/SU2_CFD/src/output/CNEMOCompOutput.cpp @@ -44,7 +44,7 @@ CNEMOCompOutput::CNEMOCompOutput(CConfig *config, unsigned short nDim) : CFlowOu turb_model = config->GetKind_Turb_Model(); lastInnerIter = curInnerIter; - gridMovement = config->GetGrid_Movement(); + gridMovement = config->GetDynamic_Grid(); nSpecies = config->GetnSpecies(); /*--- Set the default history fields if nothing is set in the config file ---*/ @@ -74,6 +74,14 @@ CNEMOCompOutput::CNEMOCompOutput(CConfig *config, unsigned short nDim) : CFlowOu nRequestedVolumeFields = requestedVolumeFields.size(); } + if (gridMovement) { + auto notFound = requestedVolumeFields.end(); + if (find(requestedVolumeFields.begin(), notFound, string("GRID_VELOCITY")) == notFound) { + requestedVolumeFields.emplace_back("GRID_VELOCITY"); + nRequestedVolumeFields ++; + } + } + stringstream ss; ss << "Zone " << config->GetiZone() << " (Comp. Fluid)"; multiZoneHeaderString = ss.str(); @@ -327,7 +335,7 @@ void CNEMOCompOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("MASSFRAC_" + std::to_string(iSpecies), "MassFrac_" + std::to_string(iSpecies), "AUXILIARY", "MassFrac_" + std::to_string(iSpecies)); // Grid velocity - if (config->GetDynamic_Grid()){ + if (gridMovement){ AddVolumeOutput("GRID_VELOCITY-X", "Grid_Velocity_x", "GRID_VELOCITY", "x-component of the grid velocity vector"); AddVolumeOutput("GRID_VELOCITY-Y", "Grid_Velocity_y", "GRID_VELOCITY", "y-component of the grid velocity vector"); if (nDim == 3 ) @@ -474,7 +482,7 @@ void CNEMOCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv break; } - if (config->GetDynamic_Grid()){ + if (gridMovement){ SetVolumeOutputValue("GRID_VELOCITY-X", iPoint, Node_Geo->GetGridVel(iPoint)[0]); SetVolumeOutputValue("GRID_VELOCITY-Y", iPoint, Node_Geo->GetGridVel(iPoint)[1]); if (nDim == 3) diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 8a8d5ad1461..e8f68c9a3af 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -526,14 +526,14 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig UpdateGridCoord(geometry[MESH_0], config); /*--- Update the dual grid. ---*/ - UpdateDualGrid(geometry[MESH_0], config); + CGeometry::UpdateGeometry(geometry, config); /*--- Check for failed deformation (negative volumes). ---*/ SetMinMaxVolume(geometry[MESH_0], config, true); /*--- The Grid Velocity is only computed if the problem is time domain ---*/ if (time_domain && !config->GetFSI_Simulation()) - ComputeGridVelocity(geometry[MESH_0], config); + ComputeGridVelocity(geometry, config); } END_SU2_OMP_PARALLEL @@ -542,14 +542,9 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig ComputeGridVelocity_FromBoundary(geometry, numerics, config); } - /*--- Update the multigrid structure. ---*/ - SU2_OMP_PARALLEL - UpdateMultiGrid(geometry, config); - END_SU2_OMP_PARALLEL - } -void CMeshSolver::UpdateGridCoord(CGeometry *geometry, CConfig *config){ +void CMeshSolver::UpdateGridCoord(CGeometry *geometry, const CConfig *config){ /*--- Update the grid coordinates using the solution of the linear system ---*/ @@ -575,17 +570,6 @@ void CMeshSolver::UpdateGridCoord(CGeometry *geometry, CConfig *config){ } -void CMeshSolver::UpdateDualGrid(CGeometry *geometry, CConfig *config){ - - /*--- After moving all nodes, update the dual mesh. Recompute the edges and - dual mesh control volumes in the domain and on the boundaries. ---*/ - - geometry->SetControlVolume(config, UPDATE); - geometry->SetBoundControlVolume(config, UPDATE); - geometry->SetMaxLength(config); - -} - void CMeshSolver::ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumerics **numerics, CConfig *config){ if (config->GetnZone() == 1) @@ -621,20 +605,30 @@ void CMeshSolver::ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumeri /*--- Solve the linear system. ---*/ Solve_System(geometry[MESH_0], config); - SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) - for (unsigned short iDim = 0; iDim < nDim; iDim++) - geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, LinSysSol(iPoint,iDim)*invVelRef); + SU2_OMP_PARALLEL { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) + for (unsigned short iDim = 0; iDim < nDim; iDim++) + geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, LinSysSol(iPoint,iDim)*invVelRef); + END_SU2_OMP_FOR + + for (auto iMGlevel = 1u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + geometry[iMGlevel]->SetRestricted_GridVelocity(geometry[iMGlevel-1], config); + } END_SU2_OMP_PARALLEL + } -void CMeshSolver::ComputeGridVelocity(CGeometry *geometry, CConfig *config){ +void CMeshSolver::ComputeGridVelocity(CGeometry **geometry, const CConfig *config) const { + + /*--- Compute the velocity of each node. ---*/ - /*--- Compute the velocity of each node in the domain of the current rank - (halo nodes are not computed as the grid velocity is later communicated). ---*/ + const bool firstOrder = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST; + const bool secondOrder = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND; + const su2double invTimeStep = 1.0 / config->GetDelta_UnstTimeND(); SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { /*--- Coordinates of the current point at n+1, n, & n-1 time levels. ---*/ @@ -642,50 +636,23 @@ void CMeshSolver::ComputeGridVelocity(CGeometry *geometry, CConfig *config){ const su2double* Disp_n = nodes->GetSolution_time_n(iPoint); const su2double* Disp_nP1 = nodes->GetSolution(iPoint); - /*--- Unsteady time step ---*/ - - su2double TimeStep = config->GetDelta_UnstTimeND(); - - /*--- Compute mesh velocity with 1st or 2nd-order approximation. ---*/ + /*--- Compute mesh velocity for this point with 1st or 2nd-order approximation. ---*/ for (unsigned short iDim = 0; iDim < nDim; iDim++) { su2double GridVel = 0.0; + if (firstOrder) + GridVel = (Disp_nP1[iDim] - Disp_n[iDim]) * invTimeStep; + else if (secondOrder) + GridVel = (1.5*Disp_nP1[iDim] - 2.0*Disp_n[iDim] + 0.5*Disp_nM1[iDim]) * invTimeStep; - if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) - GridVel = ( Disp_nP1[iDim] - Disp_n[iDim] ) / TimeStep; - if (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) - GridVel = ( 3.0*Disp_nP1[iDim] - 4.0*Disp_n[iDim] + - 1.0*Disp_nM1[iDim] ) / (2.0*TimeStep); - - /*--- Store grid velocity for this point ---*/ - - geometry->nodes->SetGridVel(iPoint, iDim, GridVel); - + geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, GridVel); } } END_SU2_OMP_FOR - /*--- The velocity was computed for nPointDomain, now we communicate it. ---*/ - geometry->InitiateComms(geometry, config, GRID_VELOCITY); - geometry->CompleteComms(geometry, config, GRID_VELOCITY); - -} - -void CMeshSolver::UpdateMultiGrid(CGeometry **geometry, CConfig *config) const{ - - /*--- Update the multigrid structure after moving the finest grid, - including computing the grid velocities on the coarser levels - when the problem is solved in unsteady conditions. ---*/ - - for (auto iMGlevel = 1u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) { - const auto iMGfine = iMGlevel-1; - geometry[iMGlevel]->SetControlVolume(config, geometry[iMGfine], UPDATE); - geometry[iMGlevel]->SetBoundControlVolume(config, geometry[iMGfine],UPDATE); - geometry[iMGlevel]->SetCoord(geometry[iMGfine]); - if (time_domain) - geometry[iMGlevel]->SetRestricted_GridVelocity(geometry[iMGfine], config); - } + for (auto iMGlevel = 1u; iMGlevel <= config->GetnMGLevels(); iMGlevel++) + geometry[iMGlevel]->SetRestricted_GridVelocity(geometry[iMGlevel-1], config); } @@ -857,10 +824,6 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * solver[MESH_0][MESH_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION); solver[MESH_0][MESH_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION); - /*--- Communicate the new coordinates at the halos ---*/ - geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, COORDINATES); - geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, COORDINATES); - /*--- Init the linear system solution. ---*/ for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { for (unsigned short iDim = 0; iDim < nDim; ++iDim) { @@ -868,23 +831,14 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * } } - /*--- Recompute the edges and dual mesh control volumes in the - domain and on the boundaries. ---*/ - UpdateDualGrid(geometry[MESH_0], config); - /*--- For time-domain problems, we need to compute the grid velocities. ---*/ if (time_domain && !config->GetFSI_Simulation()) { /*--- Update the old geometry (coordinates n and n-1) ---*/ - Restart_OldGeometry(geometry[MESH_0], config); - /*--- Once Displacement_n and Displacement_n1 are filled, - we can compute the Grid Velocity ---*/ - ComputeGridVelocity(geometry[MESH_0], config); - } + RestartOldGeometry(geometry[MESH_0], config); - /*--- Update the multigrid structure after setting up the finest grid, - including computing the grid velocities on the coarser levels - when the problem is unsteady. ---*/ - UpdateMultiGrid(geometry, config); + /*--- Once Displacement_n and Displacement_n1 are filled we can compute the Grid Velocity ---*/ + ComputeGridVelocity(geometry, config); + } /*--- Store the boundary displacements at the Bound_Disp variable. ---*/ @@ -911,7 +865,7 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * } -void CMeshSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { +void CMeshSolver::RestartOldGeometry(CGeometry *geometry, const CConfig *config) { /*--- This function is intended for dual time simulations ---*/ From d2fbc1357df49cc8c5c38303365d7afa4b021c6e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 24 May 2021 15:38:14 +0100 Subject: [PATCH 12/21] fix for turbo --- Common/src/geometry/CGeometry.cpp | 4 ---- SU2_CFD/include/solvers/CFVMFlowSolverBase.inl | 16 +++++++++++----- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index d4aab89d4c0..a7b1fe7b9fa 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -2489,10 +2489,6 @@ void CGeometry::UpdateGeometry(CGeometry **geometry_container, CConfig *config) geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, COORDINATES); geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, COORDINATES); - if (config->GetDynamic_Grid()){ - geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, GRID_VELOCITY); - geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, GRID_VELOCITY); - } geometry_container[MESH_0]->SetControlVolume(config, UPDATE); geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 374a6ec3917..67c1366fb7d 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -849,13 +849,19 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** /*--- Compute the grid velocities on the coarser levels. ---*/ if (iMesh) geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMesh-1], config); + else { + geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY); + geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY); + } - geometry[iMesh]->nodes->SetVolume_n(); - geometry[iMesh]->nodes->SetVolume_nM1(); + if (config->GetTime_Marching() != TIME_MARCHING::STEADY) { + geometry[iMesh]->nodes->SetVolume_n(); + geometry[iMesh]->nodes->SetVolume_nM1(); - if (config->GetGrid_Movement()) { - geometry[iMesh]->nodes->SetCoord_n(); - geometry[iMesh]->nodes->SetCoord_n1(); + if (config->GetGrid_Movement()) { + geometry[iMesh]->nodes->SetCoord_n(); + geometry[iMesh]->nodes->SetCoord_n1(); + } } } } From 319e51e67a11aa3dbecdc0b19b1e396e98cc79df Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 24 May 2021 17:22:09 +0100 Subject: [PATCH 13/21] fix for moving grid disc adj --- .../include/solvers/CFVMFlowSolverBase.inl | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 67c1366fb7d..5102b521a13 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -853,16 +853,6 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver ** geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY); geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY); } - - if (config->GetTime_Marching() != TIME_MARCHING::STEADY) { - geometry[iMesh]->nodes->SetVolume_n(); - geometry[iMesh]->nodes->SetVolume_nM1(); - - if (config->GetGrid_Movement()) { - geometry[iMesh]->nodes->SetCoord_n(); - geometry[iMesh]->nodes->SetCoord_n1(); - } - } } } } @@ -998,6 +988,16 @@ void CFVMFlowSolverBase::PushSolutionBackInTime(unsigned long TimeIter, bo solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n(); solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n1(); } + + if (dynamic_grid) { + geometry[iMesh]->nodes->SetVolume_n(); + geometry[iMesh]->nodes->SetVolume_nM1(); + } + + if (config->GetGrid_Movement()) { + geometry[iMesh]->nodes->SetCoord_n(); + geometry[iMesh]->nodes->SetCoord_n1(); + } } if (restart && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { From 6094cb6d753e0c9b6a83f9a0b5293d011e5cc1eb Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 24 May 2021 23:28:21 +0100 Subject: [PATCH 14/21] update tests --- TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref | 16 ++++++++-------- TestCases/hybrid_regression.py | 4 ++-- TestCases/parallel_regression.py | 8 ++++---- TestCases/serial_regression.py | 8 ++++---- 4 files changed, 18 insertions(+), 18 deletions(-) diff --git a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref index 9f9eb61496d..72cc7d70416 100644 --- a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref +++ b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref @@ -1,9 +1,9 @@ INDEX GRAD -0 -3.465631559336758e-03 -1 -1.844003066470726e-03 -2 -7.925069065625049e-04 -3 -2.742895537614673e-04 -4 -2.738092746211827e-04 -5 -7.890546914650276e-04 -6 -1.831172635608937e-03 -7 -3.431364694518416e-03 +0 -3.465632121754074e-03 +1 -1.844003360219930e-03 +2 -7.925070240684198e-04 +3 -2.742895873196515e-04 +4 -2.738093165742900e-04 +5 -7.890548338696854e-04 +6 -1.831172969426724e-03 +7 -3.431365309796604e-03 diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 48b223f8024..af6011b0eaf 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -646,7 +646,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4, 0, -3.771897, -4.163473] + fsi2d.test_vals = [4, 0, -3.743230, -4.133462] fsi2d.multizone= True fsi2d.unsteady = True test_list.append(fsi2d) @@ -665,7 +665,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.400852, -4.085149, 5.3839e-08, 121] + dyn_fsi.test_vals = [-4.355806, -4.060581, 5.3837e-08, 100] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index abee592c8eb..54728bb237a 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1175,7 +1175,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4, 0, -3.771873, -4.163491] #last 4 columns + fsi2d.test_vals = [4, 0, -3.743210, -4.133483] #last 4 columns fsi2d.su2_exec = "parallel_computation_fsi.py -f" fsi2d.timeout = 1600 fsi2d.multizone= True @@ -1200,7 +1200,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.400874, -4.085152, 0.000000, 117] + dyn_fsi.test_vals = [-4.355829, -4.060587, 5.3837e-08, 98] dyn_fsi.multizone = True dyn_fsi.unsteady = True dyn_fsi.su2_exec = "mpirun -n 2 SU2_CFD" @@ -1356,7 +1356,7 @@ def main(): pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" pywrapper_fsi2d.cfg_file = "configFSI.cfg" pywrapper_fsi2d.test_iter = 4 - pywrapper_fsi2d.test_vals = [4, 0, -3.771873, -4.163491] #last 4 columns + pywrapper_fsi2d.test_vals = [4, 0, -3.743210, -4.133483] #last 4 columns pywrapper_fsi2d.su2_exec = "mpirun -np 2 SU2_CFD.py --nZone 2 --fsi True --parallel -f" pywrapper_fsi2d.timeout = 1600 pywrapper_fsi2d.unsteady = True @@ -1382,7 +1382,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.551335, 2.295594, 0.350036, 0.093081] + pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350036, 0.093137] pywrapper_rigidMotion.su2_exec = "mpirun -np 2 python launch_flatPlate_rigidMotion.py --parallel -f" pywrapper_rigidMotion.timeout = 1600 pywrapper_rigidMotion.tol = 0.00001 diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index feb426f565e..6a715ac5478 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1319,7 +1319,7 @@ def main(): fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" fsi2d.cfg_file = "configFSI.cfg" fsi2d.test_iter = 4 - fsi2d.test_vals = [4, 0, -3.771878, -4.163489] #last 4 columns + fsi2d.test_vals = [4, 0, -3.743214, -4.133482] #last 4 columns fsi2d.su2_exec = "SU2_CFD" fsi2d.timeout = 1600 fsi2d.multizone = True @@ -1356,7 +1356,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.400854, -4.085153, 0.000000, 101] #last 4 columns + dyn_fsi.test_vals = [-4.355809, -4.060588, 5.3837e-08, 86] #last 4 columns dyn_fsi.multizone = True dyn_fsi.unsteady = True dyn_fsi.su2_exec = "SU2_CFD" @@ -1874,7 +1874,7 @@ def main(): pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" pywrapper_fsi2d.cfg_file = "configFSI.cfg" pywrapper_fsi2d.test_iter = 4 - pywrapper_fsi2d.test_vals = [4, 0, -3.771878, -4.163489] #last 4 columns + pywrapper_fsi2d.test_vals = [4, 0, -3.743214, -4.133482] #last 4 columns pywrapper_fsi2d.su2_exec = "SU2_CFD.py --nZone 2 --fsi True -f" pywrapper_fsi2d.new_output = True pywrapper_fsi2d.unsteady = True @@ -1903,7 +1903,7 @@ def main(): pywrapper_rigidMotion.cfg_dir = "py_wrapper/flatPlate_rigidMotion" pywrapper_rigidMotion.cfg_file = "flatPlate_rigidMotion_Conf.cfg" pywrapper_rigidMotion.test_iter = 5 - pywrapper_rigidMotion.test_vals = [-1.551335, 2.295594, 0.350050, 0.093081] + pywrapper_rigidMotion.test_vals = [-1.614170, 2.242953, 0.350050, 0.093137] pywrapper_rigidMotion.su2_exec = "python launch_flatPlate_rigidMotion.py -f" pywrapper_rigidMotion.new_output = True pywrapper_rigidMotion.timeout = 1600 From d5715c40a324fa63aa2dc1e08cab840cfdca659e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 25 May 2021 16:16:18 +0100 Subject: [PATCH 15/21] local indices only, no more hidden global "adjointPosition" --- Common/include/basic_types/ad_structure.hpp | 30 ++----------------- .../basic_types/datatype_structure.hpp | 10 ------- Common/include/geometry/CGeometry.hpp | 6 ++-- Common/include/geometry/dual_grid/CPoint.hpp | 2 +- Common/src/basic_types/ad_structure.cpp | 6 ---- SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp | 10 ++----- SU2_CFD/include/variables/CMeshVariable.hpp | 2 +- SU2_CFD/include/variables/CVariable.hpp | 20 +++++++------ SU2_CFD/src/variables/CMeshVariable.cpp | 5 +--- SU2_CFD/src/variables/CVariable.cpp | 10 ++----- TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg | 23 +++++++------- 11 files changed, 33 insertions(+), 91 deletions(-) diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index 620da3246f5..6ecbaaac531 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -271,14 +271,6 @@ namespace AD{ extern ExtFuncHelper* FuncHelper; - /*--- Stores the indices of the input variables (they might be overwritten) ---*/ - - extern std::vector inputValues; - - /*--- Current position inside the adjoint vector ---*/ - - extern int adjointVectorPosition; - extern bool Status; extern bool PreaccActive; @@ -297,24 +289,13 @@ namespace AD{ extern std::vector TapePositions; - extern std::vector localInputValues; - - extern std::vector localOutputValues; - extern codi::PreaccumulationHelper PreaccHelper; /*--- Reference to the tape. ---*/ - FORCEINLINE su2double::TapeType& getGlobalTape() { - return su2double::getGlobalTape(); - } + FORCEINLINE su2double::TapeType& getGlobalTape() {return su2double::getGlobalTape();} - FORCEINLINE void RegisterInput(su2double &data, bool push_index = true) { - AD::getGlobalTape().registerInput(data); - if (push_index) { - inputValues.push_back(data.getGradientData()); - } - } + FORCEINLINE void RegisterInput(su2double &data) {AD::getGlobalTape().registerInput(data);} FORCEINLINE void RegisterOutput(su2double& data) {AD::getGlobalTape().registerOutput(data);} @@ -335,7 +316,6 @@ namespace AD{ opdi::logic->prepareEvaluate(); #endif AD::getGlobalTape().evaluate(); - adjointVectorPosition = 0; } FORCEINLINE void ComputeAdjoint(unsigned short enter, unsigned short leave) { @@ -346,8 +326,6 @@ namespace AD{ #else AD::getGlobalTape().evaluate(TapePositions[enter], TapePositions[leave]); #endif - if (leave == 0) - adjointVectorPosition = 0; } FORCEINLINE void Reset() { @@ -355,10 +333,6 @@ namespace AD{ #if defined(HAVE_OPDI) opdi::logic->reset(); #endif - if (inputValues.size() != 0) { - adjointVectorPosition = 0; - inputValues.clear(); - } if (TapePositions.size() != 0) { #if defined(HAVE_OPDI) for (TapePosition& pos : TapePositions) { diff --git a/Common/include/basic_types/datatype_structure.hpp b/Common/include/basic_types/datatype_structure.hpp index 039df331200..b30d81555f1 100644 --- a/Common/include/basic_types/datatype_structure.hpp +++ b/Common/include/basic_types/datatype_structure.hpp @@ -95,19 +95,9 @@ namespace SU2_TYPE { FORCEINLINE void SetDerivative(su2double& data, const passivedouble &val) {data.setGradient(val);} -#ifdef CODI_REVERSE_TYPE - FORCEINLINE passivedouble GetSecondary(const su2double& data) { - return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]); - } - - FORCEINLINE passivedouble GetDerivative(const su2double& data) { - return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]); - } -#else // forward FORCEINLINE passivedouble GetSecondary(const su2double& data) {return data.getGradient();} FORCEINLINE passivedouble GetDerivative(const su2double& data) {return data.getGradient();} -#endif #else // passive type, no AD diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index c0d1ac313ba..91f033b5164 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1234,14 +1234,14 @@ class CGeometry { * \brief Ray Intersects Triangle (Moller and Trumbore algorithm) */ bool RayIntersectsTriangle(const su2double orig[3], const su2double dir[3], - const su2double vert0[3], const su2double vert1[3], const su2double vert2[3], - su2double *intersect); + const su2double vert0[3], const su2double vert1[3], const su2double vert2[3], + su2double *intersect); /*! * \brief Segment Intersects Triangle */ bool SegmentIntersectsTriangle(su2double point0[3], const su2double point1[3], - su2double vert0[3], su2double vert1[3], su2double vert2[3]); + su2double vert0[3], su2double vert1[3], su2double vert2[3]); /*! * \brief Segment Intersects Line (for 2D FFD Intersection) diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 56a21fc1ad9..91ff7135eb4 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -818,7 +818,7 @@ class CPoint { inline void RegisterCoordinates(unsigned long iPoint, bool input) { for (unsigned long iDim = 0; iDim < nDim; iDim++) { if(input) { - AD::RegisterInput(Coord(iPoint,iDim),false); + AD::RegisterInput(Coord(iPoint,iDim)); AD::SetIndex(AD_InputIndex(iPoint,iDim), Coord(iPoint,iDim)); } else { diff --git a/Common/src/basic_types/ad_structure.cpp b/Common/src/basic_types/ad_structure.cpp index 18342e13a90..f6defb62435 100644 --- a/Common/src/basic_types/ad_structure.cpp +++ b/Common/src/basic_types/ad_structure.cpp @@ -31,12 +31,6 @@ namespace AD { #ifdef CODI_REVERSE_TYPE /*--- Initialization of the global variables ---*/ - int adjointVectorPosition = 0; - - std::vector inputValues; - std::vector localInputValues; - std::vector localOutputValues; - TapePosition StartPosition, EndPosition; std::vector TapePositions; diff --git a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp index e81ac660eb2..750bf1dda8f 100644 --- a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp @@ -48,7 +48,6 @@ class CDiscAdjFEASolver final : public CSolver { struct SensData { unsigned short size = 0; su2double* val = nullptr; /*!< \brief Value of the variable. */ - int* AD_Idx = nullptr; /*!< \brief Derivative index in the AD tape. */ su2double* LocalSens = nullptr; /*!< \brief Local sensitivity (domain). */ su2double* GlobalSens = nullptr; /*!< \brief Global sensitivity (mpi). */ su2double* TotalSens = nullptr; /*!< \brief Total sensitivity (time domain). */ @@ -60,7 +59,6 @@ class CDiscAdjFEASolver final : public CSolver { clear(); size = n; val = new su2double[n](); - AD_Idx = new int[n](); LocalSens = new su2double[n](); GlobalSens = new su2double[n](); TotalSens = new su2double[n](); @@ -69,21 +67,17 @@ class CDiscAdjFEASolver final : public CSolver { void clear() { size = 0; delete [] val; - delete [] AD_Idx; delete [] LocalSens; delete [] GlobalSens; delete [] TotalSens; } void Register() { - for (auto i = 0u; i < size; ++i) { - AD::RegisterInput(val[i], false); - AD::SetIndex(AD_Idx[i], val[i]); - } + for (auto i = 0u; i < size; ++i) AD::RegisterInput(val[i]); } void GetDerivative() { - for (auto i = 0u; i < size; ++i) LocalSens[i] = AD::GetDerivative(AD_Idx[i]); + for (auto i = 0u; i < size; ++i) LocalSens[i] = SU2_TYPE::GetDerivative(val[i]); SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); } diff --git a/SU2_CFD/include/variables/CMeshVariable.hpp b/SU2_CFD/include/variables/CMeshVariable.hpp index c420e523206..217c51ce5f3 100644 --- a/SU2_CFD/include/variables/CMeshVariable.hpp +++ b/SU2_CFD/include/variables/CMeshVariable.hpp @@ -95,7 +95,7 @@ class CMeshVariable : public CVariable { */ inline void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) final { for (unsigned long iDim = 0; iDim < nDim; iDim++) { - adj_mesh[iDim] = AD::GetDerivative(AD_InputIndex(iPoint,iDim)); + adj_mesh[iDim] = SU2_TYPE::GetDerivative(Mesh_Coord(iPoint,iDim)); AD::ResetInput(Mesh_Coord(iPoint,iDim)); } } diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 8f049516587..1898006dbac 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -95,8 +95,6 @@ class CVariable { su2matrix AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ - su2matrix AD_Time_n_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ - su2matrix AD_Time_n1_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ unsigned long nPoint = 0; /*!< \brief Number of points in the domain. */ unsigned long nDim = 0; /*!< \brief Number of dimension of the problem. */ @@ -114,21 +112,25 @@ class CVariable { assert(false && "A base method of CVariable was used, but it should have been overridden by the derived class."); } - void RegisterContainer(bool input, su2activematrix& variable, su2matrix& ad_index) { + void RegisterContainer(bool input, su2activematrix& variable, su2matrix* ad_index = nullptr) { const auto nPoint = variable.rows(); SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads())) for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { for(unsigned long iVar=0; iVar& ad_index) { + RegisterContainer(input, variable, &ad_index); + } + public: /*--- Disable copy and assignment. ---*/ CVariable(const CVariable&) = delete; @@ -2159,13 +2161,13 @@ class CVariable { } inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const { - for (unsigned long iVar = 0; iVar < AD_Time_n_InputIndex.cols(); iVar++) - adj_sol[iVar] = AD::GetDerivative(AD_Time_n_InputIndex(iPoint,iVar)); + for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++) + adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n(iPoint,iVar)); } inline void GetAdjointSolution_time_n1(unsigned long iPoint, su2double *adj_sol) const { - for (unsigned long iVar = 0; iVar < nVar; iVar++) - adj_sol[iVar] = AD::GetDerivative(AD_Time_n1_InputIndex(iPoint,iVar)); + for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); iVar++) + adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n1(iPoint,iVar)); } /*! diff --git a/SU2_CFD/src/variables/CMeshVariable.cpp b/SU2_CFD/src/variables/CMeshVariable.cpp index aed8430db5e..143c0f1188f 100644 --- a/SU2_CFD/src/variables/CMeshVariable.cpp +++ b/SU2_CFD/src/variables/CMeshVariable.cpp @@ -38,9 +38,6 @@ CMeshVariable::CMeshVariable(unsigned long npoint, unsigned long ndim, CConfig * Mesh_Coord.resize(nPoint,nDim) = su2double(0.0); WallDistance.resize(nPoint) = su2double(1e-9); - if (config->GetDiscrete_Adjoint()) - AD_InputIndex.resize(nPoint,nDim) = -1; - /*--- Initialize the variables necessary when the problem is time domain ---*/ if (config->GetTime_Domain()) { Solution_time_n.resize(nPoint,nDim) = su2double(0.0); @@ -49,5 +46,5 @@ CMeshVariable::CMeshVariable(unsigned long npoint, unsigned long ndim, CConfig * } void CMeshVariable::Register_MeshCoord() { - RegisterContainer(true, Mesh_Coord, AD_InputIndex); + RegisterContainer(true, Mesh_Coord); } diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index ce5f9a03c4c..78e6ae7b59d 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -72,12 +72,6 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva if (!adjoint) { AD_InputIndex.resize(nPoint,nVar) = -1; AD_OutputIndex.resize(nPoint,nVar) = -1; - - if (config->GetTime_Domain()) - AD_Time_n_InputIndex.resize(nPoint,nVar) = -1; - - if (config->GetTime_Marching() != TIME_MARCHING::STEADY) - AD_Time_n1_InputIndex.resize(nPoint,nVar) = -1; } } @@ -122,9 +116,9 @@ void CVariable::RegisterSolution(bool input) { } void CVariable::RegisterSolution_time_n() { - RegisterContainer(true, Solution_time_n, AD_Time_n_InputIndex); + RegisterContainer(true, Solution_time_n); } void CVariable::RegisterSolution_time_n1() { - RegisterContainer(true, Solution_time_n1, AD_Time_n1_InputIndex); + RegisterContainer(true, Solution_time_n1); } diff --git a/TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg b/TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg index e5fb7fcd096..04cf4c83c6e 100755 --- a/TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg +++ b/TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg @@ -43,7 +43,7 @@ MARKER_DESIGNING= ( leading_edge, pressure_side, suction_side ) % Common numerics settings --------------------------------------------- % REF_DIMENSIONALIZATION= DIMENSIONAL NUM_METHOD_GRAD= GREEN_GAUSS -CFL_NUMBER= 15.0 +CFL_NUMBER= 1000.0 % % Flow numerics -------------------------------------------------------- % CONV_NUM_METHOD_FLOW= JST @@ -53,18 +53,15 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT % Linear solvers ------------------------------------------------------- % LINEAR_SOLVER= BCGSTAB LINEAR_SOLVER_PREC= ILU -LINEAR_SOLVER_ERROR= 1E-3 -LINEAR_SOLVER_ITER= 1000 -DISCADJ_LIN_SOLVER= BCGSTAB +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 25 +DISCADJ_LIN_SOLVER= SMOOTHER DISCADJ_LIN_PREC= ILU +LINEAR_SOLVER_SMOOTHER_RELAXATION= 0.7 % Multigrid -MGLEVEL= 2 -MGCYCLE= V_CYCLE -MG_PRE_SMOOTH= ( 1, 2, 3, 3 ) -MG_POST_SMOOTH= ( 0, 0, 0, 0 ) -MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 ) -MG_DAMP_RESTRICTION= 0.75 -MG_DAMP_PROLONGATION= 0.75 +MGLEVEL= 0 +NEWTON_KRYLOV= YES +QUASI_NEWTON_NUM_SAMPLES= 999 % DEFORM_LINEAR_SOLVER= CONJUGATE_GRADIENT DEFORM_LINEAR_SOLVER_PREC= ILU @@ -78,9 +75,9 @@ DEFORM_POISSONS_RATIO= 1e6 BGS_RELAXATION= FIXED_PARAMETER STAT_RELAX_PARAMETER= 1.0 % fluid -INNER_ITER= 51 +INNER_ITER= 31 CONV_STARTITER= 0 -CONV_RESIDUAL_MINVAL= -9 +CONV_RESIDUAL_MINVAL= -7 % % In\Out --------------------------------------------------------------- % MESH_FILENAME= mesh.su2 From 567064a5e2364fa992edb9a5e92aa89b1d78566c Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 25 May 2021 18:43:09 +0100 Subject: [PATCH 16/21] fixes --- Common/include/CConfig.hpp | 3 +- SU2_DOT/src/SU2_DOT.cpp | 37 ++++++------------- .../disc_adj_fsi/dyn_fsi/grad_dv.opt.ref | 16 ++++---- UnitTests/Common/simple_ad_test.cpp | 2 +- 4 files changed, 22 insertions(+), 36 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index f403284c37c..f1952f35355 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -5444,7 +5444,8 @@ class CConfig { * \param[in] val_val - Value of the design variable that we want to read. * \return Design variable step. */ - su2double GetDV_Value(unsigned short val_dv, unsigned short val_val = 0) const { return DV_Value[val_dv][val_val]; } + su2double& GetDV_Value(unsigned short val_dv, unsigned short val_val = 0) { return DV_Value[val_dv][val_val]; } + const su2double& GetDV_Value(unsigned short val_dv, unsigned short val_val = 0) const { return DV_Value[val_dv][val_val]; } /*! * \brief Set the value of the design variable step, we use this value in design problems. diff --git a/SU2_DOT/src/SU2_DOT.cpp b/SU2_DOT/src/SU2_DOT.cpp index f544b6e63fe..35aab262fff 100644 --- a/SU2_DOT/src/SU2_DOT.cpp +++ b/SU2_DOT/src/SU2_DOT.cpp @@ -688,18 +688,16 @@ void SetProjection_FD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *surface_movement, su2double** Gradient){ - su2double DV_Value, *VarCoord, Sensitivity, my_Gradient, localGradient, *Normal, Area = 0.0; - unsigned short iDV_Value = 0, iMarker, nMarker, iDim, nDim, iDV, nDV, nDV_Value; + su2double *VarCoord = nullptr, Sensitivity, my_Gradient, localGradient, *Normal, Area = 0.0; + unsigned short iDV_Value = 0, iMarker, nMarker, iDim, nDim, iDV, nDV; unsigned long iVertex, nVertex, iPoint; - int rank = SU2_MPI::GetRank(); + const int rank = SU2_MPI::GetRank(); nMarker = config->GetnMarker_All(); nDim = geometry->GetnDim(); nDV = config->GetnDV(); - VarCoord = nullptr; - /*--- Discrete adjoint gradient computation ---*/ if (rank == MASTER_NODE) @@ -712,21 +710,13 @@ void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su /*--- Register design variables as input and set them to zero * (since we want to have the derivative at alpha = 0, i.e. for the current design) ---*/ - - for (iDV = 0; iDV < nDV; iDV++){ - nDV_Value = config->GetnDV_Value(iDV); - - for (iDV_Value = 0; iDV_Value < nDV_Value; iDV_Value++){ + for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++){ - /*--- Initilization with su2double resets the index ---*/ + config->SetDV_Value(iDV, iDV_Value, 0.0); - DV_Value = 0.0; - - AD::RegisterInput(DV_Value); - - config->SetDV_Value(iDV, iDV_Value, DV_Value); + AD::RegisterInput(config->GetDV_Value(iDV, iDV_Value)); } } @@ -742,10 +732,7 @@ void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su * We need that to make sure to set the sensitivity of surface points only once * (Markers share points, so we would visit them more than once in the loop over the markers below) ---*/ - bool* visited = new bool[geometry->GetnPoint()]; - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++){ - visited[iPoint] = false; - } + vector visited(geometry->GetnPoint(), false); /*--- Initialize the derivatives of the output of the surface deformation routine * with the discrete adjoints from the CFD solution ---*/ @@ -779,18 +766,16 @@ void SetProjection_AD(CGeometry *geometry, CConfig *config, CSurfaceMovement *su } } - delete [] visited; - /*--- Compute derivatives and extract gradient ---*/ AD::ComputeAdjoint(); for (iDV = 0; iDV < nDV; iDV++){ - nDV_Value = config->GetnDV_Value(iDV); - for (iDV_Value = 0; iDV_Value < nDV_Value; iDV_Value++){ - DV_Value = config->GetDV_Value(iDV, iDV_Value); - my_Gradient = SU2_TYPE::GetDerivative(DV_Value); + for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++){ + + my_Gradient = SU2_TYPE::GetDerivative(config->GetDV_Value(iDV, iDV_Value)); + SU2_MPI::Allreduce(&my_Gradient, &localGradient, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); /*--- Angle of Attack design variable (this is different, diff --git a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref index 72cc7d70416..3bcc2d6eedf 100644 --- a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref +++ b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref @@ -1,9 +1,9 @@ INDEX GRAD -0 -3.465632121754074e-03 -1 -1.844003360219930e-03 -2 -7.925070240684198e-04 -3 -2.742895873196515e-04 -4 -2.738093165742900e-04 -5 -7.890548338696854e-04 -6 -1.831172969426724e-03 -7 -3.431365309796604e-03 +0 -3.461460672772851e-03 +1 -1.841786315630031e-03 +2 -7.915536278471477e-04 +3 -2.739622075935052e-04 +4 -2.734869089862929e-04 +5 -7.881162336341228e-04 +6 -1.828978274796512e-03 +7 -3.427219375030697e-03 diff --git a/UnitTests/Common/simple_ad_test.cpp b/UnitTests/Common/simple_ad_test.cpp index 6e95b218bb4..49d18467f43 100644 --- a/UnitTests/Common/simple_ad_test.cpp +++ b/UnitTests/Common/simple_ad_test.cpp @@ -56,5 +56,5 @@ TEST_CASE("Simple AD Test", "[AD tests]") { AD::ComputeAdjoint(); CHECK(SU2_TYPE::GetValue(y) == Approx(64)); - CHECK(SU2_TYPE::GetDerivative(y) == Approx(48)); + CHECK(SU2_TYPE::GetDerivative(x) == Approx(48)); } From ef4f23329ed9b07a91920055afe37a7d9819d094 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 26 May 2021 16:09:42 +0100 Subject: [PATCH 17/21] proper-er way to reset input --- Common/include/basic_types/ad_structure.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index 6ecbaaac531..18c430c7b2f 100644 --- a/Common/include/basic_types/ad_structure.hpp +++ b/Common/include/basic_types/ad_structure.hpp @@ -299,7 +299,7 @@ namespace AD{ FORCEINLINE void RegisterOutput(su2double& data) {AD::getGlobalTape().registerOutput(data);} - FORCEINLINE void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();} + FORCEINLINE void ResetInput(su2double &data) {data = data.getValue();} FORCEINLINE void StartRecording() {AD::getGlobalTape().setActive();} From 9df5d023626eb7c6e8089c22e2bc0ae6f8904333 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Wed, 26 May 2021 17:29:10 +0100 Subject: [PATCH 18/21] small things --- Common/src/geometry/CGeometry.cpp | 4 +--- SU2_CFD/src/drivers/CDriver.cpp | 18 ++++++++++++------ SU2_CFD/src/solvers/CSolver.cpp | 2 +- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index a7b1fe7b9fa..62e0c2b3551 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -2721,8 +2721,6 @@ void CGeometry::ComputeSurf_Curvature(CConfig *config) { su2double U[3] = {0.0}, V[3] = {0.0}, W[3] = {0.0}, Length_U, Length_V, Length_W, CosValue, Angle_Value, *K, *Angle_Defect, *Area_Vertex, *Angle_Alpha, *Angle_Beta, **NormalMeanK, MeanK, GaussK, MaxPrinK, cot_alpha, cot_beta, delta, X1, X2, X3, Y1, Y2, Y3, radius, *Buffer_Send_Coord, *Buffer_Receive_Coord, *Coord, Dist, MinDist, MaxK, MinK, SigmaK; bool *Check_Edge; - const bool fea = config->GetStructuralProblem(); - /*--- Allocate surface curvature ---*/ K = new su2double [nPoint]; for (iPoint = 0; iPoint < nPoint; iPoint++) K[iPoint] = 0.0; @@ -2986,7 +2984,7 @@ void CGeometry::ComputeSurf_Curvature(CConfig *config) { SigmaK = sqrt(SigmaK/su2double(TotalnPointDomain)); - if ((rank == MASTER_NODE) && (!fea)) + if (rank == MASTER_NODE) cout << "Max K: " << MaxK << ". Mean K: " << MeanK << ". Standard deviation K: " << SigmaK << "." << endl; Point_Critical.clear(); diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index ee98c0b25a3..677f19a7bb8 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -773,7 +773,7 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet /*--- Create the control volume structures ---*/ - if ((rank == MASTER_NODE) && (!fea)) cout << "Setting the control volume structure." << endl; + if (rank == MASTER_NODE) cout << "Setting the control volume structure." << endl; SU2_OMP_PARALLEL { geometry[MESH_0]->SetControlVolume(config, ALLOCATE); geometry[MESH_0]->SetBoundControlVolume(config, ALLOCATE); @@ -798,8 +798,10 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet /*--- Compute the surface curvature ---*/ - if ((rank == MASTER_NODE) && (!fea)) cout << "Compute the surface curvature." << endl; - geometry[MESH_0]->ComputeSurf_Curvature(config); + if (!fea) { + if (rank == MASTER_NODE) cout << "Compute the surface curvature." << endl; + geometry[MESH_0]->ComputeSurf_Curvature(config); + } /*--- Check for periodicity and disable MG if necessary. ---*/ @@ -895,13 +897,17 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet /*--- Compute the max length. ---*/ - if ((rank == MASTER_NODE) && (!fea) && (iMGlevel == MESH_0)) cout << "Finding max control volume width." << endl; - geometry[iMGlevel]->SetMaxLength(config); + if (!fea) { + if ((rank == MASTER_NODE) && (iMGlevel == MESH_0)) + cout << "Finding max control volume width." << endl; + geometry[iMGlevel]->SetMaxLength(config); + } /*--- Communicate the number of neighbors. This is needed for some centered schemes and for multigrid in parallel. ---*/ - if ((rank == MASTER_NODE) && (size > SINGLE_NODE) && (!fea) && (iMGlevel == MESH_0)) cout << "Communicating number of neighbors." << endl; + if ((rank == MASTER_NODE) && (size > SINGLE_NODE) && (iMGlevel == MESH_0)) + cout << "Communicating number of neighbors." << endl; geometry[iMGlevel]->InitiateComms(geometry[iMGlevel], config, NEIGHBORS); geometry[iMGlevel]->CompleteComms(geometry[iMGlevel], config, NEIGHBORS); } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index c36e6404f54..9f8da0fbd13 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -3129,7 +3129,7 @@ void CSolver::InterpolateRestartData(const CGeometry *geometry, const CConfig *c if (size != SINGLE_NODE && size % 2) SU2_MPI::Error("Number of ranks must be multiple of 2.", CURRENT_FUNCTION); - if (config->GetStructuralProblem() || config->GetFEMSolver()) + if (config->GetFEMSolver()) SU2_MPI::Error("Cannot interpolate the restart file for FEM problems.", CURRENT_FUNCTION); /* Challenges: From 19d3329da5b596c19a463e0ef8fbeae72e1f02e0 Mon Sep 17 00:00:00 2001 From: TobiKattmann Date: Thu, 27 May 2021 12:05:32 +0200 Subject: [PATCH 19/21] Remove errors from #1281. No auto for su2double with rhs math operation and MAXNVAR fix for flow. --- SU2_CFD/src/solvers/CHeatSolver.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index c845fa0e8f1..f8a8b862033 100644 --- a/SU2_CFD/src/solvers/CHeatSolver.cpp +++ b/SU2_CFD/src/solvers/CHeatSolver.cpp @@ -363,8 +363,8 @@ void CHeatSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe nVarFlow = solver_container[FLOW_SOL]->GetnVar(); /*--- Define some auxiliary vectors related to the primitive flow solution ---*/ - su2double Primitive_Flow_i[MAXNVAR] = {0}; - su2double Primitive_Flow_j[MAXNVAR] = {0}; + vector Primitive_Flow_i(nVarFlow, 0.0); + vector Primitive_Flow_j(nVarFlow, 0.0); for (auto iEdge = 0ul; iEdge < geometry->GetnEdge(); iEdge++) { @@ -421,7 +421,7 @@ void CHeatSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe Temp_i_Corrected = Temp_i + Project_Temp_i_Grad; Temp_j_Corrected = Temp_j + Project_Temp_j_Grad; - numerics->SetPrimitive(Primitive_Flow_i, Primitive_Flow_j); + numerics->SetPrimitive(Primitive_Flow_i.data(), Primitive_Flow_j.data()); numerics->SetTemperature(Temp_i_Corrected, Temp_j_Corrected); } @@ -571,9 +571,9 @@ void CHeatSolver::BC_Isothermal_Wall(CGeometry *geometry, CSolver **solver_conta su2double *Normal, *Coord_i, *Coord_j, Area, dist_ij, Twall, dTdn; const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - const auto laminar_viscosity = config->GetMu_ConstantND(); - const auto Prandtl_Lam = config->GetPrandtl_Lam(); - const auto thermal_diffusivity = flow ? laminar_viscosity/Prandtl_Lam : config->GetThermalDiffusivity_Solid(); + const su2double laminar_viscosity = config->GetMu_ConstantND(); + const su2double Prandtl_Lam = config->GetPrandtl_Lam(); + const su2double thermal_diffusivity = flow ? laminar_viscosity/Prandtl_Lam : config->GetThermalDiffusivity_Solid(); //su2double Prandtl_Turb = config->GetPrandtl_Turb(); //laminar_viscosity = config->GetViscosity_FreeStreamND(); // TDE check for consistency for CHT @@ -832,8 +832,8 @@ void CHeatSolver::BC_ConjugateHeat_Interface(CGeometry *geometry, CSolver **solv const bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); - const auto Temperature_Ref = config->GetTemperature_Ref(); - const auto rho_cp_solid = config->GetDensity_Solid()*config->GetSpecific_Heat_Cp(); + const su2double Temperature_Ref = config->GetTemperature_Ref(); + const su2double rho_cp_solid = config->GetDensity_Solid()*config->GetSpecific_Heat_Cp(); if (flow) { From e8cee0bd684dd45d25996bde7edd3fb0d38188b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Bl=C3=BChdorn?= Date: Thu, 27 May 2021 14:20:48 +0200 Subject: [PATCH 20/21] CoDiPack update. --- externals/codi | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/externals/codi b/externals/codi index 6a67202a388..c78918f5eab 160000 --- a/externals/codi +++ b/externals/codi @@ -1 +1 @@ -Subproject commit 6a67202a3887c8da490fdfde82bc46507de68692 +Subproject commit c78918f5eabfc3ef06c8f00a82fb9880ff2c4c75 From be27e1d5ab45ec240fd0cba610d503e2a8a2929a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Bl=C3=BChdorn?= Date: Thu, 27 May 2021 14:21:11 +0200 Subject: [PATCH 21/21] Testing RealReverseIndex. --- Common/include/code_config.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp index 17315154c34..a341541c920 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -92,7 +92,7 @@ using su2double = codi::RealReversePrimal; #elif CODI_PRIMAL_INDEX_TAPE using su2double = codi::RealReversePrimalIndex; #else -using su2double = codi::RealReverse; +using su2double = codi::RealReverseIndex; #endif #endif #elif defined(CODI_FORWARD_TYPE) // forward mode AD