From 00ae7a76147bb27375c94f02f0e90241288c2ff8 Mon Sep 17 00:00:00 2001 From: Bruno Blais Date: Tue, 26 Apr 2022 22:33:07 -0400 Subject: [PATCH] Add simple test that breaks property pool in debug --- tests/particles/particle_handler_22.cc | 256 ++++++++++++++++++ ...handler_22.with_p4est=true.mpirun=2.output | 5 + 2 files changed, 261 insertions(+) create mode 100644 tests/particles/particle_handler_22.cc create mode 100644 tests/particles/particle_handler_22.with_p4est=true.mpirun=2.output diff --git a/tests/particles/particle_handler_22.cc b/tests/particles/particle_handler_22.cc new file mode 100644 index 0000000000..5274f792dd --- /dev/null +++ b/tests/particles/particle_handler_22.cc @@ -0,0 +1,256 @@ +/* --------------------------------------------------------------------- + * + * 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 in the center + // of the domain to generate ghost particles. We generate enough particle to + // reproduce the bug. + + Point center_of_particles; + center_of_particles[0] = 0.0; + center_of_particles[1] = 0.0; + if (dim == 3) + center_of_particles[2] = 0.0; + + const double outer_radius = 0.5; + 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(); + particle_handler.exchange_ghost_particles(false); + + unsigned int n_part_after_sort = particle_handler.n_global_particles(); + } +} + +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_22.with_p4est=true.mpirun=2.output b/tests/particles/particle_handler_22.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..4c59b55710 --- /dev/null +++ b/tests/particles/particle_handler_22.with_p4est=true.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL:0::Number of particles inserted: 48 + +DEAL:1::Number of particles inserted: 48 + -- 2.39.5