diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 992a0d6e30a..bfe929f6871 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 paraviewParallelOutput -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 paraviewParallelOutput -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 paraviewParallelOutput -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 paraviewParallelOutput -s ${{matrix.testscript}} -a "--asan" - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2-asan:261011-1107 with: diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1af43985c12..02c5690f204 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2408,6 +2408,8 @@ enum class OUTPUT_TYPE { PARAVIEW_XML, /*!< \brief Paraview XML with binary data format */ SURFACE_PARAVIEW_XML, /*!< \brief Surface Paraview XML with binary data format */ PARAVIEW_MULTIBLOCK, /*!< \brief Paraview XML Multiblock */ + PARAVIEW_PARALLEL, /*!< \brief Paraview parallel XML (.pvtu), one piece per rank. */ + SURFACE_PARAVIEW_PARALLEL, /*!< \brief Surface Paraview parallel XML (.pvtu), one piece per rank. */ CGNS, /*!< \brief CGNS format. */ SURFACE_CGNS, /*!< \brief CGNS format. */ STL_ASCII, /*!< \brief STL ASCII format for surface solution output. */ @@ -2427,6 +2429,8 @@ static const MapType Output_Map = { MakePair("PARAVIEW", OUTPUT_TYPE::PARAVIEW_XML) MakePair("SURFACE_PARAVIEW", OUTPUT_TYPE::SURFACE_PARAVIEW_XML) MakePair("PARAVIEW_MULTIBLOCK", OUTPUT_TYPE::PARAVIEW_MULTIBLOCK) + MakePair("PARAVIEW_PARALLEL", OUTPUT_TYPE::PARAVIEW_PARALLEL) + MakePair("SURFACE_PARAVIEW_PARALLEL", OUTPUT_TYPE::SURFACE_PARAVIEW_PARALLEL) MakePair("MESH", OUTPUT_TYPE::MESH) MakePair("MESH_BINARY", OUTPUT_TYPE::MESH_BINARY) MakePair("RESTART_ASCII", OUTPUT_TYPE::RESTART_ASCII) @@ -2449,6 +2453,8 @@ inline bool isParaview(OUTPUT_TYPE format) { case OUTPUT_TYPE::PARAVIEW_XML: case OUTPUT_TYPE::SURFACE_PARAVIEW_XML: case OUTPUT_TYPE::PARAVIEW_MULTIBLOCK: + case OUTPUT_TYPE::PARAVIEW_PARALLEL: + case OUTPUT_TYPE::SURFACE_PARAVIEW_PARALLEL: return true; default: return false; diff --git a/SU2_CFD/include/output/filewriter/CParaviewPVTUFileWriter.hpp b/SU2_CFD/include/output/filewriter/CParaviewPVTUFileWriter.hpp new file mode 100644 index 00000000000..c7230eb1384 --- /dev/null +++ b/SU2_CFD/include/output/filewriter/CParaviewPVTUFileWriter.hpp @@ -0,0 +1,105 @@ +/*! + * \file CParaviewPVTUFileWriter.hpp + * \brief Headers of the Paraview parallel XML file writer (.pvtu and one .vtu piece per rank). + * \author A. Rausa + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include + +#include "CFileWriter.hpp" + +/*! + * \class CParaviewPVTUFileWriter + * \brief Writes the data as a Paraview parallel unstructured grid: a .pvtu file that lists one .vtu + * piece per rank. Each rank writes its own piece, so no rank gathers data of other ranks and + * ParaView can read, filter and render the pieces in parallel (unlike a single .vtu file, + * which one process has to load as a whole). + */ +class CParaviewPVTUFileWriter final : public CFileWriter { + private: + /*! + * \brief Boolean storing whether we are on a big or little endian machine. + */ + bool bigEndian; + + /*! + * \brief True to write the coordinates and fields in double precision instead of single. + */ + bool doublePrecision; + + /*! + * \brief Points of the piece of this rank: those used by its elements, wherever they belong. + */ + vector pieceData; /*!< \brief Fields of the points of the piece, field by field. */ + vector usedPoints; /*!< \brief Global indices of the points of the piece, sorted. */ + + /*! + * \brief Collect the points used by the elements of this rank and get their data from the ranks that own + * them, so that the piece of this rank is self-contained. + * \returns The number of points of the piece. + */ + unsigned long CollectExternalPoints(); + + /*! + * \brief Write the piece of this rank. + * \param[in] val_filename - Name of the piece file. + * \param[in] nPointsPiece - Number of points of the piece. + * \param[in] connInt64 - Write the connectivity and offsets as Int64 instead of Int32. + */ + void WritePiece(const string& val_filename, unsigned long nPointsPiece, bool connInt64); + + /*! + * \brief Write the .pvtu file that lists the pieces, only the master node writes it. + * \param[in] val_filename - Name of the .pvtu file. + * \param[in] pieceNames - Names of the pieces, relative to the .pvtu file. + * \param[in] connInt64 - The connectivity and offsets of the pieces are Int64 instead of Int32. + */ + void WriteIndexFile(const string& val_filename, const vector& pieceNames, bool connInt64); + + public: + /*! + * \brief File extension + */ + const static string fileExt; + + /*! + * \brief Construct a file writer using field names and the data sorter. + * \param[in] valDataSorter - The parallel sorted data to write. + * \param[in] valDoublePrecision - Write the coordinates and fields in double precision. + */ + CParaviewPVTUFileWriter(CParallelDataSorter* valDataSorter, bool valDoublePrecision = false); + + /*! + * \brief Destructor + */ + ~CParaviewPVTUFileWriter() override; + + /*! + * \brief Write sorted data to the .pvtu file and its pieces. + * \param[in] val_filename - The name of the file. + */ + void WriteData(string val_filename) override; +}; diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index f2458d8c87c..7158ff0112c 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -51,6 +51,7 @@ su2_cfd_src += files(['output/COutputFactory.cpp', 'output/filewriter/CSU2BinaryFileWriter.cpp', 'output/filewriter/CParaviewXMLFileWriter.cpp', 'output/filewriter/CParaviewVTMFileWriter.cpp', + 'output/filewriter/CParaviewPVTUFileWriter.cpp', 'output/filewriter/CSU2MeshFileWriter.cpp', 'output/filewriter/CSU2MeshBinaryFileWriter.cpp', 'output/filewriter/CCGNSFileWriter.cpp', diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp index c9ca947546c..d458d671855 100644 --- a/SU2_CFD/src/output/COutput.cpp +++ b/SU2_CFD/src/output/COutput.cpp @@ -42,6 +42,7 @@ #include "../../include/output/filewriter/CSTLFileWriter.hpp" #include "../../include/output/filewriter/CParaviewBinaryFileWriter.hpp" #include "../../include/output/filewriter/CParaviewXMLFileWriter.hpp" +#include "../../include/output/filewriter/CParaviewPVTUFileWriter.hpp" #include "../../include/output/filewriter/CParaviewVTMFileWriter.hpp" #include "../../include/output/filewriter/CTecplotFileWriter.hpp" #include "../../include/output/filewriter/CTecplotBinaryFileWriter.hpp" @@ -560,6 +561,46 @@ void COutput::WriteToFile(CConfig *config, CGeometry *geometry, OUTPUT_TYPE form break; + case OUTPUT_TYPE::PARAVIEW_PARALLEL: + + extension = CParaviewPVTUFileWriter::fileExt; + + if (fileName.empty()) + fileName = config->GetFilename(volumeFilename, "", curTimeIter); + + if (!config->GetWrt_Volume_Overwrite()) + filename_iter = config->GetFilename_Iter(fileName, curInnerIter, curOuterIter); + + /*--- Load the connectivity without sorting it: each rank then writes the elements of its own + partition, which are next to each other, so its piece needs few points of other ranks. ---*/ + + volumeDataSorter->SortConnectivity(config, geometry, false); + + LogOutputFiles("Paraview parallel"); + fileWriter = new CParaviewPVTUFileWriter(volumeDataSorter, config->GetVolume_Output_Double_Precision()); + + break; + + case OUTPUT_TYPE::SURFACE_PARAVIEW_PARALLEL: + + extension = CParaviewPVTUFileWriter::fileExt; + + if (fileName.empty()) + fileName = config->GetFilename(surfaceFilename, "", curTimeIter); + + if (!config->GetWrt_Surface_Overwrite()) + filename_iter = config->GetFilename_Iter(fileName, curInnerIter, curOuterIter); + + /*--- Load and sort the output data and connectivity. ---*/ + + surfaceDataSorter->SortConnectivity(config, geometry); + surfaceDataSorter->SortOutputData(); + + LogOutputFiles("Paraview parallel surface"); + fileWriter = new CParaviewPVTUFileWriter(surfaceDataSorter, config->GetVolume_Output_Double_Precision()); + + break; + case OUTPUT_TYPE::PARAVIEW_LEGACY_BINARY: extension = CParaviewBinaryFileWriter::fileExt; diff --git a/SU2_CFD/src/output/filewriter/CParaviewPVTUFileWriter.cpp b/SU2_CFD/src/output/filewriter/CParaviewPVTUFileWriter.cpp new file mode 100644 index 00000000000..aa7818df114 --- /dev/null +++ b/SU2_CFD/src/output/filewriter/CParaviewPVTUFileWriter.cpp @@ -0,0 +1,414 @@ +/*! + * \file CParaviewPVTUFileWriter.cpp + * \brief Filewriter class for the Paraview parallel XML format (.pvtu and one .vtu piece per rank). + * \author A. Rausa + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../../include/output/filewriter/CParaviewPVTUFileWriter.hpp" + +#include +#include +#include + +#if defined(_WIN32) +#include +#else +#include +#include +#include +#include +#endif + +namespace { + +constexpr unsigned short NCOORDS = 3; + +/*--- The element types in the order in which they are written. ---*/ +const std::array ElemTypes = {LINE, TRIANGLE, QUADRILATERAL, TETRAHEDRON, + HEXAHEDRON, PRISM, PYRAMID}; + +/*--- Write a data array in the "appended raw" format: the size in bytes of the array followed by its + values, converted to the type the file declares. Returns the number of bytes written. ---*/ +template +size_t WriteDataArray(FILE* file, const vector& values) { + const auto nBytes = static_cast(values.size() * sizeof(T)); + fwrite(&nBytes, sizeof(nBytes), 1, file); + + vector converted(values.begin(), values.end()); + fwrite(converted.data(), sizeof(T), converted.size(), file); + + return sizeof(nBytes) + nBytes; +} + +/*--- Open a new file for writing. On POSIX systems it is created with the permissions of the other output files + (owner read/write, group and others read), fopen would create it with 0666 before the umask. ---*/ +FILE* OpenFile(const string& name, const char* mode) { +#if defined(_WIN32) + return fopen(name.c_str(), mode); +#else + const int fd = open(name.c_str(), O_WRONLY | O_CREAT | O_TRUNC, S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH); + if (fd < 0) return nullptr; + FILE* file = fdopen(fd, mode); + if (!file) close(fd); + return file; +#endif +} + +/*--- Create a folder, if it does not exist yet. ---*/ +void MakeFolder(const string& name) { +#if defined(_WIN32) +#ifdef __MINGW32__ + mkdir(name.c_str()); +#else + _mkdir(name.c_str()); +#endif +#else + mkdir(name.c_str(), 0777); +#endif +} + +/*--- Size in bytes of a data array in the file, including the size written in front of it. ---*/ +size_t DataArraySize(size_t nValues, size_t typeSize) { return sizeof(uint64_t) + nValues * typeSize; } + +/*--- Name of a field without the component suffix of a vector, e.g. "Momentum_x" -> "Momentum". ---*/ +string FieldBaseName(string name) { + name.erase(remove(name.begin(), name.end(), '"'), name.end()); + return name.substr(0, name.size() - 2); +} + +} // namespace + +const string CParaviewPVTUFileWriter::fileExt = ".pvtu"; + +CParaviewPVTUFileWriter::CParaviewPVTUFileWriter(CParallelDataSorter* valDataSorter, bool valDoublePrecision) + : CFileWriter(valDataSorter, fileExt), doublePrecision(valDoublePrecision) { + /* Check for big endian. We have to swap bytes otherwise. */ + + bigEndian = false; + unsigned int i = 1; + char* c = (char*)&i; + bigEndian = *c == 0; +} + +CParaviewPVTUFileWriter::~CParaviewPVTUFileWriter() = default; + +void CParaviewPVTUFileWriter::WriteData(string val_filename) { + if (!dataSorter->GetConnectivitySorted()) { + SU2_MPI::Error("Connectivity must be sorted.", CURRENT_FUNCTION); + } + + if (bigEndian) { + SU2_MPI::Error("The Paraview parallel writer is only implemented for little endian machines.", CURRENT_FUNCTION); + } + + startTime = SU2_MPI::Wtime(); + + /*--- The pieces are written into a folder next to the .pvtu file, one file per rank. ---*/ + + const string folderName = val_filename; + + if (rank == MASTER_NODE) MakeFolder(folderName); + SU2_MPI::Barrier(SU2_MPI::GetComm()); + + const auto lastSlash = val_filename.find_last_of("/\\"); + const string baseName = (lastSlash == string::npos) ? val_filename : val_filename.substr(lastSlash + 1); + + vector pieceNames(size); + for (int i = 0; i < size; i++) pieceNames[i] = baseName + "/" + baseName + "_" + to_string(i) + ".vtu"; + + /*--- Make the data of this rank self-contained and write its piece. ---*/ + + const auto nPointsPiece = CollectExternalPoints(); + + /*--- All pieces must declare the same types, so the connectivity type is that of the largest piece. ---*/ + + unsigned long localConn = dataSorter->GetnConn(), maxConn = 0; + SU2_MPI::Allreduce(&localConn, &maxConn, 1, MPI_UNSIGNED_LONG, MPI_MAX, SU2_MPI::GetComm()); + const bool connInt64 = maxConn > static_cast(std::numeric_limits::max()); + + WritePiece(folderName + "/" + baseName + "_" + to_string(rank) + ".vtu", nPointsPiece, connInt64); + + /*--- The master node writes the file that lists the pieces. ---*/ + + WriteIndexFile(val_filename + fileExt, pieceNames, connInt64); + + stopTime = SU2_MPI::Wtime(); + usedTime = stopTime - startTime; + + su2double myFileSize = fileSize; + SU2_MPI::Allreduce(&myFileSize, &fileSize, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + bandwidth = fileSize / (1.0e6) / usedTime; +} + +unsigned long CParaviewPVTUFileWriter::CollectExternalPoints() { + const auto nFields = dataSorter->GetFieldNames().size(); + const auto firstPoint = dataSorter->GetnPointCumulative(rank); + + /*--- The piece holds exactly the points used by the elements of this rank, wherever they live. ---*/ + + usedPoints.clear(); + usedPoints.reserve(dataSorter->GetnConn()); + + for (auto type : ElemTypes) { + const auto nNodes = nPointsOfElementType(type); + for (unsigned long iElem = 0; iElem < dataSorter->GetnElem(type); iElem++) + for (unsigned short iNode = 0; iNode < nNodes; iNode++) + usedPoints.push_back(dataSorter->GetElemConnectivity(type, iElem, iNode) - 1); + } + + sort(usedPoints.begin(), usedPoints.end()); + usedPoints.erase(unique(usedPoints.begin(), usedPoints.end()), usedPoints.end()); + usedPoints.shrink_to_fit(); + + if (usedPoints.size() > static_cast(std::numeric_limits::max())) + SU2_MPI::Error("Too many points in one piece of the Paraview parallel output.", CURRENT_FUNCTION); + + /*--- Ask the rank that owns each point for its data. The points are sorted by their global index, and + the ranks own ranges of increasing global indices, so they are already grouped by rank. ---*/ + + vector nRequest(size, 0), nAnswer(size, 0), requestOffset(size + 1, 0), answerOffset(size + 1, 0); + + for (auto iPoint : usedPoints) nRequest[dataSorter->FindProcessor(iPoint)]++; + + SU2_MPI::Alltoall(nRequest.data(), 1, MPI_INT, nAnswer.data(), 1, MPI_INT, SU2_MPI::GetComm()); + + for (int i = 0; i < size; i++) { + requestOffset[i + 1] = requestOffset[i] + nRequest[i]; + answerOffset[i + 1] = answerOffset[i] + nAnswer[i]; + } + + vector answeredPoints(answerOffset[size]); + + SU2_MPI::Alltoallv(usedPoints.data(), nRequest.data(), requestOffset.data(), MPI_UNSIGNED_LONG, + answeredPoints.data(), nAnswer.data(), answerOffset.data(), MPI_UNSIGNED_LONG, + SU2_MPI::GetComm()); + + /*--- Answer with the data of the points this rank owns, one field at a time so that the MPI counts are numbers + of points. The data is passive, so the passive MPI wrapper is used also in AD builds. ---*/ + + const auto nPointsPiece = usedPoints.size(); + pieceData.resize(nPointsPiece * nFields); + + vector answerData(answeredPoints.size()); + + for (size_t iField = 0; iField < nFields; iField++) { + for (size_t i = 0; i < answeredPoints.size(); i++) + answerData[i] = dataSorter->GetData(iField, answeredPoints[i] - firstPoint); + + SelectMPIWrapper::W::Alltoallv(answerData.data(), nAnswer.data(), answerOffset.data(), MPI_DOUBLE, + &pieceData[iField * nPointsPiece], nRequest.data(), + requestOffset.data(), MPI_DOUBLE, SU2_MPI::GetComm()); + } + + return usedPoints.size(); +} + +void CParaviewPVTUFileWriter::WritePiece(const string& val_filename, unsigned long nPointsPiece, bool connInt64) { + const auto& fieldNames = dataSorter->GetFieldNames(); + const auto nFields = fieldNames.size(); + const auto nDim = dataSorter->GetnDim(); + + /*--- Local index of a point of the piece. ---*/ + + auto localIndex = [&](unsigned long iPoint) { + return static_cast(lower_bound(usedPoints.begin(), usedPoints.end(), iPoint) - usedPoints.begin()); + }; + + /*--- Connectivity and cell types of the piece. ---*/ + + const auto nElem = dataSorter->GetnElem(); + const auto nConn = dataSorter->GetnConn(); + + vector connectivity, offsets; + vector types; + connectivity.reserve(nConn); + offsets.reserve(nElem); + types.reserve(nElem); + + for (auto type : ElemTypes) { + const auto nNodes = nPointsOfElementType(type); + for (unsigned long iElem = 0; iElem < dataSorter->GetnElem(type); iElem++) { + for (unsigned short iNode = 0; iNode < nNodes; iNode++) + connectivity.push_back(localIndex(dataSorter->GetElemConnectivity(type, iElem, iNode) - 1)); + offsets.push_back(connectivity.size()); + types.push_back(static_cast(type)); + } + } + + /*--- Sizes of the arrays in the file, needed for the offsets of the appended data. ---*/ + + const size_t realSize = doublePrecision ? sizeof(double) : sizeof(float); + const string realType = doublePrecision ? "Float64" : "Float32"; + + const size_t connSize = connInt64 ? sizeof(int64_t) : sizeof(int32_t); + const string connType = connInt64 ? "Int64" : "Int32"; + + /*--- Header of the piece, the offset of each array in the appended data is written with it. ---*/ + + ostringstream header; + size_t offset = 0; + + header << "\n"; + header << "\n"; + header << "\n"; + + header << "\n\n\n"; + offset += DataArraySize(nPointsPiece * NCOORDS, realSize); + + header << "\n"; + header << "\n"; + offset += DataArraySize(connectivity.size(), connSize); + header << "\n"; + offset += DataArraySize(offsets.size(), connSize); + header << "\n"; + offset += DataArraySize(types.size(), sizeof(uint8_t)); + header << "\n"; + + /*--- The fields start after the coordinates, a field whose name ends in _x is written as a vector + together with the fields that follow it. ---*/ + + const size_t varStart = (nDim == 3) ? 3 : 2; + + header << "\n"; + for (size_t iField = varStart; iField < nFields; iField++) { + if (fieldNames[iField].find("_y") != string::npos || fieldNames[iField].find("_z") != string::npos) continue; + + const bool isVector = fieldNames[iField].find("_x") != string::npos; + string name = fieldNames[iField]; + if (isVector) { + name = FieldBaseName(name); + } else { + name.erase(remove(name.begin(), name.end(), '"'), name.end()); + } + + header << "\n"; + offset += DataArraySize(nPointsPiece * (isVector ? NCOORDS : 1), realSize); + } + header << "\n\n\n\n_"; + + /*--- Each rank writes its own file, no rank waits for another one. The arrays are written one at a time. ---*/ + + FILE* file = OpenFile(val_filename, "wb"); + if (!file) SU2_MPI::Error(string("Unable to open file ") + val_filename, CURRENT_FUNCTION); + + const auto headerString = header.str(); + fwrite(headerString.data(), sizeof(char), headerString.size(), file); + fileSize = headerString.size(); + + vector buffer; + + auto writeField = [&](size_t iField, unsigned short nComponents) { + buffer.assign(nPointsPiece * nComponents, 0.0); + for (unsigned short iComp = 0; iComp < nComponents; iComp++) { + if ((nDim == 2) && (nComponents == NCOORDS) && (iComp == 2)) continue; + const auto* field = &pieceData[(iField + iComp) * nPointsPiece]; + for (unsigned long iPoint = 0; iPoint < nPointsPiece; iPoint++) + buffer[iPoint * nComponents + iComp] = field[iPoint]; + } + fileSize += doublePrecision ? WriteDataArray(file, buffer) : WriteDataArray(file, buffer); + }; + + auto writeConn = [&](const vector& values) { + fileSize += connInt64 ? WriteDataArray(file, values) : WriteDataArray(file, values); + }; + + writeField(0, NCOORDS); + + writeConn(connectivity); + writeConn(offsets); + fileSize += WriteDataArray(file, types); + + for (size_t iField = varStart; iField < nFields; iField++) { + if (fieldNames[iField].find("_y") != string::npos || fieldNames[iField].find("_z") != string::npos) continue; + + const bool isVector = fieldNames[iField].find("_x") != string::npos; + writeField(iField, isVector ? NCOORDS : 1); + } + + const string tail = "\n\n\n"; + fwrite(tail.data(), sizeof(char), tail.size(), file); + fclose(file); + + fileSize += tail.size(); +} + +void CParaviewPVTUFileWriter::WriteIndexFile(const string& val_filename, const vector& pieceNames, + bool connInt64) { + if (rank != MASTER_NODE) return; + + const auto& fieldNames = dataSorter->GetFieldNames(); + const auto nDim = dataSorter->GetnDim(); + + const string realType = doublePrecision ? "Float64" : "Float32"; + const string connType = connInt64 ? "Int64" : "Int32"; + + ostringstream file; + + file << "\n"; + file << "\n"; + file << "\n\n\n"; + file << "\n"; + file << "\n"; + file << "\n"; + file << "\n"; + file << "\n"; + + const size_t varStart = (nDim == 3) ? 3 : 2; + + file << "\n"; + for (size_t iField = varStart; iField < fieldNames.size(); iField++) { + if (fieldNames[iField].find("_y") != string::npos || fieldNames[iField].find("_z") != string::npos) continue; + + const bool isVector = fieldNames[iField].find("_x") != string::npos; + string name = fieldNames[iField]; + if (isVector) { + name = FieldBaseName(name); + } else { + name.erase(remove(name.begin(), name.end(), '"'), name.end()); + } + + file << "\n"; + } + file << "\n"; + + for (const auto& piece : pieceNames) file << "\n"; + + file << "\n\n"; + + const auto fileString = file.str(); + + FILE* indexFile = OpenFile(val_filename, "w"); + if (!indexFile) SU2_MPI::Error(string("Unable to open file ") + val_filename, CURRENT_FUNCTION); + fwrite(fileString.data(), sizeof(char), fileString.size(), indexFile); + fclose(indexFile); + + fileSize += fileString.size(); +} diff --git a/SU2_CFD/src/output/filewriter/CSurfaceFVMDataSorter.cpp b/SU2_CFD/src/output/filewriter/CSurfaceFVMDataSorter.cpp index 55fb2296ad7..9689bffb2db 100644 --- a/SU2_CFD/src/output/filewriter/CSurfaceFVMDataSorter.cpp +++ b/SU2_CFD/src/output/filewriter/CSurfaceFVMDataSorter.cpp @@ -1167,6 +1167,8 @@ void CSurfaceFVMDataSorter::SortSurfaceConnectivity(CConfig *config, CGeometry * auto it = std::find(markerList.begin(), markerList.end(), markerTag); if (it != markerList.end()) { + /*--- Element indices restart at zero for each marker. ---*/ + std::fill(nElem_Flag, nElem_Flag + size, -1); for (int ii = 0; ii < (int)geometry->GetnElem_Bound(iMarker); ii++) { @@ -1257,6 +1259,8 @@ void CSurfaceFVMDataSorter::SortSurfaceConnectivity(CConfig *config, CGeometry * auto it = std::find(markerList.begin(), markerList.end(), markerTag); if (it != markerList.end()) { + /*--- Element indices restart at zero for each marker. ---*/ + std::fill(nElem_Flag, nElem_Flag + size, -1); for (int ii = 0; ii < (int)geometry->GetnElem_Bound(iMarker); ii++) { diff --git a/TestCases/output_writers/pvtu_output.cfg b/TestCases/output_writers/pvtu_output.cfg new file mode 100644 index 00000000000..24d977b77a8 --- /dev/null +++ b/TestCases/output_writers/pvtu_output.cfg @@ -0,0 +1,38 @@ +% Test of the Paraview parallel writer (PARAVIEW_PARALLEL and SURFACE_PARAVIEW_PARALLEL). +% The solution of vandv/rans/flatplate is written without iterating (ITER= 0) on 2 ranks. +% Each rank writes its own .vtu piece in the folders flow/ and surface_flow/, and rank 0 writes +% the .pvtu files that list the pieces and the fields. The .pvtu files are compared with +% reference files. The pieces are binary, so they are not compared here. +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +SOLVER= RANS +KIND_TURB_MODEL= SST +SST_OPTIONS= V1994m +RESTART_SOL= YES +% ----------- COMPRESSIBLE AND INCOMPRESSIBLE FREE-STREAM DEFINITION ----------% +MACH_NUMBER= 0.2 +AOA= 0.0 +FREESTREAM_TEMPERATURE= 300.0 +REYNOLDS_NUMBER= 5000000.0 +REYNOLDS_LENGTH= 1.0 +FREESTREAM_TURBULENCEINTENSITY= 0.00038729 +FREESTREAM_TURB2LAMVISCRATIO= 0.009 +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +MARKER_HEATFLUX= ( wall, 0.0 ) +MARKER_FAR= ( farfield ) +MARKER_INLET= ( inlet, 302.4, 117691.7874, 1.0, 0.0, 0.0 ) +MARKER_OUTLET= ( outlet, 114455.0 ) +MARKER_SYM= ( symmetry ) +MARKER_PLOTTING= ( wall ) +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +CONV_NUM_METHOD_FLOW= ROE +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +ITER= 0 +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +MESH_FILENAME= ../vandv/rans/flatplate/mesh_flatplate_turb_035x025.su2 +MESH_FORMAT= SU2 +SOLUTION_FILENAME= ../vandv/rans/flatplate/solution_flow +VOLUME_FILENAME= flow +SURFACE_FILENAME= surface_flow +OUTPUT_FILES= (PARAVIEW_PARALLEL, SURFACE_PARAVIEW_PARALLEL) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 89bdeb0fe11..3b48f86d124 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -2113,6 +2113,19 @@ def main(): pass_list.append(ascii_mesh_writer.run_filediff()) test_list.append(ascii_mesh_writer) + # Paraview parallel writer: each rank writes its .vtu piece, compare the .pvtu files that list them. + for tag, test_file in [("pvtu_writer", "flow.pvtu"), ("pvtu_surface_writer", "surface_flow.pvtu")]: + pvtu_writer = TestCase(tag) + pvtu_writer.cfg_dir = "output_writers" + pvtu_writer.cfg_file = "pvtu_output.cfg" + pvtu_writer.test_iter = -1 + pvtu_writer.command = TestCase.Command("mpirun -n 2", "SU2_CFD") + pvtu_writer.timeout = 1600 + pvtu_writer.reference_file = test_file + ".ref" + pvtu_writer.test_file = test_file + pass_list.append(pvtu_writer.run_filediff()) + test_list.append(pvtu_writer) + ###################################### ### RUN SU2_DEF TESTS ### ###################################### diff --git a/UnitTests/SU2_CFD/output/check_surface_marker_connectivity.py b/UnitTests/SU2_CFD/output/check_surface_marker_connectivity.py new file mode 100644 index 00000000000..e52ad53035c --- /dev/null +++ b/UnitTests/SU2_CFD/output/check_surface_marker_connectivity.py @@ -0,0 +1,207 @@ +#!/usr/bin/env python3 +"""Check that surface sorting preserves cells across distinct marker-local indices. + +Run against SU2_CFD, optionally with --launcher mpirun -np 4. Each probe runs in +an isolated temporary directory. Select --case box9 for larger MPI rank counts. +""" +import argparse +from pathlib import Path +import subprocess +import tempfile + + +CUBE_MARKERS = ("x_minus", "x_plus", "y_minus", "y_plus", "z_minus", "z_plus") + + +def canonical(nodes): + return tuple(sorted(tuple(round(v, 12) for v in p) for p in nodes)) + + +def box(n, markers): + values = [i / (n - 1) for i in range(n)] + faces = {} + for axis in range(3): + other = [i for i in range(3) if i != axis] + for side in range(2): + name = CUBE_MARKERS[2 * axis + side] + faces[name] = [] + for i in range(n - 1): + for j in range(n - 1): + nodes = [] + for a, b in ((i, j), (i + 1, j), (i + 1, j + 1), (i, j + 1)): + p = [0.0] * 3 + p[axis] = side + p[other[0]] = values[a] + p[other[1]] = values[b] + nodes.append(p) + faces[name].append(canonical(nodes)) + return set(face for name in markers for face in faces[name]) + + +def mesh_fixture(kind): + if kind == "square": + points = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + dim = 2 + volume = "9 0 1 2 3 0" + boundaries = [(3, 0), (1, 2), (0, 1), (2, 3)] + vtk = 3 + else: + points = [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)] + dim = 3 + volume = "10 0 1 2 3 0" + boundaries = [(1, 2, 3), (0, 3, 2), (0, 1, 3), (0, 2, 1)] + vtk = 5 + markers = tuple("face" + str(i) for i in range(4)) + lines = ["NDIME= " + str(dim), "NELEM= 1", volume, "NPOIN= 4"] + lines += [" ".join(map(str, p[:dim])) + " " + str(i) for i, p in enumerate(points)] + lines += ["NMARK= 4"] + for name, nodes in zip(markers, boundaries): + lines += [ + "MARKER_TAG= " + name, + "MARKER_ELEMS= 1", + str(vtk) + " " + " ".join(map(str, nodes)), + ] + return ( + "\n".join(lines) + "\n", + markers, + {canonical([points[i] for i in face]) for face in boundaries}, + ) + + +def read_surface(path): + tokens = path.read_text().split() + start = tokens.index("POINTS") + n = int(tokens[start + 1]) + offset = start + 3 + points = [ + tuple(map(float, tokens[offset + 3 * i : offset + 3 * i + 3])) for i in range(n) + ] + start = tokens.index("CELLS") + count = int(tokens[start + 1]) + offset = start + 3 + faces = [] + for _ in range(count): + width = int(tokens[offset]) + offset += 1 + ids = list(map(int, tokens[offset : offset + width])) + offset += width + faces.append(canonical([points[i] for i in ids])) + return faces + + +def run(binary, launcher, name): + markers = CUBE_MARKERS + with tempfile.TemporaryDirectory(prefix="su2-surface-marker-") as temporary: + work = Path(temporary) + if name in ("square", "tetrahedron"): + mesh, markers, expected = mesh_fixture(name) + (work / "mesh.su2").write_text(mesh) + mesh_config = "MESH_FORMAT= SU2\nMESH_FILENAME= mesh.su2\n" + else: + n = { + "single_hex": 2, + "box3": 3, + "box9": 9, + "subset": 2, + "single_marker": 2, + }[name] + if name == "subset": + markers = ("x_minus", "z_plus") + elif name == "single_marker": + markers = ("x_minus",) + expected = box(n, markers) + mesh_config = ( + "MESH_FORMAT= BOX\nMESH_FILENAME= unused.su2\nMESH_BOX_SIZE= (%s,%s,%s)\nMESH_BOX_LENGTH= (1,1,1)\n" + % (n, n, n) + ) + far = ",".join( + CUBE_MARKERS if name not in ("square", "tetrahedron") else markers + ) + cfg = ( + """SOLVER= EULER +MACH_NUMBER= 0.5 +ITER= 0 +MGLEVEL= 0 +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= NO +TIME_DISCRE_FLOW= RUNGE-KUTTA_EXPLICIT +OUTPUT_FILES= (SURFACE_PARAVIEW_ASCII) +SURFACE_FILENAME= surface +HISTORY_OUTPUT= (ITER) +SCREEN_OUTPUT= (INNER_ITER) +MARKER_FAR= (""" + + far + + """) +MARKER_PLOTTING= (""" + + ",".join(markers) + + ")\n" + + mesh_config + ) + (work / "case.cfg").write_text(cfg) + result = subprocess.run( + launcher + [str(binary), "case.cfg"], + cwd=work, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + text=True, + timeout=120, + ) + if result.returncode: + raise RuntimeError(name + " solver failed:\n" + result.stdout[-4000:]) + faces = read_surface(work / "surface.vtk") + if ( + len(faces) != len(expected) + or len(set(faces)) != len(faces) + or set(faces) != expected + ): + raise AssertionError( + "%s: wrote %s cells, expected %s; missing %s, unexpected %s" + % ( + name, + len(faces), + len(expected), + len(expected - set(faces)), + len(set(faces) - expected), + ) + ) + print( + "PASS: " + + name + + " preserves all " + + str(len(expected)) + + " distinct surface cells", + flush=True, + ) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("binary", type=Path) + parser.add_argument( + "--case", + action="append", + choices=( + "single_hex", + "box3", + "box9", + "subset", + "single_marker", + "square", + "tetrahedron", + ), + ) + parser.add_argument("--launcher", nargs=argparse.REMAINDER, default=[]) + args = parser.parse_args() + for name in args.case or ( + "single_hex", + "box3", + "subset", + "single_marker", + "square", + "tetrahedron", + ): + run(args.binary.resolve(), args.launcher, name) + + +if __name__ == "__main__": + main() diff --git a/config_template.cfg b/config_template.cfg index 0aa29398b90..95ff52a533c 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -2713,7 +2713,9 @@ MULTIZONE_ADAPT_FILENAME= YES % Files to output % Possible formats : (TECPLOT_ASCII, TECPLOT, SURFACE_TECPLOT_ASCII, % SURFACE_TECPLOT, CSV, SURFACE_CSV, PARAVIEW_ASCII, PARAVIEW_LEGACY, SURFACE_PARAVIEW_ASCII, -% SURFACE_PARAVIEW_LEGACY, PARAVIEW, SURFACE_PARAVIEW, RESTART_ASCII, RESTART, CGNS, SURFACE_CGNS, STL_ASCII, STL_BINARY) +% SURFACE_PARAVIEW_LEGACY, PARAVIEW, SURFACE_PARAVIEW, PARAVIEW_PARALLEL, SURFACE_PARAVIEW_PARALLEL, +% RESTART_ASCII, RESTART, CGNS, SURFACE_CGNS, STL_ASCII, STL_BINARY) +% PARAVIEW_PARALLEL writes a .pvtu file and a folder with one .vtu piece per MPI rank. % default : (RESTART, PARAVIEW, SURFACE_PARAVIEW) OUTPUT_FILES= (RESTART, PARAVIEW, SURFACE_PARAVIEW) %