From e7eae8caae545de7b8dc83657031cd7176ff9bde Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 20 Apr 2026 22:01:42 -0500 Subject: [PATCH 1/2] Draft --- src/coreComponents/mesh/CMakeLists.txt | 6 + .../mesh/benchmarks/CMakeLists.txt | 13 + .../benchmarks/benchmarkScatterMethods.cpp | 467 ++++++++++ .../mesh/generators/VTKMeshGenerator.cpp | 30 +- .../mesh/generators/VTKMeshGenerator.hpp | 5 + .../mesh/generators/VTKMeshScattering.cpp | 828 ++++++++++++++++++ .../mesh/generators/VTKMeshScattering.hpp | 73 ++ .../mesh/generators/VTKUtilities.cpp | 158 +--- .../mesh/generators/VTKUtilities.hpp | 11 +- .../mesh/unitTests/CMakeLists.txt | 10 + .../mesh/unitTests/testVTKMeshScattering.cpp | 302 +++++++ 11 files changed, 1743 insertions(+), 160 deletions(-) create mode 100644 src/coreComponents/mesh/benchmarks/CMakeLists.txt create mode 100644 src/coreComponents/mesh/benchmarks/benchmarkScatterMethods.cpp create mode 100644 src/coreComponents/mesh/generators/VTKMeshScattering.cpp create mode 100644 src/coreComponents/mesh/generators/VTKMeshScattering.hpp create mode 100644 src/coreComponents/mesh/unitTests/testVTKMeshScattering.cpp diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 3b7e7a825cb..c1e960b503a 100644 --- a/src/coreComponents/mesh/CMakeLists.txt +++ b/src/coreComponents/mesh/CMakeLists.txt @@ -214,6 +214,7 @@ if( ENABLE_VTK ) generators/VTKHierarchicalDataSource.hpp generators/VTKMeshGenerator.hpp generators/VTKMeshGeneratorTools.hpp + generators/VTKMeshScattering.hpp generators/VTKWellGenerator.hpp generators/VTKUtilities.hpp ) @@ -223,6 +224,7 @@ if( ENABLE_VTK ) generators/VTKHierarchicalDataSource.cpp generators/VTKMeshGenerator.cpp generators/VTKMeshGeneratorTools.cpp + generators/VTKMeshScattering.cpp generators/VTKWellGenerator.cpp generators/VTKUtilities.cpp ) @@ -268,3 +270,7 @@ if( GEOS_ENABLE_TESTS ) add_subdirectory( unitTests ) endif( ) +# Add benchmarks subdirectory +if( ENABLE_VTK ) + add_subdirectory( benchmarks ) +endif() \ No newline at end of file diff --git a/src/coreComponents/mesh/benchmarks/CMakeLists.txt b/src/coreComponents/mesh/benchmarks/CMakeLists.txt new file mode 100644 index 00000000000..3c877949bcb --- /dev/null +++ b/src/coreComponents/mesh/benchmarks/CMakeLists.txt @@ -0,0 +1,13 @@ +# +# Mesh benchmarks +# + +blt_add_executable( + NAME benchmarkScatterMethods + SOURCES benchmarkScatterMethods.cpp + OUTPUT_DIR ${BENCHMARK_OUTPUT_DIRECTORY} + DEPENDS_ON mesh + FOLDER coreComponents/mesh/benchmarks +) + +install(TARGETS benchmarkScatterMethods DESTINATION ${CMAKE_INSTALL_PREFIX}/benchmarks) \ No newline at end of file diff --git a/src/coreComponents/mesh/benchmarks/benchmarkScatterMethods.cpp b/src/coreComponents/mesh/benchmarks/benchmarkScatterMethods.cpp new file mode 100644 index 00000000000..5b35aeed5bb --- /dev/null +++ b/src/coreComponents/mesh/benchmarks/benchmarkScatterMethods.cpp @@ -0,0 +1,467 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file benchmarkScatterMethods.cpp + * @brief Benchmark comparing VTK mesh scattering methods with load balance analysis + */ + +#include "common/DataTypes.hpp" +#include "common/format/Format.hpp" +#include "common/logger/Logger.hpp" +#include "common/MpiWrapper.hpp" +#include "mesh/generators/VTKMeshScattering.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace geos; +using namespace geos::vtk; + +namespace +{ + +// ============================================================================ +// Load balance and compactness metrics +// ============================================================================ + +struct LoadBalanceStats +{ + long long totalCells; + long long minCells; + long long maxCells; + double avgCells; + double imbalanceRatio; // max / avg (1.0 = perfect) + double bboxOverlap; // sum-of-local-bbox-volumes / total-bbox-volume (1.0 = ideal) +}; + +/** + * @brief Compute bounding-box overlap ratio as a proxy for spatial compactness. + * + * Each rank's local bounding box volume is summed across all ranks and divided + * by the total mesh bounding box volume. A ratio near 1.0 means partitions + * tile the domain with little overlap; a ratio near @p nranks means every + * partition spans the entire mesh (no spatial locality). + */ +double computeBboxOverlap( vtkDataSet * localMesh, + double const totalBounds[6], + MPI_Comm comm ) +{ + bool const hasLocal = localMesh->GetNumberOfCells() > 0; + + double localBounds[6] = { 0, 0, 0, 0, 0, 0 }; + if( hasLocal ) + { + localMesh->GetBounds( localBounds ); + } + + // Only use non-degenerate dimensions (extent > epsilon). + double totalVol = 1.0; + double localVol = hasLocal ? 1.0 : 0.0; + bool hasDim = false; + for( int d = 0; d < 3; ++d ) + { + double const totalExtent = totalBounds[2 * d + 1] - totalBounds[2 * d]; + if( totalExtent > 1.0e-12 ) + { + hasDim = true; + totalVol *= totalExtent; + if( hasLocal ) + { + localVol *= std::max( localBounds[2 * d + 1] - localBounds[2 * d], 0.0 ); + } + } + } + + if( !hasDim || totalVol <= 0.0 ) + { + return 0.0; + } + + double const sumLocalVol = MpiWrapper::allReduce( localVol, MpiWrapper::Reduction::Sum, comm ); + return sumLocalVol / totalVol; +} + +LoadBalanceStats computeLoadBalance( long long localCells, + vtkDataSet * localMesh, + double const totalBounds[6], + MPI_Comm comm ) +{ + LoadBalanceStats stats{}; + int const size = MpiWrapper::commSize( comm ); + + stats.totalCells = MpiWrapper::allReduce( localCells, MpiWrapper::Reduction::Sum, comm ); + stats.minCells = MpiWrapper::allReduce( localCells, MpiWrapper::Reduction::Min, comm ); + stats.maxCells = MpiWrapper::allReduce( localCells, MpiWrapper::Reduction::Max, comm ); + stats.avgCells = static_cast< double >( stats.totalCells ) / size; + stats.imbalanceRatio = ( stats.avgCells > 0 ) + ? stats.maxCells / stats.avgCells : 0.0; + stats.bboxOverlap = computeBboxOverlap( localMesh, totalBounds, comm ); + + return stats; +} + +// ============================================================================ +// Benchmark result and output +// ============================================================================ + +struct BenchmarkResult +{ + std::string name; + double time; + LoadBalanceStats balance; + bool cellsPreserved; +}; + +/** + * @brief Save a partitioned mesh to a per-method subdirectory. + * + * Each rank writes its own .vtu piece. Only rank 0 writes the .pvtu header. + * Directory layout: //_.vtu + * //.pvtu + */ +void savePartitionedMesh( vtkUnstructuredGrid * mesh, + std::string const & dir, + std::string const & methodName, + MPI_Comm comm ) +{ + int const rank = MpiWrapper::commRank( comm ); + + // Add RankID cell data + vtkNew< vtkIntArray > rankArray; + rankArray->SetName( "RankID" ); + rankArray->SetNumberOfTuples( mesh->GetNumberOfCells() ); + rankArray->FillValue( rank ); + mesh->GetCellData()->AddArray( rankArray ); + + std::string const methodDir = GEOS_FMT( "{}/{}", dir, methodName ); + if( rank == 0 ) + { + std::filesystem::create_directories( methodDir ); + } + MpiWrapper::barrier( comm ); + + // Write using the parallel XML writer + std::string const pvtuFile = GEOS_FMT( "{}/{}.pvtu", methodDir, methodName ); + + vtkNew< vtkMPIController > controller; + vtkMPICommunicatorOpaqueComm vtkComm( &comm ); + vtkNew< vtkMPICommunicator > communicator; + communicator->InitializeExternal( &vtkComm ); + controller->SetCommunicator( communicator ); + + vtkNew< vtkXMLPUnstructuredGridWriter > writer; + writer->SetController( controller ); + writer->SetFileName( pvtuFile.c_str() ); + writer->SetInputData( mesh ); + writer->SetNumberOfPieces( MpiWrapper::commSize( comm ) ); + writer->SetStartPiece( rank ); + writer->SetEndPiece( rank ); + writer->Write(); + + controller->SetCommunicator( nullptr ); + + MpiWrapper::barrier( comm ); + + if( rank == 0 ) + { + GEOS_LOG_RANK_0( GEOS_FMT( " Saved: {}", pvtuFile ) ); + } +} + +BenchmarkResult runBenchmark( std::string const & name, + ScatterMethod method, + vtkUnstructuredGrid & mesh, + arrayView1d< int const > const & partitions, + long long totalCells, + double const totalBounds[6], + std::string const & outputDir, + MPI_Comm comm ) +{ + BenchmarkResult result; + result.name = name; + + MpiWrapper::barrier( comm ); + + auto start = std::chrono::high_resolution_clock::now(); + + vtkSmartPointer< vtkDataSet > scattered; + if( method == ScatterMethod::kdtree ) + { + // Call VTK redistribution directly to get raw kdtree metrics, + // bypassing the contiguous fallback in scatterMesh. + int const size = MpiWrapper::commSize( comm ); + vtkNew< vtkMPIController > controller; + vtkMPICommunicatorOpaqueComm vtkComm( &comm ); + vtkNew< vtkMPICommunicator > communicator; + communicator->InitializeExternal( &vtkComm ); + controller->SetCommunicator( communicator ); + vtkMultiProcessController::SetGlobalController( controller ); + + vtkNew< vtkRedistributeDataSetFilter > rdsf; + rdsf->SetInputDataObject( &mesh ); + rdsf->SetNumberOfPartitions( size ); + rdsf->Update(); + scattered = vtkDataSet::SafeDownCast( rdsf->GetOutputDataObject( 0 ) ); + } + else + { + scattered = scatterMesh( method, mesh, partitions, comm ); + } + vtkSmartPointer< vtkUnstructuredGrid > output = + vtkUnstructuredGrid::SafeDownCast( scattered ); + + auto end = std::chrono::high_resolution_clock::now(); + double localTime = std::chrono::duration< double >( end - start ).count(); + + long long localCells = output->GetNumberOfCells(); + + result.time = MpiWrapper::reduce( localTime, MpiWrapper::Reduction::Max, 0, comm ); + result.balance = computeLoadBalance( localCells, output, totalBounds, comm ); + result.cellsPreserved = ( result.balance.totalCells == totalCells ); + + GEOS_LOG_RANK_0( GEOS_FMT( "\n=== {} ===\n" + " Time: {:.3f}s\n" + " Cells: {}{}\n" + " Balance: min={} max={} avg={:.1f}" + " imbalance={:.6f}x bbox_overlap={:.3f}", + name, + result.time, + result.balance.totalCells, + result.cellsPreserved ? "" + : GEOS_FMT( " WARNING: lost {} cells!", + totalCells - result.balance.totalCells ), + result.balance.minCells, + result.balance.maxCells, + result.balance.avgCells, + result.balance.imbalanceRatio, + result.balance.bboxOverlap ) ); + + // Save partitioned mesh immediately if an output directory was provided. + if( !outputDir.empty() ) + { + savePartitionedMesh( output, outputDir, name, comm ); + } + + MpiWrapper::barrier( comm ); + return result; +} + +void printSummary( std::vector< BenchmarkResult > const & results, + long long totalCells, + MPI_Comm comm ) +{ + int const size = MpiWrapper::commSize( comm ); + + std::string const sep( 110, '=' ); + std::string const thin( 110, '-' ); + + GEOS_LOG_RANK_0( GEOS_FMT( "\n{}\nSUMMARY ({} ranks, {} cells)\n{}", + sep, size, totalCells, sep ) ); + + // Table header + GEOS_LOG_RANK_0( GEOS_FMT( "{:<14}| {:<10}| {:<12}| {:<8}| {:<10}| {:<10}| {:<12}| {:<12}", + "Method", "Time (s)", "Cells", "Status", + "Min cells", "Max cells", "Imbalance", "BBox Overlap" ) ); + GEOS_LOG_RANK_0( thin ); + + // Table rows + for( auto const & r : results ) + { + GEOS_LOG_RANK_0( GEOS_FMT( "{:<14}| {:<10.3f}| {:<12}| {:<8}| {:<10}| {:<10}| {:<12.6f}| {:<12.3f}", + r.name, r.time, + r.balance.totalCells, + r.cellsPreserved ? "OK" : "LOSS", + r.balance.minCells, + r.balance.maxCells, + r.balance.imbalanceRatio, + r.balance.bboxOverlap ) ); + } + GEOS_LOG_RANK_0( thin ); + + // Find fastest method + double fastestTime = results[0].time; + for( auto const & r : results ) + { + fastestTime = std::min( fastestTime, r.time ); + } + + // Speedup table + GEOS_LOG_RANK_0( "\nSpeedup relative to each method:" ); + for( auto const & r : results ) + { + double const speedup = r.time / fastestTime; + GEOS_LOG_RANK_0( GEOS_FMT( " {:<14}{:.2f}x slower than fastest{}", + r.name, speedup, + ( std::abs( r.time - fastestTime ) < 1e-6 ) ? " <-- fastest" : "" ) ); + } + + // Best balance and compactness + double bestImbalance = results[0].balance.imbalanceRatio; + double bestOverlap = results[0].balance.bboxOverlap; + std::string bestBalanceName = results[0].name; + std::string bestCompactName = results[0].name; + for( auto const & r : results ) + { + if( r.cellsPreserved ) + { + if( r.balance.imbalanceRatio < bestImbalance ) + { + bestImbalance = r.balance.imbalanceRatio; bestBalanceName = r.name; + } + if( r.balance.bboxOverlap < bestOverlap ) + { + bestOverlap = r.balance.bboxOverlap; bestCompactName = r.name; + } + } + } + GEOS_LOG_RANK_0( GEOS_FMT( "\nBest load balance: {} (imbalance ratio: {:.6f}x)", + bestBalanceName, bestImbalance ) ); + GEOS_LOG_RANK_0( GEOS_FMT( "Best compactness: {} (bbox overlap: {:.3f})", bestCompactName, bestOverlap ) ); + GEOS_LOG_RANK_0( sep ); +} + +} // anonymous namespace + +int main( int argc, char * argv[] ) +{ + MpiWrapper::init( &argc, &argv ); + +#ifdef GEOS_USE_CHAI + chai::ArrayManager::getInstance()->disableCallbacks(); +#endif + + MPI_Comm const comm = MPI_COMM_WORLD; + int const size = MpiWrapper::commSize( comm ); + + logger::InitializeLogger( comm ); + + if( argc < 5 ) + { + GEOS_LOG_RANK_0( GEOS_FMT( "Usage: {} [output_dir]\n" + " nx*ny*nz must equal number of MPI ranks (for Cartesian method)\n" + " output_dir (optional): save partitioned meshes (one subfolder per method)", + argv[0] ) ); + logger::FinalizeLogger(); + MpiWrapper::finalize(); + return 1; + } + + int const nx = std::atoi( argv[2] ); + int const ny = std::atoi( argv[3] ); + int const nz = std::atoi( argv[4] ); + + std::string outputDir; + if( argc > 5 ) + { + outputDir = argv[5]; + } + + GEOS_LOG_RANK_0( "=== Redistribution Performance Comparison ===" ); + GEOS_LOG_RANK_0( GEOS_FMT( "MPI Ranks: {}", size ) ); + GEOS_LOG_RANK_0( GEOS_FMT( "Grid: {} x {} x {}", nx, ny, nz ) ); + GEOS_LOG_RANK_0( GEOS_FMT( "Input file: {}", argv[1] ) ); + + // Load mesh on rank 0 + vtkSmartPointer< vtkUnstructuredGrid > mesh; + long long totalCells = 0; + double totalBounds[6] = { 0, 0, 0, 0, 0, 0 }; + + int const rank = MpiWrapper::commRank( comm ); + if( rank == 0 ) + { + vtkNew< vtkXMLUnstructuredGridReader > reader; + reader->SetFileName( argv[1] ); + reader->Update(); + mesh = reader->GetOutput(); + + if( !mesh || mesh->GetNumberOfCells() == 0 ) + { + std::cerr << "Failed to read mesh or mesh is empty\n"; + MPI_Abort( comm, 1 ); + } + + totalCells = mesh->GetNumberOfCells(); + mesh->GetBounds( totalBounds ); + } + else + { + mesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + } + + MpiWrapper::bcast( &totalCells, 1, 0, comm ); + MpiWrapper::bcast( totalBounds, 6, 0, comm ); + + GEOS_LOG_RANK_0( GEOS_FMT( "Loaded {} cells", totalCells ) ); + GEOS_LOG_RANK_0( GEOS_FMT( "Bounds: [{}, {}] x [{}, {}] x [{}, {}]", + totalBounds[0], totalBounds[1], + totalBounds[2], totalBounds[3], + totalBounds[4], totalBounds[5] ) ); + + std::vector< BenchmarkResult > results; + + array1d< int > parts( 3 ); + parts[0] = nx; parts[1] = ny; parts[2] = nz; + + auto run = [&]( std::string const & name, ScatterMethod method ) + { + GEOS_LOG_RANK_0( GEOS_FMT( "\nRunning {} scatter ({} cells, {} ranks)...", name, totalCells, size ) ); + results.push_back( runBenchmark( name, method, *mesh, parts.toViewConst(), + totalCells, totalBounds, outputDir, comm ) ); + }; + + run( "KdTree", ScatterMethod::kdtree ); + run( "Contiguous", ScatterMethod::contiguous ); + + if( nx * ny * nz == size ) + { + run( "Cartesian", ScatterMethod::cartesian ); + } + else + { + GEOS_LOG_RANK_0( GEOS_FMT( "\n=== Cartesian ===\n SKIPPED: nx*ny*nz ({}) != MPI size ({})", + nx * ny * nz, size ) ); + } + + run( "RCB", ScatterMethod::rcb ); + + // Print summary + printSummary( results, totalCells, comm ); + + logger::FinalizeLogger(); + MpiWrapper::finalize(); + return 0; +} diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp index b90bd1651e7..6d3a280921b 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.cpp @@ -80,7 +80,16 @@ VTKMeshGenerator::VTKMeshGenerator( string const & name, registerWrapper( viewKeyStruct::partitionMethodString(), &m_partitionMethod ). setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Method (library) used to partition the mesh" ); + setDescription( "Method (library) used to refine mesh partitioning" ); + + registerWrapper( viewKeyStruct::scatterMethodString(), &m_scatterMethod ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( vtk::ScatterMethod::kdtree ). + setDescription( "Method for initial mesh scatter from rank 0 to all ranks: " + "contiguous (cell ID ranges, no geometry), " + "cartesian (regular grid using -x/-y/-z partitions), " + "kdtree (VTK built-in kd-tree, default), " + "rcb (recursive coordinate bisection)" ); registerWrapper( viewKeyStruct::useGlobalIdsString(), &m_useGlobalIds ). setInputFlag( InputFlags::OPTIONAL ). @@ -136,7 +145,17 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager vtkSmartPointer< vtkMultiProcessController > controller = vtk::getController(); vtkMultiProcessController::SetGlobalController( controller ); - int const numPartZ = m_structuredIndexAttributeName.empty() ? 1 : partition.getPartitions()[2]; + array1d< int > const & partitions = partition.getPartitions(); + + if( m_scatterMethod == vtk::ScatterMethod::cartesian ) + { + int const product = partitions[0] * partitions[1] * partitions[2]; + GEOS_ERROR_IF( product != MpiWrapper::commSize( comm ), + GEOS_FMT( "scatterMethod=\"cartesian\" requires -x * -y * -z = MPI size. " + "Got {}x{}x{} = {} but MPI size is {}.", + partitions[0], partitions[1], partitions[2], + product, MpiWrapper::commSize( comm ) ) ); + } GEOS_LOG_LEVEL_RANK_0( logInfo::VTKSteps, " redistributing mesh..." ); { @@ -153,7 +172,7 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager { if( MpiWrapper::commRank() == 0 ) { - stdVector< vtkSmartPointer< vtkPartitionedDataSet > > partitions; + stdVector< vtkSmartPointer< vtkPartitionedDataSet > > vtkPartitions; vtkNew< vtkAppendFilter > appender; appender->MergePointsOn(); for( auto & [key, value] : getSubGroups()) @@ -206,11 +225,12 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager allMeshes.getMainMesh(), allMeshes.getFaceBlocks(), comm, + m_scatterMethod, + partitions.toViewConst(), m_partitionMethod, m_partitionRefinement, m_useGlobalIds, - m_structuredIndexAttributeName, - numPartZ ); + m_structuredIndexAttributeName ); m_vtkMesh = redistributedMeshes.getMainMesh(); m_faceBlockMeshes = redistributedMeshes.getFaceBlocks(); GEOS_LOG_LEVEL_RANK_0( logInfo::VTKSteps, GEOS_FMT( "{} '{}': finding neighbor ranks...", catalogName(), getName() ) ); diff --git a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp index cd968c2abd4..e5e52736f49 100644 --- a/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp +++ b/src/coreComponents/mesh/generators/VTKMeshGenerator.hpp @@ -21,6 +21,7 @@ #define GEOS_MESH_GENERATORS_VTKMESHGENERATOR_HPP #include "mesh/generators/ExternalMeshGeneratorBase.hpp" +#include "mesh/generators/VTKMeshScattering.hpp" #include "mesh/generators/VTKUtilities.hpp" #include "mesh/generators/VTKHierarchicalDataSource.hpp" #include "mesh/mpiCommunications/SpatialPartition.hpp" @@ -114,6 +115,7 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase constexpr static char const * nodesetNamesString() { return "nodesetNames"; } constexpr static char const * partitionRefinementString() { return "partitionRefinement"; } constexpr static char const * partitionMethodString() { return "partitionMethod"; } + constexpr static char const * scatterMethodString() { return "scatterMethod"; } constexpr static char const * useGlobalIdsString() { return "useGlobalIds"; } constexpr static char const * dataSourceString() { return "dataSourceName"; } }; @@ -167,6 +169,9 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase /// Method (library) used to partition the mesh vtk::PartitionMethod m_partitionMethod = vtk::PartitionMethod::parmetis; + /// Method used for initial mesh scatter from rank 0 + vtk::ScatterMethod m_scatterMethod = vtk::ScatterMethod::kdtree; + /// Lists of VTK cell ids, organized by element type, then by region vtk::CellMapType m_cellMap; diff --git a/src/coreComponents/mesh/generators/VTKMeshScattering.cpp b/src/coreComponents/mesh/generators/VTKMeshScattering.cpp new file mode 100644 index 00000000000..e5880a0d146 --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKMeshScattering.cpp @@ -0,0 +1,828 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "mesh/generators/VTKMeshScattering.hpp" + +#include "common/format/Format.hpp" +#include "common/logger/Logger.hpp" +#include "common/MpiWrapper.hpp" +#include "common/TimingMacros.hpp" +#include "LvArray/src/math.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef GEOS_USE_MPI +#include +#include +#else +#include +#endif +#include +#include + +#include +#include +#include +#include +#include + +namespace geos +{ +namespace vtk +{ + +namespace +{ + +// ============================================================================ +// Buffer serialization helpers +// ============================================================================ + +void appendBytes( stdVector< char > & buf, void const * data, int64_t n ) +{ + int64_t const pos = static_cast< int64_t >( buf.size() ); + buf.resize( static_cast< size_t >( pos + n ) ); + std::memcpy( buf.data() + pos, data, static_cast< size_t >( n ) ); +} + +template< typename T > +void appendValue( stdVector< char > & buf, T val ) +{ + appendBytes( buf, &val, sizeof( T ) ); +} + +template< typename T > +T readValue( char const * & ptr ) +{ + T val; + std::memcpy( &val, ptr, sizeof( T ) ); + ptr += sizeof( T ); + return val; +} + +// ============================================================================ +// MPI large-message helpers (handles buffers > 2 GB) +// ============================================================================ + +constexpr int64_t MPI_CHUNK_SIZE = INT64_C( 1 ) << 30; // 1 GB + +void mpiSendLarge( void const * buf, int64_t count, integer dest, integer tag, MPI_Comm comm ) +{ + char const * ptr = static_cast< char const * >( buf ); + while( count > 0 ) + { + integer const chunk = static_cast< integer >( LvArray::math::min( count, MPI_CHUNK_SIZE ) ); + MpiWrapper::send( ptr, chunk, dest, tag, comm ); + ptr += chunk; + count -= chunk; + } +} + +void mpiRecvLarge( void * buf, int64_t count, integer src, integer tag, MPI_Comm comm ) +{ + char * ptr = static_cast< char * >( buf ); + while( count > 0 ) + { + integer const chunk = static_cast< integer >( LvArray::math::min( count, MPI_CHUNK_SIZE ) ); + MPI_Recv( ptr, chunk, MPI_BYTE, src, tag, comm, MPI_STATUS_IGNORE ); + ptr += chunk; + count -= chunk; + } +} + +// ============================================================================ +// Pack / unpack vtkDataSetAttributes (cell data or point data) +// ============================================================================ + +constexpr int NUM_ATTR_TYPES = vtkDataSetAttributes::NUM_ATTRIBUTES; + +void packDataArrays( stdVector< char > & buf, vtkDataSetAttributes * attrs ) +{ + // Count only vtkDataArray entries (skip string / abstract arrays) + int32_t nArrays = 0; + for( int a = 0; a < attrs->GetNumberOfArrays(); ++a ) + { + if( attrs->GetArray( a ) != nullptr ) + { + ++nArrays; + } + } + appendValue( buf, nArrays ); + + for( int a = 0; a < attrs->GetNumberOfArrays(); ++a ) + { + vtkDataArray * arr = attrs->GetArray( a ); + if( arr == nullptr ) + { + continue; + } + char const * name = arr->GetName() ? arr->GetName() : ""; + int32_t const nameLen = static_cast< int32_t >( std::strlen( name ) ); + appendValue( buf, nameLen ); + appendBytes( buf, name, nameLen ); + + int32_t const nComp = arr->GetNumberOfComponents(); + int32_t const dataType = arr->GetDataType(); + int64_t const nTuples = arr->GetNumberOfTuples(); + int32_t const elemSize = arr->GetDataTypeSize(); + + appendValue( buf, nComp ); + appendValue( buf, dataType ); + appendValue( buf, nTuples ); + appendValue( buf, elemSize ); + + int64_t const dataBytes = nTuples * static_cast< int64_t >( nComp ) * static_cast< int64_t >( elemSize ); + appendBytes( buf, arr->GetVoidPointer( 0 ), dataBytes ); + } + + // Serialize active attribute designations (SCALARS=0 .. GLOBALIDS=5). + // For each slot store the name length + name (empty string = none). + for( int t = 0; t < NUM_ATTR_TYPES; ++t ) + { + vtkDataArray * active = attrs->GetAttribute( t ); + char const * activeName = ( active && active->GetName() ) ? active->GetName() : ""; + int32_t const nameLen = static_cast< int32_t >( std::strlen( activeName ) ); + appendValue( buf, nameLen ); + if( nameLen > 0 ) + { + appendBytes( buf, activeName, nameLen ); + } + } +} + +void unpackDataArrays( char const * & ptr, vtkDataSetAttributes * attrs ) +{ + int32_t const nArrays = readValue< int32_t >( ptr ); + + for( int a = 0; a < nArrays; ++a ) + { + int32_t const nameLen = readValue< int32_t >( ptr ); + string name( ptr, nameLen ); + ptr += nameLen; + + int32_t const nComp = readValue< int32_t >( ptr ); + int32_t const dataType = readValue< int32_t >( ptr ); + int64_t const nTuples = readValue< int64_t >( ptr ); + int32_t const elemSize = readValue< int32_t >( ptr ); + + vtkSmartPointer< vtkDataArray > arr( vtkDataArray::CreateDataArray( dataType ) ); + arr->SetName( name.c_str() ); + arr->SetNumberOfComponents( nComp ); + arr->SetNumberOfTuples( nTuples ); + + int64_t const dataBytes = nTuples * static_cast< int64_t >( nComp ) * static_cast< int64_t >( elemSize ); + std::memcpy( arr->GetVoidPointer( 0 ), ptr, static_cast< size_t >( dataBytes ) ); + ptr += dataBytes; + + attrs->AddArray( arr ); + } + + // Restore active attribute designations + for( int t = 0; t < NUM_ATTR_TYPES; ++t ) + { + int32_t const nameLen = readValue< int32_t >( ptr ); + if( nameLen > 0 ) + { + string activeName( ptr, nameLen ); + ptr += nameLen; + attrs->SetActiveAttribute( activeName.c_str(), t ); + } + } +} + +// ============================================================================ +// Pack / unpack a full vtkUnstructuredGrid + assignment vector +// ============================================================================ + +void packGrid( vtkUnstructuredGrid * grid, + stdVector< integer > const & assignment, + stdVector< char > & buf ) +{ + buf.clear(); + + int64_t const nPoints = grid->GetNumberOfPoints(); + int64_t const nCells = grid->GetNumberOfCells(); + appendValue( buf, nPoints ); + appendValue( buf, nCells ); + + // Points (always serialized as real64) + if( nPoints > 0 ) + { + vtkPoints * points = grid->GetPoints(); + if( points->GetDataType() == VTK_DOUBLE ) + { + appendBytes( buf, points->GetVoidPointer( 0 ), nPoints * 3 * sizeof( real64 ) ); + } + else + { + for( int64_t i = 0; i < nPoints; ++i ) + { + real64 p[3]; + points->GetPoint( i, p ); + appendBytes( buf, p, sizeof( p ) ); + } + } + } + + // Cell types, offsets, connectivity + if( nCells > 0 ) + { + vtkUnsignedCharArray * types = grid->GetCellTypesArray(); + appendBytes( buf, types->GetVoidPointer( 0 ), nCells * sizeof( unsigned char ) ); + + vtkCellArray * cells = grid->GetCells(); + vtkDataArray * offsets = cells->GetOffsetsArray(); + vtkDataArray * conn = cells->GetConnectivityArray(); + + int64_t const nOffsets = offsets->GetNumberOfValues(); + int64_t const connSize = conn->GetNumberOfValues(); + + appendValue( buf, connSize ); + appendBytes( buf, offsets->GetVoidPointer( 0 ), nOffsets * sizeof( vtkIdType ) ); + appendBytes( buf, conn->GetVoidPointer( 0 ), connSize * sizeof( vtkIdType ) ); + } + + // Field data arrays + packDataArrays( buf, grid->GetCellData() ); + packDataArrays( buf, grid->GetPointData() ); + + // Assignment vector + int64_t const assignSize = static_cast< int64_t >( assignment.size() ); + appendValue( buf, assignSize ); + if( assignSize > 0 ) + { + appendBytes( buf, assignment.data(), assignSize * sizeof( integer ) ); + } +} + + +std::pair< vtkSmartPointer< vtkUnstructuredGrid >, stdVector< integer > > +unpackGrid( stdVector< char > const & buf ) +{ + char const * ptr = buf.data(); + + int64_t const nPoints = readValue< int64_t >( ptr ); + int64_t const nCells = readValue< int64_t >( ptr ); + + auto grid = vtkSmartPointer< vtkUnstructuredGrid >::New(); + + // Points + if( nPoints > 0 ) + { + vtkNew< vtkPoints > points; + points->SetDataTypeToDouble(); + points->SetNumberOfPoints( nPoints ); + std::memcpy( points->GetVoidPointer( 0 ), ptr, nPoints * 3 * sizeof( real64 ) ); + ptr += nPoints * 3 * sizeof( real64 ); + grid->SetPoints( points ); + } + + // Cells + if( nCells > 0 ) + { + vtkNew< vtkUnsignedCharArray > types; + types->SetNumberOfValues( nCells ); + std::memcpy( types->GetVoidPointer( 0 ), ptr, nCells * sizeof( unsigned char ) ); + ptr += nCells * sizeof( unsigned char ); + + int64_t const connSize = readValue< int64_t >( ptr ); + + vtkNew< vtkIdTypeArray > offsets; + offsets->SetNumberOfValues( nCells + 1 ); + std::memcpy( offsets->GetVoidPointer( 0 ), ptr, ( nCells + 1 ) * sizeof( vtkIdType ) ); + ptr += ( nCells + 1 ) * sizeof( vtkIdType ); + + vtkNew< vtkIdTypeArray > conn; + conn->SetNumberOfValues( connSize ); + std::memcpy( conn->GetVoidPointer( 0 ), ptr, connSize * sizeof( vtkIdType ) ); + ptr += connSize * sizeof( vtkIdType ); + + vtkNew< vtkCellArray > cellArray; + cellArray->SetData( offsets, conn ); + grid->SetCells( types, cellArray ); + } + + // Field data + unpackDataArrays( ptr, grid->GetCellData() ); + unpackDataArrays( ptr, grid->GetPointData() ); + + // Assignment + int64_t const assignSize = readValue< int64_t >( ptr ); + stdVector< integer > assignment( assignSize ); + if( assignSize > 0 ) + { + std::memcpy( assignment.data(), ptr, assignSize * sizeof( integer ) ); + ptr += assignSize * sizeof( integer ); + } + + return { grid, std::move( assignment ) }; +} + + +// ============================================================================ +// Split cells into two subsets: those with assignment < mid (kept locally) +// and those with assignment >= mid (sent to the partner rank). +// ============================================================================ + +struct SplitResult +{ + vtkSmartPointer< vtkUnstructuredGrid > loMesh; + stdVector< integer > loAssignment; + vtkSmartPointer< vtkUnstructuredGrid > hiMesh; + stdVector< integer > hiAssignment; +}; + +SplitResult splitByMid( vtkUnstructuredGrid * mesh, + stdVector< integer > const & assignment, + integer mid ) +{ + stdVector< vtkIdType > loCells, hiCells; + stdVector< integer > loAssign, hiAssign; + loCells.reserve( assignment.size() / 2 ); + hiCells.reserve( assignment.size() / 2 ); + loAssign.reserve( assignment.size() / 2 ); + hiAssign.reserve( assignment.size() / 2 ); + + for( vtkIdType i = 0; i < static_cast< vtkIdType >( assignment.size() ); ++i ) + { + if( assignment[i] < mid ) + { + loCells.push_back( i ); + loAssign.push_back( assignment[i] ); + } + else + { + hiCells.push_back( i ); + hiAssign.push_back( assignment[i] ); + } + } + + auto extract = []( vtkUnstructuredGrid * inputMesh, + stdVector< vtkIdType > const & ids ) + -> vtkSmartPointer< vtkUnstructuredGrid > + { + if( ids.empty() ) + { + return vtkSmartPointer< vtkUnstructuredGrid >::New(); + } + vtkNew< vtkExtractCells > extractor; + extractor->SetInputData( inputMesh ); + extractor->SetCellIds( ids.data(), static_cast< vtkIdType >( ids.size() ) ); + extractor->Update(); + auto result = vtkSmartPointer< vtkUnstructuredGrid >::New(); + result->ShallowCopy( extractor->GetOutput() ); + return result; + }; + + // Extract the smaller subset first, then release the original mesh reference + // before extracting the larger one to reduce peak memory. + SplitResult result; + if( hiCells.size() <= loCells.size() ) + { + result.hiMesh = extract( mesh, hiCells ); + result.hiAssignment = std::move( hiAssign ); + result.loMesh = extract( mesh, loCells ); + result.loAssignment = std::move( loAssign ); + } + else + { + result.loMesh = extract( mesh, loCells ); + result.loAssignment = std::move( loAssign ); + result.hiMesh = extract( mesh, hiCells ); + result.hiAssignment = std::move( hiAssign ); + } + return result; +} + + +// ============================================================================ +// Centroid computation +// ============================================================================ + +stdVector< stdArray< real64, 3 > > +computeCentroids( vtkDataSet & mesh ) +{ + vtkIdType const n = mesh.GetNumberOfCells(); + stdVector< stdArray< real64, 3 > > centroids( n ); + vtkNew< vtkIdList > ptIds; + + for( vtkIdType c = 0; c < n; ++c ) + { + mesh.GetCellPoints( c, ptIds ); + real64 cx = 0.0, cy = 0.0, cz = 0.0; + vtkIdType const nPts = ptIds->GetNumberOfIds(); + for( vtkIdType i = 0; i < nPts; ++i ) + { + real64 p[3]; + mesh.GetPoint( ptIds->GetId( i ), p ); + cx += p[0]; + cy += p[1]; + cz += p[2]; + } + if( nPts > 0 ) + { + real64 const inv = 1.0 / nPts; + centroids[c] = { cx * inv, cy * inv, cz * inv }; + } + else + { + centroids[c] = { 0.0, 0.0, 0.0 }; + } + } + return centroids; +} + + +// ============================================================================ +// Cell-rank assignment: contiguous (index-based, no geometry) +// ============================================================================ + +stdVector< integer > +computeCellRanksContiguous( vtkIdType nCells, integer size ) +{ + stdVector< integer > ranks( nCells ); + vtkIdType const perRank = nCells / size; + vtkIdType const remainder = nCells % size; + + // Ranks [0, remainder) get (perRank+1) cells, the rest get perRank + vtkIdType cell = 0; + for( integer r = 0; r < size; ++r ) + { + vtkIdType const count = perRank + ( r < remainder ? 1 : 0 ); + for( vtkIdType i = 0; i < count; ++i ) + { + ranks[cell++] = r; + } + } + return ranks; +} + + +// ============================================================================ +// Cell-rank assignment: Cartesian grid partition +// ============================================================================ + +stdVector< integer > +computeCellRanksCartesian( vtkDataSet & mesh, integer nx, integer ny, integer nz ) +{ + vtkIdType const numCells = mesh.GetNumberOfCells(); + + real64 bounds[6]; + mesh.GetBounds( bounds ); + real64 const xMin = bounds[0], xMax = bounds[1]; + real64 const yMin = bounds[2], yMax = bounds[3]; + real64 const zMin = bounds[4], zMax = bounds[5]; + + GEOS_ERROR_IF( nx > 1 && xMax <= xMin, + GEOS_FMT( "computeCellRanksCartesian: nx={} but mesh has zero extent in x ([{}, {}])", nx, xMin, xMax ) ); + GEOS_ERROR_IF( ny > 1 && yMax <= yMin, + GEOS_FMT( "computeCellRanksCartesian: ny={} but mesh has zero extent in y ([{}, {}])", ny, yMin, yMax ) ); + GEOS_ERROR_IF( nz > 1 && zMax <= zMin, + GEOS_FMT( "computeCellRanksCartesian: nz={} but mesh has zero extent in z ([{}, {}])", nz, zMin, zMax ) ); + + real64 const dx = nx > 1 ? ( xMax - xMin ) / nx : 1.0; + real64 const dy = ny > 1 ? ( yMax - yMin ) / ny : 1.0; + real64 const dz = nz > 1 ? ( zMax - zMin ) / nz : 1.0; + + auto centroids = computeCentroids( mesh ); + stdVector< integer > ranks( numCells ); + + for( vtkIdType c = 0; c < numCells; ++c ) + { + integer const ix = std::clamp( static_cast< integer >( ( centroids[c][0] - xMin ) / dx ), 0, nx - 1 ); + integer const iy = std::clamp( static_cast< integer >( ( centroids[c][1] - yMin ) / dy ), 0, ny - 1 ); + integer const iz = std::clamp( static_cast< integer >( ( centroids[c][2] - zMin ) / dz ), 0, nz - 1 ); + // Rank ordering: ix + nx*(iy + ny*iz) (X-fastest, Z-slowest). + // This matches the convention GEOS uses for -x/-y/-z grid decomposition. + ranks[c] = ix + nx * ( iy + ny * iz ); + } + return ranks; +} + + +// ============================================================================ +// Cell-rank assignment: Recursive Coordinate Bisection (nth_element) +// ============================================================================ + +stdVector< integer > +computeCellRanksRCB( vtkDataSet & mesh, integer size ) +{ + auto centroids = computeCentroids( mesh ); + vtkIdType const n = mesh.GetNumberOfCells(); + + stdVector< vtkIdType > indices( n ); + std::iota( indices.begin(), indices.end(), 0 ); + + stdVector< integer > ranks( n ); + + constexpr real64 realMax = std::numeric_limits< real64 >::max(); + + std::function< void( vtkIdType, vtkIdType, integer, integer ) > bisect; + bisect = [&]( vtkIdType begin, vtkIdType end, integer rankLo, integer rankHi ) + { + if( rankHi - rankLo == 1 ) + { + for( vtkIdType i = begin; i < end; ++i ) + { + ranks[indices[i]] = rankLo; + } + return; + } + + // Find the dimension with the largest spread + real64 lo[3] = { realMax, realMax, realMax }; + real64 hi[3] = { -realMax, -realMax, -realMax }; + for( vtkIdType i = begin; i < end; ++i ) + { + auto const & c = centroids[indices[i]]; + for( integer d = 0; d < 3; ++d ) + { + lo[d] = LvArray::math::min( lo[d], c[d] ); + hi[d] = LvArray::math::max( hi[d], c[d] ); + } + } + integer bestDim = 0; + if( hi[1] - lo[1] > hi[bestDim] - lo[bestDim] ) + bestDim = 1; + if( hi[2] - lo[2] > hi[bestDim] - lo[bestDim] ) + bestDim = 2; + + // Split proportionally to balance cell counts between left/right rank groups + integer const leftParts = ( rankHi - rankLo ) / 2; + integer const totalParts = rankHi - rankLo; + vtkIdType const splitAt = begin + ( ( end - begin ) * leftParts ) / totalParts; + integer const rankMid = rankLo + leftParts; + + std::nth_element( indices.begin() + begin, + indices.begin() + splitAt, + indices.begin() + end, + [&]( vtkIdType a, vtkIdType b ) + { return centroids[a][bestDim] < centroids[b][bestDim]; } ); + + bisect( begin, splitAt, rankLo, rankMid ); + bisect( splitAt, end, rankMid, rankHi ); + }; + + bisect( 0, n, 0, size ); + return ranks; +} + + + +// ============================================================================ +// Binary-tree scatter +// +// Only rank 0 holds mesh data and the assignment vector initially. +// At each level, the "root" of each active subrange extracts the upper +// half, packs it into a raw buffer, sends it to the midpoint rank, and +// keeps only the lower half. +// ============================================================================ + +vtkSmartPointer< vtkUnstructuredGrid > +binaryTreeScatter( vtkUnstructuredGrid * inputMesh, + stdVector< integer > assignment, + MPI_Comm comm ) +{ + integer const rank = MpiWrapper::commRank( comm ); + integer const size = MpiWrapper::commSize( comm ); + + // Working mesh: only rank 0 starts with data + vtkSmartPointer< vtkUnstructuredGrid > workingMesh; + if( rank == 0 ) + { + workingMesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + workingMesh->ShallowCopy( inputMesh ); + } + + integer lo = 0; + integer hi = size; + + while( hi - lo > 1 ) + { + integer const mid = lo + ( hi - lo ) / 2; + + if( rank == lo ) + { + // Sender: always send bufSize to mid (even if mesh is empty) to avoid receiver deadlock. + stdVector< char > buffer; + if( workingMesh && workingMesh->GetNumberOfCells() > 0 ) + { + // Split into lower (keep) and upper (send) halves. + auto split = splitByMid( workingMesh, assignment, mid ); + workingMesh = nullptr; // release original before packing + + packGrid( split.hiMesh, split.hiAssignment, buffer ); + split.hiMesh = nullptr; + split.hiAssignment.clear(); + + // Keep the lower half. + workingMesh = std::move( split.loMesh ); + assignment = std::move( split.loAssignment ); + } + + int64_t bufSize = static_cast< int64_t >( buffer.size() ); + MpiWrapper::send( &bufSize, 1, mid, 0, comm ); + if( bufSize > 0 ) + { + mpiSendLarge( buffer.data(), bufSize, mid, 1, comm ); + } + + hi = mid; + } + else if( rank == mid ) + { + // Receiver: receive from rank lo + + int64_t bufSize = 0; + MPI_Recv( &bufSize, 1, MPI_INT64_T, lo, 0, comm, MPI_STATUS_IGNORE ); + + if( bufSize > 0 ) + { + stdVector< char > buffer( bufSize ); + mpiRecvLarge( buffer.data(), bufSize, lo, 1, comm ); + + auto [mesh, assign] = unpackGrid( buffer ); + workingMesh = mesh; + assignment = std::move( assign ); + } + else + { + workingMesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + assignment.clear(); + } + + lo = mid; + } + else + { + // Inactive at this level, just narrow the range + if( rank < mid ) + { + hi = mid; + } + else + { + lo = mid; + } + } + } + + if( !workingMesh ) + { + workingMesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + } + return workingMesh; +} + +} // anonymous namespace + + +// ============================================================================ +// Public API +// ============================================================================ + +vtkSmartPointer< vtkDataSet > +scatterMesh( ScatterMethod method, + vtkDataSet & mesh, + arrayView1d< integer const > cartesianPartitions, + MPI_Comm comm ) +{ + GEOS_MARK_FUNCTION; + + integer const rank = MpiWrapper::commRank( comm ); + integer const size = MpiWrapper::commSize( comm ); + + // Early exits + vtkIdType const localCells = mesh.GetNumberOfCells(); + vtkIdType const totalCells = MpiWrapper::allReduce( localCells, MpiWrapper::Reduction::Sum, comm ); + + if( totalCells == 0 ) + { + return vtkSmartPointer< vtkUnstructuredGrid >::New(); + } + if( size == 1 ) + { + auto copy = vtkSmartPointer< vtkUnstructuredGrid >::New(); + copy->ShallowCopy( &mesh ); + return copy; + } + + GEOS_ERROR_IF( rank == 0 && localCells != totalCells, + GEOS_FMT( "Rank 0 must have the complete mesh. " + "Rank 0 has {} cells but total is {}", localCells, totalCells ) ); + + // KdTree: legacy path using VTK's built-in redistribution + if( method == ScatterMethod::kdtree ) + { +#ifdef GEOS_USE_MPI + vtkNew< vtkMPIController > controller; + vtkMPICommunicatorOpaqueComm vtkComm( &comm ); + vtkNew< vtkMPICommunicator > communicator; + communicator->InitializeExternal( &vtkComm ); + controller->SetCommunicator( communicator ); +#else + vtkNew< vtkDummyController > controller; +#endif + vtkMultiProcessController::SetGlobalController( controller ); + + vtkNew< vtkRedistributeDataSetFilter > rdsf; + rdsf->SetInputDataObject( &mesh ); + rdsf->SetNumberOfPartitions( size ); + rdsf->Update(); + + vtkSmartPointer< vtkDataSet > kdResult = vtkDataSet::SafeDownCast( rdsf->GetOutputDataObject( 0 ) ); + + vtkIdType const localAfter = kdResult->GetNumberOfCells(); + vtkIdType const totalAfter = MpiWrapper::allReduce( localAfter, MpiWrapper::Reduction::Sum, comm ); + + if( totalAfter != totalCells ) + { + GEOS_WARNING_IF( rank == 0, + GEOS_FMT( "{} cells lost during kdtree scatter ({} -> {}). " + "Falling back to contiguous scatter.", + totalCells - totalAfter, totalCells, totalAfter ) ); + // Fall through to the contiguous path below. + method = ScatterMethod::contiguous; + } + else + { + return kdResult; + } + } + + // Compute cell to rank assignment (rank 0 only) + stdVector< integer > cellRanks; + + if( rank == 0 ) + { + switch( method ) + { + case ScatterMethod::contiguous: + cellRanks = computeCellRanksContiguous( totalCells, size ); + break; + case ScatterMethod::cartesian: + { + GEOS_ERROR_IF( cartesianPartitions.size() < 3, + "Cartesian method requires 3 partition values (nx, ny, nz)" ); + integer const nx = cartesianPartitions[0], ny = cartesianPartitions[1], nz = cartesianPartitions[2]; + GEOS_ERROR_IF( nx * ny * nz != size, + GEOS_FMT( "partition grid {}x{}x{} = {} does not match MPI size {}", + nx, ny, nz, nx * ny * nz, size ) ); + cellRanks = computeCellRanksCartesian( mesh, nx, ny, nz ); + break; + } + case ScatterMethod::rcb: + cellRanks = computeCellRanksRCB( mesh, size ); + break; + default: + GEOS_ERROR( GEOS_FMT( "Unknown scatter method: {}", static_cast< integer >( method ) ) ); + break; + } + } + + // Scatter via binary tree + auto * inputGrid = vtkUnstructuredGrid::SafeDownCast( &mesh ); + GEOS_ERROR_IF( rank == 0 && inputGrid == nullptr, + "input must be a vtkUnstructuredGrid" ); + + // Non-rank-0 passes a dummy; binaryTreeScatter only uses rank 0's mesh. + vtkNew< vtkUnstructuredGrid > dummyGrid; + vtkUnstructuredGrid * gridPtr = ( rank == 0 ) ? inputGrid : dummyGrid.GetPointer(); + + vtkSmartPointer< vtkUnstructuredGrid > result = binaryTreeScatter( gridPtr, std::move( cellRanks ), comm ); + + // Validate cell conservation + vtkIdType const localAfter = result->GetNumberOfCells(); + vtkIdType const totalAfter = MpiWrapper::allReduce( localAfter, MpiWrapper::Reduction::Sum, comm ); + + GEOS_WARNING_IF( rank == 0 && totalAfter != totalCells, + GEOS_FMT( "{} cells lost during scatter ({} -> {})", + totalCells - totalAfter, totalCells, totalAfter ) ); + + return result; +} + + +} // namespace vtk +} // namespace geos diff --git a/src/coreComponents/mesh/generators/VTKMeshScattering.hpp b/src/coreComponents/mesh/generators/VTKMeshScattering.hpp new file mode 100644 index 00000000000..f7b364dd621 --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKMeshScattering.hpp @@ -0,0 +1,73 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_MESH_GENERATORS_VTKMESHSCATTERING_HPP +#define GEOS_MESH_GENERATORS_VTKMESHSCATTERING_HPP + +#include "common/DataTypes.hpp" +#include "common/format/EnumStrings.hpp" +#include "common/MpiWrapper.hpp" + +#include +#include + +namespace geos +{ +namespace vtk +{ + +/** + * @brief Method used to partition and scatter the mesh across MPI ranks. + */ +enum class ScatterMethod +{ + contiguous, ///< Assign contiguous blocks of cell indices to each rank (no geometry awareness). + cartesian, ///< Partition based on a user-specified Cartesian grid (nx × ny × nz). + rcb, ///< Recursive Coordinate Bisection along the longest bounding-box axis. + kdtree ///< VTK's built-in kd-tree redistribution (vtkRedistributeDataSetFilter). Retained for backward compatibility. +}; + +ENUM_STRINGS( ScatterMethod, + "contiguous", + "cartesian", + "rcb", + "kdtree" ); + +/** + * @brief Scatter a mesh held entirely on rank 0 to all MPI ranks. + * + * The mesh is first partitioned on rank 0 using the selected method, producing + * a cell to rank assignment vector. The data is then distributed via a binary-tree + * scatter pattern. + * + * @param[in] method the partitioning strategy to use + * @param[in] mesh the input mesh (must contain all cells on rank 0; empty on other ranks) + * @param[in] cartesianPartitions additional parameters for the cartesian partitioning method: + * - For @p cartesian: must contain at least 3 values (nx, ny, nz) + * with nx*ny*nz == MPI size. + * - Ignored for other methods. + * @param[in] comm the MPI communicator + * @return the local partition of the mesh on each rank + */ +vtkSmartPointer< vtkDataSet > +scatterMesh( ScatterMethod method, + vtkDataSet & mesh, + arrayView1d< integer const > cartesianPartitions, + MPI_Comm comm ); + +} // namespace vtk +} // namespace geos + +#endif // GEOS_MESH_GENERATORS_VTKMESHSCATTERING_HPP diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 0f5bea2f589..f39699d85ff 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -48,7 +48,6 @@ #include #include #include -#include #include #include #include @@ -752,113 +751,6 @@ redistributeByCellGraph( AllMeshes & input, } -/** - * @brief Scatter the mesh by blocks (no geometric information involved, assumes rank 0 has the full mesh) - * - * @param[in] mesh a vtk grid - * @return the vtk grid redistributed - */ -vtkSmartPointer< vtkDataSet > -scatterByBlock( vtkDataSet & mesh ) -{ - GEOS_MARK_FUNCTION; - - int const rank = MpiWrapper::commRank(); - int const size = MpiWrapper::commSize(); - - // Count total cells across all ranks - vtkIdType localCells = mesh.GetNumberOfCells(); - vtkIdType totalCells = MpiWrapper::allReduce( localCells, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); - - // Handle edge cases - if( totalCells == 0 ) - { - vtkNew< vtkUnstructuredGrid > emptyMesh; - return emptyMesh; - } - - if( size == 1 ) - { - vtkNew< vtkUnstructuredGrid > copy; - copy->DeepCopy( &mesh ); - return copy; - } - - // Verify rank 0 has the complete mesh for redistribution - if( rank == 0 && localCells != totalCells ) - { - GEOS_ERROR( GEOS_FMT( "Rank 0 must have the complete mesh. Rank 0 has {} cells but total is {}", - localCells, - totalCells ) ); - } - - // Scatter cells by contiguous blocks - vtkIdType cellsPerRank = totalCells / size; - vtkIdType remainder = totalCells % size; - - // Create partitioned dataset - vtkNew< vtkPartitionedDataSet > localParts; - if( rank == 0 ) - { - // Rank 0 has the full mesh, extract cells for each rank - for( int r = 0; r < size; ++r ) - { - vtkIdType rankStart = r * cellsPerRank + std::min( (vtkIdType)r, remainder ); - vtkIdType rankEnd = rankStart + cellsPerRank + (r < remainder ? 1 : 0); - - // Validate cell range - GEOS_ERROR_IF( rankStart< 0 || rankEnd > totalCells, - GEOS_FMT( "Invalid cell range for rank {}: [{}, {}) with total cells {}", - r, - rankStart, - rankEnd, - totalCells ) ); - - if( rankEnd > rankStart ) - { - // Add cells for this rank - vtkNew< vtkExtractCells > extractor; - extractor->SetInputDataObject( &mesh ); - extractor->AddCellRange( rankStart, rankEnd - 1 ); - extractor->Update(); - vtkUnstructuredGrid * extracted = extractor->GetOutput(); - localParts->SetPartition( r, extracted ); - } - else - { - // Create empty partition for ranks with no cells - vtkNew< vtkUnstructuredGrid > emptyPartition; - localParts->SetPartition( r, emptyPartition ); - } - } - } - else - { - // Other ranks have an empty mesh, but we still need to create the - // partitioned data set structure. - localParts->SetNumberOfPartitions( size ); - for( int r = 0; r < size; ++r ) - { - vtkNew< vtkUnstructuredGrid > emptyPartition; - localParts->SetPartition( r, emptyPartition ); - } - } - - //Send cells to appropriate ranks - vtkSmartPointer< vtkUnstructuredGrid > result = vtk::redistribute( *localParts, MPI_COMM_GEOS ); - - // Final validation - vtkIdType finalLocalCells = result->GetNumberOfCells(); - vtkIdType finalTotalCells = MpiWrapper::allReduce( finalLocalCells, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); - - GEOS_ERROR_IF( finalTotalCells != totalCells, - GEOS_FMT( "Block redistribution lost cells: started with {}, ended with {}", - totalCells, - finalTotalCells ) ); - - return result; -} - /** * @brief Classify cells by dimension @@ -1399,45 +1291,6 @@ redistributeByAreaGraphAndLayer( AllMeshes & input, return AllMeshes( vtk::redistribute( *splitMesh, MPI_COMM_GEOS ), {} ); } -/** - * @brief Redistributes the mesh using a Kd-Tree - * - * @param[in] mesh a vtk grid - * @return the vtk grid redistributed - */ -vtkSmartPointer< vtkDataSet > -redistributeByKdTree( vtkDataSet & mesh ) -{ - GEOS_MARK_FUNCTION; - - // Count input cells for verification - vtkIdType localInputCells = mesh.GetNumberOfCells(); - vtkIdType globalInputCells = MpiWrapper::allReduce( localInputCells, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); - - // Use a VTK filter which employs a kd-tree partition internally - vtkNew< vtkRedistributeDataSetFilter > rdsf; - rdsf->SetInputDataObject( &mesh ); - rdsf->SetNumberOfPartitions( MpiWrapper::commSize() ); - rdsf->Update(); - - vtkSmartPointer< vtkDataSet > result = vtkDataSet::SafeDownCast( rdsf->GetOutputDataObject( 0 ) ); - - // Verify we didn't lose any cells - vtkIdType localOutputCells = result->GetNumberOfCells(); - vtkIdType globalOutputCells = MpiWrapper::allReduce( localOutputCells, MpiWrapper::Reduction::Sum, MPI_COMM_GEOS ); - - if( globalOutputCells != globalInputCells ) - { - if( MpiWrapper::commRank() == 0 ) - { - GEOS_WARNING( GEOS_FMT( "VTK KdTree redistribution lost {} elements! Falling back to block redistribution.", - globalInputCells - globalOutputCells ) ); - } - return scatterByBlock( mesh ); - } - - return result; -} stdVector< int > findNeighborRanks( stdVector< vtkBoundingBox > boundingBoxes ) @@ -1654,16 +1507,19 @@ redistributeMeshes( integer const logLevel, vtkSmartPointer< vtkDataSet > loadedMesh, stdMap< string, vtkSmartPointer< vtkDataSet > > & namesToFractures, MPI_Comm const comm, + ScatterMethod const scatterMethod, + arrayView1d< int const > partitions, PartitionMethod const method, int const partitionRefinement, int const useGlobalIds, - string const & structuredIndexAttributeName, - int const numPartZ ) + string const & structuredIndexAttributeName ) { GEOS_MARK_FUNCTION; int const numRanks = MpiWrapper::commSize( comm ); int const rank = MpiWrapper::commRank( comm ); + int const numPartZ = structuredIndexAttributeName.empty() ? 1 : partitions[2]; + stdVector< vtkSmartPointer< vtkDataSet > > fractures; for( auto & nameToFracture: namesToFractures ) { @@ -1702,8 +1558,8 @@ redistributeMeshes( integer const logLevel, // Extract 3D cells for initial redistribution vtkSmartPointer< vtkUnstructuredGrid > cells3D = extractCellsByIndices( *mesh, cells3DIndices.toViewConst() ); - // Redistribute the 3D mesh over all ranks using simple octree partitions - vtkSmartPointer< vtkDataSet > redistributed3D = redistributeByKdTree( *cells3D ); + // Scatter the 3D mesh over all ranks using the selected method + vtkSmartPointer< vtkDataSet > redistributed3D = scatterMesh( scatterMethod, *cells3D, partitions, comm ); // Check if a rank does not have a cell after the redistribution if( MpiWrapper::min( redistributed3D->GetNumberOfCells(), comm ) == 0 ) diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index ef660d7641c..af07e95592d 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.hpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.hpp @@ -23,6 +23,7 @@ #include "common/DataTypes.hpp" #include "common/MpiWrapper.hpp" #include "mesh/generators/CellBlockManager.hpp" +#include "mesh/generators/VTKMeshScattering.hpp" #include #include @@ -150,11 +151,12 @@ findNeighborRanks( stdVector< vtkBoundingBox > boundingBoxes ); * @param[in] loadedMesh the mesh that was loaded on one or several MPI ranks * @param[in] namesToFractures the fracture meshes * @param[in] comm the MPI communicator - * @param[in] method the partitioning method + * @param[in] scatterMethod the initial scatter method + * @param[in] partitions the -x/-y/-z partition counts (used for cartesian scatter and structured index) + * @param[in] method the partitioning method for refinement * @param[in] partitionRefinement number of graph partitioning refinement cycles * @param[in] useGlobalIds controls whether global id arrays from the vtk input should be used * @param[in] structuredIndexAttributeName VTK array name for structured index attribute, if present - * @param[in] numPartZ number of MPI partitions in Z direction (only if @p structuredIndexAttributeName is used) * @return the vtk grid redistributed */ AllMeshes @@ -162,11 +164,12 @@ redistributeMeshes( integer const logLevel, vtkSmartPointer< vtkDataSet > loadedMesh, stdMap< string, vtkSmartPointer< vtkDataSet > > & namesToFractures, MPI_Comm const comm, + ScatterMethod const scatterMethod, + arrayView1d< int const > partitions, PartitionMethod const method, int const partitionRefinement, int const useGlobalIds, - string const & structuredIndexAttributeName, - int const numPartZ ); + string const & structuredIndexAttributeName ); /** * @brief Collect lists of VTK cell indices organized by type and attribute value. diff --git a/src/coreComponents/mesh/unitTests/CMakeLists.txt b/src/coreComponents/mesh/unitTests/CMakeLists.txt index a365e1f0b0a..5198a317e3d 100644 --- a/src/coreComponents/mesh/unitTests/CMakeLists.txt +++ b/src/coreComponents/mesh/unitTests/CMakeLists.txt @@ -61,4 +61,14 @@ set( nranks 8 ) target_link_libraries( ${test_name} PRIVATE "stdc++fs" ) endif () endforeach() + + # VTKMeshScattering test — 4 MPI ranks + blt_add_executable( NAME testVTKMeshScattering_mpi + SOURCES testVTKMeshScattering.cpp + OUTPUT_DIR ${TEST_OUTPUT_DIRECTORY} + DEPENDS_ON ${dependencyList} ${tplDependencyList} ) + + geos_add_test( NAME testVTKMeshScattering_mpi + COMMAND testVTKMeshScattering_mpi + NUM_MPI_TASKS 4 ) endif() diff --git a/src/coreComponents/mesh/unitTests/testVTKMeshScattering.cpp b/src/coreComponents/mesh/unitTests/testVTKMeshScattering.cpp new file mode 100644 index 00000000000..100c113db2b --- /dev/null +++ b/src/coreComponents/mesh/unitTests/testVTKMeshScattering.cpp @@ -0,0 +1,302 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file testVTKMeshScattering.cpp + * @brief Unit tests for VTKMeshScattering (MPI, 4 ranks). + */ + +#include "common/DataTypes.hpp" +#include "common/MpiWrapper.hpp" +#include "mesh/generators/VTKMeshScattering.hpp" + +#include +#include +#include +#include +#include +#include + +#include + +using namespace geos; +using namespace geos::vtk; + +namespace +{ + +/** + * @brief Build a regular 4x4x4 hexahedral grid on rank 0. + * + * The mesh occupies [0, 4] x [0, 4] x [0, 4] with 64 hex cells, + * 125 points, a cell data array ("CellId") and a point data array ("PointId"). + * Returns an empty grid on non-zero ranks. + */ +vtkSmartPointer< vtkUnstructuredGrid > buildTestMesh( MPI_Comm comm ) +{ + integer const rank = MpiWrapper::commRank( comm ); + auto mesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + + if( rank != 0 ) + { + return mesh; + } + + constexpr integer N = 4; // cells per dimension + constexpr integer NP = N + 1; + + // Points: (NP)^3 = 125 points + vtkNew< vtkPoints > points; + points->SetDataTypeToDouble(); + points->Allocate( NP * NP * NP ); + + for( integer k = 0; k <= N; ++k ) + { + for( integer j = 0; j <= N; ++j ) + { + for( integer i = 0; i <= N; ++i ) + { + points->InsertNextPoint( static_cast< real64 >( i ), + static_cast< real64 >( j ), + static_cast< real64 >( k ) ); + } + } + } + mesh->SetPoints( points ); + + // Hex cells: N^3 = 64 cells + mesh->Allocate( N * N * N ); + for( integer k = 0; k < N; ++k ) + { + for( integer j = 0; j < N; ++j ) + { + for( integer i = 0; i < N; ++i ) + { + vtkIdType const base = i + j * NP + k * NP * NP; + vtkIdType hex[8] = { + base, + base + 1, + base + NP + 1, + base + NP, + base + NP * NP, + base + NP * NP + 1, + base + NP * NP + NP + 1, + base + NP * NP + NP + }; + mesh->InsertNextCell( VTK_HEXAHEDRON, 8, hex ); + } + } + } + + // Cell data: sequential cell IDs + vtkNew< vtkIdTypeArray > cellIds; + cellIds->SetName( "CellId" ); + cellIds->SetNumberOfTuples( N * N * N ); + for( vtkIdType i = 0; i < N * N * N; ++i ) + { + cellIds->SetValue( i, i ); + } + mesh->GetCellData()->AddArray( cellIds ); + + // Point data: sequential point IDs + vtkNew< vtkIdTypeArray > pointIds; + pointIds->SetName( "PointId" ); + pointIds->SetNumberOfTuples( NP * NP * NP ); + for( vtkIdType i = 0; i < NP * NP * NP; ++i ) + { + pointIds->SetValue( i, i ); + } + mesh->GetPointData()->AddArray( pointIds ); + + return mesh; +} + +/** + * @brief Helper: scatter with a given method and return the result as vtkUnstructuredGrid. + */ +vtkSmartPointer< vtkUnstructuredGrid > +scatter( ScatterMethod method, + vtkUnstructuredGrid & mesh, + arrayView1d< integer const > parts, + MPI_Comm comm ) +{ + vtkSmartPointer< vtkDataSet > result = scatterMesh( method, mesh, parts, comm ); + return vtkUnstructuredGrid::SafeDownCast( result ); +} + +} // anonymous namespace + + +class VTKMeshScatteringTest : public ::testing::Test +{ +protected: + void SetUp() override + { + comm = MPI_COMM_GEOS; + rank = MpiWrapper::commRank( comm ); + size = MpiWrapper::commSize( comm ); + mesh = buildTestMesh( comm ); + + parts.resize( 3 ); + parts[0] = 2; parts[1] = 2; parts[2] = 1; + } + + MPI_Comm comm; + integer rank; + integer size; + vtkSmartPointer< vtkUnstructuredGrid > mesh; + array1d< integer > parts; + + static constexpr vtkIdType totalCells = 64; + static constexpr vtkIdType totalPoints = 125; +}; + + +/// All cells must survive the scatter (no loss, no duplication). +TEST_F( VTKMeshScatteringTest, CellConservation ) +{ + stdVector< ScatterMethod > methods = { + ScatterMethod::contiguous, + ScatterMethod::cartesian, + ScatterMethod::rcb, + ScatterMethod::kdtree + }; + + for( auto method : methods ) + { + auto result = scatter( method, *mesh, parts.toViewConst(), comm ); + ASSERT_NE( result, nullptr ); + + vtkIdType const localCells = result->GetNumberOfCells(); + vtkIdType const globalCells = MpiWrapper::allReduce( localCells, MpiWrapper::Reduction::Sum, comm ); + + EXPECT_EQ( globalCells, totalCells ) + << "Cell conservation failed for method " << toString( method ); + } +} + + +/// With 64 cells and 4 ranks, every rank must get exactly 16 cells. +TEST_F( VTKMeshScatteringTest, LoadBalance ) +{ + stdVector< ScatterMethod > methods = { + ScatterMethod::contiguous, + ScatterMethod::cartesian, + ScatterMethod::rcb + }; + + vtkIdType const expectedPerRank = totalCells / size; + + for( auto method : methods ) + { + auto result = scatter( method, *mesh, parts.toViewConst(), comm ); + vtkIdType const localCells = result->GetNumberOfCells(); + + EXPECT_EQ( localCells, expectedPerRank ) + << "Rank " << rank << " got " << localCells << " cells for method " << toString( method ); + } +} + + +/// Cell data and point data arrays must survive the scatter. +TEST_F( VTKMeshScatteringTest, DataArrayPreservation ) +{ + auto result = scatter( ScatterMethod::rcb, *mesh, parts.toViewConst(), comm ); + + // Cell data + vtkDataArray * cellIdArr = result->GetCellData()->GetArray( "CellId" ); + ASSERT_NE( cellIdArr, nullptr ) << "CellId array lost during scatter"; + EXPECT_EQ( cellIdArr->GetNumberOfTuples(), result->GetNumberOfCells() ); + + // Point data + vtkDataArray * pointIdArr = result->GetPointData()->GetArray( "PointId" ); + ASSERT_NE( pointIdArr, nullptr ) << "PointId array lost during scatter"; + EXPECT_EQ( pointIdArr->GetNumberOfTuples(), result->GetNumberOfPoints() ); +} + + +/// Cartesian partitioning with a 2x2x1 grid should place cells in the correct spatial quadrant. +TEST_F( VTKMeshScatteringTest, CartesianSpatialCorrectness ) +{ + auto result = scatter( ScatterMethod::cartesian, *mesh, parts.toViewConst(), comm ); + + // With a 2x2x1 partition on [0,4]^3, the split is at x=2 and y=2. + // Rank = ix + 2*iy, where ix = (centroid.x < 2 ? 0 : 1), iy = (centroid.y < 2 ? 0 : 1). + for( vtkIdType c = 0; c < result->GetNumberOfCells(); ++c ) + { + real64 bounds[6]; + result->GetCell( c )->GetBounds( bounds ); + real64 const cx = ( bounds[0] + bounds[1] ) * 0.5; + real64 const cy = ( bounds[2] + bounds[3] ) * 0.5; + + integer const expectedIx = ( cx < 2.0 ) ? 0 : 1; + integer const expectedIy = ( cy < 2.0 ) ? 0 : 1; + integer const expectedRank = expectedIx + 2 * expectedIy; + + EXPECT_EQ( expectedRank, rank ) + << "Cell centroid (" << cx << ", " << cy << ") on wrong rank"; + } +} + + +/// Scatter on a single rank (size==1 early exit) returns a deep copy. +TEST_F( VTKMeshScatteringTest, SingleRankNoOp ) +{ + // Create a sub-communicator containing only rank 0. + // On other ranks we just verify we don't crash. + integer const color = ( rank == 0 ) ? 0 : MPI_UNDEFINED; + MPI_Comm singleComm = MpiWrapper::commSplit( comm, color, rank ); + + if( rank == 0 ) + { + auto result = scatter( ScatterMethod::rcb, *mesh, parts.toViewConst(), singleComm ); + EXPECT_EQ( result->GetNumberOfCells(), totalCells ); + EXPECT_EQ( result->GetNumberOfPoints(), totalPoints ); + + MpiWrapper::commFree( singleComm ); + } +} + + +/// An empty mesh on all ranks should return an empty grid without error. +TEST_F( VTKMeshScatteringTest, EmptyMesh ) +{ + auto emptyMesh = vtkSmartPointer< vtkUnstructuredGrid >::New(); + auto result = scatter( ScatterMethod::rcb, *emptyMesh, parts.toViewConst(), comm ); + + ASSERT_NE( result, nullptr ); + EXPECT_EQ( result->GetNumberOfCells(), 0 ); +} + + +int main( int argc, char * * argv ) +{ + MpiWrapper::init( &argc, &argv ); + MPI_COMM_GEOS = MpiWrapper::commDup( MPI_COMM_WORLD ); + + ::testing::InitGoogleTest( &argc, argv ); + + integer const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); + if( rank != 0 ) + { + ::testing::TestEventListeners & listeners = ::testing::UnitTest::GetInstance()->listeners(); + delete listeners.Release( listeners.default_result_printer() ); + } + + integer result = RUN_ALL_TESTS(); + + MpiWrapper::finalize(); + return result; +} From 73a367e494ac675f0f17138902e80a35d93abfb5 Mon Sep 17 00:00:00 2001 From: DENEL Bertrand Date: Mon, 27 Apr 2026 14:55:56 -0500 Subject: [PATCH 2/2] Fixed dependency --- src/coreComponents/mesh/benchmarks/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/mesh/benchmarks/CMakeLists.txt b/src/coreComponents/mesh/benchmarks/CMakeLists.txt index 3c877949bcb..160b6bf8bec 100644 --- a/src/coreComponents/mesh/benchmarks/CMakeLists.txt +++ b/src/coreComponents/mesh/benchmarks/CMakeLists.txt @@ -2,11 +2,13 @@ # Mesh benchmarks # +set( dependencyList mesh ${parallelDeps} ${tplDependencyList} ) + blt_add_executable( NAME benchmarkScatterMethods SOURCES benchmarkScatterMethods.cpp OUTPUT_DIR ${BENCHMARK_OUTPUT_DIRECTORY} - DEPENDS_ON mesh + DEPENDS_ON ${dependencyList} FOLDER coreComponents/mesh/benchmarks )