From: Rene Gassmoeller Date: Wed, 15 Jan 2025 10:57:06 +0000 (+0100) Subject: Fix division by zero in particle generator X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b42cc356d3e2e5928b6699ece7123c233482e7c0;p=dealii.git Fix division by zero in particle generator --- diff --git a/doc/news/changes/minor/20250228Gassmoeller b/doc/news/changes/minor/20250228Gassmoeller new file mode 100644 index 0000000000..134bbc7861 --- /dev/null +++ b/doc/news/changes/minor/20250228Gassmoeller @@ -0,0 +1,6 @@ +Fixed: The particle generator function 'probabilistic_locations' could +under certain conditions cause a division by zero if the user-provided +probability density function evaluated to 0 everywhere in a local domain, +but not globally. This has been fixed. +
+(Rene Gassmoeller, 2025/02/28) diff --git a/source/particles/generators.cc b/source/particles/generators.cc index c444540322..7ac588c226 100644 --- a/source/particles/generators.cc +++ b/source/particles/generators.cc @@ -381,26 +381,29 @@ namespace Particles } else { - // Compute number of particles per cell according to the ratio - // between their weight and the local weight integral - types::particle_index particles_created = 0; - - for (const auto &cell : triangulation.active_cell_iterators() | - IteratorFilters::LocallyOwnedCell()) + if (local_weight_integral > 0) { - const types::particle_index cumulative_particles_to_create = - std::llround( - static_cast(n_local_particles) * - cumulative_cell_weights[cell->active_cell_index()] / - local_weight_integral); - - // Compute particles for this cell as difference between - // number of particles that should be created including this - // cell minus the number of particles already created. - particles_per_cell[cell->active_cell_index()] = - cumulative_particles_to_create - particles_created; - particles_created += - particles_per_cell[cell->active_cell_index()]; + types::particle_index particles_created = 0; + + // Compute number of particles per cell according to the ratio + // between their weight and the local weight integral + for (const auto &cell : triangulation.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + const types::particle_index cumulative_particles_to_create = + std::llround( + static_cast(n_local_particles) * + cumulative_cell_weights[cell->active_cell_index()] / + local_weight_integral); + + // Compute particles for this cell as difference between + // number of particles that should be created including this + // cell minus the number of particles already created. + particles_per_cell[cell->active_cell_index()] = + cumulative_particles_to_create - particles_created; + particles_created += + particles_per_cell[cell->active_cell_index()]; + } } } } diff --git a/tests/particles/generators_14.cc b/tests/particles/generators_14.cc new file mode 100644 index 0000000000..9134353b97 --- /dev/null +++ b/tests/particles/generators_14.cc @@ -0,0 +1,111 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2019 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// like generators_05.cc, but in particular +// check that the function works if the integrated probability density +// function is zero on some parallel ranks. + +#include + +#include + +#include + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +test(Function &probability_density_function) +{ + { + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr, mapping); + + Particles::Generators::probabilistic_locations( + tr, probability_density_function, false, 16, particle_handler, mapping); + + deallog << "Rank: " + << std::to_string(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)) + << ". Locally owned active cells: " + << tr.n_locally_owned_active_cells() << std::endl; + + deallog << "Global particles: " << particle_handler.n_global_particles() + << std::endl; + + for (const auto &cell : tr.active_cell_iterators()) + { + if (cell->is_locally_owned()) + { + deallog << "Cell " << cell << " has " + << particle_handler.n_particles_in_cell(cell) + << " particles." << std::endl; + } + } + + for (const auto &particle : particle_handler) + { + deallog << "Particle index " << particle.get_id() << " is in cell " + << particle.get_surrounding_cell() << std::endl; + deallog << "Particle location: " << particle.get_location() + << std::endl; + } + } + + deallog << "OK" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll init; + + { + deallog.push("2d/2d"); + Functions::CutOffFunctionW1<2> probability_density_function(0.5); + test<2, 2>(probability_density_function); + deallog.pop(); + } + { + deallog.push("2d/3d"); + Functions::CutOffFunctionW1<3> probability_density_function(0.5); + test<2, 3>(probability_density_function); + } + { + deallog.pop(); + deallog.push("3d/3d"); + Functions::CutOffFunctionW1<3> probability_density_function(0.5); + test<3, 3>(probability_density_function); + deallog.pop(); + } +} diff --git a/tests/particles/generators_14.with_p4est=true.mpirun=2.output b/tests/particles/generators_14.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..fb9d7174d6 --- /dev/null +++ b/tests/particles/generators_14.with_p4est=true.mpirun=2.output @@ -0,0 +1,213 @@ + +DEAL:0:2d/2d::Rank: 0. Locally owned active cells: 8 +DEAL:0:2d/2d::Global particles: 16 +DEAL:0:2d/2d::Cell 2.0 has 10 particles. +DEAL:0:2d/2d::Cell 2.1 has 3 particles. +DEAL:0:2d/2d::Cell 2.2 has 3 particles. +DEAL:0:2d/2d::Cell 2.3 has 0 particles. +DEAL:0:2d/2d::Cell 2.4 has 0 particles. +DEAL:0:2d/2d::Cell 2.5 has 0 particles. +DEAL:0:2d/2d::Cell 2.6 has 0 particles. +DEAL:0:2d/2d::Cell 2.7 has 0 particles. +DEAL:0:2d/2d::Particle index 0 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.200844 0.125494 +DEAL:0:2d/2d::Particle index 1 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.185603 0.197844 +DEAL:0:2d/2d::Particle index 2 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.128636 0.203122 +DEAL:0:2d/2d::Particle index 3 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.0413383 0.0199530 +DEAL:0:2d/2d::Particle index 4 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.0497937 0.217265 +DEAL:0:2d/2d::Particle index 5 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.182363 0.0966280 +DEAL:0:2d/2d::Particle index 6 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.219574 0.00346588 +DEAL:0:2d/2d::Particle index 7 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.0326305 0.00256411 +DEAL:0:2d/2d::Particle index 8 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.206502 0.174084 +DEAL:0:2d/2d::Particle index 9 is in cell 2.0 +DEAL:0:2d/2d::Particle location: 0.0756403 0.0224136 +DEAL:0:2d/2d::Particle index 10 is in cell 2.1 +DEAL:0:2d/2d::Particle location: 0.470285 0.0513904 +DEAL:0:2d/2d::Particle index 11 is in cell 2.1 +DEAL:0:2d/2d::Particle location: 0.460320 0.104390 +DEAL:0:2d/2d::Particle index 12 is in cell 2.1 +DEAL:0:2d/2d::Particle location: 0.400281 0.221697 +DEAL:0:2d/2d::Particle index 13 is in cell 2.2 +DEAL:0:2d/2d::Particle location: 0.133304 0.312895 +DEAL:0:2d/2d::Particle index 14 is in cell 2.2 +DEAL:0:2d/2d::Particle location: 0.105785 0.418393 +DEAL:0:2d/2d::Particle index 15 is in cell 2.2 +DEAL:0:2d/2d::Particle location: 0.00552836 0.466790 +DEAL:0:2d/2d::OK +DEAL:0:2d/3d::Rank: 0. Locally owned active cells: 8 +DEAL:0:2d/3d::Global particles: 16 +DEAL:0:2d/3d::Cell 2.0 has 10 particles. +DEAL:0:2d/3d::Cell 2.1 has 3 particles. +DEAL:0:2d/3d::Cell 2.2 has 3 particles. +DEAL:0:2d/3d::Cell 2.3 has 0 particles. +DEAL:0:2d/3d::Cell 2.4 has 0 particles. +DEAL:0:2d/3d::Cell 2.5 has 0 particles. +DEAL:0:2d/3d::Cell 2.6 has 0 particles. +DEAL:0:2d/3d::Cell 2.7 has 0 particles. +DEAL:0:2d/3d::Particle index 0 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.200844 0.125494 0.00000 +DEAL:0:2d/3d::Particle index 1 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.197844 0.128636 0.00000 +DEAL:0:2d/3d::Particle index 2 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.0413383 0.0199530 0.00000 +DEAL:0:2d/3d::Particle index 3 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.217265 0.182363 0.00000 +DEAL:0:2d/3d::Particle index 4 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.219574 0.00346588 0.00000 +DEAL:0:2d/3d::Particle index 5 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.00256411 0.206502 0.00000 +DEAL:0:2d/3d::Particle index 6 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.0756403 0.0224136 0.00000 +DEAL:0:2d/3d::Particle index 7 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.0513904 0.210320 0.00000 +DEAL:0:2d/3d::Particle index 8 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.150281 0.221697 0.00000 +DEAL:0:2d/3d::Particle index 9 is in cell 2.0 +DEAL:0:2d/3d::Particle location: 0.0628946 0.105785 0.00000 +DEAL:0:2d/3d::Particle index 10 is in cell 2.1 +DEAL:0:2d/3d::Particle location: 0.255528 0.216790 0.00000 +DEAL:0:2d/3d::Particle index 11 is in cell 2.1 +DEAL:0:2d/3d::Particle location: 0.315964 0.238469 0.00000 +DEAL:0:2d/3d::Particle index 12 is in cell 2.1 +DEAL:0:2d/3d::Particle location: 0.270879 0.204701 0.00000 +DEAL:0:2d/3d::Particle index 13 is in cell 2.2 +DEAL:0:2d/3d::Particle location: 0.142387 0.439867 0.00000 +DEAL:0:2d/3d::Particle index 14 is in cell 2.2 +DEAL:0:2d/3d::Particle location: 0.0667096 0.407735 0.00000 +DEAL:0:2d/3d::Particle index 15 is in cell 2.2 +DEAL:0:2d/3d::Particle location: 0.196835 0.456176 0.00000 +DEAL:0:2d/3d::OK +DEAL:0:3d/3d::Rank: 0. Locally owned active cells: 32 +DEAL:0:3d/3d::Global particles: 16 +DEAL:0:3d/3d::Cell 2.0 has 8 particles. +DEAL:0:3d/3d::Cell 2.1 has 3 particles. +DEAL:0:3d/3d::Cell 2.2 has 2 particles. +DEAL:0:3d/3d::Cell 2.3 has 0 particles. +DEAL:0:3d/3d::Cell 2.4 has 3 particles. +DEAL:0:3d/3d::Cell 2.5 has 0 particles. +DEAL:0:3d/3d::Cell 2.6 has 0 particles. +DEAL:0:3d/3d::Cell 2.7 has 0 particles. +DEAL:0:3d/3d::Cell 2.8 has 0 particles. +DEAL:0:3d/3d::Cell 2.9 has 0 particles. +DEAL:0:3d/3d::Cell 2.10 has 0 particles. +DEAL:0:3d/3d::Cell 2.11 has 0 particles. +DEAL:0:3d/3d::Cell 2.12 has 0 particles. +DEAL:0:3d/3d::Cell 2.13 has 0 particles. +DEAL:0:3d/3d::Cell 2.14 has 0 particles. +DEAL:0:3d/3d::Cell 2.15 has 0 particles. +DEAL:0:3d/3d::Cell 2.16 has 0 particles. +DEAL:0:3d/3d::Cell 2.17 has 0 particles. +DEAL:0:3d/3d::Cell 2.18 has 0 particles. +DEAL:0:3d/3d::Cell 2.19 has 0 particles. +DEAL:0:3d/3d::Cell 2.20 has 0 particles. +DEAL:0:3d/3d::Cell 2.21 has 0 particles. +DEAL:0:3d/3d::Cell 2.22 has 0 particles. +DEAL:0:3d/3d::Cell 2.23 has 0 particles. +DEAL:0:3d/3d::Cell 2.24 has 0 particles. +DEAL:0:3d/3d::Cell 2.25 has 0 particles. +DEAL:0:3d/3d::Cell 2.26 has 0 particles. +DEAL:0:3d/3d::Cell 2.27 has 0 particles. +DEAL:0:3d/3d::Cell 2.28 has 0 particles. +DEAL:0:3d/3d::Cell 2.29 has 0 particles. +DEAL:0:3d/3d::Cell 2.30 has 0 particles. +DEAL:0:3d/3d::Cell 2.31 has 0 particles. +DEAL:0:3d/3d::Particle index 0 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.200844 0.125494 0.185603 +DEAL:0:3d/3d::Particle index 1 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.197844 0.128636 0.203122 +DEAL:0:3d/3d::Particle index 2 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.0413383 0.0199530 0.0497937 +DEAL:0:3d/3d::Particle index 3 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.217265 0.182363 0.0966280 +DEAL:0:3d/3d::Particle index 4 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.219574 0.00346588 0.0326305 +DEAL:0:3d/3d::Particle index 5 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.00256411 0.206502 0.174084 +DEAL:0:3d/3d::Particle index 6 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.0756403 0.0224136 0.220285 +DEAL:0:3d/3d::Particle index 7 is in cell 2.0 +DEAL:0:3d/3d::Particle location: 0.0513904 0.210320 0.104390 +DEAL:0:3d/3d::Particle index 8 is in cell 2.1 +DEAL:0:3d/3d::Particle location: 0.400281 0.221697 0.133304 +DEAL:0:3d/3d::Particle index 9 is in cell 2.1 +DEAL:0:3d/3d::Particle location: 0.312895 0.105785 0.168393 +DEAL:0:3d/3d::Particle index 10 is in cell 2.1 +DEAL:0:3d/3d::Particle location: 0.255528 0.216790 0.0928732 +DEAL:0:3d/3d::Particle index 11 is in cell 2.2 +DEAL:0:3d/3d::Particle location: 0.0659640 0.488469 0.0121682 +DEAL:0:3d/3d::Particle index 12 is in cell 2.2 +DEAL:0:3d/3d::Particle location: 0.0208788 0.454701 0.0928019 +DEAL:0:3d/3d::Particle index 13 is in cell 2.4 +DEAL:0:3d/3d::Particle location: 0.142387 0.189867 0.283471 +DEAL:0:3d/3d::Particle index 14 is in cell 2.4 +DEAL:0:3d/3d::Particle location: 0.0667096 0.157735 0.433275 +DEAL:0:3d/3d::Particle index 15 is in cell 2.4 +DEAL:0:3d/3d::Particle location: 0.196835 0.206176 0.434975 +DEAL:0:3d/3d::OK + +DEAL:1:2d/2d::Rank: 1. Locally owned active cells: 8 +DEAL:1:2d/2d::Global particles: 16 +DEAL:1:2d/2d::Cell 2.8 has 0 particles. +DEAL:1:2d/2d::Cell 2.9 has 0 particles. +DEAL:1:2d/2d::Cell 2.10 has 0 particles. +DEAL:1:2d/2d::Cell 2.11 has 0 particles. +DEAL:1:2d/2d::Cell 2.12 has 0 particles. +DEAL:1:2d/2d::Cell 2.13 has 0 particles. +DEAL:1:2d/2d::Cell 2.14 has 0 particles. +DEAL:1:2d/2d::Cell 2.15 has 0 particles. +DEAL:1:2d/2d::OK +DEAL:1:2d/3d::Rank: 1. Locally owned active cells: 8 +DEAL:1:2d/3d::Global particles: 16 +DEAL:1:2d/3d::Cell 2.8 has 0 particles. +DEAL:1:2d/3d::Cell 2.9 has 0 particles. +DEAL:1:2d/3d::Cell 2.10 has 0 particles. +DEAL:1:2d/3d::Cell 2.11 has 0 particles. +DEAL:1:2d/3d::Cell 2.12 has 0 particles. +DEAL:1:2d/3d::Cell 2.13 has 0 particles. +DEAL:1:2d/3d::Cell 2.14 has 0 particles. +DEAL:1:2d/3d::Cell 2.15 has 0 particles. +DEAL:1:2d/3d::OK +DEAL:1:3d/3d::Rank: 1. Locally owned active cells: 32 +DEAL:1:3d/3d::Global particles: 16 +DEAL:1:3d/3d::Cell 2.32 has 0 particles. +DEAL:1:3d/3d::Cell 2.33 has 0 particles. +DEAL:1:3d/3d::Cell 2.34 has 0 particles. +DEAL:1:3d/3d::Cell 2.35 has 0 particles. +DEAL:1:3d/3d::Cell 2.36 has 0 particles. +DEAL:1:3d/3d::Cell 2.37 has 0 particles. +DEAL:1:3d/3d::Cell 2.38 has 0 particles. +DEAL:1:3d/3d::Cell 2.39 has 0 particles. +DEAL:1:3d/3d::Cell 2.40 has 0 particles. +DEAL:1:3d/3d::Cell 2.41 has 0 particles. +DEAL:1:3d/3d::Cell 2.42 has 0 particles. +DEAL:1:3d/3d::Cell 2.43 has 0 particles. +DEAL:1:3d/3d::Cell 2.44 has 0 particles. +DEAL:1:3d/3d::Cell 2.45 has 0 particles. +DEAL:1:3d/3d::Cell 2.46 has 0 particles. +DEAL:1:3d/3d::Cell 2.47 has 0 particles. +DEAL:1:3d/3d::Cell 2.48 has 0 particles. +DEAL:1:3d/3d::Cell 2.49 has 0 particles. +DEAL:1:3d/3d::Cell 2.50 has 0 particles. +DEAL:1:3d/3d::Cell 2.51 has 0 particles. +DEAL:1:3d/3d::Cell 2.52 has 0 particles. +DEAL:1:3d/3d::Cell 2.53 has 0 particles. +DEAL:1:3d/3d::Cell 2.54 has 0 particles. +DEAL:1:3d/3d::Cell 2.55 has 0 particles. +DEAL:1:3d/3d::Cell 2.56 has 0 particles. +DEAL:1:3d/3d::Cell 2.57 has 0 particles. +DEAL:1:3d/3d::Cell 2.58 has 0 particles. +DEAL:1:3d/3d::Cell 2.59 has 0 particles. +DEAL:1:3d/3d::Cell 2.60 has 0 particles. +DEAL:1:3d/3d::Cell 2.61 has 0 particles. +DEAL:1:3d/3d::Cell 2.62 has 0 particles. +DEAL:1:3d/3d::Cell 2.63 has 0 particles. +DEAL:1:3d/3d::OK +