From 44dbc2b30fdeb72d6649789eff42923ec6cd635f Mon Sep 17 00:00:00 2001 From: Bruno Blais Date: Wed, 19 Jan 2022 15:55:00 -0500 Subject: [PATCH] Fix lost particle by recalculating mapping + test --- source/particles/particle_handler.cc | 12 +- tests/particles/particle_handler_21.cc | 231 ++++++++++++++++++ ...particle_handler_21.with_p4est=true.output | 5 + 3 files changed, 245 insertions(+), 3 deletions(-) create mode 100644 tests/particles/particle_handler_21.cc create mode 100644 tests/particles/particle_handler_21.with_p4est=true.output diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 1d7a9b703d..f3b6b74c30 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -1250,11 +1250,17 @@ namespace Particles reference_locations); auto particle = pic.begin(); - for (const auto &p_unit : reference_locations) + for (unsigned int p=0 ; p< reference_locations.size() ; ++p) { + auto & p_unit = reference_locations[p]; if (p_unit[0] == std::numeric_limits::infinity() || - !GeometryInfo::is_inside_unit_cell(p_unit)) - particles_out_of_cell.push_back(particle); + !GeometryInfo::is_inside_unit_cell(reference_locations[p])) + { + // Try to find the reference location again for this particle + p_unit = mapping->transform_real_to_unit_cell(cell,real_locations[p]) ; + if (!GeometryInfo::is_inside_unit_cell(p_unit)) + particles_out_of_cell.push_back(particle); + } else particle->set_reference_location(p_unit); diff --git a/tests/particles/particle_handler_21.cc b/tests/particles/particle_handler_21.cc new file mode 100644 index 0000000000..b2305e72a6 --- /dev/null +++ b/tests/particles/particle_handler_21.cc @@ -0,0 +1,231 @@ +/* --------------------------------------------------------------------- + * + * 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 +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template class VelocityField : public Function { +public: + VelocityField() : Function(dim) {} + + virtual void vector_value(const Point &point, + Vector &values) const override; +}; + +template +void VelocityField::vector_value(const Point & /*point*/, + Vector &values) const { + values[0] = 0; + values[1] = -1; + if (dim > 2) + values[2] = 0; +} + +template class ParticleTracking { +public: + ParticleTracking(); + void run(); + +private: + void generate_particles(); + + void euler_step_analytical(const double dt); + + MPI_Comm mpi_communicator; + parallel::distributed::Triangulation background_triangulation; + Particles::ParticleHandler particle_handler; + + MappingQ1 mapping; + VelocityField velocity; +}; + +template +ParticleTracking::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 void ParticleTracking::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 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 offset from the + // center of the simulation domain. We generate enough particle to reproduce + // the bug. + + Point center_of_particles; + center_of_particles[0] = 0.0; + center_of_particles[1] = -.735; + if (dim == 3) + center_of_particles[2] = 0.0; + + const double outer_radius = 0.2; + const double inner_radius = 0.01; + + parallel::distributed::Triangulation 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> properties( + particle_triangulation.n_locally_owned_active_cells(), + std::vector(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(), 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 +void ParticleTracking::euler_step_analytical(const double dt) { + const unsigned int this_mpi_rank = + Utilities::MPI::this_mpi_process(mpi_communicator); + Vector 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 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 properties = particle.get_properties(); + for (int d = 0; d < dim; ++d) + properties[d] = particle_velocity[d]; + properties[dim] = this_mpi_rank; + } +} + +template void ParticleTracking::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(); + unsigned int n_part_after_sort = particle_handler.n_global_particles(); + + deallog << "Number of particles before sort : " << n_part_before_sort + << " After : " << n_part_after_sort << std::endl; + } +} + +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(); +} diff --git a/tests/particles/particle_handler_21.with_p4est=true.output b/tests/particles/particle_handler_21.with_p4est=true.output new file mode 100644 index 0000000000..19a545fc33 --- /dev/null +++ b/tests/particles/particle_handler_21.with_p4est=true.output @@ -0,0 +1,5 @@ + +DEAL:0::Number of particles inserted: 48 +DEAL:0::Number of particles before sort : 48 After : 48 +DEAL:0::Number of particles before sort : 48 After : 48 +DEAL:0::Number of particles before sort : 48 After : 48 -- 2.39.5