Faster and parallel reading; velocity-aware GMSH preprocessing

.github/workflows/cmake.yml

@@ -29,12 +29,16 @@ jobs:
         run: |
           git submodule update --init --recursive
           mkdir build_${{ matrix.build_type}} && cd build_${{ matrix.build_type}}
-          cmake .. -DCMAKE_PREFIX_PATH=${{ runner.workspace }}/opt -DCMAKE_BUILD_TYPE=${{ matrix.build_type }}
+          cmake .. -DCMAKE_PREFIX_PATH=${{ runner.workspace }}/opt -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} -DMPIEXEC_PREFLAGS=--oversubscribe
 
       - name: Build
         working-directory: ${{ github.workspace }}/build_${{ matrix.build_type }}
         run: cmake --build .
 
+      - name: Test
+        working-directory: ${{ github.workspace }}/build_${{ matrix.build_type }}
+        run: ctest --output-on-failure
+
   pre-commit-check:
     runs-on: ubuntu-latest
     steps:

.pre-commit-config.yaml

@@ -59,5 +59,7 @@ repos:
     hooks:
     -   id: end-of-file-fixer
         name: '[GENERIC] newline eof'
+        exclude: ^tests/fixtures/
     -   id: trailing-whitespace
         name: '[GENERIC] remove trailing whitespace'
+        exclude: ^tests/fixtures/

CMakeLists.txt

@@ -13,6 +13,7 @@ set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fno-omit-frame-pointer -fsa
 set(CMAKE_LINKER_FLAGS_DEBUG "${CMAKE_LINKER_FLAGS_DEBUG} -fno-omit-frame-pointer -fsanitize=address")
 
 option(SIMMETRIX "Use simmetrix libraries" OFF)
+option(EASI "Build the velocity-aware sizing with easi (always on with SIMMETRIX)" OFF)
 set(SIMMETRIX_ROOT "" CACHE STRING "Root directory of simmetrix dev files")
 set(SIM_MPI "" CACHE STRING "MPI version used by simmetrix")
 option(NETCDF "Use netcdf" OFF)
@@ -37,11 +38,13 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O2")
 add_executable(pumgen
   ${CMAKE_CURRENT_SOURCE_DIR}/src/pumgen.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GambitReader.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GMSHParser.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GMSH2Parser.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/GMSH4Parser.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/MshInput.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/Msh4Index.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/meshreader/DistributedGMSHReader.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/input/ParallelVertexFilter.cpp
-  ${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/helper/InsphereCalculator.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/src/helper/MPIConvenience.cpp
   )
@@ -93,10 +96,6 @@ target_link_libraries(pumgen PUBLIC tinyxml2)
 if (SIMMETRIX)
     find_package(SIMMETRIX REQUIRED)
 
-    # easi
-    find_package(easi 1.0.0 REQUIRED)
-    target_link_libraries(pumgen PUBLIC easi::easi)
-
     target_sources(pumgen PUBLIC
       ${CMAKE_CURRENT_SOURCE_DIR}/src/input/AnalysisAttributes.cpp
       ${CMAKE_CURRENT_SOURCE_DIR}/src/input/EasiMeshSize.cpp
@@ -134,8 +133,40 @@ endfunction()
 
 cast_log_level_to_int(LOG_LEVEL LOG_LEVEL)
 target_compile_definitions(pumgen PUBLIC LOG_LEVEL=${LOG_LEVEL})
+
+if (SIMMETRIX)
+  set(EASI ON)
+endif()
+if (EASI)
+  find_package(easi 1.0.0 REQUIRED)
+  add_library(pumgen_sizing STATIC
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/sizing/VelocityAwareSettings.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/sizing/VelocityAwareMeshSize.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/sizing/SizeField.cpp
+  )
+  target_include_directories(pumgen_sizing PUBLIC
+    ${CMAKE_CURRENT_SOURCE_DIR}/src
+    ${CMAKE_CURRENT_SOURCE_DIR}/submodules
+  )
+  target_compile_definitions(pumgen_sizing PRIVATE LOG_LEVEL=${LOG_LEVEL})
+  target_link_libraries(pumgen_sizing PUBLIC easi::easi tinyxml2)
+  target_link_libraries(pumgen PUBLIC pumgen_sizing)
+  target_sources(pumgen PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src/sizing/VelocityCheck.cpp)
+  target_compile_definitions(pumgen PUBLIC USE_EASI)
+
+  add_executable(pumgen-sizefield ${CMAKE_CURRENT_SOURCE_DIR}/src/sizefield.cpp)
+  target_compile_definitions(pumgen-sizefield PRIVATE LOG_LEVEL=${LOG_LEVEL})
+  target_link_libraries(pumgen-sizefield PRIVATE pumgen_sizing MPI::MPI_CXX)
+  set_target_properties(pumgen-sizefield PROPERTIES INSTALL_RPATH "${INSTALL_RPATH}")
+  install(TARGETS pumgen-sizefield RUNTIME DESTINATION bin)
+endif()
 if(SIMMETRIX)
   set(INSTALL_RPATH "${SIM_PS_KRNL_LIB_DIR}:${INSTALL_RPATH}")
 endif()
 set_target_properties(pumgen PROPERTIES INSTALL_RPATH "${INSTALL_RPATH}")
 install(TARGETS pumgen RUNTIME DESTINATION bin DESTINATION ${CMAKE_INSTALL_BINDIR})
+
+include(CTest)
+if (BUILD_TESTING)
+  add_subdirectory(tests)
+endif()

REUSE.toml

@@ -0,0 +1,12 @@
+# SPDX-FileCopyrightText: 2026 SeisSol Group
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+version = 1
+
+[[annotations]]
+path = ["tests/fixtures/*.msh", "tests/fixtures/*.neu", "tests/fixtures/*.h5",
+        "tests/fixtures/*.bin"]
+precedence = "override"
+SPDX-FileCopyrightText = "2026 SeisSol Group"
+SPDX-License-Identifier = "BSD-3-Clause"

src/helper/Distributor.h

@@ -4,7 +4,9 @@
 #ifndef PUMGEN_SRC_HELPER_DISTRIBUTOR_H_
 #define PUMGEN_SRC_HELPER_DISTRIBUTOR_H_
 
+#include <algorithm>
 #include <cmath>
+#include <cstddef>
 #include <cstdlib>
 
 // no namespace for now
@@ -20,4 +22,17 @@ constexpr std::size_t getChunksum(std::size_t total, int until, int size) {
   return base + addon;
 }
 
+/**
+ * The rank whose chunk (as given by getChunksize/getChunksum) contains the item index.
+ */
+constexpr int getChunkOwner(std::size_t total, std::size_t index, int size) {
+  const std::size_t base = total / size;
+  const std::size_t rest = total % size;
+  const std::size_t largeChunks = rest * (base + 1);
+  if (index < largeChunks) {
+    return static_cast<int>(index / (base + 1));
+  }
+  return static_cast<int>(rest + (index - largeChunks) / base);
+}
+
 #endif // PUMGEN_SRC_HELPER_DISTRIBUTOR_H_

src/helper/InsphereCalculator.cpp

@@ -12,11 +12,10 @@
 
 #include "MPIConvenience.h"
 
-std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivity,
-                                      const std::vector<double>& geometry, MPI_Comm comm) {
+CellVertices::CellVertices(const std::vector<std::size_t>& connectivity,
+                           const std::vector<double>& geometry, std::size_t cellSize, MPI_Comm comm)
+    : connectivity(connectivity), geometry(geometry), cellSize(cellSize) {
   int commsize;
-  int commrank;
-
   MPI_Comm_size(comm, &commsize);
   MPI_Comm_rank(comm, &commrank);
 
@@ -30,7 +29,7 @@ std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivi
 
   std::size_t localVertices = geometry.size() / 3;
 
-  std::vector<std::size_t> vertexDist(commsize + 1);
+  vertexDist.resize(commsize + 1);
 
   MPI_Allgather(&localVertices, 1, tndm::mpi_type_t<std::size_t>(), vertexDist.data() + 1, 1,
                 tndm::mpi_type_t<std::size_t>(), comm);
@@ -39,9 +38,13 @@ std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivi
     vertexDist[i + 1] += vertexDist[i];
   }
 
-  std::vector<std::unordered_map<std::size_t, std::size_t>> outidxmap(commsize);
+  outidxmap.resize(commsize);
 
-  for (const auto& vertex : connectivity) {
+  for (std::size_t node = 0; node < connectivity.size(); ++node) {
+    if (node % cellSize >= VerticesPerCell) {
+      continue;
+    }
+    const auto vertex = connectivity[node];
     auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex);
     auto position = std::distance(vertexDist.begin(), itPosition) - 1;
     auto localVertex = vertex - vertexDist[position];
@@ -57,7 +60,7 @@ std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivi
   MPI_Alltoall(outrequests.data(), 1, tndm::mpi_type_t<std::size_t>(), inrequests.data(), 1,
                tndm::mpi_type_t<std::size_t>(), comm);
 
-  std::vector<std::size_t> outdisp(commsize + 1);
+  outdisp.resize(commsize + 1);
   std::vector<std::size_t> indisp(commsize + 1);
   for (std::size_t i = 1; i < commsize + 1; ++i) {
     outdisp[i] = outdisp[i - 1] + outrequests[i - 1];
@@ -84,7 +87,7 @@ std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivi
                     tndm::mpi_type_t<std::size_t>(), comm);
   }
 
-  std::vector<double> outvertices(3 * connectivity.size());
+  outvertices.resize(3 * outdisp[commsize]);
 
   {
     std::vector<double> invertices(3 * indisp[commsize]);
@@ -101,27 +104,36 @@ std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivi
                     outvertices.data(), outrequests.data(), outdisp.data(), vertexType, comm);
   }
 
-  std::vector<std::size_t> counter(commsize);
-  std::vector<double> inspheres(connectivity.size() / 4);
-
-  for (std::size_t i = 0; i < connectivity.size() / 4; ++i) {
-    std::array<std::array<double, 3>, 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;
-      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];
-      }
+  MPI_Type_free(&vertexType);
+}
+
+std::array<std::array<double, 3>, CellVertices::VerticesPerCell>
+CellVertices::operator()(std::size_t cell) const {
+  std::array<std::array<double, 3>, VerticesPerCell> vertices;
+  for (std::size_t j = 0; j < VerticesPerCell; ++j) {
+    auto vertex = connectivity[cell * cellSize + j];
+    auto itPosition = std::upper_bound(vertexDist.begin(), vertexDist.end(), vertex);
+    auto position = std::distance(vertexDist.begin(), itPosition) - 1;
+    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].at(localVertex) + outdisp[position];
+      vertices[j][0] = outvertices[transferidx * 3 + 0];
+      vertices[j][1] = outvertices[transferidx * 3 + 1];
+      vertices[j][2] = outvertices[transferidx * 3 + 2];
     }
+  }
+  return vertices;
+}
+
+std::vector<double> calculateInsphere(const CellVertices& cells) {
+  std::vector<double> inspheres(cells.numCells());
+
+  for (std::size_t i = 0; i < cells.numCells(); ++i) {
+    const auto vertices = cells(i);
     double a11 = vertices[1][0] - vertices[0][0];
     double a12 = vertices[1][1] - vertices[0][1];
     double a13 = vertices[1][2] - vertices[0][2];
@@ -162,7 +174,5 @@ std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivi
     inspheres[i] = gram / faces;
   }
 
-  MPI_Type_free(&vertexType);
-
   return inspheres;
 }

src/helper/InsphereCalculator.h

@@ -4,13 +4,43 @@
 #ifndef PUMGEN_SRC_HELPER_INSPHERECALCULATOR_H_
 #define PUMGEN_SRC_HELPER_INSPHERECALCULATOR_H_
 
+#include <array>
 #include <cstddef>
 #include <mpi.h>
+#include <unordered_map>
 #include <vector>
 
-// assumes a contiguous distribution of all vertices over all processes
+/**
+ * The vertex coordinates of the local cells; the coordinates of vertices owned by other ranks are
+ * fetched once. Assumes a contiguous distribution of all vertices over all processes. Each cell
+ * has cellSize nodes, the first four of which are its vertices.
+ */
+class CellVertices {
+  public:
+  static constexpr std::size_t VerticesPerCell = 4;
 
-std::vector<double> calculateInsphere(const std::vector<std::size_t>& connectivity,
-                                      const std::vector<double>& geometry, MPI_Comm comm);
+  CellVertices(const std::vector<std::size_t>& connectivity, const std::vector<double>& geometry,
+               std::size_t cellSize, MPI_Comm comm);
+
+  [[nodiscard]] std::size_t numCells() const { return connectivity.size() / cellSize; }
+
+  [[nodiscard]] std::array<std::array<double, 3>, VerticesPerCell>
+  operator()(std::size_t cell) const;
+
+  private:
+  const std::vector<std::size_t>& connectivity;
+  const std::vector<double>& geometry;
+  std::size_t cellSize;
+  int commrank = 0;
+  std::vector<std::size_t> vertexDist;
+  std::vector<std::unordered_map<std::size_t, std::size_t>> outidxmap;
+  std::vector<std::size_t> outdisp;
+  std::vector<double> outvertices;
+};
+
+/**
+ * The insphere radius of each tetrahedron.
+ */
+std::vector<double> calculateInsphere(const CellVertices& cells);
 
 #endif // PUMGEN_SRC_HELPER_INSPHERECALCULATOR_H_