From 4b8056566a681ccfdfcc194693b4e7aec7d73dea Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Mon, 30 May 2022 16:56:08 -0400 Subject: [PATCH] Prevent endless particle loop during refinement --- source/particles/particle_handler.cc | 21 ++- ...article_handler_refinement_outside_cell.cc | 158 ++++++++++++++++++ ...tside_cell.with_p4est=true.mpirun=1.output | 16 ++ 3 files changed, 194 insertions(+), 1 deletion(-) create mode 100644 tests/particles/particle_handler_refinement_outside_cell.cc create mode 100644 tests/particles/particle_handler_refinement_outside_cell.with_p4est=true.mpirun=1.output diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 9bd6af4f1d..51798fbc40 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -2394,6 +2394,8 @@ namespace Particles auto particle = loaded_particles_on_cell.begin(); for (unsigned int i = 0; i < cache->particles.size();) { + bool found_new_cell = false; + for (unsigned int child_index = 0; child_index < GeometryInfo::max_children_per_cell; ++child_index) @@ -2407,8 +2409,10 @@ namespace Particles const Point p_unit = mapping->transform_real_to_unit_cell( child, particle->get_location()); - if (GeometryInfo::is_inside_unit_cell(p_unit)) + if (GeometryInfo::is_inside_unit_cell(p_unit, + 1e-12)) { + found_new_cell = true; particle->set_reference_location(p_unit); // if the particle is not in child 0, we stored the @@ -2433,6 +2437,21 @@ namespace Particles catch (typename Mapping::ExcTransformationFailed &) {} } + + if (found_new_cell == false) + { + // If we get here, we did not find the particle in any + // child. This case may happen for particles that are at the + // boundary for strongly curved cells. We apply a tolerance + // in the call to GeometryInfo::is_inside_unit_cell to + // account for this, but if that is not enough, we still + // need to prevent an endless loop here. Delete the particle + // and move on. + signals.particle_lost(particle, + particle->get_surrounding_cell()); + cache->particles[i] = cache->particles.back(); + cache->particles.resize(cache->particles.size() - 1); + } } // clean up in case child 0 has no particle left if (cache->particles.empty()) diff --git a/tests/particles/particle_handler_refinement_outside_cell.cc b/tests/particles/particle_handler_refinement_outside_cell.cc new file mode 100644 index 0000000000..fb2e900733 --- /dev/null +++ b/tests/particles/particle_handler_refinement_outside_cell.cc @@ -0,0 +1,158 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 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. +// +// --------------------------------------------------------------------- + + + +// Check that particles are correctly transferred during grid +// refinement/coarsening. Like particle_handler_05.cc, but with +// a particle that is registered in a parent cell, but will not +// be found in any child cell. Before this fix, this test +// would lead to an endless loop. + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +void +create_regular_particle_distribution( + Particles::ParticleHandler & particle_handler, + const parallel::distributed::Triangulation &tr, + const unsigned int particles_per_direction = 3) +{ + for (unsigned int i = 0; i < particles_per_direction; ++i) + for (unsigned int j = 0; j < particles_per_direction; ++j) + { + Point position; + Point reference_position; + unsigned int id = i * particles_per_direction + j; + + position[0] = static_cast(i) / + static_cast(particles_per_direction - 1); + position[1] = static_cast(j) / + static_cast(particles_per_direction - 1); + + if (dim > 2) + for (unsigned int k = 0; k < particles_per_direction; ++k) + { + position[2] = static_cast(j) / + static_cast(particles_per_direction - 1); + } + else + { + Particles::Particle particle(position, + reference_position, + id); + + typename parallel::distributed::Triangulation:: + active_cell_iterator cell = + GridTools::find_active_cell_around_point( + tr, particle.get_location()); + + particle_handler.insert_particle(particle, cell); + } + } +} + +template +void +test() +{ + { + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tr); + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr, mapping); + + // Create a particle outside of the domain + Point position; + position[0] = 3; + Point reference_position; + reference_position[0] = 0.5; + types::particle_index id = 0; + + particle_handler.insert_particle(position, + reference_position, + id, + tr.begin_active()); + + position[0] = 1.0 + 1e-13; + id = 1; + + particle_handler.insert_particle(position, + reference_position, + id, + tr.begin_active()); + + for (const auto &particle : particle_handler) + deallog << "Before refinement particle id " << particle.get_id() + << " is in cell " << particle.get_surrounding_cell(tr) + << std::endl; + + // Check that all particles are moved to children + particle_handler.prepare_for_coarsening_and_refinement(); + tr.refine_global(1); + particle_handler.unpack_after_coarsening_and_refinement(); + + for (const auto &particle : particle_handler) + deallog << "After refinement particle id " << particle.get_id() + << " is in cell " << particle.get_surrounding_cell(tr) + << std::endl; + + // Reverse the refinement and check again + for (auto cell = tr.begin_active(); cell != tr.end(); ++cell) + cell->set_coarsen_flag(); + + particle_handler.prepare_for_coarsening_and_refinement(); + tr.execute_coarsening_and_refinement(); + particle_handler.unpack_after_coarsening_and_refinement(); + + for (const auto &particle : particle_handler) + deallog << "After coarsening particle id " << particle.get_id() + << " is in cell " << particle.get_surrounding_cell(tr) + << std::endl; + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll all; + + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/particle_handler_refinement_outside_cell.with_p4est=true.mpirun=1.output b/tests/particles/particle_handler_refinement_outside_cell.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..a82262d565 --- /dev/null +++ b/tests/particles/particle_handler_refinement_outside_cell.with_p4est=true.mpirun=1.output @@ -0,0 +1,16 @@ + +DEAL:0:2d/2d::Before refinement particle id 0 is in cell 0.0 +DEAL:0:2d/2d::Before refinement particle id 1 is in cell 0.0 +DEAL:0:2d/2d::After refinement particle id 1 is in cell 1.1 +DEAL:0:2d/2d::After coarsening particle id 1 is in cell 0.0 +DEAL:0:2d/2d::OK +DEAL:0:2d/3d::Before refinement particle id 0 is in cell 0.0 +DEAL:0:2d/3d::Before refinement particle id 1 is in cell 0.0 +DEAL:0:2d/3d::After refinement particle id 1 is in cell 1.1 +DEAL:0:2d/3d::After coarsening particle id 1 is in cell 0.0 +DEAL:0:2d/3d::OK +DEAL:0:3d/3d::Before refinement particle id 0 is in cell 0.0 +DEAL:0:3d/3d::Before refinement particle id 1 is in cell 0.0 +DEAL:0:3d/3d::After refinement particle id 1 is in cell 1.1 +DEAL:0:3d/3d::After coarsening particle id 1 is in cell 0.0 +DEAL:0:3d/3d::OK -- 2.39.5