From 15030bcef6935ecfe2d9d2e36bd912069fe4e634 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Sat, 28 Oct 2023 16:26:18 +0200 Subject: [PATCH 01/34] Add insphere calculator, add sparse alltoallv --- CMakeLists.txt | 2 + src/aux/InsphereCalculator.cpp | 148 +++++++++++++++++++++++++++++++ src/aux/InsphereCalculator.h | 12 +++ src/aux/MPIConvenience.cpp | 51 +++++++++++ src/aux/MPIConvenience.h | 14 +++ src/input/MeshData.h | 60 +++++++++++++ src/input/ParallelVertexFilter.h | 12 +-- 7 files changed, 294 insertions(+), 5 deletions(-) create mode 100644 src/aux/InsphereCalculator.cpp create mode 100644 src/aux/InsphereCalculator.h create mode 100644 src/aux/MPIConvenience.cpp create mode 100644 src/aux/MPIConvenience.h create mode 100644 src/input/MeshData.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 09dcb2b..3baa622 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,8 @@ add_executable(pumgen ${CMAKE_CURRENT_SOURCE_DIR}/src/third_party/GMSHLexer.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/third_party/GMSHParser.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/third_party/GMSH2Parser.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/aux/InsphereCalculator.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/aux/MPIConvenience.cpp ) target_include_directories(pumgen PUBLIC src diff --git a/src/aux/InsphereCalculator.cpp b/src/aux/InsphereCalculator.cpp new file mode 100644 index 0000000..6498a8f --- /dev/null +++ b/src/aux/InsphereCalculator.cpp @@ -0,0 +1,148 @@ +#include "InsphereCalculator.h" +#include "third_party/MPITraits.h" +#include +#include +#include +#include +#include + +#include "MPIConvenience.h" + +std::vector calculateInsphere(const std::vector& connectivity, + const std::vector& geometry, MPI_Comm comm) { + int commsize; + int commrank; + + MPI_Comm_size(comm, &commsize); + MPI_Comm_rank(comm, &commrank); + + MPI_Datatype vertexType; + + MPI_Type_contiguous(3, MPI_DOUBLE, &vertexType); + MPI_Type_commit(&vertexType); + + std::vector outrequests(commsize); + std::vector inrequests(commsize); + + std::size_t localVertices = geometry.size() / 3; + + std::vector vertexDist(commsize + 1); + + MPI_Allgather(&localVertices, 1, tndm::mpi_type_t(), vertexDist.data() + 1, 1, + tndm::mpi_type_t(), comm); + + for (std::size_t i = 0; i < commsize; ++i) { + vertexDist[i + 1] += vertexDist[i]; + } + + for (const auto& vertex : connectivity) { + auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); + auto position = std::distance(vertexDist.begin(), itPosition) - 1; + ++outrequests[position]; + } + + MPI_Alltoall(outrequests.data(), 1, MPI_INT, inrequests.data(), 1, MPI_INT, comm); + + std::vector outdisp(commsize + 1); + std::vector indisp(commsize + 1); + for (std::size_t i = 1; i < commsize + 1; ++i) { + outdisp[i] = outdisp[i - 1] + outrequests[i - 1]; + indisp[i] = indisp[i - 1] + inrequests[i - 1]; + } + + std::vector inidx(indisp[commsize]); + + { + std::vector outidx(connectivity.size()); + + std::vector counter(commsize); + + for (std::size_t i = 0; i < connectivity.size(); ++i) { + auto vertex = connectivity[i]; + auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); + auto position = std::distance(vertexDist.begin(), itPosition) - 1; + outidx[counter[position] + outdisp[position]] = vertex - vertexDist[position]; + ++counter[position]; + } + + // transfer indices + + sparseAlltoallv(outidx.data(), outrequests.data(), outdisp.data(), + tndm::mpi_type_t(), inidx.data(), inrequests.data(), indisp.data(), + tndm::mpi_type_t(), comm); + } + + // this is a bit inefficient, since we don't look for duplicates... TODO: maybe improve + std::vector outvertices(3 * connectivity.size()); + + { + std::vector invertices(3 * indisp[commsize]); + + for (std::size_t i = 0; i < inidx.size(); ++i) { + for (int j = 0; j < 3; ++j) { + invertices[3 * i + j] = geometry[3 * inidx[i] + j]; + } + } + + // transfer vertices + + sparseAlltoallv(invertices.data(), inrequests.data(), indisp.data(), vertexType, + outvertices.data(), outrequests.data(), outdisp.data(), vertexType, comm); + } + + std::vector counter(commsize); + std::vector inspheres(connectivity.size() / 4); + + for (std::size_t i = 0; i < connectivity.size() / 4; ++i) { + std::array vertexIds{0, 0, 0, 0}; + for (int j = 0; j < 4; ++j) { + auto vertex = connectivity[i * 4 + j]; + auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); + auto position = std::distance(vertexDist.begin(), itPosition) - 1; + vertexIds[j] = counter[position] + outdisp[position]; + ++counter[position]; + } + double a11 = outvertices[3 * vertexIds[1] + 0] - outvertices[3 * vertexIds[0] + 0]; + double a12 = outvertices[3 * vertexIds[1] + 1] - outvertices[3 * vertexIds[0] + 1]; + double a13 = outvertices[3 * vertexIds[1] + 2] - outvertices[3 * vertexIds[0] + 2]; + double a21 = outvertices[3 * vertexIds[2] + 0] - outvertices[3 * vertexIds[0] + 0]; + double a22 = outvertices[3 * vertexIds[2] + 1] - outvertices[3 * vertexIds[0] + 1]; + double a23 = outvertices[3 * vertexIds[2] + 2] - outvertices[3 * vertexIds[0] + 2]; + double a31 = outvertices[3 * vertexIds[3] + 0] - outvertices[3 * vertexIds[0] + 0]; + double a32 = outvertices[3 * vertexIds[3] + 1] - outvertices[3 * vertexIds[0] + 1]; + double a33 = outvertices[3 * vertexIds[3] + 2] - outvertices[3 * vertexIds[0] + 2]; + double b11 = outvertices[3 * vertexIds[2] + 0] - outvertices[3 * vertexIds[1] + 0]; + double b12 = outvertices[3 * vertexIds[2] + 1] - outvertices[3 * vertexIds[1] + 1]; + double b13 = outvertices[3 * vertexIds[2] + 2] - outvertices[3 * vertexIds[1] + 2]; + double b21 = outvertices[3 * vertexIds[3] + 0] - outvertices[3 * vertexIds[1] + 0]; + double b22 = outvertices[3 * vertexIds[3] + 1] - outvertices[3 * vertexIds[1] + 1]; + double b23 = outvertices[3 * vertexIds[3] + 2] - outvertices[3 * vertexIds[1] + 2]; + double det = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a13 * a22 * a31 - + a12 * a21 * a33 - a11 * a23 * a32; + double gram = std::abs(det); + + double a1a2x1 = a12 * a23 - a13 * a22; + double a1a2x2 = a13 * a21 - a11 * a23; + double a1a2x3 = a11 * a22 - a12 * a21; + double a1a3x1 = a12 * a33 - a13 * a32; + double a1a3x2 = a13 * a31 - a11 * a33; + double a1a3x3 = a11 * a32 - a12 * a31; + double a2a3x1 = a22 * a33 - a23 * a32; + double a2a3x2 = a23 * a31 - a21 * a33; + double a2a3x3 = a21 * a32 - a22 * a31; + double b1b2x1 = b12 * b23 - b13 * b22; + double b1b2x2 = b13 * b21 - b11 * b23; + double b1b2x3 = b11 * b22 - b12 * b21; + + double faces = std::sqrt(a1a2x1 * a1a2x1 + a1a2x2 * a1a2x2 + a1a2x3 * a1a2x3) + + std::sqrt(a1a3x1 * a1a3x1 + a1a3x2 * a1a3x2 + a1a3x3 * a1a3x3) + + std::sqrt(a2a3x1 * a2a3x1 + a2a3x2 * a2a3x2 + a2a3x3 * a2a3x3) + + std::sqrt(b1b2x1 * b1b2x1 + b1b2x2 * b1b2x2 + b1b2x3 * b1b2x3); + + inspheres[i] = gram / faces; + } + + MPI_Type_free(&vertexType); + + return inspheres; +} diff --git a/src/aux/InsphereCalculator.h b/src/aux/InsphereCalculator.h new file mode 100644 index 0000000..a019fb6 --- /dev/null +++ b/src/aux/InsphereCalculator.h @@ -0,0 +1,12 @@ +#ifndef PUMGEN_AUX_INSPHERE_CALCULATOR_H_ +#define PUMGEN_AUX_INSPHERE_CALCULATOR_H_ + +#include +#include + +// assumes a contiguous distribution of all vertices over all processes + +std::vector calculateInsphere(const std::vector& connectivity, + const std::vector& geometry, MPI_Comm comm); + +#endif // PUMGEN_AUX_INSPHERE_CALCULATOR_H_ diff --git a/src/aux/MPIConvenience.cpp b/src/aux/MPIConvenience.cpp new file mode 100644 index 0000000..2e8eef6 --- /dev/null +++ b/src/aux/MPIConvenience.cpp @@ -0,0 +1,51 @@ +#include "MPIConvenience.h" +#include +#include + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { + int commsize; + int commrank; + int sendtypesizePre; + int recvtypesizePre; + MPI_Comm_size(comm, &commsize); + MPI_Comm_rank(comm, &commrank); + MPI_Type_size(sendtype, &sendtypesizePre); + MPI_Type_size(sendtype, &recvtypesizePre); + + constexpr int Tag = 1000; + + std::size_t sendtypesize = sendtypesizePre; + std::size_t recvtypesize = recvtypesizePre; + + std::vector requests(commsize * 2, MPI_REQUEST_NULL); + + // (no special handling of self-to-self comm at the moment) + + for (int i = 0; i < commsize; ++i) { + if (sendsize[i] > 0) { + MPI_Isend(reinterpret_cast(sendbuf) + senddisp[i] * sendtypesize, sendsize[i], + sendtype, i, Tag, comm, &requests[i]); + } + } + for (int i = 0; i < commsize; ++i) { + if (recvsize[i] > 0) { + MPI_Irecv(reinterpret_cast(recvbuf) + recvdisp[i] * recvtypesize, recvsize[i], + recvtype, i, Tag, comm, &requests[i + commsize]); + } + } + + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); +} + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, const int* recvdisp, + MPI_Datatype recvtype, MPI_Comm comm) { + int commsize; + MPI_Comm_size(comm, &commsize); + std::vector senddispSize(senddisp, senddisp + commsize); + std::vector recvdispSize(recvdisp, recvdisp + commsize); + sparseAlltoallv(sendbuf, sendsize, senddispSize.data(), sendtype, recvbuf, recvsize, + recvdispSize.data(), recvtype, comm); +} diff --git a/src/aux/MPIConvenience.h b/src/aux/MPIConvenience.h new file mode 100644 index 0000000..128439f --- /dev/null +++ b/src/aux/MPIConvenience.h @@ -0,0 +1,14 @@ +#ifndef PUMGEN_AUX_MPI_CONVENIENCE_H_ +#define PUMGEN_AUX_MPI_CONVENIENCE_H_ + +#include + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm); + +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, const int* recvdisp, + MPI_Datatype recvtype, MPI_Comm comm); + +#endif // PUMGEN_AUX_MPI_CONVENIENCE_H_ diff --git a/src/input/MeshData.h b/src/input/MeshData.h new file mode 100644 index 0000000..526e5ca --- /dev/null +++ b/src/input/MeshData.h @@ -0,0 +1,60 @@ + +#ifndef MESH_DATA_H +#define MESH_DATA_H + +#include +#include + +/** + * Interface for mesh input + */ +class MeshData { + public: + virtual std::size_t cellCount() = 0; + virtual std::size_t vertexCount() = 0; + + virtual const std::vector& connectivity() = 0; + virtual const std::vector& geometry() = 0; + virtual const std::vector& group() = 0; + virtual const std::vector& boundary() = 0; +}; + +class FullStorageMeshData : public MeshData { + public: + virtual std::size_t cellCount() { return cellCountValue; } + virtual std::size_t vertexCount() { return vertexCountValue; } + + virtual const std::vector& connectivity() { return connectivityData; } + virtual const std::vector& geometry() { return geometryData; } + virtual const std::vector& group() { return groupData; } + virtual const std::vector& boundary() { return boundaryData; } + + protected: + std::size_t cellCountValue = 0; + std::size_t vertexCountValue = 0; + + std::vector connectivityData; + std::vector geometryData; + std::vector groupData; + std::vector boundaryData; + + int bndShift = 8; + + void setBoundary(std::size_t cell, int face, int value) { + boundaryData[cell] = static_cast(value) << (face * bndShift); + } + + void setup(std::size_t cellCount, std::size_t vertexCount) { + cellCountValue = cellCount; + vertexCountValue = vertexCount; + + connectivityData.resize(cellCount * 4); + geometryData.resize(vertexCount * 3); + groupData.resize(cellCount); + boundaryData.resize(cellCount); + + // TODO: reconsider the times 4 or 3 + } +}; + +#endif // MESH_INPUT_H diff --git a/src/input/ParallelVertexFilter.h b/src/input/ParallelVertexFilter.h index 96b1136..5a34909 100644 --- a/src/input/ParallelVertexFilter.h +++ b/src/input/ParallelVertexFilter.h @@ -30,6 +30,8 @@ #include "third_party/MPITraits.h" +#include "aux/MPIConvenience.h" + /** * Filters duplicate vertices in parallel */ @@ -200,8 +202,8 @@ class ParallelVertexFilter { sDispls[i] = sDispls[i - 1] + bucketSize[i - 1]; rDispls[i] = rDispls[i - 1] + recvSize[i - 1]; } - MPI_Alltoallv(sendVertices.data(), bucketSize.data(), sDispls.data(), vertexType, - sortVertices.data(), recvSize.data(), rDispls.data(), vertexType, m_comm); + sparseAlltoallv(sendVertices.data(), bucketSize.data(), sDispls.data(), vertexType, + sortVertices.data(), recvSize.data(), rDispls.data(), vertexType, m_comm); // Chop the last 4 bits to avoid numerical errors roundVertices.resize(numSortVertices * 3); @@ -253,9 +255,9 @@ class ParallelVertexFilter { // Send result back std::vector globalIds(numVertices); - MPI_Alltoallv(gids.data(), recvSize.data(), rDispls.data(), tndm::mpi_type_t(), - globalIds.data(), bucketSize.data(), sDispls.data(), - tndm::mpi_type_t(), m_comm); + sparseAlltoallv(gids.data(), recvSize.data(), rDispls.data(), tndm::mpi_type_t(), + globalIds.data(), bucketSize.data(), sDispls.data(), + tndm::mpi_type_t(), m_comm); // Assign the global ids to the correct vertices m_globalIds.resize(numVertices); From 0ade0857303f3b8780e7816be34e893fb33c8ce9 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Sat, 28 Oct 2023 16:28:46 +0200 Subject: [PATCH 02/34] Replace APF for SimModSuite Disable other mesh outputs temporarily --- src/input/SimModSuite.h | 207 ++++++++--- src/input/SimModSuiteApf.h | 688 +++++++++++++++++++++++++++++++++++++ src/pumgen.cpp | 184 ++-------- 3 files changed, 880 insertions(+), 199 deletions(-) create mode 100644 src/input/SimModSuiteApf.h diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index b28350a..1e834c6 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -6,8 +6,9 @@ * notice in the file 'COPYING' at the root directory of this package * and the copyright notice at https://github.com/SeisSol/PUMGen * - * @copyright 2017 Technical University of Munich + * @copyright 2017,2023 Technical University of Munich * @author Sebastian Rettenberger + * @author David Schneller */ #ifndef SIM_MOD_SUITE_H @@ -20,12 +21,6 @@ #include #include -#include -#include -#include -#include -#include - #include #include #ifdef BEFORE_SIM_18 @@ -48,12 +43,14 @@ #include "utils/path.h" #include "utils/progress.h" -#include "MeshInput.h" +#include "MeshData.h" #include "AnalysisAttributes.h" #include "EasiMeshSize.h" #include "MeshAttributes.h" +#include "ParallelVertexFilter.h" + #include #include #include @@ -72,7 +69,7 @@ pAManager SModel_attManager(pModel model); * of this class * @todo Maybe add MS_setMaxEntities to limit the number of elements */ -class SimModSuite : public MeshInput { +class SimModSuite : public FullStorageMeshData { private: EasiMeshSize easiMeshSize; pGModel m_model; @@ -174,51 +171,171 @@ class SimModSuite : public MeshInput { analyse_mesh(); } - // Convert to APF mesh - apf::Mesh* tmpMesh = apf::createMesh(m_simMesh); - gmi_register_sim(); - gmi_model* model = gmi_import_sim(m_model); + logInfo(PMU_rank()) << "Iterating over mesh to get data..."; + int parts = PM_numParts(m_simMesh); + std::size_t vertexCount = 0; + std::size_t vertexIDCount = 0; + std::size_t cellCount = 0; + + std::vector vertexStart(parts + 1); + std::vector vertexIDStart(parts + 1); + std::vector cellStart(parts + 1); + for (int i = 0; i < parts; i++) { + logInfo(PMU_rank()) << "Counting part" << i << "/" << parts; + vertexStart[i] = vertexCount; + vertexIDStart[i] = vertexIDCount; + cellStart[i] = cellCount; + { + size_t totalCount = 0; + size_t ownedCount = 0; + VIter reg_it = M_vertexIter(PM_mesh(m_simMesh, i)); + while (pVertex vertex = VIter_next(reg_it)) { + ++totalCount; + if (EN_isOwnerProc((pEntity)vertex)) { + ++ownedCount; + } + } + VIter_delete(reg_it); + + vertexCount += ownedCount; + vertexIDCount += totalCount; + } + { + size_t totalCount = 0; + size_t ownedCount = 0; + RIter reg_it = M_regionIter(PM_mesh(m_simMesh, i)); + while (pRegion reg = RIter_next(reg_it)) { + ++totalCount; + if (EN_isOwnerProc((pEntity)reg)) { + ++ownedCount; + } + } + RIter_delete(reg_it); + + cellCount += ownedCount; + assert(totalCount == ownedCount); + } + } + vertexStart.back() = vertexCount; + cellStart.back() = cellCount; + vertexIDStart.back() = vertexIDCount; - logInfo(PMU_rank()) << "Converting mesh to APF"; - m_mesh = apf::createMdsMesh(model, tmpMesh); - apf::destroyMesh(tmpMesh); + logInfo(PMU_rank()) << "Local cells:" << cellCount; + logInfo(PMU_rank()) << "Local vertices:" << vertexCount; + logInfo(PMU_rank()) << "Local vertices (with duplicates):" << vertexIDCount; - // Set the boundary conditions from the geometric model AttCase_associate(analysisCase, 0L); - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - apf::MeshIterator* it = m_mesh->begin(2); - while (apf::MeshEntity* face = m_mesh->iterate(it)) { - apf::ModelEntity* modelFace = m_mesh->toModel(face); - if (m_mesh->getModelType(modelFace) != 2) - continue; - - pGEntity simFace = reinterpret_cast(modelFace); - - pAttribute attr = GEN_attrib(simFace, "boundaryCondition"); - if (attr) { - char* image = Attribute_imageClass(attr); - int boundary = parseBoundary(image); - Sim_deleteString(image); - - m_mesh->setIntTag(face, boundaryTag, &boundary); + + ParallelVertexFilter vertexFilter; + + { + std::vector vertexData(vertexIDCount * 3); + for (int i = 0; i < parts; i++) { + logInfo(PMU_rank()) << "Processing part" << i << "/" << parts; + logInfo(PMU_rank()) << "Vertices:" << vertexStart[i] << "to" << vertexStart[i + 1]; + { + std::size_t indexID = vertexIDStart[i]; + VIter reg_it = M_vertexIter(PM_mesh(m_simMesh, i)); + while (pVertex vertex = VIter_next(reg_it)) { + // TODO: EN_isOwnerProc() + EN_setID((pEntity)vertex, indexID); + double xyz[3]; + V_coord(vertex, xyz); + vertexData[indexID * 3 + 0] = xyz[0]; + vertexData[indexID * 3 + 1] = xyz[1]; + vertexData[indexID * 3 + 2] = xyz[2]; + ++indexID; + } + VIter_delete(reg_it); + } } + vertexFilter.filter(vertexIDCount, vertexData); } - m_mesh->end(it); + vertexCount = vertexFilter.numLocalVertices(); + + logInfo(PMU_rank()) << "Local vertices (really, now):" << vertexCount; + + setup(cellCount, vertexCount); + std::copy(vertexFilter.localVertices().begin(), vertexFilter.localVertices().end(), + geometryData.begin()); + for (int i = 0; i < parts; i++) { + logInfo(PMU_rank()) << "Processing part" << i << "/" << parts; + logInfo(PMU_rank()) << "Connectivity:" << cellStart[i] << "to" << cellStart[i + 1]; + { + std::size_t index = cellStart[i]; + RIter reg_it = M_regionIter(PM_mesh(m_simMesh, i)); + pRegion reg; + while ((reg = RIter_next(reg_it))) { + if (EN_isOwnerProc((pEntity)reg)) { + pPList vertices = R_vertices(reg, 1); + void* iter = nullptr; + int j = 0; + while (pVertex vertex = (pVertex)PList_next(vertices, &iter)) { + EN_setID((pEntity)reg, index); + connectivityData[index * 4 + j] = vertexFilter.globalIds()[EN_id((pEntity)vertex)]; + ++j; + } + PList_delete(vertices); + ++index; + } + } + RIter_delete(reg_it); + } - // Set groups - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - it = m_mesh->begin(3); - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::ModelEntity* modelRegion = m_mesh->toModel(element); + logInfo(PMU_rank()) << "Groups"; + { + int groupIdx = 0; + GRIter modelRegionIter = GM_regionIter(m_model); + while (pGRegion mreg = GRIter_next(modelRegionIter)) { + RIter reg_it = M_classifiedRegionIter(PM_mesh(m_simMesh, i), mreg); + while (pRegion reg = RIter_next(reg_it)) { + if (EN_isOwnerProc((pEntity)reg)) { + groupData[EN_id((pEntity)reg)] = groupIdx; + } + } + RIter_delete(reg_it); - pGRegion simRegion = reinterpret_cast(modelRegion); - std::unordered_map::const_iterator i = groupMap.find(simRegion); - if (i == groupMap.end()) - logError() << "Mesh element with unknown region found."; + ++groupIdx; + } + GRIter_delete(modelRegionIter); + } - m_mesh->setIntTag(element, groupTag, &i->second); + logInfo(PMU_rank()) << "Boundaries"; + { + GFIter modelFaceIter = GM_faceIter(m_model); + while (pGFace mface = GFIter_next(modelFaceIter)) { + pAttribute attr = GEN_attrib(mface, "boundaryCondition"); + if (attr != nullptr) { + char* image = Attribute_imageClass(attr); + int boundary = parseBoundary(image); + FIter face_it = M_classifiedFaceIter(PM_mesh(m_simMesh, i), mface, 1); + while (pFace face = FIter_next(face_it)) { + pRegion r1 = F_region(face, 0); + pRegion r2 = F_region(face, 1); + + // suboptimal: we go over all faces here... + if (r1 != nullptr && EN_isOwnerProc((pEntity)r1)) { + for (int k = 0; k < 4; ++k) { + if (R_face(r1, k) == face) { + setBoundary(EN_id((pEntity)r1), k, boundary); + } + } + } + if (r2 != nullptr && EN_isOwnerProc((pEntity)r2)) { + for (int k = 0; k < 4; ++k) { + if (R_face(r2, k) == face) { + setBoundary(EN_id((pEntity)r2), k, boundary); + } + } + } + } + FIter_delete(face_it); + Sim_deleteString(image); + } + } + GFIter_delete(modelFaceIter); + } } - m_mesh->end(it); AttCase_unassociate(analysisCase); diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h new file mode 100644 index 0000000..b28350a --- /dev/null +++ b/src/input/SimModSuiteApf.h @@ -0,0 +1,688 @@ +/** + * @file + * This file is part of PUMGen + * + * For conditions of distribution and use, please see the copyright + * notice in the file 'COPYING' at the root directory of this package + * and the copyright notice at https://github.com/SeisSol/PUMGen + * + * @copyright 2017 Technical University of Munich + * @author Sebastian Rettenberger + */ + +#ifndef SIM_MOD_SUITE_H +#define SIM_MOD_SUITE_H + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#ifdef BEFORE_SIM_18 +#include +#include +#include +#else +#include +#include +#include +#endif +#include +#include +#ifdef PARASOLID +#include +#endif +#include + +#include "utils/logger.h" +#include "utils/path.h" +#include "utils/progress.h" + +#include "MeshInput.h" + +#include "AnalysisAttributes.h" +#include "EasiMeshSize.h" +#include "MeshAttributes.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// forward declare +pAManager SModel_attManager(pModel model); + +/** + * @todo Currently it is not supported to create more than one instance + * of this class + * @todo Maybe add MS_setMaxEntities to limit the number of elements + */ +class SimModSuite : public MeshInput { + private: + EasiMeshSize easiMeshSize; + pGModel m_model; + + pParMesh m_simMesh; + + /** Enable Simmetrix logging file */ + bool m_log; + + public: + SimModSuite(const char* modFile, const char* cadFile = 0L, const char* licenseFile = 0L, + const char* meshCaseName = "mesh", const char* analysisCaseName = "analysis", + int enforceSize = 0, const char* xmlFile = 0L, const bool analyseAR = false, + const char* logFile = 0L) { + // Init SimModSuite + SimModel_start(); + SimPartitionedMesh_start(0L, 0L); + if (logFile) { + m_log = true; + Sim_logOn(logFile); + } else + m_log = false; + Sim_readLicenseFile(licenseFile); + MS_init(); + SimDiscrete_start(0); +#ifdef PARASOLID + SimParasolid_start(1); +#endif + Sim_setMessageHandler(messageHandler); + + // Load CAD + logInfo(PMU_rank()) << "Loading model"; + + std::string smodFile = modFile; + if (cadFile != 0L) { + loadCAD(modFile, cadFile); + } else if (smodFile.substr(smodFile.find_last_of(".") + 1) == "smd") { + loadCAD(modFile, cadFile); + } else { + logError() << "unsupported model file"; + } + + // Create group lookup table + std::unordered_map groupMap; + GRIter regionIt = GM_regionIter(m_model); + while (pGRegion region = GRIter_next(regionIt)) { + groupMap[region] = static_cast(groupMap.size() + 1); + } + + // Extract cases + logInfo(PMU_rank()) << "Extracting cases"; + pACase meshCase, analysisCase; + MeshAttributes MeshAtt; + + if (xmlFile != nullptr) { + + // Read mesh Attributes from xml file + int numFaces = GM_numFaces(m_model); + MeshAtt.init(xmlFile); + AnalysisAttributes AnalysisAtt(xmlFile, numFaces); + setCases(m_model, meshCase, analysisCase, MeshAtt, AnalysisAtt, groupMap); + } else { + extractCases(m_model, meshCase, meshCaseName, analysisCase, analysisCaseName); + } + + m_simMesh = PM_new(0, m_model, PMU_size()); + + pProgress prog = Progress_new(); + Progress_setCallback(prog, progressHandler); + + // create the mesh + logInfo(PMU_rank()) << "Starting the surface mesher"; + pSurfaceMesher surfaceMesher = SurfaceMesher_new(meshCase, m_simMesh); + if (xmlFile != nullptr) { + SurfaceMesher_setSmoothing(surfaceMesher, MeshAtt.surfaceSmoothingLevel); + SurfaceMesher_setSmoothType(surfaceMesher, static_cast(MeshAtt.surfaceSmoothingType)); + SurfaceMesher_setFaceRotationLimit(surfaceMesher, MeshAtt.surfaceFaceRotationLimit); + SurfaceMesher_setSnapForDiscrete(surfaceMesher, MeshAtt.surfaceSnap); + } + SurfaceMesher_execute(surfaceMesher, prog); + SurfaceMesher_delete(surfaceMesher); + + PM_setTotalNumParts(m_simMesh, PMU_size()); + + logInfo(PMU_rank()) << "Starting the volume mesher"; + pVolumeMesher volumeMesher = VolumeMesher_new(meshCase, m_simMesh); + if (xmlFile != nullptr) { + VolumeMesher_setSmoothing(volumeMesher, MeshAtt.volumeSmoothingLevel); + VolumeMesher_setSmoothType(volumeMesher, static_cast(MeshAtt.volumeSmoothingType)); + VolumeMesher_setOptimization(volumeMesher, MeshAtt.VolumeMesherOptimization); + } + VolumeMesher_setEnforceSize(volumeMesher, enforceSize); + VolumeMesher_execute(volumeMesher, prog); + VolumeMesher_delete(volumeMesher); + + Progress_delete(prog); + + if (analyseAR) { + analyse_mesh(); + } + + // Convert to APF mesh + apf::Mesh* tmpMesh = apf::createMesh(m_simMesh); + gmi_register_sim(); + gmi_model* model = gmi_import_sim(m_model); + + logInfo(PMU_rank()) << "Converting mesh to APF"; + m_mesh = apf::createMdsMesh(model, tmpMesh); + apf::destroyMesh(tmpMesh); + + // Set the boundary conditions from the geometric model + AttCase_associate(analysisCase, 0L); + apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); + apf::MeshIterator* it = m_mesh->begin(2); + while (apf::MeshEntity* face = m_mesh->iterate(it)) { + apf::ModelEntity* modelFace = m_mesh->toModel(face); + if (m_mesh->getModelType(modelFace) != 2) + continue; + + pGEntity simFace = reinterpret_cast(modelFace); + + pAttribute attr = GEN_attrib(simFace, "boundaryCondition"); + if (attr) { + char* image = Attribute_imageClass(attr); + int boundary = parseBoundary(image); + Sim_deleteString(image); + + m_mesh->setIntTag(face, boundaryTag, &boundary); + } + } + m_mesh->end(it); + + // Set groups + apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); + it = m_mesh->begin(3); + while (apf::MeshEntity* element = m_mesh->iterate(it)) { + apf::ModelEntity* modelRegion = m_mesh->toModel(element); + + pGRegion simRegion = reinterpret_cast(modelRegion); + std::unordered_map::const_iterator i = groupMap.find(simRegion); + if (i == groupMap.end()) + logError() << "Mesh element with unknown region found."; + + m_mesh->setIntTag(element, groupTag, &i->second); + } + m_mesh->end(it); + + AttCase_unassociate(analysisCase); + + // Delete cases + MS_deleteMeshCase(meshCase); + MS_deleteMeshCase(analysisCase); + } + + virtual ~SimModSuite() { + M_release(m_simMesh); + // We cannot delete the model here because it is still + // connected to the mesh + // GM_release(m_model); + + // Finalize SimModSuite +#ifdef PARASOLID + SimParasolid_stop(1); +#endif + SimDiscrete_stop(0); + MS_exit(); + Sim_unregisterAllKeys(); + if (m_log) + Sim_logOff(); + SimPartitionedMesh_stop(); + SimModel_stop(); + } + + private: + pACase extractCase(pAManager attMngr, const char* name) { + pACase acase = AMAN_findCase(attMngr, name); + if (!acase) + logError() << "Case" << std::string(name) << "not found."; + + AttCase_setModel(acase, m_model); + + return acase; + } + + private: + /** + * @todo Make this private as soon as APF does no longer need it + */ + static unsigned int parseBoundary(const char* boundaryCondition) { + if (strcmp(boundaryCondition, "freeSurface") == 0) + return 1; + if (strcmp(boundaryCondition, "dynamicRupture") == 0) + return 3; + if (strcmp(boundaryCondition, "absorbing") == 0) + return 5; + + // check if boundaryCondition starts with pattern + std::string pattern = "BC"; + std::string sboundaryCondition(boundaryCondition, 11); + if (sboundaryCondition.find(pattern) == 0) { + std::string sNumber = sboundaryCondition.substr(pattern.length()); + return stoi(sNumber); + } + logError() << "Unknown boundary condition" << boundaryCondition; + return -1; + } + + static void messageHandler(int type, const char* msg) { + switch (type) { + case Sim_InfoMsg: + // Show sim info messages as debug messages + logDebug(PMU_rank()) << "SimModeler:" << msg; + break; + case Sim_DebugMsg: + // Ignore sim debug messages + break; + case Sim_WarningMsg: + logWarning(PMU_rank()) << "SimModeler:" << msg; + break; + case Sim_ErrorMsg: + // Use warning because error will abort the program + logWarning() << "SimModeler:" << msg; + break; + } + } + + static void progressHandler(const char* what, int level, int startVal, int endVal, int currentVal, + void* ignore) { + if (PMU_rank() != 0) + return; + if (endVal != 0) { + switch (currentVal) { + case -2: + // task is started, do nothing + logInfo(PMU_rank()) << "Progress:" << what << ", 0" + << "/" << endVal; + break; + case -1: + // end of the task + logInfo(PMU_rank()) << "Progress:" << what << ", done"; + break; + default: + logInfo(PMU_rank()) << "Progress:" << what << "," << currentVal << "/" << endVal; + break; + } + } else { + if (currentVal == -2) + logInfo(PMU_rank()) << "Progress:" << what; + } + logDebug() << what << level << startVal << endVal << currentVal; + } + + private: + void extractCases(pGModel m_model, pACase& meshCase, const char* meshCaseName, + pACase& analysisCase, const char* analysisCaseName) { + logInfo(PMU_rank()) << "Extracting cases"; + pAManager attMngr = SModel_attManager(m_model); + + MeshingOptions meshingOptions; + meshCase = MS_newMeshCase(m_model); + + pACase meshCaseFile = extractCase(attMngr, meshCaseName); + AttCase_associate(meshCaseFile, nullptr); + MS_processSimModelerMeshingAtts(meshCaseFile, meshCase, &meshingOptions); + AttCase_setModel(meshCase, m_model); + + analysisCase = extractCase(attMngr, analysisCaseName); + pPList children = AttNode_children(analysisCase); + void* iter = 0L; + while (pANode child = static_cast(PList_next(children, &iter))) + AttCase_setModel(static_cast(child), m_model); + PList_delete(children); + } + + private: + void setMeshSize(pGModel model, pACase& meshCase, MeshAttributes& MeshAtt) { + // Set mesh size on vertices, edges, surfaces and regions + const char* element_name[4] = {"vertex", "edge", "face", "region"}; + const ElementType element_type[4] = {ElementType::vertex, ElementType::edge, ElementType::face, + ElementType::region}; + + for (int element_type_id = 0; element_type_id < 4; element_type_id++) { + std::list, double>> lpair_lId_MSize = + MeshAtt.getMSizeList(element_type[element_type_id]); + for (const auto& pair_lId_MSize : lpair_lId_MSize) { + double MSize = pair_lId_MSize.second; + std::list lElements = pair_lId_MSize.first; + for (const auto& element_id : lElements) { + pGEntity entity = GM_entityByTag(model, element_type_id, element_id); + if (entity == nullptr) { + logError() << element_name[element_type_id] << "id:" << element_id + << "not found in model."; + } else { + MS_setMeshSize(meshCase, entity, 1, MSize, nullptr); + logInfo(PMU_rank()) << element_name[element_type_id] << "id:" << element_id + << ", MSize =" << MSize; + } + } + } + } + } + + void setCases(pGModel model, pACase& meshCase, pACase& analysisCase, MeshAttributes& MeshAtt, + AnalysisAttributes& AnalysisAtt, std::unordered_map groupMap) { + + logInfo(PMU_rank()) << "Setting cases"; + // ------------------------------ Set boundary conditions + // ------------------------------ + + // Create a new manager. The manager is responsible for creating + // new attributes, saving/retrieving to/from file, and to look up + // cases, AttModels etc. + pAManager attMngr = AMAN_new(); + + // Create a case. A case serves as a grouping mechanism. It will (should) + // contain all the attributes that make up a complete analysis of + // some kind, although the system has no way of verifying this. + analysisCase = AMAN_newCase(attMngr, "analysis", "", (pModel)model); + + // Now we create an attribute information node. This can be seen + // as an "attribute generator object", as we will later on create an + // actual attribute for each geometric face from this attribute + // information node. We name the attribute T1, and give it the + // information type "boundaryCondition". + + // With this code we can tag any BC, not only 1,3,5 + int numFaces = GM_numFaces(model); + std::set UniquefaceBound; + for (auto const& f : AnalysisAtt.faceBound) { + UniquefaceBound.insert(f.bcType); + } + int numBC = UniquefaceBound.size(); + std::map aBC; + + for (auto const& ufb : UniquefaceBound) { + std::string strboundaryCondition = "BC"; + std::ostringstream stringstream; + stringstream << std::setw(9) << std::setfill('0') << ufb; + std::string sNumber = stringstream.str(); + + pAttInfoVoid iBC = AMAN_newAttInfoVoid(attMngr, "BC", "boundaryCondition"); + + strboundaryCondition.append(sNumber); + AttNode_setImageClass((pANode)iBC, strboundaryCondition.c_str()); + AttCase_addNode(analysisCase, (pANode)iBC); + aBC[ufb] = AttCase_newModelAssoc(analysisCase, (pANode)iBC); + } + + for (auto const& fb : AnalysisAtt.faceBound) { + // Get the face + pGEntity face = GM_entityByTag(model, 2, fb.faceID + 1); + if (face == nullptr) { + logError() << "faceid:" << fb.faceID + 1 << "not found in model."; + } else { + // Add the face to the model association. Note that we passed + // the Attribute Information Node into the Model Association + // at the time when the Model Association was created. That prepares + // the creation of the AttributeVoid on the face as soon as the + // association process is started + if (fb.bcType != 0) { + logInfo(PMU_rank()) << "faceBound[" << fb.faceID + 1 << "] =" << fb.bcType; + AMA_addGEntity(aBC[fb.bcType], face); + } + } + } + + // ------------------------------ Set meshing parameters + // ------------------------------ + + meshCase = MS_newMeshCase(model); + + // Set global mesh size + pModelItem modelDomain = GM_domain(model); + if (MeshAtt.globalMSize > 0) { + logInfo(PMU_rank()) << "globalMSize =" << MeshAtt.globalMSize; + // ( , , <1=absolute, 2=relative>, , ) + MS_setMeshSize(meshCase, modelDomain, 1, MeshAtt.globalMSize, nullptr); + } + + // Set mesh size on vertices, edges, surfaces and regions + setMeshSize(model, meshCase, MeshAtt); + + if (MeshAtt.velocityAwareRefinementSettings.isVelocityAwareRefinementOn()) { + logInfo(PMU_rank()) << "Enabling velocity aware meshing"; + easiMeshSize = EasiMeshSize(MeshAtt.velocityAwareRefinementSettings, model, groupMap); + auto easiMeshSizeFunc = [](pSizeAttData sadata, void* userdata) { + auto* easiMeshSize = static_cast(userdata); + std::array pt{}; + int haspt = SizeAttData_point(sadata, pt.data()); + if (!haspt) { + logError() << "!haspt"; + } + return easiMeshSize->getMeshSize(pt); + }; + // set the user-defined function for isotropic size + MS_setSizeAttFunc(meshCase, "setCustomMeshSize", easiMeshSizeFunc, &easiMeshSize); + // Relative anisotropic size is set for the entire model through the + // function anisoSize + MS_setMeshSize(meshCase, GM_domain(model), MS_userDefinedType | 1, 0, "setCustomMeshSize"); + } + + if (MeshAtt.gradation > 0) { + // Set gradation relative + logInfo(PMU_rank()) << "Gradation rate =" << MeshAtt.gradation; + MS_setGlobalSizeGradationRate(meshCase, MeshAtt.gradation); + } + if (MeshAtt.vol_AspectRatio > 0) { + // Set target equivolume AspectRatio + logInfo(PMU_rank()) << "Target equivolume AspectRatio =" << MeshAtt.vol_AspectRatio; + MS_setVolumeShapeMetric(meshCase, modelDomain, ShapeMetricType_AspectRatio, + MeshAtt.vol_AspectRatio); + } + if (MeshAtt.area_AspectRatio > 0) { + // Set target equiarea AspectRatio + logInfo(PMU_rank()) << "Target equiarea AspectRatio =" << MeshAtt.area_AspectRatio; +#ifdef BEFORE_SIM_11 + MS_setSurfaceShapeMetric(meshCase, modelDomain, ShapeMetricType_AspectRatio, + MeshAtt.area_AspectRatio); +#else + MS_setSurfaceShapeMetric(meshCase, modelDomain, ShapeMetricType_AspectRatio, + ElementType_Triangle, MeshAtt.area_AspectRatio); +#endif + } + for (auto& mycube : MeshAtt.lCube) { + logInfo(PMU_rank()) << "Cube mesh refinement: " << mycube.CubeMSize; + logInfo(PMU_rank()) << "Center" << mycube.CubeCenter[0] << " " << mycube.CubeCenter[1] << " " + << mycube.CubeCenter[2]; + logInfo(PMU_rank()) << "Width" << mycube.CubeWidth[0] << " " << mycube.CubeWidth[1] << " " + << mycube.CubeWidth[2]; + logInfo(PMU_rank()) << "Height" << mycube.CubeHeight[0] << " " << mycube.CubeHeight[1] << " " + << mycube.CubeHeight[2]; + logInfo(PMU_rank()) << "Depth" << mycube.CubeDepth[0] << " " << mycube.CubeDepth[1] << " " + << mycube.CubeDepth[2]; + MS_addCubeRefinement(meshCase, mycube.CubeMSize, &mycube.CubeCenter[0], &mycube.CubeWidth[0], + &mycube.CubeHeight[0], &mycube.CubeDepth[0]); + } + + if (MeshAtt.MeshSizePropagationDistance > 0.0) { + for (auto& iElem : MeshAtt.lFaceIdMeshSizePropagation) { + pGEntity face = GM_entityByTag(model, 2, iElem); + if (face == nullptr) { + logError() << "MeshSizeProp faceid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "MeshSizeProp faceid:" << iElem + << ", distance =" << MeshAtt.MeshSizePropagationDistance + << ", scaling factor =" << MeshAtt.MeshSizePropagationScalingFactor; +#ifdef BEFORE_SIM_15 + MS_setMeshSizePropagation(meshCase, face, 1, MeshAtt.MeshSizePropagationDistance, + MeshAtt.MeshSizePropagationScalingFactor); +#else + MS_setMeshSizePropagation(meshCase, face, 2, 1, MeshAtt.MeshSizePropagationDistance, + MeshAtt.MeshSizePropagationScalingFactor); +#endif + } + } + } + for (auto& iElem : MeshAtt.lFaceIdUseDiscreteMesh) { + pGEntity face = GM_entityByTag(model, 2, iElem); + if (face == nullptr) { + logError() << "UseDiscreteMesh; faceid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "UseDiscreteMesh; faceid, noModification:" << iElem + << MeshAtt.UseDiscreteMesh_noModification; + MS_useDiscreteGeometryMesh(meshCase, face, MeshAtt.UseDiscreteMesh_noModification); + // MS_limitSurfaceMeshModification(meshCase,face,UseDiscreteMesh_noModification); + } + } + for (auto& iElem : MeshAtt.lFaceIdNoMesh) { + pGEntity face = GM_entityByTag(model, 2, iElem); + if (face == nullptr) { + logError() << "No Mesh; faceid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "No Mesh; faceid:" << iElem; + MS_setNoMesh(meshCase, face, 1); + } + } + for (auto& iElem : MeshAtt.lRegionIdNoMesh) { + pGEntity region = GM_entityByTag(model, 3, iElem); + if (region == nullptr) { + logError() << "No Mesh; regionid:" << iElem << "not found in model."; + } else { + logInfo(PMU_rank()) << "No Mesh; regionid:" << iElem; + MS_setNoMesh(meshCase, region, 1); + } + } + } + + private: + void loadCAD(const char* modFile, const char* cadFile) { + // Load CAD + std::string sCadFile; + if (cadFile) + sCadFile = cadFile; + else { + sCadFile = modFile; + utils::StringUtils::replaceLast(sCadFile, ".smd", "_nat.x_t"); + } + pNativeModel nativeModel = 0L; +#ifdef PARASOLID + if (utils::Path(sCadFile).exists()) + nativeModel = ParasolidNM_createFromFile(sCadFile.c_str(), 0); +#endif + m_model = GM_load(modFile, nativeModel, 0L); + nativeModel = GM_nativeModel(m_model); + + if (nativeModel) + NM_release(nativeModel); + +#ifdef BEFORE_SIM_18 + // check for model errors + pPList modelErrors = PList_new(); + if (!GM_isValid(m_model, 1, modelErrors)) { + void* iter = 0L; + while (pSimError err = static_cast(PList_next(modelErrors, &iter))) { + logInfo(PMU_rank()) << " Error code: " << SimError_code(err) << std::endl; + logInfo(PMU_rank()) << " Error string: " << SimError_toString(err) << std::endl; + } + logError() << "Input model is not valid"; + } + PList_delete(modelErrors); +#else + // check for model infos + pPList modelInfos = PList_new(); + if (!GM_isValid(m_model, 1, modelInfos)) { + void* iter = 0L; + while (pSimInfo info = static_cast(PList_next(modelInfos, &iter))) { + logInfo(PMU_rank()) << " Info code: " << SimInfo_code(info) << std::endl; + logInfo(PMU_rank()) << " Info string: " << SimInfo_toString(info) << std::endl; + } + logError() << "Input model is not valid"; + } + PList_delete(modelInfos); +#endif + + // check for self-intersecting CAD mesh + pProgress prog = Progress_new(); + Progress_setCallback(prog, progressHandler); + pMesh mesh = DM_getMesh(static_cast(m_model)); + pPList entities = PList_new(); + if (MS_checkMeshIntersections(mesh, entities, prog) > 0) { + + const char* element_name[4] = {"vertex", "edge", "face", "region"}; + + for (int i = 0; i < PList_size(entities); i += 2) { + pEntity firstEnt = (pEntity)PList_item(entities, i); + pEntity secondEnt = (pEntity)PList_item(entities, i + 1); + logInfo(PMU_rank()) << "Self-intersection between" << element_name[EN_whatInType(firstEnt)] + << GEN_tag(EN_whatIn(firstEnt)) << "and" + << element_name[EN_whatInType(secondEnt)] + << GEN_tag(EN_whatIn(secondEnt)); + double centroid[3], centroid2[3]; + EN_centroid(firstEnt, ¢roid[0]); + EN_centroid(secondEnt, ¢roid2[0]); + logInfo(PMU_rank()) << "centroids" << centroid[0] << centroid[1] << centroid[2] << "and " + << centroid2[0] << centroid2[1] << centroid2[2] << std::flush; + } + logError() << PList_size(entities) / 2 << "Self-intersection(s) detected in CAD mesh"; + } + PList_delete(entities); + Progress_delete(prog); + } + + private: + void analyse_mesh() { + int num_bins = 8; + double AR[8] = {0, 2, 4, 6, 10, 20, 40, 100}; + long int AR_vol_bins[num_bins]; + long int skew_area_bins[num_bins]; + memset(AR_vol_bins, 0, num_bins * sizeof(long int)); + memset(skew_area_bins, 0, num_bins * sizeof(long int)); + + // Accumulate equivolume AspectRatio data + RIter reg_it; + int num_partMeshes = PM_numParts(m_simMesh); + pRegion reg; + double maxAR = 0, AR_global = 0, elAR; + int kAR; + for (int i = 0; i < num_partMeshes; i++) { + reg_it = M_regionIter(PM_mesh(m_simMesh, i)); + while ((reg = RIter_next(reg_it))) { + kAR = num_bins; + for (int k = 1; k < num_bins; k++) { + elAR = R_aspectRatio(reg); + if (elAR < AR[k]) { + kAR = k; + break; + } + } + AR_vol_bins[kAR - 1]++; + maxAR = std::max(maxAR, elAR); + } + } + RIter_delete(reg_it); + + // Print the statistics + logInfo(PMU_rank()) << "AR statistics:"; + MPI_Allreduce(&maxAR, &AR_global, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + logInfo(PMU_rank()) << "AR max:" << AR_global; + logInfo(PMU_rank()) << "AR (target: < ~10):"; + long int bin_global; + for (int i = 0; i < num_bins - 1; i++) { + MPI_Allreduce(&AR_vol_bins[i], &bin_global, 1, MPI_LONG, MPI_SUM, MPI_COMM_WORLD); + logInfo(PMU_rank()) << std::fixed << std::setprecision(2) << "[" << AR[i] << "," << AR[i + 1] + << "):" << bin_global; + } + MPI_Allreduce(&AR_vol_bins[num_bins - 1], &bin_global, 1, MPI_LONG, MPI_SUM, MPI_COMM_WORLD); + logInfo(PMU_rank()) << std::fixed << std::setprecision(2) << "[" << AR[num_bins - 1] + << ",inf):" << bin_global; + } +}; + +#endif // SIM_MOD_SUITE_H diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 45dd8ba..e202192 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -49,6 +49,8 @@ #include "meshreader/ParallelGambitReader.h" #include "third_party/GMSH2Parser.h" +#include "aux/InsphereCalculator.h" + template static TT _checkH5Err(TT&& status, const char* file, int line) { if (status < 0) { logError() << utils::nospace << "An HDF5 error occurred (" << file << ": " << line << ")"; @@ -78,7 +80,7 @@ template static void iterate(apf::Mesh* mesh, int dim, F&& function constexpr std::size_t NoSecondDim = 0; template -static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf::Mesh* mesh, +static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, void* mesh, int meshdim, hid_t h5memtype, hid_t h5outtype, std::size_t chunk, std::size_t localSize, std::size_t globalSize, bool reduceInts, int filterEnable, std::size_t filterChunksize) { @@ -142,8 +144,6 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: std::size_t written = 0; std::size_t index = 0; - apf::MeshIterator* it = mesh->begin(meshdim); - hsize_t nullstart[2] = {0, 0}; for (std::size_t i = 0; i < rounds; ++i) { @@ -151,15 +151,7 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: start[0] = offset + written; count[0] = std::min(localSize - written, bufferSize); - while (apf::MeshEntity* element = mesh->iterate(it)) { - if (handler(element, data.begin() + index * SecondSize)) { - ++index; - } - - if (index >= count[0]) { - break; - } - } + std::copy_n(handler.begin() + written * SecondSize, count[0] * SecondSize, data.begin()); checkH5Err(H5Sselect_hyperslab(h5memspace, H5S_SELECT_SET, nullstart, nullptr, count, nullptr)); @@ -170,8 +162,6 @@ static void writeH5Data(F&& handler, hid_t h5file, const std::string& name, apf: written += count[0]; } - mesh->end(it); - if (filterEnable > 0) { checkH5Err(H5Pclose(h5filter)); } @@ -306,24 +296,23 @@ int main(int argc, char* argv[]) { } // Create/read the mesh - MeshInput* meshInput = 0L; - apf::Mesh2* mesh = 0L; + MeshData* meshInput = nullptr; switch (args.getArgument("source", 0)) { case 0: logInfo(rank) << "Using Gambit mesh"; - meshInput = new SerialMeshFile(inputFile); + // meshInput = new SerialMeshFile(inputFile); break; case 1: logInfo(rank) << "Using Fidap mesh"; - meshInput = new SerialMeshFile(inputFile); + // meshInput = new SerialMeshFile(inputFile); break; case 2: logInfo(rank) << "Using GMSH mesh format 2 (msh2) mesh"; - meshInput = new SerialMeshFile>(inputFile); + // meshInput = new SerialMeshFile>(inputFile); break; case 3: logInfo(rank) << "Using GMSH mesh format 4 (msh4) mesh"; - meshInput = new SerialMeshFile>(inputFile); + // meshInput = new SerialMeshFile>(inputFile); break; case 4: #ifdef USE_NETCDF @@ -335,7 +324,7 @@ int main(int argc, char* argv[]) { break; case 5: logInfo(rank) << "Using APF native format"; - meshInput = new ApfNative(inputFile, args.getArgument("model", 0L)); + // meshInput = new ApfNative(inputFile, args.getArgument("model", 0L)); break; case 6: #ifdef USE_SIMMOD @@ -355,67 +344,24 @@ int main(int argc, char* argv[]) { logError() << "Unknown source"; } - mesh = meshInput->getMesh(); - logInfo(rank) << "Parsed mesh successfully, writing output..."; - // Check mesh - if (alignMdsMatches(mesh)) - logWarning() << "Fixed misaligned matches"; - mesh->verify(); - - // Dump mesh for later usage? - const char* dumpFile = args.getArgument("dump", 0L); - if (dumpFile) { - logInfo(PCU_Comm_Self()) << "Writing native APF mesh"; - mesh->writeNative(dumpFile); - - const char* modelFile = args.getArgument("model", 0L); - if (modelFile) - gmi_write_dmg(mesh->getModel(), modelFile); - } - - // Dump VTK mesh - const char* vtkPrefix = args.getArgument("vtk", 0L); - if (vtkPrefix) { - logInfo(PCU_Comm_Self()) << "Writing VTK mesh"; - apf::writeVtkFiles(vtkPrefix, mesh); - } - - apf::Sharing* sharing = apf::getSharing(mesh); - - // oriented at the apf::countOwned method ... But with size_t - auto countOwnedLong = [sharing, mesh](int dim) { - std::size_t counter = 0; - - apf::MeshIterator* it = mesh->begin(dim); - while (apf::MeshEntity* element = mesh->iterate(it)) { - if (sharing->isOwned(element)) { - ++counter; - } - } - mesh->end(it); - - return counter; - }; + void* mesh = nullptr; // TODO(David): replace by apf::countOwned again once extended // Get local/global size - std::size_t localSize[2] = {countOwnedLong(3), countOwnedLong(0)}; + std::size_t localSize[2] = {meshInput->cellCount(), meshInput->vertexCount()}; std::size_t globalSize[2] = {localSize[0], localSize[1]}; MPI_Allreduce(MPI_IN_PLACE, globalSize, 2, tndm::mpi_type_t(), MPI_SUM, MPI_COMM_WORLD); - logInfo(rank) << "Mesh size:" << globalSize[0]; + logInfo(rank) << "Total cell count:" << globalSize[0]; + logInfo(rank) << "Total vertex count:" << globalSize[1]; - // Compute min insphere radius - double min = std::numeric_limits::max(); - apf::MeshIterator* it = mesh->begin(3); - while (apf::MeshEntity* element = mesh->iterate(it)) { - min = std::min(min, ma::getInsphere(mesh, element)); - } - mesh->end(it); + auto inspheres = + calculateInsphere(meshInput->connectivity(), meshInput->geometry(), MPI_COMM_WORLD); + double min = *std::min_element(inspheres.begin(), inspheres.end()); MPI_Reduce((rank == 0 ? MPI_IN_PLACE : &min), &min, 1, MPI_DOUBLE, MPI_MIN, 0, MPI_COMM_WORLD); logInfo(rank) << "Minimum insphere found:" << min; @@ -425,10 +371,6 @@ int main(int argc, char* argv[]) { offsets[0] -= localSize[0]; offsets[1] -= localSize[1]; - // Create numbering for the vertices/elements - apf::GlobalNumbering* vertexNum = apf::makeGlobal(apf::numberOwnedNodes(mesh, "vertices")); - apf::synchronize(vertexNum); - // Create the H5 file hid_t h5falist = H5Pcreate(H5P_FILE_ACCESS); checkH5Err(h5falist); @@ -445,99 +387,35 @@ int main(int argc, char* argv[]) { // Write cells std::size_t connectSize = 8; logInfo(rank) << "Writing cells"; - writeH5Data( - [&](auto& element, auto&& data) { - apf::NewArray vn; - apf::getElementNumbers(vertexNum, element, vn); - - for (int i = 0; i < 4; i++) { - *(data + i) = vn[i]; - } - return true; - }, - h5file, "connect", mesh, 3, H5T_NATIVE_ULONG, H5T_STD_U64LE, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->connectivity(), h5file, "connect", mesh, 3, + H5T_NATIVE_ULONG, H5T_STD_U64LE, chunksize, localSize[0], + globalSize[0], reduceInts, filterEnable, filterChunksize); // Vertices logInfo(rank) << "Writing vertices"; - writeH5Data( - [&](auto& element, auto&& data) { - if (!sharing->isOwned(element)) { - return false; - } - - /*long gid = apf::getNumber(vertexNum, apf::Node(element, 0)); - - if (gid != static_cast(offsets[1] + index)) { - logError() << "Global vertex numbering is incorrect"; - }*/ - - apf::Vector3 point; - mesh->getPoint(element, 0, point); - double geometry[3]; - point.toArray(geometry); - - *(data + 0) = geometry[0]; - *(data + 1) = geometry[1]; - *(data + 2) = geometry[2]; - - return true; - }, - h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, H5T_IEEE_F64LE, chunksize, localSize[1], - globalSize[1], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->geometry(), h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, + H5T_IEEE_F64LE, chunksize, localSize[1], globalSize[1], reduceInts, + filterEnable, filterChunksize); // Group information - apf::MeshTag* groupTag = mesh->findTag("group"); std::size_t groupSize = 4; - if (groupTag) { + if (true) { logInfo(rank) << "Writing group information"; - writeH5Data( - [&](auto& element, auto&& data) { - int group; - mesh->getIntTag(element, groupTag, &group); - *data = group; - return true; - }, - h5file, "group", mesh, 3, H5T_NATIVE_INT, H5T_STD_I32LE, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->group(), h5file, "group", mesh, 3, H5T_NATIVE_INT, + H5T_STD_I32LE, chunksize, localSize[0], globalSize[0], reduceInts, + filterEnable, filterChunksize); } else { logInfo() << "No group information found in mesh"; } // Write boundary condition logInfo(rank) << "Writing boundary condition"; - apf::MeshTag* boundaryTag = mesh->findTag("boundary condition"); - assert(boundaryTag); auto i32limit = std::numeric_limits::max(); auto i64limit = std::numeric_limits::max(); - writeH5Data( - [&](auto& element, auto&& data) { - long boundary = 0; - apf::Downward faces; - mesh->getDownward(element, 2, faces); - - for (int i = 0; i < 4; i++) { - if (mesh->hasTag(faces[i], boundaryTag)) { - int b; - mesh->getIntTag(faces[i], boundaryTag, &b); - - if (b <= 0 || (b > i32limit && boundaryType == 0) || - (b > i64limit && boundaryType == 1)) { - logError() << "Cannot handle boundary condition" << b; - } - - // NOTE: this fails for boundary values per face larger than 2**31 - 1 due to sign - // extension - boundary |= static_cast(b) << (i * faceOffset); - } - } - - *data = boundary; - return true; - }, - h5file, "boundary", mesh, 3, H5T_NATIVE_LONG, boundaryDatatype, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->boundary(), h5file, "boundary", mesh, 3, + H5T_NATIVE_LONG, boundaryDatatype, chunksize, localSize[0], + globalSize[0], reduceInts, filterEnable, filterChunksize); // Writing XDMF file if (rank == 0) { @@ -590,8 +468,6 @@ int main(int argc, char* argv[]) { checkH5Err(H5Fclose(h5file)); - delete sharing; - delete mesh; delete meshInput; logInfo(rank) << "Finished successfully"; From 151c9219a1facb2b8f70bf23a5ac3402a507c1c9 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Sat, 28 Oct 2023 16:35:33 +0200 Subject: [PATCH 03/34] Allow larger vertex arrays --- src/aux/InsphereCalculator.cpp | 4 ++-- src/input/ParallelVertexFilter.h | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/aux/InsphereCalculator.cpp b/src/aux/InsphereCalculator.cpp index 6498a8f..f2910a6 100644 --- a/src/aux/InsphereCalculator.cpp +++ b/src/aux/InsphereCalculator.cpp @@ -43,8 +43,8 @@ std::vector calculateInsphere(const std::vector& connectivi MPI_Alltoall(outrequests.data(), 1, MPI_INT, inrequests.data(), 1, MPI_INT, comm); - std::vector outdisp(commsize + 1); - std::vector indisp(commsize + 1); + std::vector outdisp(commsize + 1); + std::vector indisp(commsize + 1); for (std::size_t i = 1; i < commsize + 1; ++i) { outdisp[i] = outdisp[i - 1] + outrequests[i - 1]; indisp[i] = indisp[i - 1] + inrequests[i - 1]; diff --git a/src/input/ParallelVertexFilter.h b/src/input/ParallelVertexFilter.h index 5a34909..360c5b8 100644 --- a/src/input/ParallelVertexFilter.h +++ b/src/input/ParallelVertexFilter.h @@ -194,8 +194,8 @@ class ParallelVertexFilter { // Allocate buffer for the vertices and exchange them std::vector sortVertices(3 * numSortVertices); - std::vector sDispls(m_numProcs); - std::vector rDispls(m_numProcs); + std::vector sDispls(m_numProcs); + std::vector rDispls(m_numProcs); sDispls[0] = 0; rDispls[0] = 0; for (int i = 1; i < m_numProcs; i++) { From cef998a6885df61f8776834281fa97167d5151c3 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Sat, 28 Oct 2023 16:36:33 +0200 Subject: [PATCH 04/34] Make Serial Mesh files work without APF --- src/aux/Distributor.h | 12 ++ src/input/SerialMeshFile.h | 86 ++++----------- src/input/SerialMeshFileApf.h | 151 ++++++++++++++++++++++++++ src/meshreader/GMSHBuilder.h | 2 +- src/meshreader/ParallelGMSHReader.h | 8 +- src/meshreader/ParallelGambitReader.h | 21 +--- src/pumgen.cpp | 8 +- 7 files changed, 196 insertions(+), 92 deletions(-) create mode 100644 src/aux/Distributor.h create mode 100644 src/input/SerialMeshFileApf.h diff --git a/src/aux/Distributor.h b/src/aux/Distributor.h new file mode 100644 index 0000000..b5918a9 --- /dev/null +++ b/src/aux/Distributor.h @@ -0,0 +1,12 @@ +#ifndef PUMGEN_AUX_DISTRIBUTOR_H_ +#define PUMGEN_AUX_DISTRIBUTOR_H_ + +#include + +// no namespace for now + +constexpr std::size_t getChunksize(std::size_t total, int rank, int size) { + return (total / size) + (total % size > rank); +} + +#endif // PUMGEN_AUX_DISTRIBUTOR_H_ diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 5377314..3332710 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -18,18 +18,20 @@ #endif // PARALLEL #include "ApfConvertWrapper.h" +#include "aux/Distributor.h" #include #include #include #include #include -#include "MeshInput.h" +#include "MeshData.h" +#include "utils/logger.h" /** * Read a mesh from a serial file */ -template class SerialMeshFile : public MeshInput { +template class SerialMeshFile : public FullStorageMeshData { private: #ifdef PARALLEL MPI_Comm m_comm; @@ -73,77 +75,27 @@ template class SerialMeshFile : public MeshInput { const std::size_t nVertices = m_meshReader.nVertices(); const std::size_t nElements = m_meshReader.nElements(); - std::size_t nLocalVertices = (nVertices + m_nProcs - 1) / m_nProcs; - std::size_t nLocalElements = (nElements + m_nProcs - 1) / m_nProcs; - if (m_rank == m_nProcs - 1) { - nLocalVertices = nVertices - (m_nProcs - 1) * nLocalVertices; - nLocalElements = nElements - (m_nProcs - 1) * nLocalElements; - } + const std::size_t nLocalVertices = getChunksize(nVertices, m_rank, m_nProcs); + const std::size_t nLocalElements = getChunksize(nElements, m_rank, m_nProcs); - gmi_register_null(); - gmi_model* model = gmi_load(".null"); - m_mesh = apf::makeEmptyMdsMesh(model, 3, false); - - // Create elements - apf::GlobalToVert vertMap; - { - std::vector elements; - { - std::vector elementsInt(nLocalElements * 4); - m_meshReader.readElements(elementsInt.data()); - elements = std::vector(elementsInt.begin(), elementsInt.end()); - } - logInfo(m_rank) << "Create APF connectivity"; - apf::construct(m_mesh, elements.data(), nLocalElements, apf::Mesh::TET, vertMap); - } + setup(nLocalElements, nLocalVertices); - apf::alignMdsRemotes(m_mesh); - apf::deriveMdsModel(m_mesh); + logInfo(m_rank) << "Read vertex coordinates"; + m_meshReader.readVertices(geometryData.data()); - // Set vertices - { - std::vector vertices(nLocalVertices * 3); - m_meshReader.readVertices(vertices.data()); - logInfo(m_rank) << "Set coordinates in APF"; - apf::setCoords(m_mesh, vertices.data(), nLocalVertices, vertMap); - } + logInfo(m_rank) << "Read cell vertices"; + m_meshReader.readElements(connectivityData.data()); - // Set boundaries - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - { - std::vector boundaries(nLocalElements * 4); - m_meshReader.readBoundaries(boundaries.data()); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::Adjacent adjacent; - m_mesh->getAdjacent(element, 2, adjacent); - - for (int j = 0; j < 4; j++) { - if (!boundaries[i * 4 + j]) { - continue; - } - - m_mesh->setIntTag(adjacent[j], boundaryTag, &boundaries[i * 4 + j]); - } - - i++; - } - m_mesh->end(it); - } + logInfo(m_rank) << "Read cell groups"; + m_meshReader.readGroups(groupData.data()); - // Set groups - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - { - std::vector groups(nLocalElements); - m_meshReader.readGroups(groups.data()); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - m_mesh->setIntTag(element, groupTag, &groups[i]); - i++; + logInfo(m_rank) << "Read boundary conditions"; + std::vector preBoundaryData(nLocalElements * 4); + m_meshReader.readBoundaries(preBoundaryData.data()); + for (std::size_t i = 0; i < nLocalElements; ++i) { + for (int j = 0; j < 4; ++j) { + setBoundary(i, j, preBoundaryData[4 * i + j]); } - m_mesh->end(it); } } }; diff --git a/src/input/SerialMeshFileApf.h b/src/input/SerialMeshFileApf.h new file mode 100644 index 0000000..10ed374 --- /dev/null +++ b/src/input/SerialMeshFileApf.h @@ -0,0 +1,151 @@ +/** + * @file + * This file is part of PUMGen + * + * For conditions of distribution and use, please see the copyright + * notice in the file 'COPYING' at the root directory of this package + * and the copyright notice at https://github.com/SeisSol/PUMGen + * + * @copyright 2017 Technical University of Munich + * @author Sebastian Rettenberger + */ + +#ifndef SERIAL_MESH_FILE_H +#define SERIAL_MESH_FILE_H + +#ifdef PARALLEL +#include +#endif // PARALLEL + +#include "ApfConvertWrapper.h" +#include "aux/Distributor.h" +#include +#include +#include +#include +#include + +#include "MeshInput.h" + +/** + * Read a mesh from a serial file + */ +template class SerialMeshFile : public MeshInput { + private: +#ifdef PARALLEL + MPI_Comm m_comm; +#endif // PARALLEL + + int m_rank; + int m_nProcs; + + T m_meshReader; + + public: +#ifdef PARALLEL + SerialMeshFile(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) + : m_comm(comm), m_meshReader(comm) { + init(); + open(meshFile); + } +#else // PARALLEL + SerialMeshFile(const char* meshFile) { + init(); + open(meshFile); + } +#endif // PARALLEL + + private: + /** + * Sets some parameters (called from the constructor) + */ + void init() { +#ifdef PARALLEL + MPI_Comm_rank(m_comm, &m_rank); + MPI_Comm_size(m_comm, &m_nProcs); +#else // PARALLLEL + m_rank = 0; + m_nProcs = 1; +#endif // PARALLEL + } + + void open(const char* meshFile) { + m_meshReader.open(meshFile); + + const std::size_t nVertices = m_meshReader.nVertices(); + const std::size_t nElements = m_meshReader.nElements(); + const std::size_t nLocalVertices = getChunksize(nVertices, m_rank, m_nProcs); + const std::size_t nLocalElements = getChunksize(nElements, m_rank, m_nProcs); + + gmi_register_null(); + gmi_model* model = gmi_load(".null"); + m_mesh = apf::makeEmptyMdsMesh(model, 3, false); + + // Create elements + logInfo(m_rank) << "Handle cell connectivity"; + apf::GlobalToVert vertMap; + { + std::vector elements; + { + std::vector elementsInt(nLocalElements * 4); + m_meshReader.readElements(elementsInt.data()); + elements = std::vector(elementsInt.begin(), elementsInt.end()); + } + logInfo(m_rank) << "Create APF connectivity"; + apf::construct(m_mesh, elements.data(), nLocalElements, apf::Mesh::TET, vertMap); + } + + apf::alignMdsRemotes(m_mesh); + apf::deriveMdsModel(m_mesh); + + // Set vertices + logInfo(m_rank) << "Handle vertex coordinates"; + { + std::vector vertices(nLocalVertices * 3); + m_meshReader.readVertices(vertices.data()); + apf::setCoords(m_mesh, vertices.data(), nLocalVertices, vertMap); + } + + // Set boundaries + logInfo(m_rank) << "Handle boundary conditions"; + apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); + { + std::vector boundaries(nLocalElements * 4); + m_meshReader.readBoundaries(boundaries.data()); + apf::MeshIterator* it = m_mesh->begin(3); + std::size_t i = 0; + while (apf::MeshEntity* element = m_mesh->iterate(it)) { + apf::Adjacent adjacent; + m_mesh->getAdjacent(element, 2, adjacent); + + for (int j = 0; j < 4; j++) { + if (!boundaries[i * 4 + j]) { + continue; + } + + m_mesh->setIntTag(adjacent[j], boundaryTag, &boundaries[i * 4 + j]); + } + + i++; + } + m_mesh->end(it); + } + + // Set groups + logInfo(m_rank) << "Handle cell groups"; + apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); + { + std::vector groups(nLocalElements); + m_meshReader.readGroups(groups.data()); + apf::MeshIterator* it = m_mesh->begin(3); + std::size_t i = 0; + while (apf::MeshEntity* element = m_mesh->iterate(it)) { + m_mesh->setIntTag(element, groupTag, &groups[i]); + i++; + } + m_mesh->end(it); + } + } +}; + +#endif // SERIAL_MESH_FILE_H diff --git a/src/meshreader/GMSHBuilder.h b/src/meshreader/GMSHBuilder.h index 1224920..4bc11be 100644 --- a/src/meshreader/GMSHBuilder.h +++ b/src/meshreader/GMSHBuilder.h @@ -36,7 +36,7 @@ template class GMSHBuilder : public tndm::GMSHMeshBuilder { std::vector bcs; void setNumVertices(std::size_t numVertices) { vertices.resize(numVertices); } - void setVertex(long id, std::array const& x) { + void setVertex(long id, std::array const& x) { for (std::size_t i = 0; i < D; ++i) { vertices[id][i] = x[i]; } diff --git a/src/meshreader/ParallelGMSHReader.h b/src/meshreader/ParallelGMSHReader.h index 66b2902..eed6677 100644 --- a/src/meshreader/ParallelGMSHReader.h +++ b/src/meshreader/ParallelGMSHReader.h @@ -11,6 +11,8 @@ #include #include +#include "aux/Distributor.h" + namespace puml { template class ParallelGMSHReader { @@ -40,7 +42,7 @@ template class ParallelGMSHReader { } MPI_Bcast(&nVertices_, 1, tndm::mpi_type_t(), 0, comm_); - MPI_Bcast(&nElements_, 1, tndm::mpi_type_t(), 0, comm_); + MPI_Bcast(&nElements_, 1, tndm::mpi_type_t(), 0, comm_); } [[nodiscard]] std::size_t nVertices() const { return nVertices_; } @@ -140,13 +142,11 @@ template class ParallelGMSHReader { MPI_Comm_rank(comm_, &rank); MPI_Comm_size(comm_, &procs); - const std::size_t chunkSize = (numElements + procs - 1) / procs; auto sendcounts = std::vector(procs); auto displs = std::vector(procs + 1); - std::fill(sendcounts.begin(), sendcounts.end(), numPerElement * chunkSize); - sendcounts.back() = numPerElement * (numElements - (procs - 1u) * chunkSize); displs[0] = 0; for (int i = 0; i < procs; ++i) { + sendcounts[i] = numPerElement * getChunksize(numElements, i, procs); displs[i + 1] = displs[i] + sendcounts[i]; } diff --git a/src/meshreader/ParallelGambitReader.h b/src/meshreader/ParallelGambitReader.h index 2bfc970..b5b3d8c 100644 --- a/src/meshreader/ParallelGambitReader.h +++ b/src/meshreader/ParallelGambitReader.h @@ -17,6 +17,7 @@ #include "GambitReader.h" #include "ParallelMeshReader.h" +#include "aux/Distributor.h" #include "third_party/MPITraits.h" @@ -37,7 +38,7 @@ class ParallelGambitReader : public ParallelMeshReader { * This is a collective operation. */ void readGroups(int* groups) { - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; + std::size_t chunkSize = getChunksize(nElements(), m_rank, m_nProcs); if (m_rank == 0) { std::size_t maxChunkSize = chunkSize; @@ -50,9 +51,7 @@ class ParallelGambitReader : public ParallelMeshReader { for (int i = 0; i < m_nProcs; i++) { logInfo() << "Reading group information part" << (i + 1) << "of" << m_nProcs; - if (i == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } + chunkSize = getChunksize(nElements(), i, m_nProcs); m_serialReader.readGroups(i * maxChunkSize, chunkSize, map.data()); @@ -92,10 +91,6 @@ class ParallelGambitReader : public ParallelMeshReader { MPI_Waitall((m_nProcs - 1) * 2, requests.data(), MPI_STATUSES_IGNORE); } else { - if (m_rank == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } - // Allocate enough memory std::vector buf(chunkSize * 2); @@ -126,7 +121,7 @@ class ParallelGambitReader : public ParallelMeshReader { * @todo Only tetrahedral meshes are supported */ void readBoundaries(int* boundaries) { - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; + std::size_t chunkSize = getChunksize(nBoundaries(), m_rank, m_nProcs); if (m_rank == 0) { std::size_t maxChunkSize = chunkSize; @@ -140,9 +135,7 @@ class ParallelGambitReader : public ParallelMeshReader { for (std::size_t i = 0; i < nChunks; i++) { logInfo() << "Reading boundary conditions part" << (i + 1) << "of" << nChunks; - if (i == nChunks - 1) { - chunkSize = nBoundaries() - (nChunks - 1) * chunkSize; - } + chunkSize = getChunksize(nBoundaries(), i, m_nProcs); m_serialReader.readBoundaries(i * maxChunkSize, chunkSize, faces.data()); @@ -188,10 +181,6 @@ class ParallelGambitReader : public ParallelMeshReader { } MPI_Waitall(m_nProcs - 1, requests.data(), MPI_STATUSES_IGNORE); } else { - if (m_rank == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } - // Allocate enough memory std::vector buf(chunkSize * 2); diff --git a/src/pumgen.cpp b/src/pumgen.cpp index e202192..b3b7fc3 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -300,19 +300,19 @@ int main(int argc, char* argv[]) { switch (args.getArgument("source", 0)) { case 0: logInfo(rank) << "Using Gambit mesh"; - // meshInput = new SerialMeshFile(inputFile); + meshInput = new SerialMeshFile(inputFile); break; case 1: logInfo(rank) << "Using Fidap mesh"; - // meshInput = new SerialMeshFile(inputFile); + meshInput = new SerialMeshFile(inputFile); break; case 2: logInfo(rank) << "Using GMSH mesh format 2 (msh2) mesh"; - // meshInput = new SerialMeshFile>(inputFile); + meshInput = new SerialMeshFile>(inputFile); break; case 3: logInfo(rank) << "Using GMSH mesh format 4 (msh4) mesh"; - // meshInput = new SerialMeshFile>(inputFile); + meshInput = new SerialMeshFile>(inputFile); break; case 4: #ifdef USE_NETCDF From 86c0f5227693c9191bca5df71a7be8fee738544c Mon Sep 17 00:00:00 2001 From: David Schneller Date: Sat, 28 Oct 2023 17:00:17 +0200 Subject: [PATCH 05/34] Improve insphere calculation --- src/aux/InsphereCalculator.cpp | 70 +++++++++++++++++++++------------- 1 file changed, 43 insertions(+), 27 deletions(-) diff --git a/src/aux/InsphereCalculator.cpp b/src/aux/InsphereCalculator.cpp index f2910a6..2845565 100644 --- a/src/aux/InsphereCalculator.cpp +++ b/src/aux/InsphereCalculator.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include "MPIConvenience.h" @@ -35,10 +36,19 @@ std::vector calculateInsphere(const std::vector& connectivi vertexDist[i + 1] += vertexDist[i]; } + std::vector> outidxmap(commsize); + for (const auto& vertex : connectivity) { auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); auto position = std::distance(vertexDist.begin(), itPosition) - 1; - ++outrequests[position]; + auto localVertex = vertex - vertexDist[position]; + + // only transfer each vertex coordinate once, and only if it's not already on the same rank + if (position != commrank && + outidxmap[position].find(localVertex) == outidxmap[position].end()) { + outidxmap[position][localVertex] = outrequests[position]; + ++outrequests[position]; + } } MPI_Alltoall(outrequests.data(), 1, MPI_INT, inrequests.data(), 1, MPI_INT, comm); @@ -53,16 +63,14 @@ std::vector calculateInsphere(const std::vector& connectivi std::vector inidx(indisp[commsize]); { - std::vector outidx(connectivity.size()); + std::vector outidx(outdisp[commsize]); std::vector counter(commsize); - for (std::size_t i = 0; i < connectivity.size(); ++i) { - auto vertex = connectivity[i]; - auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); - auto position = std::distance(vertexDist.begin(), itPosition) - 1; - outidx[counter[position] + outdisp[position]] = vertex - vertexDist[position]; - ++counter[position]; + for (std::size_t i = 0; i < commsize; ++i) { + for (const auto& [localVertex, j] : outidxmap[i]) { + outidx[j + outdisp[i]] = localVertex; + } } // transfer indices @@ -72,7 +80,6 @@ std::vector calculateInsphere(const std::vector& connectivi tndm::mpi_type_t(), comm); } - // this is a bit inefficient, since we don't look for duplicates... TODO: maybe improve std::vector outvertices(3 * connectivity.size()); { @@ -94,29 +101,38 @@ std::vector calculateInsphere(const std::vector& connectivi std::vector inspheres(connectivity.size() / 4); for (std::size_t i = 0; i < connectivity.size() / 4; ++i) { - std::array vertexIds{0, 0, 0, 0}; + std::array, 4> vertices; for (int j = 0; j < 4; ++j) { auto vertex = connectivity[i * 4 + j]; auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex); auto position = std::distance(vertexDist.begin(), itPosition) - 1; - vertexIds[j] = counter[position] + outdisp[position]; - ++counter[position]; + auto localVertex = vertex - vertexDist[position]; + if (position == commrank) { + vertices[j][0] = geometry[localVertex * 3 + 0]; + vertices[j][1] = geometry[localVertex * 3 + 1]; + vertices[j][2] = geometry[localVertex * 3 + 2]; + } else { + auto transferidx = outidxmap[position][localVertex] + outdisp[position]; + vertices[j][0] = outvertices[transferidx * 3 + 0]; + vertices[j][1] = outvertices[transferidx * 3 + 1]; + vertices[j][2] = outvertices[transferidx * 3 + 2]; + } } - double a11 = outvertices[3 * vertexIds[1] + 0] - outvertices[3 * vertexIds[0] + 0]; - double a12 = outvertices[3 * vertexIds[1] + 1] - outvertices[3 * vertexIds[0] + 1]; - double a13 = outvertices[3 * vertexIds[1] + 2] - outvertices[3 * vertexIds[0] + 2]; - double a21 = outvertices[3 * vertexIds[2] + 0] - outvertices[3 * vertexIds[0] + 0]; - double a22 = outvertices[3 * vertexIds[2] + 1] - outvertices[3 * vertexIds[0] + 1]; - double a23 = outvertices[3 * vertexIds[2] + 2] - outvertices[3 * vertexIds[0] + 2]; - double a31 = outvertices[3 * vertexIds[3] + 0] - outvertices[3 * vertexIds[0] + 0]; - double a32 = outvertices[3 * vertexIds[3] + 1] - outvertices[3 * vertexIds[0] + 1]; - double a33 = outvertices[3 * vertexIds[3] + 2] - outvertices[3 * vertexIds[0] + 2]; - double b11 = outvertices[3 * vertexIds[2] + 0] - outvertices[3 * vertexIds[1] + 0]; - double b12 = outvertices[3 * vertexIds[2] + 1] - outvertices[3 * vertexIds[1] + 1]; - double b13 = outvertices[3 * vertexIds[2] + 2] - outvertices[3 * vertexIds[1] + 2]; - double b21 = outvertices[3 * vertexIds[3] + 0] - outvertices[3 * vertexIds[1] + 0]; - double b22 = outvertices[3 * vertexIds[3] + 1] - outvertices[3 * vertexIds[1] + 1]; - double b23 = outvertices[3 * vertexIds[3] + 2] - outvertices[3 * vertexIds[1] + 2]; + double a11 = vertices[1][0] - vertices[0][0]; + double a12 = vertices[1][1] - vertices[0][1]; + double a13 = vertices[1][2] - vertices[0][2]; + double a21 = vertices[2][0] - vertices[0][0]; + double a22 = vertices[2][1] - vertices[0][1]; + double a23 = vertices[2][2] - vertices[0][2]; + double a31 = vertices[3][0] - vertices[0][0]; + double a32 = vertices[3][1] - vertices[0][1]; + double a33 = vertices[3][2] - vertices[0][2]; + double b11 = vertices[2][0] - vertices[1][0]; + double b12 = vertices[2][1] - vertices[1][1]; + double b13 = vertices[2][2] - vertices[1][2]; + double b21 = vertices[3][0] - vertices[1][0]; + double b22 = vertices[3][1] - vertices[1][1]; + double b23 = vertices[3][2] - vertices[1][2]; double det = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a13 * a22 * a31 - a12 * a21 * a33 - a11 * a23 * a32; double gram = std::abs(det); From be03f73895445ce06ed85c3f5706a8cf0f4f7826 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 31 Oct 2023 10:52:01 +0100 Subject: [PATCH 06/34] Add cstddef includes --- src/aux/InsphereCalculator.h | 1 + src/aux/MPIConvenience.h | 1 + src/input/SerialMeshFile.h | 1 + src/input/SimModSuite.h | 1 + src/pumgen.cpp | 1 + 5 files changed, 5 insertions(+) diff --git a/src/aux/InsphereCalculator.h b/src/aux/InsphereCalculator.h index a019fb6..2b0bbb6 100644 --- a/src/aux/InsphereCalculator.h +++ b/src/aux/InsphereCalculator.h @@ -3,6 +3,7 @@ #include #include +#include // assumes a contiguous distribution of all vertices over all processes diff --git a/src/aux/MPIConvenience.h b/src/aux/MPIConvenience.h index 128439f..f5aa221 100644 --- a/src/aux/MPIConvenience.h +++ b/src/aux/MPIConvenience.h @@ -2,6 +2,7 @@ #define PUMGEN_AUX_MPI_CONVENIENCE_H_ #include +#include void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, MPI_Datatype sendtype, void* recvbuf, const int* recvsize, diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 3332710..067fc03 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -24,6 +24,7 @@ #include #include #include +#include #include "MeshData.h" #include "utils/logger.h" diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index 1e834c6..7e0b945 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -60,6 +60,7 @@ #include #include #include +#include // forward declare pAManager SModel_attManager(pModel model); diff --git a/src/pumgen.cpp b/src/pumgen.cpp index b3b7fc3..7ad6767 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -16,6 +16,7 @@ #include #include +#include #include #include #include From a6cd259b2537cd275cb8518c2df8bd6df80f20b8 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Thu, 16 Nov 2023 10:38:16 +0100 Subject: [PATCH 07/34] Fix boundaries --- src/input/MeshData.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input/MeshData.h b/src/input/MeshData.h index 526e5ca..d800eb7 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -41,7 +41,7 @@ class FullStorageMeshData : public MeshData { int bndShift = 8; void setBoundary(std::size_t cell, int face, int value) { - boundaryData[cell] = static_cast(value) << (face * bndShift); + boundaryData[cell] |= static_cast(value) << (face * bndShift); } void setup(std::size_t cellCount, std::size_t vertexCount) { From c5d91943ff543aec92839e84a080d4be9d4a5613 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 26 Dec 2023 21:03:29 +0100 Subject: [PATCH 08/34] Remove APF serial mesh file --- src/input/SerialMeshFile.h | 5 -- src/input/SerialMeshFileApf.h | 151 ---------------------------------- 2 files changed, 156 deletions(-) delete mode 100644 src/input/SerialMeshFileApf.h diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 067fc03..662d130 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -17,12 +17,7 @@ #include #endif // PARALLEL -#include "ApfConvertWrapper.h" #include "aux/Distributor.h" -#include -#include -#include -#include #include #include diff --git a/src/input/SerialMeshFileApf.h b/src/input/SerialMeshFileApf.h deleted file mode 100644 index 10ed374..0000000 --- a/src/input/SerialMeshFileApf.h +++ /dev/null @@ -1,151 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#ifndef SERIAL_MESH_FILE_H -#define SERIAL_MESH_FILE_H - -#ifdef PARALLEL -#include -#endif // PARALLEL - -#include "ApfConvertWrapper.h" -#include "aux/Distributor.h" -#include -#include -#include -#include -#include - -#include "MeshInput.h" - -/** - * Read a mesh from a serial file - */ -template class SerialMeshFile : public MeshInput { - private: -#ifdef PARALLEL - MPI_Comm m_comm; -#endif // PARALLEL - - int m_rank; - int m_nProcs; - - T m_meshReader; - - public: -#ifdef PARALLEL - SerialMeshFile(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) - : m_comm(comm), m_meshReader(comm) { - init(); - open(meshFile); - } -#else // PARALLEL - SerialMeshFile(const char* meshFile) { - init(); - open(meshFile); - } -#endif // PARALLEL - - private: - /** - * Sets some parameters (called from the constructor) - */ - void init() { -#ifdef PARALLEL - MPI_Comm_rank(m_comm, &m_rank); - MPI_Comm_size(m_comm, &m_nProcs); -#else // PARALLLEL - m_rank = 0; - m_nProcs = 1; -#endif // PARALLEL - } - - void open(const char* meshFile) { - m_meshReader.open(meshFile); - - const std::size_t nVertices = m_meshReader.nVertices(); - const std::size_t nElements = m_meshReader.nElements(); - const std::size_t nLocalVertices = getChunksize(nVertices, m_rank, m_nProcs); - const std::size_t nLocalElements = getChunksize(nElements, m_rank, m_nProcs); - - gmi_register_null(); - gmi_model* model = gmi_load(".null"); - m_mesh = apf::makeEmptyMdsMesh(model, 3, false); - - // Create elements - logInfo(m_rank) << "Handle cell connectivity"; - apf::GlobalToVert vertMap; - { - std::vector elements; - { - std::vector elementsInt(nLocalElements * 4); - m_meshReader.readElements(elementsInt.data()); - elements = std::vector(elementsInt.begin(), elementsInt.end()); - } - logInfo(m_rank) << "Create APF connectivity"; - apf::construct(m_mesh, elements.data(), nLocalElements, apf::Mesh::TET, vertMap); - } - - apf::alignMdsRemotes(m_mesh); - apf::deriveMdsModel(m_mesh); - - // Set vertices - logInfo(m_rank) << "Handle vertex coordinates"; - { - std::vector vertices(nLocalVertices * 3); - m_meshReader.readVertices(vertices.data()); - apf::setCoords(m_mesh, vertices.data(), nLocalVertices, vertMap); - } - - // Set boundaries - logInfo(m_rank) << "Handle boundary conditions"; - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - { - std::vector boundaries(nLocalElements * 4); - m_meshReader.readBoundaries(boundaries.data()); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::Adjacent adjacent; - m_mesh->getAdjacent(element, 2, adjacent); - - for (int j = 0; j < 4; j++) { - if (!boundaries[i * 4 + j]) { - continue; - } - - m_mesh->setIntTag(adjacent[j], boundaryTag, &boundaries[i * 4 + j]); - } - - i++; - } - m_mesh->end(it); - } - - // Set groups - logInfo(m_rank) << "Handle cell groups"; - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - { - std::vector groups(nLocalElements); - m_meshReader.readGroups(groups.data()); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - m_mesh->setIntTag(element, groupTag, &groups[i]); - i++; - } - m_mesh->end(it); - } - } -}; - -#endif // SERIAL_MESH_FILE_H From ba84b61a7b58993ed83e68e697327be94de52d20 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 26 Dec 2023 21:13:22 +0100 Subject: [PATCH 09/34] Make clang-format happy --- src/meshreader/GMSHBuilder.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/meshreader/GMSHBuilder.h b/src/meshreader/GMSHBuilder.h index 4bc11be..f4c35c9 100644 --- a/src/meshreader/GMSHBuilder.h +++ b/src/meshreader/GMSHBuilder.h @@ -10,6 +10,9 @@ namespace puml { template struct GMSHSimplexType {}; + +// clang format seems to mess up the following lines from time to time, hence we disable it +// clang-format off template <> struct GMSHSimplexType<0u> { static constexpr long type = 15; }; // 15 -> point @@ -22,6 +25,7 @@ template <> struct GMSHSimplexType<2u> { template <> struct GMSHSimplexType<3u> { static constexpr long type = 4; }; // 4 -> tetrahedron +// clang-format on template class GMSHBuilder : public tndm::GMSHMeshBuilder { public: From a15910062c17f7105a0cdeda65be84f439f4ced0 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 01:29:06 +0100 Subject: [PATCH 10/34] Adjust for building without SCOREC --- src/input/SerialMeshFile.h | 2 +- src/input/SimModSuite.h | 2 +- src/input/SimModSuiteApf.h | 2 +- src/pumgen.cpp | 53 ++++++++++++++++++++++++-------------- 4 files changed, 37 insertions(+), 22 deletions(-) diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 662d130..17d548a 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -18,8 +18,8 @@ #endif // PARALLEL #include "aux/Distributor.h" -#include #include +#include #include "MeshData.h" #include "utils/logger.h" diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index 7e0b945..633d2d5 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -53,6 +53,7 @@ #include #include +#include #include #include #include @@ -60,7 +61,6 @@ #include #include #include -#include // forward declare pAManager SModel_attManager(pModel model); diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h index b28350a..41f63b5 100644 --- a/src/input/SimModSuiteApf.h +++ b/src/input/SimModSuiteApf.h @@ -72,7 +72,7 @@ pAManager SModel_attManager(pModel model); * of this class * @todo Maybe add MS_setMaxEntities to limit the number of elements */ -class SimModSuite : public MeshInput { +class SimModSuite : public ApfMeshInput { private: EasiMeshSize easiMeshSize; pGModel m_model; diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 7ad6767..9f860d9 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -10,8 +10,6 @@ * @author Sebastian Rettenberger */ -#include -#include #include #include @@ -26,10 +24,16 @@ #include +#ifdef USE_SCOREC #include #include // #include #include +#include "input/ApfNative.h" +#ifdef USE_SIMMOD +#include "input/SimModSuiteApf.h" +#endif +#endif #include #include "third_party/MPITraits.h" @@ -40,7 +44,6 @@ #ifdef USE_NETCDF #include "input/NetCDFMesh.h" #endif // USE_NETCDF -#include "input/ApfNative.h" #ifdef USE_SIMMOD #include "input/SimModSuite.h" #endif // USE_SIMMOD @@ -70,14 +73,6 @@ static int ilog(std::size_t value, int expbase = 1) { return count; } -template static void iterate(apf::Mesh* mesh, int dim, F&& function) { - apf::MeshIterator* it = mesh->begin(dim); - while (apf::MeshEntity* element = mesh->iterate(it)) { - std::invoke(function, element); - } - mesh->end(it); -} - constexpr std::size_t NoSecondDim = 0; template @@ -184,15 +179,15 @@ int main(int argc, char* argv[]) { MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &processes); +#ifdef USE_SCOREC PCU_Comm_Init(); +#endif // Parse command line arguments utils::Args args; - const char* source[] = {"gambit", "fidap", "msh2", "msh4", "netcdf", "apf", "simmodsuite"}; + const char* source[] = {"gambit", "fidap", "msh2", "msh4", + "netcdf", "apf", "simmodsuite", "simmodsuite-apf"}; args.addEnumOption("source", source, 's', "Mesh source (default: gambit)", false); - args.addOption("dump", 'd', "Dump APF mesh before partitioning it", utils::Args::Required, false); - args.addOption("model", 0, "Dump/Load a specific model file", utils::Args::Required, false); - args.addOption("vtk", 0, "Dump mesh to VTK files", utils::Args::Required, false); const char* filters[] = {"none", "scaleoffset", "deflate1", "deflate2", "deflate3", "deflate4", "deflate5", "deflate6", @@ -325,7 +320,12 @@ int main(int argc, char* argv[]) { break; case 5: logInfo(rank) << "Using APF native format"; - // meshInput = new ApfNative(inputFile, args.getArgument("model", 0L)); +#ifdef USE_SCOREC + meshInput = ApfNative(inputFile, args.getArgument("input", 0L)).generate(); +#else + logError() << "This version of PUMgen has been compiled without SCOREC. Hence, the APF format " + "is not available."; +#endif break; case 6: #ifdef USE_SIMMOD @@ -338,8 +338,23 @@ int main(int argc, char* argv[]) { args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), args.isSet("analyseAR"), args.getArgument("sim_log", 0L)); #else // USE_SIMMOD - logError() << "SimModSuite is not supported in this version"; + logError() << "SimModSuite is not supported in this version."; #endif // USE_SIMMOD + case 7: +#ifdef USE_SIMMOD +#ifdef USE_SCOREC + meshInput = new SimModSuiteApf( + inputFile, args.getArgument("cad", 0L), + args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), + args.getArgument("analysis", "analysis"), + args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), + args.isSet("analyseAR"), args.getArgument("sim_log", 0L)).generate(); +#else + logError() << "This version of PUMgen has been compiled without SCOREC. Hence, this reader for the SimModSuite is not available here."; +#endif +#else + logError() << "SimModSuite is not supported in this version."; +#endif break; default: logError() << "Unknown source"; @@ -349,8 +364,6 @@ int main(int argc, char* argv[]) { void* mesh = nullptr; - // TODO(David): replace by apf::countOwned again once extended - // Get local/global size std::size_t localSize[2] = {meshInput->cellCount(), meshInput->vertexCount()}; std::size_t globalSize[2] = {localSize[0], localSize[1]}; @@ -473,7 +486,9 @@ int main(int argc, char* argv[]) { logInfo(rank) << "Finished successfully"; +#ifdef USE_SCOREC PCU_Comm_Free(); +#endif MPI_Finalize(); return 0; From f06da888683394f1547d7fe6ea052ab4babc359a Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 01:30:01 +0100 Subject: [PATCH 11/34] Remove FIDAP support --- src/meshreader/FidapReader.cpp | 24 -- src/meshreader/FidapReader.h | 414 --------------------------- src/meshreader/ParallelFidapReader.h | 237 --------------- src/pumgen.cpp | 39 ++- 4 files changed, 19 insertions(+), 695 deletions(-) delete mode 100644 src/meshreader/FidapReader.cpp delete mode 100644 src/meshreader/FidapReader.h delete mode 100644 src/meshreader/ParallelFidapReader.h diff --git a/src/meshreader/FidapReader.cpp b/src/meshreader/FidapReader.cpp deleted file mode 100644 index fe39208..0000000 --- a/src/meshreader/FidapReader.cpp +++ /dev/null @@ -1,24 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#include "FidapReader.h" - -const char* FidapReader::FIDAP_FILE_ID = "** FIDAP NEUTRAL FILE"; -const char* FidapReader::NODAL_COORDINATES = "NODAL COORDINATES"; -const char* FidapReader::BOUNDARY_CONDITIONS = "BOUNDARY CONDITIONS"; -const char* FidapReader::ELEMENT_GROUPS = "ELEMENT GROUPS"; - -const char* FidapReader::ZONE_GROUP = "GROUP:"; -const char* FidapReader::ZONE_ELEMENTS = "ELEMENTS:"; -const char* FidapReader::ZONE_NODES = "NODES:"; -const char* FidapReader::ZONE_GEOMETRY = "GEOMETRY:"; -const char* FidapReader::ZONE_TYPE = "TYPE:"; diff --git a/src/meshreader/FidapReader.h b/src/meshreader/FidapReader.h deleted file mode 100644 index cbc39cc..0000000 --- a/src/meshreader/FidapReader.h +++ /dev/null @@ -1,414 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#ifndef FIDAP_READER_H -#define FIDAP_READER_H - -#include -#include -#include -#include -#include -#include - -#include "utils/stringutils.h" - -#include "MeshReader.h" - -struct FidapBoundaryFace { - /** The vertices of the face */ - std::array vertices; - /** The type of the boundary */ - int type; -}; - -/** Defines a face by three vertices */ -struct FaceVertex { - FaceVertex() {} - - FaceVertex(std::size_t v1, std::size_t v2, std::size_t v3) { - vertices[0] = v1; - vertices[1] = v2; - vertices[2] = v3; - std::sort(vertices.begin(), vertices.end()); - } - - FaceVertex(const std::array& v) { - std::copy(v.begin(), v.end(), vertices.begin()); - std::sort(vertices.begin(), vertices.end()); - } - - bool operator<(const FaceVertex& other) const { - return vertices[0] < other.vertices[0] || - (vertices[0] == other.vertices[0] && vertices[1] < other.vertices[1]) || - (vertices[0] == other.vertices[0] && vertices[1] == other.vertices[1] && - vertices[2] < other.vertices[2]); - } - - std::array vertices; -}; - -/** Defines a face by element and face side */ -struct FaceElement { - FaceElement() { - element = -1; - side = -1; - } - - FaceElement(std::size_t e, int s) { - element = e; - side = s; - } - - std::size_t element; - int side; -}; - -class FidapReader : public MeshReader { - private: - struct ElementSection : FileSection { - /** Group id for elements in this section */ - int group; - }; - - struct BoundarySection : FileSection { - /** Type of the boundary */ - int type; - }; - - // Section descriptions - FileSection m_vertices; - std::vector m_elements; - std::vector m_boundaries; - - public: - void open(const char* meshFile) { - MeshReader::open(meshFile); - - std::string line; - - // Read header information - getline(m_mesh, line); - utils::StringUtils::trim(line); - if (line != FIDAP_FILE_ID) - logError() << "Not an Fidap mesh file:" << meshFile; - - std::string name; - getline(m_mesh, name); // Internal name - - getline(m_mesh, line); // Skip line: - // VERSION: x.y - getline(m_mesh, line); // Date - getline(m_mesh, line); // Empty line - - std::size_t nElements, nZones, t, dimensions; - m_mesh >> m_vertices.nLines; - m_mesh >> nElements; - m_mesh >> nZones; - m_mesh >> t; // I don't know what this variable contains - m_mesh >> dimensions; - if (dimensions != 3) - logError() << "Fidap file does not contain a 3 dimensional mesh"; - getline(m_mesh, line); // Read newline - - getline(m_mesh, line); // Empty line - getline(m_mesh, line); // Unknown line - getline(m_mesh, line); // Empty line - getline(m_mesh, line); // Unknown line - getline(m_mesh, line); // Empty line - getline(m_mesh, line); // Unknown line - - // Find seek positions - // Vertices - getline(m_mesh, line); - if (line.find(NODAL_COORDINATES) == std::string::npos) - logError() << "Invalid Fidap format: Coordinates expected, found" << line; - m_vertices.seekPosition = m_mesh.tellg(); - getline(m_mesh, line); - m_vertices.lineSize = line.length() + 1; - m_mesh.seekg(m_vertices.seekPosition + m_vertices.nLines * m_vertices.lineSize); - - // Boundary conditions (currently ignored) - getline(m_mesh, line); - if (line.find(BOUNDARY_CONDITIONS) == std::string::npos) - logError() << "Invalid Fidap format: Boundary conditions expected, found" << line; - getline(m_mesh, line); - - // Element groups - getline(m_mesh, line); - if (line.find(ELEMENT_GROUPS) == std::string::npos) - logError() << "Invalid Fidap format: Element groups expected, found" << line; - - for (std::size_t i = 0; i < nZones; i++) { - FileSection sec; - - // Group (unused) - m_mesh >> name; - if (name.find(ZONE_GROUP) == std::string::npos) - logError() << "Invalid Fidap format: Expected group, found" << name; - unsigned int id; - m_mesh >> id; - // Elements - m_mesh >> name; - if (name.find(ZONE_ELEMENTS) == std::string::npos) - logError() << "Invalid Fidap format: Expected elements, found" << name; - m_mesh >> sec.nLines; - // Nodes - m_mesh >> name; - if (name.find(ZONE_NODES) == std::string::npos) - logError() << "Invalid Fidap format: Expected nodes, found" << name; - unsigned int nodeType; - m_mesh >> nodeType; - // Geometry (unused) - m_mesh >> name; - if (name.find(ZONE_GEOMETRY) == std::string::npos) - logError() << "Invalid Fidap format: Expected geometry, found" << name; - unsigned int geoType; - m_mesh >> geoType; - // Type (unused) - m_mesh >> name; - if (name.find(ZONE_TYPE) == std::string::npos) - logError() << "Invalid Fidap format: Expected type, found" << name; - unsigned int type; - m_mesh >> type; - getline(m_mesh, line); // Read newline - - // Entity name - m_mesh >> name; - m_mesh >> name; - std::string entityName; - getline(m_mesh, entityName); - utils::StringUtils::trim(entityName); - - sec.seekPosition = m_mesh.tellg(); - getline(m_mesh, line); - sec.lineSize = line.size() + 1; - - switch (nodeType) { - case 3: { - BoundarySection boundary; - boundary.seekPosition = sec.seekPosition; - boundary.nLines = sec.nLines; - boundary.lineSize = sec.lineSize; - - boundary.type = convertType(entityName); - - m_boundaries.push_back(boundary); - break; - } - case 4: { - ElementSection element; - element.seekPosition = sec.seekPosition; - element.nLines = sec.nLines; - element.lineSize = sec.lineSize; - - element.group = m_elements.size() + 1; - - m_elements.push_back(element); - break; - } - default: - logError() << "Unsupported node type" << nodeType << "found"; - } - - m_mesh.seekg(sec.seekPosition + sec.nLines * sec.lineSize); - } - } - - std::size_t nVertices() const { return m_vertices.nLines; } - - std::size_t nElements() const { - std::size_t count = 0; - for (const auto& element : m_elements) { - count += element.nLines; - } - - return count; - } - - std::size_t nBoundaries() const { - std::size_t count = 0; - for (const auto& boundary : m_boundaries) { - count += boundary.nLines; - } - - return count; - } - - void readVertices(std::size_t start, std::size_t count, double* vertices) { - m_mesh.clear(); - - m_mesh.seekg(m_vertices.seekPosition + start * m_vertices.lineSize + m_vertices.lineSize - - 3 * COORDINATE_SIZE - 1); - - for (std::size_t i = 0; i < count; i++) { - for (int j = 0; j < 3; j++) - m_mesh >> vertices[i * 3 + j]; - - m_mesh.seekg(m_vertices.lineSize - 3 * COORDINATE_SIZE, std::fstream::cur); - } - } - - void readElements(std::size_t start, std::size_t count, std::size_t* elements) { - m_mesh.clear(); - - // Get the boundary, were we should start reading - std::vector::const_iterator section; - for (section = m_elements.begin(); section != m_elements.end() && section->nLines < start; - section++) - start -= section->nLines; - - m_mesh.seekg(section->seekPosition + start * section->lineSize + section->lineSize - - 4 * VERTEX_ID_SIZE - 1); - - for (std::size_t i = 0; i < count; i++) { - if (start >= section->nLines) { - // Are we starting a new section in this iteration? - start = 0; - section++; - - m_mesh.seekg(section->seekPosition + section->lineSize - 4 * VERTEX_ID_SIZE - 1); - } - - for (int j = 0; j < 4; j++) { - m_mesh >> elements[i * 4 + j]; - elements[i * 4 + j]--; - } - - m_mesh.seekg(section->lineSize - 4 * VERTEX_ID_SIZE, std::fstream::cur); - - start++; // Line in the current section - } - } - - /** - * Read group ids from start to start+count. The group ids are sorted - * in the same order as the elements. - */ - void readGroups(std::size_t start, std::size_t count, int* groups) { - // Get the boundary, were we should start reading - std::vector::const_iterator section; - for (section = m_elements.begin(); section != m_elements.end() && section->nLines < start; - section++) - start -= section->nLines; - - for (std::size_t i = 0; i < count; i++) { - if (start >= section->nLines) { - // Are we starting a new section in this iteration? - start = 0; - section++; - } - - groups[i] = section->group; - - start++; // "Line" in the current section - } - } - - /** - * Read all groups - */ - void readGroups(int* groups) { readGroups(0, nElements(), groups); } - - void readBoundaries(std::size_t start, std::size_t count, FidapBoundaryFace* boundaries) { - m_mesh.clear(); - - // Get the boundary, were we should start reading - std::vector::const_iterator section; - for (section = m_boundaries.begin(); section != m_boundaries.end() && section->nLines < start; - section++) - start -= section->nLines; - - m_mesh.seekg(section->seekPosition + start * section->lineSize + section->lineSize - - 3 * VERTEX_ID_SIZE - 1); - - for (std::size_t i = 0; i < count; i++) { - if (start >= section->nLines) { - // Are we starting a new section in this iteration? - start = 0; - section++; - - m_mesh.seekg(section->seekPosition + section->lineSize - 3 * VERTEX_ID_SIZE - 1); - } - - for (int j = 0; j < 3; j++) { - m_mesh >> boundaries[i].vertices[j]; - boundaries[i].vertices[j]--; - } - boundaries[i].type = section->type; - - // Seek to next position - m_mesh.seekg(section->lineSize - 3 * VERTEX_ID_SIZE, std::fstream::cur); - - start++; // Line in the current section - } - } - - // TODO void readBoundaries(std::size_t* boundaries) {} - - public: - /** - * Creates a map from vertex definition of a face to a element/side definition - * - * @param elements Elements (4 vertices per element) - * @param n Number of elements - * @param[out] map The map - */ - static void createFaceMap(const std::size_t* elements, std::size_t n, - std::map& map) { - for (std::size_t i = 0; i < n; i++) { - FaceVertex v1(elements[i * 4], elements[i * 4 + 1], elements[i * 4 + 2]); - FaceElement e1(i, 0); - map[v1] = e1; - - FaceVertex v2(elements[i * 4], elements[i * 4 + 1], elements[i * 4 + 3]); - FaceElement e2(i, 1); - map[v2] = e2; - - FaceVertex v3(elements[i * 4 + 1], elements[i * 4 + 2], elements[i * 4 + 3]); - FaceElement e3(i, 2); - map[v3] = e3; - - FaceVertex v4(elements[i * 4], elements[i * 4 + 2], elements[i * 4 + 3]); - FaceElement e4(i, 3); - map[v4] = e4; - } - } - - private: - static int convertType(const std::string& name) { - std::istringstream ss(name.substr(2)); - - int type; - ss >> type; - - return type % 100; - } - - /** Number of characters required to store a coordinate */ - static const size_t COORDINATE_SIZE = 20ul; - /** Number of characters required to store a vertex id */ - static const size_t VERTEX_ID_SIZE = 8ul; - - static const char* FIDAP_FILE_ID; - static const char* NODAL_COORDINATES; - static const char* BOUNDARY_CONDITIONS; - static const char* ELEMENT_GROUPS; - - static const char* ZONE_GROUP; - static const char* ZONE_ELEMENTS; - static const char* ZONE_NODES; - static const char* ZONE_GEOMETRY; - static const char* ZONE_TYPE; -}; - -#endif // FIDAP_READER_H diff --git a/src/meshreader/ParallelFidapReader.h b/src/meshreader/ParallelFidapReader.h deleted file mode 100644 index dea856a..0000000 --- a/src/meshreader/ParallelFidapReader.h +++ /dev/null @@ -1,237 +0,0 @@ -/** - * @file - * This file is part of PUMGen - * - * For conditions of distribution and use, please see the copyright - * notice in the file 'COPYING' at the root directory of this package - * and the copyright notice at https://github.com/SeisSol/PUMGen - * - * @copyright 2017 Technical University of Munich - * @author Sebastian Rettenberger - */ - -#ifndef PARALLEL_FIDAP_READER -#define PARALLEL_FIDAP_READER - -#include - -#include - -#include "FidapReader.h" -#include "ParallelMeshReader.h" -#include "third_party/MPITraits.h" - -class ParallelFidapReader : public ParallelMeshReader { - private: - /** Face map required for boundaries */ - std::map m_faceMap; - - public: - ParallelFidapReader(MPI_Comm comm = MPI_COMM_WORLD) : ParallelMeshReader(comm) {} - - ParallelFidapReader(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) - : ParallelMeshReader(meshFile, comm) {} - - void readElements(std::size_t* elements) { - ParallelMeshReader::readElements(elements); - - // Create the face map - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; - const std::size_t maxChunkSize = chunkSize; - if (m_rank == m_nProcs - 1) { - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - } - - FidapReader::createFaceMap(elements, chunkSize, m_faceMap); - for (auto& face : m_faceMap) { - face.second.element += m_rank * maxChunkSize; - } - } - - /** - * @copydoc ParallelGambitReader::readGroups - */ - void readGroups(int* groups) { - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; - - if (m_rank == 0) { - // Allocate second buffer so we can read and send in parallel - std::vector tempGroups(chunkSize); - int* groups2 = tempGroups.data(); - if (m_nProcs % 2 == 0) - // swap once so we have the correct buffer at the end - swap(groups, groups2); - - MPI_Request request = MPI_REQUEST_NULL; - - for (int i = 1; i < m_nProcs - 1; i++) { - logInfo() << "Reading group information part" << i << "of" << m_nProcs; - m_serialReader.readGroups(i * chunkSize, chunkSize, groups); - MPI_Wait(&request, MPI_STATUS_IGNORE); - MPI_Isend(groups, chunkSize, MPI_INT, i, 0, m_comm, &request); - swap(groups, groups2); - } - - if (m_nProcs > 1) { - // Read last one - const std::size_t lastChunkSize = nElements() - (m_nProcs - 1) * chunkSize; - logInfo() << "Reading group information part" << (m_nProcs - 1) << "of" << m_nProcs; - m_serialReader.readGroups((m_nProcs - 1) * chunkSize, lastChunkSize, groups); - MPI_Wait(&request, MPI_STATUS_IGNORE); - MPI_Isend(groups, lastChunkSize, MPI_INT, m_nProcs - 1, 0, m_comm, &request); - swap(groups, groups2); - } - - // Finally read the first part - logInfo() << "Reading group information part" << m_nProcs << "of" << m_nProcs; - m_serialReader.readGroups(0, chunkSize, groups); - MPI_Wait(&request, MPI_STATUS_IGNORE); - } else { - if (m_rank == m_nProcs - 1) - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - - MPI_Recv(groups, chunkSize, MPI_INT, 0, 0, m_comm, MPI_STATUS_IGNORE); - } - } - - /** - * @copydoc ParallelGambitReader::readBoundaries - */ - void readBoundaries(int* boundaries) { - // Create the distributed face map - int vertexChunk = (nVertices() + m_nProcs - 1) / m_nProcs; - PCU_Comm_Begin(); - for (std::map::iterator i = m_faceMap.begin(); i != m_faceMap.end();) { - int proc = i->first.vertices[1] / vertexChunk; - if (proc == m_rank) - i++; - else { - PCU_COMM_PACK(proc, i->first); - PCU_COMM_PACK(proc, i->second); - m_faceMap.erase(i++); - } - } - PCU_Comm_Send(); - while (PCU_Comm_Receive()) { - FaceVertex v; - PCU_COMM_UNPACK(v); - FaceElement e; - PCU_COMM_UNPACK(e); - m_faceMap[v] = e; - } - - // Distribute the faces - std::vector facePos; - std::vector faceType; - - std::size_t chunkSize = (nElements() + m_nProcs - 1) / m_nProcs; - const std::size_t maxChunkSize = chunkSize; - - if (m_rank == 0) { - std::vector faces(maxChunkSize); - - std::vector> aggregator(m_nProcs - 1); - std::vector sizes(m_nProcs - 1); - std::vector requests((m_nProcs - 1) * 2, MPI_REQUEST_NULL); - - const std::size_t nChunks = (nBoundaries() + chunkSize - 1) / chunkSize; - for (std::size_t i = 0; i < nChunks; i++) { - logInfo() << "Reading boundary conditions part" << (i + 1) << "of" << nChunks; - - if (i == nChunks - 1) { - chunkSize = nBoundaries() - (nChunks - 1) * chunkSize; - } - - m_serialReader.readBoundaries(i * maxChunkSize, chunkSize, faces.data()); - - // Wait for all sending from last iteration - MPI_Waitall((m_nProcs - 1) * 2, requests.data(), MPI_STATUSES_IGNORE); - for (int j = 0; j < m_nProcs - 1; j++) { - aggregator[j].clear(); - } - - // Sort boundary conditions into the corresponding aggregator - for (std::size_t j = 0; j < chunkSize; j++) { - FaceVertex v(faces[j].vertices); - - if (v.vertices[1] < vertexChunk) { - facePos.push_back(m_faceMap.at(v)); - faceType.push_back(faces[j].type); - } else { - // Face for another processor - // Serials the struct to make sending easier - int proc = v.vertices[1] / vertexChunk; - aggregator[proc - 1].insert(aggregator[proc - 1].end(), v.vertices.begin(), - v.vertices.end()); - aggregator[proc - 1].push_back(faces[j].type); - } - } - - // Send send aggregated values - for (int j = 0; j < m_nProcs - 1; j++) { - if (aggregator[j].empty()) - continue; - - sizes[j] = aggregator[j].size() / 4; // 3 vertices + face type - MPI_Isend(&sizes[j], 1, tndm::mpi_type_t(), j + 1, 0, m_comm, - &requests[j * 2]); - MPI_Isend(&aggregator[j][0], aggregator[j].size(), MPI_INT, j + 1, 0, m_comm, - &requests[j * 2 + 1]); - } - } - - MPI_Waitall((m_nProcs - 1) * 2, requests.data(), MPI_STATUSES_IGNORE); - - // Send finish signal to all other processes - std::fill(sizes.begin(), sizes.end(), 0); - for (int i = 0; i < m_nProcs - 1; i++) { - MPI_Isend(&sizes[i], 1, tndm::mpi_type_t(), i + 1, 0, m_comm, &requests[i]); - } - MPI_Waitall(m_nProcs - 1, requests.data(), MPI_STATUSES_IGNORE); - } else { - if (m_rank == m_nProcs - 1) - chunkSize = nElements() - (m_nProcs - 1) * chunkSize; - - // Allocate enough memory - std::vector buf(chunkSize * 4); - - while (true) { - std::size_t size; - MPI_Recv(&size, 1, tndm::mpi_type_t(), 0, 0, m_comm, MPI_STATUS_IGNORE); - if (size == 0) - // Finished - break; - - MPI_Recv(buf.data(), size * 4, tndm::mpi_type_t(), 0, 0, m_comm, - MPI_STATUS_IGNORE); - - for (std::size_t i = 0; i < size * 4; i += 4) { - FaceVertex v(buf[i], buf[i + 1], buf[i + 2]); - facePos.push_back(m_faceMap.at(v)); - faceType.push_back(buf[i + 3]); - } - } - } - - // Distribute the faces to the final ranks - PCU_Comm_Begin(); - for (std::size_t i = 0; i < facePos.size(); i++) { - if (facePos[i].element / maxChunkSize == static_cast(m_rank)) - boundaries[(facePos[i].element % maxChunkSize) * 4 + facePos[i].side] = faceType[i]; - else { - PCU_COMM_PACK(facePos[i].element / maxChunkSize, facePos[i]); - PCU_COMM_PACK(facePos[i].element / maxChunkSize, faceType[i]); - } - } - PCU_Comm_Send(); - while (PCU_Comm_Receive()) { - FaceElement e; - PCU_COMM_UNPACK(e); - int type; - PCU_COMM_UNPACK(type); - boundaries[(e.element % maxChunkSize) * 4 + e.side] = type; - } - } -}; - -#endif // PARALLEL_FIDAP_READER diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 9f860d9..f53ef46 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -28,8 +28,8 @@ #include #include // #include -#include #include "input/ApfNative.h" +#include #ifdef USE_SIMMOD #include "input/SimModSuiteApf.h" #endif @@ -48,7 +48,6 @@ #include "input/SimModSuite.h" #endif // USE_SIMMOD #include "meshreader/GMSH4Parser.h" -#include "meshreader/ParallelFidapReader.h" #include "meshreader/ParallelGMSHReader.h" #include "meshreader/ParallelGambitReader.h" #include "third_party/GMSH2Parser.h" @@ -185,7 +184,7 @@ int main(int argc, char* argv[]) { // Parse command line arguments utils::Args args; - const char* source[] = {"gambit", "fidap", "msh2", "msh4", + const char* source[] = {"gambit", "msh2", "msh4", "netcdf", "apf", "simmodsuite", "simmodsuite-apf"}; args.addEnumOption("source", source, 's', "Mesh source (default: gambit)", false); @@ -299,18 +298,14 @@ int main(int argc, char* argv[]) { meshInput = new SerialMeshFile(inputFile); break; case 1: - logInfo(rank) << "Using Fidap mesh"; - meshInput = new SerialMeshFile(inputFile); - break; - case 2: logInfo(rank) << "Using GMSH mesh format 2 (msh2) mesh"; meshInput = new SerialMeshFile>(inputFile); break; - case 3: + case 2: logInfo(rank) << "Using GMSH mesh format 4 (msh4) mesh"; meshInput = new SerialMeshFile>(inputFile); break; - case 4: + case 3: #ifdef USE_NETCDF logInfo(rank) << "Using netCDF mesh"; meshInput = new NetCDFMesh(inputFile); @@ -318,7 +313,7 @@ int main(int argc, char* argv[]) { logError() << "netCDF is not supported in this version"; #endif // USE_NETCDF break; - case 5: + case 4: logInfo(rank) << "Using APF native format"; #ifdef USE_SCOREC meshInput = ApfNative(inputFile, args.getArgument("input", 0L)).generate(); @@ -327,7 +322,7 @@ int main(int argc, char* argv[]) { "is not available."; #endif break; - case 6: + case 5: #ifdef USE_SIMMOD logInfo(rank) << "Using SimModSuite"; @@ -340,24 +335,28 @@ int main(int argc, char* argv[]) { #else // USE_SIMMOD logError() << "SimModSuite is not supported in this version."; #endif // USE_SIMMOD - case 7: + case 6: #ifdef USE_SIMMOD #ifdef USE_SCOREC - meshInput = new SimModSuiteApf( - inputFile, args.getArgument("cad", 0L), - args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), - args.getArgument("analysis", "analysis"), - args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), - args.isSet("analyseAR"), args.getArgument("sim_log", 0L)).generate(); + meshInput = + new SimModSuiteApf(inputFile, args.getArgument("cad", 0L), + args.getArgument("license", 0L), + args.getArgument("mesh", "mesh"), + args.getArgument("analysis", "analysis"), + args.getArgument("enforce-size", 0), + args.getArgument("xml", 0L), args.isSet("analyseAR"), + args.getArgument("sim_log", 0L)) + .generate(); #else - logError() << "This version of PUMgen has been compiled without SCOREC. Hence, this reader for the SimModSuite is not available here."; + logError() << "This version of PUMgen has been compiled without SCOREC. Hence, this reader for " + "the SimModSuite is not available here."; #endif #else logError() << "SimModSuite is not supported in this version."; #endif break; default: - logError() << "Unknown source"; + logError() << "Unknown source."; } logInfo(rank) << "Parsed mesh successfully, writing output..."; From 12135d6b43c2402565d10a94b4b366120a2a410c Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 01:30:23 +0100 Subject: [PATCH 12/34] More SCOREC removal --- CMakeLists.txt | 10 +++++++--- cmake/FindSIMMETRIX.cmake | 9 +++++++-- src/input/ApfNative.h | 2 +- src/input/MeshData.h | 4 ++++ src/input/MeshInput.h | 11 ++++++++--- 5 files changed, 27 insertions(+), 9 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3baa622..848ca7a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,6 +13,7 @@ set(SIMMETRIX_ROOT "" CACHE STRING "Root directory of simmetrix dev files") set(SIM_MPI "mpich3" CACHE STRING "MPI version used by simmetrix") option(NETCDF "Use netcdf" OFF) option(PARASOLID "enable support for parasolid files" OFF) +option(SCOREC "Use PUMI/APF" OFF) set(LOG_LEVEL "info" CACHE STRING "Log level for the code") set_property(CACHE LOG_LEVEL PROPERTY STRINGS "debug" "info" "warning" "error") @@ -76,9 +77,12 @@ target_link_libraries(pumgen PUBLIC MPI::MPI_CXX) target_compile_definitions(pumgen PUBLIC PARALLEL) -find_package(APF REQUIRED) -target_include_directories(pumgen PUBLIC ${APF_INCLUDE_DIR}) -target_link_libraries(pumgen PUBLIC ${APF_LIBRARIES}) +if (SCOREC) + find_package(APF REQUIRED) + target_include_directories(pumgen PUBLIC ${APF_INCLUDE_DIR}) + target_link_libraries(pumgen PUBLIC ${APF_LIBRARIES}) + target_compile_definitions(pumgen PUBLIC USE_SCOREC) +endif() add_library(tinyxml2 submodules/tinyxml2/tinyxml2.cpp) target_link_libraries(pumgen PUBLIC tinyxml2) diff --git a/cmake/FindSIMMETRIX.cmake b/cmake/FindSIMMETRIX.cmake index 55601a0..8b77a97 100644 --- a/cmake/FindSIMMETRIX.cmake +++ b/cmake/FindSIMMETRIX.cmake @@ -30,8 +30,6 @@ find_library(SIM_PS_KRNL_LIB pskernel ${SIM_LIB_HINT}/psKrnl) get_filename_component(SIM_PS_KRNL_LIB_DIR ${SIM_PS_KRNL_LIB} DIRECTORY) list(APPEND SIMMETRIX_LIBRARIES - "${GMI_SIM_LIB}" - "${APF_SIM_LIB}" "${SIM_DISCRETE_LIB}" "${SIM_EXPORT_LIB}" "${SIM_MESHING_LIB}" @@ -42,6 +40,13 @@ list(APPEND SIMMETRIX_LIBRARIES "${SIM_PS_KRNL_LIB}" ) +if (SCOREC) +list(APPEND SIMMETRIX_LIBRARIES + "${GMI_SIM_LIB}" + "${APF_SIM_LIB}" +) +endif() + if (PARASOLID) list(APPEND SIMMETRIX_LIBRARIES "${SIM_PARASOLID_LIB}" diff --git a/src/input/ApfNative.h b/src/input/ApfNative.h index 46ab3c2..af59ccc 100644 --- a/src/input/ApfNative.h +++ b/src/input/ApfNative.h @@ -21,7 +21,7 @@ #include "MeshInput.h" -class ApfNative : public MeshInput { +class ApfNative : public ApfMeshInput { public: ApfNative(const char* mesh, const char* model = 0L) { if (model) diff --git a/src/input/MeshData.h b/src/input/MeshData.h index d800eb7..d14b699 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -5,6 +5,10 @@ #include #include +// TODO: to save space, the connectivity, geometry, group, and boundary arrays can also be +// constructed on the fly (for that, change the MeshData data structure to work on an iterator +// instead) + /** * Interface for mesh input */ diff --git a/src/input/MeshInput.h b/src/input/MeshInput.h index b30a3e2..c7a2e10 100644 --- a/src/input/MeshInput.h +++ b/src/input/MeshInput.h @@ -10,6 +10,7 @@ * @author Sebastian Rettenberger */ +#include "MeshData.h" #include #ifndef MESH_INTPUT_H @@ -18,14 +19,18 @@ /** * Interface for mesh input */ -class MeshInput { +class ApfMeshInput { protected: apf::Mesh2* m_mesh; public: - virtual ~MeshInput() {} - apf::Mesh2* getMesh() { return m_mesh; } + + MeshData* generate() { + FullStorageMeshData* data = new FullStorageMeshData(); + + return data; + } }; #endif // MESH_INPUT_H From fe74e88339fe1011a4321b56893ba58a687fcea5 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 01:30:31 +0100 Subject: [PATCH 13/34] Larger MPI --- src/aux/InsphereCalculator.cpp | 7 +++--- src/aux/InsphereCalculator.h | 2 +- src/aux/MPIConvenience.cpp | 37 ++++++++++++++++++++++++-------- src/aux/MPIConvenience.h | 6 +++++- src/input/ParallelVertexFilter.h | 7 +++--- 5 files changed, 42 insertions(+), 17 deletions(-) diff --git a/src/aux/InsphereCalculator.cpp b/src/aux/InsphereCalculator.cpp index 2845565..9f0735b 100644 --- a/src/aux/InsphereCalculator.cpp +++ b/src/aux/InsphereCalculator.cpp @@ -22,8 +22,8 @@ std::vector calculateInsphere(const std::vector& connectivi MPI_Type_contiguous(3, MPI_DOUBLE, &vertexType); MPI_Type_commit(&vertexType); - std::vector outrequests(commsize); - std::vector inrequests(commsize); + std::vector outrequests(commsize); + std::vector inrequests(commsize); std::size_t localVertices = geometry.size() / 3; @@ -51,7 +51,8 @@ std::vector calculateInsphere(const std::vector& connectivi } } - MPI_Alltoall(outrequests.data(), 1, MPI_INT, inrequests.data(), 1, MPI_INT, comm); + MPI_Alltoall(outrequests.data(), 1, tndm::mpi_type_t(), inrequests.data(), 1, + tndm::mpi_type_t(), comm); std::vector outdisp(commsize + 1); std::vector indisp(commsize + 1); diff --git a/src/aux/InsphereCalculator.h b/src/aux/InsphereCalculator.h index 2b0bbb6..acd4b04 100644 --- a/src/aux/InsphereCalculator.h +++ b/src/aux/InsphereCalculator.h @@ -1,9 +1,9 @@ #ifndef PUMGEN_AUX_INSPHERE_CALCULATOR_H_ #define PUMGEN_AUX_INSPHERE_CALCULATOR_H_ +#include #include #include -#include // assumes a contiguous distribution of all vertices over all processes diff --git a/src/aux/MPIConvenience.cpp b/src/aux/MPIConvenience.cpp index 2e8eef6..15d473f 100644 --- a/src/aux/MPIConvenience.cpp +++ b/src/aux/MPIConvenience.cpp @@ -1,9 +1,10 @@ #include "MPIConvenience.h" +#include #include #include -void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, - MPI_Datatype sendtype, void* recvbuf, const int* recvsize, +void sparseAlltoallv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const std::size_t* recvsize, const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { int commsize; int commrank; @@ -19,26 +20,44 @@ void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t std::size_t sendtypesize = sendtypesizePre; std::size_t recvtypesize = recvtypesizePre; - std::vector requests(commsize * 2, MPI_REQUEST_NULL); + std::vector requests; + requests.reserve(commsize * 2); // (no special handling of self-to-self comm at the moment) + constexpr std::size_t DataIncrement = std::numeric_limits::max(); + for (int i = 0; i < commsize; ++i) { - if (sendsize[i] > 0) { - MPI_Isend(reinterpret_cast(sendbuf) + senddisp[i] * sendtypesize, sendsize[i], - sendtype, i, Tag, comm, &requests[i]); + for (std::size_t position = 0; position < sendsize[i]; position += DataIncrement) { + int localSize = static_cast(std::min(sendsize[i] - position, DataIncrement)); + requests.push_back(MPI_REQUEST_NULL); + MPI_Isend(reinterpret_cast(sendbuf) + senddisp[i] * sendtypesize + position, + localSize, sendtype, i, Tag, comm, &requests.back()); } } for (int i = 0; i < commsize; ++i) { - if (recvsize[i] > 0) { - MPI_Irecv(reinterpret_cast(recvbuf) + recvdisp[i] * recvtypesize, recvsize[i], - recvtype, i, Tag, comm, &requests[i + commsize]); + for (std::size_t position = 0; position < recvsize[i]; position += DataIncrement) { + int localSize = static_cast(std::min(recvsize[i] - position, DataIncrement)); + requests.push_back(MPI_REQUEST_NULL); + MPI_Irecv(reinterpret_cast(recvbuf) + recvdisp[i] * recvtypesize + position, localSize, + recvtype, i, Tag, comm, &requests.back()); } } MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); } +void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const int* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { + int commsize; + MPI_Comm_size(comm, &commsize); + std::vector sendsizeSize(sendsize, sendsize + commsize); + std::vector recvsizeSize(recvsize, recvsize + commsize); + sparseAlltoallv(sendbuf, sendsizeSize.data(), senddisp, sendtype, recvbuf, recvsizeSize.data(), + recvdisp, recvtype, comm); +} + void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddisp, MPI_Datatype sendtype, void* recvbuf, const int* recvsize, const int* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { diff --git a/src/aux/MPIConvenience.h b/src/aux/MPIConvenience.h index f5aa221..a1a36bd 100644 --- a/src/aux/MPIConvenience.h +++ b/src/aux/MPIConvenience.h @@ -1,8 +1,12 @@ #ifndef PUMGEN_AUX_MPI_CONVENIENCE_H_ #define PUMGEN_AUX_MPI_CONVENIENCE_H_ -#include #include +#include + +void sparseAlltoallv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, const std::size_t* recvsize, + const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm); void sparseAlltoallv(const void* sendbuf, const int* sendsize, const std::size_t* senddisp, MPI_Datatype sendtype, void* recvbuf, const int* recvsize, diff --git a/src/input/ParallelVertexFilter.h b/src/input/ParallelVertexFilter.h index 360c5b8..144058e 100644 --- a/src/input/ParallelVertexFilter.h +++ b/src/input/ParallelVertexFilter.h @@ -164,15 +164,16 @@ class ParallelVertexFilter { } // Determine the (local and total) bucket size - std::vector bucketSize(m_numProcs); + std::vector bucketSize(m_numProcs); for (std::size_t i = 0; i < numVertices; i++) { ++bucketSize[bucket[i]]; } // Tell all processes what we are going to send them - std::vector recvSize(m_numProcs); + std::vector recvSize(m_numProcs); - MPI_Alltoall(bucketSize.data(), 1, MPI_INT, recvSize.data(), 1, MPI_INT, m_comm); + MPI_Alltoall(bucketSize.data(), 1, tndm::mpi_type_t(), recvSize.data(), 1, + tndm::mpi_type_t(), m_comm); std::size_t numSortVertices = 0; #ifdef _OPENMP From 9725ed8ee58c4fb3ac6a7180b042ad7cbac2248f Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 01:30:48 +0100 Subject: [PATCH 14/34] Clang-format --- src/pumgen.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pumgen.cpp b/src/pumgen.cpp index f53ef46..ee02fe7 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -184,8 +184,8 @@ int main(int argc, char* argv[]) { // Parse command line arguments utils::Args args; - const char* source[] = {"gambit", "msh2", "msh4", - "netcdf", "apf", "simmodsuite", "simmodsuite-apf"}; + const char* source[] = {"gambit", "msh2", "msh4", "netcdf", + "apf", "simmodsuite", "simmodsuite-apf"}; args.addEnumOption("source", source, 's', "Mesh source (default: gambit)", false); const char* filters[] = {"none", "scaleoffset", "deflate1", "deflate2", From f5c838fe339ae8e148a8e8f14770da5fcb804b1b Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 20:15:42 +0100 Subject: [PATCH 15/34] Remove FIDAP source file --- CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 848ca7a..f762fe9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -32,7 +32,6 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O2") add_executable(pumgen ${CMAKE_CURRENT_SOURCE_DIR}/src/pumgen.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/FidapReader.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GambitReader.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GMSH4Parser.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/input/ParallelVertexFilter.cpp From bd3398cf6b521d532faae3b9e27d7f1c3ab50c07 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 27 Dec 2023 23:21:58 +0100 Subject: [PATCH 16/34] Update NetCDF reader --- src/input/ApfConvertWrapper.h | 31 --------- src/input/NetCDFMesh.h | 116 ++++++++++------------------------ 2 files changed, 33 insertions(+), 114 deletions(-) delete mode 100644 src/input/ApfConvertWrapper.h diff --git a/src/input/ApfConvertWrapper.h b/src/input/ApfConvertWrapper.h deleted file mode 100644 index dc503c7..0000000 --- a/src/input/ApfConvertWrapper.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef APF_CONVERT_WRAPPER_H -#define APF_CONVERT_WRAPPER_H - -#include -#include - -// This is a small fix, because some time after PUMI (i.e. SCOREC/core) v2.2.7, -// the element ID was switched from a mere int to a named type apf::GiD, a typedef to long. -// Thus, PUMgen did not compile anymore, but to my knowledge, there was no flag to check the -// version. (of course, I could have missed some) To make PUMgen work with newer versions of PUMI, -// we therefore read out the type of the only function that actually was affected by the switch to -// apf::GiD and use that one in PUMgen. - -// inspired by https://stackoverflow.com/a/70954691 - -namespace internal { -template struct GetSecondFunctionArgument {}; - -template -struct GetSecondFunctionArgument { - using Type = Arg2; -}; - -using ElementIDConstPtr = GetSecondFunctionArgument::Type; -using ElementID = std::remove_const_t>; - -} // namespace internal - -using ElementID = internal::ElementID; - -#endif diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index 934b6fb..f241cd4 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -18,12 +18,7 @@ #include #include -#include "ApfConvertWrapper.h" -#include -#include -#include -#include - +#include "MeshData.h" #include "utils/logger.h" #include "MeshInput.h" @@ -33,7 +28,7 @@ /** * Read PUMGen generated mesh files */ -class NetCDFMesh : public MeshInput { +class NetCDFMesh : public FullStorageMeshData { public: NetCDFMesh(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) { int rank = 0; @@ -41,10 +36,6 @@ class NetCDFMesh : public MeshInput { MPI_Comm_rank(comm, &rank); MPI_Comm_size(comm, &nProcs); - gmi_register_null(); - gmi_model* model = gmi_load(".null"); - m_mesh = apf::makeEmptyMdsMesh(model, 3, false); - int ncFile; checkNcError(nc_open_par(meshFile, NC_MPIIO, comm, MPI_INFO_NULL, &ncFile)); @@ -71,14 +62,8 @@ class NetCDFMesh : public MeshInput { checkNcError(nc_open_par(meshFile, NC_MPIIO, commIO, MPI_INFO_NULL, &ncFile)); } - PCU_Switch_Comm(commIO); - std::size_t nElements = 0; std::size_t nVertices = 0; - std::vector elements; - std::vector vertices; - std::vector boundaries; - std::vector groups; if (nLocalPart > 0) { std::vector partitions(nLocalPart); @@ -161,94 +146,59 @@ class NetCDFMesh : public MeshInput { // Copy to the buffer std::vector elementsLocal(nElements * 4); - elements.resize(nElements * 4); - vertices.resize(nVertices * 3); + std::vector verticesLocal(nVertices * 3); - boundaries.resize(nElements * 4); - groups.resize(nElements); + std::size_t vertexOffset = 0; + for (std::size_t i = 0; i < nLocalPart; i++) { + std::copy_n(partitions[i].vertices(), partitions[i].nVertices() * 3, + geometryData.begin() + vertexOffset * 3); + vertexOffset += partitions[i].nVertices(); + } + + logInfo(rank) << "Running vertex filter"; + ParallelVertexFilter filter(commIO); + filter.filter(nVertices, verticesLocal); + + nVertices = filter.numLocalVertices(); + + setup(nElements, nVertices); + + std::copy_n(filter.localVertices().begin(), nVertices * 3, geometryData.begin()); std::size_t elementOffset = 0; - std::size_t vertexOffset = 0; for (std::size_t i = 0; i < nLocalPart; i++) { #ifdef _OPENMP -#pragma omp parallel schedule(static) +#pragma omp parallel for schedule(static) #endif for (std::size_t j = 0; j < partitions[i].nElements() * 4; j++) { elementsLocal[elementOffset * 4 + j] = partitions[i].elements()[j] + vertexOffset; } - std::copy_n(partitions[i].vertices(), partitions[i].nVertices() * 3, - vertices.begin() + vertexOffset * 3); - partitions[i].convertBoundary(); - std::copy_n(partitions[i].boundaries(), partitions[i].nElements() * 4, - boundaries.begin() + elementOffset * 4); + +#ifdef _OPENMP +#pragma omp parallel for schedule(static) +#endif + for (std::size_t j = 0; j < partitions[i].nElements(); ++j) { + for (int k = 0; k < 4; ++k) { + setBoundary(elementOffset + j, k, partitions[i].boundaries()[j * 4 + k]); + } + } + std::copy_n(partitions[i].groups(), partitions[i].nElements(), - groups.begin() + elementOffset); + groupData.begin() + elementOffset); elementOffset += partitions[i].nElements(); - vertexOffset += partitions[i].nVertices(); } - logInfo(rank) << "Running vertex filter"; - ParallelVertexFilter filter(commIO); - filter.filter(nVertices, vertices); - - // Create filtered vertex list - - nVertices = filter.numLocalVertices(); - vertices.resize(nVertices * 3); - std::copy_n(filter.localVertices(), nVertices * 3, vertices.begin()); - logInfo(rank) << "Converting local to global vertex identifier"; #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel for #endif for (std::size_t i = 0; i < nElements * 4; i++) { - elements[i] = filter.globalIds()[elementsLocal[i]]; - } - } - - logInfo(rank) << "Constructing the mesh"; - apf::GlobalToVert vertMap; - apf::construct(m_mesh, elements.data(), nElements, apf::Mesh::TET, vertMap); - - apf::alignMdsRemotes(m_mesh); - apf::deriveMdsModel(m_mesh); - - logInfo(rank) << "Set coordinates in APF"; - apf::setCoords(m_mesh, vertices.data(), nVertices, vertMap); - - // Set boundaries - apf::MeshTag* boundaryTag = m_mesh->createIntTag("boundary condition", 1); - apf::MeshIterator* it = m_mesh->begin(3); - std::size_t i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - apf::Adjacent adjacent; - m_mesh->getAdjacent(element, 2, adjacent); - - for (int j = 0; j < 4; j++) { - if (!boundaries[i * 4 + j]) - continue; - - m_mesh->setIntTag(adjacent[j], boundaryTag, &boundaries[i * 4 + j]); + connectivityData[i] = filter.globalIds()[elementsLocal[i]]; } - - i++; } - m_mesh->end(it); - - // Set groups - apf::MeshTag* groupTag = m_mesh->createIntTag("group", 1); - it = m_mesh->begin(3); - i = 0; - while (apf::MeshEntity* element = m_mesh->iterate(it)) { - m_mesh->setIntTag(element, groupTag, &groups[i]); - i++; - } - m_mesh->end(it); - - PCU_Switch_Comm(MPI_COMM_WORLD); } private: From 72678b4066a47e565e0da2992fd7acafefe3d202 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 13:17:39 +0100 Subject: [PATCH 17/34] Include fix --- src/input/NetCDFMesh.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index f241cd4..2a4a030 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -21,7 +21,7 @@ #include "MeshData.h" #include "utils/logger.h" -#include "MeshInput.h" +#include "MeshData.h" #include "NetCDFPartition.h" #include "ParallelVertexFilter.h" @@ -85,7 +85,7 @@ class NetCDFMesh : public FullStorageMeshData { bool useGroups = true; if (nc_inq_varid(ncFile, "element_group", &ncVarElemGroup) != NC_NOERR) { useGroups = false; - logWarning() << "No group found, using group 0 for all elements"; + logWarning(rank) << "No group found, using group 0 for all elements"; } else { collectiveAccess(ncFile, ncVarElemGroup); } From 17e66c78b8c4f0b813a9aa0b941aebf3ad367330 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 13:18:10 +0100 Subject: [PATCH 18/34] Add APF to new data converter --- src/input/MeshInput.h | 94 +++++++++++++++++++++++++++++++++++++++--- src/input/NetCDFMesh.h | 2 +- src/pumgen.cpp | 19 ++++----- 3 files changed, 99 insertions(+), 16 deletions(-) diff --git a/src/input/MeshInput.h b/src/input/MeshInput.h index c7a2e10..55c627c 100644 --- a/src/input/MeshInput.h +++ b/src/input/MeshInput.h @@ -11,7 +11,10 @@ */ #include "MeshData.h" +#include "utils/logger.h" #include +#include +#include #ifndef MESH_INTPUT_H #define MESH_INTPUT_H @@ -19,17 +22,98 @@ /** * Interface for mesh input */ -class ApfMeshInput { +class ApfMeshInput : public FullStorageMeshData { + private: + template static void iterate(apf::Mesh* mesh, int dim, F&& function) { + apf::MeshIterator* it = mesh->begin(dim); + std::size_t index = 0; + while (apf::MeshEntity* element = mesh->iterate(it)) { + std::invoke(function, index, element); + ++index; + } + mesh->end(it); + } + protected: - apf::Mesh2* m_mesh; + apf::Mesh2* m_mesh = nullptr; public: apf::Mesh2* getMesh() { return m_mesh; } - MeshData* generate() { - FullStorageMeshData* data = new FullStorageMeshData(); + ~ApfMeshInput() { + delete m_mesh; + m_mesh = nullptr; + } + + void generate() { + if (alignMdsMatches(m_mesh)) { + logWarning() << "Fixed misaligned matches"; + } + m_mesh->verify(); + apf::GlobalNumbering* vertexNum = apf::makeGlobal(apf::numberOwnedNodes(m_mesh, "vertices")); + apf::synchronize(vertexNum); + + apf::Sharing* sharing = apf::getSharing(m_mesh); + + setup(apf::countOwned(m_mesh, 3), apf::countOwned(m_mesh, 0)); + + // connectivity + iterate(m_mesh, 3, [&](auto index, auto* element) { + apf::NewArray vn; + apf::getElementNumbers(vertexNum, element, vn); + + for (int i = 0; i < 4; i++) { + connectivityData[index * 4 + i] = vn[i]; + } + }); + + std::size_t vertexIndex = 0; + + // geometry + iterate(m_mesh, 0, [&](auto index, auto* element) { + if (!sharing->isOwned(element)) { + return; + } + + apf::Vector3 point; + m_mesh->getPoint(element, 0, point); + double geometry[3]; + point.toArray(geometry); + + for (int i = 0; i < 3; ++i) { + geometryData[vertexIndex * 3 + i] = geometry[i]; + } + ++vertexIndex; + }); + + // groups + iterate(m_mesh, 3, [&](auto index, auto* element) { + int group; + mesh->getIntTag(element, groupTag, &group); + groupData[index] = group; + }); + + // boundary + iterate(m_mesh, 3, [&](auto index, auto* element) { + apf::Downward faces; + mesh->getDownward(element, 2, faces); + + for (int i = 0; i < 4; i++) { + if (mesh->hasTag(faces[i], boundaryTag)) { + int b; + mesh->getIntTag(faces[i], boundaryTag, &b); + + if (b <= 0 || (b > i32limit && boundaryType == 0) || + (b > i64limit && boundaryType == 1)) { + logError() << "Cannot handle boundary condition" << b; + } + + setBoundary(inex, i, b); + } + } + }); - return data; + delete sharing; } }; diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index 2a4a030..f742c95 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -193,7 +193,7 @@ class NetCDFMesh : public FullStorageMeshData { logInfo(rank) << "Converting local to global vertex identifier"; #ifdef _OPENMP -#pragma omp parallel for +#pragma omp parallel for #endif for (std::size_t i = 0; i < nElements * 4; i++) { connectivityData[i] = filter.globalIds()[elementsLocal[i]]; diff --git a/src/pumgen.cpp b/src/pumgen.cpp index ee02fe7..69f6929 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -316,7 +316,8 @@ int main(int argc, char* argv[]) { case 4: logInfo(rank) << "Using APF native format"; #ifdef USE_SCOREC - meshInput = ApfNative(inputFile, args.getArgument("input", 0L)).generate(); + meshInput = ApfNative(inputFile, args.getArgument("input", 0L)); + (dynamic_cast(meshInput))->generate(); #else logError() << "This version of PUMgen has been compiled without SCOREC. Hence, the APF format " "is not available."; @@ -338,15 +339,13 @@ int main(int argc, char* argv[]) { case 6: #ifdef USE_SIMMOD #ifdef USE_SCOREC - meshInput = - new SimModSuiteApf(inputFile, args.getArgument("cad", 0L), - args.getArgument("license", 0L), - args.getArgument("mesh", "mesh"), - args.getArgument("analysis", "analysis"), - args.getArgument("enforce-size", 0), - args.getArgument("xml", 0L), args.isSet("analyseAR"), - args.getArgument("sim_log", 0L)) - .generate(); + meshInput = new SimModSuiteApf( + inputFile, args.getArgument("cad", 0L), + args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), + args.getArgument("analysis", "analysis"), + args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), + args.isSet("analyseAR"), args.getArgument("sim_log", 0L)); + (dynamic_cast(meshInput))->generate(); #else logError() << "This version of PUMgen has been compiled without SCOREC. Hence, this reader for " "the SimModSuite is not available here."; From c7367e0fdff88ea0b72d35f43c638e3f31840ecf Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 13:41:33 +0100 Subject: [PATCH 19/34] Allow Simmod with and without APF --- src/input/MeshInput.h | 17 +++++++---------- src/input/SimModSuiteApf.h | 18 +++++++++--------- src/pumgen.cpp | 3 ++- 3 files changed, 18 insertions(+), 20 deletions(-) diff --git a/src/input/MeshInput.h b/src/input/MeshInput.h index 55c627c..95a5f0f 100644 --- a/src/input/MeshInput.h +++ b/src/input/MeshInput.h @@ -87,28 +87,25 @@ class ApfMeshInput : public FullStorageMeshData { }); // groups + apf::MeshTag* groupTag = m_mesh->findTag("group"); iterate(m_mesh, 3, [&](auto index, auto* element) { int group; - mesh->getIntTag(element, groupTag, &group); + m_mesh->getIntTag(element, groupTag, &group); groupData[index] = group; }); // boundary + apf::MeshTag* boundaryTag = m_mesh->findTag("boundary condition"); iterate(m_mesh, 3, [&](auto index, auto* element) { apf::Downward faces; - mesh->getDownward(element, 2, faces); + m_mesh->getDownward(element, 2, faces); for (int i = 0; i < 4; i++) { - if (mesh->hasTag(faces[i], boundaryTag)) { + if (m_mesh->hasTag(faces[i], boundaryTag)) { int b; - mesh->getIntTag(faces[i], boundaryTag, &b); + m_mesh->getIntTag(faces[i], boundaryTag, &b); - if (b <= 0 || (b > i32limit && boundaryType == 0) || - (b > i64limit && boundaryType == 1)) { - logError() << "Cannot handle boundary condition" << b; - } - - setBoundary(inex, i, b); + setBoundary(index, i, b); } } }); diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h index 41f63b5..c49cbf9 100644 --- a/src/input/SimModSuiteApf.h +++ b/src/input/SimModSuiteApf.h @@ -10,8 +10,8 @@ * @author Sebastian Rettenberger */ -#ifndef SIM_MOD_SUITE_H -#define SIM_MOD_SUITE_H +#ifndef SIM_MOD_SUITE_APF_H +#define SIM_MOD_SUITE_APF_H #include @@ -72,7 +72,7 @@ pAManager SModel_attManager(pModel model); * of this class * @todo Maybe add MS_setMaxEntities to limit the number of elements */ -class SimModSuite : public ApfMeshInput { +class SimModSuiteApf : public ApfMeshInput { private: EasiMeshSize easiMeshSize; pGModel m_model; @@ -83,10 +83,10 @@ class SimModSuite : public ApfMeshInput { bool m_log; public: - SimModSuite(const char* modFile, const char* cadFile = 0L, const char* licenseFile = 0L, - const char* meshCaseName = "mesh", const char* analysisCaseName = "analysis", - int enforceSize = 0, const char* xmlFile = 0L, const bool analyseAR = false, - const char* logFile = 0L) { + SimModSuiteApf(const char* modFile, const char* cadFile = 0L, const char* licenseFile = 0L, + const char* meshCaseName = "mesh", const char* analysisCaseName = "analysis", + int enforceSize = 0, const char* xmlFile = 0L, const bool analyseAR = false, + const char* logFile = 0L) { // Init SimModSuite SimModel_start(); SimPartitionedMesh_start(0L, 0L); @@ -227,7 +227,7 @@ class SimModSuite : public ApfMeshInput { MS_deleteMeshCase(analysisCase); } - virtual ~SimModSuite() { + virtual ~SimModSuiteApf() { M_release(m_simMesh); // We cannot delete the model here because it is still // connected to the mesh @@ -685,4 +685,4 @@ class SimModSuite : public ApfMeshInput { } }; -#endif // SIM_MOD_SUITE_H +#endif // SIM_MOD_SUITE_APF_H diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 69f6929..264fae0 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -25,6 +25,7 @@ #include #ifdef USE_SCOREC +#include #include #include // #include @@ -316,7 +317,7 @@ int main(int argc, char* argv[]) { case 4: logInfo(rank) << "Using APF native format"; #ifdef USE_SCOREC - meshInput = ApfNative(inputFile, args.getArgument("input", 0L)); + meshInput = new ApfNative(inputFile, args.getArgument("input", 0L)); (dynamic_cast(meshInput))->generate(); #else logError() << "This version of PUMgen has been compiled without SCOREC. Hence, the APF format " From 4fd6b8270b16b59b3816f958b2e914518c755919 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 13:41:40 +0100 Subject: [PATCH 20/34] Add mesh data check --- src/input/MeshData.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/input/MeshData.h b/src/input/MeshData.h index d14b699..a47a513 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -3,8 +3,11 @@ #define MESH_DATA_H #include +#include #include +#include "utils/logger.h" + // TODO: to save space, the connectivity, geometry, group, and boundary arrays can also be // constructed on the fly (for that, change the MeshData data structure to work on an iterator // instead) @@ -14,6 +17,8 @@ */ class MeshData { public: + virtual ~MeshData() = default; + virtual std::size_t cellCount() = 0; virtual std::size_t vertexCount() = 0; @@ -45,6 +50,12 @@ class FullStorageMeshData : public MeshData { int bndShift = 8; void setBoundary(std::size_t cell, int face, int value) { + constexpr auto i32limit = std::numeric_limits::max(); + constexpr auto i64limit = std::numeric_limits::max(); + if (value < 0 || (value > i32limit && bndShift == 8) || (value > i64limit && bndShift == 16)) { + logError() << "Cannot handle boundary condition" << value; + } + boundaryData[cell] |= static_cast(value) << (face * bndShift); } From 4ed6a1a7f0583afb7de61b931ccceac393d8ffd3 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 13:42:37 +0100 Subject: [PATCH 21/34] MPI_Init_thread --- src/pumgen.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 264fae0..7849b68 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -173,8 +173,9 @@ static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, int main(int argc, char* argv[]) { int rank = 0; int processes = 1; + int mpithreadstate = MPI_THREAD_SINGLE; - MPI_Init(&argc, &argv); + MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &mpithreadstate); MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &processes); From 56138043f38ecdf983b654143c45e3357a379a36 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 14:13:36 +0100 Subject: [PATCH 22/34] Fixing call simmodeler twice --- src/input/SimModSuite.h | 12 ++++++------ src/pumgen.cpp | 3 +++ 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index 633d2d5..d46ff76 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -189,14 +189,14 @@ class SimModSuite : public FullStorageMeshData { { size_t totalCount = 0; size_t ownedCount = 0; - VIter reg_it = M_vertexIter(PM_mesh(m_simMesh, i)); - while (pVertex vertex = VIter_next(reg_it)) { + VIter vert_it = M_vertexIter(PM_mesh(m_simMesh, i)); + while (pVertex vertex = VIter_next(vert_it)) { ++totalCount; if (EN_isOwnerProc((pEntity)vertex)) { ++ownedCount; } } - VIter_delete(reg_it); + VIter_delete(vert_it); vertexCount += ownedCount; vertexIDCount += totalCount; @@ -238,7 +238,6 @@ class SimModSuite : public FullStorageMeshData { std::size_t indexID = vertexIDStart[i]; VIter reg_it = M_vertexIter(PM_mesh(m_simMesh, i)); while (pVertex vertex = VIter_next(reg_it)) { - // TODO: EN_isOwnerProc() EN_setID((pEntity)vertex, indexID); double xyz[3]; V_coord(vertex, xyz); @@ -254,7 +253,7 @@ class SimModSuite : public FullStorageMeshData { } vertexCount = vertexFilter.numLocalVertices(); - logInfo(PMU_rank()) << "Local vertices (really, now):" << vertexCount; + logInfo(PMU_rank()) << "Local vertices (after filtering):" << vertexCount; setup(cellCount, vertexCount); std::copy(vertexFilter.localVertices().begin(), vertexFilter.localVertices().end(), @@ -268,11 +267,12 @@ class SimModSuite : public FullStorageMeshData { pRegion reg; while ((reg = RIter_next(reg_it))) { if (EN_isOwnerProc((pEntity)reg)) { + EN_setID((pEntity)reg, index); + pPList vertices = R_vertices(reg, 1); void* iter = nullptr; int j = 0; while (pVertex vertex = (pVertex)PList_next(vertices, &iter)) { - EN_setID((pEntity)reg, index); connectivityData[index * 4 + j] = vertexFilter.globalIds()[EN_id((pEntity)vertex)]; ++j; } diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 7849b68..d50ed7f 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -338,9 +338,12 @@ int main(int argc, char* argv[]) { #else // USE_SIMMOD logError() << "SimModSuite is not supported in this version."; #endif // USE_SIMMOD + break; case 6: #ifdef USE_SIMMOD #ifdef USE_SCOREC + logInfo(rank) << "Using SimModSuite (with APF)"; + meshInput = new SimModSuiteApf( inputFile, args.getArgument("cad", 0L), args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), From 66bd465ca736bfe4ad928bbdcd1e8a47aabacca9 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 14:27:58 +0100 Subject: [PATCH 23/34] Fix virtual destructors --- src/input/MeshData.h | 2 ++ src/input/MeshInput.h | 2 +- src/input/NetCDFMesh.h | 2 ++ src/input/SerialMeshFile.h | 3 +++ src/pumgen.cpp | 2 +- 5 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/input/MeshData.h b/src/input/MeshData.h index a47a513..3f56f42 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -30,6 +30,8 @@ class MeshData { class FullStorageMeshData : public MeshData { public: + virtual ~FullStorageMeshData() = default; + virtual std::size_t cellCount() { return cellCountValue; } virtual std::size_t vertexCount() { return vertexCountValue; } diff --git a/src/input/MeshInput.h b/src/input/MeshInput.h index 95a5f0f..c50cd02 100644 --- a/src/input/MeshInput.h +++ b/src/input/MeshInput.h @@ -40,7 +40,7 @@ class ApfMeshInput : public FullStorageMeshData { public: apf::Mesh2* getMesh() { return m_mesh; } - ~ApfMeshInput() { + virtual ~ApfMeshInput() { delete m_mesh; m_mesh = nullptr; } diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index f742c95..47d1d36 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -30,6 +30,8 @@ */ class NetCDFMesh : public FullStorageMeshData { public: + virtual ~NetCDFMesh() = default; + NetCDFMesh(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) { int rank = 0; int nProcs = 1; diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 17d548a..962548a 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -28,6 +28,9 @@ * Read a mesh from a serial file */ template class SerialMeshFile : public FullStorageMeshData { + public: + virtual ~SerialMeshFile() = default; + private: #ifdef PARALLEL MPI_Comm m_comm; diff --git a/src/pumgen.cpp b/src/pumgen.cpp index d50ed7f..2e207de 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -342,7 +342,7 @@ int main(int argc, char* argv[]) { case 6: #ifdef USE_SIMMOD #ifdef USE_SCOREC - logInfo(rank) << "Using SimModSuite (with APF)"; + logInfo(rank) << "Using SimModSuite with APF (deprecated)"; meshInput = new SimModSuiteApf( inputFile, args.getArgument("cad", 0L), From 59b404b8f4d68a0e49665ea3b497689a6cfa1d9b Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 14:50:55 +0100 Subject: [PATCH 24/34] Propagate boundary offset --- src/input/ApfNative.h | 3 ++- src/input/MeshData.h | 2 ++ src/input/MeshInput.h | 2 ++ src/input/NetCDFMesh.h | 3 ++- src/input/SerialMeshFile.h | 6 +++--- src/input/SimModSuite.h | 9 +++++---- src/input/SimModSuiteApf.h | 9 +++++---- src/pumgen.cpp | 16 +++++++++------- 8 files changed, 30 insertions(+), 20 deletions(-) diff --git a/src/input/ApfNative.h b/src/input/ApfNative.h index af59ccc..a9794bf 100644 --- a/src/input/ApfNative.h +++ b/src/input/ApfNative.h @@ -23,7 +23,8 @@ class ApfNative : public ApfMeshInput { public: - ApfNative(const char* mesh, const char* model = 0L) { + ApfNative(const char* mesh, int boundarySize, const char* model = 0L) + : ApfMeshInput(boundarySize) { if (model) gmi_register_mesh(); else { diff --git a/src/input/MeshData.h b/src/input/MeshData.h index 3f56f42..201cc54 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -41,6 +41,8 @@ class FullStorageMeshData : public MeshData { virtual const std::vector& boundary() { return boundaryData; } protected: + FullStorageMeshData(int boundarySize) : bndShift(boundarySize) {} + std::size_t cellCountValue = 0; std::size_t vertexCountValue = 0; diff --git a/src/input/MeshInput.h b/src/input/MeshInput.h index c50cd02..fae39a2 100644 --- a/src/input/MeshInput.h +++ b/src/input/MeshInput.h @@ -37,6 +37,8 @@ class ApfMeshInput : public FullStorageMeshData { protected: apf::Mesh2* m_mesh = nullptr; + ApfMeshInput(int boundarySize) : FullStorageMeshData(boundarySize) {} + public: apf::Mesh2* getMesh() { return m_mesh; } diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index 47d1d36..358b5ed 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -32,7 +32,8 @@ class NetCDFMesh : public FullStorageMeshData { public: virtual ~NetCDFMesh() = default; - NetCDFMesh(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) { + NetCDFMesh(const char* meshFile, int boundarySize, MPI_Comm comm = MPI_COMM_WORLD) + : FullStorageMeshData(boundarySize) { int rank = 0; int nProcs = 1; MPI_Comm_rank(comm, &rank); diff --git a/src/input/SerialMeshFile.h b/src/input/SerialMeshFile.h index 962548a..362b26e 100644 --- a/src/input/SerialMeshFile.h +++ b/src/input/SerialMeshFile.h @@ -43,13 +43,13 @@ template class SerialMeshFile : public FullStorageMeshData { public: #ifdef PARALLEL - SerialMeshFile(const char* meshFile, MPI_Comm comm = MPI_COMM_WORLD) - : m_comm(comm), m_meshReader(comm) { + SerialMeshFile(const char* meshFile, int boundarySize, MPI_Comm comm = MPI_COMM_WORLD) + : FullStorageMeshData(boundarySize), m_comm(comm), m_meshReader(comm) { init(); open(meshFile); } #else // PARALLEL - SerialMeshFile(const char* meshFile) { + SerialMeshFile(const char* meshFile, int boundarySize) : FullStorageMeshData(boundarySize) { init(); open(meshFile); } diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index d46ff76..416448a 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -81,10 +81,11 @@ class SimModSuite : public FullStorageMeshData { bool m_log; public: - SimModSuite(const char* modFile, const char* cadFile = 0L, const char* licenseFile = 0L, - const char* meshCaseName = "mesh", const char* analysisCaseName = "analysis", - int enforceSize = 0, const char* xmlFile = 0L, const bool analyseAR = false, - const char* logFile = 0L) { + SimModSuite(const char* modFile, int boundarySize, const char* cadFile = 0L, + const char* licenseFile = 0L, const char* meshCaseName = "mesh", + const char* analysisCaseName = "analysis", int enforceSize = 0, + const char* xmlFile = 0L, const bool analyseAR = false, const char* logFile = 0L) + : FullStorageMeshData(boundarySize) { // Init SimModSuite SimModel_start(); SimPartitionedMesh_start(0L, 0L); diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h index c49cbf9..ad07f30 100644 --- a/src/input/SimModSuiteApf.h +++ b/src/input/SimModSuiteApf.h @@ -83,10 +83,11 @@ class SimModSuiteApf : public ApfMeshInput { bool m_log; public: - SimModSuiteApf(const char* modFile, const char* cadFile = 0L, const char* licenseFile = 0L, - const char* meshCaseName = "mesh", const char* analysisCaseName = "analysis", - int enforceSize = 0, const char* xmlFile = 0L, const bool analyseAR = false, - const char* logFile = 0L) { + SimModSuiteApf(const char* modFile, int boundarySize, const char* cadFile = 0L, + const char* licenseFile = 0L, const char* meshCaseName = "mesh", + const char* analysisCaseName = "analysis", int enforceSize = 0, + const char* xmlFile = 0L, const bool analyseAR = false, const char* logFile = 0L) + : ApfMeshInput(boundarySize) { // Init SimModSuite SimModel_start(); SimPartitionedMesh_start(0L, 0L); diff --git a/src/pumgen.cpp b/src/pumgen.cpp index 2e207de..9f72f4d 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -297,20 +297,22 @@ int main(int argc, char* argv[]) { switch (args.getArgument("source", 0)) { case 0: logInfo(rank) << "Using Gambit mesh"; - meshInput = new SerialMeshFile(inputFile); + meshInput = new SerialMeshFile(inputFile, faceOffset); break; case 1: logInfo(rank) << "Using GMSH mesh format 2 (msh2) mesh"; - meshInput = new SerialMeshFile>(inputFile); + meshInput = + new SerialMeshFile>(inputFile, faceOffset); break; case 2: logInfo(rank) << "Using GMSH mesh format 4 (msh4) mesh"; - meshInput = new SerialMeshFile>(inputFile); + meshInput = + new SerialMeshFile>(inputFile, faceOffset); break; case 3: #ifdef USE_NETCDF logInfo(rank) << "Using netCDF mesh"; - meshInput = new NetCDFMesh(inputFile); + meshInput = new NetCDFMesh(inputFile, faceOffset); #else // USE_NETCDF logError() << "netCDF is not supported in this version"; #endif // USE_NETCDF @@ -318,7 +320,7 @@ int main(int argc, char* argv[]) { case 4: logInfo(rank) << "Using APF native format"; #ifdef USE_SCOREC - meshInput = new ApfNative(inputFile, args.getArgument("input", 0L)); + meshInput = new ApfNative(inputFile, faceOffset, args.getArgument("input", 0L)); (dynamic_cast(meshInput))->generate(); #else logError() << "This version of PUMgen has been compiled without SCOREC. Hence, the APF format " @@ -330,7 +332,7 @@ int main(int argc, char* argv[]) { logInfo(rank) << "Using SimModSuite"; meshInput = new SimModSuite( - inputFile, args.getArgument("cad", 0L), + inputFile, faceOffset, args.getArgument("cad", 0L), args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), args.getArgument("analysis", "analysis"), args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), @@ -345,7 +347,7 @@ int main(int argc, char* argv[]) { logInfo(rank) << "Using SimModSuite with APF (deprecated)"; meshInput = new SimModSuiteApf( - inputFile, args.getArgument("cad", 0L), + inputFile, faceOffset, args.getArgument("cad", 0L), args.getArgument("license", 0L), args.getArgument("mesh", "mesh"), args.getArgument("analysis", "analysis"), args.getArgument("enforce-size", 0), args.getArgument("xml", 0L), From 73d3149bb3882ccc4459296d8333ac72f22c8646 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Fri, 29 Dec 2023 15:06:01 +0100 Subject: [PATCH 25/34] Avoid looking for unused libraries --- cmake/FindSIMMETRIX.cmake | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/cmake/FindSIMMETRIX.cmake b/cmake/FindSIMMETRIX.cmake index 8b77a97..e8fdbef 100644 --- a/cmake/FindSIMMETRIX.cmake +++ b/cmake/FindSIMMETRIX.cmake @@ -1,18 +1,24 @@ include(FindPackageHandleStandardArgs) -find_path(GMI_SIM_INCLUDE_DIR gmi_sim.h) -find_path(APF_SIM_INCLUDE_DIR apfSIM.h) +if (SCOREC) + find_path(GMI_SIM_INCLUDE_DIR gmi_sim.h) + find_path(APF_SIM_INCLUDE_DIR apfSIM.h) + + list(APPEND SIMMETRIX_INCLUDE_DIR + ${GMI_SIM_INCLUDE_DIR} + ${APF_SIM_INCLUDE_DIR} + ) + + find_library(GMI_SIM_LIB gmi_sim) + find_library(APF_SIM_LIB apf_sim) +endif() + find_path(MESH_SIM_INCLUDE_DIR MeshSim.h) list(APPEND SIMMETRIX_INCLUDE_DIR - ${GMI_SIM_INCLUDE_DIR} - ${APF_SIM_INCLUDE_DIR} ${MESH_SIM_INCLUDE_DIR} ) -find_library(GMI_SIM_LIB gmi_sim) -find_library(APF_SIM_LIB apf_sim) - set(SIM_LIB_HINT ${SIMMETRIX_ROOT}/lib/x64_rhel7_gcc48) find_library(SIM_DISCRETE_LIB SimDiscrete ${SIM_LIB_HINT}) From 5dd700da89e8481c443297a05805ed1ddd0528a4 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Mon, 19 Feb 2024 11:01:35 +0100 Subject: [PATCH 26/34] Remove group if --- src/pumgen.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/pumgen.cpp b/src/pumgen.cpp index c321ec2..d485c06 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -419,13 +419,10 @@ int main(int argc, char* argv[]) { // Group information std::size_t groupBytesPerData = 4; - logInfo(rank) << "Writing group information"; - writeH5Data(meshInput->group(), h5file, "group", mesh, 3, H5T_NATIVE_INT, - H5T_STD_I32LE, chunksize, localSize[0], globalSize[0], reduceInts, - filterEnable, filterChunksize); - } else { - logInfo() << "No group information found in mesh"; - } + logInfo(rank) << "Writing group information"; + writeH5Data(meshInput->group(), h5file, "group", mesh, 3, H5T_NATIVE_INT, + H5T_STD_I32LE, chunksize, localSize[0], globalSize[0], reduceInts, + filterEnable, filterChunksize); // Write boundary condition logInfo(rank) << "Writing boundary condition"; From a41c3734da1a52a3feb7321106501a925155a58c Mon Sep 17 00:00:00 2001 From: David Schneller Date: Wed, 6 Mar 2024 10:05:08 +0100 Subject: [PATCH 27/34] Fix typos in the Netcdf reader --- src/input/NetCDFMesh.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/input/NetCDFMesh.h b/src/input/NetCDFMesh.h index 358b5ed..fac05f7 100644 --- a/src/input/NetCDFMesh.h +++ b/src/input/NetCDFMesh.h @@ -154,7 +154,7 @@ class NetCDFMesh : public FullStorageMeshData { std::size_t vertexOffset = 0; for (std::size_t i = 0; i < nLocalPart; i++) { std::copy_n(partitions[i].vertices(), partitions[i].nVertices() * 3, - geometryData.begin() + vertexOffset * 3); + verticesLocal.begin() + vertexOffset * 3); vertexOffset += partitions[i].nVertices(); } @@ -169,6 +169,7 @@ class NetCDFMesh : public FullStorageMeshData { std::copy_n(filter.localVertices().begin(), nVertices * 3, geometryData.begin()); std::size_t elementOffset = 0; + vertexOffset = 0; for (std::size_t i = 0; i < nLocalPart; i++) { #ifdef _OPENMP #pragma omp parallel for schedule(static) @@ -192,6 +193,7 @@ class NetCDFMesh : public FullStorageMeshData { groupData.begin() + elementOffset); elementOffset += partitions[i].nElements(); + vertexOffset += partitions[i].nVertices(); } logInfo(rank) << "Converting local to global vertex identifier"; From 5e96ddaa5e88be701a19c0fdf08327d7195076ec Mon Sep 17 00:00:00 2001 From: David Schneller Date: Thu, 25 Apr 2024 15:20:14 +0200 Subject: [PATCH 28/34] Implement Large-range MPI for scatter, more boundary conditions --- CMakeLists.txt | 2 +- src/aux/MPIConvenience.cpp | 80 ++++++++++++++++++++-------- src/aux/MPIConvenience.h | 4 ++ src/input/MeshData.h | 23 +++++--- src/meshreader/ParallelGMSHReader.h | 8 +-- src/pumgen.cpp | 81 ++++++++++++++++++----------- 6 files changed, 134 insertions(+), 64 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a40d19..0efd978 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,7 +25,7 @@ set(CMAKE_CXX_EXTENSIONS OFF) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "RelWithDebInfo") # MinSizeRel is useless for us if(NOT CMAKE_BUILD_TYPE) - set(CMAKE_BUILD_TYPE Release) + set(CMAKE_BUILD_TYPE Release CACHE STRING "" FORCE) message(STATUS "Set build type to Release as none was supplied.") endif() set(CMAKE_CXX_FLAGS_RELEASE "-O2") diff --git a/src/aux/MPIConvenience.cpp b/src/aux/MPIConvenience.cpp index 15d473f..4dd0022 100644 --- a/src/aux/MPIConvenience.cpp +++ b/src/aux/MPIConvenience.cpp @@ -3,6 +3,41 @@ #include #include +namespace { +constexpr std::size_t DataIncrement = std::numeric_limits::max(); +constexpr int Tag = 1000; + +std::vector largeIsend(const void* sendbuf, std::size_t sendsize, std::size_t senddisp, + MPI_Datatype sendtype, int dest, int tag, MPI_Comm comm) { + int sendtypesizePre; + std::vector requests; + MPI_Type_size(sendtype, &sendtypesizePre); + std::size_t sendtypesize = sendtypesizePre; + for (std::size_t position = 0; position < sendsize; position += DataIncrement) { + int localSize = static_cast(std::min(sendsize - position, DataIncrement)); + requests.push_back(MPI_REQUEST_NULL); + MPI_Isend(reinterpret_cast(sendbuf) + senddisp * sendtypesize + position, + localSize, sendtype, dest, tag, comm, &requests.back()); + } + return requests; +} + +std::vector largeIrecv(void* recvbuf, std::size_t recvsize, std::size_t recvdisp, + MPI_Datatype recvtype, int dest, int tag, MPI_Comm comm) { + int recvtypesizePre; + std::vector requests; + MPI_Type_size(recvtype, &recvtypesizePre); + std::size_t recvtypesize = recvtypesizePre; + for (std::size_t position = 0; position < recvsize; position += DataIncrement) { + int localSize = static_cast(std::min(recvsize - position, DataIncrement)); + requests.push_back(MPI_REQUEST_NULL); + MPI_Irecv(reinterpret_cast(recvbuf) + recvdisp * recvtypesize + position, localSize, + recvtype, dest, tag, comm, &requests.back()); + } + return requests; +} +} // namespace + void sparseAlltoallv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, MPI_Datatype sendtype, void* recvbuf, const std::size_t* recvsize, const std::size_t* recvdisp, MPI_Datatype recvtype, MPI_Comm comm) { @@ -12,36 +47,19 @@ void sparseAlltoallv(const void* sendbuf, const std::size_t* sendsize, const std int recvtypesizePre; MPI_Comm_size(comm, &commsize); MPI_Comm_rank(comm, &commrank); - MPI_Type_size(sendtype, &sendtypesizePre); - MPI_Type_size(sendtype, &recvtypesizePre); - - constexpr int Tag = 1000; - - std::size_t sendtypesize = sendtypesizePre; - std::size_t recvtypesize = recvtypesizePre; std::vector requests; requests.reserve(commsize * 2); // (no special handling of self-to-self comm at the moment) - constexpr std::size_t DataIncrement = std::numeric_limits::max(); - for (int i = 0; i < commsize; ++i) { - for (std::size_t position = 0; position < sendsize[i]; position += DataIncrement) { - int localSize = static_cast(std::min(sendsize[i] - position, DataIncrement)); - requests.push_back(MPI_REQUEST_NULL); - MPI_Isend(reinterpret_cast(sendbuf) + senddisp[i] * sendtypesize + position, - localSize, sendtype, i, Tag, comm, &requests.back()); - } + auto localRequests = largeIsend(sendbuf, sendsize[i], senddisp[i], sendtype, i, Tag, comm); + requests.insert(requests.end(), localRequests.begin(), localRequests.end()); } for (int i = 0; i < commsize; ++i) { - for (std::size_t position = 0; position < recvsize[i]; position += DataIncrement) { - int localSize = static_cast(std::min(recvsize[i] - position, DataIncrement)); - requests.push_back(MPI_REQUEST_NULL); - MPI_Irecv(reinterpret_cast(recvbuf) + recvdisp[i] * recvtypesize + position, localSize, - recvtype, i, Tag, comm, &requests.back()); - } + auto localRequests = largeIrecv(recvbuf, recvsize[i], recvdisp[i], recvtype, i, Tag, comm); + requests.insert(requests.end(), localRequests.begin(), localRequests.end()); } MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); @@ -68,3 +86,23 @@ void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddi sparseAlltoallv(sendbuf, sendsize, senddispSize.data(), sendtype, recvbuf, recvsize, recvdispSize.data(), recvtype, comm); } + +void largeScatterv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, std::size_t recvsize, + MPI_Datatype recvtype, int root, MPI_Comm comm) { + std::vector requests; + int commrank; + int commsize; + MPI_Comm_size(comm, &commsize); + MPI_Comm_rank(comm, &commrank); + if (commrank == root) { + requests.reserve(commsize + 1); + for (int i = 0; i < commsize; ++i) { + auto sendRequests = largeIsend(sendbuf, sendsize[i], senddisp[i], sendtype, i, Tag, comm); + requests.insert(requests.end(), sendRequests.begin(), sendRequests.end()); + } + } + auto recvRequests = largeIrecv(recvbuf, recvsize, 0, recvtype, root, Tag, comm); + requests.insert(requests.end(), recvRequests.begin(), recvRequests.end()); + MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); +} diff --git a/src/aux/MPIConvenience.h b/src/aux/MPIConvenience.h index a1a36bd..b757e24 100644 --- a/src/aux/MPIConvenience.h +++ b/src/aux/MPIConvenience.h @@ -16,4 +16,8 @@ void sparseAlltoallv(const void* sendbuf, const int* sendsize, const int* senddi MPI_Datatype sendtype, void* recvbuf, const int* recvsize, const int* recvdisp, MPI_Datatype recvtype, MPI_Comm comm); +void largeScatterv(const void* sendbuf, const std::size_t* sendsize, const std::size_t* senddisp, + MPI_Datatype sendtype, void* recvbuf, std::size_t recvsize, + MPI_Datatype recvtype, int root, MPI_Comm comm); + #endif // PUMGEN_AUX_MPI_CONVENIENCE_H_ diff --git a/src/input/MeshData.h b/src/input/MeshData.h index 201cc54..230c209 100644 --- a/src/input/MeshData.h +++ b/src/input/MeshData.h @@ -54,13 +54,18 @@ class FullStorageMeshData : public MeshData { int bndShift = 8; void setBoundary(std::size_t cell, int face, int value) { - constexpr auto i32limit = std::numeric_limits::max(); - constexpr auto i64limit = std::numeric_limits::max(); - if (value < 0 || (value > i32limit && bndShift == 8) || (value > i64limit && bndShift == 16)) { - logError() << "Cannot handle boundary condition" << value; + if (bndShift < 0) { + boundaryData[cell * 4 + face] = value; + } else { + constexpr auto i32limit = std::numeric_limits::max(); + constexpr auto i64limit = std::numeric_limits::max(); + if (value < 0 || (value > i32limit && bndShift == 8) || + (value > i64limit && bndShift == 16)) { + logError() << "Cannot handle boundary condition" << value; + } + + boundaryData[cell] |= static_cast(value) << (face * bndShift); } - - boundaryData[cell] |= static_cast(value) << (face * bndShift); } void setup(std::size_t cellCount, std::size_t vertexCount) { @@ -70,7 +75,11 @@ class FullStorageMeshData : public MeshData { connectivityData.resize(cellCount * 4); geometryData.resize(vertexCount * 3); groupData.resize(cellCount); - boundaryData.resize(cellCount); + if (bndShift < 0) { + boundaryData.resize(cellCount * 4); + } else { + boundaryData.resize(cellCount); + } // TODO: reconsider the times 4 or 3 } diff --git a/src/meshreader/ParallelGMSHReader.h b/src/meshreader/ParallelGMSHReader.h index eed6677..47433b1 100644 --- a/src/meshreader/ParallelGMSHReader.h +++ b/src/meshreader/ParallelGMSHReader.h @@ -142,8 +142,8 @@ template class ParallelGMSHReader { MPI_Comm_rank(comm_, &rank); MPI_Comm_size(comm_, &procs); - auto sendcounts = std::vector(procs); - auto displs = std::vector(procs + 1); + auto sendcounts = std::vector(procs); + auto displs = std::vector(procs + 1); displs[0] = 0; for (int i = 0; i < procs; ++i) { sendcounts[i] = numPerElement * getChunksize(numElements, i, procs); @@ -151,8 +151,8 @@ template class ParallelGMSHReader { } auto recvcount = sendcounts[rank]; - MPI_Scatterv(sendbuf, sendcounts.data(), displs.data(), tndm::mpi_type_t(), recvbuf, - recvcount, tndm::mpi_type_t(), 0, comm_); + largeScatterv(sendbuf, sendcounts.data(), displs.data(), tndm::mpi_type_t(), recvbuf, + recvcount, tndm::mpi_type_t(), 0, comm_); } MPI_Comm comm_; diff --git a/src/pumgen.cpp b/src/pumgen.cpp index d485c06..e825338 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -75,27 +75,27 @@ static int ilog(std::size_t value, int expbase = 1) { constexpr std::size_t NoSecondDim = 0; -template +template static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, void* mesh, int meshdim, hid_t h5memtype, hid_t h5outtype, std::size_t chunk, std::size_t localSize, std::size_t globalSize, bool reduceInts, - int filterEnable, std::size_t filterChunksize) { - constexpr std::size_t SecondSize = std::max(SecondDim, static_cast(1)); - constexpr std::size_t Dimensions = SecondDim == 0 ? 1 : 2; - const std::size_t chunkSize = chunk / SecondSize / sizeof(T); + int filterEnable, std::size_t filterChunksize, std::size_t secondDim) { + const std::size_t secondSize = std::max(secondDim, static_cast(1)); + const std::size_t dimensions = secondDim == 0 ? 1 : 2; + const std::size_t chunkSize = chunk / secondSize / sizeof(T); const std::size_t bufferSize = std::min(localSize, chunkSize); - std::vector data(SecondSize * bufferSize); + std::vector data(secondSize * bufferSize); std::size_t rounds = (localSize + chunkSize - 1) / chunkSize; MPI_Allreduce(MPI_IN_PLACE, &rounds, 1, tndm::mpi_type_t(), MPI_MAX, MPI_COMM_WORLD); - hsize_t globalSizes[2] = {globalSize, SecondDim}; - hid_t h5space = H5Screate_simple(Dimensions, globalSizes, nullptr); + hsize_t globalSizes[2] = {globalSize, secondDim}; + hid_t h5space = H5Screate_simple(dimensions, globalSizes, nullptr); - hsize_t sizes[2] = {bufferSize, SecondDim}; - hid_t h5memspace = H5Screate_simple(Dimensions, sizes, nullptr); + hsize_t sizes[2] = {bufferSize, secondDim}; + hid_t h5memspace = H5Screate_simple(dimensions, sizes, nullptr); std::size_t offset = localSize; @@ -104,7 +104,7 @@ static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, offset -= localSize; hsize_t start[2] = {offset, 0}; - hsize_t count[2] = {bufferSize, SecondDim}; + hsize_t count[2] = {bufferSize, secondDim}; hid_t h5dxlist = H5Pcreate(H5P_DATASET_XFER); checkH5Err(h5dxlist); @@ -124,7 +124,7 @@ static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, hid_t h5filter = H5P_DEFAULT; if (filterEnable > 0) { h5filter = checkH5Err(H5Pcreate(H5P_DATASET_CREATE)); - hsize_t chunk[2] = {std::min(filterChunksize, bufferSize), SecondDim}; + hsize_t chunk[2] = {std::min(filterChunksize, bufferSize), secondDim}; checkH5Err(H5Pset_chunk(h5filter, 2, chunk)); if (filterEnable == 1 && std::is_integral_v) { checkH5Err(H5Pset_scaleoffset(h5filter, H5Z_SO_INT, H5Z_SO_INT_MINBITS_DEFAULT)); @@ -147,7 +147,7 @@ static void writeH5Data(const F& handler, hid_t h5file, const std::string& name, start[0] = offset + written; count[0] = std::min(localSize - written, bufferSize); - std::copy_n(handler.begin() + written * SecondSize, count[0] * SecondSize, data.begin()); + std::copy_n(handler.begin() + written * secondSize, count[0] * secondSize, data.begin()); checkH5Err(H5Sselect_hyperslab(h5memspace, H5S_SELECT_SET, nullstart, nullptr, count, nullptr)); @@ -202,9 +202,9 @@ int main(int argc, char* argv[]) { utils::Args::Required, false); args.addOption("chunksize", 0, "Chunksize for writing (default=1 GiB).", utils::Args::Required, false); - const char* boundarytypes[] = {"int32", "int64"}; + const char* boundarytypes[] = {"int32", "int64", "i32", "i64", "i32x4"}; args.addEnumOption("boundarytype", boundarytypes, 0, - "Type for writing out boundary data (default: int32).", false); + "Type for writing out boundary data (default: i32).", false); args.addOption("license", 'l', "License file (only used by SimModSuite)", utils::Args::Required, false); @@ -272,14 +272,26 @@ int main(int argc, char* argv[]) { int boundaryType = args.getArgument("boundarytype", 0); hid_t boundaryDatatype; int faceOffset; - if (boundaryType == 0) { + int secondShape = NoSecondDim; + int boundaryFormatAttr = 0; + if (boundaryType == 0 || boundaryType == 2) { boundaryDatatype = H5T_STD_I32LE; faceOffset = 8; - logInfo(rank) << "Using 32-bit integer boundary type conditions."; - } else if (boundaryType == 1) { + secondShape = NoSecondDim; + boundaryFormatAttr = 1; + logInfo(rank) << "Using 32-bit integer boundary type conditions, or 8 bit per face (i32)."; + } else if (boundaryType == 1 || boundaryType == 3) { boundaryDatatype = H5T_STD_I64LE; faceOffset = 16; - logInfo(rank) << "Using 64-bit integer boundary type conditions."; + secondShape = NoSecondDim; + boundaryFormatAttr = 2; + logInfo(rank) << "Using 64-bit integer boundary type conditions, or 16 bit per face (i64)."; + } else if (boundaryType == 4) { + boundaryDatatype = H5T_STD_I32LE; + secondShape = 4; + faceOffset = -1; + boundaryFormatAttr = 0; + logInfo(rank) << "Using 32-bit integer per boundary face (i32x4)."; } std::string xdmfFile = outputFile; @@ -406,31 +418,38 @@ int main(int argc, char* argv[]) { // Write cells std::size_t connectBytesPerData = 8; logInfo(rank) << "Writing cells"; - writeH5Data(meshInput->connectivity(), h5file, "connect", mesh, 3, - H5T_NATIVE_ULONG, H5T_STD_U64LE, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->connectivity(), h5file, "connect", mesh, 3, + H5T_NATIVE_ULONG, H5T_STD_U64LE, chunksize, localSize[0], + globalSize[0], reduceInts, filterEnable, filterChunksize, 4); // Vertices logInfo(rank) << "Writing vertices"; - writeH5Data(meshInput->geometry(), h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, - H5T_IEEE_F64LE, chunksize, localSize[1], globalSize[1], reduceInts, - filterEnable, filterChunksize); + writeH5Data(meshInput->geometry(), h5file, "geometry", mesh, 0, H5T_IEEE_F64LE, + H5T_IEEE_F64LE, chunksize, localSize[1], globalSize[1], reduceInts, + filterEnable, filterChunksize, 3); // Group information std::size_t groupBytesPerData = 4; logInfo(rank) << "Writing group information"; - writeH5Data(meshInput->group(), h5file, "group", mesh, 3, H5T_NATIVE_INT, - H5T_STD_I32LE, chunksize, localSize[0], globalSize[0], reduceInts, - filterEnable, filterChunksize); + writeH5Data(meshInput->group(), h5file, "group", mesh, 3, H5T_NATIVE_INT, H5T_STD_I32LE, + chunksize, localSize[0], globalSize[0], reduceInts, filterEnable, + filterChunksize, NoSecondDim); // Write boundary condition logInfo(rank) << "Writing boundary condition"; auto i32limit = std::numeric_limits::max(); auto i64limit = std::numeric_limits::max(); - writeH5Data(meshInput->boundary(), h5file, "boundary", mesh, 3, - H5T_NATIVE_LONG, boundaryDatatype, chunksize, localSize[0], - globalSize[0], reduceInts, filterEnable, filterChunksize); + writeH5Data(meshInput->boundary(), h5file, "boundary", mesh, 3, H5T_NATIVE_LONG, + boundaryDatatype, chunksize, localSize[0], globalSize[0], reduceInts, + filterEnable, filterChunksize, secondShape); + + hid_t attrSpace = checkH5Err(H5Screate(H5S_SCALAR)); + hid_t attrBoundary = checkH5Err( + H5Acreate(h5file, "boundary-format", H5T_NATIVE_INT, attrSpace, H5P_DEFAULT, H5P_DEFAULT)); + checkH5Err(H5Awrite(attrBoundary, H5T_NATIVE_INT, &boundaryFormatAttr)); + checkH5Err(H5Aclose(attrBoundary)); + checkH5Err(H5Sclose(attrSpace)); // Writing XDMF file if (rank == 0) { From 81317a2ded094f00e4c3a6a4e8ca8430298e5bb2 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Thu, 25 Apr 2024 15:28:33 +0200 Subject: [PATCH 29/34] Add missing header --- src/meshreader/ParallelGMSHReader.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/meshreader/ParallelGMSHReader.h b/src/meshreader/ParallelGMSHReader.h index 47433b1..99388b8 100644 --- a/src/meshreader/ParallelGMSHReader.h +++ b/src/meshreader/ParallelGMSHReader.h @@ -12,6 +12,7 @@ #include #include "aux/Distributor.h" +#include "aux/MPIConvenience.h" namespace puml { From 88d9533784c833871f24c91d65180c5487bad42b Mon Sep 17 00:00:00 2001 From: David Schneller Date: Thu, 25 Apr 2024 15:39:30 +0200 Subject: [PATCH 30/34] Fix position increment --- src/aux/MPIConvenience.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aux/MPIConvenience.cpp b/src/aux/MPIConvenience.cpp index 4dd0022..168a650 100644 --- a/src/aux/MPIConvenience.cpp +++ b/src/aux/MPIConvenience.cpp @@ -16,7 +16,7 @@ std::vector largeIsend(const void* sendbuf, std::size_t sendsize, s for (std::size_t position = 0; position < sendsize; position += DataIncrement) { int localSize = static_cast(std::min(sendsize - position, DataIncrement)); requests.push_back(MPI_REQUEST_NULL); - MPI_Isend(reinterpret_cast(sendbuf) + senddisp * sendtypesize + position, + MPI_Isend(reinterpret_cast(sendbuf) + (senddisp + position) * sendtypesize, localSize, sendtype, dest, tag, comm, &requests.back()); } return requests; @@ -31,7 +31,7 @@ std::vector largeIrecv(void* recvbuf, std::size_t recvsize, std::si for (std::size_t position = 0; position < recvsize; position += DataIncrement) { int localSize = static_cast(std::min(recvsize - position, DataIncrement)); requests.push_back(MPI_REQUEST_NULL); - MPI_Irecv(reinterpret_cast(recvbuf) + recvdisp * recvtypesize + position, localSize, + MPI_Irecv(reinterpret_cast(recvbuf) + (recvdisp + position) * recvtypesize, localSize, recvtype, dest, tag, comm, &requests.back()); } return requests; From dd32ce738c9cf6e455a10fa29fcfa984fbe8c99a Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 7 May 2024 10:28:51 +0200 Subject: [PATCH 31/34] Update clang-format pipeline --- .github/workflows/cmake.yml | 3 ++- src/input/MeshAttributes.cpp | 5 ++--- src/input/SimModSuite.h | 3 +-- src/input/SimModSuiteApf.h | 3 +-- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index dfec676..70325cf 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -87,6 +87,7 @@ jobs: - name: clang-format run: | sudo apt-get update - sudo apt-get install -qq python3 clang-format + sudo apt-get install -qq python3 python3-pip + pip3 install clang-format git submodule update --init ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format --exclude src/third_party -r src diff --git a/src/input/MeshAttributes.cpp b/src/input/MeshAttributes.cpp index de77884..da316bf 100644 --- a/src/input/MeshAttributes.cpp +++ b/src/input/MeshAttributes.cpp @@ -210,9 +210,8 @@ void MeshAttributes::set_velocity_aware_meshing() { logInfo(PMU_rank()) << "Adding velocity aware refinement region targeting" << targetedFrequency << "Hz, centered at x =" << cuboid.center[0] << "y=" << cuboid.center[1] << "z=" << cuboid.center[2] - << "with half sizes" - << "x =" << cuboid.halfSize[0] << "y =" << cuboid.halfSize[1] - << "z =" << cuboid.halfSize[2]; + << "with half sizes" << "x =" << cuboid.halfSize[0] + << "y =" << cuboid.halfSize[1] << "z =" << cuboid.halfSize[2]; if (std::abs(cuboid.rotationZ) > 0.0) { logInfo(PMU_rank()) << "rotated around z axis by " << cuboid.rotationZ << "degree(s) counterclockwise from x axis."; diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index 416448a..6de183f 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -426,8 +426,7 @@ class SimModSuite : public FullStorageMeshData { switch (currentVal) { case -2: // task is started, do nothing - logInfo(PMU_rank()) << "Progress:" << what << ", 0" - << "/" << endVal; + logInfo(PMU_rank()) << "Progress:" << what << ", 0" << "/" << endVal; break; case -1: // end of the task diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h index ad07f30..3b40e56 100644 --- a/src/input/SimModSuiteApf.h +++ b/src/input/SimModSuiteApf.h @@ -308,8 +308,7 @@ class SimModSuiteApf : public ApfMeshInput { switch (currentVal) { case -2: // task is started, do nothing - logInfo(PMU_rank()) << "Progress:" << what << ", 0" - << "/" << endVal; + logInfo(PMU_rank()) << "Progress:" << what << ", 0" << "/" << endVal; break; case -1: // end of the task From 990652340be814a4126d6e49702b6d5ba8e5df0b Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 7 May 2024 11:48:55 +0200 Subject: [PATCH 32/34] Add boundary format as attribute --- src/input/AnalysisAttributes.h | 2 +- src/input/MeshAttributes.h | 2 +- src/pumgen.cpp | 28 +++++++++++++++++++--------- 3 files changed, 21 insertions(+), 11 deletions(-) diff --git a/src/input/AnalysisAttributes.h b/src/input/AnalysisAttributes.h index e7d52b0..ed91eee 100644 --- a/src/input/AnalysisAttributes.h +++ b/src/input/AnalysisAttributes.h @@ -12,7 +12,7 @@ using namespace tinyxml2; struct faceBoundary { - faceBoundary(int faceID, int bcType) : faceID(faceID), bcType(bcType){}; + faceBoundary(int faceID, int bcType) : faceID(faceID), bcType(bcType) {}; int faceID; int bcType; }; diff --git a/src/input/MeshAttributes.h b/src/input/MeshAttributes.h index 591e99c..86fe757 100644 --- a/src/input/MeshAttributes.h +++ b/src/input/MeshAttributes.h @@ -33,7 +33,7 @@ struct VelocityRefinementCube { VelocityRefinementCube(SimpleCuboid cuboid, double targetedFrequency, int bypassFindRegionAndUseGroup) : cuboid(cuboid), targetedFrequency(targetedFrequency), - bypassFindRegionAndUseGroup(bypassFindRegionAndUseGroup){}; + bypassFindRegionAndUseGroup(bypassFindRegionAndUseGroup) {}; SimpleCuboid cuboid; double targetedFrequency; diff --git a/src/pumgen.cpp b/src/pumgen.cpp index e825338..98a5a97 100644 --- a/src/pumgen.cpp +++ b/src/pumgen.cpp @@ -10,6 +10,7 @@ * @author Sebastian Rettenberger */ +#include #include #include @@ -273,24 +274,30 @@ int main(int argc, char* argv[]) { hid_t boundaryDatatype; int faceOffset; int secondShape = NoSecondDim; - int boundaryFormatAttr = 0; + std::string boundaryFormatAttr = ""; + int boundaryPrecision = 0; + std::string secondDimBoundary = ""; if (boundaryType == 0 || boundaryType == 2) { boundaryDatatype = H5T_STD_I32LE; faceOffset = 8; secondShape = NoSecondDim; - boundaryFormatAttr = 1; + boundaryFormatAttr = "i32"; + boundaryPrecision = 4; logInfo(rank) << "Using 32-bit integer boundary type conditions, or 8 bit per face (i32)."; } else if (boundaryType == 1 || boundaryType == 3) { boundaryDatatype = H5T_STD_I64LE; faceOffset = 16; secondShape = NoSecondDim; - boundaryFormatAttr = 2; + boundaryFormatAttr = "i64"; + boundaryPrecision = 8; logInfo(rank) << "Using 64-bit integer boundary type conditions, or 16 bit per face (i64)."; } else if (boundaryType == 4) { boundaryDatatype = H5T_STD_I32LE; secondShape = 4; faceOffset = -1; - boundaryFormatAttr = 0; + boundaryFormatAttr = "i32x4"; + boundaryPrecision = 4; + secondDimBoundary = " 4"; logInfo(rank) << "Using 32-bit integer per boundary face (i32x4)."; } @@ -445,11 +452,14 @@ int main(int argc, char* argv[]) { filterEnable, filterChunksize, secondShape); hid_t attrSpace = checkH5Err(H5Screate(H5S_SCALAR)); + hid_t attrType = checkH5Err(H5Tcopy(H5T_C_S1)); + checkH5Err(H5Tset_size(attrType, boundaryFormatAttr.size())); hid_t attrBoundary = checkH5Err( - H5Acreate(h5file, "boundary-format", H5T_NATIVE_INT, attrSpace, H5P_DEFAULT, H5P_DEFAULT)); - checkH5Err(H5Awrite(attrBoundary, H5T_NATIVE_INT, &boundaryFormatAttr)); + H5Acreate(h5file, "boundary-format", attrType, attrSpace, H5P_DEFAULT, H5P_DEFAULT)); + checkH5Err(H5Awrite(attrBoundary, attrType, boundaryFormatAttr.data())); checkH5Err(H5Aclose(attrBoundary)); checkH5Err(H5Sclose(attrSpace)); + checkH5Err(H5Tclose(attrType)); // Writing XDMF file if (rank == 0) { @@ -491,9 +501,9 @@ int main(int argc, char* argv[]) { << globalSize[0] << "\">" << basename << ":/group" << std::endl << " " << std::endl << " " << std::endl - << " " << basename - << ":/boundary" << std::endl + << " " + << basename << ":/boundary" << std::endl << " " << std::endl << " " << std::endl << " " << std::endl From 5e9758a80e527b5301a817fe5d129e660cd14c84 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 7 May 2024 11:53:57 +0200 Subject: [PATCH 33/34] Formatting again, fix clang-format version --- .github/workflows/cmake.yml | 2 +- src/input/MeshAttributes.cpp | 5 +++-- src/input/SimModSuite.h | 3 ++- src/input/SimModSuiteApf.h | 3 ++- 4 files changed, 8 insertions(+), 5 deletions(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 70325cf..a336d4d 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -88,6 +88,6 @@ jobs: run: | sudo apt-get update sudo apt-get install -qq python3 python3-pip - pip3 install clang-format + pip3 install clang-format==18.1.5 git submodule update --init ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format --exclude src/third_party -r src diff --git a/src/input/MeshAttributes.cpp b/src/input/MeshAttributes.cpp index da316bf..de77884 100644 --- a/src/input/MeshAttributes.cpp +++ b/src/input/MeshAttributes.cpp @@ -210,8 +210,9 @@ void MeshAttributes::set_velocity_aware_meshing() { logInfo(PMU_rank()) << "Adding velocity aware refinement region targeting" << targetedFrequency << "Hz, centered at x =" << cuboid.center[0] << "y=" << cuboid.center[1] << "z=" << cuboid.center[2] - << "with half sizes" << "x =" << cuboid.halfSize[0] - << "y =" << cuboid.halfSize[1] << "z =" << cuboid.halfSize[2]; + << "with half sizes" + << "x =" << cuboid.halfSize[0] << "y =" << cuboid.halfSize[1] + << "z =" << cuboid.halfSize[2]; if (std::abs(cuboid.rotationZ) > 0.0) { logInfo(PMU_rank()) << "rotated around z axis by " << cuboid.rotationZ << "degree(s) counterclockwise from x axis."; diff --git a/src/input/SimModSuite.h b/src/input/SimModSuite.h index 6de183f..416448a 100644 --- a/src/input/SimModSuite.h +++ b/src/input/SimModSuite.h @@ -426,7 +426,8 @@ class SimModSuite : public FullStorageMeshData { switch (currentVal) { case -2: // task is started, do nothing - logInfo(PMU_rank()) << "Progress:" << what << ", 0" << "/" << endVal; + logInfo(PMU_rank()) << "Progress:" << what << ", 0" + << "/" << endVal; break; case -1: // end of the task diff --git a/src/input/SimModSuiteApf.h b/src/input/SimModSuiteApf.h index 3b40e56..ad07f30 100644 --- a/src/input/SimModSuiteApf.h +++ b/src/input/SimModSuiteApf.h @@ -308,7 +308,8 @@ class SimModSuiteApf : public ApfMeshInput { switch (currentVal) { case -2: // task is started, do nothing - logInfo(PMU_rank()) << "Progress:" << what << ", 0" << "/" << endVal; + logInfo(PMU_rank()) << "Progress:" << what << ", 0" + << "/" << endVal; break; case -1: // end of the task From 8fb55e3087270789fc8c174c793127b921d10861 Mon Sep 17 00:00:00 2001 From: David Schneller Date: Tue, 7 May 2024 12:01:39 +0200 Subject: [PATCH 34/34] Remove old CI pipelines --- .travis.yml | 37 ----------------------- azure-pipelines-build-template.yml | 26 ---------------- azure-pipelines-clang-format-template.yml | 13 -------- azure-pipelines.yml | 3 -- 4 files changed, 79 deletions(-) delete mode 100644 .travis.yml delete mode 100644 azure-pipelines-build-template.yml delete mode 100644 azure-pipelines-clang-format-template.yml delete mode 100644 azure-pipelines.yml diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index cb0066d..0000000 --- a/.travis.yml +++ /dev/null @@ -1,37 +0,0 @@ -sudo: required - -notifications: - email: - on_success: always - on_failure: always - -language: cpp - -compiler: gcc - -dist: bionic - -before_install: - - sudo add-apt-repository -y ppa:ubuntu-toolchain-r/test - - sudo apt-get update -qq - -install: - - sudo apt-get install -qq g++-5 openmpi-bin openmpi-common libopenmpi-dev hdf5-tools libhdf5-openmpi-100 libhdf5-openmpi-dev libnetcdf-dev clang-format-10 - - cd .. -# Download simmetrix libraries -# check license -#build pumi - - git clone https://github.com/SCOREC/core.git - - mkdir pumi - - mkdir core/build - - cd core/build - - cmake .. -DCMAKE_C_COMPILER="`which mpicc`" -DCMAKE_CXX_COMPILER="`which mpiCC`" -DCMAKE_C_FLAGS="-O2 -g -Wall" -DCMAKE_CXX_FLAGS="-O2 -g -Wall" -DCMAKE_INSTALL_PREFIX=../../pumi - - make -j2 && make install - - make -j2 && make install -#Go back to PUMGen - - cd ../../PUMGen - -script: - - ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format-10 --exclude src/third_party -r src - - mkdir build && cd build - - CC=mpicc CXX=mpiCC cmake .. -DCMAKE_PREFIX_PATH=../pumi && make -j2 diff --git a/azure-pipelines-build-template.yml b/azure-pipelines-build-template.yml deleted file mode 100644 index bf89ddb..0000000 --- a/azure-pipelines-build-template.yml +++ /dev/null @@ -1,26 +0,0 @@ -jobs: - - job: - displayName: pumgen-build - pool: - vmImage: 'ubuntu-22.04' - steps: - - bash: | - set -euo pipefail - export IFS=$'\n\t' - sudo apt-get update - sudo apt-get install -qq g++ openmpi-bin openmpi-common libopenmpi-dev hdf5-tools libhdf5-openmpi-103 libhdf5-openmpi-dev libnetcdf-dev - cd .. - # build pumi - git clone --recursive --branch v2.2.7 https://github.com/SCOREC/core.git - cd core - git --no-pager log --pretty=oneline -1 - mkdir build - cd build - cmake .. -DCMAKE_C_COMPILER="`which mpicc`" -DCMAKE_CXX_COMPILER="`which mpiCC`" -DCMAKE_C_FLAGS="-O2 -g -Wall" -DCMAKE_CXX_FLAGS="-O2 -g -Wall" -DSCOREC_CXX_WARNINGS=OFF -DCMAKE_INSTALL_PREFIX=../../pumi - make -j $(nproc) && make install - #Go back to PUMGen and build it - cd ../../s - git submodule update --init - mkdir build && cd build - CC=mpicc CXX=mpiCC cmake .. -DCMAKE_PREFIX_PATH=../pumi && make -j $(nproc) - diff --git a/azure-pipelines-clang-format-template.yml b/azure-pipelines-clang-format-template.yml deleted file mode 100644 index 9cb6996..0000000 --- a/azure-pipelines-clang-format-template.yml +++ /dev/null @@ -1,13 +0,0 @@ -jobs: - - job: - displayName: pumgen-clang-format - pool: - vmImage: 'ubuntu-22.04' - steps: - - bash: | - set -euo pipefail - export IFS=$'\n\t' - sudo apt-get update - sudo apt-get install -qq python3 clang-format-14 - git submodule update --init - ./submodules/run-clang-format/run-clang-format.py --clang-format-executable clang-format-14 --exclude src/third_party -r src diff --git a/azure-pipelines.yml b/azure-pipelines.yml deleted file mode 100644 index 997c5bf..0000000 --- a/azure-pipelines.yml +++ /dev/null @@ -1,3 +0,0 @@ -jobs: - - template: azure-pipelines-clang-format-template.yml - - template: azure-pipelines-build-template.yml