diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index c8a9819435f3..f1952f353559 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. @@ -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/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp index 620da3246f53..18c430c7b2f7 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,28 +289,17 @@ 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);} - FORCEINLINE void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();} + FORCEINLINE void ResetInput(su2double &data) {data = data.getValue();} FORCEINLINE void StartRecording() {AD::getGlobalTape().setActive();} @@ -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 039df3312006..b30d81555f12 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/code_config.hpp b/Common/include/code_config.hpp index 1524b720affd..c1600b310f42 100644 --- a/Common/include/code_config.hpp +++ b/Common/include/code_config.hpp @@ -123,9 +123,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 diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 61dd15714a18..bdde61d4290f 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) @@ -1250,16 +1250,15 @@ 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. * \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 417b88974fff..487b0836292e 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -563,7 +563,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(); @@ -831,52 +831,23 @@ class CPoint { } /*! - * \brief Set the adjoint values of the coordinates. + * \brief Register coordinates of a point. * \param[in] iPoint - Index of the point. - * \param[in] adj_sol - The adjoint values of the coordinates. + * \param[in] input - Register as input or output. */ - 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 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)); + inline void RegisterCoordinates(unsigned long iPoint, bool input) { + for (unsigned long iDim = 0; iDim < nDim; iDim++) { + if(input) { + AD::RegisterInput(Coord(iPoint,iDim)); + 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 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); - /*! * \brief Set wall roughnesses according to stored closest wall information. * \param[in] roughness - Mapping [rank][zone][marker] -> roughness @@ -892,4 +863,5 @@ class CPoint { } } } + }; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index abe3f2359927..8a51da7a7286 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/Common/src/basic_types/ad_structure.cpp b/Common/src/basic_types/ad_structure.cpp index 18342e13a905..f6defb624350 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/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index 3f1c90667921..07634cd18501 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -2476,18 +2476,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 } @@ -2496,10 +2490,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); @@ -2732,8 +2722,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; @@ -2997,7 +2985,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/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index 6fd84b1afde8..19ec490c29e3 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; } @@ -153,17 +153,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 613b5a7fb35e..750bf1dda8f8 100644 --- a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp @@ -48,8 +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. */ - 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). */ @@ -61,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,28 +66,18 @@ class CDiscAdjFEASolver final : public CSolver { void clear() { size = 0; - localIdx = false; delete [] val; - delete [] AD_Idx; delete [] LocalSens; delete [] GlobalSens; 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]); } 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] = SU2_TYPE::GetDerivative(val[i]); SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); } diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index fa533937244b..5102b521a13c 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,22 @@ 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); - } + for (iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) { - /*--- Recompute the edges and dual mesh control volumes in the - domain and on the boundaries. ---*/ - - geometry[MESH_0]->SetControlVolume(config, UPDATE); - geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); - geometry[MESH_0]->SetMaxLength(config); - - /*--- Update the multigrid structure after setting up the finest grid, - including computing the grid velocities on the coarser levels. ---*/ - - 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); + /*--- 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]->SetMaxLength(config); } } @@ -973,7 +965,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); } @@ -996,26 +988,36 @@ 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 && (TimeIter == config->GetRestart_Iter()) && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { - /*--- Load an additional restart file for a 2nd-order restart ---*/ + if (restart && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) { + + /*--- 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); + 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 5a4a6d927f03..27fd06dd39c7 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 946c4e36502b..ccf0bbf74089 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/include/variables/CFEAVariable.hpp b/SU2_CFD/include/variables/CFEAVariable.hpp index 5dc269577178..975f4b24f562 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/CMeshVariable.hpp b/SU2_CFD/include/variables/CMeshVariable.hpp index 7b917f51ffc0..217c51ce5f33 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++) + inline void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) final { + for (unsigned long iDim = 0; iDim < nDim; 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 db37de808482..1898006dbac1 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,6 +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 = 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; @@ -1926,7 +1943,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. @@ -2024,7 +2041,7 @@ class CVariable { /*! * \brief A virtual member. Recover the value of the adjoint of the mesh coordinates. */ - inline virtual void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) const { } + inline virtual void GetAdjoint_MeshCoord(unsigned long iPoint, su2double *adj_mesh) { } /*! * \brief A virtual member. Get the value of the displacement imposed at the boundary. @@ -2112,84 +2129,47 @@ 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++) + for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); 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)); - } - /*! * \brief Set the sensitivity at the node * \param[in] iDim - spacial component diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index fc2b099129d3..c52d5afe0e05 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -949,7 +949,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); } } diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index c789bf70f381..677f19a7bb8d 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 @@ -235,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; @@ -250,9 +239,7 @@ 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; - + 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], @@ -262,7 +249,9 @@ CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunica } 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); } @@ -330,40 +319,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"); @@ -436,8 +413,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]; } } @@ -446,20 +422,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]; } @@ -480,10 +453,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; @@ -497,11 +469,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; @@ -513,18 +482,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; @@ -658,7 +622,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++) { @@ -756,7 +720,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. ---*/ @@ -809,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); @@ -834,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. ---*/ @@ -892,6 +858,7 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet if (config->GetnMGLevels() != requestedMGlevels) { delete geometry[iMGlevel]; + geometry[iMGlevel] = nullptr; break; } @@ -902,16 +869,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(); + } } } @@ -928,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); } @@ -964,7 +937,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); @@ -1040,7 +1013,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); @@ -1050,16 +1023,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 ---*/ @@ -1137,170 +1107,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 = config->GetRestart_Iter()-1; - else val_iter = 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 ((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 (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 (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 = config->GetRestart_Iter()-1; - else val_iter = 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, @@ -1359,7 +1215,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; @@ -1546,9 +1402,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; } @@ -1573,9 +1427,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 : @@ -1593,9 +1445,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++) @@ -1712,9 +1562,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; } @@ -1729,18 +1577,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; } @@ -1839,9 +1683,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 : @@ -1850,9 +1692,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; } @@ -1904,9 +1744,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; } @@ -1914,9 +1752,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; } @@ -1986,9 +1822,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; } @@ -2002,9 +1836,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++) { @@ -2015,9 +1847,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; } @@ -2067,9 +1897,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++) { @@ -2077,9 +1905,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; } @@ -2122,9 +1948,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; } } @@ -2148,17 +1972,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 : @@ -2171,9 +1991,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; } @@ -2200,18 +2018,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; } @@ -2273,25 +2087,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; } @@ -2454,7 +2262,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 50e3bbd96e47..314f21daace5 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -239,12 +239,12 @@ 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); } - 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 fd7322330242..2b9f488d4931 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], @@ -164,14 +164,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/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index 6fe780efba46..eb1a1890a4f2 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -173,10 +173,10 @@ void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geom solver[iZone][iInst][MESH_0][FEA_SOL]->RegisterVariables(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(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 6e17f208f397..87449a46cc4d 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]); } @@ -456,11 +456,10 @@ 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()) { - /*--- 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]); @@ -504,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); } @@ -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/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp index 7eee95121dea..be710ab4b2ed 100644 --- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp @@ -188,7 +188,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 ---*/ @@ -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 077aa5b988f9..5f053a80275f 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 b1b23e74c42d..229d61f866c5 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 c3de5fbf3a5f..291127f4e08b 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/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 485549393321..3f0e634230ee 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,82 +229,48 @@ 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); } -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++){ - - /*--- Extract the adjoint solution ---*/ - - if(multizone) { - direct_solver->GetNodes()->GetAdjointSolution_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); - } + if (!config->GetMultizone_Problem()) nodes->Set_OldSolution(); - /*--- Store 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); nodes->SetSolution(iPoint,Solution); - } - /*--- Solution for acceleration (u'') and velocity (u') at time n ---*/ - - if (dynamic){ - - /*--- NOW: The solution at time n ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++){ - - /*--- Extract the adjoint solution at time n ---*/ + if (CrossTerm) return; - if(multizone) { - direct_solver->GetNodes()->GetAdjointSolution_time_n_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - } + /*--- Extract and store the adjoint solution at time n (including accel. and velocity) ---*/ - /*--- Store the adjoint solution at time n ---*/ - - if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); + if (config->GetTime_Domain()) { + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + su2double Solution[MAXNVAR] = {0.0}; + direct_solver->GetNodes()->GetAdjointSolution_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)); @@ -387,10 +339,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 +367,18 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo /*--- Extract the geometric sensitivities ---*/ + const bool time_domain = config->GetTime_Domain(); + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - su2double *Coord = geometry->nodes->GetCoord(iPoint); + auto Coord = geometry->nodes->GetCoord(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); + AD::ResetInput(Coord[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/CDiscAdjMeshSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp index 0c81be0c3793..98a93919dc0f 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 ce0dc2fff6f8..84ce1b3906c0 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,25 +292,19 @@ 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){ +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}; @@ -319,19 +312,15 @@ 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 ---*/ - if(config->GetMultizone_Problem()) { - direct_solver->GetNodes()->GetAdjointSolution_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); - } + su2double Solution[MAXNVAR] = {0.0}; + direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution); /*--- Relax and store the adjoint solution, compute the residuals. ---*/ @@ -352,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++) { @@ -373,15 +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++) { - - if(config->GetMultizone_Problem()) { - direct_solver->GetNodes()->GetAdjointSolution_time_n_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); - } - - if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); + su2double Solution[MAXNVAR] = {0.0}; + direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); + nodes->Set_Solution_time_n(iPoint,Solution); } END_SU2_OMP_FOR } @@ -390,15 +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++) { - - if(config->GetMultizone_Problem()) { - direct_solver->GetNodes()->GetAdjointSolution_time_n1_LocalIndex(iPoint,Solution); - } - else { - direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); - } - - if (!CrossTerm) nodes->Set_Solution_time_n1(iPoint,Solution); + su2double Solution[MAXNVAR] = {0.0}; + direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); + nodes->Set_Solution_time_n1(iPoint,Solution); } END_SU2_OMP_FOR } @@ -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 } @@ -518,17 +498,7 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve 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. ---*/ - + 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/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 30307741607f..7763d01f1b14 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) { diff --git a/SU2_CFD/src/solvers/CHeatSolver.cpp b/SU2_CFD/src/solvers/CHeatSolver.cpp index c845fa0e8f11..f8a8b8620333 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) { diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 06af9b8408c6..e8f68c9a3afb 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 @@ -519,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 @@ -535,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 ---*/ @@ -568,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) @@ -593,10 +584,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 @@ -606,27 +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++) { - su2double val_vel = LinSysSol(iPoint, iDim); - - /*--- Non-dimensionalize velocity ---*/ - val_vel /= config->GetVelocity_Ref(); + 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 - geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, val_vel); - } + 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. ---*/ @@ -634,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); } @@ -830,8 +805,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. ---*/ @@ -851,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) { @@ -862,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){ + /*--- 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. ---*/ @@ -905,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 ---*/ @@ -928,8 +888,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) { diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index c36e6404f54b..9f8da0fbd13b 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: diff --git a/SU2_CFD/src/variables/CMeshVariable.cpp b/SU2_CFD/src/variables/CMeshVariable.cpp index 95aa08ab1aad..143c0f1188fc 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; @@ -42,14 +39,12 @@ CMeshVariable::CMeshVariable(unsigned long npoint, unsigned long ndim, CConfig * WallDistance.resize(nPoint) = su2double(1e-9); /*--- 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); } diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 5882fd2189b3..78e6ae7b59d6 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -65,19 +65,13 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva if (config->GetTime_Marching() != TIME_MARCHING::STEADY) Solution_time_n1.resize(nPoint,nVar) = su2double(0.0); - if (config->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; - - 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; } } @@ -117,31 +111,14 @@ 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) { - 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; iVarGetnMarker_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/Airfoil_2d/config.cfg b/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg index 012ad702857a..818f359cdac9 100755 --- a/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg +++ b/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg @@ -7,7 +7,7 @@ MARKER_ZONE_INTERFACE= (pressure_side,pressure_side_s, suction_side,suction_side CONSERVATIVE_INTERPOLATION= NO OUTER_ITER= 9 -OUTPUT_WRT_FREQ= 99999 +OUTPUT_WRT_FREQ= 5 MESH_FILENAME= mesh.su2 diff --git a/TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg b/TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg index e5fb7fcd096f..04cf4c83c6e9 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 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 9f9eb61496da..3bcc2d6eedf1 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.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/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index fb06b8d5a6a2..af6011b0eaf1 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.743230, -4.133462] 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.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 b0c92bcb386f..54728bb237a0 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.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.398550, -4.086739, 0.000000, 117.000000] + 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.768496, -4.159961] #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 98fd0feb4995..6a715ac54781 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.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.398530, -4.086741, 0.000000, 101.000000] #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.768501, -4.159959] #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 diff --git a/UnitTests/Common/simple_ad_test.cpp b/UnitTests/Common/simple_ad_test.cpp index 6e95b218bb44..49d18467f43c 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)); }