diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 7c07dd8a0cc..b9efe9edff4 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,
@@ -2914,7 +2914,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..942895279b9 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 dedicated 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..f6c99d4a455 100644
--- a/Common/src/geometry/CMultiGridGeometry.cpp
+++ b/Common/src/geometry/CMultiGridGeometry.cpp
@@ -33,14 +33,24 @@
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) {
@@ -67,15 +77,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,42 +110,93 @@ 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 one of them is 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;
+ 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]]++;
+ }
+ }
}
}
@@ -132,7 +204,8 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
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 ---*/
@@ -149,34 +222,45 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
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();
+ /*--- 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);
+ 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 ---*/
+ /*--- 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 ---*/
+ 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 ---*/
+ if (SetBoundAgglomeration(CVPoint, marker_seed, fine_grid, config)) {
+ /*--- We set the value of the parent ---*/
- fine_grid->nodes->SetParent_CV(CVPoint, Index_CoarseCV);
+ 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. ---*/
+ /*--- 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);
+ if (fine_grid->nodes->GetAgglomerate_Indirect(CVPoint))
+ nodes->SetAgglomerate_Indirect(Index_CoarseCV, true);
- /*--- We set the value of the child ---*/
+ /*--- We set the value of the child ---*/
- nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint);
- nChildren++;
+ nodes->SetChildren_CV(Index_CoarseCV, nChildren, CVPoint);
+ nChildren++;
+ /*--- Apply maxAgglomSize limit for 3D internal boundary face nodes. ---*/
+ if (nChildren == maxAgglomSize) break;
+ }
}
}
}
@@ -223,7 +307,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())) {
@@ -271,6 +355,7 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
MGQueue_InnerCV.Update(CVPoint, fine_grid);
}
+ if (nChildren == maxAgglomSize) break;
}
/*--- Identify the indirect neighbors ---*/
@@ -282,6 +367,8 @@ 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))) {
@@ -343,13 +430,54 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
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));
+
+ /*--- 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);
- auto nChildren = nodes->GetnChildren_CV(iCoarsePoint_Complete);
+ /*--- 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);
@@ -369,6 +497,16 @@ CMultiGridGeometry::CMultiGridGeometry(CGeometry* fine_grid, CConfig* config, un
nodes->ResetPoints();
#ifdef HAVE_MPI
+ /*--- Reset halo point parents before MPI agglomeration.
+ 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. ---*/
+
+ for (auto iPoint = fine_grid->GetnPointDomain(); iPoint < fine_grid->GetnPoint(); iPoint++) {
+ fine_grid->nodes->SetParent_CV(iPoint, std::numeric_limits::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 +514,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 +562,88 @@ 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(), std::numeric_limits::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 != std::numeric_limits::max()) continue;
+
+ /*--- Skip if received parent is invalid (sending rank didn't agglomerate this point) ---*/
+ if (Parent_Remote[iVertex] == std::numeric_limits::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] = std::numeric_limits::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;
}
}
+
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 nSkipped = 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. ---*/
+ /*--- Skip if this halo point was already agglomerated (first-wins policy) ---*/
+ auto existing_parent = fine_grid->nodes->GetParent_CV(iPoint_Fine);
+ if (existing_parent != std::numeric_limits::max()) {
+ nSkipped++;
+ continue; // First-wins: keep existing assignment
+ }
+
+ /*--- Skip if parent mapping is invalid (sender didn't agglomerate) ---*/
+ const auto iPoint_Coarse = Parent_Local[iVertex];
+ if (iPoint_Coarse == std::numeric_limits::max()) {
+ nSkipped++;
+ continue;
+ }
+
+ /*--- 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++;
}
}
}
+
#endif // HAVE_MPI
/*--- Update the number of points after the MPI agglomeration ---*/
@@ -463,25 +651,27 @@ 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();
+ 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);
- const su2double CFL = factor * config->GetCFL(iMesh - 1) / Coeff;
+ /*--- Note: CFL at the coarse levels have a large impact on convergence,
+ this should be rewritten to use adaptive CFL. ---*/
+ const su2double Coeff = 1.5;
+ const su2double CFL = config->GetCFL(iMesh - 1) / Coeff;
config->SetCFL(iMesh, CFL);
}
const su2double ratio = su2double(Global_nPointFine) / su2double(Global_nPointCoarse);
+ cout << "********** ratio = " << ratio << endl;
+ // lower value leads to more levels being accepted.
- if (((nDim == 2) && (ratio < 2.5)) || ((nDim == 3) && (ratio < 2.5))) {
+ 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 +693,65 @@ 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;
+ }
+
+ return (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 +761,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 +775,54 @@ 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;
+
+ // we also allow agglomeration if the seed has 2 markers, one of them is the same as the candidate, and the
+ // other is a send-receive marker.
+ if ((marker_seed.size() == 2) && ((copy_marker[0] == marker_seed[0]) || (copy_marker[0] == marker_seed[1]))) {
+ // check that the other marker is a send-receive marker
+ unsigned short other_marker = (copy_marker[0] == marker_seed[0]) ? marker_seed[1] : marker_seed[0];
+ if (config->GetMarker_All_KindBC(other_marker) == SEND_RECEIVE) {
+ agglomerate_CV = true;
}
}
}
/*--- 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 ---*/
-
- 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))) {
+ /*--- 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 (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 +842,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 +873,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 +904,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 +922,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 +944,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 +979,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 +995,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 +1019,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 +1035,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 +1062,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 +1076,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 +1103,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 +1195,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 +1234,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..fbc568bd7fc 100644
--- a/SU2_CFD/include/integration/CMultiGridIntegration.hpp
+++ b/SU2_CFD/include/integration/CMultiGridIntegration.hpp
@@ -164,8 +164,6 @@ class CMultiGridIntegration final : public CIntegration {
* \param[in] geo_fine - Geometrical definition of the fine grid.
* \param[in] geo_coarse - Geometrical definition of the coarse grid.
* \param[in] config - Definition of the particular problem.
- * \param[in] iMesh - Index of the mesh in multigrid computations.
- * \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);
@@ -181,4 +179,28 @@ class CMultiGridIntegration final : public CIntegration {
void Adjoint_Setup(CGeometry ****geometry, CSolver *****solver_container, CConfig **config,
unsigned short RunTime_EqSystem, unsigned long Iteration, unsigned short iZone);
+ /*!
+ * \brief Compute adaptive CFL for multigrid coarse levels.
+ * \param[in] config - Problem configuration.
+ * \param[in] solver_coarse - Coarse grid solver.
+ * \param[in] geometry_coarse - Coarse grid geometry.
+ * \param[in] iMesh - Current multigrid level.
+ * \param[in] CFL_fine - Fine grid CFL value (passive).
+ * \param[in] CFL_coarse_current - Current coarse grid CFL value (passive).
+ * \return New CFL value for the coarse grid.
+ */
+ passivedouble computeMultigridCFL(CConfig* config, CSolver* solver_coarse, CGeometry* geometry_coarse,
+ unsigned short iMesh, passivedouble CFL_fine, passivedouble CFL_coarse_current);
+
+ /*--- CFL adaptation state variables ---*/
+ static constexpr int MAX_MG_LEVELS = 10;
+ su2double current_avg[MAX_MG_LEVELS] = {};
+ su2double prev_avg[MAX_MG_LEVELS] = {};
+ su2double last_res[MAX_MG_LEVELS] = {};
+ bool last_was_increase[MAX_MG_LEVELS] = {};
+ int oscillation_count[MAX_MG_LEVELS] = {};
+ unsigned long last_check_iter[MAX_MG_LEVELS] = {};
+ unsigned long last_update_iter[MAX_MG_LEVELS] = {};
+ unsigned long last_reset_iter = std::numeric_limits::max();
+
};
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp b/SU2_CFD/include/solvers/CFVMFlowSolverBase.hpp
index f42ef9a731e..55f9ff3c32a 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..697bef04d7b 100644
--- a/SU2_CFD/src/integration/CMultiGridIntegration.cpp
+++ b/SU2_CFD/src/integration/CMultiGridIntegration.cpp
@@ -28,6 +28,211 @@
#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) continue;
+
+ 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
+
+passivedouble CMultiGridIntegration::computeMultigridCFL(CConfig* config, CSolver* solver_coarse, CGeometry* geometry_coarse,
+ unsigned short iMesh, passivedouble CFL_fine, passivedouble CFL_coarse_current) {
+ passivedouble 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;
+
+ /*--- Only master thread performs CFL adaptation to ensure deterministic behavior ---*/
+ /*--- All adaptive CFL state and computation must be done by a single thread ---*/
+ passivedouble 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 ---*/
+ 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;
+
+ /*--- 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());
+ }
+
+ /*--- 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 ---*/
+ passivedouble new_coeff = current_coeff;
+ const passivedouble 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 (wrapped to prevent tape recording) ---*/
+ const bool wasActive = AD::BeginPassive();
+ config->SetCFL(iMesh+1, CFL_coarse_new);
+ AD::EndPassive(wasActive);
+ }
+ END_SU2_OMP_MASTER
+ /*--- Implicit barrier at end of master region ensures all threads see updated CFL ---*/
+
+ return CFL_coarse_new;
+}
CMultiGridIntegration::CMultiGridIntegration() : CIntegration() { }
@@ -64,8 +269,6 @@ 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. ---*/
-
SU2_OMP_PARALLEL_(if(solver_container[iZone][iInst][MESH_0][Solver_Position]->GetHasHybridParallel()))
{
@@ -170,6 +373,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++) {
@@ -201,8 +407,11 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
Space_Integration(geometry_fine, solver_container_fine, numerics_fine, config, iMesh, iRKStep, RunTime_EqSystem);
- /*--- Time integration, update solution using the old solution plus the solution increment ---*/
-
+ /*--- Time integration, update solution using the old solution plus the solution increment.
+ Must be called from ALL threads cooperatively: CSysSolve::Solve uses internal
+ SU2_OMP_BARRIER / SU2_OMP_FOR constructs that require participation of all threads.
+ Wrapping in SU2_OMP_MASTER would park non-master threads at the outer barrier,
+ starving the inner barriers and causing data races / deadlocks. ---*/
Time_Integration(geometry_fine, solver_container_fine, config, iRKStep, RunTime_EqSystem);
/*--- Send-Receive boundary conditions, and postprocessing ---*/
@@ -257,6 +466,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 +475,53 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
if (iMesh == config->GetnMGLevels()-2)
nextRecurseParam = 0;
+ /*--- Synchronize prior to recursive calls ---*/
+ 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);
+ /*--- Synchronize before reading CFL to avoid race with concurrent writes from computeMultigridCFL ---*/
+ SU2_OMP_BARRIER
+
+ /*--- Read current CFL values in master thread to ensure thread-safe access.
+ Extract passive values to avoid AD tape recording in parallel region. ---*/
+ // passivedouble CFL_fine_passive = 0.0;
+ // passivedouble CFL_coarse_current_passive = 0.0;
+
+ // SU2_OMP_MASTER
+ // {
+ // CFL_fine_passive = SU2_TYPE::GetValue(config->GetCFL(iMesh));
+ // CFL_coarse_current_passive = SU2_TYPE::GetValue(config->GetCFL(iMesh+1));
+ // }
+ // END_SU2_OMP_MASTER
+ /*--- Implicit barrier here ensures all threads see the CFL values before proceeding ---*/
+
+ /*--- Compute adaptive CFL for coarse grid using passive values (no tape recording) ---*/
+
+ //nijso: just put a fixed value to see if the omp problem is really here...
+ passivedouble CFL_coarse_new = 1.0; //computeMultigridCFL(config, solver_coarse, geometry_coarse, iMesh, CFL_fine_passive, CFL_coarse_current_passive);
+
+ /*--- Explicit barrier to ensure all threads see the updated CFL from computeMultigridCFL
+ before proceeding. This prevents data races where threads at different recursion depths
+ might read stale CFL values during nested MultiGrid_Cycle calls. ---*/
+ SU2_OMP_BARRIER
+
+ /*--- 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
+
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 +531,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);
@@ -295,11 +548,13 @@ void CMultiGridIntegration::MultiGrid_Cycle(CGeometry ****geometry,
Space_Integration(geometry_fine, solver_container_fine, numerics_fine, config, iMesh, iRKStep, RunTime_EqSystem);
+ /*--- Must be called cooperatively from all threads (see presmooth comment). ---*/
Time_Integration(geometry_fine, solver_container_fine, config, iRKStep, RunTime_EqSystem);
solver_fine->Postprocessing(geometry_fine, solver_container_fine, config, iMesh);
}
+
}
}
@@ -307,9 +562,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 +569,47 @@ 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;
+ /*--- 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 ---*/
@@ -370,9 +628,9 @@ void CMultiGridIntegration::GetProlongated_Correction(unsigned short RunTime_EqS
sol_coarse->CompleteComms(geo_coarse, config, MPI_QUANTITIES::SOLUTION_OLD);
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 +645,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 +657,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 +680,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 +715,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 +744,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 +758,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 +766,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 +783,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 +795,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 +806,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();
@@ -606,6 +859,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 +871,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 f5bea2d55b8..f70a9df21ea 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.510301, 6.055876, 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.589038, 6.055876, 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
@@ -1442,7 +1442,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)
@@ -1514,7 +1514,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.260964, -0.001386, 0.172980]
+ 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 55fdf623942..a445aa7006d 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.510318, 6.056300, 0.812738, 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.591046, 6.056300, 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.260960, 0.000771, 0.172982]
+ 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 424c2a560c3..91962d2b9a8 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 b2dfe5ffd92..f117b092c2a 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.060602, -0.001097]
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.290816, 0.051445]
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.124663, 31.721714]
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.165195, 0.238368]
tutorial_design_inv_naca0012.no_restart = True
test_list.append(tutorial_design_inv_naca0012)