diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 494c95429221..3a5f64a4fdf6 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -996,6 +996,8 @@ class CConfig {
bool DeadLoad; /*!< \brief Application of dead loads to the FE analysis */
bool PseudoStatic; /*!< \brief Application of dead loads to the FE analysis */
bool SteadyRestart; /*!< \brief Restart from a steady state for FSI problems. */
+ bool Turbo_MultiPsgs; /*!< \brief Restart multi passage simulation from a single passage solution for turbomachinery application. */
+ unsigned short nPassages; /*!< \brief Number of user defined passages. */
su2double Newmark_beta, /*!< \brief Parameter alpha for Newmark method. */
Newmark_gamma; /*!< \brief Parameter delta for Newmark method. */
unsigned short nIntCoeffs; /*!< \brief Number of integration coeffs for structural calculations. */
@@ -5219,6 +5221,18 @@ class CConfig {
*/
void SetnBlades(unsigned short val_iZone, su2double nblades) { nBlades[val_iZone] = nblades;}
+ /*!
+ * \brief Get the number of user defined passages to be modelled in computation.
+ * \param[in] nPassages - user defined number of passages
+ */
+ unsigned short GetnPassages(void) const { return nPassages; }
+
+ /*!
+ * \brief Identifies if we want to restart multi passage simulation from a single passage solution for turbomachinery case.
+ * \return TRUE if we restart from multi passage simulation, FALSE otherwise.
+ */
+ bool GetTurbo_MultiPsgs(void) const { return Turbo_MultiPsgs; }
+
/*!
* \brief Verify if there is any Giles Boundary Condition option specified from config file.
* \return boolean.
diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp
index 582679ca6afa..ff77217cbc40 100644
--- a/Common/include/geometry/CGeometry.hpp
+++ b/Common/include/geometry/CGeometry.hpp
@@ -317,6 +317,54 @@ class CGeometry {
SU2_MPI::Request* req_PeriodicSend{nullptr}; /*!< \brief Data structure for periodic send requests. */
SU2_MPI::Request* req_PeriodicRecv{nullptr}; /*!< \brief Data structure for periodic recv requests. */
+ /*--- Data structures for mesh Full annulus duplication. ---*/
+
+ int nPassages{1}; /*!< \brief the user defined number of passages. Used only for mesh duplication.*/
+ int nPsgs_FullAnnu{0}; /*!< \brief the number of passages for full annulus. Used only for mesh duplication.*/
+ bool Fullannus = false; /*!< \brief whether it's full annulus. Used only for mesh duplication.*/
+ unsigned long Global_nPoint_OldPsg{
+ 0}; /*!< \brief the number of points in the original one passage. Used only for mesh duplication.*/
+ su2double PitchAngle; /*!< \brief the pitch angle for one passage. Used only for mesh duplication.*/
+ int* nPerPointOnRank_donor = new int[size]; /*!< \brief vector contains the number of periodic points on each node.
+ Used only for mesh duplication. */
+ int* nPerPointOnRank_target = new int[size]; /*!< \brief vector contains the number of periodic points on each node.
+ Used only for mesh duplication. */
+
+ unsigned long nPerPointAll_donor = 0;
+ unsigned long nPerPointAll_target = 0;
+ unsigned long nPsgPoint_RemPer =
+ 0; /*!< \brief the number of points in one passage after removing donor points. Used only for mesh duplication.*/
+ unsigned long* PerPointsMatched_donor{nullptr}; /*!< \brief Global index of the matched points on the target of one
+ pair periodic BC. Used only for mesh duplication.*/
+ unsigned long* PerPointsGlbIndx_donor{nullptr}; /*!< \brief Global index of the points on the donor of one pair
+ periodic BC. Used only for mesh duplication.*/
+ su2double* PerPointsCoord_donor{nullptr}; /*!< \brief Coordinates of the points on the donor of one pair periodic BC.
+ Used only for mesh duplication.*/
+ unsigned long* PsgPointsOldGlbIndx_RemPer{nullptr}; /*!< \brief new Global index of the points in one passage after
+ removing donor points. Used only for mesh duplication.*/
+ su2double* PsgPointsCoord_RemPer{nullptr}; /*!< \brief coordinates of the points in one passage after removing donor
+ points. Used only for mesh duplication.*/
+ unordered_map
+ PerPoint_Match_Map; /*!< \brief Global index map for matched points on periodic boundaries. Used only for mesh
+ duplication.*/
+ unordered_map
+ Donor_Indx_Map; /*!< \brief Global to donor periodic boudnry index map. Used only for mesh duplication.*/
+ unordered_map Psg_Old2New_Map; /*!< \brief Old to new global index map for points in one
+ passage. Used only for mesh duplication.*/
+ unordered_map Psg_Old2New_Map2; /*!< \brief Old to new global index map for points in one
+ passage. Used only for mesh duplication.*/
+
+ // to be removed
+ unsigned long* nElem_Bound_FullAnnu{nullptr};
+ unsigned long Global_nPoint_Fullannulus = 0; /*!< \brief the new global number of points for full annulus*/
+ unsigned long nPoint_FullAnnu = 0; /*!< \brief the new local number of points for full annulus*/
+ unsigned long* LocalPointsGlbIndx_FullAnnu{nullptr};
+ su2double* LocalPointsCoord_FullAnnu{nullptr};
+ unsigned long numberOfLocalElements_FullAnnu =
+ 0; /*!< \brief Number of local elements on this rank for full annulus. */
+ vector localVolumeElementConnectivity_FullAnnu; /*!< \brief Vector containing the element connectivity
+ from the mesh file for the local elements. */
+
/*--- Mesh quality metrics. ---*/
vector Orthogonality; /*!< \brief Measure of dual CV orthogonality angle (0 to 90 deg., 90 being best). */
@@ -1871,4 +1919,86 @@ class CGeometry {
* \return A pointer to the reference node coordinate vector.
*/
inline virtual const su2double* GetStreamwise_Periodic_RefNode() const { return nullptr; }
+
+ /*!
+ * \brief Set the number of passages.
+ * \return The number of passages.
+ */
+ // inline void SetGlobal_nPoint_OldPsg(unsigned long val_nPoint) { Global_nPoint_OldPsg = val_nPoint; }
+
+ /*!
+ * \brief Retrieve the bool Fullannus.
+ * \return bool Fullannus.
+ */
+ inline bool GetFullannus() const { return Fullannus; }
+
+ /*!
+ * \brief Retrieve the old number of points before full annulus transformation.
+ * \return The old number of points before full annulus transformation.
+ */
+ inline unsigned long GetGlobal_nPoint_OldPsg() const { return Global_nPoint_OldPsg; }
+
+ /*!
+ * \brief Retrieve the new number of points in one passage after removing one periodic boundary.
+ * \return The new number of points in one passage after removing one periodic boundary.
+ */
+ inline unsigned long GetGlobal_nPsgPoint_RemPer() const { return nPsgPoint_RemPer; }
+
+ /*!
+ * \brief Retrieve the number of passages.
+ * \return The number of passages.
+ */
+ inline int GetnPassages() const { return nPassages; }
+
+ /*!
+ * \brief Retrieve the number of passages for full annulus.
+ * \return The number of passages.
+ */
+ inline int GetnPassages_FullAnnu() const { return nPsgs_FullAnnu; }
+
+ /*!
+ * \brief Retrieve the pitch angle for duplicating passages.
+ * \return The pitch angle for duplicating passages.
+ */
+ inline su2double GetPitchangle() const { return PitchAngle; }
+
+ /*!
+ * \brief Get the matched point's global index for periodic BC.
+ * \param[in] val_ipoint - Global point.
+ * \return Local index that correspond with the global index, -1 if not found on the current rank (process).
+ */
+ inline void Set_PsgOld2New_Indx_2(unsigned long val_old_indx, unsigned long val_new_indx) {
+ Psg_Old2New_Map2[val_old_indx] = val_new_indx;
+ }
+
+ /*!
+ * \brief Get the matched point's global index for periodic BC.
+ * \param[in] val_ipoint - Global point.
+ * \return Local index that correspond with the global index, -1 if not found on the current rank (process).
+ */
+ inline void Set_Donor_Indx(unsigned long val_donor_indx, unsigned long val_new_indx) {
+ Donor_Indx_Map[val_donor_indx] = val_new_indx;
+ }
+
+ /*!
+ * \brief Get the matched point's global index for periodic BC.
+ * \param[in] val_ipoint - Global point.
+ * \return Local index that correspond with the global index, -1 if not found on the current rank (process).
+ */
+ inline unsigned long Get_PsgOld2New_Indx_2(unsigned long val_old_indx) {
+ unsigned long val_new_indx = Psg_Old2New_Map2[val_old_indx];
+ // return Psg_Old2New_Map[val_old_indx];
+ return val_new_indx;
+ }
+
+ /*!
+ * \brief Get the matched point's global index for periodic BC.
+ * \param[in] val_ipoint - Global point.
+ * \return Local index that correspond with the global index, -1 if not found on the current rank (process).
+ */
+ inline unsigned long Get_Donor_Indx(unsigned long val_donor_indx) {
+ unsigned long val_new_indx = Donor_Indx_Map[val_donor_indx];
+ // return Psg_Old2New_Map[val_old_indx];
+ return val_new_indx;
+ }
};
diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp
index 10e781cb85e9..42d2596b835e 100644
--- a/Common/include/geometry/CPhysicalGeometry.hpp
+++ b/Common/include/geometry/CPhysicalGeometry.hpp
@@ -328,6 +328,41 @@ class CPhysicalGeometry final : public CGeometry {
*/
void LoadUnpartitionedSurfaceElements(CConfig* config, CMeshReaderFVM* mesh);
+ /*!
+ * \brief Match the periodic boundary, useful for full annulus mesh duplication.
+ * \param[in] config - definition of the particular problem.
+ * \param[in] val_periodic - index of the periodic pair that is going to be treated.
+ */
+ void MatchPeriodicPoints(CConfig* config, unsigned short val_periodic);
+
+ /*!
+ * \brief Reset the rotating angle for periodic boudary conditions, useful for full annulus mesh duplication.
+ * \param[in] config - definition of the particular problem.
+ * \param[in] val_periodic - index of the periodic pair that is going to be treated.
+ */
+ void ReSetPeriodicAngle_Z(CConfig* config);
+
+ /*!
+ * \brief duplicate the points, and redistribute points on each node.
+ * \param[in] config - definition of the particular problem.
+ * \param[in] mesh - mesh reader object containing the current zone data.
+ */
+ void DuplicatePoints(CConfig* config, CMeshReaderFVM* mesh);
+
+ /*!
+ * \brief duplicate the volume elems, and redistribute volume elems on each node.
+ * \param[in] config - definition of the particular problem.
+ * \param[in] mesh - mesh reader object containing the current zone data.
+ */
+ void DuplicateVolumeElems(CConfig* config, CMeshReaderFVM* mesh);
+
+ /*!
+ * \brief duplicate the surface elems, and redistribute surface elems on each node.
+ * \param[in] config - definition of the particular problem.
+ * \param[in] mesh - mesh reader object containing the current zone data.
+ */
+ void DuplicateSurfaceElems(CConfig* config, CMeshReaderFVM* mesh);
+
/*!
* \brief Prepares the grid point adjacency based on a linearly partitioned mesh object needed by ParMETIS for graph
* partitioning in parallel. \param[in] config - Definition of the particular problem.
diff --git a/Common/include/parallelization/mpi_structure.hpp b/Common/include/parallelization/mpi_structure.hpp
index b212d770e8f6..0aee9610dc67 100644
--- a/Common/include/parallelization/mpi_structure.hpp
+++ b/Common/include/parallelization/mpi_structure.hpp
@@ -218,6 +218,16 @@ class CBaseMPIWrapper {
MPI_Scatter(sendbuf, sendcnt, sendtype, recvbuf, recvcnt, recvtype, root, comm);
}
+ static inline void Scatterv(const void* sendbuf, const int* sendcnts, const int* displs, Datatype sendtype,
+ void* recvbuf, int recvcnt, Datatype recvtype, int root, Comm comm) {
+ MPI_Scatterv(sendbuf, sendcnts, displs, sendtype, recvbuf, recvcnt, recvtype, root, comm);
+ }
+
+ static inline void Gatherv(const void* sendbuf, int sendcnt, Datatype sendtype, void* recvbuf, const int* recvcounts,
+ const int* displs, Datatype recvtype, int root, Comm comm) {
+ MPI_Gatherv(sendbuf, sendcnt, sendtype, recvbuf, recvcounts, displs, recvtype, root, comm);
+ }
+
static inline void Allgather(const void* sendbuf, int sendcnt, Datatype sendtype, void* recvbuf, int recvcnt,
Datatype recvtype, Comm comm) {
MPI_Allgather(sendbuf, sendcnt, sendtype, recvbuf, recvcnt, recvtype, comm);
@@ -425,6 +435,18 @@ class CMediMPIWrapper : public CBaseMPIWrapper {
convertComm(comm));
}
+ static inline void Scatterv(const void* sendbuf, const int* sendcnts, const int* displs, Datatype sendtype,
+ void* recvbuf, int recvcnt, Datatype recvtype, int root, Comm comm) {
+ AMPI_Scatterv(sendbuf, sendcnts, displs, convertDatatype(sendtype), recvbuf, recvcnt, convertDatatype(recvtype),
+ root, comm);
+ }
+
+ static inline void Gatherv(const void* sendbuf, int sendcnt, Datatype sendtype, void* recvbuf, const int* recvcounts,
+ const int* displs, Datatype recvtype, int root, Comm comm) {
+ MPI_Gatherv(sendbuf, sendcnt, convertDatatype(sendtype), recvbuf, recvcounts, displs, convertDatatype(recvtype),
+ root, comm);
+ }
+
static inline void Allgather(const void* sendbuf, int sendcnt, Datatype sendtype, void* recvbuf, int recvcnt,
Datatype recvtype, Comm comm) {
AMPI_Allgather(sendbuf, sendcnt, convertDatatype(sendtype), recvbuf, recvcnt, convertDatatype(recvtype),
@@ -572,6 +594,16 @@ class CBaseMPIWrapper {
CopyData(sendbuf, recvbuf, sendcnt, sendtype);
}
+ static inline void Scatterv(const void* sendbuf, const int* sendcnts, const int* displs, Datatype sendtype,
+ void* recvbuf, int recvcnt, Datatype recvtype, int root, Comm comm) {
+ CopyData(sendbuf, recvbuf, sendcnts[0], sendtype);
+ }
+
+ static inline void Gatherv(const void* sendbuf, int sendcnt, Datatype sendtype, void* recvbuf, const int* recvcounts,
+ const int* displs, Datatype recvtype, int root, Comm comm) {
+ CopyData(sendbuf, recvbuf, sendcnt, sendtype);
+ }
+
static inline void Allgatherv(const void* sendbuf, int sendcnt, Datatype sendtype, void* recvbuf, const int* recvcnt,
const int* displs, Datatype recvtype, Comm comm) {
CopyData(sendbuf, recvbuf, sendcnt, sendtype, displs[0]);
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 5ff351cfee5c..3f34895ae76e 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -1631,6 +1631,10 @@ void CConfig::SetConfig_Options() {
addTurboPerfOption("MARKER_TURBOMACHINERY", nMarker_Turbomachinery, Marker_TurboBoundIn, Marker_TurboBoundOut);
/*!\brief NUM_SPANWISE_SECTIONS \n DESCRIPTION: Integer number of spanwise sections to compute 3D turbo BC and Performance for turbomachinery */
addUnsignedShortOption("NUM_SPANWISE_SECTIONS", nSpanWiseSections_User, 1);
+ /*!\brief PassageNumber\n DESCRIPTION: User defined number of passages. DEFAULT: 1 \ingroup Config*/
+ addUnsignedShortOption("PASSAGE_NUMBER", nPassages, 1);
+ /*!\brief TURBOMACHINERY_MULTIPASSAGE_RESTART \n DESCRIPTION: Restart multi passage simulation from single passage results \n Options: NO, YES \ingroup Config */
+ addBoolOption("TURBOMACHINERY_MULTIPASSAGE_RESTART", Turbo_MultiPsgs, false);
/*!\brief SPANWISE_KIND \n DESCRIPTION: type of algorithm to identify the span-wise sections at the turbo boundaries.
\n OPTIONS: see \link SpanWise_Map \endlink \n Default: AUTOMATIC */
addEnumOption("SPANWISE_KIND", Kind_SpanWise, SpanWise_Map, AUTOMATIC);
diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp
index 29f10a492b30..b706666a0b82 100644
--- a/Common/src/geometry/CPhysicalGeometry.cpp
+++ b/Common/src/geometry/CPhysicalGeometry.cpp
@@ -188,6 +188,23 @@ CPhysicalGeometry::CPhysicalGeometry(CGeometry* geometry, CConfig* config) : CGe
nDim = geometry->GetnDim();
nZone = geometry->GetnZone();
+ if (config->GetTurbo_MultiPsgs()) {
+ /*--- Set important info regarding full annulus transformation in the new geometry class ---*/
+ nPassages = config->GetnPassages();
+ PitchAngle = geometry->GetPitchangle();
+ Global_nPoint_OldPsg = geometry->GetGlobal_nPoint_OldPsg();
+ nPsgPoint_RemPer = geometry->GetGlobal_nPsgPoint_RemPer();
+ Fullannus = geometry->GetFullannus();
+ unsigned long new_indx;
+ for (unsigned long iPoint = 0; iPoint < Global_nPoint_OldPsg; iPoint++) {
+ new_indx = geometry->Get_PsgOld2New_Indx_2(iPoint);
+ Set_PsgOld2New_Indx_2(iPoint, new_indx);
+ if (new_indx == -1) {
+ Set_Donor_Indx(iPoint, geometry->Get_Donor_Indx(iPoint));
+ }
+ }
+ }
+
/*--- Recompute the linear partitioning offsets. ---*/
PrepareOffsets(geometry->GetGlobal_nPoint());
@@ -3521,6 +3538,49 @@ void CPhysicalGeometry::Read_Mesh_FVM(CConfig* config, const string& val_mesh_fi
LoadLinearlyPartitionedVolumeElements(config, MeshFVM);
LoadUnpartitionedSurfaceElements(config, MeshFVM);
+ if (config->GetTurbo_MultiPsgs()) {
+ /*--- Match the periodic boundaries. ---*/
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Required] "
+ << " Matching periodic points. " << endl;
+
+ nPassages = config->GetnPassages();
+
+ for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
+ MatchPeriodicPoints(config, iPeriodic);
+ }
+
+ /*--- Reinitialize counters for local/global points & elements ---*/
+
+ Global_nPoint_OldPsg = Global_nPointDomain;
+
+ Global_nPoint = 0;
+ Global_nPointDomain = 0;
+ Global_nElem = 0;
+ Global_nElemDomain = 0;
+ nelem_edge = 0;
+ Global_nelem_edge = 0;
+ nelem_triangle = 0;
+ Global_nelem_triangle = 0;
+ nelem_quad = 0;
+ Global_nelem_quad = 0;
+ nelem_tetra = 0;
+ Global_nelem_tetra = 0;
+ nelem_hexa = 0;
+ Global_nelem_hexa = 0;
+ nelem_prism = 0;
+ Global_nelem_prism = 0;
+ nelem_pyramid = 0;
+ Global_nelem_pyramid = 0;
+
+ /*--- Mesh duplication ---*/
+ DuplicatePoints(config, MeshFVM);
+ DuplicateVolumeElems(config, MeshFVM);
+ DuplicateSurfaceElems(config, MeshFVM);
+
+ // ReSetPeriodicAngle_Z(config);
+ }
+
/*--- Prepare the nodal adjacency structures for ParMETIS. ---*/
PrepareAdjacency(config);
@@ -3770,6 +3830,1309 @@ void CPhysicalGeometry::LoadUnpartitionedSurfaceElements(CConfig* config, CMeshR
}
}
+void CPhysicalGeometry::MatchPeriodicPoints(CConfig* config, unsigned short val_periodic) {
+ /*---We take the advantage that up to now all the marker elements
+ are stored in the master node, thus spare some coding effort of MPI
+ communication, and do most of the matching work in serial.
+ Here we try to align with the MatchPeriodic() function ---*/
+
+ int* displs_donor = new int[size];
+ int* displs_target = new int[size];
+ vector PerPoints_donor; /*!< \brief Vector to store the donor points (as global index). */
+ vector PerPoints_target; /*!< \brief Vector to store the target points (as global index). */
+ vector::iterator it;
+
+ unsigned long Global_Index;
+ unsigned short iDim;
+ unsigned long iNode, iPoint, jPoint;
+ unsigned short iMarker;
+ unsigned short iPeriodic, nPeriodic;
+ int PeriodicFound = 0;
+
+ if (rank == MASTER_NODE) {
+ /*--- Get the number of periodic boundaries ---*/
+ nPeriodic = config->GetnMarker_Periodic();
+
+ /*--- Send an initial message to the console. ---*/
+ cout << "Mesh duplication pre: Matching the periodic boundary points for marker pair ";
+ cout << val_periodic << "." << endl;
+
+ /*--- Find all points on the periodic boundary pair and load into our list. ---*/
+ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+ if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
+ iPeriodic = config->GetMarker_All_PerBound(iMarker);
+
+ if (iPeriodic == val_periodic) {
+ PeriodicFound = 1;
+
+ /*--- Loop over each point on the periodic marker and find the matching point from the donor marker. ---*/
+ for (unsigned long iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) {
+ unsigned short NODES_PER_ELEMENT = 0;
+ unsigned short Elem_Type = bound[iMarker][iElem]->GetVTK_Type();
+ switch (Elem_Type) {
+ case LINE:
+ NODES_PER_ELEMENT = N_POINTS_LINE;
+ break;
+ case TRIANGLE:
+ NODES_PER_ELEMENT = N_POINTS_TRIANGLE;
+ break;
+ case QUADRILATERAL:
+ NODES_PER_ELEMENT = N_POINTS_QUADRILATERAL;
+ break;
+ default:
+ SU2_MPI::Error("Unrecognized element type.", CURRENT_FUNCTION);
+ break;
+ }
+
+ for (iNode = 0; iNode < NODES_PER_ELEMENT; iNode++) {
+ /*--- Get the global index of the current point. ---*/
+ Global_Index = bound[iMarker][iElem]->GetNode(iNode);
+ PerPoints_donor.push_back(Global_Index);
+ }
+ }
+ }
+
+ if (iPeriodic == val_periodic + nPeriodic / 2) {
+ PeriodicFound = 1;
+
+ /*--- Loop over each point on the periodic marker and find the matching point from the donor marker. ---*/
+ for (unsigned long iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) {
+ unsigned short NODES_PER_ELEMENT = 0;
+ unsigned short Elem_Type = bound[iMarker][iElem]->GetVTK_Type();
+ switch (Elem_Type) {
+ case LINE:
+ NODES_PER_ELEMENT = N_POINTS_LINE;
+ break;
+ case TRIANGLE:
+ NODES_PER_ELEMENT = N_POINTS_TRIANGLE;
+ break;
+ case QUADRILATERAL:
+ NODES_PER_ELEMENT = N_POINTS_QUADRILATERAL;
+ break;
+ default:
+ SU2_MPI::Error("Unrecognized element type.", CURRENT_FUNCTION);
+ break;
+ }
+
+ for (iNode = 0; iNode < NODES_PER_ELEMENT; iNode++) {
+ /*--- Get the global index of the current point. ---*/
+ Global_Index = bound[iMarker][iElem]->GetNode(iNode);
+ PerPoints_target.push_back(Global_Index);
+ }
+ }
+ }
+ }
+ }
+
+ if (PeriodicFound == 1) {
+ for (int iRank = 0; iRank < size; iRank++) {
+ nPerPointOnRank_donor[iRank] = 0;
+ nPerPointOnRank_target[iRank] = 0;
+ }
+
+ /*--- Post-process the PerPoints vector. ---*/
+ sort(PerPoints_donor.begin(), PerPoints_donor.end());
+ it = unique(PerPoints_donor.begin(), PerPoints_donor.end());
+ PerPoints_donor.resize(it - PerPoints_donor.begin());
+
+ sort(PerPoints_target.begin(), PerPoints_target.end());
+ it = unique(PerPoints_target.begin(), PerPoints_target.end());
+ PerPoints_target.resize(it - PerPoints_target.begin());
+
+ /*--- Get a linear partitioner to track the partition counts. ---*/
+ CLinearPartitioner pointPartitioner(Global_nPointDomain, 0);
+
+ /*--- prepare the index vectors for MPI communication ---*/
+
+ displs_donor[0] = 0;
+ for (int iRank = 0; iRank < size; iRank++) {
+ int count = 0;
+ for (iPoint = 0; iPoint < PerPoints_donor.size(); iPoint++) {
+ if (PerPoints_donor[iPoint] >= pointPartitioner.GetFirstIndexOnRank(iRank) &&
+ PerPoints_donor[iPoint] < pointPartitioner.GetLastIndexOnRank(iRank)) {
+ count++;
+ }
+ }
+ nPerPointOnRank_donor[iRank] = count;
+ if (iRank != size - 1) displs_donor[iRank + 1] = displs_donor[iRank] + nPerPointOnRank_donor[iRank];
+ }
+
+ displs_target[0] = 0;
+ for (int iRank = 0; iRank < size; iRank++) {
+ int count = 0;
+ for (iPoint = 0; iPoint < PerPoints_target.size(); iPoint++) {
+ if (PerPoints_target[iPoint] >= pointPartitioner.GetFirstIndexOnRank(iRank) &&
+ PerPoints_target[iPoint] < pointPartitioner.GetLastIndexOnRank(iRank)) {
+ count++;
+ }
+ }
+ nPerPointOnRank_target[iRank] = count;
+ if (iRank != size - 1) displs_target[iRank + 1] = displs_target[iRank] + nPerPointOnRank_target[iRank];
+ }
+
+ /*
+ for (int iRank = 0; iRank < size; iRank++) {
+ cout<<"(Master Node) nPerPoint_donor on Rank "<< iRank <<" : "<GetGlobalIndex(jPoint)) {
+ su2double* Coord_j = nodes->GetCoord(jPoint);
+ for (iDim = 0; iDim < nDim; iDim++) Buffer_Send_Coord_donor[iPoint * nDim + iDim] = Coord_j[iDim];
+ }
+ }
+ }
+
+ for (iPoint = 0; iPoint < nPerPointOnRank_target[rank]; iPoint++) {
+ for (jPoint = 0; jPoint < nPoint; jPoint++) {
+ if (PerPointsOnRank_target[iPoint] == nodes->GetGlobalIndex(jPoint)) {
+ su2double* Coord_j = nodes->GetCoord(jPoint);
+ for (iDim = 0; iDim < nDim; iDim++) Buffer_Send_Coord_target[iPoint * nDim + iDim] = Coord_j[iDim];
+ }
+ }
+ }
+
+ int* CoordCountOnRank_donor = new int[size];
+ int* CoordCountOnRank_target = new int[size];
+ int* displs_coord_donor = new int[size];
+ int* displs_coord_target = new int[size];
+
+ for (int iRank = 0; iRank < size; iRank++) {
+ CoordCountOnRank_donor[iRank] = nPerPointOnRank_donor[iRank] * nDim;
+ CoordCountOnRank_target[iRank] = nPerPointOnRank_target[iRank] * nDim;
+ displs_coord_donor[iRank] = displs_donor[iRank] * nDim;
+ displs_coord_target[iRank] = displs_target[iRank] * nDim;
+ }
+
+ /*--- Gather the coordinates of all periodic points. ---*/
+ SU2_MPI::Gatherv(Buffer_Send_Coord_donor, nPerPointOnRank_donor[rank] * nDim, MPI_DOUBLE, Buffer_Recv_Coord_donor,
+ CoordCountOnRank_donor, displs_coord_donor, MPI_DOUBLE, MASTER_NODE, SU2_MPI::GetComm());
+ SU2_MPI::Gatherv(Buffer_Send_Coord_target, nPerPointOnRank_target[rank] * nDim, MPI_DOUBLE,
+ Buffer_Recv_Coord_target, CoordCountOnRank_target, displs_coord_target, MPI_DOUBLE, MASTER_NODE,
+ SU2_MPI::GetComm());
+
+ /*---Match the points based on rotation.---*/
+ string Marker_Tag;
+ su2double translation[3] = {0.0, 0.0, 0.0}, dx, dy, dz;
+ const su2double *center, *angles, *trans;
+ su2double Coord_i[3], Coord_j[3];
+ su2double rotCoord[3] = {0.0, 0.0, 0.0};
+ su2double rotMatrix[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+ su2double Theta, Phi, Psi, cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi;
+ su2double dist;
+ unsigned long pPoint = 0;
+ unsigned long nPointMatch = 0;
+ su2double mindist;
+ /*--- Tolerance for distance-based match to report warning. ---*/
+ su2double epsilon = 1e-6;
+
+ /*---Global vectors store the donor and matched point's global index---*/
+ PerPointsMatched_donor = new unsigned long[nPerPointAll_donor];
+ PerPointsGlbIndx_donor = new unsigned long[nPerPointAll_donor];
+ PerPointsCoord_donor = new su2double[nPerPointAll_donor * nDim];
+
+ /*--- Again, the master node takes care of matching.---*/
+ if (rank == MASTER_NODE) {
+ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+ if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
+ /*--- Find the donor boundary. ---*/
+ iPeriodic = config->GetMarker_All_PerBound(iMarker);
+ if (iPeriodic == val_periodic) {
+ /*--- Retrieve the supplied periodic information. ---*/
+ Marker_Tag = config->GetMarker_All_TagBound(iMarker);
+ center = config->GetPeriodicRotCenter(Marker_Tag);
+ angles = config->GetPeriodicRotAngles(Marker_Tag);
+ trans = config->GetPeriodicTranslation(Marker_Tag);
+
+ /*--- Store the pitch angle for duplication ---*/
+ /*--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ---*/
+ /*--- Note that we assume the rotating axis is Z-axis
+ And we reverse the duplication direction, so as from target to donor.
+ This is very important, otherwise the renumbering will goes wrong.
+ Later we will delete the donor side ---*/
+ /*--- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ---*/
+
+ // PitchAngle is used for grids duplication.
+ PitchAngle = -angles[2];
+ nPsgs_FullAnnu = abs(int(2 * PI_NUMBER / PitchAngle));
+
+ /*--- Store (center+trans) as it is constant and will be added. ---*/
+ translation[0] = center[0] + trans[0];
+ translation[1] = center[1] + trans[1];
+ translation[2] = center[2] + trans[2];
+
+ /*--- Store angles separately for clarity. Compute sines/cosines. ---*/
+ Theta = angles[0];
+ Phi = angles[1];
+ Psi = angles[2];
+ cosTheta = cos(Theta);
+ cosPhi = cos(Phi);
+ cosPsi = cos(Psi);
+ sinTheta = sin(Theta);
+ sinPhi = sin(Phi);
+ sinPsi = sin(Psi);
+
+ /*--- Compute the rotation matrix. Note that the implicit
+ ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+ rotMatrix[0][0] = cosPhi * cosPsi;
+ rotMatrix[1][0] = cosPhi * sinPsi;
+ rotMatrix[2][0] = -sinPhi;
+
+ rotMatrix[0][1] = sinTheta * sinPhi * cosPsi - cosTheta * sinPsi;
+ rotMatrix[1][1] = sinTheta * sinPhi * sinPsi + cosTheta * cosPsi;
+ rotMatrix[2][1] = sinTheta * cosPhi;
+
+ rotMatrix[0][2] = cosTheta * sinPhi * cosPsi + sinTheta * sinPsi;
+ rotMatrix[1][2] = cosTheta * sinPhi * sinPsi - sinTheta * cosPsi;
+ rotMatrix[2][2] = cosTheta * cosPhi;
+
+ /*--- loop through the points on the donor boundary ---*/
+ for (iPoint = 0; iPoint < nPerPointAll_donor; iPoint++) {
+ for (iDim = 0; iDim < nDim; iDim++) {
+ Coord_i[iDim] = Buffer_Recv_Coord_donor[iPoint * nDim + iDim];
+ }
+
+ /*--- Get the position vector from rotation center to point. ---*/
+ dx = Coord_i[0] - center[0];
+ dy = Coord_i[1] - center[1];
+ if (nDim == 3)
+ dz = Coord_i[2] - center[2];
+ else
+ dz = 0.0;
+
+ /*--- Compute transformed point coordinates. ---*/
+ rotCoord[0] = (rotMatrix[0][0] * dx + rotMatrix[0][1] * dy + rotMatrix[0][2] * dz + translation[0]);
+ rotCoord[1] = (rotMatrix[1][0] * dx + rotMatrix[1][1] * dy + rotMatrix[1][2] * dz + translation[1]);
+ rotCoord[2] = (rotMatrix[2][0] * dx + rotMatrix[2][1] * dy + rotMatrix[2][2] * dz + translation[2]);
+
+ /*--- Check if the point lies on the axis of rotation. If it does,
+ the rotated coordinate and the original coordinate are the same. ---*/
+
+ mindist = 1E6;
+ pPoint = 0;
+
+ for (jPoint = 0; jPoint < nPerPointAll_target; jPoint++) {
+ /*--- Compute the distance between the candidate periodic
+ point and the transformed coordinates of the owned point. ---*/
+ dist = 0.0;
+ for (iDim = 0; iDim < nDim; iDim++) {
+ Coord_j[iDim] = Buffer_Recv_Coord_target[jPoint * nDim + iDim];
+ dist += pow(Coord_j[iDim] - rotCoord[iDim], 2.0);
+ }
+ dist = sqrt(dist);
+
+ if (dist < mindist) {
+ /*--- We have found an intermediate match. Store the
+ data for this point before continuing the search. ---*/
+ mindist = dist;
+ pPoint = jPoint;
+ }
+ }
+
+ if (mindist > epsilon) {
+ cout << " Bad match for point " << iPoint << ".\tNearest";
+ cout << " donor distance: " << scientific << mindist << " Matched pPoint " << pPoint << "." << endl;
+ }
+
+ /*--- Store the global index for the best matched point. ---*/
+ PerPointsMatched_donor[iPoint] = PerPoints_target[pPoint];
+ }
+ }
+ }
+ }
+ }
+
+ if (rank == MASTER_NODE) {
+ cout << "[Multi-Passages Check] "
+ << " Matched " << nPerPointAll_target << endl;
+ cout << "[Multi-Passages Check] "
+ << " Match pair built. " << endl;
+ }
+
+ if (rank == MASTER_NODE) {
+ for (iPoint = 0; iPoint < PerPoints_donor.size(); iPoint++) {
+ PerPointsGlbIndx_donor[iPoint] = PerPoints_donor[iPoint];
+ for (iDim = 0; iDim < nDim; iDim++) {
+ PerPointsCoord_donor[iPoint * nDim + iDim] = Buffer_Recv_Coord_donor[iPoint * nDim + iDim];
+ }
+ }
+ }
+
+ /*--- Broadcast the final matched vector. ---*/
+ SU2_MPI::Bcast(PerPointsMatched_donor, nPerPointAll_donor, MPI_UNSIGNED_LONG, MASTER_NODE, SU2_MPI::GetComm());
+ SU2_MPI::Bcast(PerPointsGlbIndx_donor, nPerPointAll_donor, MPI_UNSIGNED_LONG, MASTER_NODE, SU2_MPI::GetComm());
+ SU2_MPI::Bcast(PerPointsCoord_donor, nPerPointAll_donor * nDim, MPI_DOUBLE, MASTER_NODE, SU2_MPI::GetComm());
+ for (iPoint = 0; iPoint < nPerPointAll_donor; iPoint++) {
+ PerPoint_Match_Map[PerPointsGlbIndx_donor[iPoint]] = PerPointsMatched_donor[iPoint];
+ Donor_Indx_Map[PerPointsGlbIndx_donor[iPoint]] = iPoint;
+ }
+
+ delete[] PerPointsOnRank_donor;
+ delete[] PerPointsOnRank_target;
+ delete[] Buffer_Send_Coord_donor;
+ delete[] Buffer_Recv_Coord_donor;
+ delete[] Buffer_Send_Coord_target;
+ delete[] Buffer_Recv_Coord_target;
+
+ // delete[] PerPointsMatched_donor;
+ // delete[] PerPointsGlbIndx_donor;
+ // decltype(PerPoint_Match_Map)().swap(PerPoint_Match_Map);
+ }
+
+ /*--- broadcast the total number of passages for full annulus ---*/
+ SU2_MPI::Bcast(&nPsgs_FullAnnu, 1, MPI_INT, MASTER_NODE, SU2_MPI::GetComm());
+ SU2_MPI::Bcast(&PitchAngle, 1, MPI_DOUBLE, MASTER_NODE, SU2_MPI::GetComm());
+
+ if (nPassages == nPsgs_FullAnnu) Fullannus = true;
+
+ delete[] displs_donor;
+ delete[] displs_target;
+ // decltype(PerPoints_donor)().swap(PerPoints_donor);
+ // decltype(PerPoints_target)().swap(PerPoints_target);
+}
+
+void CPhysicalGeometry::ReSetPeriodicAngle_Z(CConfig* config) {
+ unsigned short iMarker;
+ unsigned short val_periodic, nPeriodic;
+ const su2double *angles_donor, *angles_target;
+ string Donor_Tag, Target_Tag;
+
+ /*--- Evaluate the number of periodic boundary conditions ---*/
+ nPeriodic = config->GetnMarker_Periodic();
+
+ char* donor_buf = new char[MAX_STRING_SIZE];
+ char* target_buf = new char[MAX_STRING_SIZE];
+
+ /*--- First we find the string tag of 2 periodic boundaries on master core---*/
+ if (rank == MASTER_NODE) {
+ for (unsigned short iPeriodic = 1; iPeriodic <= nPeriodic / 2; iPeriodic++) {
+ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+ if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
+ val_periodic = config->GetMarker_All_PerBound(iMarker);
+ if (iPeriodic == val_periodic) {
+ Donor_Tag = config->GetMarker_All_TagBound(iMarker);
+ SPRINTF(&donor_buf[0], "%s", Donor_Tag.c_str());
+ }
+ if (iPeriodic == val_periodic + nPeriodic / 2) {
+ Target_Tag = config->GetMarker_All_TagBound(iMarker);
+ SPRINTF(&target_buf[0], "%s", Target_Tag.c_str());
+ }
+ }
+ }
+ }
+ }
+ cout << "rank " << rank << " Before Bcast Donor_Tag " << Donor_Tag << endl;
+ /*--- Then we broadcast this tag info to all cores ---*/
+ SU2_MPI::Bcast(donor_buf, MAX_STRING_SIZE, MPI_CHAR, MASTER_NODE, SU2_MPI::GetComm());
+ SU2_MPI::Bcast(target_buf, MAX_STRING_SIZE, MPI_CHAR, MASTER_NODE, SU2_MPI::GetComm());
+ cout << "rank " << rank << " After Bcast Donor_Tag " << Donor_Tag << endl;
+ // angles_donor = config->GetPeriodicRotAngles(Donor_Tag);
+ // angles_target = config->GetPeriodicRotAngles(Target_Tag);
+ /*--- Store the new pitch angle for user defined number of passages ---*/
+ /*--- Note that we assume Z-axis to be the rotating axis ---*/
+ // config->SetPeriodicRotAngle_Z(Donor_Tag, 2, angles_donor[2]*nPassages);
+ // config->SetPeriodicRotAngle_Z(Target_Tag, 2, angles_target[2]*nPassages);
+}
+
+/*--- duplicate the points to form full annulus grids ---*/
+void CPhysicalGeometry::DuplicatePoints(CConfig* config, CMeshReaderFVM* mesh) {
+ /*--- We need to rebuild the linear partition and redistribute points
+ on each node accordingly. Casue it is called later e.g. in DistributeColoring() ---*/
+
+ /*--- first we remove one marker(donor) of the periodic pair
+ and reorder the points on each node seperately ---*/
+
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Check] "
+ << " Remove one periodic boundary. Renumber points" << endl;
+
+ unsigned long nPoint_RemPer = nPoint - nPerPointOnRank_donor[rank];
+ CLinearPartitioner pointPartitioner(mesh->GetNumberOfGlobalPoints(), 0);
+ unsigned long GlobalIndex_Old = pointPartitioner.GetFirstIndexOnRank(rank);
+ // unsigned long GlobalIndex_New = GlobalIndex_Old;
+ /*
+ for (int iRank = 0; iRank < rank; iRank++){
+ GlobalIndex_New -= nPerPointOnRank_donor[iRank];
+ }*/
+
+ /*--- Get the linearly partitioned coordinates from the mesh object. ---*/
+ const auto& gridCoords = mesh->GetLocalPointCoordinates();
+
+ su2double* Coord_RemPer = new su2double[nPoint_RemPer * nDim];
+ unsigned long* OldGlobalIndex_RemPer = new unsigned long[nPoint_RemPer];
+
+ /*--- remove the donor points and renumber the point global index on each node ---*/
+ int iPoint_RemPer = 0;
+ unsigned long PerPointsStart = 0;
+
+ for (int iRank = 0; iRank < rank; iRank++) {
+ PerPointsStart += nPerPointOnRank_donor[iRank];
+ }
+
+ for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
+ bool PerPoint_Found = false;
+
+ for (unsigned long jPoint = 0; jPoint < nPerPointOnRank_donor[rank]; jPoint++) {
+ if (GlobalIndex_Old == PerPointsGlbIndx_donor[PerPointsStart + jPoint]) {
+ PerPoint_Found = true;
+ break;
+ }
+ }
+
+ if (!PerPoint_Found) {
+ for (unsigned short iDim = 0; iDim < nDim; ++iDim) {
+ Coord_RemPer[iPoint_RemPer * nDim + iDim] = gridCoords[iDim][iPoint];
+ }
+ OldGlobalIndex_RemPer[iPoint_RemPer] = GlobalIndex_Old;
+ //++GlobalIndex_New;
+ ++iPoint_RemPer;
+ }
+
+ ++GlobalIndex_Old;
+ }
+
+ /*--- We then collect all the points of one passage after removing
+ the donor periodic points, and distribute the vector to all nodes ---*/
+
+ /*--- prepare recvcounts and displs for mpi communication ---*/
+ int* displs_idx = new int[size];
+ int* displs_coord = new int[size];
+ for (int iRank = 0; iRank < size; iRank++) {
+ displs_idx[iRank] = 0;
+ displs_coord[iRank] = 0;
+ }
+ int* nPoint_RemPer_OnRank = new int[size];
+ int* CoordCount_RemPer_OnRank = new int[size];
+
+ for (int iRank = 0; iRank < size; iRank++) {
+ nPoint_RemPer_OnRank[iRank] = pointPartitioner.GetSizeOnRank(iRank) - nPerPointOnRank_donor[iRank];
+ CoordCount_RemPer_OnRank[iRank] = nPoint_RemPer_OnRank[iRank] * nDim;
+ }
+
+ for (int iRank = 0; iRank < size; iRank++) {
+ if (iRank == 0) {
+ displs_idx[0] = 0;
+ displs_coord[0] = 0;
+ } else {
+ displs_idx[iRank] = displs_idx[iRank - 1] + nPoint_RemPer_OnRank[iRank - 1];
+ displs_coord[iRank] = displs_coord[iRank - 1] + CoordCount_RemPer_OnRank[iRank - 1];
+ }
+ }
+
+ nPsgPoint_RemPer = mesh->GetNumberOfGlobalPoints() - nPerPointAll_donor;
+
+ /*--- prepare the recvbuf on all nodes to store the point information of one passage ---*/
+ PsgPointsOldGlbIndx_RemPer = new unsigned long[nPsgPoint_RemPer];
+ PsgPointsCoord_RemPer = new su2double[nPsgPoint_RemPer * nDim];
+
+ /*--- Gather the coordinates of all periodic points. ---*/
+
+ SU2_MPI::Allgatherv(OldGlobalIndex_RemPer, nPoint_RemPer, MPI_UNSIGNED_LONG, PsgPointsOldGlbIndx_RemPer,
+ nPoint_RemPer_OnRank, displs_idx, MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
+ SU2_MPI::Allgatherv(Coord_RemPer, nPoint_RemPer * nDim, MPI_DOUBLE, PsgPointsCoord_RemPer, CoordCount_RemPer_OnRank,
+ displs_coord, MPI_DOUBLE, SU2_MPI::GetComm());
+
+ for (unsigned long iPoint = 0; iPoint < nPoint_RemPer; iPoint++) {
+ // cout<<"rank "<SetGlobalIndex(iPoint, GlobalIndex_FullAnnu);
+
+ /*--- find which passage this iPoint locates ---*/
+ int nPsg = floor(GlobalIndex_FullAnnu / nPsgPoint_RemPer);
+
+ /*--- find the index inside one passage ---*/
+ unsigned long IndxInPsg = int(GlobalIndex_FullAnnu % nPsgPoint_RemPer);
+
+ /*--- coordinates before rotating ---*/
+ for (unsigned short iDim = 0; iDim < nDim; ++iDim) oriCoord[iDim] = PsgPointsCoord_RemPer[IndxInPsg * nDim + iDim];
+
+ /*--- Treatment for non-full-annulus case ---*/
+
+ /*--- For non-full-annulus case, an additional periodic
+ boundary shall be added at the last duplication, since
+ we removed one periodic boundary for each passage ---*/
+
+ /*--- Calculate the start index for the additonal periodic boundary point on this rank ---*/
+ if (!FirstAddPerPoint && (nPsg == nPassages)) {
+ FirstAddPerPoint = true;
+ IndxInPerDonor = nPerPointAll_donor - (nPoint - iPoint);
+ }
+ /*--- Retrieve the coordinates for points on the additional boundary ---*/
+ if (nPsg == nPassages) {
+ for (unsigned short iDim = 0; iDim < nDim; ++iDim)
+ oriCoord[iDim] = PerPointsCoord_donor[IndxInPerDonor * nDim + iDim];
+
+ // For points on the additional boundary, the rotating angle is PitchAngle*(nPssg-1)
+ nPsg--;
+
+ IndxInPerDonor++;
+ }
+
+ /*--- Store angles separately for clarity. Compute sines/cosines. ---*/
+ Theta = 0;
+ Phi = 0;
+ Psi = PitchAngle * nPsg;
+ cosTheta = cos(Theta);
+ cosPhi = cos(Phi);
+ cosPsi = cos(Psi);
+ sinTheta = sin(Theta);
+ sinPhi = sin(Phi);
+ sinPsi = sin(Psi);
+
+ /*--- Compute the rotation matrix. Note that the implicit
+ ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+ rotMatrix[0][0] = cosPhi * cosPsi;
+ rotMatrix[1][0] = cosPhi * sinPsi;
+ rotMatrix[2][0] = -sinPhi;
+
+ rotMatrix[0][1] = sinTheta * sinPhi * cosPsi - cosTheta * sinPsi;
+ rotMatrix[1][1] = sinTheta * sinPhi * sinPsi + cosTheta * cosPsi;
+ rotMatrix[2][1] = sinTheta * cosPhi;
+
+ rotMatrix[0][2] = cosTheta * sinPhi * cosPsi + sinTheta * sinPsi;
+ rotMatrix[1][2] = cosTheta * sinPhi * sinPsi - sinTheta * cosPsi;
+ rotMatrix[2][2] = cosTheta * cosPhi;
+
+ /*--- Compute transformed point coordinates. ---*/
+ rotCoord[0] = (rotMatrix[0][0] * oriCoord[0] + rotMatrix[0][1] * oriCoord[1] + rotMatrix[0][2] * oriCoord[2]);
+ rotCoord[1] = (rotMatrix[1][0] * oriCoord[0] + rotMatrix[1][1] * oriCoord[1] + rotMatrix[1][2] * oriCoord[2]);
+ rotCoord[2] = (rotMatrix[2][0] * oriCoord[0] + rotMatrix[2][1] * oriCoord[1] + rotMatrix[2][2] * oriCoord[2]);
+
+ /*--- Store new coordinates after rotating. ---*/
+ for (unsigned short iDim = 0; iDim < nDim; ++iDim) {
+ nodes->SetCoord(iPoint, iDim, rotCoord[iDim]);
+ LocalPointsCoord_FullAnnu[iPoint * nDim + iDim] = rotCoord[iDim];
+ }
+
+ GlobalIndex_FullAnnu++;
+ }
+
+ /*---last step: we make a map from the old index to new index
+ in one passage including the removed periodic points ---*/
+ for (unsigned long iPoint = 0; iPoint < nPsgPoint_RemPer; iPoint++) {
+ Psg_Old2New_Map[PsgPointsOldGlbIndx_RemPer[iPoint]] = iPoint;
+ Psg_Old2New_Map2[PsgPointsOldGlbIndx_RemPer[iPoint]] = iPoint;
+ }
+
+ for (unsigned long iPoint = 0; iPoint < nPerPointAll_donor; iPoint++) {
+ unsigned long newGlbIndx_donor = Psg_Old2New_Map[PerPointsMatched_donor[iPoint]] + nPsgPoint_RemPer;
+ Psg_Old2New_Map[PerPointsGlbIndx_donor[iPoint]] = newGlbIndx_donor;
+ Psg_Old2New_Map2[PerPointsGlbIndx_donor[iPoint]] = -1;
+ }
+
+ delete[] OldGlobalIndex_RemPer;
+ delete[] Coord_RemPer;
+ delete[] displs_idx;
+ delete[] displs_coord;
+ delete[] nPoint_RemPer_OnRank;
+ delete[] CoordCount_RemPer_OnRank;
+}
+
+/*--- duplicate the elements to form full annulus grids ---*/
+void CPhysicalGeometry::DuplicateVolumeElems(CConfig* config, CMeshReaderFVM* mesh) {
+ /*--- First, collect all elements info and distribute to all nodes---*/
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Check] "
+ << " Gather volume element info" << endl;
+ /*--- Gather the nElem info and distribute to all nodes ---*/
+ int* nElem_OnRank = new int[size];
+ int* displs = new int[size];
+ int* recvcounts = new int[size];
+ for (int iRank = 0; iRank < size; iRank++) {
+ displs[iRank] = iRank;
+ recvcounts[iRank] = 1;
+ }
+ int nElem_local = mesh->GetNumberOfLocalElements();
+
+ SU2_MPI::Allgatherv(&nElem_local, 1, MPI_INT, nElem_OnRank, recvcounts, displs, MPI_INT, SU2_MPI::GetComm());
+
+ auto& PsgConnElems_local = mesh->GetLocalVolumeElementConnectivity();
+
+ unsigned long nPsgElemGlobal = mesh->GetNumberOfGlobalElements();
+
+ /*--- There will be repeated elems in the PsgConnElems_Glb container,
+ compute the total number first. .*/
+ unsigned long nElem_collected = 0;
+ for (int iRank = 0; iRank < size; iRank++) nElem_collected += nElem_OnRank[iRank];
+
+ unsigned long* PsgConnElems_Glb = new unsigned long[nElem_collected * SU2_CONN_SIZE];
+
+ /*--- Gather the connectivity info and distribute to all nodes ---*/
+ int* recvcountss = new int[size];
+ int* displss = new int[size];
+ displss[0] = 0;
+ for (int iRank = 1; iRank < size; iRank++) {
+ displss[iRank] = 0;
+ for (int jRank = 0; jRank < iRank; jRank++) {
+ displss[iRank] += nElem_OnRank[jRank] * SU2_CONN_SIZE;
+ }
+ }
+ for (int iRank = 0; iRank < size; iRank++) {
+ recvcountss[iRank] = nElem_OnRank[iRank] * SU2_CONN_SIZE;
+ }
+ SU2_MPI::Allgatherv(&PsgConnElems_local[0], nElem * SU2_CONN_SIZE, MPI_UNSIGNED_LONG, PsgConnElems_Glb, recvcountss,
+ displss, MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
+
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Check] "
+ << " Duplicate and Renumber new volume elements in parallel" << endl;
+
+ /*--- Now we have the gathered element info on all nodes.
+ Start duplicating and distributing full annulus elements to all nodes.
+ Note that there is no repeated elements after duplication.
+
+ Secifically, loop through all base elements in one passage,
+ then, for each base element, loop through its series due to duplications.
+ Decide if anyone in each element series contains a point owned by this node(rank)
+ with the help of linear partition. If so, store all these elements locally.---*/
+
+ numberOfLocalElements_FullAnnu = 0;
+
+ /* Get a partitioner to help with linear partitioning. */
+ CLinearPartitioner pointPartitioner_fullannulus(Global_nPoint, 0);
+ const unsigned long localpoint_start = pointPartitioner_fullannulus.GetFirstIndexOnRank(rank);
+ const unsigned long localpoint_end = pointPartitioner_fullannulus.GetLastIndexOnRank(rank);
+ vector element_checked;
+ unsigned long nPoint_Check = nPsgPoint_RemPer * nPassages;
+ int counter = 0;
+ /*---First we loop through elems in PsgConnElems_Glb
+ recompute the local number of elements and the element index---*/
+ for (unsigned long iElem = 0; iElem < nElem_collected; iElem++) {
+ // if (rank == MASTER_NODE)
+ // cout<<"[ Multi-Passages Check] "<<" iElem "<(PsgConnElems_Glb[iElem * SU2_CONN_SIZE + 1]);
+ auto connectivity = &PsgConnElems_Glb[iElem * SU2_CONN_SIZE + SU2_CONN_SKIP];
+ bool isOwned = false;
+ unsigned long Global_Index_Elem_FullAnnu, Indx_FullAnnu;
+ int iPsg = 0;
+
+ // Avoid repeatition
+ if (find(element_checked.begin(), element_checked.end(), Global_Index_Elem) != element_checked.end()) continue;
+
+ element_checked.push_back(Global_Index_Elem);
+
+ switch (vtk_type) {
+ case TRIANGLE:
+ /* Check whether any of the points reside in our linear partition. */
+ for (iPsg = 0; iPsg < nPassages; iPsg++) {
+ isOwned = false;
+ for (unsigned short i = 0; i < N_POINTS_TRIANGLE; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ if ((Indx_FullAnnu >= localpoint_start) && (Indx_FullAnnu < localpoint_end)) {
+ isOwned = true;
+ break;
+ }
+ }
+ /* If so, we need to store the element locally. */
+ if (isOwned) {
+ Global_Index_Elem_FullAnnu = Global_Index_Elem + nPsgElemGlobal * iPsg;
+ localVolumeElementConnectivity_FullAnnu.push_back(Global_Index_Elem_FullAnnu);
+ localVolumeElementConnectivity_FullAnnu.push_back(vtk_type);
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ // Indx_FullAnnu = Psg_Old2New_Map[PerPoint_Match_Map[connectivity[i]]];
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ localVolumeElementConnectivity_FullAnnu.push_back(Indx_FullAnnu);
+ }
+ numberOfLocalElements_FullAnnu++;
+ }
+ }
+
+ break;
+
+ case QUADRILATERAL:
+ /* Check whether any of the points reside in our linear partition. */
+ for (iPsg = 0; iPsg < nPassages; iPsg++) {
+ isOwned = false;
+ for (unsigned short i = 0; i < N_POINTS_QUADRILATERAL; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ if ((Indx_FullAnnu >= localpoint_start) && (Indx_FullAnnu < localpoint_end)) {
+ isOwned = true;
+ break;
+ }
+ }
+ /* If so, we need to store the element locally. */
+ if (isOwned) {
+ Global_Index_Elem_FullAnnu = Global_Index_Elem + nPsgElemGlobal * iPsg;
+ localVolumeElementConnectivity_FullAnnu.push_back(Global_Index_Elem_FullAnnu);
+ localVolumeElementConnectivity_FullAnnu.push_back(vtk_type);
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ // Indx_FullAnnu = Psg_Old2New_Map[PerPoint_Match_Map[connectivity[i]]];
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ localVolumeElementConnectivity_FullAnnu.push_back(Indx_FullAnnu);
+ }
+ numberOfLocalElements_FullAnnu++;
+ }
+ }
+
+ break;
+
+ case TETRAHEDRON:
+ /* Check whether any of the points reside in our linear partition. */
+ for (iPsg = 0; iPsg < nPassages; iPsg++) {
+ isOwned = false;
+ for (unsigned short i = 0; i < N_POINTS_TETRAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ if ((Indx_FullAnnu >= localpoint_start) && (Indx_FullAnnu < localpoint_end)) {
+ isOwned = true;
+ break;
+ }
+ }
+ /* If so, we need to store the element locally. */
+ if (isOwned) {
+ Global_Index_Elem_FullAnnu = Global_Index_Elem + nPsgElemGlobal * iPsg;
+ localVolumeElementConnectivity_FullAnnu.push_back(Global_Index_Elem_FullAnnu);
+ localVolumeElementConnectivity_FullAnnu.push_back(vtk_type);
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ // Indx_FullAnnu = Psg_Old2New_Map[PerPoint_Match_Map[connectivity[i]]];
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ localVolumeElementConnectivity_FullAnnu.push_back(Indx_FullAnnu);
+ }
+ numberOfLocalElements_FullAnnu++;
+ }
+ }
+
+ break;
+
+ case HEXAHEDRON:
+ /* Check whether any of the points reside in our linear partition. */
+ for (iPsg = 0; iPsg < nPassages; iPsg++) {
+ isOwned = false;
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ if ((Indx_FullAnnu >= localpoint_start) && (Indx_FullAnnu < localpoint_end)) {
+ isOwned = true;
+ break;
+ }
+ }
+ /* If so, we need to store the element locally. */
+ if (isOwned) {
+ Global_Index_Elem_FullAnnu = Global_Index_Elem + nPsgElemGlobal * iPsg;
+ localVolumeElementConnectivity_FullAnnu.push_back(Global_Index_Elem_FullAnnu);
+ localVolumeElementConnectivity_FullAnnu.push_back(vtk_type);
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ localVolumeElementConnectivity_FullAnnu.push_back(Indx_FullAnnu);
+ }
+ numberOfLocalElements_FullAnnu++;
+ }
+ }
+
+ break;
+
+ case PRISM:
+ /* Check whether any of the points reside in our linear partition. */
+ for (iPsg = 0; iPsg < nPassages; iPsg++) {
+ isOwned = false;
+ for (unsigned short i = 0; i < N_POINTS_PRISM; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ if ((Indx_FullAnnu >= localpoint_start) && (Indx_FullAnnu < localpoint_end)) {
+ isOwned = true;
+ break;
+ }
+ }
+ /* If so, we need to store the element locally. */
+ if (isOwned) {
+ Global_Index_Elem_FullAnnu = Global_Index_Elem + nPsgElemGlobal * iPsg;
+ localVolumeElementConnectivity_FullAnnu.push_back(Global_Index_Elem_FullAnnu);
+ localVolumeElementConnectivity_FullAnnu.push_back(vtk_type);
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ // Indx_FullAnnu = Psg_Old2New_Map[PerPoint_Match_Map[connectivity[i]]];
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ localVolumeElementConnectivity_FullAnnu.push_back(Indx_FullAnnu);
+ }
+ numberOfLocalElements_FullAnnu++;
+ }
+ }
+
+ break;
+
+ case PYRAMID:
+ /* Check whether any of the points reside in our linear partition. */
+ for (iPsg = 0; iPsg < nPassages; iPsg++) {
+ isOwned = false;
+ for (unsigned short i = 0; i < N_POINTS_PYRAMID; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ if ((Indx_FullAnnu >= localpoint_start) && (Indx_FullAnnu < localpoint_end)) {
+ isOwned = true;
+ break;
+ }
+ }
+ /* If so, we need to store the element locally. */
+ if (isOwned) {
+ Global_Index_Elem_FullAnnu = Global_Index_Elem + nPsgElemGlobal * iPsg;
+ localVolumeElementConnectivity_FullAnnu.push_back(Global_Index_Elem_FullAnnu);
+ localVolumeElementConnectivity_FullAnnu.push_back(vtk_type);
+ for (unsigned short i = 0; i < N_POINTS_HEXAHEDRON; i++) {
+ Indx_FullAnnu = Psg_Old2New_Map[connectivity[i]] + iPsg * nPsgPoint_RemPer;
+ // The removed periodic points on the last duplicated passage locates on the first passage
+ if (Indx_FullAnnu >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case we have the following
+ // Indx_FullAnnu = Psg_Old2New_Map[PerPoint_Match_Map[connectivity[i]]];
+ Indx_FullAnnu = Indx_FullAnnu - nPoint_Check;
+ } else {
+ // for non full annulus case we have the following
+ Indx_FullAnnu = Donor_Indx_Map[connectivity[i]] + nPoint_Check;
+ }
+ }
+ localVolumeElementConnectivity_FullAnnu.push_back(Indx_FullAnnu);
+ }
+ numberOfLocalElements_FullAnnu++;
+ }
+ }
+
+ break;
+
+ default:
+ cout << " rank " << rank << " iElem " << iElem << " VTK_Type " << vtk_type << endl;
+ SU2_MPI::Error("Element type not supported!", CURRENT_FUNCTION);
+ break;
+ }
+ }
+
+ /*--- Now we know the local number of elements.
+ We begin to store the local elements to SU2's data structure.
+ We keep in alignment with the LoadLinearlyPartitionedVolumeElements(). ---*/
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Check] "
+ << " Build new connectivities inside element." << endl;
+ nElem = numberOfLocalElements_FullAnnu;
+ Global_nElem = nPsgElemGlobal * nPassages;
+ Global_nElemDomain = nPsgElemGlobal * nPassages;
+ elem = new CPrimalGrid*[nElem]();
+ Global_to_Local_Elem.clear();
+
+ for (unsigned long iElem = 0; iElem < nElem; iElem++) {
+ const auto Global_Index_Elem = localVolumeElementConnectivity_FullAnnu[iElem * SU2_CONN_SIZE + 0];
+ Global_to_Local_Elem[Global_Index_Elem] = iElem;
+ const auto vtk_type = static_cast(localVolumeElementConnectivity_FullAnnu[iElem * SU2_CONN_SIZE + 1]);
+ auto connectivity = &localVolumeElementConnectivity_FullAnnu[iElem * SU2_CONN_SIZE + SU2_CONN_SKIP];
+
+ switch (vtk_type) {
+ case TRIANGLE:
+ elem[iElem] = new CTriangle(connectivity[0], connectivity[1], connectivity[2]);
+ nelem_triangle++;
+ break;
+
+ case QUADRILATERAL:
+ elem[iElem] = new CQuadrilateral(connectivity[0], connectivity[1], connectivity[2], connectivity[3]);
+ nelem_quad++;
+ break;
+
+ case TETRAHEDRON:
+ elem[iElem] = new CTetrahedron(connectivity[0], connectivity[1], connectivity[2], connectivity[3]);
+ nelem_tetra++;
+ break;
+
+ case HEXAHEDRON:
+ elem[iElem] = new CHexahedron(connectivity[0], connectivity[1], connectivity[2], connectivity[3],
+ connectivity[4], connectivity[5], connectivity[6], connectivity[7]);
+ nelem_hexa++;
+ break;
+
+ case PRISM:
+ elem[iElem] = new CPrism(connectivity[0], connectivity[1], connectivity[2], connectivity[3], connectivity[4],
+ connectivity[5]);
+ nelem_prism++;
+ break;
+
+ case PYRAMID:
+ elem[iElem] = new CPyramid(connectivity[0], connectivity[1], connectivity[2], connectivity[3], connectivity[4]);
+ nelem_pyramid++;
+ break;
+
+ default:
+ SU2_MPI::Error("Element type not supported!", CURRENT_FUNCTION);
+ break;
+ }
+ }
+
+ /*--- Reduce the global counts of all element types found in
+ the CGNS grid with all ranks. ---*/
+
+ auto reduce = [](unsigned long p, unsigned long& t) {
+ SU2_MPI::Allreduce(&p, &t, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+ };
+ reduce(nelem_triangle, Global_nelem_triangle);
+ reduce(nelem_quad, Global_nelem_quad);
+ reduce(nelem_hexa, Global_nelem_hexa);
+ reduce(nelem_tetra, Global_nelem_tetra);
+ reduce(nelem_prism, Global_nelem_prism);
+ reduce(nelem_pyramid, Global_nelem_pyramid);
+
+ decltype(element_checked)().swap(element_checked);
+ delete[] nElem_OnRank;
+ delete[] displs;
+ delete[] recvcounts;
+ delete[] recvcountss;
+ delete[] displss;
+ delete[] PsgConnElems_Glb;
+}
+
+/*--- duplicate the elements to form full annulus grids ---*/
+void CPhysicalGeometry::DuplicateSurfaceElems(CConfig* config, CMeshReaderFVM* mesh) {
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Check] "
+ << " Duplicate and Renumber new surface elements" << endl;
+
+ if (rank == MASTER_NODE)
+ cout << "[Multi-Passages Check] "
+ << " Build new connectivities inside surface element." << endl;
+
+ /*--- We align with LoadUnpartitionedSurfaceElements(),
+ the master node takes care of loading all markers and surface elements---*/
+ if (rank == MASTER_NODE) {
+ const vector& sectionNames = mesh->GetMarkerNames();
+ if (nPassages == nPsgs_FullAnnu) {
+ // remove the periodic pair
+ bound = new CPrimalGrid**[nMarker - 2];
+ nElem_Bound = new unsigned long[nMarker - 2];
+ /*--- We leave the nMarker_Max unchanged.
+ Redundant memory won't change the progress ---*/
+ Tag_to_Marker = new string[config->GetnMarker_Max() - 2];
+ } else {
+ bound = new CPrimalGrid**[nMarker];
+ nElem_Bound = new unsigned long[nMarker];
+ Tag_to_Marker = new string[config->GetnMarker_Max()];
+ }
+
+ int npe, vtk_type;
+ vector connectivity(N_POINTS_HEXAHEDRON);
+ unsigned long iElem = 0;
+ int iMarker_New = 0;
+ unsigned long nPoint_Check = nPsgPoint_RemPer * nPassages;
+ unsigned short val_periodic;
+
+ /*--- Loop over all markers ---*/
+ for (int iMarker = 0; iMarker < nMarker; iMarker++) {
+ // Skip the periodic pair for full annulus case
+ // if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY)
+ if ((config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) && (nPassages == nPsgs_FullAnnu)) continue;
+
+ // Identify donor and target periodic oundary for non full annulus case
+ bool IsDonor = false;
+ bool IsTarget = false;
+ if ((config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) && (nPassages < nPsgs_FullAnnu)) {
+ for (unsigned short iPeriodic = 1; iPeriodic <= config->GetnMarker_Periodic() / 2; iPeriodic++) {
+ val_periodic = config->GetMarker_All_PerBound(iMarker);
+ if (iPeriodic == val_periodic) IsDonor = true;
+ if (iPeriodic + config->GetnMarker_Periodic() / 2 == val_periodic) IsTarget = true;
+ }
+ }
+
+ /*--- Reinitialize some counter variables ---*/
+ nelem_edge_bound = 0;
+ nelem_triangle_bound = 0;
+ nelem_quad_bound = 0;
+ iElem = 0;
+
+ /*--- Get the string name for this marker. ---*/
+ string Marker_Tag = sectionNames[iMarker];
+
+ /* Get the marker info and surface connectivity from the mesh object. */
+ const unsigned long surfElems = mesh->GetNumberOfSurfaceElementsForMarker(iMarker);
+ const vector& connElems = mesh->GetSurfaceElementConnectivityForMarker(iMarker);
+ if (config->GetMarker_All_KindBC(iMarker) == PERIODIC_BOUNDARY) {
+ nElem_Bound[iMarker_New] = surfElems;
+ } else {
+ nElem_Bound[iMarker_New] = surfElems * nPassages;
+ }
+ /*--- Report the number and name of the marker to the console. ---*/
+
+ cout << "[Full Annulus] " << nElem_Bound[iMarker] << " boundary elements in index ";
+ cout << iMarker << " (Marker = " << Marker_Tag << ")." << endl;
+
+ /*--- Reinstantiate the list of elements in the data structure. ---*/
+ bound[iMarker_New] = new CPrimalGrid*[nElem_Bound[iMarker_New]];
+
+ for (int iPsg = 0; iPsg < nPassages; iPsg++) {
+ /*--- Skip donor periodic boundary until the last loop ---*/
+ if (IsDonor && iPsg != nPassages - 1) {
+ continue;
+ }
+ /*--- Skip target periodic boundary unless the first loop ---*/
+ if (IsTarget && iPsg != 0) {
+ continue;
+ }
+
+ for (unsigned long jElem = 0; jElem < surfElems; jElem++) {
+ vtk_type = (int)connElems[jElem * SU2_CONN_SIZE + 1];
+
+ /*--- Store the loop size more easily. ---*/
+ npe = (int)(SU2_CONN_SIZE - SU2_CONN_SKIP);
+
+ /*--- Store the nodes for this element more clearly. ---*/
+ for (int j = 0; j < npe; j++) {
+ unsigned long nn = jElem * SU2_CONN_SIZE + SU2_CONN_SKIP + j;
+ connectivity[j] = Psg_Old2New_Map[connElems[nn]] + iPsg * nPsgPoint_RemPer;
+
+ if (connectivity[j] >= nPoint_Check) {
+ if (nPassages == nPsgs_FullAnnu) {
+ // for full annulus case
+ // connectivity[j] = Psg_Old2New_Map[PerPoint_Match_Map[connElems[nn]]];
+ connectivity[j] = connectivity[j] - nPoint_Check;
+ } else {
+ // for non full annulus case
+ connectivity[j] = Donor_Indx_Map[connElems[nn]] + nPoint_Check;
+ }
+ }
+ }
+
+ /*--- Instantiate the boundary element object. ---*/
+ switch (vtk_type) {
+ case LINE:
+ bound[iMarker_New][iElem] = new CLine(connectivity[0], connectivity[1]);
+ iElem++;
+ nelem_edge_bound++;
+ break;
+ case TRIANGLE:
+ bound[iMarker_New][iElem] = new CTriangle(connectivity[0], connectivity[1], connectivity[2]);
+ iElem++;
+ nelem_triangle_bound++;
+ break;
+ case QUADRILATERAL:
+ bound[iMarker_New][iElem] =
+ new CQuadrilateral(connectivity[0], connectivity[1], connectivity[2], connectivity[3]);
+ iElem++;
+ nelem_quad_bound++;
+ break;
+ }
+ }
+ }
+
+ /*--- Update config file lists in order to store the boundary
+ information for this marker in the correct place. ---*/
+
+ Tag_to_Marker[config->GetMarker_CfgFile_TagBound(Marker_Tag)] = Marker_Tag;
+ config->SetMarker_All_TagBound(iMarker_New, Marker_Tag);
+ config->SetMarker_All_KindBC(iMarker_New, config->GetMarker_CfgFile_KindBC(Marker_Tag));
+ config->SetMarker_All_Monitoring(iMarker_New, config->GetMarker_CfgFile_Monitoring(Marker_Tag));
+ config->SetMarker_All_GeoEval(iMarker_New, config->GetMarker_CfgFile_GeoEval(Marker_Tag));
+ config->SetMarker_All_Designing(iMarker_New, config->GetMarker_CfgFile_Designing(Marker_Tag));
+ config->SetMarker_All_Plotting(iMarker_New, config->GetMarker_CfgFile_Plotting(Marker_Tag));
+ config->SetMarker_All_Analyze(iMarker_New, config->GetMarker_CfgFile_Analyze(Marker_Tag));
+ config->SetMarker_All_ZoneInterface(iMarker_New, config->GetMarker_CfgFile_ZoneInterface(Marker_Tag));
+ config->SetMarker_All_DV(iMarker_New, config->GetMarker_CfgFile_DV(Marker_Tag));
+ config->SetMarker_All_Moving(iMarker_New, config->GetMarker_CfgFile_Moving(Marker_Tag));
+ config->SetMarker_All_Deform_Mesh(iMarker_New, config->GetMarker_CfgFile_Deform_Mesh(Marker_Tag));
+ config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker_New,
+ config->GetMarker_CfgFile_Deform_Mesh_Sym_Plane(Marker_Tag));
+ config->SetMarker_All_Fluid_Load(iMarker_New, config->GetMarker_CfgFile_Fluid_Load(Marker_Tag));
+ config->SetMarker_All_PyCustom(iMarker_New, config->GetMarker_CfgFile_PyCustom(Marker_Tag));
+ config->SetMarker_All_PerBound(iMarker_New, config->GetMarker_CfgFile_PerBound(Marker_Tag));
+ config->SetMarker_All_SendRecv(iMarker_New, NONE);
+ config->SetMarker_All_Turbomachinery(iMarker_New, config->GetMarker_CfgFile_Turbomachinery(Marker_Tag));
+ config->SetMarker_All_TurbomachineryFlag(iMarker_New, config->GetMarker_CfgFile_TurbomachineryFlag(Marker_Tag));
+ config->SetMarker_All_MixingPlaneInterface(iMarker_New,
+ config->GetMarker_CfgFile_MixingPlaneInterface(Marker_Tag));
+
+ iMarker_New++;
+ }
+
+ nMarker = iMarker_New;
+ config->SetnMarker_All(nMarker);
+ }
+}
+
void CPhysicalGeometry::PrepareAdjacency(const CConfig* config) {
#ifdef HAVE_MPI
#ifdef HAVE_PARMETIS
@@ -6755,7 +8118,7 @@ void CPhysicalGeometry::MatchPeriodic(const CConfig* config, unsigned short val_
Theta = angles[0];
Phi = angles[1];
- Psi = angles[2];
+ Psi = angles[2] * nPassages;
cosTheta = cos(Theta);
cosPhi = cos(Phi);
cosPsi = cos(Psi);
diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
index eb3edb511dbd..89738f786e05 100644
--- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
+++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl
@@ -882,75 +882,176 @@ void CFVMFlowSolverBase::LoadRestart_impl(CGeometry **geometry, CSolver **
/*--- Load data from the restart into correct containers. ---*/
unsigned long counter = 0;
- for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) {
- /*--- Retrieve local index. If this node from the restart file lives
- on the current processor, we will load and instantiate the vars. ---*/
+ if (config->GetTurbo_MultiPsgs()){
+
+ long PsgIndx_RemPer;
+ unsigned long GlobalIndx_FullAnnu;
+ unsigned long nPsgPoint_RemPer = geometry[MESH_0]->GetGlobal_nPsgPoint_RemPer();
+ //bool Fullannulus = false;
+ //if (geometry[MESH_0]->GetnPassages() == geometry[MESH_0]->GetnPassages_FullAnnu())
+ //Fullannulus = true;
+
+ su2double Theta, Phi, Psi, cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi;
+ su2double oriVel[3] = {0.0, 0.0, 0.0};
+ su2double rotVel[3] = {0.0, 0.0, 0.0};
+ su2double rotMatrix[3][3] = {{1.0,0.0,0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}};
+ su2double PitchAngle = geometry[MESH_0]->GetPitchangle();
+
+
+ for (int iPsg = 0; iPsg < geometry[MESH_0]->GetnPassages(); iPsg++) {
+
+
+ /*--- Store angles separately for clarity. Compute sines/cosines. ---*/
+ Theta = 0; Phi = 0; Psi = PitchAngle*iPsg;
+ cosTheta = cos(Theta); cosPhi = cos(Phi); cosPsi = cos(Psi);
+ sinTheta = sin(Theta); sinPhi = sin(Phi); sinPsi = sin(Psi);
+
+ /*--- Compute the rotation matrix. Note that the implicit
+ ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+ rotMatrix[0][0] = cosPhi*cosPsi;
+ rotMatrix[1][0] = cosPhi*sinPsi;
+ rotMatrix[2][0] = -sinPhi;
+
+ rotMatrix[0][1] = sinTheta*sinPhi*cosPsi - cosTheta*sinPsi;
+ rotMatrix[1][1] = sinTheta*sinPhi*sinPsi + cosTheta*cosPsi;
+ rotMatrix[2][1] = sinTheta*cosPhi;
+
+ rotMatrix[0][2] = cosTheta*sinPhi*cosPsi + sinTheta*sinPsi;
+ rotMatrix[1][2] = cosTheta*sinPhi*sinPsi - sinTheta*cosPsi;
+ rotMatrix[2][2] = cosTheta*cosPhi;
+
+
+ for (unsigned long iPoint = 0; iPoint < geometry[MESH_0]->GetGlobal_nPoint_OldPsg(); iPoint++ ) {
+
+ PsgIndx_RemPer = geometry[MESH_0]->Get_PsgOld2New_Indx_2(iPoint);
+
+ /*--- Skip the removed periodic points ---*/
+ if (PsgIndx_RemPer == -1){
+ if (!geometry[MESH_0]->GetFullannus() && (iPsg == geometry[MESH_0]->GetnPassages()-1)){
+ PsgIndx_RemPer = geometry[MESH_0]->Get_Donor_Indx(iPoint) + nPsgPoint_RemPer;
+
+ }else{
+ continue;
+ }
+ }
+
+
+ GlobalIndx_FullAnnu = PsgIndx_RemPer + nPsgPoint_RemPer*iPsg;
+
+ /*--- Retrieve local index. If this node from the restart file lives
+ on the current processor, we will load and instantiate the vars. ---*/
+ auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(GlobalIndx_FullAnnu);
+
+ if (iPoint_Local > -1) {
+
+ /*--- We need to store this point's data, so jump to the correct
+ offset in the buffer of data from the restart file and load it. ---*/
+ auto index = counter*Restart_Vars[1] + skipVars;
+
+ /*--- We need to transferm the velocity before
+ loading the restart solution . ---*/
+
+ /*--- velocities before rotating ---*/
+ for (unsigned short iDim = 0; iDim < nDim; ++iDim)
+ oriVel[iDim] = Restart_Data[index+1+iDim];
+
+ /*--- Compute transformed velocities. ---*/
+ rotVel[0] = (rotMatrix[0][0]*oriVel[0] + rotMatrix[0][1]*oriVel[1] + rotMatrix[0][2]*oriVel[2]);
+ rotVel[1] = (rotMatrix[1][0]*oriVel[0] + rotMatrix[1][1]*oriVel[1] + rotMatrix[1][2]*oriVel[2]);
+ rotVel[2] = (rotMatrix[2][0]*oriVel[0] + rotMatrix[2][1]*oriVel[1] + rotMatrix[2][2]*oriVel[2]);
+
+ /*--- rewrite the velocity ---*/
+ for (unsigned short iDim = 0; iDim < nDim; ++iDim)
+ Restart_Data[index+1+iDim] = rotVel[iDim];
+
+
+ if (SolutionRestart == nullptr) {
+ for (auto iVar = 0; iVar < nVar_Restart; iVar++)
+ nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
+ }
+ else {
+ /*--- Used as buffer, allows defaults for nVar > nVar_Restart. ---*/
+ for (auto iVar = 0; iVar < nVar_Restart; iVar++)
+ SolutionRestart[iVar] = Restart_Data[index+iVar];
+ nodes->SetSolution(iPoint_Local, SolutionRestart);
+ }
+
+ counter++;
+ }
+ }
+ }
+ }else{
- const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
+ for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) {
- if (iPoint_Local > -1) {
+ /*--- Retrieve local index. If this node from the restart file lives
+ on the current processor, we will load and instantiate the vars. ---*/
- /*--- We need to store this point's data, so jump to the correct
- offset in the buffer of data from the restart file and load it. ---*/
+ const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
- auto index = counter * Restart_Vars[1] + skipVars;
+ if (iPoint_Local > -1) {
- if (SolutionRestart == nullptr) {
- for (auto iVar = 0u; iVar < nVar_Restart; iVar++)
- nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
- }
- else {
- /*--- Used as buffer, allows defaults for nVar > nVar_Restart. ---*/
- for (auto iVar = 0u; iVar < nVar_Restart; iVar++)
- SolutionRestart[iVar] = Restart_Data[index + iVar];
- nodes->SetSolution(iPoint_Local, SolutionRestart);
- }
+ /*--- We need to store this point's data, so jump to the correct
+ offset in the buffer of data from the restart file and load it. ---*/
- /*--- For dynamic meshes, read in and store the
- grid coordinates and grid velocities for each node. ---*/
+ auto index = counter * Restart_Vars[1] + skipVars;
- if (dynamic_grid && update_geo) {
+ if (SolutionRestart == nullptr) {
+ for (auto iVar = 0u; iVar < nVar_Restart; iVar++)
+ nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
+ }
+ else {
+ /*--- Used as buffer, allows defaults for nVar > nVar_Restart. ---*/
+ for (auto iVar = 0u; iVar < nVar_Restart; iVar++)
+ SolutionRestart[iVar] = Restart_Data[index + iVar];
+ nodes->SetSolution(iPoint_Local, SolutionRestart);
+ }
- /*--- Read in the next 2 or 3 variables which are the grid velocities ---*/
- /*--- If we are restarting the solution from a previously computed static calculation (no grid movement) ---*/
- /*--- the grid velocities are set to 0. This is useful for FSI computations ---*/
+ /*--- For dynamic meshes, read in and store the
+ grid coordinates and grid velocities for each node. ---*/
- /*--- Rewind the index to retrieve the Coords. ---*/
- index = counter * Restart_Vars[1];
- const auto* Coord = &Restart_Data[index];
-
- su2double GridVel[MAXNDIM] = {0.0};
- if (!steady_restart) {
- /*--- Move the index forward to get the grid velocities. ---*/
- index += skipVars + nVar_Restart + config->GetnTurbVar();
- for (auto iDim = 0u; iDim < nDim; iDim++) { GridVel[iDim] = Restart_Data[index+iDim]; }
- }
+ if (dynamic_grid && update_geo) {
- for (auto iDim = 0u; iDim < nDim; iDim++) {
- geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]);
- geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, GridVel[iDim]);
+ /*--- Read in the next 2 or 3 variables which are the grid velocities ---*/
+ /*--- If we are restarting the solution from a previously computed static calculation (no grid movement) ---*/
+ /*--- the grid velocities are set to 0. This is useful for FSI computations ---*/
+
+ /*--- Rewind the index to retrieve the Coords. ---*/
+ index = counter * Restart_Vars[1];
+ const auto* Coord = &Restart_Data[index];
+
+ su2double GridVel[MAXNDIM] = {0.0};
+ if (!steady_restart) {
+ /*--- Move the index forward to get the grid velocities. ---*/
+ index += skipVars + nVar_Restart + config->GetnTurbVar();
+ for (auto iDim = 0u; iDim < nDim; iDim++) { GridVel[iDim] = Restart_Data[index+iDim]; }
+ }
+
+ for (auto iDim = 0u; iDim < nDim; iDim++) {
+ geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]);
+ geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, GridVel[iDim]);
+ }
}
- }
- /*--- For static FSI problems, grid_movement is 0 but we need to read in and store the
- grid coordinates for each node (but not the grid velocities, as there are none). ---*/
+ /*--- For static FSI problems, grid_movement is 0 but we need to read in and store the
+ grid coordinates for each node (but not the grid velocities, as there are none). ---*/
- if (static_fsi && update_geo) {
- /*--- Rewind the index to retrieve the Coords. ---*/
- index = counter*Restart_Vars[1];
- const auto* Coord = &Restart_Data[index];
+ if (static_fsi && update_geo) {
+ /*--- Rewind the index to retrieve the Coords. ---*/
+ index = counter*Restart_Vars[1];
+ const auto* Coord = &Restart_Data[index];
- for (auto iDim = 0u; iDim < nDim; iDim++) {
- geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]);
+ for (auto iDim = 0u; iDim < nDim; iDim++) {
+ geometry[MESH_0]->nodes->SetCoord(iPoint_Local, iDim, Coord[iDim]);
+ }
}
- }
- /*--- Increment the overall counter for how many points have been loaded. ---*/
- counter++;
+ /*--- Increment the overall counter for how many points have been loaded. ---*/
+ counter++;
+ }
}
}
-
/*--- Detect a wrong solution file ---*/
if (counter != nPointDomain) {
diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp
index 1f2ca053b409..eaef5416d408 100644
--- a/SU2_CFD/src/solvers/CSolver.cpp
+++ b/SU2_CFD/src/solvers/CSolver.cpp
@@ -461,7 +461,7 @@ void CSolver::InitiatePeriodicComms(CGeometry *geometry,
su2double Theta = angles[0];
su2double Phi = angles[1];
- su2double Psi = angles[2];
+ su2double Psi = angles[2]*geometry->GetnPassages();
/*--- Compute the rotation matrix. Note that the implicit
ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
@@ -2858,37 +2858,92 @@ void CSolver::Read_SU2_Restart_ASCII(CGeometry *geometry, const CConfig *config,
Restart_Data.resize(Restart_Vars[1]*geometry->GetnPointDomain());
- /*--- Read all lines in the restart file and extract data. ---*/
+ //bool FullAnnu_Transform = 1;
+ long PsgIndx_RemPer;
+ unsigned long GlobalIndx_FullAnnu;
+ unsigned long nPsgPoint_RemPer = geometry->GetGlobal_nPsgPoint_RemPer();
+ //bool Fullannulus = false;
+
+ //if (geometry->GetnPassages() == geometry->GetnPassages_FullAnnu())
+ // Fullannulus = true;
+
+ if (config->GetTurbo_MultiPsgs()){
+
+ for (int iPsg = 0; iPsg < geometry->GetnPassages(); iPsg++) {
+
+ /*--- read restart file for nPassages times ---*/
+ restart_file.close();
+ restart_file.open(val_filename.data(), ios::in);
+ getline (restart_file, text_line);
+
+ for (unsigned long iPoint = 0; iPoint < geometry->GetGlobal_nPoint_OldPsg(); iPoint++ ) {
+
+ getline (restart_file, text_line);
+ vector point_line = PrintingToolbox::split(text_line, delimiter);
+
+ PsgIndx_RemPer = geometry->Get_PsgOld2New_Indx_2(iPoint);
+
+ /*--- Skip the removed periodic points ---*/
+ if (PsgIndx_RemPer == -1){
+ if (!geometry->GetFullannus() && (iPsg == geometry->GetnPassages()-1)){
+ PsgIndx_RemPer = geometry->Get_Donor_Indx(iPoint) + nPsgPoint_RemPer;
+
+ }else{
+ continue;
+ }
+ }
- for (iPoint_Global = 0; iPoint_Global < geometry->GetGlobal_nPointDomain(); iPoint_Global++) {
- if (!getline (restart_file, text_line)) break;
+ GlobalIndx_FullAnnu = PsgIndx_RemPer + nPsgPoint_RemPer*iPsg;
+
+ /*--- Retrieve local index. If this node from the restart file lives
+ on the current processor, we will load and instantiate the vars. ---*/
+ iPoint_Local = geometry->GetGlobal_to_Local_Point(GlobalIndx_FullAnnu);
- /*--- Retrieve local index. If this node from the restart file lives
- on the current processor, we will load and instantiate the vars. ---*/
+ if (iPoint_Local > -1) {
- iPoint_Local = geometry->GetGlobal_to_Local_Point(iPoint_Global);
+ for (iVar = 0; iVar < Restart_Vars[1]; iVar++)
+ Restart_Data[counter*Restart_Vars[1] + iVar] = SU2_TYPE::GetValue(PrintingToolbox::stod(point_line[iVar+1]));
- if (iPoint_Local > -1) {
+ counter++;
+ }
+ }
- vector point_line = PrintingToolbox::split(text_line, delimiter);
+ }
+ }else{
- /*--- Store the solution (starting with node coordinates) --*/
+ /*--- Read all lines in the restart file and extract data. ---*/
- for (iVar = 0; iVar < Restart_Vars[1]; iVar++)
- Restart_Data[counter*Restart_Vars[1] + iVar] = SU2_TYPE::GetValue(PrintingToolbox::stod(point_line[iVar+1]));
+ for (iPoint_Global = 0; iPoint_Global < geometry->GetGlobal_nPointDomain(); iPoint_Global++) {
- /*--- Increment our local point counter. ---*/
+ if (!getline (restart_file, text_line)) break;
- counter++;
+ /*--- Retrieve local index. If this node from the restart file lives
+ on the current processor, we will load and instantiate the vars. ---*/
- }
- }
+ iPoint_Local = geometry->GetGlobal_to_Local_Point(iPoint_Global);
- if (iPoint_Global != geometry->GetGlobal_nPointDomain())
- SU2_MPI::Error("The solution file does not match the mesh, currently only binary files can be interpolated.",
- CURRENT_FUNCTION);
+ if (iPoint_Local > -1) {
+
+ vector point_line = PrintingToolbox::split(text_line, delimiter);
+
+ /*--- Store the solution (starting with node coordinates) --*/
+ for (iVar = 0; iVar < Restart_Vars[1]; iVar++)
+ Restart_Data[counter*Restart_Vars[1] + iVar] = SU2_TYPE::GetValue(PrintingToolbox::stod(point_line[iVar+1]));
+
+ /*--- Increment our local point counter. ---*/
+
+ counter++;
+
+ }
+ }
+
+
+ if (iPoint_Global != geometry->GetGlobal_nPointDomain())
+ SU2_MPI::Error("The solution file does not match the mesh, currently only binary files can be interpolated.",
+ CURRENT_FUNCTION);
+ }
}
void CSolver::Read_SU2_Restart_Binary(CGeometry *geometry, const CConfig *config, string val_filename) {
diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp
index d0744a7e2236..05395b1d507b 100644
--- a/SU2_CFD/src/solvers/CTurbSolver.cpp
+++ b/SU2_CFD/src/solvers/CTurbSolver.cpp
@@ -136,24 +136,70 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig*
/*--- Load data from the restart into correct containers. ---*/
unsigned long counter = 0;
- for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) {
- /*--- Retrieve local index. If this node from the restart file lives
- on the current processor, we will load and instantiate the vars. ---*/
- const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
+ unsigned long nPsgPoint_RemPer = geometry[MESH_0]->GetGlobal_nPsgPoint_RemPer();
+ //bool Fullannulus = false;
+ //if (geometry[MESH_0]->GetnPassages() == geometry[MESH_0]->GetnPassages_FullAnnu())
+ //Fullannulus = true;
+
+ if (config->GetTurbo_MultiPsgs()){
+ for (int iPsg = 0; iPsg < geometry[MESH_0]->GetnPassages(); iPsg++) {
+ for (unsigned long iPoint = 0; iPoint < geometry[MESH_0]->GetGlobal_nPoint_OldPsg(); iPoint++ ) {
+
+ long PsgIndx_RemPer = geometry[MESH_0]->Get_PsgOld2New_Indx_2(iPoint);
+
+ /*--- Skip the removed periodic points ---*/
+ if (PsgIndx_RemPer == -1){
+ if (!geometry[MESH_0]->GetFullannus() && (iPsg == geometry[MESH_0]->GetnPassages()-1)){
+ PsgIndx_RemPer = geometry[MESH_0]->Get_Donor_Indx(iPoint) + nPsgPoint_RemPer;
+
+ }else{
+ continue;
+ }
+ }
+
+ unsigned long GlobalIndx_FullAnnu = PsgIndx_RemPer + nPsgPoint_RemPer*iPsg;
+
+ /*--- Retrieve local index. If this node from the restart file lives
+ on the current processor, we will load and instantiate the vars. ---*/
+ auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(GlobalIndx_FullAnnu);
+
+ if (iPoint_Local > -1) {
+
+ /*--- We need to store this point's data, so jump to the correct
+ offset in the buffer of data from the restart file and load it. ---*/
+ const auto index = counter*Restart_Vars[1] + skipVars;
+ for (auto iVar = 0; iVar < nVar; ++iVar)
+ nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index+iVar]);
+
+ /*--- Increment the overall counter for how many points have been loaded. ---*/
+ counter++;
+ }
+
+ //cout<<"iPoint_Local "< -1) {
- /*--- We need to store this point's data, so jump to the correct
- offset in the buffer of data from the restart file and load it. ---*/
+ for (auto iPoint_Global = 0ul; iPoint_Global < geometry[MESH_0]->GetGlobal_nPointDomain(); iPoint_Global++) {
+ /*--- Retrieve local index. If this node from the restart file lives
+ on the current processor, we will load and instantiate the vars. ---*/
- const auto index = counter * Restart_Vars[1] + skipVars;
- for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]);
+ const auto iPoint_Local = geometry[MESH_0]->GetGlobal_to_Local_Point(iPoint_Global);
- /*--- Increment the overall counter for how many points have been loaded. ---*/
- counter++;
+ if (iPoint_Local > -1) {
+ /*--- We need to store this point's data, so jump to the correct
+ offset in the buffer of data from the restart file and load it. ---*/
+
+ const auto index = counter * Restart_Vars[1] + skipVars;
+ for (auto iVar = 0u; iVar < nVar; iVar++) nodes->SetSolution(iPoint_Local, iVar, Restart_Data[index + iVar]);
+
+ /*--- Increment the overall counter for how many points have been loaded. ---*/
+ counter++;
+ }
}
}
-
/*--- Detect a wrong solution file ---*/
if (counter != nPointDomain) {