Let particle managers move particles around within cells to reduce particle clustering

benchmarks/particle_distribution/active_redistribution_benchmark.prm

@@ -0,0 +1,118 @@
+# This parameter file demonstrates the "Active redistribution within cells"
+# parameter of ASPECT's particle manager. The parameter enables a particle manager
+# to move particles around within a cell so that the particles are distributed more
+# evenly throughout the cell.
+
+set Dimension                              = 2
+set Use years instead of seconds           = false
+set End time                               = 6
+set Nonlinear solver scheme                = single Advection, no Stokes
+set Output directory                       = output-active-redistribution
+set Resume computation                     = false
+set Maximum time step                      = 0.03
+
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set Y periodic = true
+    set X extent = 0.5
+    set Y extent = 1
+    set Y repetitions = 2
+  end
+end
+
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  =
+    set Function expression = 1
+  end
+end
+
+subsection Boundary temperature model
+  set List of model names = box
+  set Fixed temperature boundary indicators   = left, right
+
+  subsection Box
+    set Left temperature   = 0
+    set Right temperature  = 0
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0.0
+  end
+end
+
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1e-5
+    set Thermal expansion coefficient = 0
+    set Viscosity                     = 0
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement              = 2
+  set Initial global refinement                = 1
+  set Time steps between mesh refinement       = 0
+
+  set Strategy = minimum refinement function
+
+  subsection Minimum refinement function
+    set Coordinate system   = cartesian
+    set Variable names      = x,y
+    set Function expression = if ( y<=4, 2, 4)
+  end
+
+end
+
+
+subsection Prescribed Stokes solution
+  set Model name = function
+  subsection Velocity function
+    set Variable names      = x,y,t
+    set Function constants  = velConstant=1
+    set Function expression = 0; if (y<=0.5, velConstant, velConstant/5)
+  end
+end
+
+
+subsection Particles
+  set Particle generator name =  uniform box
+  set Active redistribution within cells = true
+  set Maximum distribution score per cell = 0.01
+  subsection Generator
+    subsection Uniform box
+      set Number of particles = 500
+      set Maximum y = 0.5
+    end
+  end
+
+end
+
+subsection Postprocess
+  set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, visualization, particles, particle distribution score,particle distribution statistics
+
+  subsection Visualization
+    set Time between graphical output = .001
+    set Interpolate output = false
+  end
+
+  subsection Particles
+    set Time between data output = 0.001
+    set Data output format       = vtu
+  end
+
+end

doc/modules/changes/20260608_bakerdunn

@@ -0,0 +1,5 @@
+New: The particle manager class can now redistribute particles
+within cells if it calculates that the particles within that
+cell are too clustered.
+<br>
+(Jarett Baker-Dunn, 2026/06/08)

doc/sphinx/parameters/Particles.md

@@ -5,6 +5,24 @@
 ## **Subsection:** Particles
 
 
+::::{dropdown} __Parameter:__ {ref}`Active redistribution histogram granularity<parameters:Particles/Active_20redistribution_20histogram_20granularity>`
+:name: parameters:Particles/Active_20redistribution_20histogram_20granularity
+**Default value:** 2
+
+**Pattern:** [Integer range 2...2147483647 (inclusive)]
+
+**Documentation:** The number of subdivisions of each cell in each spatial dimension when evaluating cells for particle clustering when actively redistributing particles.
+::::
+
+::::{dropdown} __Parameter:__ {ref}`Active redistribution within cells<parameters:Particles/Active_20redistribution_20within_20cells>`
+:name: parameters:Particles/Active_20redistribution_20within_20cells
+**Default value:** false
+
+**Pattern:** [Bool]
+
+**Documentation:** Sometimes particles will cluster in one part of the cell, which can make the particles less useful for visualization. When this parameter is set to true, the particle manager will use a histogram to quantify the level of particle clustering in each cell. If the particles are sufficiently clustered, the particle manager will redistribute the particles in the cell so that they are more evenly spread throughout the cell.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Addition histogram granularity<parameters:Particles/Addition_20histogram_20granularity>`
 :name: parameters:Particles/Addition_20histogram_20granularity
 **Default value:** 3
@@ -36,11 +54,20 @@
 :name: parameters:Particles/Bandwidth
 **Default value:** 0.3
 
-**Pattern:** [Double 0.3...MAX_DOUBLE (inclusive)]
+**Pattern:** [Double 0.1...0.9 (inclusive)]
 
 **Documentation:** The bandwidth value is used to scale the kernel function when generating the point density function of particles. The bandwidth is measured as a fraction of the cells extent in one spatial dimension. For example, the default bandwidth of 0.3 represents a size equal to 30 percent of the cells size in one spatial dimension.
 ::::
 
+::::{dropdown} __Parameter:__ {ref}`Distribution attempts max<parameters:Particles/Distribution_20attempts_20max>`
+:name: parameters:Particles/Distribution_20attempts_20max
+**Default value:** 10
+
+**Pattern:** [Integer range 1...100 (inclusive)]
+
+**Documentation:** The maximum amount of times to move a particle in order to reduce a cell&rsquo;s distribution score below the Maximum distribution score per cell before giving up. The particle manager will stop moving particles if the distribution score is brought below the Maximum distribution score per cell even if the number of particles moved is below this number.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Integration scheme<parameters:Particles/Integration_20scheme>`
 :name: parameters:Particles/Integration_20scheme
 **Default value:** rk2
@@ -149,6 +176,15 @@ The following properties are available:
 **Documentation:** Strategy that is used to balance the computational load across processors for adaptive meshes.
 ::::
 
+::::{dropdown} __Parameter:__ {ref}`Maximum distribution score per cell<parameters:Particles/Maximum_20distribution_20score_20per_20cell>`
+:name: parameters:Particles/Maximum_20distribution_20score_20per_20cell
+**Default value:** 0.10
+
+**Pattern:** [Double 0.001...1 (inclusive)]
+
+**Documentation:** The maximum allowable distribution score a cell can have before its particles are redistributed inside the cell. The distribution score is measured using a histogram based method which is described in more detail in the documentation for the Particle Distribution Score plugin. Paticles are only redistributed within cells if the &lsquo;Active redistribution within cells&lsquo; parameter is set to true.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Maximum particles per cell<parameters:Particles/Maximum_20particles_20per_20cell>`
 :name: parameters:Particles/Maximum_20particles_20per_20cell
 **Default value:** 100

doc/sphinx/parameters/Particles_202.md

@@ -5,6 +5,24 @@
 ## **Subsection:** Particles 2
 
 
+::::{dropdown} __Parameter:__ {ref}`Active redistribution histogram granularity<parameters:Particles_202/Active_20redistribution_20histogram_20granularity>`
+:name: parameters:Particles_202/Active_20redistribution_20histogram_20granularity
+**Default value:** 2
+
+**Pattern:** [Integer range 2...2147483647 (inclusive)]
+
+**Documentation:** The number of subdivisions of each cell in each spatial dimension when evaluating cells for particle clustering when actively redistributing particles.
+::::
+
+::::{dropdown} __Parameter:__ {ref}`Active redistribution within cells<parameters:Particles_202/Active_20redistribution_20within_20cells>`
+:name: parameters:Particles_202/Active_20redistribution_20within_20cells
+**Default value:** false
+
+**Pattern:** [Bool]
+
+**Documentation:** Sometimes particles will cluster in one part of the cell, which can make the particles less useful for visualization. When this parameter is set to true, the particle manager will use a histogram to quantify the level of particle clustering in each cell. If the particles are sufficiently clustered, the particle manager will redistribute the particles in the cell so that they are more evenly spread throughout the cell.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Addition histogram granularity<parameters:Particles_202/Addition_20histogram_20granularity>`
 :name: parameters:Particles_202/Addition_20histogram_20granularity
 **Default value:** 3
@@ -36,11 +54,20 @@
 :name: parameters:Particles_202/Bandwidth
 **Default value:** 0.3
 
-**Pattern:** [Double 0.3...MAX_DOUBLE (inclusive)]
+**Pattern:** [Double 0.1...0.9 (inclusive)]
 
 **Documentation:** The bandwidth value is used to scale the kernel function when generating the point density function of particles. The bandwidth is measured as a fraction of the cells extent in one spatial dimension. For example, the default bandwidth of 0.3 represents a size equal to 30 percent of the cells size in one spatial dimension.
 ::::
 
+::::{dropdown} __Parameter:__ {ref}`Distribution attempts max<parameters:Particles_202/Distribution_20attempts_20max>`
+:name: parameters:Particles_202/Distribution_20attempts_20max
+**Default value:** 10
+
+**Pattern:** [Integer range 1...100 (inclusive)]
+
+**Documentation:** The maximum amount of times to move a particle in order to reduce a cell&rsquo;s distribution score below the Maximum distribution score per cell before giving up. The particle manager will stop moving particles if the distribution score is brought below the Maximum distribution score per cell even if the number of particles moved is below this number.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Integration scheme<parameters:Particles_202/Integration_20scheme>`
 :name: parameters:Particles_202/Integration_20scheme
 **Default value:** rk2
@@ -149,6 +176,15 @@ The following properties are available:
 **Documentation:** Strategy that is used to balance the computational load across processors for adaptive meshes.
 ::::
 
+::::{dropdown} __Parameter:__ {ref}`Maximum distribution score per cell<parameters:Particles_202/Maximum_20distribution_20score_20per_20cell>`
+:name: parameters:Particles_202/Maximum_20distribution_20score_20per_20cell
+**Default value:** 0.10
+
+**Pattern:** [Double 0.001...1 (inclusive)]
+
+**Documentation:** The maximum allowable distribution score a cell can have before its particles are redistributed inside the cell. The distribution score is measured using a histogram based method which is described in more detail in the documentation for the Particle Distribution Score plugin. Paticles are only redistributed within cells if the &lsquo;Active redistribution within cells&lsquo; parameter is set to true.
+::::
+
 ::::{dropdown} __Parameter:__ {ref}`Maximum particles per cell<parameters:Particles_202/Maximum_20particles_20per_20cell>`
 :name: parameters:Particles_202/Maximum_20particles_20per_20cell
 **Default value:** 100

include/aspect/particle/histogram.h

@@ -0,0 +1,85 @@
+/*
+  Copyright (C) 2025 - by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file LICENSE.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _aspect_particle_histogram_h
+#define _aspect_particle_histogram_h
+
+#include <aspect/simulator_access.h>
+#include <deal.II/particles/property_pool.h>
+#include <deal.II/particles/particle_handler.h>
+
+
+namespace aspect
+{
+  namespace Particle
+  {
+    /**
+     * This class computes a spatial histogram describing the distribution of
+     * particles in a cell. The number of spatial bins is equal to granularity^dimensions.
+     */
+    template <int dim>
+    class ParticleHistogram
+    {
+      public:
+        const unsigned int granularity;
+
+        /**
+         * This is the constructor for ParticleHistogram. ParticleHistogram sorts
+         * the particles in a cell into different histogram bins depending on the
+         * position of the particle in the cell. These histogram bins are compared to
+         * a virtual histogram representing a uniform distribution of particles in the
+         * cell and the difference between the real and virtual histograms is reported
+         * as a score representing the degree of particle clustering in the cell.
+         *
+         * @param granularity determines the number spatial bins into which ParticleHistogram
+         * divides each cell. The real number of bins is equal to granularity^dim, where dim
+         * is the number of spatial dimensions, either 2 or 3.
+         */
+        ParticleHistogram(const unsigned int granularity);
+
+        /**
+         * Calculates and returns a score between 0 and 1 representing the level of clustering of the
+         * particle range within the cell, where 0 represents a completely even distribution
+         * and 1 represents the most clustered distribution that the histogram is able to
+         * measure. The precision of the score depends on the ParticleHistogram's granularity
+         * parameter.
+         * @param particle_range is the particle range for which to evaluate the clustering
+         * score
+         * @param n_particles_in_cell is the number of particles in the particle range
+         * @return A double ranging from 0 to 1 representing the degree to which the particles in
+         * the provided particle range are clustered together, with 1 being the most clustered
+         */
+        double evaluate_particle_range(const typename Particles::ParticleHandler<dim>::particle_iterator_range &particle_range,
+                                       const unsigned int n_particles_in_cell);
+
+        /**
+         * Sorts all of the particles within the cell into a deal.II table based on their position.
+         * @param cell The cell for which to compute the particle distribution.
+         * @param bucket_width The size (relative to the size of the cell) of each bucket in the table.
+         * @return The table with the particle information.
+         */
+        Table<dim,unsigned int>
+        sort_particles_into_buckets(const typename Particles::ParticleHandler<dim>::particle_iterator_range &particle_range,
+                                    const double bucket_width) const;
+    };
+  }
+}
+
+#endif