Ghost particle property pool breakage

source/particles/particle_handler.cc

@@ -287,6 +287,10 @@ namespace Particles
         // into a new one (keeping alive the possibly large vectors with
         // particles on cells) into a new container.
         particle_container sorted_particles;
+
+        // note that this call updates owned_particles_end, so that
+        // particle_container_owned_end() below already points to the
+        // new container
         reset_particle_container(sorted_particles);
 
         // iterate over cells and insert the entries in the new order
@@ -295,9 +299,12 @@ namespace Particles
             if (cells_to_particle_cache[cell->active_cell_index()] !=
                 particles.end())
               {
+                // before we move the sorted_particles into particles
+                // particle_container_ghost_end() still points to the
+                // old particles container. Therefore this condition looks quirky.
                 typename particle_container::iterator insert_position =
                   cell->is_locally_owned() ? particle_container_owned_end() :
-                                             particle_container_ghost_end();
+                                             --sorted_particles.end();
                 typename particle_container::iterator new_entry =
                   sorted_particles.insert(
                     insert_position, typename particle_container::value_type());

tests/particles/particle_handler_22.cc

@@ -0,0 +1,256 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2020 - 2021 by the deal.II authors
+ *
+ * This file is part of the deal.II library.
+ *
+ * The deal.II library is free software; you can use it, redistribute
+ * it, and/or modify it under the terms of the GNU Lesser General
+ * Public License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ * The full text of the license can be found in the file LICENSE.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+
+ * This test case is an extremely simplified version of step-68.
+ * A ball made of 48 particles is generated. This ball is moved down slightly.
+ * The particles remain in the simulation domain and they should not disappear.
+ * At the moment of the creation of this test, a bug in the particle_handler
+ would
+ * make one particle disappear and the number would decrease from 48 to 47 at
+ the
+ * second time step.
+ */
+
+// Include files
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/discrete_time.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/parameter_acceptor.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/particles/data_out.h>
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+class VelocityField : public Function<dim>
+{
+public:
+  VelocityField()
+    : Function<dim>(dim)
+  {}
+
+  virtual void
+  vector_value(const Point<dim> &point, Vector<double> &values) const override;
+};
+
+template <int dim>
+void
+VelocityField<dim>::vector_value(const Point<dim> & /*point*/,
+                                 Vector<double> &values) const
+{
+  values[0] = 0;
+  values[1] = -1;
+  if (dim > 2)
+    values[2] = 0;
+}
+
+template <int dim>
+class ParticleTracking
+{
+public:
+  ParticleTracking();
+  void
+  run();
+
+private:
+  void
+  generate_particles();
+
+  void
+  euler_step_analytical(const double dt);
+
+  MPI_Comm                                  mpi_communicator;
+  parallel::distributed::Triangulation<dim> background_triangulation;
+  Particles::ParticleHandler<dim>           particle_handler;
+
+  MappingQ1<dim>     mapping;
+  VelocityField<dim> velocity;
+};
+
+template <int dim>
+ParticleTracking<dim>::ParticleTracking()
+  : mpi_communicator(MPI_COMM_WORLD)
+  , background_triangulation(mpi_communicator)
+{}
+
+// This function generates the tracer particles and the background
+// triangulation on which these particles evolve.
+
+template <int dim>
+void
+ParticleTracking<dim>::generate_particles()
+{
+  // We create an hyper ball triangulation which we globally refine. The bug
+  // that this test tries to reproduce only occured in curevd geometry.
+  Point<dim> center_of_triangulation;
+  center_of_triangulation[0] = 0.;
+  center_of_triangulation[1] = 0.;
+  if (dim == 3)
+    center_of_triangulation[2] = 0.;
+
+  GridGenerator::hyper_ball(background_triangulation,
+                            center_of_triangulation,
+                            1);
+  background_triangulation.refine_global(3);
+
+  // This initializes the background triangulation where the particles are
+  // living and the number of properties of the particles.
+  particle_handler.initialize(background_triangulation, mapping, 1 + dim);
+
+  // We create a particle triangulation which is solely used to generate
+  // the points which will be used to insert the particles. This
+  // triangulation is a hyper shell which is in the center
+  // of the domain to generate ghost particles. We generate enough particle to
+  // reproduce the bug.
+
+  Point<dim> center_of_particles;
+  center_of_particles[0] = 0.0;
+  center_of_particles[1] = 0.0;
+  if (dim == 3)
+    center_of_particles[2] = 0.0;
+
+  const double outer_radius = 0.5;
+  const double inner_radius = 0.01;
+
+  parallel::distributed::Triangulation<dim> particle_triangulation(
+    mpi_communicator);
+
+  GridGenerator::hyper_shell(
+    particle_triangulation, center_of_particles, inner_radius, outer_radius, 6);
+  particle_triangulation.refine_global(1);
+
+  // We generate the necessary bounding boxes for the particles generator.
+  // These bounding boxes are required to quickly identify in which
+  // process's subdomain the inserted particle lies, and which cell owns it.
+  const auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
+    background_triangulation, IteratorFilters::LocallyOwnedCell());
+  const auto global_bounding_boxes =
+    Utilities::MPI::all_gather(mpi_communicator, my_bounding_box);
+
+  // We generate an empty vector of properties. We will attribute the
+  // properties to the particles once they are generated.
+  std::vector<std::vector<double>> properties(
+    particle_triangulation.n_locally_owned_active_cells(),
+    std::vector<double>(dim + 1, 0.));
+
+  // We generate the particles at the position of a single
+  // point quadrature. Consequently, one particle will be generated
+  // at the centroid of each cell.
+  Particles::Generators::quadrature_points(particle_triangulation,
+                                           QMidpoint<dim>(),
+                                           global_bounding_boxes,
+                                           particle_handler,
+                                           mapping,
+                                           properties);
+
+  deallog << "Number of particles inserted: "
+          << particle_handler.n_global_particles() << std::endl;
+}
+
+// We integrate the particle trajectories using a first order explicit Euler
+// scheme.
+template <int dim>
+void
+ParticleTracking<dim>::euler_step_analytical(const double dt)
+{
+  const unsigned int this_mpi_rank =
+    Utilities::MPI::this_mpi_process(mpi_communicator);
+  Vector<double> particle_velocity(dim);
+
+  // Looping over all particles in the domain using a
+  // particle iterator
+  for (auto &particle : particle_handler)
+    {
+      // We calculate the velocity of the particles using their current
+      // location.
+      Point<dim> particle_location = particle.get_location();
+      velocity.vector_value(particle_location, particle_velocity);
+
+      // This updates the position of the particles and sets the old position
+      // equal to the new position of the particle.
+      for (int d = 0; d < dim; ++d)
+        particle_location[d] += particle_velocity[d] * dt;
+
+      particle.set_location(particle_location);
+
+      // We store the processor id (a scalar) and the particle velocity (a
+      // vector) in the particle properties.
+      ArrayView<double> properties = particle.get_properties();
+      for (int d = 0; d < dim; ++d)
+        properties[d] = particle_velocity[d];
+      properties[dim] = this_mpi_rank;
+    }
+}
+
+template <int dim>
+void
+ParticleTracking<dim>::run()
+{
+  DiscreteTime discrete_time(0, 0.015, 0.005);
+
+  generate_particles();
+
+  // The particles are advected by looping over time.
+  while (!discrete_time.is_at_end())
+    {
+      discrete_time.advance_time();
+      velocity.set_time(discrete_time.get_previous_time());
+
+      euler_step_analytical(discrete_time.get_previous_step_size());
+
+      unsigned int n_part_before_sort = particle_handler.n_global_particles();
+
+      particle_handler.sort_particles_into_subdomains_and_cells();
+      particle_handler.exchange_ghost_particles(false);
+
+      unsigned int n_part_after_sort = particle_handler.n_global_particles();
+    }
+}
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  MPILogInitAll all;
+
+  ParticleTracking<3> particle_advection_problem;
+  particle_advection_problem.run();
+}

tests/particles/particle_handler_22.with_p4est=true.mpirun=2.output

@@ -0,0 +1,5 @@
+
+DEAL:0::Number of particles inserted: 48
+
+DEAL:1::Number of particles inserted: 48
+