From: Rene Gassmoeller Date: Mon, 16 Sep 2019 21:23:34 +0000 (-0700) Subject: Add random particle generation X-Git-Tag: v9.2.0-rc1~1098^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=86e06890c0fe954e32c61d5c51464a05c54d504e;p=dealii.git Add random particle generation --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 1fee0117ab..50f8639f93 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -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} +} diff --git a/include/deal.II/particles/generators.h b/include/deal.II/particles/generators.h index 77018c454e..272967c781 100644 --- a/include/deal.II/particles/generators.h +++ b/include/deal.II/particles/generators.h @@ -16,13 +16,18 @@ #ifndef dealii_particles_particle_generator_h #define dealii_particles_particle_generator_h +#include + #include #include #include +#include #include +#include + DEAL_II_NAMESPACE_OPEN #ifdef DEAL_II_WITH_P4EST @@ -63,6 +68,100 @@ namespace Particles const Mapping & mapping = StaticMappingQ1::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 + Particle + random_particle_in_cell( + const typename Triangulation::active_cell_iterator &cell, + const types::particle_index id, + std::mt19937 & random_number_generator, + const Mapping &mapping = + StaticMappingQ1::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 + void + probabilistic_locations( + const Triangulation &triangulation, + const Function & probability_density_function, + const bool random_cell_selection, + const types::particle_index n_particles_to_create, + ParticleHandler & particle_handler, + const Mapping & mapping = + StaticMappingQ1::mapping, + const unsigned int random_number_seed = 5432); + } // namespace Generators } // namespace Particles diff --git a/source/particles/generators.cc b/source/particles/generators.cc index cda50b6aa2..45ff913e44 100644 --- a/source/particles/generators.cc +++ b/source/particles/generators.cc @@ -13,8 +13,14 @@ // // --------------------------------------------------------------------- +#include +#include +#include #include +#include +#include + #include DEAL_II_NAMESPACE_OPEN @@ -25,6 +31,79 @@ namespace Particles { namespace Generators { + namespace + { + /** + * Exception + */ + template + DeclException1( + ProbabilityFunctionNegative, + Point, + << "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 + std::vector + compute_local_cumulative_cell_weights( + const Triangulation &triangulation, + const Mapping & mapping, + const Function & probability_density_function) + { + std::vector cumulative_cell_weights( + triangulation.n_active_cells()); + double cumulative_weight = 0.0; + + // Evaluate function at all cell midpoints + const QMidpoint 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 alibi_finite_element; + FEValues 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( + 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 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 *>( - &triangulation)) + if (const auto tria = + dynamic_cast *>( + &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 + Particle + random_particle_in_cell( + const typename Triangulation::active_cell_iterator &cell, + const types::particle_index id, + std::mt19937 & random_number_generator, + const Mapping &mapping) + { + // Uniform distribution on the interval [0,1]. This + // will be used to generate random particle locations. + std::uniform_real_distribution uniform_distribution_01(0, 1); + + const BoundingBox cell_bounding_box(cell->bounding_box()); + const std::pair, Point> 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 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 p_unit = + mapping.transform_real_to_unit_cell(cell, particle_position); + if (GeometryInfo::is_inside_unit_cell(p_unit)) + { + // Generate the particle + return Particle(particle_position, p_unit, id); + } + } + catch (typename Mapping::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(); + } + + + + template + void + probabilistic_locations( + const Triangulation &triangulation, + const Function & probability_density_function, + const bool random_cell_selection, + const types::particle_index n_particles_to_create, + ParticleHandler & particle_handler, + const Mapping & mapping, + const unsigned int random_number_seed) + { + unsigned int combined_seed = random_number_seed; + if (const auto tria = + dynamic_cast *>( + &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 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 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 *>( + &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::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 *>( + &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(n_particles_to_create) * + local_start_weight / global_weight_integral); + n_local_particles = + std::llround(static_cast(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 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::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(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::active_cell_iterator, + Particle> + 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 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 diff --git a/source/particles/generators.inst.in b/source/particles/generators.inst.in index e347d022f3..bef4d7dc82 100644 --- a/source/particles/generators.inst.in +++ b/source/particles/generators.inst.in @@ -30,6 +30,27 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) ParticleHandler &particle_handler, const Mapping &mapping); + + template Particle + 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 &mapping); + + template void + probabilistic_locations( + const Triangulation + & triangulation, + const Function &probability_density_function, + const bool random_cell_selection, + const types::particle_index n_particles_to_create, + ParticleHandler + &particle_handler, + const Mapping &mapping, + const unsigned int random_number_seed); \} \} #endif