Skip to content

Make Particle distribution score/statistics postprocessors work with multiple particle managers #6908

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
JarettBakerDunn wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Make Particle distribution score/statistics postprocessors work with multiple particle managers #6908

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b93b94e
Make particle_distribution_score work correctly with multiple particl…
JarettBakerDunn Feb 27, 2026
fbcffc4
make particle_distribution_statistics work with multiple particle man…
JarettBakerDunn Mar 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions include/aspect/postprocess/particle_distribution_score.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,7 @@ namespace aspect
*/
Table<dim,unsigned int>
sort_particles_into_buckets(const typename Triangulation<dim>::active_cell_iterator &cell,
const unsigned int particle_manager_index,
const double bucket_width) const;
};
}
Expand Down
292 changes: 162 additions & 130 deletions source/postprocess/particle_distribution_score.cc
View file
Open in desktop

Large diffs are not rendered by default.

173 changes: 102 additions & 71 deletions source/postprocess/particle_distribution_statistics.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,38 +37,31 @@ namespace aspect
std::pair<std::string,std::string>
ParticleDistributionStatistics<dim>::execute (TableHandler &statistics)
{
unsigned int cells_with_particles = 0;
double standard_deviation_sum = 0;
double standard_deviation_min = std::numeric_limits<double>::max();
double standard_deviation_max = std::numeric_limits<double>::min();
double function_min_sum = 0;
double function_max_sum = 0;
double function_min_min = std::numeric_limits<double>::max();
double function_max_max = std::numeric_limits<double>::min();
// These need to be vectors to account for multiple particle managers.
std::vector<unsigned int> cells_with_particles(this->n_particle_managers(),0);
std::vector<double> standard_deviation_sums(this->n_particle_managers(),0);
std::vector<double> standard_deviation_mins(this->n_particle_managers(),std::numeric_limits<double>::max());
std::vector<double> standard_deviation_maxs(this->n_particle_managers(),0);
std::vector<double> function_min_sums(this->n_particle_managers(),0);
std::vector<double> function_max_sums(this->n_particle_managers(),0);
std::vector<double> function_min_mins(this->n_particle_managers(),std::numeric_limits<double>::max());
std::vector<double> function_max_maxs(this->n_particle_managers(),std::numeric_limits<double>::min());


for (const auto &cell : this->get_dof_handler().active_cell_iterators())
{
if (cell->is_locally_owned())
{
unsigned int particles_in_cell = 0;
for (unsigned int particle_manager_index = 0; particle_manager_index < this->n_particle_managers(); ++particle_manager_index)
particles_in_cell += this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);

if (particles_in_cell > 0)
{
++cells_with_particles;
unsigned int particles_in_cell = this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);

if (KDE_per_particle == false)
if (particles_in_cell > 0)
{
/*
Loop through every particle manager here, since the ParticlePDF class operates
on a single particle_iterator_range. The ParticlePDF class is written to operate
on a single particle_iterator_range at a time so that it can be called directly
by the particle_manager class to assist in particle load balancing.
*/
for (unsigned int particle_manager_index = 0; particle_manager_index < this->n_particle_managers(); ++particle_manager_index)
++cells_with_particles[particle_manager_index];
if (KDE_per_particle == false)
{

const auto &particle_handler = this->get_particle_manager(particle_manager_index).get_particle_handler();
Particle::ParticlePDF<dim> pdf(granularity,bandwidth,kernel_function);
const std::vector<typename Particles::ParticleHandler<dim>::particle_iterator_range> particle_ranges_to_sum_over = get_neighboring_particle_ranges(cell,particle_handler);
Expand All @@ -80,19 +73,16 @@ namespace aspect
cell);
pdf.compute_statistical_values();

standard_deviation_min = std::min(standard_deviation_min, pdf.get_standard_deviation());
standard_deviation_max = std::max(standard_deviation_max, pdf.get_standard_deviation());
standard_deviation_sum += pdf.get_standard_deviation();
standard_deviation_mins[particle_manager_index] = std::min(standard_deviation_mins[particle_manager_index], pdf.get_standard_deviation());
standard_deviation_maxs[particle_manager_index] = std::max(standard_deviation_maxs[particle_manager_index], pdf.get_standard_deviation());
standard_deviation_sums[particle_manager_index] += pdf.get_standard_deviation();

function_min_sum += pdf.get_min();
function_max_sum += pdf.get_max();
function_min_min = std::min(function_min_min, pdf.get_min());
function_max_max = std::max(function_max_max, pdf.get_max());
function_min_sums[particle_manager_index] += pdf.get_min();
function_max_sums[particle_manager_index] += pdf.get_max();
function_min_mins[particle_manager_index] = std::min(function_min_mins[particle_manager_index], pdf.get_min());
function_max_maxs[particle_manager_index] = std::max(function_max_maxs[particle_manager_index], pdf.get_max());
}
}
else
{
for (unsigned int particle_manager_index = 0; particle_manager_index < this->n_particle_managers(); ++particle_manager_index)
else
{
const auto &particle_handler = this->get_particle_manager(particle_manager_index).get_particle_handler();
Particle::ParticlePDF<dim> pdf(bandwidth,kernel_function);
Expand All @@ -105,52 +95,93 @@ namespace aspect
cell);
pdf.compute_statistical_values();

standard_deviation_min = std::min(standard_deviation_min, pdf.get_standard_deviation());
standard_deviation_max = std::max(standard_deviation_max, pdf.get_standard_deviation());
standard_deviation_sum += pdf.get_standard_deviation();
standard_deviation_mins[particle_manager_index] = std::min(standard_deviation_mins[particle_manager_index], pdf.get_standard_deviation());
standard_deviation_maxs[particle_manager_index] = std::max(standard_deviation_maxs[particle_manager_index], pdf.get_standard_deviation());
standard_deviation_sums[particle_manager_index] += pdf.get_standard_deviation();

function_min_sum += pdf.get_min();
function_max_sum += pdf.get_max();
function_min_min = std::min(function_min_min, pdf.get_min());
function_max_max = std::max(function_max_max, pdf.get_max());
function_min_sums[particle_manager_index] += pdf.get_min();
function_max_sums[particle_manager_index] += pdf.get_max();
function_min_mins[particle_manager_index] = std::min(function_min_mins[particle_manager_index], pdf.get_min());
function_max_maxs[particle_manager_index] = std::max(function_max_maxs[particle_manager_index], pdf.get_max());
}
}
}
}
}

// Get final values from all processors
const double global_standard_deviation_max = Utilities::MPI::max (standard_deviation_max, this->get_mpi_communicator());
const double global_standard_deviation_min = Utilities::MPI::min (standard_deviation_min, this->get_mpi_communicator());
const double global_cells_with_particles = Utilities::MPI::sum (cells_with_particles, this->get_mpi_communicator());
const double global_standard_deviation_sum = Utilities::MPI::sum (standard_deviation_sum, this->get_mpi_communicator());
const double global_standard_deviation_mean = global_standard_deviation_sum/global_cells_with_particles;
const double global_function_min_min = Utilities::MPI::min(function_min_min,this->get_mpi_communicator());
const double global_function_max_max = Utilities::MPI::min(function_max_max,this->get_mpi_communicator());

// Get the average of the functions max and min values
const double global_function_min_sum = Utilities::MPI::sum (function_min_sum, this->get_mpi_communicator());
const double global_function_max_sum = Utilities::MPI::sum (function_max_sum, this->get_mpi_communicator());
const double global_function_min_mean = global_function_min_sum/global_cells_with_particles;
const double global_function_max_mean = global_function_max_sum/global_cells_with_particles;

// Write to statistics file
statistics.add_value ("Minimum PDF standard deviation: ", global_standard_deviation_min);
statistics.add_value ("Mean of PDF standard deviation: ", global_standard_deviation_mean);
statistics.add_value ("Maximum PDF standard deviation: ", global_standard_deviation_max);
statistics.add_value ("Mean of PDF minimum values: ", global_function_min_mean);
statistics.add_value ("Mean PDF maximum values: ", global_function_max_mean);
statistics.add_value ("Minimum of PDF minimum values: ", global_function_min_min);
statistics.add_value ("Maximum of PDF maximum values: ", global_function_max_max);


std::ostringstream output;
output << global_standard_deviation_min << "/" << global_standard_deviation_mean << "/" << global_standard_deviation_max << ", "
<< global_function_min_mean << "/" << global_function_max_mean << ", " << global_function_min_min << "/" << global_function_max_max;


return std::pair<std::string, std::string> ("Particle Distribution Stats (stddev min/mean/max, mean min/max, absolute min/max):",
output.str());
if (this->n_particle_managers()==1)
{
// Get final values from all processors
const double global_standard_deviation_max = Utilities::MPI::max (standard_deviation_maxs[0], this->get_mpi_communicator());
const double global_standard_deviation_min = Utilities::MPI::min (standard_deviation_mins[0], this->get_mpi_communicator());
const double global_cells_with_particles = Utilities::MPI::sum (cells_with_particles[0], this->get_mpi_communicator());
const double global_standard_deviation_sum = Utilities::MPI::sum (standard_deviation_sums[0], this->get_mpi_communicator());
const double global_standard_deviation_mean = global_standard_deviation_sum/global_cells_with_particles;
const double global_function_min_min = Utilities::MPI::min(function_min_mins[0],this->get_mpi_communicator());
const double global_function_max_max = Utilities::MPI::min(function_max_maxs[0],this->get_mpi_communicator());

// Get the average of the functions max and min values
const double global_function_min_sum = Utilities::MPI::sum (function_min_sums[0], this->get_mpi_communicator());
const double global_function_max_sum = Utilities::MPI::sum (function_max_sums[0], this->get_mpi_communicator());
const double global_function_min_mean = global_function_min_sum/global_cells_with_particles;
const double global_function_max_mean = global_function_max_sum/global_cells_with_particles;

// Write to statistics file
statistics.add_value ("Minimum PDF standard deviation: ", global_standard_deviation_min);
statistics.add_value ("Mean of PDF standard deviation: ", global_standard_deviation_mean);
statistics.add_value ("Maximum PDF standard deviation: ", global_standard_deviation_max);
statistics.add_value ("Mean of PDF minimum values: ", global_function_min_mean);
statistics.add_value ("Mean PDF maximum values: ", global_function_max_mean);
statistics.add_value ("Minimum of PDF minimum values: ", global_function_min_min);
statistics.add_value ("Maximum of PDF maximum values: ", global_function_max_max);


std::ostringstream output;
output << global_standard_deviation_min << "/" << global_standard_deviation_mean << "/" << global_standard_deviation_max << ", "
<< global_function_min_mean << "/" << global_function_max_mean << ", " << global_function_min_min << "/" << global_function_max_max;


return std::pair<std::string, std::string> ("Particle Distribution Stats (stddev min/mean/max, mean min/max, absolute min/max):",
output.str());
}
else
{
std::ostringstream output;
for (unsigned int particle_manager_index = 0; particle_manager_index < this->n_particle_managers(); ++particle_manager_index)
{
// Get final values from all processors
const double global_standard_deviation_max = Utilities::MPI::max (standard_deviation_maxs[particle_manager_index], this->get_mpi_communicator());
const double global_standard_deviation_min = Utilities::MPI::min (standard_deviation_mins[particle_manager_index], this->get_mpi_communicator());
const double global_cells_with_particles = Utilities::MPI::sum (cells_with_particles[particle_manager_index], this->get_mpi_communicator());
const double global_standard_deviation_sum = Utilities::MPI::sum (standard_deviation_sums[particle_manager_index], this->get_mpi_communicator());
const double global_standard_deviation_mean = global_standard_deviation_sum/global_cells_with_particles;
const double global_function_min_min = Utilities::MPI::min(function_min_mins[particle_manager_index],this->get_mpi_communicator());
const double global_function_max_max = Utilities::MPI::min(function_max_maxs[particle_manager_index],this->get_mpi_communicator());

// Get the average of the functions max and min values
const double global_function_min_sum = Utilities::MPI::sum (function_min_sums[particle_manager_index], this->get_mpi_communicator());
const double global_function_max_sum = Utilities::MPI::sum (function_max_sums[particle_manager_index], this->get_mpi_communicator());
const double global_function_min_mean = global_function_min_sum/global_cells_with_particles;
const double global_function_max_mean = global_function_max_sum/global_cells_with_particles;

// Write to statistics file
std::string particle_manager_index_prefix = "Particle Manager " + std::to_string(particle_manager_index) + " ";

statistics.add_value (particle_manager_index_prefix+"Minimum PDF standard deviation: ", global_standard_deviation_min);
statistics.add_value (particle_manager_index_prefix+"Mean of PDF standard deviation: ", global_standard_deviation_mean);
statistics.add_value (particle_manager_index_prefix+"Maximum PDF standard deviation: ", global_standard_deviation_max);
statistics.add_value (particle_manager_index_prefix+"Mean of PDF minimum values: ", global_function_min_mean);
statistics.add_value (particle_manager_index_prefix+"Mean PDF maximum values: ", global_function_max_mean);
statistics.add_value (particle_manager_index_prefix+"Minimum of PDF minimum values: ", global_function_min_min);
statistics.add_value (particle_manager_index_prefix+"Maximum of PDF maximum values: ", global_function_max_max);

output << "Particle Manager index " << particle_manager_index << ": " << global_standard_deviation_min << "/" << global_standard_deviation_mean << "/" << global_standard_deviation_max << ", "
<< global_function_min_mean << "/" << global_function_max_mean << ", " << global_function_min_min << "/" << global_function_max_max << ", ";

}
return std::pair<std::string, std::string> ("Particle Distribution Stats (stddev min/mean/max, mean min/max, absolute min/max): ",
output.str());
}
}


Expand Down
131 changes: 131 additions & 0 deletions tests/particle_score_statistics_multiple_particle_managers.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
# This parameter file is used to test running the particle distribution score and statistics
# postprocessors when they are used on models with multiple particle managers.

set Dimension = 2
set Use years instead of seconds = false
set End time = 1
set Nonlinear solver scheme = single Advection, no Stokes
set Resume computation = false
set Maximum time step = 0.03

subsection Geometry model
set Model name = box

subsection Box
set Y periodic = false
set X extent = 0.5
set Y extent = 5
set Y repetitions = 10
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 Particles
set Number of particle systems = 2
set Particle generator name = reference cell
set Maximum particles per cell = 20
set Minimum particles per cell = 16
set Load balancing strategy = remove and add particles
set Particle removal algorithm = random
set Particle addition algorithm = random
subsection Generator
subsection Reference cell
set Number of particles per cell per direction = 4
end
end
end

subsection Particles 2
set Particle generator name = reference cell
set Maximum particles per cell = 12
set Minimum particles per cell = 9
set Load balancing strategy = remove and add particles
set Particle removal algorithm = random
set Particle addition algorithm = random
subsection Generator
subsection Reference cell
set Number of particles per cell per direction = 3
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


subsection Prescribed Stokes solution
set Model name = function
subsection Velocity function
set Variable names = x,y,t
set Function constants = velConstant=-0.1
set Function expression = 0; velConstant
end
end
Loading
Loading