diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 992a0d6e30ad..960f52ac1d7f 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -235,7 +235,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107 with: @@ -284,7 +284,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tapetests" + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} -a "--tapetests" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:261011-1107 with: @@ -334,7 +334,7 @@ jobs: PMIX_MCA_gds: hash with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--tsan" + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} -a "--tsan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-tsan:261011-1107 with: @@ -379,7 +379,7 @@ jobs: uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107 with: # -t -c - args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} -a "--asan" + args: -b ${{github.ref}} -t develop -c fix_periodic_axis -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107 with: diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index ca8a3a26aabd..1c4f29dbd969 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -344,8 +344,14 @@ class CGeometry { a particular vertex to be sent in periodic comms. */ *Local_Marker_PeriodicRecv{nullptr}; /*!< \brief Data structure holding the local index of the periodic marker for a particular vertex to be received in periodic comms. */ - su2double* bufD_PeriodicRecv{nullptr}; /*!< \brief Data structure for su2double periodic receive. */ - su2double* bufD_PeriodicSend{nullptr}; /*!< \brief Data structure for su2double periodic send. */ + vector Local_Copy_PeriodicSend; /*!< \brief For points on a rotation axis, which are their own periodic + match, the number of times the rotation is applied to the data + sent (one send per other copy of the control volume). Zero for + all other points. */ + vector Local_Copy_PeriodicRecv; /*!< \brief Same as Local_Copy_PeriodicSend, for the data received. */ + bool PeriodicAxisPoints{false}; /*!< \brief Whether this rank has periodic points on a rotation axis. */ + su2double* bufD_PeriodicRecv{nullptr}; /*!< \brief Data structure for su2double periodic receive. */ + su2double* bufD_PeriodicSend{nullptr}; /*!< \brief Data structure for su2double periodic send. */ unsigned short* bufS_PeriodicRecv{nullptr}; /*!< \brief Data structure for unsigned long periodic receive. */ unsigned short* bufS_PeriodicSend{nullptr}; /*!< \brief Data structure for unsigned long periodic send. */ SU2_MPI::Request* req_PeriodicSend{nullptr}; /*!< \brief Data structure for periodic send requests. */ @@ -519,6 +525,12 @@ class CGeometry { */ void AllocatePeriodicComms(unsigned short val_countPerPeriodicPoint); + /*! \brief Sum partial scalar or vector geometry data over periodic copies (called by the master). */ + void SumPeriodicGeometry(const CConfig* config, su2activematrix& values, int vectorIndex); + + /*! \brief Refresh donor volumes after a mesh update. */ + void UpdatePeriodicVolumes(const CConfig* config); + /*! * \brief Routine to launch non-blocking recvs only for all periodic communication with neighboring partitions. * \note This routine is called by any class that has loaded data into the generic communication buffers. diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp index acafc2b7c025..9ed01302a3e9 100644 --- a/Common/src/geometry/CGeometry.cpp +++ b/Common/src/geometry/CGeometry.cpp @@ -1194,6 +1194,56 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { using PointMarkerSet = std::unordered_set; std::vector Points_Send_All(size, PointMarkerSet(0, pairHash)); + /*--- Points on a rotation axis are their own periodic match. The other copies of their control + volume are obtained by applying the rotation 1, 2, ... N-1 times (N = 360 deg / angle), instead + of once in each direction. Therefore, the vertex on the first marker of the pair sends N-1 times + to itself, and the vertex on the second marker does not send. ---*/ + + const unsigned long nMarker_All = config->GetnMarker_All(); + PeriodicAxisPoints = false; + + /*--- Number of copies of the domain around the rotation axis of a periodic marker (N). ---*/ + + auto nCopiesAroundAxis = [&](unsigned short jMarker) { + /*--- Tolerance, in radians, to accept the rotation angle as an integer fraction of 360 degrees. ---*/ + constexpr passivedouble angleTol = 1e-4; + + /*--- The trace of a rotation matrix is 1 + 2 cos(angle). ---*/ + const auto markerTag = config->GetMarker_All_TagBound(jMarker); + const su2double* angles = config->GetPeriodicRotAngles(markerTag); + su2double rotMatrix[MAXNDIM][MAXNDIM]; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotMatrix); + su2double cosAngle = 0.5 * (rotMatrix[0][0] + rotMatrix[1][1] + rotMatrix[2][2] - 1.0); + if (nDim == 2) cosAngle = cos(angles[2]); + cosAngle = min(cosAngle, 1.0); + cosAngle = max(cosAngle, -1.0); + const su2double angle = acos(cosAngle); + + const int nCopy = (angle > angleTol) ? SU2_TYPE::Int(2 * PI_NUMBER / angle + 0.5) : 0; + + if (fabs(nCopy * angle - 2 * PI_NUMBER) > angleTol) { + SU2_MPI::Error( + "Periodic points on the rotation axis require an angle that divides 360 degrees (marker " + markerTag + ").", + CURRENT_FUNCTION); + } + return nCopy; + }; + + /*--- Number of sends of a vertex (1 for the points that are not on a rotation axis), + and whether it is on a rotation axis. ---*/ + + auto nSendOfVertex = [&](unsigned short jMarker, unsigned long jVertex, bool& onAxis) -> unsigned long { + const auto* vert = geometry->vertex[jMarker][jVertex]; + onAxis = (static_cast(vert->GetDonorProcessor()) == rank) && + (static_cast(vert->GetDonorPoint()) == vert->GetNode()); + if (!onAxis) return 1; + if (config->GetMarker_All_PerBound(jMarker) > config->GetnMarker_Periodic() / 2) return 0; + + PeriodicAxisPoints = true; + return nCopiesAroundAxis(jMarker) - 1; + }; + bool onAxis = false; + /*--- Loop through all of our periodic markers and track our sends with each rank. ---*/ @@ -1213,8 +1263,11 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { iRank = static_cast(geometry->vertex[iMarker][iVertex]->GetDonorProcessor()); - /*--- Store the (point, marker) pair in the set for the destination rank. ---*/ - Points_Send_All[iRank].insert(std::make_pair(iPoint, static_cast(iMarker))); + /*--- Store the (point, marker) pair in the set for the destination rank, + once for each time the vertex sends (see above). ---*/ + const auto nCopy = nSendOfVertex(iMarker, iVertex, onAxis); + for (auto iCopy = 0ul; iCopy < nCopy; iCopy++) + Points_Send_All[iRank].insert(std::make_pair(iPoint, iMarker + iCopy * nMarker_All)); } } } @@ -1305,6 +1358,9 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { Local_Marker_PeriodicRecv = new unsigned long[nPoint_PeriodicRecv[nPeriodicRecv]]; for (iRecv = 0; iRecv < nPoint_PeriodicRecv[nPeriodicRecv]; iRecv++) Local_Marker_PeriodicRecv[iRecv] = 0; + Local_Copy_PeriodicSend.assign(nPoint_PeriodicSend[nPeriodicSend], 0); + Local_Copy_PeriodicRecv.assign(nPoint_PeriodicRecv[nPeriodicRecv], 0); + /*--- We allocate the buffers for communicating values in a later step once we know the maximum packet size that we need to communicate. This memory is deallocated and reallocated automatically in the case that @@ -1325,11 +1381,12 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { req_PeriodicRecv = new SU2_MPI::Request[nPeriodicRecv]; } - /*--- Allocate arrays for sending the periodic point index and marker - index to the recv rank so that it can store the local values. Therefore, - the recv rank can quickly loop through the buffers to unpack the data. ---*/ + /*--- Allocate arrays for sending the periodic point index, the marker index, + and the copy index (for points on a rotation axis) to the recv rank so that + it can store the local values. Therefore, the recv rank can quickly loop + through the buffers to unpack the data. ---*/ - unsigned short nPackets = 2; + const unsigned short nPackets = 3; auto* idSend = new unsigned long[nPoint_PeriodicSend[nPeriodicSend] * nPackets]; for (iSend = 0; iSend < nPoint_PeriodicSend[nPeriodicSend] * nPackets; iSend++) idSend[iSend] = 0; @@ -1366,17 +1423,24 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { if (iRank == destRank) { /*--- Check if we have already added this (point, marker) pair for this destination rank. Use the result if insert(), which is - a pair whose second element is success. ---*/ - const auto pointMarkerPair = std::make_pair(iPoint, static_cast(iMarker)); - const auto insertResult = Points_Send_All[destRank].insert(pointMarkerPair); - if (insertResult.second) { - Local_Point_PeriodicSend[ii] = iPoint; - Local_Marker_PeriodicSend[ii] = static_cast(iMarker); - jj = ii * nPackets; - idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); - jj++; - idSend[jj] = static_cast(iPeriodic); - ii++; + a pair whose second element is success. Points on a rotation + axis send more than once (see above). ---*/ + const auto nCopy = nSendOfVertex(iMarker, iVertex, onAxis); + for (auto iCopy = 0ul; iCopy < nCopy; iCopy++) { + const auto pointMarkerPair = std::make_pair(iPoint, iMarker + iCopy * nMarker_All); + const auto insertResult = Points_Send_All[destRank].insert(pointMarkerPair); + if (insertResult.second) { + Local_Point_PeriodicSend[ii] = iPoint; + Local_Marker_PeriodicSend[ii] = static_cast(iMarker); + Local_Copy_PeriodicSend[ii] = onAxis ? iCopy + 1 : 0; + jj = ii * nPackets; + idSend[jj] = geometry->vertex[iMarker][iVertex]->GetDonorPoint(); + jj++; + idSend[jj] = static_cast(iPeriodic); + jj++; + idSend[jj] = Local_Copy_PeriodicSend[ii]; + ii++; + } } } } @@ -1485,10 +1549,83 @@ void CGeometry::PreprocessPeriodicComms(CGeometry* geometry, CConfig* config) { ii++; Local_Marker_PeriodicRecv[iRecv] = idRecv[ii]; ii++; + Local_Copy_PeriodicRecv[iRecv] = idRecv[ii]; + ii++; } delete[] idSend; delete[] idRecv; + ComputeModifiedSymmetryNormals(config); +} + +void CGeometry::SumPeriodicGeometry(const CConfig* config, su2activematrix& values, int vectorIndex) { + if (!nPeriodicSend && !nPeriodicRecv) return; + const auto count = values.cols(); + AllocatePeriodicComms(count); + for (auto iPair = 1u; iPair <= config->GetnMarker_Periodic() / 2; ++iPair) { + PostPeriodicRecvs(this, config, COMM_TYPE::DOUBLE, count); + for (auto iMessage = 0; iMessage < nPeriodicSend; ++iMessage) { + for (auto iSend = nPoint_PeriodicSend[iMessage]; iSend < nPoint_PeriodicSend[iMessage + 1]; ++iSend) { + const auto iPoint = Local_Point_PeriodicSend[iSend]; + auto* buffer = bufD_PeriodicSend + count * iSend; + for (auto iVar = 0u; iVar < count; ++iVar) buffer[iVar] = values(iPoint, iVar); + if (vectorIndex >= 0) { + const auto* angles = + config->GetPeriodicRotAngles(config->GetMarker_All_TagBound(Local_Marker_PeriodicSend[iSend])); + su2double rotation[3][3], vector[3] = {}; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotation); + /*--- Axis points contribute the rotated normal of each additional copy. ---*/ + const auto nCopies = max(1ul, Local_Copy_PeriodicSend[iSend]); + for (auto iCopy = 0ul; iCopy < nCopies; ++iCopy) { + std::copy(buffer + vectorIndex, buffer + vectorIndex + nDim, vector); + for (auto iDim = 0u; iDim < nDim; ++iDim) { + buffer[vectorIndex + iDim] = 0; + for (auto jDim = 0u; jDim < nDim; ++jDim) + buffer[vectorIndex + iDim] += rotation[iDim][jDim] * vector[jDim]; + } + } + } + } +#ifdef HAVE_MPI + PostPeriodicSends(this, config, COMM_TYPE::DOUBLE, count, iMessage); +#else + const auto message = PeriodicRecv2Neighbor[rank]; + const auto start = count * nPoint_PeriodicSend[iMessage]; + const auto end = count * nPoint_PeriodicSend[iMessage + 1]; + std::copy(bufD_PeriodicSend + start, bufD_PeriodicSend + end, + bufD_PeriodicRecv + count * nPoint_PeriodicRecv[message]); +#endif + } + for (auto iMessage = 0; iMessage < nPeriodicRecv; ++iMessage) { + auto message = iMessage; +#ifdef HAVE_MPI + SU2_MPI::Status status; + int index; + SU2_MPI::Waitany(nPeriodicRecv, req_PeriodicRecv, &index, &status); + message = PeriodicRecv2Neighbor[status.MPI_SOURCE]; +#endif + for (auto iRecv = nPoint_PeriodicRecv[message]; iRecv < nPoint_PeriodicRecv[message + 1]; ++iRecv) { + const auto periodic = Local_Marker_PeriodicRecv[iRecv]; + if (periodic != iPair && periodic != iPair + config->GetnMarker_Periodic() / 2) continue; + const auto iPoint = Local_Point_PeriodicRecv[iRecv]; + for (auto iVar = 0u; iVar < count; ++iVar) values(iPoint, iVar) += bufD_PeriodicRecv[count * iRecv + iVar]; + } + } + SU2_MPI::Waitall(nPeriodicSend, req_PeriodicSend, MPI_STATUS_IGNORE); + } +} + +void CGeometry::UpdatePeriodicVolumes(const CConfig* config) { + if (!config->GetnMarker_Periodic() || (!nPeriodicSend && !nPeriodicRecv)) return; + AllocatePeriodicComms(1); + BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { + su2activematrix volumes(nPoint, 1); + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) volumes(iPoint, 0) = nodes->GetVolume(iPoint); + SumPeriodicGeometry(config, volumes, -1); + for (auto iPoint = 0ul; iPoint < nPointDomain; ++iPoint) + nodes->SetPeriodicVolume(iPoint, volumes(iPoint, 0) - nodes->GetVolume(iPoint)); + } + END_SU2_OMP_SAFE_GLOBAL_ACCESS } void CGeometry::AllocatePeriodicComms(unsigned short countPerPeriodicPoint) { @@ -2734,6 +2871,7 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) geometry_container[MESH_0]->SetControlVolume(config, UPDATE); geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry_container[MESH_0]->SetMaxLength(config); + geometry_container[MESH_0]->UpdatePeriodicVolumes(config); for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) { /*--- Update the control volume structures ---*/ @@ -2741,6 +2879,7 @@ void CGeometry::UpdateGeometry(CGeometry** geometry_container, CConfig* config) geometry_container[iMesh]->SetControlVolume(geometry_container[iMesh - 1], UPDATE); geometry_container[iMesh]->SetBoundControlVolume(geometry_container[iMesh - 1], config, UPDATE); geometry_container[iMesh]->SetCoord(geometry_container[iMesh - 1]); + geometry_container[iMesh]->UpdatePeriodicVolumes(config); } /*--- Compute the global surface areas for all markers. ---*/ @@ -2877,6 +3016,46 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) { } } + std::vector>> periodicNormals(nMarker); + std::vector periodicPoints(nPoint, false); + if (nPeriodicRecv) + for (auto i = 0; i < nPoint_PeriodicRecv[nPeriodicRecv]; ++i) periodicPoints[Local_Point_PeriodicRecv[i]] = true; + + /*--- A slip wall meeting a periodic face needs the normal of the complete + * wall patch. Keep markers separate to preserve distinct corner constraints. ---*/ + if (config->GetnMarker_Periodic() && (nPeriodicSend || nPeriodicRecv)) { + for (auto iCfgMarker = 0u; iCfgMarker < config->GetnMarker_CfgFile(); ++iCfgMarker) { + const auto tag = config->GetMarker_CfgFile_TagBound(iCfgMarker); + const auto kind = config->GetMarker_CfgFile_KindBC(tag); + if (kind != SYMMETRY_PLANE && kind != EULER_WALL) continue; + su2activematrix normals(nPoint, nDim); + normals = su2double(0); + for (const auto iMarker : symMarkers) { + if (config->GetMarker_All_TagBound(iMarker) != tag) continue; + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) { + const auto iPoint = vertex[iMarker][iVertex]->GetNode(); + for (auto iDim = 0u; iDim < nDim; ++iDim) normals(iPoint, iDim) = vertex[iMarker][iVertex]->GetNormal(iDim); + } + } + SumPeriodicGeometry(config, normals, 0); + for (const auto iMarker : symMarkers) { + if (config->GetMarker_All_TagBound(iMarker) != tag) continue; + for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; ++iVertex) { + const auto iPoint = vertex[iMarker][iVertex]->GetNode(); + if (!periodicPoints[iPoint]) continue; + const auto area = GeometryToolbox::Norm(nDim, normals[iPoint]); + if (area <= MIN_AREA) continue; + auto& normal = symmetryNormals[iMarker][iVertex]; + auto& areaNormal = periodicNormals[iMarker][iVertex]; + for (auto iDim = 0u; iDim < nDim; ++iDim) { + areaNormal[iDim] = normals(iPoint, iDim); + normal[iDim] = areaNormal[iDim] / area; + } + } + } + } + } + /*--- Merge the normals of curved markers to be independent of how surfaces are divided. ---*/ std::unordered_map> mergedNormals; @@ -2893,7 +3072,11 @@ void CGeometry::ComputeModifiedSymmetryNormals(const CConfig* config) { if (count < 2) continue; std::array normal = {}; - vertex[iMarker][iVertex]->GetNormal(normal.data()); + const auto corrected = periodicNormals[iMarker].find(iVertex); + if (corrected != periodicNormals[iMarker].end()) + normal = corrected->second; + else + vertex[iMarker][iVertex]->GetNormal(normal.data()); auto result = mergedNormals.emplace(iPoint, normal); const auto inserted = result.second; diff --git a/Common/src/geometry/CMultiGridGeometry.cpp b/Common/src/geometry/CMultiGridGeometry.cpp index b8caa00e9c94..9b1977e20996 100644 --- a/Common/src/geometry/CMultiGridGeometry.cpp +++ b/Common/src/geometry/CMultiGridGeometry.cpp @@ -1333,6 +1333,8 @@ void CMultiGridGeometry::SetBoundControlVolume(const CGeometry* fine_grid, const } END_SU2_OMP_FOR + /*--- Allocate on the whole team before the master computes periodic slip normals. ---*/ + if (nPeriodicSend || nPeriodicRecv) AllocatePeriodicComms(nDim); SU2_OMP_SAFE_GLOBAL_ACCESS(ComputeModifiedSymmetryNormals(config);) } diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 6928c2d922b5..18af5b585ca1 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -6984,6 +6984,8 @@ void CPhysicalGeometry::SetBoundControlVolume(const CConfig* config, unsigned sh } END_SU2_OMP_FOR + /*--- Allocate on the whole team before the master computes periodic slip normals. ---*/ + if (nPeriodicSend || nPeriodicRecv) AllocatePeriodicComms(nDim); SU2_OMP_SAFE_GLOBAL_ACCESS(ComputeModifiedSymmetryNormals(config);) } diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index ca6001b7f9f9..4e7677ea040c 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -45,6 +45,7 @@ dual mesh control volumes in the domain and on the boundaries. ---*/ geometry->SetControlVolume(config, UPDATE); geometry->SetBoundControlVolume(config, UPDATE); geometry->SetMaxLength(config); + geometry->UpdatePeriodicVolumes(config); } void CVolumetricMovement::UpdateMultiGrid(CGeometry** geometry, CConfig* config) { @@ -58,6 +59,7 @@ including computing the grid velocities on the coarser levels. ---*/ geometry[iMGlevel]->SetControlVolume(geometry[iMGfine], UPDATE); geometry[iMGlevel]->SetBoundControlVolume(geometry[iMGfine], config, UPDATE); geometry[iMGlevel]->SetCoord(geometry[iMGfine]); + geometry[iMGlevel]->UpdatePeriodicVolumes(config); if (config->GetGrid_Movement()) geometry[iMGlevel]->SetRestricted_GridVelocity(geometry[iMGfine]); } } diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 407e1c9bc7a3..4424f784a1a9 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -337,6 +337,25 @@ namespace PeriodicCommHelpers { break; } } + + /*--- Replaces a rotation matrix by the matrix that applies the same rotation nRot times. ---*/ + void applyRotationNTimes(unsigned long nRot, su2double rotMatrix[3][3]) { + su2double rotOnce[3][3], rotPrev[3][3]; + for (auto iDim = 0u; iDim < 3; iDim++) + for (auto jDim = 0u; jDim < 3; jDim++) + rotOnce[iDim][jDim] = rotMatrix[iDim][jDim]; + + for (auto iRot = 1ul; iRot < nRot; iRot++) { + for (auto iDim = 0u; iDim < 3; iDim++) + for (auto jDim = 0u; jDim < 3; jDim++) + rotPrev[iDim][jDim] = rotMatrix[iDim][jDim]; + + for (auto iDim = 0u; iDim < 3; iDim++) + for (auto jDim = 0u; jDim < 3; jDim++) + rotMatrix[iDim][jDim] = rotOnce[iDim][0]*rotPrev[0][jDim] + rotOnce[iDim][1]*rotPrev[1][jDim] + + rotOnce[iDim][2]*rotPrev[2][jDim]; + } + } } bool CSolver::HasPeriodicRotation(const CConfig& config) const { @@ -544,6 +563,32 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix3D); } + /*--- Points on a rotation axis are their own periodic match, they send once for every + other copy of their control volume, i.e. with the rotation applied more than once. ---*/ + + const auto nRot = geometry->Local_Copy_PeriodicSend[msg_offset + iSend]; + + if (nRot > 1) { + if (nDim==2) { + const su2double PsiTotal = nRot*Psi; + GeometryToolbox::RotationMatrix(PsiTotal, rotMatrix2D); + } else { + PeriodicCommHelpers::applyRotationNTimes(nRot, rotMatrix3D); + } + } + + /*--- The neighbors on periodic faces are skipped where edges must be counted once, because + the periodic match also has these edges. From the second copy of a point on a rotation + axis onwards, the edges on the face of this marker (not along the axis) are new. ---*/ + + const auto donorMarker = (nRot > 1) ? config->GetMarker_Periodic_Donor(Marker_Tag) : 0; + + auto sharedEdge = [&](unsigned long jPoint) { + if (!geometry->nodes->GetPeriodicBoundary(jPoint)) return false; + return (nRot < 2) || (geometry->nodes->GetVertex(jPoint, iPeriodic) < 0) || + (geometry->nodes->GetVertex(jPoint, donorMarker) >= 0); + }; + /*--- Whether the vector components of the solution are rotated for this marker. ---*/ const bool rotation = rotate_periodic && GeometryToolbox::HasRotation(angles); @@ -577,7 +622,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, that we avoid double counting neighbors on both sides. If not, increment the count of neighbors for the donor. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) + if (!sharedEdge(jPoint)) nNeighbor++; } @@ -695,7 +740,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + if (!sharedEdge(jPoint)) { /*--- Solution differences ---*/ @@ -761,7 +806,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid halos and boundary points so that we don't duplicate edges on both sides of the periodic BC. ---*/ - if (geometry->nodes->GetPeriodicBoundary(jPoint)) continue; + if (sharedEdge(jPoint)) continue; /*--- Use density instead of pressure for incomp. flows. ---*/ @@ -904,7 +949,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry, /*--- Avoid periodic boundary points so that we do not duplicate edges on both sides of the periodic BC. ---*/ - if (!geometry->nodes->GetPeriodicBoundary(jPoint)) { + if (!sharedEdge(jPoint)) { /*--- Get coordinates for the neighbor point. ---*/ @@ -1194,7 +1239,10 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, nRecv = (geometry->nPoint_PeriodicRecv[jRecv+1] - geometry->nPoint_PeriodicRecv[jRecv]); - SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + /*--- Points on a rotation axis receive more than once, in that case + the loop is not shared by the threads (one chunk). ---*/ + + SU2_OMP_FOR_STAT(geometry->PeriodicAxisPoints ? nRecv + 1 : size_t(OMP_MIN_SIZE)) for (iRecv = 0; iRecv < nRecv; iRecv++) { /*--- Get the local index for this communicated data. ---*/ @@ -1202,6 +1250,10 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, iPoint = geometry->Local_Point_PeriodicRecv[msg_offset + iRecv]; iPeriodic = geometry->Local_Marker_PeriodicRecv[msg_offset + iRecv]; + /*--- For points on a rotation axis, which copy of the control volume this is (0 otherwise). ---*/ + + const auto iCopy = geometry->Local_Copy_PeriodicRecv[msg_offset + iRecv]; + /*--- While all periodic face data was accumulated, we only store the values for the current pair of periodic faces. This is slightly inefficient when we have multiple pairs of periodic faces, but @@ -1273,6 +1325,18 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, Jacobian.AddBlock2Diag(iPoint, Jacobian_i); + /*--- A point on a rotation axis accumulates the residual and the diagonal block of all + the copies of its control volume, the blocks of its neighbors (which are the same in + all copies, up to the rotation) need the same factor, iCopy+1 after this copy. ---*/ + + if (iCopy > 0) { + const passivedouble factor = (iCopy + 1.0) / iCopy; + for (auto jPoint : geometry->nodes->GetPoints(iPoint)) { + auto* block = Jacobian.GetBlock(iPoint, jPoint); + for (auto iEntry = 0; iEntry < nVar*nVar; iEntry++) block[iEntry] *= factor; + } + } + if (iPeriodic == val_periodic_index + nPeriodic/2) { for (iVar = 0; iVar < nVar; iVar++) { LinSysRes(iPoint, iVar) = 0.0; @@ -1307,6 +1371,22 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, } + /*--- A point on a rotation axis takes the average over all the copies of its control + volume (running average over the copies received), which removes the velocity + normal to the axis: a vector that is the same in all the copies is along the axis. ---*/ + + else if (iCopy > 0) { + + for (auto iField = 0u; iField < nVar; iField++) { + const su2double average = base_nodes->GetSolution(iPoint, iField) + + (bufDRecv[buf_offset] - base_nodes->GetSolution(iPoint, iField)) / su2double(iCopy + 1); + base_nodes->SetSolution(iPoint, iField, average); + base_nodes->SetSolution_Old(iPoint, iField, average); + buf_offset++; + } + + } + break; case PERIODIC_LAPLACIAN: @@ -1314,8 +1394,8 @@ void CSolver::CompletePeriodicComms(CGeometry *geometry, /*--- Adjust the undivided Laplacian. The accumulation was with a subtraction before communicating, so now just add. ---*/ - for (iVar = 0; iVar < nVar; iVar++) - base_nodes->AddUnd_Lapl(iPoint, iVar, bufDRecv[buf_offset+iVar]); + for (auto iField = 0u; iField < nVar; iField++) + base_nodes->AddUnd_Lapl(iPoint, iField, bufDRecv[buf_offset+iField]); break; diff --git a/TestCases/navierstokes/periodic3D_axis/config.cfg b/TestCases/navierstokes/periodic3D_axis/config.cfg new file mode 100644 index 000000000000..fd1725ad14ac --- /dev/null +++ b/TestCases/navierstokes/periodic3D_axis/config.cfg @@ -0,0 +1,69 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Laminar flow in a pipe, 30 degree sector (one cell in the % +% circumferential direction) with nodes on the rotation axis % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +RESTART_SOL= NO + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +% The initial flow is along the axis (z) +MACH_NUMBER= 0.2 +AOA= 90.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_PRESSURE= 99000.0 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 80 +REYNOLDS_LENGTH= 1.0 +REF_DIMENSIONALIZATION= DIMENSIONAL +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_PERIODIC= ( per1, per2, 0,0,0, 0,0,30, 0,0,0 ) +MARKER_HEATFLUX= ( wall, 0.0 ) +INLET_TYPE= TOTAL_CONDITIONS +MARKER_INLET= ( inlet, 300.0, 101000.0, 0.0, 0.0, 1.0 ) +MARKER_OUTLET= ( outlet, 100000.0 ) +% +MARKER_MONITORING= ( wall ) +MARKER_ANALYZE= ( outlet ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 10.0 +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1e-6 +LINEAR_SOLVER_ITER= 20 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= NO + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD= RMS_DENSITY +CONV_RESIDUAL_MINVAL= -12 +ITER= 5000 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FORMAT= SU2 +MESH_FILENAME= wedge.su2 +OUTPUT_WRT_FREQ= 9999 +SCREEN_OUTPUT= ( INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-Z, RMS_ENERGY, SURFACE_MASSFLOW ) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 89bdeb0fe118..8756c29efcc2 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -476,6 +476,14 @@ def main(): periodic2d.test_vals = [-3.266248, -0.617784, -0.620252, 2.216476, -1.007761, -1035.000000] test_list.append(periodic2d) + # 3D rotational periodic pipe sector with nodes on the rotation axis + periodic3d_axis = TestCase('periodic3d_axis') + periodic3d_axis.cfg_dir = "navierstokes/periodic3D_axis" + periodic3d_axis.cfg_file = "config.cfg" + periodic3d_axis.test_iter = 100 + periodic3d_axis.test_vals = [-3.353165, 0.157501, 2.051458, -9.594400, -9.594400] + test_list.append(periodic3d_axis) + ########################## ### Compressible RANS ### ########################## diff --git a/UnitTests/Common/geometry/CGeometry_test.cpp b/UnitTests/Common/geometry/CGeometry_test.cpp index 22556b16055d..e2ff2695a208 100644 --- a/UnitTests/Common/geometry/CGeometry_test.cpp +++ b/UnitTests/Common/geometry/CGeometry_test.cpp @@ -26,7 +26,9 @@ */ #include "catch.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include "../../UnitQuadTestCase.hpp" +#include "../../../Common/include/geometry/CMultiGridGeometry.hpp" std::unique_ptr TestCase; @@ -140,3 +142,112 @@ TEST_CASE("Set bound control volume", "[Geometry]") { CHECK(TestCase->geometry->vertex[3][2]->GetNormal()[1] == -0.0625); CHECK(TestCase->geometry->vertex[5][3]->GetNormal()[2] == 0.03125); } + +TEST_CASE("Periodic slip-wall normal", "[Periodic]") { + const bool multigrid = GENERATE(false, true); + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, "MARKER_EULER= (y_minus,y_plus)\nMARKER_CUSTOM= (z_minus,z_plus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 30,0,0, 0,0,0)"); + if (multigrid) field.AddOption("MGLEVEL= 1"); + field.InitConfig(); + field.InitGeometry(true); + /*--- Bend the BOX into an annular sector about x: y is radius, x is azimuth. ---*/ + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + const su2double angle = x[0] * PI_NUMBER / 6, radius = 1 + x[1], axial = x[2]; + const su2double coordinate[] = {axial, -radius * sin(angle), radius * cos(angle)}; + field.geometry->nodes->SetCoord(i, coordinate); + } + field.geometry->SetControlVolume(field.config.get(), UPDATE); + field.geometry->SetBoundControlVolume(field.config.get(), UPDATE); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + std::unique_ptr coarse; + CGeometry* geometry = field.geometry.get(); + if (multigrid) { + coarse.reset(new CMultiGridGeometry(geometry, field.config.get(), 1)); + coarse->SetPoint_Connectivity(geometry); + coarse->SetEdges(); + coarse->SetVertex(geometry, field.config.get()); + coarse->SetControlVolume(geometry, ALLOCATE); + coarse->SetBoundControlVolume(geometry, field.config.get(), ALLOCATE); + coarse->SetCoord(geometry); + coarse->SetMGLevel(1); + coarse->MatchPeriodic(field.config.get(), 1); + coarse->PreprocessPeriodicComms(coarse.get(), field.config.get()); + geometry = coarse.get(); + } + su2double error = 0; + unsigned long checked = 0; + for (auto marker = 0u; marker < geometry->GetnMarker(); ++marker) { + if (field.config->GetMarker_All_TagBound(marker) != "y_plus") continue; + for (auto v = 0ul; v < geometry->GetnVertex(marker); ++v) { + const auto i = geometry->vertex[marker][v]->GetNode(); + const auto* x = geometry->nodes->GetCoord(i); + if (fabs(x[1]) > 1e-12 || fabs(x[2] - 2) > 1e-12 || !geometry->nodes->GetDomain(i)) continue; + const auto it = geometry->symmetryNormals[marker].find(v); + const auto* normal = + it == geometry->symmetryNormals[marker].end() ? geometry->vertex[marker][v]->GetNormal() : it->second.data(); + error = std::max(error, fabs(normal[1]) / GeometryToolbox::Norm(3, normal)); + ++checked; + } + } + unsigned long total = 0; + SU2_MPI::Allreduce(&checked, &total, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(total > 0); + CHECK(error < 1e-12); +} + +TEST_CASE("Periodic volume refresh after mesh deformation", "[Periodic]") { + UnitQuadTestCase field; + const auto start = field.config_options.find("MARKER_HEATFLUX="); + const auto end = field.config_options.find("VISCOSITY_MODEL="); + field.config_options.replace(start, end - start, "MARKER_CUSTOM= (y_minus,y_plus,z_plus,z_minus)\n"); + field.AddOption("MARKER_PERIODIC= (x_minus,x_plus, 0,0,0, 0,0,0, 1,0,0)"); + field.AddOption("DEFORM_MESH= YES"); + field.InitConfig(); + field.InitGeometry(true); + field.geometry->MatchPeriodic(field.config.get(), 1); + field.geometry->PreprocessPeriodicComms(field.geometry.get(), field.config.get()); + field.InitSolver(); + std::vector> original(field.geometry->GetnPoint()); + for (auto i = 0ul; i < original.size(); ++i) + for (auto d = 0u; d < 3; ++d) original[i][d] = field.geometry->nodes->GetCoord(i, d); + CGeometry* meshes[] = {field.geometry.get()}; + for (const auto amplitude : {0.05, -0.02, 0.0}) { + for (auto i = 0ul; i < original.size(); ++i) { + const auto& x = original[i]; + field.geometry->nodes->SetCoord(i, 0, x[0] + amplitude * sin(2 * PI_NUMBER * x[0]) * sin(PI_NUMBER * x[2])); + } + SU2_OMP_PARALLEL { CGeometry::UpdateGeometry(meshes, field.config.get()); } + END_SU2_OMP_PARALLEL + su2double expected = 0, stored = 0; + for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + if (fabs(x[1] - 0.5) > 1e-12 || fabs(x[2] - 0.5) > 1e-12) continue; + if (fabs(x[0] - 1) < 1e-12) expected = field.geometry->nodes->GetVolume(i); + if (fabs(x[0]) < 1e-12) stored = field.geometry->nodes->GetPeriodicVolume(i); + } + su2double totals[2] = {stored, expected}, global[2] = {}; + SU2_MPI::Allreduce(totals, global, 2, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + REQUIRE(global[1] > 0); + CHECK(global[0] == Approx(global[1]).margin(1e-12)); + auto* solver = field.solver[FLOW_SOL]; + auto* nodes = solver->GetNodes(); + for (auto i = 0ul; i < field.geometry->GetnPoint(); ++i) + for (auto v = 0u; v < solver->GetnPrimVarGrad(); ++v) + nodes->SetPrimitive(i, v, 2 + field.geometry->nodes->GetCoord(i, 1)); + SU2_OMP_PARALLEL { solver->SetPrimitive_Gradient_GG(field.geometry.get(), field.config.get()); } + END_SU2_OMP_PARALLEL + su2double error = 0; + for (auto i = 0ul; i < field.geometry->GetnPointDomain(); ++i) { + const auto* x = field.geometry->nodes->GetCoord(i); + if (fabs(x[1] - 0.5) > 1e-12 || fabs(x[2] - 0.5) > 1e-12) continue; + if (fabs(x[0]) > 1e-12 && fabs(x[0] - 1) > 1e-12) continue; + error = std::max(error, fabs(nodes->GetGradient_Primitive()(i, 0, 1) - 1)); + } + CHECK(error < 1e-12); + } +} diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index d8da75dc0960..b6e8bfdd6bc3 100644 --- a/UnitTests/UnitQuadTestCase.hpp +++ b/UnitTests/UnitQuadTestCase.hpp @@ -83,10 +83,13 @@ struct UnitQuadTestCase { /*! * \brief Initialize the geometry */ - void InitGeometry() { + void InitGeometry(bool partition = false) { cout.rdbuf(nullptr); { - auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), 0, 1)); + const auto rank = partition ? SU2_MPI::GetRank() : 0; + const auto size = partition ? SU2_MPI::GetSize() : 1; + auto aux_geometry = std::unique_ptr(new CPhysicalGeometry(config.get(), rank, size)); + if (partition) aux_geometry->SetColorGrid_Parallel(config.get()); geometry = std::unique_ptr(new CPhysicalGeometry(aux_geometry.get(), config.get())); } geometry->SetSendReceive(config.get());