AMReXParticlesReader: handle inexact division of grids among MPI ranks

- this commit fixes the grid-process assignment logic robust
when the number of grids cannot be exactly divided among
the number of MPI ranks (pieces in this code)
- this commit also adds MPI awareness to an existing unit test.
- the unit test TestAMReXParticlesReader is run with MPI, when enabled,
on 3 nodes. It passes only when the sum of the number of points across
all nodes is equal to a known value.
- without the change in vtkAMReXParticlesReader, the unit test will faill with MPI

(cherry picked from commit de94dced)

Documentation/release/dev/fix-AMReXParticlesReader-in-MPI.md

+## Fix vtkAMReXParticlesReader when using MPI
+
+The `vtkAMReXParticlesReader` can now function correctly in parallel mode, with MPI.
+Earlier, a bug caused incorrect output data when the number of grids were not
+exactly divided by the number of MPI processes.

IO/AMR/Testing/Cxx/CMakeLists.txt

@@ -6,6 +6,15 @@ ExternalData_Expand_Arguments(VTKData _
   "DATA{${_vtk_build_TEST_INPUT_DATA_DIRECTORY}/Data/AMReX/NonZeroOrigin/,RECURSE:,REGEX:.*}"
   )
 
+if (TARGET VTK::ParallelMPI)
+  set(TestAMReXParticlesReader_NUMPROCS 3)
+  vtk_add_test_mpi(vtkIOAMRCxxTests-MPI tests
+    TESTING_DATA
+    TestAMReXParticlesReader.cxx,NO_VALID,NO_OUTPUT
+    )
+  vtk_test_cxx_executable(vtkIOAMRCxxTests-MPI tests)
+endif()
+
 vtk_add_test_cxx(vtkIOAMRCxxTests tests
   NO_VALID NO_OUTPUT
   TestAMRFlashReader.cxx

IO/AMR/Testing/Cxx/TestAMReXParticlesReader.cxx

 // SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 // SPDX-License-Identifier: BSD-3-Clause
+
+#include <cstdlib>
+#if VTK_MODULE_ENABLE_VTK_ParallelMPI
+#include "vtkMPIController.h"
+#else
+#include "vtkDummyController.h"
+#endif
+
 #include "vtkAMReXParticlesReader.h"
 #include "vtkDataArraySelection.h"
-#include "vtkIdTypeArray.h"
+#include "vtkLogger.h"
 #include "vtkMultiBlockDataSet.h"
 #include "vtkMultiPieceDataSet.h"
 #include "vtkNew.h"
@@ -16,24 +24,40 @@
     if (!(x))                                                                                      \
     {                                                                                              \
       cerr << "FAILED: " << msg << endl;                                                           \
+      controller->Finalize();                                                                      \
       return EXIT_FAILURE;                                                                         \
     }                                                                                              \
   } while (false)
 
 int Validate(vtkMultiBlockDataSet* mb)
 {
+  auto* controller = vtkMultiProcessController::GetGlobalController();
   ensure(mb != nullptr, "expecting vtkMultiBlockDataSet.");
   ensure(mb->GetNumberOfBlocks() == 1, "expecting num-blocks == num-levels == 1");
 
   auto mp = vtkMultiPieceDataSet::SafeDownCast(mb->GetBlock(0));
   ensure(mp != nullptr, "expecting level is maintained in a vtkMultiPieceDataSet.");
   ensure(mp->GetNumberOfPieces() == 8, "expecting 8 datasets in level 0");
-  for (int cc = 0; cc < 8; ++cc)
+  vtkIdType numberOfPointsPerProcess = 0;
+  for (unsigned int cc = 0; cc < mp->GetNumberOfPieces(); ++cc)
   {
-    auto pd = vtkPolyData::SafeDownCast(mp->GetPiece(cc));
-    ensure(pd != nullptr, "expecting polydata for index " << cc);
-    ensure(pd->GetNumberOfPoints() > 0, "expecting non-null points.");
-    ensure(pd->GetPointData()->GetArray("density") != nullptr, "missing density");
+    if (auto pd = vtkPolyData::SafeDownCast(mp->GetPiece(cc)))
+    {
+      ensure(pd != nullptr, "expecting polydata for index " << cc);
+      numberOfPointsPerProcess += pd->GetNumberOfPoints();
+      ensure(numberOfPointsPerProcess > 0, "expecting non-null points.");
+      ensure(pd->GetPointData()->GetArray("density") != nullptr, "missing density");
+    }
+  }
+  vtkIdType totalNumberOfPoints = 0;
+  controller->AllReduce(
+    &numberOfPointsPerProcess, &totalNumberOfPoints, 1, vtkCommunicator::SUM_OP);
+  if (totalNumberOfPoints != 9776)
+  {
+    vtkLog(ERROR, << "# points per process: " << numberOfPointsPerProcess);
+    vtkLog(ERROR, << "Expected total # points: 9776");
+    vtkLog(ERROR, << "Got total # points: " << totalNumberOfPoints);
+    return EXIT_FAILURE;
   }
 
   return EXIT_SUCCESS;
@@ -41,6 +65,16 @@ int Validate(vtkMultiBlockDataSet* mb)
 
 int TestAMReXParticlesReader(int argc, char* argv[])
 {
+#if VTK_MODULE_ENABLE_VTK_ParallelMPI
+  vtkNew<vtkMPIController> controller;
+#else
+  vtkNew<vtkDummyController> controller;
+#endif
+  controller->Initialize(&argc, &argv);
+  const int processId = controller->GetLocalProcessId();
+  const int numberOfProcesses = controller->GetNumberOfProcesses();
+  vtkLogger::SetThreadName("processId=" + std::to_string(processId));
+  vtkMultiProcessController::SetGlobalController(controller);
   // Test 3D
   {
     char* fname = vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/AMReX/MFIX-Exa/plt00000");
@@ -53,9 +87,10 @@ int TestAMReXParticlesReader(int argc, char* argv[])
     reader->UpdateInformation();
     ensure(reader->GetPointDataArraySelection()->ArrayIsEnabled("proc") == 0,
       "`proc` should be disabled.");
-    reader->Update();
+    reader->UpdatePiece(processId, numberOfProcesses, 0);
     if (Validate(reader->GetOutput()) == EXIT_FAILURE)
     {
+      controller->Finalize();
       return EXIT_FAILURE;
     }
   }
@@ -76,5 +111,6 @@ int TestAMReXParticlesReader(int argc, char* argv[])
     ensure(bds[4] == bds[5], "expecting 2D dataset");
   }
 
+  controller->Finalize();
   return EXIT_SUCCESS;
 }

IO/AMR/vtk.module

@@ -24,3 +24,6 @@ TEST_DEPENDS
   VTK::IOXML
   VTK::TestingCore
   VTK::TestingRendering
+TEST_OPTIONAL_DEPENDS
+  VTK::ParallelMPI
+  VTK::ParallelCore

IO/AMR/vtkAMReXParticlesReader.cxx

@@ -831,7 +831,8 @@ bool vtkAMReXParticlesReader::ReadLevel(
   const int quotient = num_grids / num_pieces;
   const int remainder = num_grids % num_pieces;
 
-  const int start_grid_idx = (piece_idx * quotient) + ((piece_idx < remainder) ? 1 : 0);
+  const int start_grid_idx =
+    (piece_idx * quotient) + ((piece_idx < remainder) ? piece_idx : remainder);
   const int grids_count = quotient + ((piece_idx < remainder) ? 1 : 0);
 
   levelDS->SetNumberOfPieces(num_grids);