From: Rene Gassmoeller Date: Fri, 13 Oct 2017 21:52:43 +0000 (-0600) Subject: Add first particle handler test. Simplify update X-Git-Tag: v9.0.0-rc1~763^2~9 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3fade51fa288f98b732bbc3f7e57b239ef1bcec7;p=dealii.git Add first particle handler test. Simplify update --- diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index 098ceeea37..094e3b8798 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -97,6 +97,16 @@ namespace Particles */ void clear_particles(); + /** + * Updates all internally cached numbers. Note that all functions that + * modify internal data structures and act on multiple particles will + * call this function automatically (e.g. insert_particles), while + * functions that act on single particles will not call this function + * (e.g. insert_particle). This is done because the update is + * expensive compared to single operations. + */ + void update_cache(); + /** * Return an iterator to the first particle. */ @@ -352,32 +362,6 @@ namespace Particles */ unsigned int data_offset; - /** - * Calculates the number of particles in the global model domain. - */ - void - update_n_global_particles(); - - /** - * Calculates and stores the number of particles in the cell that - * contains the most particles in the global model (stored in the - * member variable global_max_particles_per_cell). This variable is a - * state variable, because it is needed to serialize and deserialize - * the particle data correctly in parallel (it determines the size of - * the data chunks per cell that are stored and read). Before accessing - * the variable this function has to be called, unless the state was - * read from another source (e.g. after resuming from a checkpoint). - */ - void - update_global_max_particles_per_cell(); - - /** - * Calculates the next free particle index in the global model domain. - * This equals one plus the highest particle index currently active. - */ - void - update_next_free_particle_index(); - /** * Transfer particles that have crossed subdomain boundaries to other * processors. diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 99a9f72267..13839888cc 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -103,6 +103,36 @@ namespace Particles + template + void + ParticleHandler::update_cache() + { + + types::particle_index locally_highest_index = 0; + unsigned int local_max_particles_per_cell = 0; + unsigned int current_particles_per_cell = 0; + typename Triangulation::active_cell_iterator current_cell; + + for (particle_iterator particle = begin(); particle != end(); ++particle) + { + locally_highest_index = std::max(locally_highest_index,particle->get_id()); + + if (particle->get_surrounding_cell(*triangulation) != current_cell) + current_particles_per_cell = 0; + + ++current_particles_per_cell; + local_max_particles_per_cell = std::max(local_max_particles_per_cell, + current_particles_per_cell); + } + + global_number_of_particles = dealii::Utilities::MPI::sum (particles.size(), triangulation->get_communicator()); + next_free_particle_index = dealii::Utilities::MPI::max (locally_highest_index, triangulation->get_communicator()) + 1; + global_max_particles_per_cell = dealii::Utilities::MPI::max(local_max_particles_per_cell,triangulation->get_communicator()); + } + + + + template typename ParticleHandler::particle_iterator ParticleHandler::begin() const @@ -202,6 +232,7 @@ namespace Particles ParticleHandler::insert_particles(const std::multimap > &new_particles) { particles.insert(new_particles.begin(),new_particles.end()); + update_cache(); } @@ -233,15 +264,6 @@ namespace Particles - template - void - ParticleHandler::update_n_global_particles() - { - global_number_of_particles = dealii::Utilities::MPI::sum (particles.size(), triangulation->get_communicator()); - } - - - template unsigned int ParticleHandler::n_particles_in_cell(const typename Triangulation::active_cell_iterator &cell) const @@ -260,37 +282,6 @@ namespace Particles - template - void - ParticleHandler::update_next_free_particle_index() - { - types::particle_index locally_highest_index = 0; - for (particle_iterator particle = begin(); particle != end(); ++particle) - locally_highest_index = std::max(locally_highest_index,particle->get_id()); - - next_free_particle_index = dealii::Utilities::MPI::max (locally_highest_index, triangulation->get_communicator()) + 1; - } - - - - template - void - ParticleHandler::update_global_max_particles_per_cell() - { - unsigned int local_max_particles_per_cell(0); - typename Triangulation::active_cell_iterator cell = triangulation->begin_active(); - for (; cell!=triangulation->end(); ++cell) - if (cell->is_locally_owned()) - { - local_max_particles_per_cell = std::max(local_max_particles_per_cell, - n_particles_in_cell(cell)); - } - - global_max_particles_per_cell = dealii::Utilities::MPI::max(local_max_particles_per_cell,triangulation->get_communicator()); - } - - - template PropertyPool & ParticleHandler::get_property_pool () const @@ -535,6 +526,7 @@ namespace Particles remove_particle(particles_out_of_cell[i]); particles.insert(sorted_particles_map.begin(),sorted_particles_map.end()); + update_cache(); } @@ -764,7 +756,7 @@ namespace Particles // Only save and load particles if there are any, we might get here for // example if somebody created a ParticleHandler but generated 0 particles. - update_global_max_particles_per_cell(); + update_cache(); if (global_max_particles_per_cell > 0) { @@ -864,7 +856,7 @@ namespace Particles // Reset offset and update global number of particles. The number // can change because of discarded or newly generated particles data_offset = numbers::invalid_unsigned_int; - update_n_global_particles(); + update_cache(); } } diff --git a/tests/particles/particle_02.cc b/tests/particles/particle_02.cc index e4601726a4..922d9bb7be 100644 --- a/tests/particles/particle_02.cc +++ b/tests/particles/particle_02.cc @@ -33,7 +33,7 @@ void test () reference_position(0) = 0.2; reference_position(1) = 0.4; - const Particles::types::particle_index index(7); + const types::particle_index index(7); Particles::Particle<2> particle(position,reference_position,index); diff --git a/tests/particles/particle_03.cc b/tests/particles/particle_03.cc index aaa6a882bf..86d7686866 100644 --- a/tests/particles/particle_03.cc +++ b/tests/particles/particle_03.cc @@ -37,7 +37,7 @@ void test () reference_position(0) = 0.2; reference_position(1) = 0.4; - const Particles::types::particle_index index(7); + const types::particle_index index(7); std::vector properties = {0.15,0.45,0.75}; diff --git a/tests/particles/particle_handler_01.cc b/tests/particles/particle_handler_01.cc new file mode 100644 index 0000000000..24557329e5 --- /dev/null +++ b/tests/particles/particle_handler_01.cc @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + + +// check the creation and destruction of particle within the particle handler class. + +#include "../tests.h" +#include + +#include +#include +#include +#include + +template +void test () +{ + { + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tr); + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr,mapping); + + Particles::Particle particle; + + Point position; + position(0) = 0.3; + if (spacedim>1) + position(1) = 0.5; + if (spacedim>2) + position(2) = 0.7; + + Point reference_position; + reference_position(0) = 0.2; + if (dim>1) + reference_position(1) = 0.4; + if (dim>2) + reference_position(1) = 0.6; + + particle.set_location(position); + particle.set_reference_location(reference_position); + deallog << "Particle location: " << particle.get_location() << std::endl; + + + std::pair::active_cell_iterator, Point > cell_position = + GridTools::find_active_cell_around_point(mapping, tr, particle.get_location()); + + particle_handler.insert_particle(particle,cell_position.first); + particle_handler.update_cache(); + + deallog << "Particle number: " << particle_handler.n_global_particles() << std::endl; + + for (auto particle = particle_handler.begin(); particle != particle_handler.end(); ++particle) + { + deallog << "Particle location: " << particle->get_location() << std::endl; + deallog << "Particle reference location: " << particle->get_reference_location() << std::endl; + } + } + + deallog << "OK" << std::endl; +} + + + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + + initlog(); + + test<2,2>(); + test<2,3>(); + + test<3,3>(); +} diff --git a/tests/particles/particle_handler_01.output b/tests/particles/particle_handler_01.output new file mode 100644 index 0000000000..2cf7a27522 --- /dev/null +++ b/tests/particles/particle_handler_01.output @@ -0,0 +1,16 @@ + +DEAL::Particle location: 0.300000 0.500000 +DEAL::Particle number: 1 +DEAL::Particle location: 0.300000 0.500000 +DEAL::Particle reference location: 0.200000 0.400000 +DEAL::OK +DEAL::Particle location: 0.300000 0.500000 0.700000 +DEAL::Particle number: 1 +DEAL::Particle location: 0.300000 0.500000 0.700000 +DEAL::Particle reference location: 0.200000 0.400000 +DEAL::OK +DEAL::Particle location: 0.300000 0.500000 0.700000 +DEAL::Particle number: 1 +DEAL::Particle location: 0.300000 0.500000 0.700000 +DEAL::Particle reference location: 0.200000 0.600000 0.00000 +DEAL::OK diff --git a/tests/particles/particle_iterator_01.cc b/tests/particles/particle_iterator_01.cc index 67ccd4597f..2eac94b1c6 100644 --- a/tests/particles/particle_iterator_01.cc +++ b/tests/particles/particle_iterator_01.cc @@ -44,7 +44,7 @@ void test () if (dim>1) reference_position(1) = 0.4; - const Particles::types::particle_index index(7); + const types::particle_index index(7); std::vector properties = {0.15,0.45,0.75}; @@ -52,9 +52,9 @@ void test () particle.set_property_pool(pool); particle.set_properties(ArrayView(&properties[0],properties.size())); - std::multimap > map; + std::multimap > map; - Particles::types::LevelInd level_index = std::make_pair(0,0); + types::LevelInd level_index = std::make_pair(0,0); map.insert(std::make_pair(level_index,particle)); particle.get_properties()[0] = 0.05;