From: Rene Gassmoeller Date: Wed, 10 Nov 2021 16:52:16 +0000 (-0500) Subject: Optimize particle generation X-Git-Tag: v9.4.0-rc1~837^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cf37c2d0d7b974a69fcfc6b87108e16314dd357d;p=dealii.git Optimize particle generation --- diff --git a/include/deal.II/particles/generators.h b/include/deal.II/particles/generators.h index 603f82b09a..5059bf415d 100644 --- a/include/deal.II/particles/generators.h +++ b/include/deal.II/particles/generators.h @@ -111,6 +111,23 @@ namespace Particles (ReferenceCells::get_hypercube() .template get_default_linear_mapping())); + /** + * A function that generates one particle at a random location in cell @p cell and with + * index @p id. This version of the function above immediately inserts the generated + * particle into the @p particle_handler and returns a iterator to it instead of + * a particle object. This avoids unnecessary copies of the particle. + */ + template + ParticleIterator + random_particle_in_cell_insert( + const typename Triangulation::active_cell_iterator &cell, + const types::particle_index id, + std::mt19937 & random_number_generator, + ParticleHandler &particle_handler, + const Mapping & mapping = + (ReferenceCells::get_hypercube() + .template get_default_linear_mapping())); + /** * A function that generates particles randomly in the domain with a * particle density diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index e5753282cf..ea42cb7618 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -151,6 +151,23 @@ namespace Particles void clear_particles(); + /** + * This function can be used to preemptively reserve memory for particle + * data. Calling this function before inserting particles will reduce + * memory allocations and therefore increase the performance. Calling + * this function is optional; if memory is not already allocated it will + * be allocated automatically during the insertion. It is recommended to + * use this function if you know the number of particles that will be + * inserted, but cannot use one of the collective particle insertion + * functions. + * + * @param n_particles Number of particles to reserve memory for. Note that + * this is the total number of particles to be stored, not the number of + * particles to be newly inserted. + */ + void + reserve(std::size_t n_particles); + /** * Update all internally cached numbers. Note that all functions that * modify internal data structures and act on multiple particles will @@ -270,6 +287,20 @@ namespace Particles const Particle &particle, const typename Triangulation::active_cell_iterator &cell); + /** + * Insert a particle into the collection of particles given all the + * properties necessary for a particle. This function is used internally to + * efficiently generate particles without the detour through a Particle + * object. + */ + particle_iterator + insert_particle( + const Point & position, + const Point & reference_position, + const types::particle_index particle_index, + const typename Triangulation::active_cell_iterator &cell, + const ArrayView &properties = {}); + /** * Insert a number of particles into the collection of particles. * This function involves a copy of the particles and their properties. @@ -864,20 +895,6 @@ namespace Particles const void *& data, const typename Triangulation::active_cell_iterator &cell); - /** - * Insert a particle into the collection of particles given all the - * properties necessary for a particle. This function is used internally to - * efficiently generate particles without the detour through a Particle - * object. - */ - particle_iterator - insert_particle( - const Point & position, - const Point & reference_position, - const types::particle_index particle_index, - const typename Triangulation::active_cell_iterator &cell, - const ArrayView &properties = {}); - /** * Perform the local insertion operation into the particle container. This * function is used in the higher-level functions inserting particles. diff --git a/source/particles/generators.cc b/source/particles/generators.cc index 32750d09fd..c6d96a77ac 100644 --- a/source/particles/generators.cc +++ b/source/particles/generators.cc @@ -106,6 +106,68 @@ namespace Particles return cumulative_cell_weights; } + + + + // This function generates a random position in the given cell and + // returns the position and its coordinates in the unit cell. It first + // tries to generate a random and uniformly distributed point in the + // real space, but if that fails (e.g. because the cell has a bad aspect + // ratio) it reverts to generating a random point in the unit cell. + template + std::pair, Point> + random_location_in_cell( + const typename Triangulation::active_cell_iterator &cell, + const Mapping &mapping, + std::mt19937 & random_number_generator) + { + // 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 + // or we exceed the maximum number of attempts. + const unsigned int n_attempts = 100; + Point position; + Point position_unit; + for (unsigned int i = 0; i < n_attempts; ++i) + { + for (unsigned int d = 0; d < spacedim; ++d) + { + position[d] = uniform_distribution_01(random_number_generator) * + (cell_bounds.second[d] - cell_bounds.first[d]) + + cell_bounds.first[d]; + } + + try + { + position_unit = + mapping.transform_real_to_unit_cell(cell, position); + + if (GeometryInfo::is_inside_unit_cell(position_unit)) + return std::make_pair(position, position_unit); + } + catch (typename Mapping::ExcTransformationFailed &) + { + // The point is not in this cell. Do nothing, just try again. + } + } + + // If the above algorithm has not worked (e.g. because of badly + // deformed cells), retry generating particles + // randomly within the reference cell. This is not generating a + // uniform distribution in real space, but will always succeed. + for (unsigned int d = 0; d < dim; ++d) + position_unit[d] = uniform_distribution_01(random_number_generator); + + position = mapping.transform_unit_to_real_cell(cell, position_unit); + + return std::make_pair(position, position_unit); + } } // namespace template @@ -117,15 +179,16 @@ namespace Particles const Mapping & mapping) { types::particle_index particle_index = 0; + types::particle_index n_particles_to_generate = + triangulation.n_active_cells() * particle_reference_locations.size(); #ifdef DEAL_II_WITH_MPI if (const auto tria = dynamic_cast *>( &triangulation)) { - const types::particle_index n_particles_to_generate = - tria->n_locally_owned_active_cells() * - particle_reference_locations.size(); + n_particles_to_generate = tria->n_locally_owned_active_cells() * + particle_reference_locations.size(); // The local particle start index is the number of all particles // generated on lower MPI ranks. @@ -139,6 +202,9 @@ namespace Particles } #endif + particle_handler.reserve(particle_handler.n_locally_owned_particles() + + n_particles_to_generate); + for (const auto &cell : triangulation.active_cell_iterators()) { if (cell->is_locally_owned()) @@ -150,10 +216,10 @@ namespace Particles mapping.transform_unit_to_real_cell(cell, reference_location); - const Particle particle(position_real, - reference_location, - particle_index); - particle_handler.insert_particle(particle, cell); + particle_handler.insert_particle(position_real, + reference_location, + particle_index, + cell); ++particle_index; } } @@ -172,54 +238,31 @@ namespace Particles 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(); + const auto position_and_reference_position = + random_location_in_cell(cell, mapping, random_number_generator); + return Particle(position_and_reference_position.first, + position_and_reference_position.second, + id); + } + + + + template + ParticleIterator + random_particle_in_cell_insert( + const typename Triangulation::active_cell_iterator &cell, + const types::particle_index id, + std::mt19937 & random_number_generator, + ParticleHandler &particle_handler, + const Mapping & mapping) + { + const auto position_and_reference_position = + random_location_in_cell(cell, mapping, random_number_generator); + return particle_handler.insert_particle( + position_and_reference_position.first, + position_and_reference_position.second, + id, + cell); } @@ -377,13 +420,10 @@ namespace Particles // Now generate as many particles per cell as determined above { + particle_handler.reserve(particle_handler.n_locally_owned_particles() + + n_local_particles); 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()) { @@ -391,19 +431,17 @@ namespace Particles 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)); + random_particle_in_cell_insert(cell, + current_particle_index, + random_number_generator, + particle_handler, + mapping); + ++current_particle_index; } } - particle_handler.insert_particles(particles); + particle_handler.update_cached_numbers(); } } @@ -461,6 +499,8 @@ namespace Particles const std::vector> &particle_reference_locations = quadrature.get_points(); std::vector> points_to_generate; + points_to_generate.reserve(particle_reference_locations.size() * + triangulation.n_active_cells()); // Loop through cells and gather gauss points for (const auto &cell : triangulation.active_cell_iterators()) diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index bcfa91ea94..bdc3b69301 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -221,6 +221,15 @@ namespace Particles + template + void + ParticleHandler::reserve(std::size_t n_particles) + { + property_pool->reserve(n_particles); + } + + + template void ParticleHandler::reset_particle_container( @@ -668,6 +677,7 @@ namespace Particles typename Triangulation::active_cell_iterator, Particle> &new_particles) { + reserve(n_locally_owned_particles() + new_particles.size()); for (const auto &cell_and_particle : new_particles) insert_particle(cell_and_particle.second, cell_and_particle.first); @@ -684,6 +694,7 @@ namespace Particles Assert(triangulation != nullptr, ExcInternalError()); update_cached_numbers(); + reserve(n_locally_owned_particles() + positions.size()); // Determine the starting particle index of this process, which // is the highest currently existing particle index plus the sum