]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add random particle generation
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Mon, 16 Sep 2019 21:23:34 +0000 (14:23 -0700)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Mon, 16 Sep 2019 21:23:34 +0000 (14:23 -0700)
doc/doxygen/references.bib
include/deal.II/particles/generators.h
source/particles/generators.cc
source/particles/generators.inst.in

index 1fee0117abeab0b226c29a627cb4dd0c95b64a10..50f8639f93998ad4da9740803734d982b4105d67 100644 (file)
@@ -358,3 +358,13 @@ MRREVIEWER = {Jose Luis Gracia},
    publisher = {SIAM},
    doi       = {10.1137/S0036142902404923}
 }
+
+@article{GLHPW2018,
+  title={Flexible and Scalable Particle-in-Cell Methods With Adaptive Mesh Refinement for Geodynamic Computations},
+  author={Gassm{\"o}ller, Rene and Lokavarapu, Harsha and Heien, Eric and Puckett, Elbridge Gerry and Bangerth, Wolfgang},
+  journal={Geochemistry, Geophysics, Geosystems},
+  volume={19},
+  number={9},
+  pages={3596--3604},
+  year={2018}
+}
index 77018c454e9454f7f0c971acdeefff4ebb135cc4..272967c78111cea623bf3c8ff545e24318c090dd 100644 (file)
 #ifndef dealii_particles_particle_generator_h
 #define dealii_particles_particle_generator_h
 
+#include <deal.II/base/function.h>
+
 #include <deal.II/distributed/tria.h>
 
 #include <deal.II/fe/mapping.h>
 #include <deal.II/fe/mapping_q1.h>
 
+#include <deal.II/particles/particle.h>
 #include <deal.II/particles/particle_handler.h>
 
+#include <random>
+
 DEAL_II_NAMESPACE_OPEN
 
 #ifdef DEAL_II_WITH_P4EST
@@ -63,6 +68,100 @@ namespace Particles
       const Mapping<dim, spacedim> &      mapping =
         StaticMappingQ1<dim, spacedim>::mapping);
 
+    /**
+     * A function that generates one particle at a random location in cell @p cell and with
+     * index @p id. The function expects a random number generator to avoid the expensive generation
+     * and destruction of a generator for every particle and optionally takes
+     * into account a mapping for the cell. The algorithm implemented in the
+     * function is described in @cite GLHPW2018. In short, the algorithm
+     * generates
+     * random locations within the bounding box of the @p cell. It then inverts the mapping
+     * to check if the generated particle is within the cell itself. This makes
+     * sure the algorithm produces statistically random locations even for
+     * nonlinear mappings and distorted cells. However, if the ratio between
+     * bounding box and cell volume becomes very large -- i.e. the cells become
+     * strongly deformed, for example a pencil shaped cell that lies diagonally
+     * in the domain -- then the algorithm can become very inefficient.
+     * Therefore, it only tries to find a location ni the cell a fixed number of
+     * times before throwing an error message.
+     *
+     * @param[in] cell The cell in which a particle is generated.
+     *
+     * @param[in] id The particle index that will be assigned to the new
+     * particle.
+     *
+     * @param[in,out] random_number_generator A random number generator that
+     * will be used for the creation of th particle.
+     *
+     * @param[in] mapping An optional mapping object that is used to map
+     * reference location in the unit cell to the real cell. If no mapping is
+     * provided a MappingQ1 is assumed.
+     */
+    template <int dim, int spacedim = dim>
+    Particle<dim, spacedim>
+    random_particle_in_cell(
+      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+      const types::particle_index                                        id,
+      std::mt19937 &                random_number_generator,
+      const Mapping<dim, spacedim> &mapping =
+        StaticMappingQ1<dim, spacedim>::mapping);
+
+    /**
+     * A function that generates particles randomly in the domain with a
+     * particle density
+     * according to a provided probability density function @p probability_density_function.
+     * The total number of particles that is added to the @p particle_handler object is
+     * @p n_particles_to_create. An optional @p mapping argument
+     * can be used to map from @p particle_reference_locations to the real particle locations.
+     * The function can compute the number of particles per cell either
+     * deterministically by computing the integral of the probability density
+     * function for each cell and creating
+     * particles accordingly (if option @p random_cell_selection set to false), or it can
+     * select cells randomly based on the probability density function and the
+     * cell size
+     * (if option @p random_cell_selection set to true). In either case the position of
+     * individual particles inside the cell is computed randomly.
+     *
+     * The algorithm implemented in the function is described in @cite
+     * GLHPW2018.
+     *
+     * @param[in] triangulation The triangulation associated with the @p particle_handler.
+     *
+     * @param[in] probability_density_function A function with non-negative
+     * values that determines the probability density of a particle to be
+     * generated in this location. The function does not need to be normalized.
+     *
+     * @param[in] random_cell_selection A bool that determines, how the number
+     * of particles per cell is computed (see the description above).
+     *
+     * @param[in] n_particles_to_create The number of particles that will be
+     * created by this function.
+     *
+     * @param[in,out] particle_handler The particle handler that will take
+     * ownership of the generated particles.
+     *
+     * @param[in] mapping An optional mapping object that is used to map
+     * reference location in the unit cell to the real cells of the
+     * triangulation. If no mapping is provided a MappingQ1 is assumed.
+     *
+     * @param[in] random_number_seed An optional seed that determines the
+     * initial state of the random number generator. Use the same number to get
+     * a reproducible particle distribution, or a changing number (e.g. based on
+     * system time) to generate different particle distributions for each call
+     * to this function.
+     */
+    template <int dim, int spacedim = dim>
+    void
+    probabilistic_locations(
+      const Triangulation<dim, spacedim> &triangulation,
+      const Function<spacedim> &          probability_density_function,
+      const bool                          random_cell_selection,
+      const types::particle_index         n_particles_to_create,
+      ParticleHandler<dim, spacedim> &    particle_handler,
+      const Mapping<dim, spacedim> &      mapping =
+        StaticMappingQ1<dim, spacedim>::mapping,
+      const unsigned int random_number_seed = 5432);
+
   } // namespace Generators
 } // namespace Particles
 
index cda50b6aa25e453b032ace8b8bca84c1bb464b5a..45ff913e44c2008b0b1cad0228f152f1ca662d2d 100644 (file)
 //
 // ---------------------------------------------------------------------
 
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/geometry_info.h>
+#include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/signaling_nan.h>
 
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
 #include <deal.II/particles/generators.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -25,6 +31,79 @@ namespace Particles
 {
   namespace Generators
   {
+    namespace
+    {
+      /**
+       * Exception
+       */
+      template <int dim>
+      DeclException1(
+        ProbabilityFunctionNegative,
+        Point<dim>,
+        << "Your probability density function in the particle generator "
+           "returned a negative value for the following position: "
+        << arg1 << ". Please check your function expression.");
+
+
+      // This function integrates the given probability density
+      // function over all cells of the triangulation. For each
+      // cell it stores the cumulative sum (of all previous
+      // cells including the current cell) in a vector that is then
+      // returned. Therefore the returned vector has as many entries
+      // as active cells, the first entry being the integral over the
+      // first cell, and the last entry the integral over the whole
+      // locally owned domain. Cells that are not locally owned
+      // simply store the same value as the cell before (equivalent
+      // to assuming a probability density function value of 0).
+      template <int dim, int spacedim = dim>
+      std::vector<double>
+      compute_local_cumulative_cell_weights(
+        const Triangulation<dim, spacedim> &triangulation,
+        const Mapping<dim, spacedim> &      mapping,
+        const Function<spacedim> &          probability_density_function)
+      {
+        std::vector<double> cumulative_cell_weights(
+          triangulation.n_active_cells());
+        double cumulative_weight = 0.0;
+
+        // Evaluate function at all cell midpoints
+        const QMidpoint<dim> quadrature_formula;
+
+        // In the simplest case we do not even need a FEValues object, because
+        // using cell->center() and cell->measure() would be equivalent. This
+        // fails however for higher-order mappings.
+        FE_Nothing<dim, spacedim> alibi_finite_element;
+        FEValues<dim, spacedim>   fe_values(mapping,
+                                          alibi_finite_element,
+                                          quadrature_formula,
+                                          update_quadrature_points |
+                                            update_JxW_values);
+
+        // compute the integral weight
+        for (const auto &cell : triangulation.active_cell_iterators())
+          {
+            if (cell->is_locally_owned())
+              {
+                fe_values.reinit(cell);
+                const double quadrature_point_weight =
+                  probability_density_function.value(
+                    fe_values.get_quadrature_points()[0]);
+
+                AssertThrow(quadrature_point_weight >= 0.0,
+                            ProbabilityFunctionNegative<dim>(
+                              quadrature_formula.point(0)));
+
+                // get_cell_weight makes sure to return positive values
+                cumulative_weight += quadrature_point_weight * fe_values.JxW(0);
+              }
+            cumulative_cell_weights[cell->active_cell_index()] =
+              cumulative_weight;
+          }
+
+        return cumulative_cell_weights;
+      }
+    } // namespace
+
     template <int dim, int spacedim>
     void
     regular_reference_locations(
@@ -35,9 +114,9 @@ namespace Particles
     {
       types::particle_index particle_index = 0;
 
-      if (const auto tria = dynamic_cast<
-            const parallel::distributed::Triangulation<dim, spacedim> *>(
-            &triangulation))
+      if (const auto tria =
+            dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+              &triangulation))
         {
           const types::particle_index n_particles_to_generate =
             tria->n_locally_owned_active_cells() *
@@ -75,6 +154,241 @@ namespace Particles
 
       particle_handler.update_cached_numbers();
     }
+
+
+
+    template <int dim, int spacedim = dim>
+    Particle<dim, spacedim>
+    random_particle_in_cell(
+      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+      const types::particle_index                                        id,
+      std::mt19937 &                random_number_generator,
+      const Mapping<dim, spacedim> &mapping)
+    {
+      // Uniform distribution on the interval [0,1]. This
+      // will be used to generate random particle locations.
+      std::uniform_real_distribution<double> uniform_distribution_01(0, 1);
+
+      const BoundingBox<spacedim> cell_bounding_box(cell->bounding_box());
+      const std::pair<Point<spacedim>, Point<spacedim>> cell_bounds(
+        cell_bounding_box.get_boundary_points());
+
+      // Generate random points in these bounds until one is within the cell
+      unsigned int       iteration          = 0;
+      const unsigned int maximum_iterations = 100;
+      Point<spacedim>    particle_position;
+      while (iteration < maximum_iterations)
+        {
+          for (unsigned int d = 0; d < spacedim; ++d)
+            {
+              particle_position[d] =
+                uniform_distribution_01(random_number_generator) *
+                  (cell_bounds.second[d] - cell_bounds.first[d]) +
+                cell_bounds.first[d];
+            }
+          try
+            {
+              const Point<dim> p_unit =
+                mapping.transform_real_to_unit_cell(cell, particle_position);
+              if (GeometryInfo<dim>::is_inside_unit_cell(p_unit))
+                {
+                  // Generate the particle
+                  return Particle<dim, spacedim>(particle_position, p_unit, id);
+                }
+            }
+          catch (typename Mapping<dim>::ExcTransformationFailed &)
+            {
+              // The point is not in this cell. Do nothing, just try again.
+            }
+          ++iteration;
+        }
+      AssertThrow(
+        iteration < maximum_iterations,
+        ExcMessage(
+          "Couldn't generate a particle position within the maximum number of tries. "
+          "The ratio between the bounding box volume in which the particle is "
+          "generated and the actual cell volume is approximately: " +
+          std::to_string(
+            cell->measure() /
+            (cell_bounds.second - cell_bounds.first).norm_square())));
+
+      return Particle<dim, spacedim>();
+    }
+
+
+
+    template <int dim, int spacedim>
+    void
+    probabilistic_locations(
+      const Triangulation<dim, spacedim> &triangulation,
+      const Function<spacedim> &          probability_density_function,
+      const bool                          random_cell_selection,
+      const types::particle_index         n_particles_to_create,
+      ParticleHandler<dim, spacedim> &    particle_handler,
+      const Mapping<dim, spacedim> &      mapping,
+      const unsigned int                  random_number_seed)
+    {
+      unsigned int combined_seed = random_number_seed;
+      if (const auto tria =
+            dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+              &triangulation))
+        {
+          const unsigned int my_rank =
+            Utilities::MPI::this_mpi_process(tria->get_communicator());
+          combined_seed += my_rank;
+        }
+      std::mt19937 random_number_generator(combined_seed);
+
+      types::particle_index start_particle_id(0);
+      types::particle_index n_local_particles(0);
+
+      std::vector<types::particle_index> particles_per_cell(
+        triangulation.n_active_cells(), 0);
+
+      // First determine how many particles to generate in which cell
+      {
+        // Get the local accumulated probabilities for every cell
+        const std::vector<double> cumulative_cell_weights =
+          compute_local_cumulative_cell_weights(triangulation,
+                                                mapping,
+                                                probability_density_function);
+
+        // Sum the local integrals over all nodes
+        double local_weight_integral = cumulative_cell_weights.back();
+        double global_weight_integral;
+
+        if (const auto tria =
+              dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+                &triangulation))
+          {
+            global_weight_integral =
+              Utilities::MPI::sum(local_weight_integral,
+                                  tria->get_communicator());
+          }
+        else
+          {
+            global_weight_integral = local_weight_integral;
+          }
+
+        AssertThrow(global_weight_integral > std::numeric_limits<double>::min(),
+                    ExcMessage(
+                      "The integral of the user prescribed probability "
+                      "density function over the domain equals zero, "
+                      "deal.II has no way to determine the cell of "
+                      "generated particles. Please ensure that the "
+                      "provided function is positive in at least a "
+                      "part of the domain; also check the syntax of "
+                      "the function."));
+
+        // Determine the starting weight of this process, which is the sum of
+        // the weights of all processes with a lower rank
+        double local_start_weight = 0.0;
+
+        if (const auto tria =
+              dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
+                &triangulation))
+          {
+            MPI_Exscan(&local_weight_integral,
+                       &local_start_weight,
+                       1,
+                       MPI_DOUBLE,
+                       MPI_SUM,
+                       tria->get_communicator());
+          }
+
+        // Calculate start id and number of local particles
+        start_particle_id =
+          std::llround(static_cast<double>(n_particles_to_create) *
+                       local_start_weight / global_weight_integral);
+        n_local_particles =
+          std::llround(static_cast<double>(n_particles_to_create) *
+                       local_weight_integral / global_weight_integral);
+
+        if (random_cell_selection)
+          {
+            // Uniform distribution on the interval [0,local_weight_integral).
+            // This will be used to randomly select cells for all local
+            // particles.
+            std::uniform_real_distribution<double> uniform_distribution(
+              0.0, local_weight_integral);
+
+            // Loop over all particles to create locally and pick their cells
+            for (types::particle_index current_particle_index = 0;
+                 current_particle_index < n_local_particles;
+                 ++current_particle_index)
+              {
+                // Draw the random number that determines the cell of the
+                // particle
+                const double random_weight =
+                  uniform_distribution(random_number_generator);
+
+                const std::vector<double>::const_iterator selected_cell =
+                  std::lower_bound(cumulative_cell_weights.begin(),
+                                   cumulative_cell_weights.end(),
+                                   random_weight);
+                const unsigned int cell_index =
+                  std::distance(cumulative_cell_weights.begin(), selected_cell);
+
+                ++particles_per_cell[cell_index];
+              }
+          }
+        else
+          {
+            // Compute number of particles per cell according to the ratio
+            // between their weight and the local weight integral
+            types::particle_index particles_created = 0;
+
+            for (const auto &cell : triangulation.active_cell_iterators())
+              if (cell->is_locally_owned())
+                {
+                  const types::particle_index cumulative_particles_to_create =
+                    std::llround(
+                      static_cast<double>(n_local_particles) *
+                      cumulative_cell_weights[cell->active_cell_index()] /
+                      local_weight_integral);
+
+                  // Compute particles for this cell as difference between
+                  // number of particles that should be created including this
+                  // cell minus the number of particles already created.
+                  particles_per_cell[cell->active_cell_index()] =
+                    cumulative_particles_to_create - particles_created;
+                  particles_created +=
+                    particles_per_cell[cell->active_cell_index()];
+                }
+          }
+      }
+
+      // Now generate as many particles per cell as determined above
+      {
+        unsigned int current_particle_index = start_particle_id;
+
+        std::multimap<
+          typename Triangulation<dim, spacedim>::active_cell_iterator,
+          Particle<dim, spacedim>>
+          particles;
+
+        for (const auto &cell : triangulation.active_cell_iterators())
+          if (cell->is_locally_owned())
+            {
+              for (unsigned int i = 0;
+                   i < particles_per_cell[cell->active_cell_index()];
+                   ++i)
+                {
+                  Particle<dim, spacedim> particle =
+                    random_particle_in_cell(cell,
+                                            current_particle_index,
+                                            random_number_generator,
+                                            mapping);
+                  particles.emplace_hint(particles.end(),
+                                         cell,
+                                         std::move(particle));
+                  ++current_particle_index;
+                }
+            }
+
+        particle_handler.insert_particles(particles);
+      }
+    }
   } // namespace Generators
 } // namespace Particles
 
index e347d022f33837592d55f7fa16b1eb6723941f82..bef4d7dc825c8f69ba846bdcc5faba1454801b0a 100644 (file)
@@ -30,6 +30,27 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
           ParticleHandler<deal_II_dimension, deal_II_space_dimension>
             &particle_handler,
           const Mapping<deal_II_dimension, deal_II_space_dimension> &mapping);
+
+        template Particle<deal_II_dimension, deal_II_space_dimension>
+        random_particle_in_cell(
+          const typename Triangulation<
+            deal_II_dimension,
+            deal_II_space_dimension>::active_cell_iterator &cell,
+          const types::particle_index                       id,
+          std::mt19937 &random_number_generator,
+          const Mapping<deal_II_dimension, deal_II_space_dimension> &mapping);
+
+        template void
+        probabilistic_locations<deal_II_dimension, deal_II_space_dimension>(
+          const Triangulation<deal_II_dimension, deal_II_space_dimension>
+            &                                      triangulation,
+          const Function<deal_II_space_dimension> &probability_density_function,
+          const bool                               random_cell_selection,
+          const types::particle_index              n_particles_to_create,
+          ParticleHandler<deal_II_dimension, deal_II_space_dimension>
+            &particle_handler,
+          const Mapping<deal_II_dimension, deal_II_space_dimension> &mapping,
+          const unsigned int random_number_seed);
       \}
     \}
 #endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.