diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index dc96c21ad3b..b15c00e333b 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -488,7 +488,7 @@ class CConfig { unsigned short **DegreeFFDBox; /*!< \brief Degree of the FFD boxes. */ string *FFDTag; /*!< \brief Parameters of the design variable. */ string *TagFFDBox; /*!< \brief Tag of the FFD box. */ - unsigned short GeometryMode; /*!< \brief Gemoetry mode (analysis or gradient computation). */ + unsigned short GeometryMode; /*!< \brief Geometry mode (analysis or gradient computation). */ unsigned short MGCycle; /*!< \brief Kind of multigrid cycle. */ unsigned short FinestMesh; /*!< \brief Finest mesh for the full multigrid approach. */ unsigned short nFFD_Fix_IDir, @@ -2912,7 +2912,7 @@ class CConfig { unsigned short GetFinestMesh(void) const { return FinestMesh; } /*! - * \brief Get the kind of multigrid (V or W). + * \brief Get the kind of multigrid (V, W or FULLMG). * \note This variable is used in a recursive way to perform the different kind of cycles * \return 0 or 1 depending of we are dealing with a V or W cycle. */ diff --git a/Common/include/geometry/CMultiGridGeometry.hpp b/Common/include/geometry/CMultiGridGeometry.hpp index 00faa8126c6..84fe5723e15 100644 --- a/Common/include/geometry/CMultiGridGeometry.hpp +++ b/Common/include/geometry/CMultiGridGeometry.hpp @@ -28,28 +28,29 @@ #pragma once #include "CGeometry.hpp" +class CMultiGridQueue; /*! * \class CMultiGridGeometry - * \brief Class for defining the multigrid geometry, the main delicated part is the + * \brief Class for defining the multigrid geometry, the main delegated part is the * agglomeration stage, which is done in the declaration. * \author F. Palacios */ class CMultiGridGeometry final : public CGeometry { private: /*! - * \brief Determine if a CVPoint van be agglomerated, if it have the same marker point as the seed. + * \brief Determine if a CVPoint can be agglomerated, if it has the same marker point as the seed. * \param[in] CVPoint - Control volume to be agglomerated. * \param[in] marker_seed - Marker of the seed. * \param[in] fine_grid - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. * \return TRUE or FALSE depending if the control volume can be agglomerated. */ - bool SetBoundAgglomeration(unsigned long CVPoint, short marker_seed, const CGeometry* fine_grid, + bool SetBoundAgglomeration(unsigned long CVPoint, vector marker_seed, const CGeometry* fine_grid, const CConfig* config) const; /*! - * \brief Determine if a can be agglomerated using geometrical criteria. + * \brief Determine if a Point can be agglomerated using geometrical criteria. * \param[in] iPoint - Seed point. * \param[in] fine_grid - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. @@ -57,7 +58,7 @@ class CMultiGridGeometry final : public CGeometry { bool GeometricalCheck(unsigned long iPoint, const CGeometry* fine_grid, const CConfig* config) const; /*! - * \brief Determine if a CVPoint van be agglomerated, if it have the same marker point as the seed. + * \brief Determine if a CVPoint can be agglomerated, if it has the same marker point as the seed. * \param[out] Suitable_Indirect_Neighbors - List of Indirect Neighbours that can be agglomerated. * \param[in] iPoint - Seed point. * \param[in] Index_CoarseCV - Index of agglomerated point. @@ -66,6 +67,15 @@ class CMultiGridGeometry final : public CGeometry { void SetSuitableNeighbors(vector& Suitable_Indirect_Neighbors, unsigned long iPoint, unsigned long Index_CoarseCV, const CGeometry* fine_grid) const; + /*! + * \brief Compute local curvature at a boundary vertex on Euler wall. + * \param[in] fine_grid - Fine grid geometry. + * \param[in] iPoint - Point index. + * \param[in] iMarker - Marker index. + * \return Maximum angle (in degrees) between this vertex normal and adjacent vertex normals. + */ + su2double ComputeLocalCurvature(const CGeometry* fine_grid, unsigned long iPoint, unsigned short iMarker) const; + public: /*--- This is to suppress Woverloaded-virtual, omitting it has no negative impact. ---*/ using CGeometry::SetBoundControlVolume; diff --git a/Common/include/geometry/CMultiGridQueue.hpp b/Common/include/geometry/CMultiGridQueue.hpp index 296d4e8f73d..c16e23520e4 100644 --- a/Common/include/geometry/CMultiGridQueue.hpp +++ b/Common/include/geometry/CMultiGridQueue.hpp @@ -93,7 +93,7 @@ class CMultiGridQueue { void IncrPriorityCV(unsigned long incrPoint); /*! - * \brief Increase the priority of the CV. + * \brief Reduce the priority of the CV. * \param[in] redPoint - Index of the control volume. */ void RedPriorityCV(unsigned long redPoint); diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index a52d111237c..9937e7faef4 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -29,20 +29,36 @@ #include "../../include/geometry/CMultiGridQueue.hpp" #include "../../include/toolboxes/printing_toolbox.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include +#include +#include +#include +#include +#include +#include CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, unsigned short iMesh) : CGeometry() { nDim = fine_grid->GetnDim(); // Write the number of dimensions of the coarse grid. - /*--- Create a queue system to do the agglomeration + /*--- Maximum agglomeration size in 2D is 4 nodes, in 3D is 8 nodes. ---*/ + const short int maxAgglomSize = (nDim == 2) ? 4 : 8; + + /*--- Inherit boundary properties from fine grid ---*/ + boundIsStraight = fine_grid->boundIsStraight; + + /*--- Agglomeration Scheme II (Nishikawa, Diskin, Thomas) + Create a queue system to do the agglomeration 1st) More than two markers ---> Vertices (never agglomerate) 2nd) Two markers ---> Edges (agglomerate if same BC, never agglomerate if different BC) 3rd) One marker ---> Surface (always agglomerate) 4th) No marker ---> Internal Volume (always agglomerate) ---*/ + // Note that for MPI, we introduce interfaces and we can choose to have agglomeration over + // the interface or not. Nishikawa chooses not to agglomerate over interfaces. + /*--- Set a marker to indicate indirect agglomeration, for quads and hexs, - i.e. consider up to neighbors of neighbors of neighbors. + i.e. consider up to neighbors of neighbors. For other levels this information is propagated down during their construction. ---*/ - if (iMesh == MESH_1) { for (auto iPoint = 0ul; iPoint < fine_grid->GetnPoint(); iPoint++) fine_grid->nodes->SetAgglomerate_Indirect(iPoint, false); @@ -59,7 +75,6 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } /*--- Create the coarse grid structure using as baseline the fine grid ---*/ - CMultiGridQueue MGQueue_InnerCV(fine_grid->GetnPoint()); vector Suitable_Indirect_Neighbors; @@ -67,15 +82,26 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un unsigned long Index_CoarseCV = 0; - /*--- The first step is the boundary agglomeration. ---*/ + /*--- Statistics for Euler wall agglomeration ---*/ + map euler_wall_agglomerated, euler_wall_rejected_curvature, + euler_wall_rejected_straight; + for (unsigned short iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) { + euler_wall_agglomerated[iMarker] = 0; + euler_wall_rejected_curvature[iMarker] = 0; + euler_wall_rejected_straight[iMarker] = 0; + } + } + /*--- STEP 1: The first step is the boundary agglomeration. ---*/ for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { for (auto iVertex = 0ul; iVertex < fine_grid->GetnVertex(iMarker); iVertex++) { const auto iPoint = fine_grid->vertex[iMarker][iVertex]->GetNode(); /*--- If the element has not been previously agglomerated and it - belongs to this physical domain, and it meets the geometrical - criteria, the agglomeration is studied. ---*/ + belongs to this physical domain, and it meets the geometrical + criteria, the agglomeration is studied. ---*/ + vector marker_seed; if ((!fine_grid->nodes->GetAgglomerate(iPoint)) && (fine_grid->nodes->GetDomain(iPoint)) && (GeometricalCheck(iPoint, fine_grid, config))) { @@ -89,94 +115,146 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- We add the seed point (child) to the parent control volume ---*/ nodes->SetChildren_CV(Index_CoarseCV, 0, iPoint); - bool agglomerate_seed = true; + bool agglomerate_seed = false; auto counter = 0; unsigned short copy_marker[3] = {}; - const auto marker_seed = iMarker; + marker_seed.push_back(iMarker); /*--- For a particular point in the fine grid we save all the markers that are in that point ---*/ - for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker() && counter < 3; jMarker++) { + for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker(); jMarker++) { + const string Marker_Tag = config->GetMarker_All_TagBound(iMarker); if (fine_grid->nodes->GetVertex(iPoint, jMarker) != -1) { copy_marker[counter] = jMarker; counter++; + + if (jMarker != iMarker) { + marker_seed.push_back(jMarker); + } } } - /*--- To aglomerate a vertex it must have only one physical bc!! - This can be improved. If there is only a marker, it is a good + /*--- To agglomerate a vertex it must have only one physical bc. + This can be improved. If there is only one marker, it is a good candidate for agglomeration ---*/ + /*--- 1 BC, so either an edge in 2D or the interior of a plane in 3D ---*/ + /*--- Valley -> Valley : conditionally allowed when both points are on the same marker. ---*/ + /*--- ! Note that in the case of MPI SEND_RECEIVE markers, we might need other conditions ---*/ if (counter == 1) { + // The seed/parent is one valley, so we set this part to true + // if the child is on the same valley, we set it to true as well. agglomerate_seed = true; - - /*--- Euler walls can be curved and agglomerating them leads to difficulties ---*/ - if (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL) agglomerate_seed = false; + /*--- Euler walls: check curvature-based agglomeration criterion ---*/ + if (config->GetMarker_All_KindBC(marker_seed[0]) == EULER_WALL) { + /*--- Allow agglomeration if marker is straight OR local curvature is small ---*/ + if (!boundIsStraight[marker_seed[0]]) { + /*--- Compute local curvature at this point ---*/ + su2double local_curvature = ComputeLocalCurvature(fine_grid, iPoint, marker_seed[0]); + // limit to 45 degrees + if (local_curvature >= 45.0) { + agglomerate_seed = false; // High curvature: do not agglomerate + euler_wall_rejected_curvature[marker_seed[0]]++; + } else { + euler_wall_agglomerated[marker_seed[0]]++; + } + } else { + /*--- Straight wall: agglomerate ---*/ + euler_wall_agglomerated[marker_seed[0]]++; + } + } + /*--- Note that if the (single) marker is a SEND_RECEIVE, then the node is actually an interior point. + In that case it can only be agglomerated with another interior point. ---*/ + if (config->GetMarker_All_KindBC(marker_seed[0]) == SEND_RECEIVE) { + agglomerate_seed = true; + } } - /*--- If there are two markers, we will agglomerate if any of the - markers is SEND_RECEIVE ---*/ + /*--- Note that in 2D, this is a corner and we do not agglomerate unless they both are SEND_RECEIVE. ---*/ + /*--- In 3D, we agglomerate if the 2 markers are the same. ---*/ if (counter == 2) { - agglomerate_seed = (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) || - (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE); - - /* --- Euler walls can also not be agglomerated when the point has 2 markers ---*/ - if ((config->GetMarker_All_KindBC(copy_marker[0]) == EULER_WALL) || - (config->GetMarker_All_KindBC(copy_marker[1]) == EULER_WALL)) { - agglomerate_seed = false; + /*--- agglomerate if one of the 2 markers are MPI markers. ---*/ + if (nDim == 2) { + // agglomerate_seed = ((config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) || + // (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE)); + } + /*--- agglomerate if both markers are the same. ---*/ + if (nDim == 3) agglomerate_seed = (copy_marker[0] == copy_marker[1]); + + /*--- Euler walls: check curvature-based agglomeration criterion for both markers ---*/ + // only in 3d because in 2d it's a corner + bool euler_wall_rejected_here = false; + for (unsigned short i = 0; i < 2; i++) { + if ((nDim == 3) && (config->GetMarker_All_KindBC(copy_marker[i]) == EULER_WALL)) { + if (!boundIsStraight[copy_marker[i]]) { + /*--- Compute local curvature at this point ---*/ + su2double local_curvature = ComputeLocalCurvature(fine_grid, iPoint, copy_marker[i]); + // limit to 45 degrees + if (local_curvature >= 45.0) { + agglomerate_seed = false; // High curvature: do not agglomerate + euler_wall_rejected_curvature[copy_marker[i]]++; + euler_wall_rejected_here = true; + } + } + /*--- Track agglomeration if not rejected ---*/ + if (agglomerate_seed && !euler_wall_rejected_here) { + euler_wall_agglomerated[copy_marker[i]]++; + } + } } } /*--- If there are more than 2 markers, the aglomeration will be discarded ---*/ - if (counter > 2) agglomerate_seed = false; - /*--- If the seed can be agglomerated, we try to agglomerate more points ---*/ - + /*--- If the seed (parent) can be agglomerated, we try to agglomerate connected childs to the parent ---*/ + /*--- Note that in 2D we allow a maximum of 4 nodes to be agglomerated ---*/ if (agglomerate_seed) { /*--- Now we do a sweep over all the nodes that surround the seed point ---*/ - for (auto CVPoint : fine_grid->nodes->GetPoints(iPoint)) { /*--- The new point can be agglomerated ---*/ - if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) { /*--- We set the value of the parent ---*/ - fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); /*--- We set the value of the child ---*/ - nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); nChildren++; + /*--- In 2D, we agglomerate exactly 2 nodes if the nodes are on the line edge. ---*/ + if ((nDim == 2) && (counter == 1)) break; + /*--- In 3D, we agglomerate exactly 2 nodes if the nodes are on the surface edge. ---*/ + if ((nDim == 3) && (counter == 2)) break; + /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes (counter==1 in 3D). ---*/ + if (nChildren == maxAgglomSize) break; } } - Suitable_Indirect_Neighbors.clear(); - - if (fine_grid->nodes->GetAgglomerate_Indirect(iPoint)) - SetSuitableNeighbors(Suitable_Indirect_Neighbors, iPoint, Index_CoarseCV, fine_grid); - - /*--- Now we do a sweep over all the indirect nodes that can be added ---*/ - - for (auto CVPoint : Suitable_Indirect_Neighbors) { - /*--- The new point can be agglomerated ---*/ - - if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) { - /*--- We set the value of the parent ---*/ - - fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); - - /*--- We set the indirect agglomeration information of the corse point - based on its children in the fine grid. ---*/ - - if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint)) - nodes->SetAgglomerate_Indirect(Index_CoarseCV, true); - - /*--- We set the value of the child ---*/ - - nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); - nChildren++; + /*--- Only take into account indirect neighbors for 3D faces, not 2D. ---*/ + if (nDim == 3) { + Suitable_Indirect_Neighbors.clear(); + + if (fine_grid->nodes->GetAgglomerate_Indirect(iPoint)) + SetSuitableNeighbors(Suitable_Indirect_Neighbors, iPoint, Index_CoarseCV, fine_grid); + + /*--- Now we do a sweep over all the indirect nodes that can be added ---*/ + for (auto CVPoint : Suitable_Indirect_Neighbors) { + /*--- The new point can be agglomerated ---*/ + if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) { + /*--- We set the value of the parent ---*/ + fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); + + /*--- We set the indirect agglomeration information of the corse point + based on its children in the fine grid. ---*/ + if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint)) + nodes->SetAgglomerate_Indirect(Index_CoarseCV, true); + + /*--- We set the value of the child ---*/ + nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint); + nChildren++; + /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes. ---*/ + if (nChildren == maxAgglomSize) break; + } } } } @@ -223,8 +301,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } - /*--- Agglomerate the domain points. ---*/ - + /*--- STEP 2: Agglomerate the domain points. ---*/ auto iteration = 0ul; while (!MGQueue_InnerCV.EmptyQueue() && (iteration < fine_grid->GetnPoint())) { const auto iPoint = MGQueue_InnerCV.NextCV(); @@ -271,6 +348,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un MGQueue_InnerCV.Update(CVPoint, fine_grid); } + if (nChildren == maxAgglomSize) break; } /*--- Identify the indirect neighbors ---*/ @@ -282,11 +360,11 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un /*--- Now we do a sweep over all the indirect nodes that can be added ---*/ for (auto CVPoint : Suitable_Indirect_Neighbors) { + // if we have reached the maximum, get out. + if (nChildren == maxAgglomSize) break; /*--- The new point can be agglomerated ---*/ - if ((!fine_grid->nodes->GetAgglomerate(CVPoint)) && (fine_grid->nodes->GetDomain(CVPoint))) { /*--- We set the value of the parent ---*/ - fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV); /*--- We set the indirect agglomeration information ---*/ @@ -340,16 +418,56 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un for (auto iCoarsePoint = 0ul; iCoarsePoint < nPointDomain; iCoarsePoint++) { if (nodes->GetnPoint(iCoarsePoint) == 1) { /*--- Find the neighbor of the isolated point. This neighbor is the right control volume ---*/ - const auto iCoarsePoint_Complete = nodes->GetPoint(iCoarsePoint, 0); - /*--- Add the children to the connected control volume (and modify its parent indexing). - Identify the child CV from the finest grid and add it to the correct control volume. - Set the parent CV of iFinePoint. Instead of using the original one - (iCoarsePoint), use the new one (iCoarsePoint_Complete) ---*/ + /*--- Check if merging would exceed the maximum agglomeration size ---*/ + auto nChildren_Target = nodes->GetnChildren_CV(iCoarsePoint_Complete); + auto nChildren_Isolated = nodes->GetnChildren_CV(iCoarsePoint); + auto nChildren_Total = nChildren_Target + nChildren_Isolated; + + /*--- If the total would exceed maxAgglomSize, try to redistribute children to neighbors ---*/ + if (nChildren_Total > maxAgglomSize) { + /*--- Find neighbors of the target coarse point that have room ---*/ + unsigned short nChildrenToRedistribute = nChildren_Total - maxAgglomSize; + + for (auto jCoarsePoint : nodes->GetPoints(iCoarsePoint_Complete)) { + if (nChildrenToRedistribute == 0) break; + + auto nChildren_Neighbor = nodes->GetnChildren_CV(jCoarsePoint); + if (nChildren_Neighbor < maxAgglomSize) { + unsigned short nCanTransfer = + min(nChildrenToRedistribute, static_cast(maxAgglomSize - nChildren_Neighbor)); - auto nChildren = nodes->GetnChildren_CV(iCoarsePoint_Complete); + /*--- Transfer children from target to neighbor ---*/ + for (unsigned short iTransfer = 0; iTransfer < nCanTransfer; iTransfer++) { + /*--- Take from the end of the target's children list ---*/ + auto nChildren_Current = nodes->GetnChildren_CV(iCoarsePoint_Complete); + if (nChildren_Current > 0) { + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint_Complete, nChildren_Current - 1); + /*--- Add to neighbor ---*/ + auto nChildren_Neighbor_Current = nodes->GetnChildren_CV(jCoarsePoint); + nodes->SetChildren_CV(jCoarsePoint, nChildren_Neighbor_Current, iFinePoint); + nodes->SetnChildren_CV(jCoarsePoint, nChildren_Neighbor_Current + 1); + + /*--- Update parent ---*/ + fine_grid->nodes->SetParent_CV(iFinePoint, jCoarsePoint); + + /*--- Remove from target (by reducing count) ---*/ + nodes->SetnChildren_CV(iCoarsePoint_Complete, nChildren_Current - 1); + + nChildrenToRedistribute--; + } + } + } + } + + /*--- Update the target's child count after redistribution ---*/ + nChildren_Target = nodes->GetnChildren_CV(iCoarsePoint_Complete); + } + + /*--- Add the isolated point's children to the target control volume ---*/ + auto nChildren = nChildren_Target; for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { const auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); nodes->SetChildren_CV(iCoarsePoint_Complete, nChildren, iFinePoint); @@ -358,7 +476,6 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } /*--- Update the number of children control volumes ---*/ - nodes->SetnChildren_CV(iCoarsePoint_Complete, nChildren); nodes->SetnChildren_CV(iCoarsePoint, 0); } @@ -369,6 +486,16 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un nodes->ResetPoints(); #ifdef HAVE_MPI + /*--- Reset halo point parents before MPI agglomeration. + This is critical for multi-level multigrid: when creating level N from level N-1, + the fine grid (level N-1) already has Parent_CV set from when it was created from level N-2. + Those parent indices point to level N, but when creating level N+1, they would be + incorrectly interpreted as level N+1 indices. Resetting ensures clean agglomeration. ---*/ + + for (auto iPoint = fine_grid->GetnPointDomain(); iPoint < fine_grid->GetnPoint(); iPoint++) { + fine_grid->nodes->SetParent_CV(iPoint, ULONG_MAX); + } + /*--- Dealing with MPI parallelization, the objective is that the received nodes must be agglomerated in the same way as the donor (send) nodes. Send the node agglomeration information of the donor (parent and children). The agglomerated halos of this rank are set according to the rank where @@ -376,8 +503,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && (config->GetMarker_All_SendRecv(iMarker) > 0)) { - const auto MarkerS = iMarker; - const auto MarkerR = iMarker + 1; + const auto MarkerS = iMarker; // sending marker + const auto MarkerR = iMarker + 1; // receiving marker const auto send_to = config->GetMarker_All_SendRecv(MarkerS) - 1; const auto receive_from = abs(config->GetMarker_All_SendRecv(MarkerR)) - 1; @@ -424,38 +551,151 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un vector Parent_Local(nVertexR); vector Children_Local(nVertexR); + /*--- First pass: Determine which parents will actually be used (have non-skipped children). + This prevents creating orphaned halo CVs that have coordinates (0,0,0). ---*/ + vector parent_used(Aux_Parent.size(), false); + vector parent_local_index(Aux_Parent.size(), ULONG_MAX); + for (auto iVertex = 0ul; iVertex < nVertexR; iVertex++) { - /*--- We use the same sorting as in the donor domain, i.e. the local parents - are numbered according to their order in the remote rank. ---*/ + const auto iPoint_Fine = fine_grid->vertex[MarkerR][iVertex]->GetNode(); + auto existing_parent = fine_grid->nodes->GetParent_CV(iPoint_Fine); + + /*--- Skip if already agglomerated (first-wins policy) ---*/ + if (existing_parent != ULONG_MAX) continue; + + /*--- Find which parent this vertex maps to ---*/ + for (auto jVertex = 0ul; jVertex < Aux_Parent.size(); jVertex++) { + if (Parent_Remote[iVertex] == Aux_Parent[jVertex]) { + parent_used[jVertex] = true; + break; + } + } + } + /*--- Assign local indices only to used parents ---*/ + unsigned long nUsedParents = 0; + for (auto jVertex = 0ul; jVertex < Aux_Parent.size(); jVertex++) { + if (parent_used[jVertex]) { + parent_local_index[jVertex] = Index_CoarseCV + nUsedParents; + nUsedParents++; + } + } + + /*--- Now map each received vertex to its local parent ---*/ + for (auto iVertex = 0ul; iVertex < nVertexR; iVertex++) { + Parent_Local[iVertex] = ULONG_MAX; for (auto jVertex = 0ul; jVertex < Aux_Parent.size(); jVertex++) { if (Parent_Remote[iVertex] == Aux_Parent[jVertex]) { - Parent_Local[iVertex] = jVertex + Index_CoarseCV; + Parent_Local[iVertex] = parent_local_index[jVertex]; break; } } + + /*--- Validate that parent mapping was found (only matters if not skipped later) ---*/ + if (Parent_Local[iVertex] == ULONG_MAX) { + SU2_MPI::Error(string("MPI agglomeration failed to map parent index ") + to_string(Parent_Remote[iVertex]) + + string(" for vertex ") + to_string(iVertex), + CURRENT_FUNCTION); + } + Children_Local[iVertex] = fine_grid->vertex[MarkerR][iVertex]->GetNode(); } - Index_CoarseCV += Aux_Parent.size(); + /*--- Only increment by the number of parents that will actually be used ---*/ + Index_CoarseCV += nUsedParents; - vector nChildren_MPI(Index_CoarseCV, 0); + /*--- Debug counters ---*/ + unsigned long nConflicts = 0, nSkipped = 0, nOutOfBounds = 0, nSuccess = 0; /*--- Create the final structure ---*/ for (auto iVertex = 0ul; iVertex < nVertexR; iVertex++) { const auto iPoint_Coarse = Parent_Local[iVertex]; const auto iPoint_Fine = Children_Local[iVertex]; - /*--- Be careful, it is possible that a node changes the agglomeration configuration, - the priority is always when receiving the information. ---*/ + /*--- Debug: Check for out-of-bounds access ---*/ + if (iPoint_Coarse >= Index_CoarseCV) { + cout << "ERROR [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR + << "]: Out-of-bounds coarse CV index " << iPoint_Coarse << " >= " << Index_CoarseCV << " (vertex " + << iVertex << ", fine point " << iPoint_Fine << ")" << endl; + nOutOfBounds++; + continue; + } + + /*--- Solution 1: Skip if this halo point was already agglomerated ---*/ + auto existing_parent = fine_grid->nodes->GetParent_CV(iPoint_Fine); + if (existing_parent != ULONG_MAX) { + if (existing_parent != iPoint_Coarse) { + /*--- Conflict detected: different parent from different interface ---*/ + nConflicts++; + + /*--- Only print detailed info for first few conflicts or if suspicious ---*/ + if (nConflicts <= 5 || existing_parent < nPointDomain) { + cout << "INFO [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR << "]: Halo point " + << iPoint_Fine << " already agglomerated to parent " << existing_parent + << (existing_parent < nPointDomain ? " (DOMAIN CV!)" : " (halo CV)") << ", skipping reassignment to " + << iPoint_Coarse << " (from rank " << receive_from << ")" << endl; + } + } else { + /*--- Same parent from different interface (duplicate) - just skip silently ---*/ + nSkipped++; + } + continue; // First-wins: keep existing assignment + } + + /*--- First assignment for this halo point - proceed with agglomeration ---*/ + + /*--- Critical fix: Append to existing children, don't overwrite ---*/ + auto existing_children_count = nodes->GetnChildren_CV(iPoint_Coarse); + fine_grid->nodes->SetParent_CV(iPoint_Fine, iPoint_Coarse); - nodes->SetChildren_CV(iPoint_Coarse, nChildren_MPI[iPoint_Coarse], iPoint_Fine); - nChildren_MPI[iPoint_Coarse]++; - nodes->SetnChildren_CV(iPoint_Coarse, nChildren_MPI[iPoint_Coarse]); + nodes->SetChildren_CV(iPoint_Coarse, existing_children_count, iPoint_Fine); + nodes->SetnChildren_CV(iPoint_Coarse, existing_children_count + 1); nodes->SetDomain(iPoint_Coarse, false); + nSuccess++; + } + + /*--- Debug: Report statistics for this marker pair ---*/ + cout << "MPI Agglomeration [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR << " (rank " << send_to + << " <-> " << receive_from << ")]: " << nSuccess << " assigned, " << nSkipped << " duplicates, " + << nConflicts << " conflicts"; + if (nOutOfBounds > 0) { + cout << ", " << nOutOfBounds << " OUT-OF-BOUNDS (CRITICAL!)"; + } + cout << endl; + + if (nConflicts > 5) { + cout << " Note: Only first 5 conflicts shown in detail, total conflicts = " << nConflicts << endl; + } + + /*--- Debug: Validate buffer size assumption ---*/ + if (nVertexS != nVertexR) { + cout << "WARNING [Rank " << rank << ", Marker " << MarkerS << "/" << MarkerR + << "]: Asymmetric interface - nVertexS=" << nVertexS << " != nVertexR=" << nVertexR << endl; } } } + + /*--- Post-process: Check for orphaned coarse CVs (should be none with new logic) ---*/ + + if (size > SINGLE_NODE) { + unsigned long nOrphaned = 0; + + /*--- Count orphaned CVs for reporting ---*/ + for (auto iCoarse = nPointDomain; iCoarse < Index_CoarseCV; iCoarse++) { + if (nodes->GetnChildren_CV(iCoarse) == 0) { + nOrphaned++; + /*--- This shouldn't happen with the new parent prefiltering logic ---*/ + cout << "WARNING [Rank " << rank << "]: Orphaned halo CV " << iCoarse + << " detected (should not occur with current logic)" << endl; + } + } + + if (nOrphaned > 0) { + cout << "WARNING [Rank " << rank << "]: " << nOrphaned + << " orphaned halo coarse CVs found - this indicates a logic error!" << endl; + } + } + #endif // HAVE_MPI /*--- Update the number of points after the MPI agglomeration ---*/ @@ -463,25 +703,28 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un nPoint = Index_CoarseCV; /*--- Console output with the summary of the agglomeration ---*/ - - unsigned long nPointFine = fine_grid->GetnPoint(); + // nijso: do not include halo points in the count + unsigned long nPointFine = fine_grid->GetnPointDomain(); unsigned long Global_nPointCoarse, Global_nPointFine; - SU2_MPI::Allreduce(&nPoint, &Global_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&nPointDomain, &Global_nPointCoarse, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); SU2_MPI::Allreduce(&nPointFine, &Global_nPointFine, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); SetGlobal_nPointDomain(Global_nPointCoarse); if (iMesh != MESH_0) { - const su2double factor = 1.5; - const su2double Coeff = pow(su2double(Global_nPointFine) / Global_nPointCoarse, 1.0 / nDim); + /*--- Note: CFL at the coarse levels have a large impact on convergence, + this should be rewritten to use adaptive CFL. ---*/ + const su2double factor = 1.0; + const su2double Coeff = 1.5; const su2double CFL = factor * config->GetCFL(iMesh - 1) / Coeff; config->SetCFL(iMesh, CFL); } const su2double ratio = su2double(Global_nPointFine) / su2double(Global_nPointCoarse); - - if (((nDim == 2) && (ratio < 2.5)) || ((nDim == 3) && (ratio < 2.5))) { + cout << "********** ratio = " << ratio << endl; + // lower value leads to more levels being accepted. + if (((nDim == 2) && (ratio < 1.5)) || ((nDim == 3) && (ratio < 1.5))) { config->SetMGLevels(iMesh - 1); } else if (rank == MASTER_NODE) { PrintingToolbox::CTablePrinter MGTable(&std::cout); @@ -503,11 +746,68 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un } } + /*--- Output Euler wall agglomeration statistics ---*/ + if (rank == MASTER_NODE) { + /*--- Gather global statistics for Euler walls ---*/ + bool has_euler_walls = false; + for (unsigned short iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) { + has_euler_walls = true; + break; + } + } + + if (has_euler_walls) { + cout << endl; + cout << "Euler Wall Agglomeration Statistics (45° curvature threshold):" << endl; + cout << "----------------------------------------------------------------" << endl; + + for (unsigned short iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) { + string marker_name = config->GetMarker_All_TagBound(iMarker); + unsigned long agglomerated = euler_wall_agglomerated[iMarker]; + unsigned long rejected = euler_wall_rejected_curvature[iMarker]; + unsigned long total = agglomerated + rejected; + + if (total > 0) { + su2double accept_rate = 100.0 * su2double(agglomerated) / su2double(total); + cout << " Marker: " << marker_name << endl; + cout << " Seeds agglomerated: " << agglomerated << " (" << std::setprecision(1) << std::fixed + << accept_rate << "%)" << endl; + cout << " Seeds rejected (>45° curv): " << rejected << " (" << std::setprecision(1) << std::fixed + << (100.0 - accept_rate) << "%)" << endl; + } + } + } + cout << "----------------------------------------------------------------" << endl; + } + } + edgeColorGroupSize = config->GetEdgeColoringGroupSize(); } -bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, short marker_seed, const CGeometry* fine_grid, - const CConfig* config) const { +bool CMultiGridGeometry::GeometricalCheck(unsigned long iPoint, const CGeometry* fine_grid, + const CConfig* config) const { + su2double max_dimension = 1.2; + + /*--- Evaluate the total size of the element ---*/ + + bool Volume = true; + su2double ratio = pow(fine_grid->nodes->GetVolume(iPoint), 1.0 / su2double(nDim)) * max_dimension; + su2double limit = pow(config->GetDomainVolume(), 1.0 / su2double(nDim)); + if (ratio > limit) { + Volume = false; + cout << "Volume limit reached!" << endl; + } + /*--- Evaluate the stretching of the element ---*/ + + bool Stretching = true; + + return (Stretching && Volume); +} + +bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, vector marker_seed, + const CGeometry* fine_grid, const CConfig* config) const { bool agglomerate_CV = false; /*--- Basic condition, the point has not been previously agglomerated, it belongs to the domain, @@ -517,12 +817,13 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, short mark (GeometricalCheck(CVPoint, fine_grid, config))) { /*--- If the point belongs to a boundary, its type must be compatible with the seed marker. ---*/ + int counter = 0; + unsigned short copy_marker[3] = {}; + if (fine_grid->nodes->GetBoundary(CVPoint)) { /*--- Identify the markers of the vertex that we want to agglomerate ---*/ // count number of markers on the agglomeration candidate - int counter = 0; - unsigned short copy_marker[3] = {}; for (auto jMarker = 0u; jMarker < fine_grid->GetnMarker() && counter < 3; jMarker++) { if (fine_grid->nodes->GetVertex(CVPoint, jMarker) != -1) { copy_marker[counter] = jMarker; @@ -530,94 +831,61 @@ bool CMultiGridGeometry::SetBoundAgglomeration(unsigned long CVPoint, short mark } } - /*--- The basic condition is that the aglomerated vertex must have the same physical marker, + /*--- The basic condition is that the agglomerated vertex must have the same physical marker, but eventually a send-receive condition ---*/ - /*--- Only one marker in the vertex that is going to be aglomerated ---*/ + /*--- Only one marker in the vertex that is going to be agglomerated ---*/ + /*--- Valley -> Valley: only if of the same type---*/ if (counter == 1) { /*--- We agglomerate if there is only one marker and it is the same marker as the seed marker ---*/ - // note that this should be the same marker id, not just the same marker type - if (copy_marker[0] == marker_seed) agglomerate_CV = true; - - /*--- If there is only one marker, but the marker is the SEND_RECEIVE ---*/ - - if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) { - agglomerate_CV = true; - } - - if ((config->GetMarker_All_KindBC(marker_seed) == SYMMETRY_PLANE) || - (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL)) { - if (config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) { - agglomerate_CV = false; - } - } + // So this is the case when in 2D we are on an edge, and in 3D we are in the interior of a surface. + // note that this should be the same marker id, not just the same marker type. + if ((marker_seed.size() == 1) && (copy_marker[0] == marker_seed[0])) agglomerate_CV = true; + + // /*--- If seed has 2 markers and one is SEND_RECEIVE, do not agglomerate with physical boundary ---*/ + // if ((marker_seed.size() == 2) && (config->GetMarker_All_KindBC(marker_seed[0]) == SEND_RECEIVE)) { + // if (copy_marker[0] == marker_seed[1]) { + // agglomerate_CV = false; + // } + // } + // if ((marker_seed.size() == 2) && (config->GetMarker_All_KindBC(marker_seed[1]) == SEND_RECEIVE)) { + // if (copy_marker[0] == marker_seed[0]) { + // agglomerate_CV = false; + // } + // } + + /*--- Note: If there is only one marker, but the marker is the SEND_RECEIVE, then the point is actually an + interior point and we do not agglomerate. ---*/ } /*--- If there are two markers in the vertex that is going to be aglomerated ---*/ if (counter == 2) { - /*--- First we verify that the seed is a physical boundary ---*/ - - if (config->GetMarker_All_KindBC(marker_seed) != SEND_RECEIVE) { - /*--- Then we check that one of the markers is equal to the seed marker, and the other is send/receive ---*/ + /*--- In 2D this is a corner and we do not agglomerate ---*/ + if (nDim == 2) { + agglomerate_CV = false; + } + /*--- Both markers have to be the same. ---*/ - if (((copy_marker[0] == marker_seed) && (config->GetMarker_All_KindBC(copy_marker[1]) == SEND_RECEIVE)) || - ((config->GetMarker_All_KindBC(copy_marker[0]) == SEND_RECEIVE) && (copy_marker[1] == marker_seed))) { + if (marker_seed.size() == 2) { + if (((copy_marker[0] == marker_seed[0]) && (copy_marker[1] == marker_seed[1])) || + ((copy_marker[0] == marker_seed[1]) && (copy_marker[1] == marker_seed[0]))) { agglomerate_CV = true; } } } } - /*--- If the element belongs to the domain, it is always agglomerated. ---*/ + /*--- If the element belongs to the domain, it is never agglomerated with a boundary node. ---*/ else { - agglomerate_CV = true; - - // actually, for symmetry (and possibly other cells) we only agglomerate cells that are on the marker - // at this point, the seed was on the boundary and the CV was not. so we check if the seed is a symmetry - if ((config->GetMarker_All_KindBC(marker_seed) == SYMMETRY_PLANE) || - (config->GetMarker_All_KindBC(marker_seed) == EULER_WALL)) { - agglomerate_CV = false; - } + agglomerate_CV = false; } } return agglomerate_CV; } -bool CMultiGridGeometry::GeometricalCheck(unsigned long iPoint, const CGeometry* fine_grid, - const CConfig* config) const { - su2double max_dimension = 1.2; - - /*--- Evaluate the total size of the element ---*/ - - bool Volume = true; - su2double ratio = pow(fine_grid->nodes->GetVolume(iPoint), 1.0 / su2double(nDim)) * max_dimension; - su2double limit = pow(config->GetDomainVolume(), 1.0 / su2double(nDim)); - if (ratio > limit) Volume = false; - - /*--- Evaluate the stretching of the element ---*/ - - bool Stretching = true; - - /* unsigned short iNode, iDim; - unsigned long jPoint; - su2double *Coord_i = fine_grid->nodes->GetCoord(iPoint); - su2double max_dist = 0.0 ; su2double min_dist = 1E20; - for (iNode = 0; iNode < fine_grid->nodes->GetnPoint(iPoint); iNode ++) { - jPoint = fine_grid->nodes->GetPoint(iPoint, iNode); - su2double *Coord_j = fine_grid->nodes->GetCoord(jPoint); - su2double distance = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - distance += (Coord_j[iDim]-Coord_i[iDim])*(Coord_j[iDim]-Coord_i[iDim]); - distance = sqrt(distance); - max_dist = max(distance, max_dist); - min_dist = min(distance, min_dist); - } - if ( max_dist/min_dist > 100.0 ) Stretching = false;*/ - - return (Stretching && Volume); -} +/*--- ---*/ void CMultiGridGeometry::SetSuitableNeighbors(vector& Suitable_Indirect_Neighbors, unsigned long iPoint, unsigned long Index_CoarseCV, const CGeometry* fine_grid) const { @@ -637,8 +905,8 @@ void CMultiGridGeometry::SetSuitableNeighbors(vector& Suitable_In auto end = First_Neighbor_Points.end(); if (find(First_Neighbor_Points.begin(), end, kPoint) == end) { - Second_Neighbor_Points.push_back(kPoint); - Second_Origin_Points.push_back(jPoint); + Second_Neighbor_Points.push_back(kPoint); // neighbor of a neighbor, not connected to original ipoint + Second_Origin_Points.push_back(jPoint); // the neighbor that is connected to ipoint } } } @@ -668,53 +936,11 @@ void CMultiGridGeometry::SetSuitableNeighbors(vector& Suitable_In sort(Suitable_Second_Neighbors.begin(), Suitable_Second_Neighbors.end()); auto it1 = unique(Suitable_Second_Neighbors.begin(), Suitable_Second_Neighbors.end()); Suitable_Second_Neighbors.resize(it1 - Suitable_Second_Neighbors.begin()); - - /*--- Create a list with the third neighbors, without first, second, and seed neighbors. ---*/ - - /// TODO: This repeats the process above but I doubt it catches any more points. - - vector Third_Neighbor_Points, Third_Origin_Points; - - for (auto kPoint : Suitable_Second_Neighbors) { - for (auto lPoint : fine_grid->nodes->GetPoints(kPoint)) { - /*--- Check that the third neighbor does not belong to the first neighbors or the seed ---*/ - - auto end1 = First_Neighbor_Points.end(); - if (find(First_Neighbor_Points.begin(), end1, lPoint) != end1) continue; - - /*--- Check that the third neighbor does not belong to the second neighbors ---*/ - - auto end2 = Suitable_Second_Neighbors.end(); - if (find(Suitable_Second_Neighbors.begin(), end2, lPoint) != end2) continue; - - Third_Neighbor_Points.push_back(lPoint); - Third_Origin_Points.push_back(kPoint); - } - } - - /*--- Identify those third neighbors that are repeated (candidate to be added). ---*/ - - for (auto iNeighbor = 0ul; iNeighbor < Third_Neighbor_Points.size(); iNeighbor++) { - for (auto jNeighbor = iNeighbor + 1; jNeighbor < Third_Neighbor_Points.size(); jNeighbor++) { - /*--- Repeated third neighbor with different origin ---*/ - - if ((Third_Neighbor_Points[iNeighbor] == Third_Neighbor_Points[jNeighbor]) && - (Third_Origin_Points[iNeighbor] != Third_Origin_Points[jNeighbor])) { - Suitable_Indirect_Neighbors.push_back(Third_Neighbor_Points[iNeighbor]); - } - } - } - - /*--- Remove duplicates from the final list of Suitable Indirect Neighbors. ---*/ - - sort(Suitable_Indirect_Neighbors.begin(), Suitable_Indirect_Neighbors.end()); - auto it2 = unique(Suitable_Indirect_Neighbors.begin(), Suitable_Indirect_Neighbors.end()); - Suitable_Indirect_Neighbors.resize(it2 - Suitable_Indirect_Neighbors.begin()); } void CMultiGridGeometry::SetPoint_Connectivity(const CGeometry* fine_grid) { /*--- Temporary, CPoint (nodes) then compresses this structure. ---*/ - vector > points(nPoint); + vector> points(nPoint); for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) { /*--- For each child CV (of the fine grid), ---*/ @@ -741,17 +967,14 @@ void CMultiGridGeometry::SetPoint_Connectivity(const CGeometry* fine_grid) { } void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* config) { - unsigned long iVertex, iFinePoint, iCoarsePoint; - unsigned short iMarker, iMarker_Tag, iChildren; - nMarker = fine_grid->GetnMarker(); unsigned short nMarker_Max = config->GetnMarker_Max(); /*--- If any children node belong to the boundary then the entire control volume will belong to the boundary ---*/ - for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) - for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { - iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); if (fine_grid->nodes->GetBoundary(iFinePoint)) { nodes->SetBoundary(iCoarsePoint, nMarker); break; @@ -762,20 +985,20 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co nVertex = new unsigned long[nMarker]; Tag_to_Marker = new string[nMarker_Max]; - for (iMarker_Tag = 0; iMarker_Tag < nMarker_Max; iMarker_Tag++) + for (auto iMarker_Tag = 0u; iMarker_Tag < nMarker_Max; iMarker_Tag++) Tag_to_Marker[iMarker_Tag] = fine_grid->GetMarker_Tag(iMarker_Tag); /*--- Compute the number of vertices to do the dimensionalization ---*/ - for (iMarker = 0; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0; + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0; - for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) { + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) { if (nodes->GetBoundary(iCoarsePoint)) { - for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { - iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); - for (iMarker = 0; iMarker < nMarker; iMarker++) { + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { if ((fine_grid->nodes->GetVertex(iFinePoint, iMarker) != -1) && (nodes->GetVertex(iCoarsePoint, iMarker) == -1)) { - iVertex = nVertex[iMarker]; + auto iVertex = nVertex[iMarker]; nodes->SetVertex(iCoarsePoint, iVertex, iMarker); nVertex[iMarker]++; } @@ -784,25 +1007,25 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co } } - for (iMarker = 0; iMarker < nMarker; iMarker++) { + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { vertex[iMarker] = new CVertex*[fine_grid->GetnVertex(iMarker) + 1]; nVertex[iMarker] = 0; } - for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) if (nodes->GetBoundary(iCoarsePoint)) - for (iMarker = 0; iMarker < nMarker; iMarker++) nodes->SetVertex(iCoarsePoint, -1, iMarker); + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) nodes->SetVertex(iCoarsePoint, -1, iMarker); - for (iMarker = 0; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0; + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) nVertex[iMarker] = 0; - for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) { + for (auto iCoarsePoint = 0ul; iCoarsePoint < nPoint; iCoarsePoint++) { if (nodes->GetBoundary(iCoarsePoint)) { - for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { - iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); - for (iMarker = 0; iMarker < fine_grid->GetnMarker(); iMarker++) { + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + for (auto iMarker = 0u; iMarker < fine_grid->GetnMarker(); iMarker++) { if ((fine_grid->nodes->GetVertex(iFinePoint, iMarker) != -1) && (nodes->GetVertex(iCoarsePoint, iMarker) == -1)) { - iVertex = nVertex[iMarker]; + auto iVertex = nVertex[iMarker]; vertex[iMarker][iVertex] = new CVertex(iCoarsePoint, nDim); nodes->SetVertex(iCoarsePoint, iVertex, iMarker); @@ -819,15 +1042,13 @@ void CMultiGridGeometry::SetVertex(const CGeometry* fine_grid, const CConfig* co } void CMultiGridGeometry::MatchActuator_Disk(const CConfig* config) { - unsigned short iMarker; - unsigned long iVertex, iPoint; int iProcessor = size; - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) == ACTDISK_INLET) || (config->GetMarker_All_KindBC(iMarker) == ACTDISK_OUTLET)) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - iPoint = vertex[iMarker][iVertex]->GetNode(); + for (auto iVertex = 0u; iVertex < nVertex[iMarker]; iVertex++) { + auto iPoint = vertex[iMarker][iVertex]->GetNode(); if (nodes->GetDomain(iPoint)) { vertex[iMarker][iVertex]->SetDonorPoint(iPoint, nodes->GetGlobalIndex(iPoint), iVertex, iMarker, iProcessor); } @@ -837,20 +1058,18 @@ void CMultiGridGeometry::MatchActuator_Disk(const CConfig* config) { } void CMultiGridGeometry::MatchPeriodic(const CConfig* config, unsigned short val_periodic) { - unsigned short iMarker, iPeriodic, nPeriodic; - unsigned long iVertex, iPoint; int iProcessor = rank; /*--- Evaluate the number of periodic boundary conditions ---*/ - nPeriodic = config->GetnMarker_Periodic(); + auto nPeriodic = config->GetnMarker_Periodic(); - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) { - iPeriodic = config->GetMarker_All_PerBound(iMarker); + auto iPeriodic = config->GetMarker_All_PerBound(iMarker); if ((iPeriodic == val_periodic) || (iPeriodic == val_periodic + nPeriodic / 2)) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - iPoint = vertex[iMarker][iVertex]->GetNode(); + for (auto iVertex = 0u; iVertex < nVertex[iMarker]; iVertex++) { + auto iPoint = vertex[iMarker][iVertex]->GetNode(); if (nodes->GetDomain(iPoint)) { vertex[iMarker][iVertex]->SetDonorPoint(iPoint, nodes->GetGlobalIndex(iPoint), iVertex, iMarker, iProcessor); @@ -863,18 +1082,12 @@ void CMultiGridGeometry::MatchPeriodic(const CConfig* config, unsigned short val void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned short action) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { - unsigned long iFinePoint, iCoarsePoint, iEdge, iParent; - long FineEdge, CoarseEdge; - unsigned short iChildren; - bool change_face_orientation; - su2double Coarse_Volume, Area; - /*--- Compute the area of the coarse volume ---*/ - for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) { + for (auto iCoarsePoint = 0u; iCoarsePoint < nPoint; iCoarsePoint++) { nodes->SetVolume(iCoarsePoint, 0.0); - Coarse_Volume = 0.0; - for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { - iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + su2double Coarse_Volume = 0.0; + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); Coarse_Volume += fine_grid->nodes->GetVolume(iFinePoint); } nodes->SetVolume(iCoarsePoint, Coarse_Volume); @@ -885,19 +1098,19 @@ void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned s edges->SetZeroValues(); } - for (iCoarsePoint = 0; iCoarsePoint < nPoint; iCoarsePoint++) - for (iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { - iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + for (auto iCoarsePoint = 0u; iCoarsePoint < nPoint; iCoarsePoint++) + for (auto iChildren = 0u; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); for (auto iFinePoint_Neighbor : fine_grid->nodes->GetPoints(iFinePoint)) { - iParent = fine_grid->nodes->GetParent_CV(iFinePoint_Neighbor); + auto iParent = fine_grid->nodes->GetParent_CV(iFinePoint_Neighbor); if ((iParent != iCoarsePoint) && (iParent < iCoarsePoint)) { - FineEdge = fine_grid->FindEdge(iFinePoint, iFinePoint_Neighbor); + auto FineEdge = fine_grid->FindEdge(iFinePoint, iFinePoint_Neighbor); - change_face_orientation = false; + bool change_face_orientation = false; if (iFinePoint < iFinePoint_Neighbor) change_face_orientation = true; - CoarseEdge = FindEdge(iParent, iCoarsePoint); + auto CoarseEdge = FindEdge(iParent, iCoarsePoint); const auto Normal = fine_grid->edges->GetNormal(FineEdge); @@ -912,9 +1125,9 @@ void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned s /*--- Check if there is a normal with null area ---*/ - for (iEdge = 0; iEdge < nEdge; iEdge++) { + for (auto iEdge = 0u; iEdge < nEdge; iEdge++) { const auto NormalFace = edges->GetNormal(iEdge); - Area = GeometryToolbox::Norm(nDim, NormalFace); + su2double Area = GeometryToolbox::Norm(nDim, NormalFace); if (Area == 0.0) { su2double DefaultNormal[3] = {EPS * EPS}; edges->SetNormal(iEdge, DefaultNormal); @@ -926,24 +1139,23 @@ void CMultiGridGeometry::SetControlVolume(const CGeometry* fine_grid, unsigned s void CMultiGridGeometry::SetBoundControlVolume(const CGeometry* fine_grid, const CConfig* config, unsigned short action) { - unsigned long iCoarsePoint, iFinePoint, FineVertex, iVertex; - su2double Normal[MAXNDIM] = {0.0}, Area, *NormalFace = nullptr; + su2double Normal[MAXNDIM] = {0.0}, *NormalFace = nullptr; if (action != ALLOCATE) { SU2_OMP_FOR_DYN(1) - for (auto iMarker = 0; iMarker < nMarker; iMarker++) - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) vertex[iMarker][iVertex]->SetZeroValues(); + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) vertex[iMarker][iVertex]->SetZeroValues(); END_SU2_OMP_FOR } SU2_OMP_FOR_DYN(1) - for (auto iMarker = 0; iMarker < nMarker; iMarker++) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - iCoarsePoint = vertex[iMarker][iVertex]->GetNode(); + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) { + auto iCoarsePoint = vertex[iMarker][iVertex]->GetNode(); for (auto iChildren = 0; iChildren < nodes->GetnChildren_CV(iCoarsePoint); iChildren++) { - iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); + auto iFinePoint = nodes->GetChildren_CV(iCoarsePoint, iChildren); if (fine_grid->nodes->GetVertex(iFinePoint, iMarker) != -1) { - FineVertex = fine_grid->nodes->GetVertex(iFinePoint, iMarker); + auto FineVertex = fine_grid->nodes->GetVertex(iFinePoint, iMarker); fine_grid->vertex[iMarker][FineVertex]->GetNormal(Normal); vertex[iMarker][iVertex]->AddNormal(Normal); } @@ -954,10 +1166,10 @@ void CMultiGridGeometry::SetBoundControlVolume(const CGeometry* fine_grid, const /*--- Check if there is a normal with null area ---*/ SU2_OMP_FOR_DYN(1) - for (auto iMarker = 0; iMarker < nMarker; iMarker++) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) { NormalFace = vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, NormalFace); + su2double Area = GeometryToolbox::Norm(nDim, NormalFace); if (Area == 0.0) for (auto iDim = 0; iDim < nDim; iDim++) NormalFace[iDim] = EPS * EPS; } @@ -1046,36 +1258,32 @@ void CMultiGridGeometry::SetRestricted_GridVelocity(const CGeometry* fine_grid) } void CMultiGridGeometry::FindNormal_Neighbor(const CConfig* config) { - unsigned short iMarker, iDim; - unsigned long iPoint, iVertex; - - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE && config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY && config->GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) { - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - iPoint = vertex[iMarker][iVertex]->GetNode(); + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) { + auto iPoint = vertex[iMarker][iVertex]->GetNode(); /*--- If the node belong to the domain ---*/ if (nodes->GetDomain(iPoint)) { /*--- Compute closest normal neighbor ---*/ - su2double cos_max, scalar_prod, norm_vect, norm_Normal, cos_alpha, diff_coord; unsigned long Point_Normal = 0; su2double* Normal = vertex[iMarker][iVertex]->GetNormal(); - cos_max = -1.0; + su2double cos_max = -1.0; for (auto jPoint : nodes->GetPoints(iPoint)) { - scalar_prod = 0.0; - norm_vect = 0.0; - norm_Normal = 0.0; - for (iDim = 0; iDim < nDim; iDim++) { - diff_coord = nodes->GetCoord(jPoint, iDim) - nodes->GetCoord(iPoint, iDim); + su2double scalar_prod = 0.0; + su2double norm_vect = 0.0; + su2double norm_Normal = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) { + su2double diff_coord = nodes->GetCoord(jPoint, iDim) - nodes->GetCoord(iPoint, iDim); scalar_prod += diff_coord * Normal[iDim]; norm_vect += diff_coord * diff_coord; norm_Normal += Normal[iDim] * Normal[iDim]; } norm_vect = sqrt(norm_vect); norm_Normal = sqrt(norm_Normal); - cos_alpha = scalar_prod / (norm_vect * norm_Normal); + su2double cos_alpha = scalar_prod / (norm_vect * norm_Normal); /*--- Get maximum cosine (not minimum because normals are oriented inwards) ---*/ if (cos_alpha >= cos_max) { @@ -1089,3 +1297,67 @@ void CMultiGridGeometry::FindNormal_Neighbor(const CConfig* config) { } } } + +su2double CMultiGridGeometry::ComputeLocalCurvature(const CGeometry* fine_grid, unsigned long iPoint, + unsigned short iMarker) const { + /*--- Compute local curvature (maximum angle between adjacent face normals) at a boundary vertex. + This is used to determine if agglomeration is safe based on a curvature threshold. ---*/ + + /*--- Get the vertex index for this point on this marker ---*/ + long iVertex = fine_grid->nodes->GetVertex(iPoint, iMarker); + if (iVertex < 0) return 0.0; // Point not on this marker + + /*--- Get the normal at this vertex ---*/ + su2double Normal_i[MAXNDIM] = {0.0}; + fine_grid->vertex[iMarker][iVertex]->GetNormal(Normal_i); + su2double Area_i = GeometryToolbox::Norm(int(nDim), Normal_i); + + if (Area_i < EPS) return 0.0; // Skip degenerate vertices + + /*--- Normalize the normal ---*/ + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Normal_i[iDim] /= Area_i; + } + + /*--- Find maximum angle with neighboring vertices on the same marker ---*/ + su2double max_angle = 0.0; + + /*--- Loop over edges connected to this point ---*/ + for (unsigned short iEdge = 0; iEdge < fine_grid->nodes->GetnPoint(iPoint); iEdge++) { + unsigned long jPoint = fine_grid->nodes->GetPoint(iPoint, iEdge); + + /*--- Check if neighbor is also on this marker ---*/ + long jVertex = fine_grid->nodes->GetVertex(jPoint, iMarker); + if (jVertex < 0) continue; // Not on this marker + + /*--- Get normal at neighbor vertex ---*/ + su2double Normal_j[MAXNDIM] = {0.0}; + fine_grid->vertex[iMarker][jVertex]->GetNormal(Normal_j); + su2double Area_j = GeometryToolbox::Norm(int(nDim), Normal_j); + + if (Area_j < EPS) continue; // Skip degenerate neighbor + + /*--- Normalize the neighbor normal ---*/ + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Normal_j[iDim] /= Area_j; + } + + /*--- Compute dot product: cos(angle) = n_i · n_j ---*/ + su2double dot_product = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + dot_product += Normal_i[iDim] * Normal_j[iDim]; + } + + /*--- Clamp to [-1, 1] to avoid numerical issues with acos ---*/ + dot_product = max(-1.0, min(1.0, dot_product)); + + /*--- Compute angle in degrees ---*/ + su2double angle_rad = acos(dot_product); + su2double angle_deg = angle_rad * 180.0 / PI_NUMBER; + + /*--- Track maximum angle ---*/ + max_angle = max(max_angle, angle_deg); + } + + return max_angle; +} diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index 361792eac5e..b5bc0dfdf35 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -25,6 +25,7 @@ * License along with SU2. If not, see . */ +#include #include "../../../include/geometry/dual_grid/CPoint.hpp" #include "../../../include/CConfig.hpp" #include "../../../include/parallelization/omp_structure.hpp" @@ -73,7 +74,7 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) { /*--- Multigrid structures. ---*/ if (config->GetnMGLevels() > 0) { - Parent_CV.resize(npoint) = 0; + Parent_CV.resize(npoint) = ULONG_MAX; // Initialize to ULONG_MAX to indicate unagglomerated Agglomerate.resize(npoint) = false; Agglomerate_Indirect.resize(npoint) = false; /*--- The finest grid does not have children CV's. ---*/ diff --git a/SU2_CFD/include/integration/CMultiGridIntegration.hpp b/SU2_CFD/include/integration/CMultiGridIntegration.hpp index eb2b7ded797..39e99f27bd5 100644 --- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp +++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp @@ -168,7 +168,7 @@ class CMultiGridIntegration final : public CIntegration { * \param[in] InclSharedDomain - Include the shared domain in the interpolation. */ void SetRestricted_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config); + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, unsigned short iMesh = 999); /*! * \brief Initialize the adjoint solution using the primal problem. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp index b61d3d61f33..26c4f8be33c 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp @@ -578,7 +578,7 @@ class CFVMFlowSolverBase : public CSolver { } - /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is diferent from 0. + /*--- Recompute the unsteady time step for the dual time strategy if the unsteady CFL is different from 0. * This is only done once because in dual time the time step cannot be variable. ---*/ if (dual_time && (Iteration == config->GetRestart_Iter()) && (config->GetUnst_CFL() != 0.0) && (iMesh == MESH_0)) { diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 961ce04832f..77c093f6c10 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -810,7 +810,7 @@ void CDriver::InitializeGeometryFVM(CConfig *config, CGeometry **&geometry) { geometry[MESH_0]->SetMGLevel(MESH_0); if ((config->GetnMGLevels() != 0) && (rank == MASTER_NODE)) - cout << "Setting the multigrid structure." << endl; + cout << "Setting the multigrid structure. NMG="<< config->GetnMGLevels() << endl; /*--- Loop over all the new grid ---*/ diff --git a/SU2_CFD/src/integration/CMultiGridIntegration.cpp b/SU2_CFD/src/integration/CMultiGridIntegration.cpp index b616f0e6649..2acca70f839 100644 --- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp +++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp @@ -28,6 +28,60 @@ #include "../../include/integration/CMultiGridIntegration.hpp" #include "../../../Common/include/parallelization/omp_structure.hpp" +namespace { +/*! + * \brief Helper function to enforce Euler wall BC by projecting momentum to tangent plane. + * \param[in] geo_coarse - Coarse grid geometry. + * \param[in] config - Problem configuration. + * \param[in,out] sol_coarse - Coarse grid solver (to access and modify solution/correction). + * \param[in] use_solution_old - If true, project Solution_Old (corrections); if false, project Solution. + */ +void ProjectEulerWallToTangentPlane(CGeometry* geo_coarse, const CConfig* config, CSolver* sol_coarse, + bool use_solution_old) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) == EULER_WALL) { + + SU2_OMP_FOR_STAT(32) + for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + const auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + + if (!geo_coarse->nodes->GetDomain(Point_Coarse)) continue; + + /*--- Get coarse grid normal ---*/ + su2double Normal[3] = {0.0}; + geo_coarse->vertex[iMarker][iVertex]->GetNormal(Normal); + const auto nDim = geo_coarse->GetnDim(); + su2double Area = GeometryToolbox::Norm(nDim, Normal); + + if (Area < EPS) continue; + + /*--- Normalize normal vector ---*/ + su2double UnitNormal[3] = {0.0}; + for (auto iDim = 0u; iDim < nDim; iDim++) { + UnitNormal[iDim] = Normal[iDim] / Area; + } + + /*--- Get current solution or correction ---*/ + su2double* solution_coarse = use_solution_old ? + sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse) : + sol_coarse->GetNodes()->GetSolution(Point_Coarse); + + /*--- Compute normal component of momentum (v·n or correction·n) ---*/ + su2double momentum_n = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) { + momentum_n += solution_coarse[iDim + 1] * UnitNormal[iDim]; + } + + /*--- Project to tangent plane: solution_coarse -= (solution_coarse·n)n ---*/ + for (auto iDim = 0u; iDim < nDim; iDim++) { + solution_coarse[iDim + 1] -= momentum_n * UnitNormal[iDim]; + } + } + END_SU2_OMP_FOR + } + } +} +} // anonymous namespace CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { } @@ -65,6 +119,17 @@ void CMultiGridIntegration::MultiGrid_Iteration(CGeometry ****geometry, const unsigned short Solver_Position = config[iZone]->GetContainerPosition(RunTime_EqSystem); /*--- Start an OpenMP parallel region covering the entire MG iteration, if the solver supports it. ---*/ + /*--- Set thread affinity for deterministic execution ---*/ + +#ifdef _OPENMP + /*--- Store original thread affinity settings ---*/ + static bool affinity_set = false; + if (!affinity_set) { + // This ensures threads are bound consistently across runs + omp_set_schedule(omp_sched_static, 0); + affinity_set = true; + } +#endif SU2_OMP_PARALLEL_(if(solver_container[iZone][iInst][MESH_0][Solver_Position]->GetHasHybridParallel())) { @@ -170,6 +235,9 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, for (unsigned short iPreSmooth = 0; iPreSmooth < config->GetMG_PreSmooth(iMesh); iPreSmooth++) { + /*--- Synchronize before each smoothing iteration ---*/ + SU2_OMP_BARRIER + /*--- Time and space integration ---*/ for (unsigned short iRKStep = 0; iRKStep < iRKLimit; iRKStep++) { @@ -211,6 +279,13 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, } + /*--- MPI sync after RK stage to ensure halos have updated solution for next smoothing iteration ---*/ + solver_fine->InitiateComms(geometry_fine, config, MPI_QUANTITIES::SOLUTION); + solver_fine->CompleteComms(geometry_fine, config, MPI_QUANTITIES::SOLUTION); + + /*--- OpenMP memory fence to ensure all threads see consistent state ---*/ + SU2_OMP_BARRIER + } /*--- Compute Forcing Term $P_(k+1) = I^(k+1)_k(P_k+F_k(u_k))-F_(k+1)(I^(k+1)_k u_k)$ and update solution for multigrid ---*/ @@ -240,7 +315,7 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, /*--- Compute $r_(k+1) = F_(k+1)(I^(k+1)_k u_k)$ ---*/ - SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + SetRestricted_Solution(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config, iMesh); solver_coarse->Preprocessing(geometry_coarse, solver_container_coarse, config, iMesh+1, NO_RK_ITER, RunTime_EqSystem, false); @@ -257,6 +332,8 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, } /*--- Recursive call to MultiGrid_Cycle (this routine). ---*/ + /*--- Execute multigrid cycles sequentially to ensure deterministic recursion order ---*/ + /*--- This prevents accumulation of floating-point variations across recursive calls ---*/ for (unsigned short imu = 0; imu <= RecursiveParam; imu++) { @@ -264,14 +341,215 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, if (iMesh == config->GetnMGLevels()-2) nextRecurseParam = 0; + /*--- Synchronize prior to recursive calls for determinism ---*/ + SU2_OMP_BARRIER + MultiGrid_Cycle(geometry, solver_container, numerics_container, config_container, iMesh+1, nextRecurseParam, RunTime_EqSystem, iZone, iInst); + + /*--- Synchronize after each multigrid cycle ---*/ + SU2_OMP_BARRIER } /*--- Compute prolongated solution, and smooth the correction $u^(new)_k = u_k + Smooth(I^k_(k+1)(u_(k+1)-I^(k+1)_k u_k))$ ---*/ GetProlongated_Correction(RunTime_EqSystem, solver_fine, solver_coarse, geometry_fine, geometry_coarse, config); + /*--- Get current CFL values ---*/ + su2double CFL_fine = config->GetCFL(iMesh); + su2double CFL_coarse_current = config->GetCFL(iMesh+1); + su2double current_coeff = CFL_coarse_current / CFL_fine; + + /*--- Adaptive CFL using Exponential Moving Average (EMA) ---*/ + /*--- All operations performed in master thread for determinism ---*/ + constexpr int AVG_WINDOW = 5; + constexpr int MAX_MG_LEVELS = 10; + + /*--- Use static storage shared across all threads for deterministic CFL adaptation ---*/ + /*--- Each MPI rank maintains its own independent adaptive CFL state ---*/ + /*--- Note: MPI ranks do not share memory, so static is safe for MPI parallelism ---*/ + static su2double current_avg[MAX_MG_LEVELS] = {}; + static su2double prev_avg[MAX_MG_LEVELS] = {}; + static su2double last_res[MAX_MG_LEVELS] = {}; + static bool last_was_increase[MAX_MG_LEVELS] = {}; + static int oscillation_count[MAX_MG_LEVELS] = {}; + static unsigned long last_check_iter[MAX_MG_LEVELS] = {}; + static unsigned long last_update_iter[MAX_MG_LEVELS] = {}; + + /*--- Only master thread performs CFL adaptation to ensure deterministic behavior ---*/ + /*--- All adaptive CFL state and computation must be done by a single thread ---*/ + su2double CFL_coarse_new = CFL_coarse_current; // Default: keep current value + + SU2_OMP_MASTER + { + /*--- Get global iteration count first ---*/ + unsigned long current_iter; + if (config->GetTime_Domain()) + current_iter = config->GetTimeIter(); + else + current_iter = config->GetInnerIter(); + + /*--- Reset state at the beginning of a new solve (iter 0 or 1) ---*/ + /*--- This ensures deterministic behavior across multiple runs ---*/ + static unsigned long last_reset_iter = ULONG_MAX; + if (current_iter <= 1 && last_reset_iter != current_iter) { + for (int i = 0; i < MAX_MG_LEVELS; i++) { + current_avg[i] = 0.0; + prev_avg[i] = 0.0; + last_res[i] = 0.0; + last_was_increase[i] = false; + oscillation_count[i] = 0; + last_check_iter[i] = 0; + last_update_iter[i] = 0; + } + last_reset_iter = current_iter; + } + + unsigned short lvl = min(iMesh, (unsigned short)(MAX_MG_LEVELS - 1)); + unsigned long iter = current_iter; + + /*--- Get sum of all RMS residuals for all variables (local to this rank) ---*/ + su2double rms_res_coarse_local = 0.0; + for (unsigned short iVar = 0; iVar < solver_coarse->GetnVar(); iVar++) { + rms_res_coarse_local += solver_coarse->GetRes_RMS(iVar); + } + + /*--- MPI synchronization: ensure all ranks use the same global residual value ---*/ + /*--- This is critical for consistent CFL adaptation across all ranks ---*/ + su2double rms_res_coarse = rms_res_coarse_local; + +#ifdef HAVE_MPI + /*--- For coarse grids, residuals are not globally reduced by default ---*/ + /*--- We need to synchronize them for consistent adaptive CFL decisions ---*/ + if (geometry_coarse->GetMGLevel() > 0) { + su2double rms_global_sum = 0.0; + SU2_MPI::Allreduce(&rms_res_coarse_local, &rms_global_sum, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + rms_res_coarse = rms_global_sum / static_cast(SU2_MPI::GetSize()); + } +#endif + + /*--- Flip-flop detection: detect oscillating residuals (once per outer iteration) ---*/ + bool oscillation_detected = false; + if (iter != last_check_iter[lvl]) { + last_check_iter[lvl] = iter; + + if (last_res[lvl] > EPS) { + bool current_is_increase = (rms_res_coarse > last_res[lvl]); + if (current_is_increase != last_was_increase[lvl]) { + /*--- Direction changed, increment oscillation counter ---*/ + oscillation_count[lvl]++; + if (oscillation_count[lvl] >= 4) { + /*--- Detected 4 consecutive direction changes = oscillation ---*/ + oscillation_detected = true; + oscillation_count[lvl] = 0; // Reset counter after detecting + } + } else { + /*--- Same direction, reset counter ---*/ + oscillation_count[lvl] = 0; + } + last_was_increase[lvl] = current_is_increase; + } + last_res[lvl] = rms_res_coarse; + } + + /*--- Update exponential moving average ---*/ + if (current_avg[lvl] < EPS) { + current_avg[lvl] = rms_res_coarse; // Initialize with first value + } else { + current_avg[lvl] = (current_avg[lvl] * (AVG_WINDOW - 1) + rms_res_coarse) / AVG_WINDOW; + } + + /*--- Check if we should compare and adapt CFL ---*/ + su2double new_coeff = current_coeff; + const su2double MIN_REDUCTION_FACTOR = 0.98; // Require at least 2% reduction + const int UPDATE_INTERVAL = 5; // Update reference every N iterations + + /*--- Initialize prev_avg on first use ---*/ + if (prev_avg[lvl] < EPS) { + prev_avg[lvl] = current_avg[lvl]; + } + + /*--- Periodically update prev_avg to allow ratio to reflect accumulated decrease ---*/ + bool should_update = (iter - last_update_iter[lvl] >= UPDATE_INTERVAL); + + /*--- Asymmetric adaptation for robustness ---*/ + if (prev_avg[lvl] > EPS) { + su2double ratio = current_avg[lvl] / prev_avg[lvl]; + bool sufficient_decrease = (ratio < MIN_REDUCTION_FACTOR); + bool increasing_trend = (ratio >= 1.0); + + if (increasing_trend) { + /*--- Residual increasing: reduce CFL immediately for robustness ---*/ + new_coeff = current_coeff * 0.90; + /*--- Update reference since we're reacting immediately ---*/ + prev_avg[lvl] = current_avg[lvl]; + last_update_iter[lvl] = iter; + } else if (sufficient_decrease && should_update) { + /*--- Residual decreasing sufficiently: increase CFL ---*/ + new_coeff = current_coeff * 1.05; + /*--- Update reference only when we actually increase CFL ---*/ + prev_avg[lvl] = current_avg[lvl]; + last_update_iter[lvl] = iter; + } + } + + /*--- CFL reduction for oscillation detection ---*/ + if (oscillation_detected) { + new_coeff = current_coeff * 0.75; + /*--- Update reference after oscillation response ---*/ + prev_avg[lvl] = current_avg[lvl]; + last_update_iter[lvl] = iter; + } + + /*--- Clamp coefficient between 0.5 and 1.0 ---*/ + new_coeff = max(0.5, min(1.0, new_coeff)); + + /*--- Update coarse grid CFL ---*/ + CFL_coarse_new = max(0.5 * CFL_fine, min(CFL_fine, CFL_fine * new_coeff)); + +#ifdef HAVE_MPI + /*--- Ensure all ranks use the same CFL value (broadcast from rank 0) ---*/ + SU2_MPI::Bcast(&CFL_coarse_new, 1, MPI_DOUBLE, 0, SU2_MPI::GetComm()); +#endif + + /*--- Update the shared config object ---*/ + config->SetCFL(iMesh+1, CFL_coarse_new); + } + END_SU2_OMP_MASTER + /*--- Implicit barrier at end of master region ensures all threads see updated CFL ---*/ + + /*--- Update LocalCFL at each coarse grid point ---*/ + SU2_OMP_FOR_STAT(roundUpDiv(geometry_coarse->GetnPoint(), omp_get_num_threads())) + for (auto iPoint = 0ul; iPoint < geometry_coarse->GetnPoint(); iPoint++) { + solver_coarse->GetNodes()->SetLocalCFL(iPoint, CFL_coarse_new); + } + END_SU2_OMP_FOR + + /*--- Output monitoring information periodically ---*/ +// +////#ifdef HAVE_MPI +//// /*--- Synchronize all ranks before output to ensure consistent state ---*/ +//// SU2_MPI::Barrier(SU2_MPI::GetComm()); +////#endif +/* + if (SU2_MPI::GetRank() == 0 && iter % 1 == 0) { + bool cfl_increased = (new_coeff < current_coeff); // Lower coeff = higher CFL + bool cfl_decreased = (new_coeff > current_coeff); // Higher coeff = lower CFL + string action = "unchanged"; + if (cfl_increased) action = "increased"; + else if (cfl_decreased) action = "decreased"; + + su2double percent_reduction = (1.0 - ratio_for_display) * 100.0; + + cout << " MG Level " << iMesh+1 + << ", iter = " << iter + << ", percent reduction = " << percent_reduction + << ", action = " << action + << ", oscillation = " << (oscillation_detected ? "yes" : "no") + << ", coeff = " << new_coeff << " (was " << current_coeff << ")" + << ", CFL = " << CFL_coarse_new << endl; + } +*/ SmoothProlongated_Correction(RunTime_EqSystem, solver_fine, geometry_fine, config->GetMG_CorrecSmooth(iMesh), 1.25, config); SetProlongated_Correction(solver_fine, geometry_fine, config, iMesh); @@ -281,6 +559,9 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, for (unsigned short iPostSmooth = 0; iPostSmooth < config->GetMG_PostSmooth(iMesh); iPostSmooth++) { + /*--- Synchronize before each post-smoothing iteration ---*/ + SU2_OMP_BARRIER + for (unsigned short iRKStep = 0; iRKStep < iRKLimit; iRKStep++) { solver_fine->Preprocessing(geometry_fine, solver_container_fine, config, iMesh, iRKStep, RunTime_EqSystem, false); @@ -300,6 +581,11 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, solver_fine->Postprocessing(geometry_fine, solver_container_fine, config, iMesh); } + + /*--- MPI sync after RK stage to ensure halos have updated solution for next smoothing iteration ---*/ + solver_fine->InitiateComms(geometry_fine, config, MPI_QUANTITIES::SOLUTION); + solver_fine->CompleteComms(geometry_fine, config, MPI_QUANTITIES::SOLUTION); + } } @@ -307,9 +593,6 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry, void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { - unsigned long Point_Fine, Point_Coarse, iVertex; - unsigned short iMarker, iChildren, iVar; - su2double Area_Parent, Area_Children; const su2double *Solution_Fine = nullptr, *Solution_Coarse = nullptr; const unsigned short nVar = sol_coarse->GetnVar(); @@ -317,41 +600,50 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS auto *Solution = new su2double[nVar]; SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { - Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); + su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); - for (iVar = 0; iVar < nVar; iVar++) Solution[iVar] = 0.0; + for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] = 0.0; - for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - Area_Children = geo_fine->nodes->GetVolume(Point_Fine); + /*--- Accumulate children contributions with stable ordering ---*/ + /*--- Process all children in sequential order to ensure deterministic FP summation ---*/ + auto nChildren = geo_coarse->nodes->GetnChildren_CV(Point_Coarse); + for (auto iChildren = 0u; iChildren < nChildren; iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + su2double Area_Children = geo_fine->nodes->GetVolume(Point_Fine); Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine); - for (iVar = 0; iVar < nVar; iVar++) - Solution[iVar] -= Solution_Fine[iVar]*Area_Children/Area_Parent; + su2double weight = Area_Children / Area_Parent; + for (auto iVar = 0u; iVar < nVar; iVar++) + Solution[iVar] -= Solution_Fine[iVar] * weight; } Solution_Coarse = sol_coarse->GetNodes()->GetSolution(Point_Coarse); - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] += Solution_Coarse[iVar]; - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) sol_coarse->GetNodes()->SetSolution_Old(Point_Coarse,Solution); } END_SU2_OMP_FOR delete [] Solution; + /*--- Ensure all threads complete before MPI communication ---*/ + SU2_OMP_BARRIER + + /*--- Enforce Euler wall BC on corrections by projecting to tangent plane ---*/ + ProjectEulerWallToTangentPlane(geo_coarse, config, sol_coarse, true); + /*--- Remove any contributions from no-slip walls. ---*/ - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetViscous_Wall(iMarker)) { SU2_OMP_FOR_STAT(32) - for (iVertex = 0; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { - - Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); /*--- For dirichlet boundary condtions, set the correction to zero. Note that Solution_Old stores the correction not the actual value ---*/ @@ -369,10 +661,13 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD); sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD); + /*--- Ensure MPI completion visible to all threads ---*/ + SU2_OMP_BARRIER + SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { - for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); sol_fine->LinSysRes.SetBlock(Point_Fine, sol_coarse->GetNodes()->GetSolution_Old(Point_Coarse)); } } @@ -387,13 +682,11 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ if (val_nSmooth == 0) return; const su2double *Residual_Old, *Residual_Sum, *Residual_j; - unsigned short iVar, iSmooth, iMarker, iNeigh; - unsigned long iPoint, jPoint, iVertex; const unsigned short nVar = solver->GetnVar(); SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); iPoint++) { Residual_Old = solver->LinSysRes.GetBlock(iPoint); solver->GetNodes()->SetResidual_Old(iPoint,Residual_Old); } @@ -401,17 +694,19 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ /*--- Jacobi iterations. ---*/ - for (iSmooth = 0; iSmooth < val_nSmooth; iSmooth++) { + for (auto iSmooth = 0u; iSmooth < val_nSmooth; iSmooth++) { /*--- Loop over all mesh points (sum the residuals of direct neighbors). ---*/ + /*--- Use static scheduling to ensure deterministic iteration order for reproducibility ---*/ SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) - for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) { + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { solver->GetNodes()->SetResidualSumZero(iPoint); - for (iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { - jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); + /*--- Sum neighbor contributions in deterministic order ---*/ + for (auto iNeigh = 0u; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) { + auto jPoint = geometry->nodes->GetPoint(iPoint, iNeigh); Residual_j = solver->LinSysRes.GetBlock(jPoint); solver->GetNodes()->AddResidual_Sum(iPoint, Residual_j); } @@ -422,28 +717,28 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ /*--- Loop over all mesh points (update residuals with the neighbor averages). ---*/ SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPoint(), omp_get_num_threads())) - for (iPoint = 0; iPoint < geometry->GetnPoint(); ++iPoint) { + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); ++iPoint) { su2double factor = 1.0/(1.0+val_smooth_coeff*su2double(geometry->nodes->GetnPoint(iPoint))); Residual_Sum = solver->GetNodes()->GetResidual_Sum(iPoint); Residual_Old = solver->GetNodes()->GetResidual_Old(iPoint); - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) solver->LinSysRes(iPoint,iVar) = (Residual_Old[iVar] + val_smooth_coeff*Residual_Sum[iVar])*factor; } END_SU2_OMP_FOR /*--- Restore original residuals (without average) at boundary points. ---*/ - for (iMarker = 0; iMarker < geometry->GetnMarker(); iMarker++) { + for (auto iMarker = 0u; iMarker < geometry->GetnMarker(); iMarker++) { if ((config->GetMarker_All_KindBC(iMarker) != INTERNAL_BOUNDARY) && (config->GetMarker_All_KindBC(iMarker) != NEARFIELD_BOUNDARY) && (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { SU2_OMP_FOR_STAT(32) - for (iVertex = 0; iVertex < geometry->GetnVertex(iMarker); iVertex++) { - iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); Residual_Old = solver->GetNodes()->GetResidual_Old(iPoint); solver->LinSysRes.SetBlock(iPoint, Residual_Old); } @@ -457,22 +752,22 @@ void CMultiGridIntegration::SmoothProlongated_Correction(unsigned short RunTime_ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeometry *geo_fine, CConfig *config, unsigned short iMesh) { - unsigned long Point_Fine; - unsigned short iVar; su2double *Solution_Fine, *Residual_Fine; const unsigned short nVar = sol_fine->GetnVar(); const su2double factor = config->GetDamp_Correc_Prolong(); SU2_OMP_FOR_STAT(roundUpDiv(geo_fine->GetnPointDomain(), omp_get_num_threads())) - for (Point_Fine = 0; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) { + for (auto Point_Fine = 0ul; Point_Fine < geo_fine->GetnPointDomain(); Point_Fine++) { Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); Solution_Fine = sol_fine->GetNodes()->GetSolution(Point_Fine); - for (iVar = 0; iVar < nVar; iVar++) { + for (auto iVar = 0u; iVar < nVar; iVar++) { /*--- Prevent a fine grid divergence due to a coarse grid divergence ---*/ if (Residual_Fine[iVar] != Residual_Fine[iVar]) Residual_Fine[iVar] = 0.0; - Solution_Fine[iVar] += factor*Residual_Fine[iVar]; + + su2double correction = factor * Residual_Fine[iVar]; + Solution_Fine[iVar] += correction; } } END_SU2_OMP_FOR @@ -486,13 +781,11 @@ void CMultiGridIntegration::SetProlongated_Correction(CSolver *sol_fine, CGeomet void CMultiGridIntegration::SetProlongated_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { - unsigned long Point_Fine, Point_Coarse; - unsigned short iChildren; SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { - for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); sol_fine->GetNodes()->SetSolution(Point_Fine, sol_coarse->GetNodes()->GetSolution(Point_Coarse)); } } @@ -502,8 +795,6 @@ void CMultiGridIntegration::SetProlongated_Solution(unsigned short RunTime_EqSys void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, unsigned short iMesh) { - unsigned long Point_Fine, Point_Coarse, iVertex; - unsigned short iMarker, iVar, iChildren; const su2double *Residual_Fine; const unsigned short nVar = sol_coarse->GetnVar(); @@ -512,16 +803,15 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar auto *Residual = new su2double[nVar]; SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { sol_coarse->GetNodes()->SetRes_TruncErrorZero(Point_Coarse); - for (iVar = 0; iVar < nVar; iVar++) Residual[iVar] = 0.0; - - for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + for (auto iVar = 0u; iVar < nVar; iVar++) Residual[iVar] = 0.0; + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + auto Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); Residual_Fine = sol_fine->LinSysRes.GetBlock(Point_Fine); - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) Residual[iVar] += factor*Residual_Fine[iVar]; } sol_coarse->GetNodes()->AddRes_TruncError(Point_Coarse, Residual); @@ -530,11 +820,11 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar delete [] Residual; - for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) { if (config->GetViscous_Wall(iMarker)) { SU2_OMP_FOR_STAT(32) - for (iVertex = 0; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { - Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); + for (auto iVertex = 0ul; iVertex < geo_coarse->nVertex[iMarker]; iVertex++) { + auto Point_Coarse = geo_coarse->vertex[iMarker][iVertex]->GetNode(); sol_coarse->GetNodes()->SetVel_ResTruncError_Zero(Point_Coarse); } END_SU2_OMP_FOR @@ -542,7 +832,7 @@ void CMultiGridIntegration::SetForcing_Term(CSolver *sol_fine, CSolver *sol_coar } SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPointDomain(), omp_get_num_threads())) - for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPointDomain(); Point_Coarse++) { sol_coarse->GetNodes()->SubtractRes_TruncError(Point_Coarse, sol_coarse->LinSysRes.GetBlock(Point_Coarse)); } END_SU2_OMP_FOR @@ -553,7 +843,7 @@ void CMultiGridIntegration::SetResidual_Term(CGeometry *geometry, CSolver *solve AD::StartNoSharedReading(); SU2_OMP_FOR_STAT(roundUpDiv(geometry->GetnPointDomain(), omp_get_num_threads())) - for (unsigned long iPoint = 0; iPoint < geometry->GetnPointDomain(); iPoint++) + for (auto iPoint = 0ul; iPoint < geometry->GetnPointDomain(); iPoint++) solver->LinSysRes.AddBlock(iPoint, solver->GetNodes()->GetResTruncError(iPoint)); END_SU2_OMP_FOR AD::EndNoSharedReading(); @@ -561,7 +851,7 @@ void CMultiGridIntegration::SetResidual_Term(CGeometry *geometry, CSolver *solve } void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, - CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { + CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config, unsigned short iMesh) { const unsigned short Solver_Position = config->GetContainerPosition(RunTime_EqSystem); const bool grid_movement = config->GetGrid_Movement(); @@ -606,6 +896,9 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst } } + /*--- Enforce Euler wall BC by projecting velocity to tangent plane ---*/ + ProjectEulerWallToTangentPlane(geo_coarse, config, sol_coarse, false); + /*--- MPI the new interpolated solution ---*/ sol_coarse->InitiateComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION); @@ -615,39 +908,36 @@ void CMultiGridIntegration::SetRestricted_Solution(unsigned short RunTime_EqSyst void CMultiGridIntegration::SetRestricted_Gradient(unsigned short RunTime_EqSystem, CSolver *sol_fine, CSolver *sol_coarse, CGeometry *geo_fine, CGeometry *geo_coarse, CConfig *config) { - unsigned long Point_Fine, Point_Coarse; - unsigned short iVar, iDim, iChildren; - su2double Area_Parent, Area_Children; const unsigned short nDim = geo_coarse->GetnDim(); const unsigned short nVar = sol_coarse->GetnVar(); auto **Gradient = new su2double* [nVar]; - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) Gradient[iVar] = new su2double [nDim]; SU2_OMP_FOR_STAT(roundUpDiv(geo_coarse->GetnPoint(), omp_get_num_threads())) - for (Point_Coarse = 0; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { - Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); + for (auto Point_Coarse = 0ul; Point_Coarse < geo_coarse->GetnPoint(); Point_Coarse++) { + su2double Area_Parent = geo_coarse->nodes->GetVolume(Point_Coarse); - for (iVar = 0; iVar < nVar; iVar++) - for (iDim = 0; iDim < nDim; iDim++) + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) Gradient[iVar][iDim] = 0.0; - for (iChildren = 0; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { - Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); - Area_Children = geo_fine->nodes->GetVolume(Point_Fine); + for (auto iChildren = 0u; iChildren < geo_coarse->nodes->GetnChildren_CV(Point_Coarse); iChildren++) { + unsigned long Point_Fine = geo_coarse->nodes->GetChildren_CV(Point_Coarse, iChildren); + su2double Area_Children = geo_fine->nodes->GetVolume(Point_Fine); auto Gradient_fine = sol_fine->GetNodes()->GetGradient(Point_Fine); - for (iVar = 0; iVar < nVar; iVar++) - for (iDim = 0; iDim < nDim; iDim++) + for (auto iVar = 0u; iVar < nVar; iVar++) + for (auto iDim = 0u; iDim < nDim; iDim++) Gradient[iVar][iDim] += Gradient_fine[iVar][iDim]*Area_Children/Area_Parent; } sol_coarse->GetNodes()->SetGradient(Point_Coarse,Gradient); } END_SU2_OMP_FOR - for (iVar = 0; iVar < nVar; iVar++) + for (auto iVar = 0u; iVar < nVar; iVar++) delete [] Gradient[iVar]; delete [] Gradient; diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp index 6477fa6f361..7113ba92a51 100644 --- a/SU2_CFD/src/solvers/CIncNSSolver.cpp +++ b/SU2_CFD/src/solvers/CIncNSSolver.cpp @@ -739,6 +739,7 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); const auto Point_Normal = geometry->vertex[iMarker][iVertex]->GetNormal_Neighbor(); + /*--- On the finest mesh compute also on halo nodes to avoid communication of tau wall. ---*/ if ((!geometry->nodes->GetDomain(iPoint)) && !(MGLevel==MESH_0)) continue; diff --git a/SU2_CFD/src/variables/CAdjEulerVariable.cpp b/SU2_CFD/src/variables/CAdjEulerVariable.cpp index 565c80f8c77..b2ac36a5b8b 100644 --- a/SU2_CFD/src/variables/CAdjEulerVariable.cpp +++ b/SU2_CFD/src/variables/CAdjEulerVariable.cpp @@ -95,6 +95,9 @@ CAdjEulerVariable::CAdjEulerVariable(su2double psirho, const su2double *phi, su2 if (config->GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE) HB_Source.resize(nPoint,nVar) = su2double(0.0); + /*--- Allocate LocalCFL for multigrid support ---*/ + LocalCFL.resize(nPoint) = su2double(0.0); + if (config->GetMultizone_Problem()) Set_BGSSolution_k(); diff --git a/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref b/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref index 8b2887947f5..9a1f6b7d1e7 100644 --- a/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref +++ b/TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref @@ -1,39 +1,39 @@ VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP" - 0 , 292.459 , 0.001 - 1 , -8318.5 , 0.001 - 2 , -16158.4 , 0.001 - 3 , -21277.2 , 0.001 - 4 , -23366.9 , 0.001 - 5 , -22727.6 , 0.001 - 6 , -19872.3 , 0.001 - 7 , -15278.1 , 0.001 - 8 , -9283.99 , 0.001 - 9 , -2179.87 , 0.001 - 10 , 5520.63 , 0.001 - 11 , 12726.0 , 0.001 - 12 , 17538.0 , 0.001 - 13 , 17441.1 , 0.001 - 14 , 10335.6 , 0.001 - 15 , -2686.37 , 0.001 - 16 , -10522.5 , 0.001 - 17 , 24712.2 , 0.001 - 18 , 166438.0 , 0.001 - 19 , -15618.6 , 0.001 - 20 , -14178.7 , 0.001 - 21 , -12765.5 , 0.001 - 22 , -12007.2 , 0.001 - 23 , -13597.5 , 0.001 - 24 , -19002.9 , 0.001 - 25 , -28729.6 , 0.001 - 26 , -41946.6 , 0.001 - 27 , -56289.3 , 0.001 - 28 , -67832.1 , 0.001 - 29 , -71484.0 , 0.001 - 30 , -62334.5 , 0.001 - 31 , -38478.2 , 0.001 - 32 , -4757.34 , 0.001 - 33 , 26448.5 , 0.001 - 34 , 45049.5 , 0.001 - 35 , 60960.9 , 0.001 - 36 , 83515.9 , 0.001 - 37 , 8837.4 , 0.001 + 0 , 17852.4 , 0.001 + 1 , 17762.2 , 0.001 + 2 , 13023.3 , 0.001 + 3 , 7235.23 , 0.001 + 4 , 2025.24 , 0.001 + 5 , -2173.83 , 0.001 + 6 , -5516.67 , 0.001 + 7 , -8317.99 , 0.001 + 8 , -10814.3 , 0.001 + 9 , -13092.8 , 0.001 + 10 , -15126.6 , 0.001 + 11 , -16823.4 , 0.001 + 12 , -17994.0 , 0.001 + 13 , -18196.9 , 0.001 + 14 , -16453.7 , 0.001 + 15 , -10589.4 , 0.001 + 16 , 4413.36 , 0.001 + 17 , 30522.9 , 0.001 + 18 , 5834.86 , 0.001 + 19 , -18131.9 , 0.001 + 20 , -21995.1 , 0.001 + 21 , -20862.7 , 0.001 + 22 , -18622.4 , 0.001 + 23 , -18164.8 , 0.001 + 24 , -20842.9 , 0.001 + 25 , -26432.9 , 0.001 + 26 , -33388.6 , 0.001 + 27 , -39214.0 , 0.001 + 28 , -40910.6 , 0.001 + 29 , -35624.7 , 0.001 + 30 , -21761.1 , 0.001 + 31 , -718.652 , 0.001 + 32 , 21422.9 , 0.001 + 33 , 34942.8 , 0.001 + 34 , 34300.9 , 0.001 + 35 , 28342.4 , 0.001 + 36 , 23237.9 , 0.001 + 37 , -20230.5 , 0.001 diff --git a/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref b/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref index 1a899cb0668..64fffa595ba 100644 --- a/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref +++ b/TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref @@ -1,4 +1,4 @@ VARIABLES="VARIABLE" , "DRAG" , "EFFICIENCY" , "FORCE_X" , "FORCE_Y" , "FORCE_Z" , "LIFT" , "MOMENT_X" , "MOMENT_Y" , "MOMENT_Z" , "SIDEFORCE" - 0 , 0.24008251 , -117.3444057 , 0.2742430499 , -1.56293638 , 0.0 , -1.568547024 , 0.0 , 0.0 , 1.189018284 , 0.0 - 1 , 0.4064005433 , -189.77779 , 0.4586830911 , -2.391641926 , 0.0 , -2.401078899 , 0.0 , 0.0 , 1.030793484 , 0.0 - 2 , 0.5421052294 , -249.5397676 , 0.6095878335 , -3.086770277 , 0.0 , -3.099333798 , 0.0 , 0.0 , 0.6218682473 , 0.0 + 0 , 0.2000855462 , -50.31379686 , 0.2356436894 , -1.627423951 , 0.0 , -1.632177209 , 0.0 , 0.0 , 0.987584412 , 0.0 + 1 , 0.3003955622 , -70.66636166 , 0.348808693 , -2.215465435 , 0.0 , -2.222547435 , 0.0 , 0.0 , 0.9429020878 , 0.0 + 2 , 0.4098512262 , -88.09224197 , 0.4674108675 , -2.633450055 , 0.0 , -2.643019879 , 0.0 , 0.0 , 0.6725005135 , 0.0 diff --git a/TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref b/TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref index eca2a7ba183..9e7d23efcaa 100644 --- a/TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref +++ b/TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref @@ -1,5 +1,5 @@ VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP" - 0 , 0.00765473 , 0.0001 - 1 , 0.00497838 , 0.0001 - 2 , 0.0024697 , 0.0001 - 3 , 0.00090216 , 0.0001 + 0 , 0.00770575 , 0.0001 + 1 , 0.00508228 , 0.0001 + 2 , 0.00250931 , 0.0001 + 3 , 0.000896287 , 0.0001 diff --git a/TestCases/euler/channel/inv_channel_RK.cfg b/TestCases/euler/channel/inv_channel_RK.cfg index a719b413d73..bf374054e77 100644 --- a/TestCases/euler/channel/inv_channel_RK.cfg +++ b/TestCases/euler/channel/inv_channel_RK.cfg @@ -44,24 +44,24 @@ MARKER_MONITORING= ( upper_wall, lower_wall ) % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES -CFL_NUMBER= 2.5 +CFL_NUMBER= 1.0 CFL_ADAPT= NO CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 ) RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) ITER= 110 LINEAR_SOLVER= BCGSTAB -LINEAR_SOLVER_ERROR= 1E-6 -LINEAR_SOLVER_ITER= 5 +LINEAR_SOLVER_ERROR= 1E-1 +LINEAR_SOLVER_ITER= 10 % -------------------------- MULTIGRID PARAMETERS -----------------------------% % MGLEVEL= 3 MGCYCLE= W_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.9 -MG_DAMP_PROLONGATION= 0.9 +MG_PRE_SMOOTH= ( 4, 4, 4, 4 ) +MG_POST_SMOOTH= ( 4, 4, 4, 4 ) +MG_CORRECTION_SMOOTH= ( 4, 4, 4, 4 ) +MG_DAMP_RESTRICTION= 0.5 +MG_DAMP_PROLONGATION= 0.5 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py index 21f86015e12..9e84aec251a 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -66,7 +66,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [-3.254503, -3.495599, 0.052373, 0.000000] + discadj_arina2k.test_vals = [-3.229402, -3.466984, 0.043984, 0.000000] test_list.append(discadj_arina2k) #################################### @@ -99,7 +99,7 @@ def main(): discadj_incomp_NACA0012.cfg_dir = "disc_adj_incomp_euler/naca0012" discadj_incomp_NACA0012.cfg_file = "incomp_NACA0012_disc.cfg" discadj_incomp_NACA0012.test_iter = 20 - discadj_incomp_NACA0012.test_vals = [20.000000, -4.091640, -2.655563, 0.000000] + discadj_incomp_NACA0012.test_vals = [20.000000, -3.977753, -2.562520, 0.000000] test_list.append(discadj_incomp_NACA0012) ##################################### @@ -187,7 +187,7 @@ def main(): discadj_pitchingNACA0012.cfg_dir = "disc_adj_euler/naca0012_pitching" discadj_pitchingNACA0012.cfg_file = "inv_NACA0012_pitching.cfg" discadj_pitchingNACA0012.test_iter = 4 - discadj_pitchingNACA0012.test_vals = [-1.220333, -1.646832, -0.007539, 0.000013] + discadj_pitchingNACA0012.test_vals = [-1.222085, -1.652435, -0.004121, -0.000003] discadj_pitchingNACA0012.unsteady = True discadj_pitchingNACA0012.enabled_with_tsan = False test_list.append(discadj_pitchingNACA0012) diff --git a/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref b/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref index fce3c63c676..383fee127a9 100644 --- a/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref +++ b/TestCases/multiple_ffd/naca0012/of_grad_cd.dat.ref @@ -1,3 +1,3 @@ VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP" - 0 , 0.0790825 , 0.001 - 1 , -0.112974 , 0.001 + 0 , 0.0645833 , 0.001 + 1 , -0.118221 , 0.001 diff --git a/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref b/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref index 27fc6f627ae..48fe6fd25d9 100644 --- a/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref +++ b/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref @@ -1,3 +1,3 @@ VARIABLES="VARIABLE" , "DRAG" , "EFFICIENCY" , "FORCE_X" , "FORCE_Y" , "FORCE_Z" , "LIFT" , "MOMENT_X" , "MOMENT_Y" , "MOMENT_Z" , "SIDEFORCE" - 0 , 0.05922508178 , -0.6082054907 , 0.05650610011 , 0.1252552372 , 0.0 , 0.1239927558 , 0.0 , 0.0 , 0.00502894879 , 0.0 - 1 , -0.07750209216 , 8.684127966 , -0.08326907905 , 0.2634518171 , 0.0 , 0.2652056281 , 0.0 , 0.0 , -0.006643227386 , 0.0 + 0 , 0.06218459594 , 0.1261488787 , 0.05968184702 , 0.1153777146 , 0.0 , 0.1140483052 , 0.0 , 0.0 , 0.01868706266 , 0.0 + 1 , -0.1116496871 , 6.293726662 , -0.1174209952 , 0.2632773471 , 0.0 , 0.2657762197 , 0.0 , 0.0 , 0.04765040929 , 0.0 diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 0f78e0d8314..e19aa7a87bc 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -220,7 +220,7 @@ def main(): channel.cfg_dir = "euler/channel" channel.cfg_file = "inv_channel_RK.cfg" channel.test_iter = 20 - channel.test_vals = [-2.904401, 2.536139, 0.020924, 0.042340] + channel.test_vals = [-2.370377, 3.168384, 0.124246, 0.142878] test_list.append(channel) # NACA0012 @@ -228,7 +228,7 @@ def main(): naca0012.cfg_dir = "euler/naca0012" naca0012.cfg_file = "inv_NACA0012_Roe.cfg" naca0012.test_iter = 20 - naca0012.test_vals = [-4.607257, -4.138750, 0.327682, 0.022685] + naca0012.test_vals = [-5.024272, -4.487585, 0.332525, 0.023058] test_list.append(naca0012) # Supersonic wedge @@ -236,7 +236,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-1.387215, 4.281574, -0.245443, 0.043264] + wedge.test_vals = [-1.308216, 4.363572, -0.237802, 0.041962] test_list.append(wedge) # ONERA M6 Wing @@ -244,7 +244,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [-11.500303, -10.971884, 0.280800, 0.008623] + oneram6.test_vals = [-11.457168, -10.921423, 0.280800, 0.008623] oneram6.timeout = 3200 test_list.append(oneram6) @@ -253,7 +253,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-3.840915, 1.693979, 0.301144, 0.019489] + fixedCL_naca0012.test_vals = [-3.811341, 1.726020, 0.301248, 0.019495] test_list.append(fixedCL_naca0012) # Polar sweep of the inviscid NACA0012 @@ -262,7 +262,7 @@ def main(): polar_naca0012.cfg_file = "inv_NACA0012.cfg" polar_naca0012.polar = True polar_naca0012.test_iter = 10 - polar_naca0012.test_vals = [-1.096894, 4.372054, 0.002074, 0.031515] + polar_naca0012.test_vals = [-1.113609, 4.361071, 0.000427, 0.023153] polar_naca0012.test_vals_aarch64 = [-1.083394, 4.386134, 0.001588, 0.033513] polar_naca0012.command = TestCase.Command(exec = "compute_polar.py", param = "-i 11") # flaky test on arm64 @@ -318,7 +318,7 @@ def main(): flatplate_udobj.cfg_dir = "user_defined_functions" flatplate_udobj.cfg_file = "lam_flatplate.cfg" flatplate_udobj.test_iter = 20 - flatplate_udobj.test_vals = [-6.760101, -1.283906, -0.745653, 0.000587, -0.000038, 0.000977, -0.001015, 596.450000, 299.550000, 296.900000, 21.318000, 0.586640, 36.553000, 2.188800] + flatplate_udobj.test_vals = [-7.488320, -2.009145, 0.001084, 0.036247, 2.361500, -2.325300, 0.000000, 0.000000] test_list.append(flatplate_udobj) # Probe performance test (11 probes, ADT path) @@ -335,7 +335,7 @@ def main(): cylinder.cfg_dir = "navierstokes/cylinder" cylinder.cfg_file = "lam_cylinder.cfg" cylinder.test_iter = 25 - cylinder.test_vals = [-8.422012, -2.930518, -0.003394, 1.608420, 0.000000] + cylinder.test_vals = [-8.484215, -2.995693, -0.007163, 1.601414, 0.000000] test_list.append(cylinder) # Laminar cylinder (low Mach correction) @@ -343,7 +343,7 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.841606, -1.379533, -1.266782, 76.118168, 0.000000] + cylinder_lowmach.test_vals = [-6.490398, -1.028528, -0.149830, -141.498221, 0.000000] test_list.append(cylinder_lowmach) # 2D Poiseuille flow (body force driven with periodic inlet / outlet) @@ -360,7 +360,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.007496, -7.227093, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.011853, -7.672985, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.007498, -7.226926, -0.000000, 2.089953] test_list.append(poiseuille_profile) @@ -373,7 +373,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.004689, -5.265730, 0.809462, 0.062011, 0.000000] + rae2822_sa.test_vals = [-1.813393, -5.141451, 0.685984, 0.032019, 0.000000] test_list.append(rae2822_sa) # RAE2822 SST @@ -381,7 +381,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-0.510305, 5.386832, 0.813794, 0.062425, 0.000000] + rae2822_sst.test_vals = [-0.509972, 5.229105, 0.601704, 0.013220, 0.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -389,7 +389,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.644127, 5.386832, 0.813793, 0.062425] + rae2822_sst_sust.test_vals = [-2.473983, 5.229105, 0.601704, 0.013220] test_list.append(rae2822_sst_sust) # Flat plate @@ -397,7 +397,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.297192, -6.731227, -0.187632, 0.057700] + turb_flatplate.test_vals = [-3.966255, -6.830370, -0.198132, 0.022553] test_list.append(turb_flatplate) # Flat plate (compressible) with species inlet @@ -405,7 +405,7 @@ def main(): turb_flatplate_species.cfg_dir = "rans/flatplate" turb_flatplate_species.cfg_file = "turb_SA_flatplate_species.cfg" turb_flatplate_species.test_iter = 20 - turb_flatplate_species.test_vals = [-4.249462, -0.634904, -1.716285, 1.223210, -3.307930, 9.000000, -6.633666, 5.000000, -6.986787, 10.000000, -6.255670, 0.999903, 0.9999033] + turb_flatplate_species.test_vals = [-3.946196, -1.084969, -1.496278, 1.527583, -3.293946, 9.000000, -6.133018, 5.000000, -6.620066, 10.000000, -5.669045, 0.999913, 0.999913] test_list.append(turb_flatplate_species) # Flat plate SST compressibility correction Wilcox @@ -577,7 +577,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.574959, -5.057159, 0.830274, -0.009283, 0.078035] + turb_naca0012_sst_restart_mg.test_vals = [-6.604004, -5.057159, 0.830274, -0.008167, 0.077889] turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 test_list.append(turb_naca0012_sst_restart_mg) @@ -591,7 +591,7 @@ def main(): inc_euler_naca0012.cfg_dir = "incomp_euler/naca0012" inc_euler_naca0012.cfg_file = "incomp_NACA0012.cfg" inc_euler_naca0012.test_iter = 20 - inc_euler_naca0012.test_vals = [-7.154146, -6.507095, 0.532001, 0.008466] + inc_euler_naca0012.test_vals = [-7.156010, -6.592360, 0.531999, 0.008466] test_list.append(inc_euler_naca0012) # C-D nozzle with pressure inlet and mass flow outlet @@ -599,7 +599,7 @@ def main(): inc_nozzle.cfg_dir = "incomp_euler/nozzle" inc_nozzle.cfg_file = "inv_nozzle.cfg" inc_nozzle.test_iter = 20 - inc_nozzle.test_vals = [-6.407323, -5.715668, -0.003225, 0.126214] + inc_nozzle.test_vals = [-4.692695, -4.147623, 0.005125, 0.159938] test_list.append(inc_nozzle) # Laminar wall mounted cylinder, Euler walls, cylinder wall diagonally split @@ -619,7 +619,7 @@ def main(): inc_lam_cylinder.cfg_dir = "incomp_navierstokes/cylinder" inc_lam_cylinder.cfg_file = "incomp_cylinder.cfg" inc_lam_cylinder.test_iter = 10 - inc_lam_cylinder.test_vals = [-4.004072, -3.194881, -0.076553, 7.780048] + inc_lam_cylinder.test_vals = [-3.995539, -3.270342, -0.100583, 6.567311] test_list.append(inc_lam_cylinder) # Laminar sphere, Re=1. Last column: Cd=24/Re @@ -627,7 +627,7 @@ def main(): inc_lam_sphere.cfg_dir = "incomp_navierstokes/sphere" inc_lam_sphere.cfg_file = "sphere.cfg" inc_lam_sphere.test_iter = 5 - inc_lam_sphere.test_vals = [-8.342048, -9.328063, 0.121003, 25.782687] + inc_lam_sphere.test_vals = [-8.373846, -9.382816, 0.121003, 25.782686] test_list.append(inc_lam_sphere) # Buoyancy-driven cavity @@ -765,7 +765,7 @@ def main(): turbmod_sa_bsl_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_bsl_rae2822.cfg_file = "turb_SA_BSL_RAE2822.cfg" turbmod_sa_bsl_rae2822.test_iter = 20 - turbmod_sa_bsl_rae2822.test_vals = [-2.004689, 0.742307, 0.497308, -5.265730, 0.809462, 0.062011] + turbmod_sa_bsl_rae2822.test_vals = [-1.813393, 0.818366, 0.644951, -5.141451, 0.685984, 0.032019] test_list.append(turbmod_sa_bsl_rae2822) # SA Negative @@ -782,7 +782,7 @@ def main(): turbmod_sa_comp_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_comp_rae2822.cfg_file = "turb_SA_COMP_RAE2822.cfg" turbmod_sa_comp_rae2822.test_iter = 20 - turbmod_sa_comp_rae2822.test_vals = [-2.004688, 0.742305, 0.497309, -5.266066, 0.809466, 0.062026] + turbmod_sa_comp_rae2822.test_vals = [-1.813396, 0.818361, 0.644949, -5.141555, 0.685988, 0.032021] test_list.append(turbmod_sa_comp_rae2822) # SA Edwards @@ -790,7 +790,7 @@ def main(): turbmod_sa_edw_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_edw_rae2822.cfg_file = "turb_SA_EDW_RAE2822.cfg" turbmod_sa_edw_rae2822.test_iter = 20 - turbmod_sa_edw_rae2822.test_vals = [-2.004687, 0.742306, 0.497310, -5.290787, 0.809485, 0.062034] + turbmod_sa_edw_rae2822.test_vals = [-1.813370, 0.818386, 0.644973, -5.177642, 0.685990, 0.032017] test_list.append(turbmod_sa_edw_rae2822) # SA Compressibility and Edwards @@ -798,7 +798,7 @@ def main(): turbmod_sa_comp_edw_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_comp_edw_rae2822.cfg_file = "turb_SA_COMP_EDW_RAE2822.cfg" turbmod_sa_comp_edw_rae2822.test_iter = 20 - turbmod_sa_comp_edw_rae2822.test_vals = [-2.004686, 0.742307, 0.497311, -5.290776, 0.809487, 0.062044] + turbmod_sa_comp_edw_rae2822.test_vals = [-1.813375, 0.818380, 0.644969, -5.177681, 0.685994, 0.032020] test_list.append(turbmod_sa_comp_edw_rae2822) # SA QCR @@ -806,7 +806,7 @@ def main(): turbmod_sa_qcr_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_qcr_rae2822.cfg_file = "turb_SA_QCR_RAE2822.cfg" turbmod_sa_qcr_rae2822.test_iter = 20 - turbmod_sa_qcr_rae2822.test_vals = [-2.004794, 0.742354, 0.497315, -5.265910, 0.807835, 0.062022] + turbmod_sa_qcr_rae2822.test_vals = [-1.813960, 0.818091, 0.644509, -5.141963, 0.685849, 0.032055] test_list.append(turbmod_sa_qcr_rae2822) ############################ @@ -830,7 +830,7 @@ def main(): contadj_naca0012.cfg_dir = "cont_adj_euler/naca0012" contadj_naca0012.cfg_file = "inv_NACA0012.cfg" contadj_naca0012.test_iter = 5 - contadj_naca0012.test_vals = [-9.662585, -14.998832, -0.726250, 0.020280] + contadj_naca0012.test_vals = [-9.586529, -15.036640, -0.726250, 0.020280] contadj_naca0012.test_vals_aarch64 = [-9.662546, -14.998818, -0.726250, 0.020280] test_list.append(contadj_naca0012) @@ -839,7 +839,7 @@ def main(): contadj_oneram6.cfg_dir = "cont_adj_euler/oneram6" contadj_oneram6.cfg_file = "inv_ONERAM6.cfg" contadj_oneram6.test_iter = 10 - contadj_oneram6.test_vals = [-12.032190, -12.587083, -1.086100, 0.007556] + contadj_oneram6.test_vals = [-11.981564, -12.514942, -1.086100, 0.007556] test_list.append(contadj_oneram6) # Inviscid WEDGE: tests averaged outflow total pressure adjoint @@ -855,7 +855,7 @@ def main(): contadj_fixed_CL_naca0012.cfg_dir = "fixed_cl/naca0012" contadj_fixed_CL_naca0012.cfg_file = "inv_NACA0012_ContAdj.cfg" contadj_fixed_CL_naca0012.test_iter = 100 - contadj_fixed_CL_naca0012.test_vals = [0.748407, -4.810872, -0.520110, -0.000291] + contadj_fixed_CL_naca0012.test_vals = [0.837601, -4.799651, -0.371420, -0.000601] test_list.append(contadj_fixed_CL_naca0012) ################################### @@ -867,7 +867,7 @@ def main(): contadj_ns_cylinder.cfg_dir = "cont_adj_navierstokes/cylinder" contadj_ns_cylinder.cfg_file = "lam_cylinder.cfg" contadj_ns_cylinder.test_iter = 20 - contadj_ns_cylinder.test_vals = [-3.651430, -9.113079, 2.056700, -0.000000] + contadj_ns_cylinder.test_vals = [-3.600568, -9.041782, 2.056700, -0.000000] test_list.append(contadj_ns_cylinder) # Adjoint laminar naca0012 subsonic @@ -911,7 +911,7 @@ def main(): contadj_rans_rae2822.cfg_dir = "cont_adj_rans/rae2822" contadj_rans_rae2822.cfg_file = "turb_SA_RAE2822.cfg" contadj_rans_rae2822.test_iter = 20 - contadj_rans_rae2822.test_vals = [-5.372407, -10.874841, -0.212470, 0.005448] + contadj_rans_rae2822.test_vals = [-5.385804, -10.889762, -0.212470, 0.005448] test_list.append(contadj_rans_rae2822) ############################# @@ -999,7 +999,7 @@ def main(): cavity.cfg_dir = "moving_wall/cavity" cavity.cfg_file = "lam_cavity.cfg" cavity.test_iter = 25 - cavity.test_vals = [-5.610923, -0.146741, 1.115860, 1.490430] + cavity.test_vals = [-5.348408, 0.117949, -0.454913, -5.607752] test_list.append(cavity) # Spinning cylinder @@ -1007,7 +1007,7 @@ def main(): spinning_cylinder.cfg_dir = "moving_wall/spinning_cylinder" spinning_cylinder.cfg_file = "spinning_cylinder.cfg" spinning_cylinder.test_iter = 25 - spinning_cylinder.test_vals = [-7.806025, -2.364860, 1.685247, 1.518292] + spinning_cylinder.test_vals = [-7.710894, -2.258032, 1.684241, 1.582301] test_list.append(spinning_cylinder) ###################################### @@ -1028,7 +1028,7 @@ def main(): sine_gust.cfg_dir = "gust" sine_gust.cfg_file = "inv_gust_NACA0012.cfg" sine_gust.test_iter = 5 - sine_gust.test_vals = [-1.977498, 3.481817, -0.010773, -0.008068] + sine_gust.test_vals = [-1.977318, 3.481877, -0.000142, -0.007861] sine_gust.unsteady = True test_list.append(sine_gust) @@ -1037,7 +1037,7 @@ def main(): aeroelastic.cfg_dir = "aeroelastic" aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" aeroelastic.test_iter = 2 - aeroelastic.test_vals = [0.075023, 0.027483, -0.001643, -0.000126] + aeroelastic.test_vals = [0.069157, 0.027462, -0.001829, 0.000025] aeroelastic.unsteady = True test_list.append(aeroelastic) @@ -1077,7 +1077,7 @@ def main(): edge_VW.cfg_dir = "nicf/edge" edge_VW.cfg_file = "edge_VW.cfg" edge_VW.test_iter = 50 - edge_VW.test_vals = [-9.057409, -2.833203, -0.000009, 0.000000] + edge_VW.test_vals = [-2.770557, 3.430361, -0.000009, 0.000000] test_list.append(edge_VW) # Rarefaction shock wave edge_PPR @@ -1085,7 +1085,7 @@ def main(): edge_PPR.cfg_dir = "nicf/edge" edge_PPR.cfg_file = "edge_PPR.cfg" edge_PPR.test_iter = 50 - edge_PPR.test_vals = [-9.781896, -3.630892, -0.000034, 0.000000] + edge_PPR.test_vals = [-10.876684, -4.693540, -0.000034, 0.000000] test_list.append(edge_PPR) # Rarefaction Q1D nozzle, include CoolProp fluid model @@ -1434,7 +1434,7 @@ def main(): pywrapper_naca0012 = TestCase('pywrapper_naca0012') pywrapper_naca0012.cfg_dir = "euler/naca0012" pywrapper_naca0012.cfg_file = "inv_NACA0012_Roe.cfg" - pywrapper_naca0012.test_iter = 100 + pywrapper_naca0012.test_iter = 80 pywrapper_naca0012.test_vals = [-9.249009, -8.546597, 0.335769, 0.023275] pywrapper_naca0012.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f") test_list.append(pywrapper_naca0012) @@ -1506,7 +1506,7 @@ def main(): pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" pywrapper_unsteadyCHT.cfg_file = "unsteady_CHT_FlatPlate_Conf.cfg" pywrapper_unsteadyCHT.test_iter = 5 - pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.264219, -0.001386, 0.172978] + pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.263483, -0.007487, 0.100530] pywrapper_unsteadyCHT.command = TestCase.Command("mpirun -np 2", "python", "launch_unsteady_CHT_FlatPlate.py --parallel -f") pywrapper_unsteadyCHT.unsteady = True test_list.append(pywrapper_unsteadyCHT) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index d648bb852c7..76005d82f45 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -93,7 +93,7 @@ def main(): channel.cfg_dir = "euler/channel" channel.cfg_file = "inv_channel_RK.cfg" channel.test_iter = 10 - channel.test_vals = [-2.691660, 2.781413, -0.009401, 0.011862] + channel.test_vals = [-2.262915, 3.282305, 0.084026, 0.157104] test_list.append(channel) # NACA0012 @@ -101,7 +101,7 @@ def main(): naca0012.cfg_dir = "euler/naca0012" naca0012.cfg_file = "inv_NACA0012_Roe.cfg" naca0012.test_iter = 20 - naca0012.test_vals = [-4.766184, -4.287722, 0.326688, 0.022661] + naca0012.test_vals = [-5.063957, -4.497026, 0.333597, 0.023126] test_list.append(naca0012) # Supersonic wedge @@ -109,7 +109,7 @@ def main(): wedge.cfg_dir = "euler/wedge" wedge.cfg_file = "inv_wedge_HLLC.cfg" wedge.test_iter = 20 - wedge.test_vals = [-1.379426, 4.288828, -0.245341, 0.043244] + wedge.test_vals = [-1.315151, 4.355264, -0.239167, 0.042188] test_list.append(wedge) # ONERA M6 Wing @@ -117,7 +117,7 @@ def main(): oneram6.cfg_dir = "euler/oneram6" oneram6.cfg_file = "inv_ONERAM6.cfg" oneram6.test_iter = 10 - oneram6.test_vals = [-11.498143, -10.969216, 0.280800, 0.008623] + oneram6.test_vals = [-11.428953, -10.889491, 0.280800, 0.008623] oneram6.timeout = 9600 test_list.append(oneram6) @@ -126,7 +126,7 @@ def main(): fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" fixedCL_naca0012.cfg_file = "inv_NACA0012.cfg" fixedCL_naca0012.test_iter = 10 - fixedCL_naca0012.test_vals = [-3.837516, 1.700577, 0.301169, 0.019490] + fixedCL_naca0012.test_vals = [-3.751143, 1.788901, 0.301332, 0.019499] test_list.append(fixedCL_naca0012) # Polar sweep of the inviscid NACA0012 @@ -135,7 +135,7 @@ def main(): polar_naca0012.cfg_file = "inv_NACA0012.cfg" polar_naca0012.polar = True polar_naca0012.test_iter = 10 - polar_naca0012.test_vals = [-1.077848, 4.386916, -0.000333, 0.029678] + polar_naca0012.test_vals = [-1.138311, 4.325266, 0.002828, 0.019184] polar_naca0012.test_vals_aarch64 = [-1.063447, 4.401847, 0.000291, 0.031696] polar_naca0012.command = TestCase.Command(exec = "compute_polar.py", param = "-n 1 -i 11") # flaky test on arm64 @@ -166,7 +166,7 @@ def main(): flatplate.cfg_dir = "navierstokes/flatplate" flatplate.cfg_file = "lam_flatplate.cfg" flatplate.test_iter = 20 - flatplate.test_vals = [-5.097199, 0.382306, 0.001326, 0.027904, 2.361500, -2.333600, 0.000000, 0.000000] + flatplate.test_vals = [-5.110155, 0.370560, 0.001251, 0.025767, 2.361500, -2.335800, 0.000000, 0.000000] test_list.append(flatplate) # Laminar cylinder (steady) @@ -174,7 +174,7 @@ def main(): cylinder.cfg_dir = "navierstokes/cylinder" cylinder.cfg_file = "lam_cylinder.cfg" cylinder.test_iter = 25 - cylinder.test_vals = [-8.363897, -2.882485, -0.017777, 1.607978, 0.000000] + cylinder.test_vals = [-8.643645, -3.150468, -0.002769, 1.602063, 0.000000] test_list.append(cylinder) # Laminar cylinder (low Mach correction) @@ -182,7 +182,7 @@ def main(): cylinder_lowmach.cfg_dir = "navierstokes/cylinder" cylinder_lowmach.cfg_file = "cylinder_lowmach.cfg" cylinder_lowmach.test_iter = 25 - cylinder_lowmach.test_vals = [-6.830989, -1.368842, -0.143868, 73.962350, 0.000000] + cylinder_lowmach.test_vals = [-6.504996, -1.043200, -0.119206, -136.159598, 0.000000] test_list.append(cylinder_lowmach) # 2D Poiseuille flow (body force driven with periodic inlet / outlet) @@ -198,7 +198,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.009004, -7.262034, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.012568, -7.696994, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns test_list.append(poiseuille_profile) @@ -218,7 +218,7 @@ def main(): rae2822_sa.cfg_dir = "rans/rae2822" rae2822_sa.cfg_file = "turb_SA_RAE2822.cfg" rae2822_sa.test_iter = 20 - rae2822_sa.test_vals = [-2.020123, -5.269264, 0.807147, 0.060494, 0.000000] + rae2822_sa.test_vals = [-1.806841, -5.123806, 0.685516, 0.029584, 0.000000] test_list.append(rae2822_sa) # RAE2822 SST @@ -226,7 +226,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-0.510321, 5.387978, 0.812737, 0.061080, 0.000000] + rae2822_sst.test_vals = [-0.510040, 5.230415, 0.612550, 0.009921, 0.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -234,7 +234,7 @@ def main(): rae2822_sst_sust.cfg_dir = "rans/rae2822" rae2822_sst_sust.cfg_file = "turb_SST_SUST_RAE2822.cfg" rae2822_sst_sust.test_iter = 20 - rae2822_sst_sust.test_vals = [-2.645780, 5.387978, 0.812737, 0.061080] + rae2822_sst_sust.test_vals = [-2.457541, 5.230415, 0.612550, 0.009921] test_list.append(rae2822_sst_sust) # Flat plate @@ -242,7 +242,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-4.316128, -6.738720, -0.187461, 0.057469] + turb_flatplate.test_vals = [-4.166987, -6.847620, -0.199569, 0.021651] test_list.append(turb_flatplate) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SST @@ -359,7 +359,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 50 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.576830, -5.081440, 0.810957, -0.008626, 0.077900] + turb_naca0012_sst_restart_mg.test_vals = [-6.609250, -5.081439, 0.810958, -0.008778, 0.077715] turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 test_list.append(turb_naca0012_sst_restart_mg) @@ -380,7 +380,7 @@ def main(): inc_euler_naca0012.cfg_dir = "incomp_euler/naca0012" inc_euler_naca0012.cfg_file = "incomp_NACA0012.cfg" inc_euler_naca0012.test_iter = 20 - inc_euler_naca0012.test_vals = [-7.140809, -6.485990, 0.531993, 0.008466] + inc_euler_naca0012.test_vals = [-7.127273, -6.574250, 0.531996, 0.008466] test_list.append(inc_euler_naca0012) # C-D nozzle with pressure inlet and mass flow outlet @@ -407,7 +407,7 @@ def main(): inc_lam_cylinder.cfg_dir = "incomp_navierstokes/cylinder" inc_lam_cylinder.cfg_file = "incomp_cylinder.cfg" inc_lam_cylinder.test_iter = 10 - inc_lam_cylinder.test_vals = [-4.004277, -3.227956, 0.003852, 7.626578] + inc_lam_cylinder.test_vals = [-4.077094, -3.309090, 0.006632, 5.905895] test_list.append(inc_lam_cylinder) # Buoyancy-driven cavity @@ -586,7 +586,7 @@ def main(): contadj_naca0012.cfg_dir = "cont_adj_euler/naca0012" contadj_naca0012.cfg_file = "inv_NACA0012.cfg" contadj_naca0012.test_iter = 5 - contadj_naca0012.test_vals = [-9.748339, -15.067997, -0.726250, 0.020280] + contadj_naca0012.test_vals = [-9.584244, -15.043117, -0.726250, 0.020280] contadj_naca0012.tol = 0.001 test_list.append(contadj_naca0012) @@ -595,7 +595,7 @@ def main(): contadj_oneram6.cfg_dir = "cont_adj_euler/oneram6" contadj_oneram6.cfg_file = "inv_ONERAM6.cfg" contadj_oneram6.test_iter = 10 - contadj_oneram6.test_vals = [-12.034680, -12.592674, -1.086100, 0.007556] + contadj_oneram6.test_vals = [-11.970149, -12.503534, -1.086100, 0.007556] test_list.append(contadj_oneram6) # Inviscid WEDGE: tests averaged outflow total pressure adjoint @@ -611,7 +611,7 @@ def main(): contadj_fixedCL_naca0012.cfg_dir = "fixed_cl/naca0012" contadj_fixedCL_naca0012.cfg_file = "inv_NACA0012_ContAdj.cfg" contadj_fixedCL_naca0012.test_iter = 100 - contadj_fixedCL_naca0012.test_vals = [0.754936, -4.793625, -0.524550, -0.000227] + contadj_fixedCL_naca0012.test_vals = [0.888169, -4.738760, -0.372890, -0.000614] test_list.append(contadj_fixedCL_naca0012) ################################### @@ -630,7 +630,7 @@ def main(): contadj_ns_cylinder.cfg_dir = "cont_adj_navierstokes/cylinder" contadj_ns_cylinder.cfg_file = "lam_cylinder.cfg" contadj_ns_cylinder.test_iter = 20 - contadj_ns_cylinder.test_vals = [-3.665842, -9.132048, 2.056700, -0.000000] + contadj_ns_cylinder.test_vals = [-3.600653, -9.043319, 2.056700, -0.000000] test_list.append(contadj_ns_cylinder) # Adjoint laminar naca0012 subsonic @@ -674,7 +674,7 @@ def main(): contadj_rans_rae2822.cfg_dir = "cont_adj_rans/rae2822" contadj_rans_rae2822.cfg_file = "turb_SA_RAE2822.cfg" contadj_rans_rae2822.test_iter = 20 - contadj_rans_rae2822.test_vals = [-5.369688, -10.872211, -0.212470, 0.005448] + contadj_rans_rae2822.test_vals = [-5.387481, -10.891118, -0.212470, 0.005448] test_list.append(contadj_rans_rae2822) ############################# @@ -760,7 +760,7 @@ def main(): cavity.cfg_dir = "moving_wall/cavity" cavity.cfg_file = "lam_cavity.cfg" cavity.test_iter = 25 - cavity.test_vals = [-5.627870, -0.164404, 0.054706, 2.545833] + cavity.test_vals = [-5.437082, 0.028919, -0.157082, -1.62953] test_list.append(cavity) # Spinning cylinder @@ -768,7 +768,7 @@ def main(): spinning_cylinder.cfg_dir = "moving_wall/spinning_cylinder" spinning_cylinder.cfg_file = "spinning_cylinder.cfg" spinning_cylinder.test_iter = 25 - spinning_cylinder.test_vals = [-7.894513, -2.469083, 1.703823, 1.670412] + spinning_cylinder.test_vals = [-7.690469, -2.232131, 1.517075, 1.520828] test_list.append(spinning_cylinder) ###################################### @@ -789,7 +789,7 @@ def main(): sine_gust.cfg_dir = "gust" sine_gust.cfg_file = "inv_gust_NACA0012.cfg" sine_gust.test_iter = 5 - sine_gust.test_vals = [-1.977498, 3.481817, -0.010657, -0.007976] + sine_gust.test_vals = [-1.977319, 3.481878, 0.001450, -0.007578] sine_gust.unsteady = True test_list.append(sine_gust) @@ -798,7 +798,7 @@ def main(): aeroelastic.cfg_dir = "aeroelastic" aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" aeroelastic.test_iter = 2 - aeroelastic.test_vals = [0.074208, 0.027599, -0.001641, -0.000128] + aeroelastic.test_vals = [0.069270, 0.027429, -0.001828, 0.000019] aeroelastic.unsteady = True test_list.append(aeroelastic) @@ -825,7 +825,7 @@ def main(): unst_deforming_naca0012.cfg_dir = "disc_adj_euler/naca0012_pitching_def" unst_deforming_naca0012.cfg_file = "inv_NACA0012_pitching_deform.cfg" unst_deforming_naca0012.test_iter = 5 - unst_deforming_naca0012.test_vals = [-3.665152, -3.793306, -3.716483, -3.148336] + unst_deforming_naca0012.test_vals = [-3.658963, -3.761517, -3.698530, -3.138323] unst_deforming_naca0012.unsteady = True test_list.append(unst_deforming_naca0012) @@ -846,7 +846,7 @@ def main(): edge_VW.cfg_dir = "nicf/edge" edge_VW.cfg_file = "edge_VW.cfg" edge_VW.test_iter = 20 - edge_VW.test_vals = [-0.772250, 5.429879, -0.000470, 0.000000] + edge_VW.test_vals = [-0.557647, 5.644222, -0.006682, 0.000000] test_list.append(edge_VW) # Rarefaction shock wave edge_PPR @@ -854,7 +854,7 @@ def main(): edge_PPR.cfg_dir = "nicf/edge" edge_PPR.cfg_file = "edge_PPR.cfg" edge_PPR.test_iter = 20 - edge_PPR.test_vals = [-2.126694, 4.066051, -0.000013, 0.000000] + edge_PPR.test_vals = [-1.689532, 4.503673, 0.000353, 0.000000] test_list.append(edge_PPR) @@ -1092,7 +1092,7 @@ def main(): airfoilRBF.cfg_file = "config.cfg" airfoilRBF.test_iter = 1 - airfoilRBF.test_vals = [1.000000, 0.086004, -3.365759] + airfoilRBF.test_vals = [1.000000, 0.120839, -3.389864] airfoilRBF.tol = 0.0001 airfoilRBF.multizone = True test_list.append(airfoilRBF) @@ -1499,7 +1499,7 @@ def main(): opt_multiobj1surf_py.cfg_dir = "optimization_euler/multiobjective_wedge" opt_multiobj1surf_py.cfg_file = "inv_wedge_ROE_multiobj_1surf.cfg" opt_multiobj1surf_py.test_iter = 1 - opt_multiobj1surf_py.test_vals = [1.000000, 1.000000, 36.122920, 4.611743] + opt_multiobj1surf_py.test_vals = [1.000000, 1.000000, 36.117840, 3.758531] opt_multiobj1surf_py.command = TestCase.Command(exec = "shape_optimization.py", param = "-g CONTINUOUS_ADJOINT -f") opt_multiobj1surf_py.timeout = 1600 opt_multiobj1surf_py.tol = 0.00001 @@ -1512,7 +1512,7 @@ def main(): opt_2surf1obj_py.cfg_dir = "optimization_euler/multiobjective_wedge" opt_2surf1obj_py.cfg_file = "inv_wedge_ROE_2surf_1obj.cfg" opt_2surf1obj_py.test_iter = 1 - opt_2surf1obj_py.test_vals = [1.000000, 1.000000, 2.005099, 0.000383] + opt_2surf1obj_py.test_vals = [1.000000, 1.000000, 2.005079, 0.000312] opt_2surf1obj_py.command = TestCase.Command(exec = "shape_optimization.py", param = "-g CONTINUOUS_ADJOINT -f") opt_2surf1obj_py.timeout = 1600 opt_2surf1obj_py.tol = 0.00001 @@ -1529,7 +1529,7 @@ def main(): pywrapper_naca0012.cfg_dir = "euler/naca0012" pywrapper_naca0012.cfg_file = "inv_NACA0012_Roe.cfg" pywrapper_naca0012.test_iter = 20 - pywrapper_naca0012.test_vals = [-4.766184, -4.287722, 0.326688, 0.022661] + pywrapper_naca0012.test_vals = [-5.063957, -4.497026, 0.333597, 0.023126] pywrapper_naca0012.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_naca0012.timeout = 1600 pywrapper_naca0012.tol = 0.00001 @@ -1570,7 +1570,7 @@ def main(): pywrapper_aeroelastic.cfg_dir = "aeroelastic" pywrapper_aeroelastic.cfg_file = "aeroelastic_NACA64A010.cfg" pywrapper_aeroelastic.test_iter = 2 - pywrapper_aeroelastic.test_vals = [0.074208, 0.027599, -0.001641, -0.000128] + pywrapper_aeroelastic.test_vals = [0.069270, 0.027429, -0.001828, 0.000019] pywrapper_aeroelastic.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_aeroelastic.timeout = 1600 pywrapper_aeroelastic.tol = 0.00001 @@ -1599,7 +1599,7 @@ def main(): pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" pywrapper_unsteadyCHT.cfg_file = "unsteady_CHT_FlatPlate_Conf.cfg" pywrapper_unsteadyCHT.test_iter = 5 - pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.264215, 0.000771, 0.172980] + pywrapper_unsteadyCHT.test_vals = [-1.614167, 2.263544, -0.004139, 0.099061] pywrapper_unsteadyCHT.command = TestCase.Command(exec = "python", param = "launch_unsteady_CHT_FlatPlate.py -f") pywrapper_unsteadyCHT.timeout = 1600 pywrapper_unsteadyCHT.tol = 0.00001 @@ -1627,7 +1627,7 @@ def main(): pywrapper_custom_inlet.cfg_dir = "py_wrapper/custom_inlet" pywrapper_custom_inlet.cfg_file = "lam_flatplate.cfg" pywrapper_custom_inlet.test_iter = 20 - pywrapper_custom_inlet.test_vals = [-4.123206, -1.543215, -3.735006, 1.339481, -0.793478, 0.161210, -0.007009, 0.513560, -0.520570] + pywrapper_custom_inlet.test_vals = [-4.122866, -1.542850, -3.664440, 1.339854, -0.793567, 0.159301, -0.007466, 0.514290, -0.521750] pywrapper_custom_inlet.command = TestCase.Command(exec = "python", param = "run.py") pywrapper_custom_inlet.timeout = 1600 pywrapper_custom_inlet.tol = 0.0001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 50381c04724..2d811a77c34 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -76,7 +76,7 @@ def main(): discadj_arina2k.cfg_dir = "disc_adj_euler/arina2k" discadj_arina2k.cfg_file = "Arina2KRS.cfg" discadj_arina2k.test_iter = 20 - discadj_arina2k.test_vals = [-3.254490, -3.495569, 0.052370, 0.000000] + discadj_arina2k.test_vals = [-3.229386, -3.466956, 0.043983, 0.000000] test_list.append(discadj_arina2k) ####################################################### @@ -108,7 +108,7 @@ def main(): discadj_incomp_NACA0012.cfg_dir = "disc_adj_incomp_euler/naca0012" discadj_incomp_NACA0012.cfg_file = "incomp_NACA0012_disc.cfg" discadj_incomp_NACA0012.test_iter = 20 - discadj_incomp_NACA0012.test_vals = [20.000000, -4.091644, -2.655563, 0.000000] + discadj_incomp_NACA0012.test_vals = [20.000000, -3.977751, -2.562520, 0.000000] test_list.append(discadj_incomp_NACA0012) ##################################### @@ -158,7 +158,7 @@ def main(): discadj_pitchingNACA0012.cfg_dir = "disc_adj_euler/naca0012_pitching" discadj_pitchingNACA0012.cfg_file = "inv_NACA0012_pitching.cfg" discadj_pitchingNACA0012.test_iter = 4 - discadj_pitchingNACA0012.test_vals = [-1.219518, -1.646199, -0.007607, 0.000013] + discadj_pitchingNACA0012.test_vals = [1.221535, -1.651330, -0.004201, -0.000004] discadj_pitchingNACA0012.unsteady = True test_list.append(discadj_pitchingNACA0012) @@ -167,7 +167,7 @@ def main(): unst_deforming_naca0012.cfg_dir = "disc_adj_euler/naca0012_pitching_def" unst_deforming_naca0012.cfg_file = "inv_NACA0012_pitching_deform_ad.cfg" unst_deforming_naca0012.test_iter = 4 - unst_deforming_naca0012.test_vals = [-1.960419, -1.844186, 2970.700000, 0.000004] + unst_deforming_naca0012.test_vals = [-1.965271, -1.847913, 2944.900000, 0.000000] unst_deforming_naca0012.unsteady = True test_list.append(unst_deforming_naca0012) diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index b02e28e7063..283dc8d606d 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -175,7 +175,7 @@ def main(): tutorial_inv_bump.cfg_dir = "../Tutorials/compressible_flow/Inviscid_Bump" tutorial_inv_bump.cfg_file = "inv_channel.cfg" tutorial_inv_bump.test_iter = 0 - tutorial_inv_bump.test_vals = [-1.437425, 4.075857, 0.035200, 0.019230] + tutorial_inv_bump.test_vals = [-1.437425, 4.075857, 0.062355, 0.000402] test_list.append(tutorial_inv_bump) # Inviscid Wedge @@ -183,7 +183,7 @@ def main(): tutorial_inv_wedge.cfg_dir = "../Tutorials/compressible_flow/Inviscid_Wedge" tutorial_inv_wedge.cfg_file = "inv_wedge_HLLC.cfg" tutorial_inv_wedge.test_iter = 0 - tutorial_inv_wedge.test_vals = [-0.481460, 5.253008, -0.281099, 0.049535] + tutorial_inv_wedge.test_vals = [-0.481460, 5.253008, -0.290711, 0.051397] tutorial_inv_wedge.no_restart = True test_list.append(tutorial_inv_wedge) @@ -192,7 +192,7 @@ def main(): tutorial_inv_onera.cfg_dir = "../Tutorials/compressible_flow/Inviscid_ONERAM6" tutorial_inv_onera.cfg_file = "inv_ONERAM6.cfg" tutorial_inv_onera.test_iter = 0 - tutorial_inv_onera.test_vals = [-5.204928, -4.597762, 0.294332, 0.115223] + tutorial_inv_onera.test_vals = [-5.204928, -4.597762, 0.261695, 0.083515] tutorial_inv_onera.no_restart = True test_list.append(tutorial_inv_onera) @@ -201,7 +201,7 @@ def main(): tutorial_lam_cylinder.cfg_dir = "../Tutorials/compressible_flow/Laminar_Cylinder" tutorial_lam_cylinder.cfg_file = "lam_cylinder.cfg" tutorial_lam_cylinder.test_iter = 0 - tutorial_lam_cylinder.test_vals = [-6.162141, -0.699617, 0.126007, 69.619462] + tutorial_lam_cylinder.test_vals = [-6.162141, -0.699617, 0.014753, 31.991973] tutorial_lam_cylinder.no_restart = True test_list.append(tutorial_lam_cylinder) @@ -323,7 +323,7 @@ def main(): tutorial_design_inv_naca0012.cfg_dir = "../Tutorials/design/Inviscid_2D_Unconstrained_NACA0012" tutorial_design_inv_naca0012.cfg_file = "inv_NACA0012_basic.cfg" tutorial_design_inv_naca0012.test_iter = 0 - tutorial_design_inv_naca0012.test_vals = [-3.585391, -2.989014, 0.169337, 0.235131] + tutorial_design_inv_naca0012.test_vals = [-3.585391, -2.989014, 0.166575, 0.241919] tutorial_design_inv_naca0012.no_restart = True test_list.append(tutorial_design_inv_naca0012)