diff --git a/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp b/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp index a2f13b88748..e2697ac1599 100644 --- a/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp +++ b/src/coreComponents/integrationTests/meshTests/testVTKImport.cpp @@ -51,11 +51,17 @@ using namespace geos::dataRepository; template< class V > -void TestMeshImport( string const & meshFilePath, V const & validate, string const fractureName="" ) +void TestMeshImport( string const & meshFilePath, V const & validate, string const fractureName="", string const scatterMethod="" ) { // Automatically use global IDs when fractures are present string const useGlobalIdsStr = fractureName.empty() ? "0" : "1"; + string scatterAttr; + if( !scatterMethod.empty() ) + { + scatterAttr = GEOS_FMT( "scatterMethod=\"{}\"", scatterMethod ); + } + string const pattern = R"xml( )xml"; string const meshNode = GEOS_FMT( pattern, meshFilePath, useGlobalIdsStr, + scatterAttr, fractureName.empty() ? "" : "faceBlocks=\"{" + fractureName + "}\"" ); xmlWrapper::xmlDocument xmlDocument; diff --git a/src/coreComponents/mesh/CMakeLists.txt b/src/coreComponents/mesh/CMakeLists.txt index 65db05207f1..133179ffac7 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/VTKSuperCellPartitioning.hpp generators/VTKUtilities.hpp @@ -224,6 +225,7 @@ if( ENABLE_VTK ) generators/VTKHierarchicalDataSource.cpp generators/VTKMeshGenerator.cpp generators/VTKMeshGeneratorTools.cpp + generators/VTKMeshScattering.cpp generators/VTKWellGenerator.cpp generators/VTKSuperCellPartitioning.cpp generators/VTKUtilities.cpp @@ -270,3 +272,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..160b6bf8bec --- /dev/null +++ b/src/coreComponents/mesh/benchmarks/CMakeLists.txt @@ -0,0 +1,15 @@ +# +# Mesh benchmarks +# + +set( dependencyList mesh ${parallelDeps} ${tplDependencyList} ) + +blt_add_executable( + NAME benchmarkScatterMethods + SOURCES benchmarkScatterMethods.cpp + OUTPUT_DIR ${BENCHMARK_OUTPUT_DIRECTORY} + DEPENDS_ON ${dependencyList} + 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 4a881414782..4ecf47fc41a 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), " + "rcb (recursive coordinate bisection), " + "kdtree (VTK built-in kd-tree, default; automatically falls back to rcb when fractures are present)" ); registerWrapper( viewKeyStruct::partitionFractureWeightString(), &m_partitionFractureWeight ). setInputFlag( InputFlags::OPTIONAL ). @@ -140,7 +149,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..." ); { @@ -157,7 +176,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()) @@ -210,12 +229,13 @@ void VTKMeshGenerator::fillCellBlockManager( CellBlockManager & cellBlockManager allMeshes.getMainMesh(), allMeshes.getFaceBlocks(), comm, + m_scatterMethod, + partitions.toViewConst(), m_partitionMethod, m_partitionRefinement, m_partitionFractureWeight, 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 e31ada4bb1e..5a50ef21348 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" @@ -115,6 +116,7 @@ class VTKMeshGenerator : public ExternalMeshGeneratorBase constexpr static char const * partitionRefinementString() { return "partitionRefinement"; } constexpr static char const * partitionMethodString() { return "partitionMethod"; } constexpr static char const * partitionFractureWeightString() { return "partitionFractureWeight"; } + constexpr static char const * scatterMethodString() { return "scatterMethod"; } constexpr static char const * useGlobalIdsString() { return "useGlobalIds"; } constexpr static char const * dataSourceString() { return "dataSourceName"; } }; @@ -171,6 +173,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..e09c0e72e78 --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKMeshScattering.cpp @@ -0,0 +1,843 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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; +} + +} // anonymous namespace + +// ============================================================================ +// Public API +// ============================================================================ + +// ---------------------------------------------------------------------------- +// 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 > +scatterByRankAssignment( 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; +} + +// ---------------------------------------------------------------------------- +// Per-cell rank assignment dispatcher (rank 0 only) +// ---------------------------------------------------------------------------- +stdVector< integer > +computeCellRanks( ScatterMethod method, + vtkDataSet & mesh, + arrayView1d< integer const > cartesianPartitions, + integer numRanks ) +{ + GEOS_ERROR_IF( method == ScatterMethod::kdtree, + "computeCellRanks: kdtree method does not produce a per-cell rank assignment; " + "use scatterMesh() for kdtree." ); + + vtkIdType const numCells = mesh.GetNumberOfCells(); + stdVector< integer > cellRanks; + switch( method ) + { + case ScatterMethod::contiguous: + cellRanks = computeCellRanksContiguous( numCells, numRanks ); + 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 != numRanks, + GEOS_FMT( "partition grid {}x{}x{} = {} does not match MPI size {}", + nx, ny, nz, nx * ny * nz, numRanks ) ); + cellRanks = computeCellRanksCartesian( mesh, nx, ny, nz ); + break; + } + case ScatterMethod::rcb: + cellRanks = computeCellRanksRCB( mesh, numRanks ); + break; + default: + GEOS_ERROR( GEOS_FMT( "Unknown scatter method: {}", static_cast< integer >( method ) ) ); + break; + } + return cellRanks; +} + + +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 ) + { + cellRanks = computeCellRanks( method, mesh, cartesianPartitions, size ); + } + + // 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; scatterByRankAssignment only uses rank 0's mesh. + vtkNew< vtkUnstructuredGrid > dummyGrid; + vtkUnstructuredGrid * gridPtr = ( rank == 0 ) ? inputGrid : dummyGrid.GetPointer(); + + vtkSmartPointer< vtkUnstructuredGrid > result = scatterByRankAssignment( 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..ca0aa650056 --- /dev/null +++ b/src/coreComponents/mesh/generators/VTKMeshScattering.hpp @@ -0,0 +1,109 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 +#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 Compute the per-cell destination rank on rank 0 for a given scatter method. + * + * This is the assignment step of @ref scatterMesh, exposed so other code (e.g. the + * super-cell aware scatter) can build its own assignment on top of the same algorithms. + * + * @param[in] method the partitioning strategy. @p kdtree is not supported here + * (it does not produce a separable per-cell assignment). + * @param[in] mesh the input mesh (must contain all cells on rank 0). + * @param[in] cartesianPartitions {nx, ny, nz}, only used for @p cartesian. + * @param[in] numRanks total number of MPI ranks. + * @return rank-0 only: @c cellRanks[i] is the destination rank for cell @p i; empty on other ranks. + */ +stdVector< integer > +computeCellRanks( ScatterMethod method, + vtkDataSet & mesh, + arrayView1d< integer const > cartesianPartitions, + integer numRanks ); + +/** + * @brief Ship cells from rank 0 to the ranks specified by @p assignment via a binary-tree exchange. + * + * This is the shipping step of @ref scatterMesh, exposed so other code can drive it with a + * custom @p assignment (e.g. atom-aware partitioning that keeps fracture-connected cells together). + * + * @param[in] inputMesh rank 0: source mesh with all cells. Other ranks: ignored (may be empty). + * @param[in] assignment rank 0 only: @c assignment[i] is the destination rank for cell @p i. + * @param[in] comm the MPI communicator. + * @return the local subset of the mesh assigned to this rank. + */ +vtkSmartPointer< vtkUnstructuredGrid > +scatterByRankAssignment( vtkUnstructuredGrid * inputMesh, + stdVector< integer > assignment, + MPI_Comm comm ); + +/** + * @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/VTKSuperCellPartitioning.cpp b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp index 679cdc95f11..c58c4cafd55 100644 --- a/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp +++ b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.cpp @@ -260,245 +260,117 @@ SuperCellInfo reconstructSuperCellInfo( vtkSmartPointer< vtkUnstructuredGrid > m vtkSmartPointer< vtkDataSet > redistributeBySuperCellBlocks( vtkSmartPointer< vtkUnstructuredGrid > cells3D, MPI_Comm comm, - InitialDistributionStrategy strategy ) + ScatterMethod scatterMethod, + arrayView1d< integer const > cartesianPartitions ) { GEOS_MARK_FUNCTION; int const rank = MpiWrapper::commRank( comm ); int const numRanks = MpiWrapper::commSize( comm ); - vtkSmartPointer< vtkPartitionedDataSet > partitionedMesh = - vtkSmartPointer< vtkPartitionedDataSet >::New(); + if( scatterMethod == ScatterMethod::kdtree ) + { + GEOS_LOG_RANK_0( "scatterMethod=kdtree is not supported with fractures (cannot preserve " + "super-cell atomicity). Automatically falling back to rcb." ); + scatterMethod = ScatterMethod::rcb; + } + + // Per-cell rank assignment is computed on rank 0 (empty elsewhere) then consumed by + // scatterByRankAssignment. Super-cell atomicity is enforced by assigning one rank + // per super-cell and propagating it to every cell of that super-cell. + stdVector< integer > cellRanks; if( rank == 0 ) { vtkIdTypeArray * superCellIdArray = vtkIdTypeArray::SafeDownCast( cells3D->GetCellData()->GetArray( "SuperCellId" ) ); - GEOS_ERROR_IF( !superCellIdArray, "SuperCellId array not found" ); vtkIdType const numCells = cells3D->GetNumberOfCells(); - // Build super-cell metadata - stdMap< vtkIdType, stdVector< vtkIdType > > superCellToLocalCells; - + // Group cells by super-cell ID. The iteration order of the map gives a stable + // 0..numSuperCells-1 indexing used by the virtual "atom mesh" below. + stdMap< vtkIdType, stdVector< vtkIdType > > superCellToCells; for( vtkIdType i = 0; i < numCells; ++i ) { - vtkIdType scId = superCellIdArray->GetValue( i ); - superCellToLocalCells.get_inserted( scId ).push_back( i ); + superCellToCells.get_inserted( superCellIdArray->GetValue( i ) ).push_back( i ); } - vtkIdType numSuperCells = superCellToLocalCells.size(); + vtkIdType const numSuperCells = superCellToCells.size(); - GEOS_LOG_RANK_0( GEOS_FMT( "Initial distribution: {} super-cells across {} ranks ({} ordering)", + GEOS_LOG_RANK_0( GEOS_FMT( "Initial distribution: {} super-cells across {} ranks (method={})", numSuperCells, numRanks, - (strategy == InitialDistributionStrategy::MORTON ? "MORTON" : "BLOCK") ) ); - - // ----------------------------------------------------------------------- - // Step 1: Build super-cell list and optionally sort by spatial locality - // ----------------------------------------------------------------------- - struct SuperCellDistributionInfo - { - vtkIdType scId; - stdVector< vtkIdType > cellIndices; - stdArray< double, 3 > centroid; // Only computed for Morton strategy - }; - - stdVector< SuperCellDistributionInfo > superCells; - superCells.reserve( numSuperCells ); - - if( strategy == InitialDistributionStrategy::MORTON ) - { - vtkPoints * meshPoints = cells3D->GetPoints(); - - // Compute super-cell centroids directly - for( auto const & [scId, cellIndices] : superCellToLocalCells ) - { - stdArray< double, 3 > centroid = {0.0, 0.0, 0.0}; - vtkIdType totalPoints = 0; - - // Accumulate all points from all cells in this super-cell - for( vtkIdType cellIdx : cellIndices ) - { - vtkIdType npts; - const vtkIdType * pts; - cells3D->GetCellPoints( cellIdx, npts, pts ); + EnumStrings< ScatterMethod >::toString( scatterMethod ) ) ); - for( vtkIdType i = 0; i < npts; ++i ) - { - double pt[3]; - meshPoints->GetPoint( pts[i], pt ); - centroid[0] += pt[0]; - centroid[1] += pt[1]; - centroid[2] += pt[2]; - } - totalPoints += npts; - } - - // Average over all points - centroid[0] /= totalPoints; - centroid[1] /= totalPoints; - centroid[2] /= totalPoints; + // Build a virtual mesh with one VTK_VERTEX cell per super-cell, placed at the + // super-cell's centroid (average over all points of all member cells). This lets + // computeCellRanks() see super-cells as atomic and partition them as such. + vtkNew< vtkPoints > atomPoints; + atomPoints->SetNumberOfPoints( numSuperCells ); - superCells.push_back( SuperCellDistributionInfo{ scId, cellIndices, centroid } ); - } + vtkNew< vtkUnstructuredGrid > atomMesh; + atomMesh->SetPoints( atomPoints ); + atomMesh->Allocate( numSuperCells ); - // Find bounding box - double minCoord[3] = {std::numeric_limits< double >::max(), - std::numeric_limits< double >::max(), - std::numeric_limits< double >::max()}; - double maxCoord[3] = {std::numeric_limits< double >::lowest(), - std::numeric_limits< double >::lowest(), - std::numeric_limits< double >::lowest()}; + stdVector< integer > cellToAtom( numCells ); + vtkPoints * meshPoints = cells3D->GetPoints(); - for( auto const & sc : superCells ) - { - for( int d = 0; d < 3; ++d ) - { - minCoord[d] = std::min( minCoord[d], sc.centroid[d] ); - maxCoord[d] = std::max( maxCoord[d], sc.centroid[d] ); - } - } + vtkIdType atomIdx = 0; + for( auto const & kv : superCellToCells ) + { + stdVector< vtkIdType > const & cellIndices = kv.second; - // Morton encoding - auto computeMorton = []( stdArray< double, 3 > const & centroid, - double bounds_min[3], - double bounds_max[3] ) -> uint64_t + real64 cx = 0.0, cy = 0.0, cz = 0.0; + vtkIdType totalPoints = 0; + for( vtkIdType cellIdx : cellIndices ) { - auto normalize = [&]( double val, int dim ) -> uint32_t + vtkIdType npts; + vtkIdType const * pts; + cells3D->GetCellPoints( cellIdx, npts, pts ); + for( vtkIdType i = 0; i < npts; ++i ) { - double range = bounds_max[dim] - bounds_min[dim]; - if( range < 1e-10 ) - return 0; - double norm = (val - bounds_min[dim]) / range; - norm = std::max( 0.0, std::min( 1.0, norm ) ); - return static_cast< uint32_t >( norm * ((1u << 21) - 1) ); - }; - - uint32_t x = normalize( centroid[0], 0 ); - uint32_t y = normalize( centroid[1], 1 ); - uint32_t z = normalize( centroid[2], 2 ); - - uint64_t code = 0; - for( int i = 0; i < 21; ++i ) - { - code |= ((x & (1u << i)) ? (1ull << (3*i)) : 0); - code |= ((y & (1u << i)) ? (1ull << (3*i + 1)) : 0); - code |= ((z & (1u << i)) ? (1ull << (3*i + 2)) : 0); + real64 p[3]; + meshPoints->GetPoint( pts[i], p ); + cx += p[0]; + cy += p[1]; + cz += p[2]; } - return code; - }; - - // Sort by Morton code - std::sort( superCells.begin(), superCells.end(), - [&]( SuperCellDistributionInfo const & a, SuperCellDistributionInfo const & b ) - { - return computeMorton( a.centroid, minCoord, maxCoord ) < - computeMorton( b.centroid, minCoord, maxCoord ); - } ); - - } - else - { - // BLOCK: Simple ordering by super-cell ID (no centroid computation needed) - for( auto const & [scId, cellIndices] : superCellToLocalCells ) - { - superCells.push_back( SuperCellDistributionInfo{ scId, cellIndices, {0.0, 0.0, 0.0} } ); + totalPoints += npts; + cellToAtom[cellIdx] = static_cast< integer >( atomIdx ); } - // Sort by super-cell ID for deterministic partitioning - std::sort( superCells.begin(), superCells.end(), - []( SuperCellDistributionInfo const & a, SuperCellDistributionInfo const & b ) - { - return a.scId < b.scId; - } ); - } - - // ----------------------------------------------------------------------- - // Step 2: Assign super-cells to ranks in contiguous blocks - // ----------------------------------------------------------------------- - array1d< int64_t > cellPartitions( numCells ); - stdVector< vtkIdType > cellsPerRank( numRanks, 0 ); + real64 const inv = (totalPoints > 0) ? 1.0 / totalPoints : 0.0; + atomPoints->SetPoint( atomIdx, cx * inv, cy * inv, cz * inv ); - vtkIdType superCellsPerRank = (numSuperCells + numRanks - 1) / numRanks; - - for( vtkIdType scIdx = 0; scIdx < numSuperCells; ++scIdx ) - { - int targetRank = std::min( static_cast< int >(scIdx / superCellsPerRank), numRanks - 1 ); + vtkIdType const pid = atomIdx; + atomMesh->InsertNextCell( VTK_VERTEX, 1, &pid ); - // All cells in this super-cell go to the same rank - for( vtkIdType cellIdx : superCells[scIdx].cellIndices ) - { - cellPartitions[cellIdx] = targetRank; - cellsPerRank[targetRank]++; - } + ++atomIdx; } - // ----------------------------------------------------------------------- - // Step 3: Build partitions - // ----------------------------------------------------------------------- - partitionedMesh->SetNumberOfPartitions( numRanks ); - - for( int r = 0; r < numRanks; ++r ) - { - vtkNew< vtkIdList > cellsForRank; - - for( vtkIdType i = 0; i < numCells; ++i ) - { - if( cellPartitions[i] == r ) - { - cellsForRank->InsertNextId( i ); - } - } - - if( cellsForRank->GetNumberOfIds() == 0 ) - { - vtkNew< vtkUnstructuredGrid > emptyPart; - partitionedMesh->SetPartition( r, emptyPart ); - continue; - } - - vtkNew< vtkExtractCells > extractor; - extractor->SetInputData( cells3D ); - extractor->SetCellList( cellsForRank ); - extractor->Update(); + // One rank per super-cell, then expand to per-cell so atomicity is preserved. + stdVector< integer > const atomRanks = + computeCellRanks( scatterMethod, *atomMesh, cartesianPartitions, numRanks ); - vtkSmartPointer< vtkUnstructuredGrid > partition = - vtkUnstructuredGrid::SafeDownCast( extractor->GetOutput() ); - - partitionedMesh->SetPartition( r, partition ); - } - } - else - { - // Other ranks: create empty partitioned dataset - partitionedMesh->SetNumberOfPartitions( numRanks ); - for( int r = 0; r < numRanks; ++r ) + cellRanks.resize( numCells ); + for( vtkIdType c = 0; c < numCells; ++c ) { - vtkNew< vtkUnstructuredGrid > emptyPart; - partitionedMesh->SetPartition( r, emptyPart ); + cellRanks[c] = atomRanks[ cellToAtom[c] ]; } } - // ----------------------------------------------------------------------- - // Step 4: Redistribute using VTK - // ----------------------------------------------------------------------- - vtkSmartPointer< vtkDataSet > result = vtk::redistribute( *partitionedMesh, comm ); + // All ranks ship cells + vtkSmartPointer< vtkUnstructuredGrid > result = + scatterByRankAssignment( cells3D.Get(), std::move( cellRanks ), comm ); - partitionedMesh = nullptr; if( rank == 0 ) { cells3D = nullptr; } - vtkIdType localCells = result->GetNumberOfCells(); - - // Verify SuperCellId array survived redistribution - if( localCells > 0 ) + // SuperCellId must survive redistribution + if( result->GetNumberOfCells() > 0 ) { - vtkIdTypeArray * resultSuperCellIdArray = - vtkIdTypeArray::SafeDownCast( result->GetCellData()->GetArray( "SuperCellId" ) ); - - GEOS_ERROR_IF( !resultSuperCellIdArray, + GEOS_ERROR_IF( !vtkIdTypeArray::SafeDownCast( result->GetCellData()->GetArray( "SuperCellId" ) ), GEOS_FMT( "Rank {}: SuperCellId array lost during redistribution", rank ) ); } diff --git a/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp index 16f2db45fa3..727a12b2d5c 100644 --- a/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp +++ b/src/coreComponents/mesh/generators/VTKSuperCellPartitioning.hpp @@ -20,6 +20,7 @@ #include "common/MpiWrapper.hpp" #include "common/TimingMacros.hpp" #include "mesh/generators/VTKMeshGeneratorTools.hpp" +#include "mesh/generators/VTKMeshScattering.hpp" #include "mesh/generators/ParMETISInterface.hpp" #include "LvArray/src/ArrayOfArrays.hpp" @@ -33,15 +34,6 @@ namespace geos namespace vtk { -/** - * @brief Initial distribution strategy for super-cells - */ -enum class InitialDistributionStrategy -{ - MORTON, ///< Morton Z-curve ordering for spatial locality - BLOCK ///< Contiguous block distribution -}; - /** * @brief Super-cell metadata for constrained partitioning * @@ -92,17 +84,18 @@ SuperCellInfo reconstructSuperCellInfo( vtkSmartPointer< vtkUnstructuredGrid > m * @brief Initial redistribution preserving super-cell integrity * * Distributes super-cells across ranks without graph partitioning. - * Faster than ParMETIS for initial scatter. * - * @param cells3D Input mesh (non-empty on rank 0 only) - * @param comm MPI communicator - * @param strategy Distribution strategy (MORTON or BLOCK) + * @param cells3D Input mesh (non-empty on rank 0 only) + * @param comm MPI communicator + * @param scatterMethod Partitioning method. @c kdtree is rejected because it cannot preserve atomicity + * @param cartesianPartitions {nx, ny, nz}, only used for @c cartesian. * @return Redistributed mesh with preserved SuperCellId array */ vtkSmartPointer< vtkDataSet > redistributeBySuperCellBlocks( vtkSmartPointer< vtkUnstructuredGrid > cells3D, MPI_Comm comm, - InitialDistributionStrategy strategy = InitialDistributionStrategy::MORTON ); + ScatterMethod scatterMethod, + arrayView1d< integer const > cartesianPartitions ); /** * @brief Build super-cell adjacency graph diff --git a/src/coreComponents/mesh/generators/VTKUtilities.cpp b/src/coreComponents/mesh/generators/VTKUtilities.cpp index 78f6b905404..4ef60a1ac3a 100644 --- a/src/coreComponents/mesh/generators/VTKUtilities.cpp +++ b/src/coreComponents/mesh/generators/VTKUtilities.cpp @@ -49,7 +49,6 @@ #include #include #include -#include #include #include #include @@ -845,113 +844,6 @@ redistributeBySuperCellGraph( } -/** - * @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 @@ -1779,45 +1671,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 ) @@ -2034,17 +1887,20 @@ 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 partitionFractureWeight, 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 ) { @@ -2181,12 +2037,11 @@ redistributeMeshes( integer const logLevel, { if( hasSuperCells ) { - redistributed3D = redistributeBySuperCellBlocks( cells3D, comm ); + redistributed3D = redistributeBySuperCellBlocks( cells3D, comm, scatterMethod, partitions ); } else { - GEOS_LOG_RANK_0( "Initial redistribution using KD-tree..." ); - redistributed3D = redistributeByKdTree( *cells3D ); + redistributed3D = scatterMesh( scatterMethod, *cells3D, partitions, comm ); if( MpiWrapper::min( redistributed3D->GetNumberOfCells(), comm ) == 0 ) { diff --git a/src/coreComponents/mesh/generators/VTKUtilities.hpp b/src/coreComponents/mesh/generators/VTKUtilities.hpp index a0181fbfd8b..de963f0467d 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 @@ -151,12 +152,13 @@ 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] partitionFractureWeight additional weight to fracture-connected super-cells during partitioning * @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 @@ -164,12 +166,13 @@ 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 partitionFractureWeight, 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..df8e74b7d86 --- /dev/null +++ b/src/coreComponents/mesh/unitTests/testVTKMeshScattering.cpp @@ -0,0 +1,371 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "mesh/generators/VTKSuperCellPartitioning.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 ); +} + + +/// Super-cell atomicity: every super-cell must end up on a single rank, regardless of method. +/// Tags pairs of adjacent cells with the same SuperCellId, so the 64-cell grid becomes +/// 32 atoms of 2 cells each. After redistribution, each SuperCellId must be owned by exactly +/// one rank. +TEST_F( VTKMeshScatteringTest, SuperCellAtomicity ) +{ + static constexpr vtkIdType cellsPerSuperCell = 2; + static constexpr vtkIdType numSuperCells = totalCells / cellsPerSuperCell; + + if( rank == 0 ) + { + vtkNew< vtkIdTypeArray > scIds; + scIds->SetName( "SuperCellId" ); + scIds->SetNumberOfTuples( totalCells ); + for( vtkIdType c = 0; c < totalCells; ++c ) + { + scIds->SetValue( c, c / cellsPerSuperCell ); + } + mesh->GetCellData()->AddArray( scIds ); + } + + stdVector< ScatterMethod > const methods = { + ScatterMethod::contiguous, + ScatterMethod::cartesian, + ScatterMethod::rcb + }; + + for( auto method : methods ) + { + vtkSmartPointer< vtkDataSet > result = + redistributeBySuperCellBlocks( mesh, comm, method, parts.toViewConst() ); + ASSERT_NE( result, nullptr ); + + // Cell conservation + 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 ); + + // Per-rank ownership: ownership[s] = 1 if this rank holds any cell of SuperCellId s. + array1d< integer > ownership( numSuperCells ); + vtkIdTypeArray * scArr = + vtkIdTypeArray::SafeDownCast( result->GetCellData()->GetArray( "SuperCellId" ) ); + ASSERT_NE( scArr, nullptr ) + << "SuperCellId array lost during redistribution for method " << toString( method ); + + for( vtkIdType c = 0; c < localCells; ++c ) + { + vtkIdType const s = scArr->GetValue( c ); + ASSERT_GE( s, 0 ); + ASSERT_LT( s, numSuperCells ); + ownership[s] = 1; + } + + array1d< integer > totalOwners( numSuperCells ); + MpiWrapper::allReduce( ownership, totalOwners, MpiWrapper::Reduction::Sum, comm ); + + for( vtkIdType s = 0; s < numSuperCells; ++s ) + { + EXPECT_EQ( totalOwners[s], 1 ) + << "Super-cell " << s << " ended up on " << totalOwners[s] + << " ranks (expected exactly 1) for method " << toString( method ); + } + } +} + + + +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; +}