Make Particle distribution score/statistics postprocessors work with multiple particle managers by JarettBakerDunn · Pull Request #6908 · geodynamics/aspect

JarettBakerDunn · 2026-03-19T16:17:49Z

Before this pull request the Particle distribution score and Particle distribution statistics postprocessors would calculate values across multiple particle managers (if there were multiple particle managers used in the model). This would make the results of the postprocessors inaccurate since the particle managers take care of particles separately from one another.

Now each postprocessor behaves as before if there is only one particle manager, but if there are more than one the postprocessors will calculate and report their statistics separately for each particle manager.

…e systems

…agers

gassmoeller

Yes, the ideas are all correct, but I would suggest modifications to the implementations, see below. Let me know when I should take another look.

gassmoeller · 2026-03-30T13:19:44Z

-                      for (unsigned int y=0; y<granularity; ++y)
+                      const double particles_in_cell_double = static_cast<double>(particles_in_cell);
+                      const double granularity_double = static_cast<double>(granularity);
+                      const double ideal_n_particles_per_bucket = particles_in_cell_double/(Utilities::fixed_power<dim>(granularity_double));


you can compute granularity_double and ideal_n_particles_per_bucket outside of the loop over all cells, that will save compute time.

I agree for granularity_double, but I think that ideal_n_particles_per_bucket depends on the number of particles in each cell which could vary depending on the cell. For example in 2D with a granularity of 2 we would have 4 buckets per cell. Then a cell with 4 particle would want 1 particle per bucket, while a cell with 8 particles would want 2 particles per bucket? Or were you thinking of calculating the ideal particles per bucket from the minimum or maximum # of allowed particles per cell?

Correct, I wasnt paying attention.

gassmoeller · 2026-03-30T13:22:20Z

+                      The "buckets_ideal" table contains doubles because the ideal
+                      number of particles per bucket will not always be an integer.
+                      */
+                      Table<dim,double> buckets_ideal;


I didnt catch this before, but it seems you dont need buckets_ideal in this postprocessor at all. You fill all entries of buckets_ideal with ideal_n_particles_per_bucket and then later use the entries. That means you can just use ideal_n_particles_per_bucket in line 104 and 111 and remove buckets_ideal altogether.

gassmoeller · 2026-03-30T13:23:42Z

-      // get final values for min and max score from all processors
-      const double global_max_score = Utilities::MPI::max (local_max_score, this->get_mpi_communicator());
-      const double global_min_score = Utilities::MPI::min (local_min_score, this->get_mpi_communicator());
+      if (this->n_particle_managers()==1)


I dont think you want the if condition here. See my suggestion below.

gassmoeller · 2026-03-30T13:31:54Z

-      std::ostringstream output;
-      output << global_min_score << '/' << mean_and_standard_deviation.first << '/' << global_max_score << '/' << mean_and_standard_deviation.second;
+          return std::pair<std::string, std::string> ("Particle Distribution Score (min/avg/max/stdev): ",
+                                                      output.str());


I would generalize this to the case of 1 particle handler as well, that will simplify your code. I.e. remove the whole if branch above, and change this to:

Suggested change

output.str());

std::ostringstream output;

for (unsigned int particle_manager_index = 0; particle_manager_index < this->n_particle_managers(); ++particle_manager_index)

{

// Calculate the mean and standard deviation of cell scores across all processors

const std::pair<double, double> mean_and_standard_deviation =

Utilities::MPI::mean_and_standard_deviation(cell_scores[particle_manager_index].begin(),

cell_scores[particle_manager_index].end(),

this->get_mpi_communicator());

// get final values for min and max score from all processors

const double global_max_score = Utilities::MPI::max (local_max_scores[particle_manager_index], this->get_mpi_communicator());

const double global_min_score = Utilities::MPI::min (local_min_scores[particle_manager_index], this->get_mpi_communicator());

// write to statistics file for all particle managers

const std::string particle_manager_index_prefix = (particle_manager_index == 0) ? "" : "Particles " + std::to_string(particle_manager_index+1) + ": ";

statistics.add_value (particle_manager_index_prefix+"Minimal particle distribution score: ", global_min_score);

statistics.add_value (particle_manager_index_prefix+"Average particle distribution score: ", mean_and_standard_deviation.first);

statistics.add_value (particle_manager_index_prefix+"Maximal particle distribution score: ", global_max_score);

statistics.add_value (particle_manager_index_prefix+"Cell Score Standard Deviation: ", mean_and_standard_deviation.second);

// write screen output for the first particle manager

if (particle_manager_index == 0)

output << global_min_score << '/' << mean_and_standard_deviation.first << '/' << global_max_score << '/' << mean_and_standard_deviation.second;

}

return std::pair<std::string, std::string> ("Particle distribution score min/avg/max/stdev: ",

output.str());

Also leave the old output capitalization and wording of the return string, it is more consistent with our other postprocessors. (It is the particle distribution statistics postprocessor that is not following conventions and should be changed, but not in this PR).

Thanks for letting me know, I will open another PR correcting the wording of particle distribution statistics after this.

gassmoeller · 2026-03-30T13:35:50Z

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


Same comment as above. Remove the if-branch. Modify the else-branch to output the statistics file and the screen output as before for particle manager 0. Do not output screen output for particle manager >0, and prefix the statistics columns for particle manager >0 with "Particles x:" where x is particle manager index +1 (this is so that it is consistent with the numbering we use in the input file).

gassmoeller · 2026-03-30T13:36:11Z

     Writing particle output:                                                           output-benchmark_particle_pdf_removal/particles/particles-00000
-     Particle distribution score min/avg/max/stdev:                                     0/0/0/0
-     Particle Distribution Stats (stddev min/mean/max, mean min/max, absolute min/max): 0/0/2.22507e-308, 0.0592697/0.0592697, 0.0454545/0.1
+     Particle Distribution Score (min/avg/max/stdev):                                   0/0/0/0


with my suggested changed above, you shouldnt see any test changes.

gassmoeller · 2026-03-30T13:37:04Z

This file needs to be named screen-output, not screen-output.txt. It is a leftover from very early in ASPECT development.

gassmoeller

Close, but the test results are still changing. Check what happens if you include my suggestions.

gassmoeller · 2026-03-31T07:29:49Z

-
-                  double actual_error_squared = 0;
-                  for (unsigned int x=0; x<granularity; ++x)
+                  unsigned int particles_in_cell = this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);


Suggested change

unsigned int particles_in_cell = this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);

const unsigned int particles_in_cell = this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);

gassmoeller · 2026-03-31T07:31:52Z

+                      TableIndices<dim> bucket_sizes;
+                      for (unsigned int i=0; i<dim; ++i)
+                        bucket_sizes[i] = granularity;


Are bucket_sizes still used? If not, they can be removed.

gassmoeller · 2026-03-31T07:39:34Z

+      // These need to be vectors to account for multiple particle managers.
+      std::vector<double> local_min_scores(this->n_particle_managers(),std::numeric_limits<double>::max());
+      std::vector<double> local_max_scores(this->n_particle_managers(),0);
+      std::vector<std::vector<double>> cell_scores(this->n_particle_managers());


While we are touching this line anyway: Can you add the following:

Suggested change

std::vector<std::vector<double>> cell_scores(this->n_particle_managers());

std::vector<std::vector<double>> cell_scores(this->n_particle_managers());

for (auto &scores: cell_scores)

scores.reserve(this->get_triangulation().n_active_cells());

This will avoid unnecessary memory allocation and copy operations.

gassmoeller · 2026-03-31T07:40:20Z

-                      for (unsigned int y=0; y<granularity; ++y)
+                      const double particles_in_cell_double = static_cast<double>(particles_in_cell);
+                      const double granularity_double = static_cast<double>(granularity);
+                      const double ideal_n_particles_per_bucket = particles_in_cell_double/(Utilities::fixed_power<dim>(granularity_double));


Correct, I wasnt paying attention.

gassmoeller · 2026-03-31T07:48:15Z

+          if (particle_manager_index == 0)
+            {
+              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 << ", ";


Suggested change

<< global_function_min_mean << "/" << global_function_max_mean << ", " << global_function_min_min << "/" << global_function_max_max << ", ";

<< global_function_min_mean << "/" << global_function_max_mean << ", " << global_function_min_min << "/" << global_function_max_max;

gassmoeller · 2026-03-31T07:48:43Z

+                     << 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): ",


Suggested change

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

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

JarettBakerDunn · 2026-03-31T21:27:19Z

I think this is ready for another review.

gassmoeller · 2026-04-03T07:14:51Z

/rebuild

gassmoeller

Thanks! Except for my one question about the test output this is ready. Please have a look and let me know what you find.

gassmoeller · 2026-04-03T07:30:10Z

+     Writing graphical output:                                                           output-benchmark_particle_pdf_removal/solution/solution-00000
+     Writing particle output:                                                            output-benchmark_particle_pdf_removal/particles/particles-00000
+     Particle distribution score min/avg/max/stdev:                                      0/0/0/0
+     Particle Distribution Stats (stddev min/mean/max, mean min/max, absolute min/max):  0/0/0, 0.0592697/0.0592697, 0.0454545/0.1, 


I dont understand why there is a trailing comma in this line. I dont see it written in the code, but the tester seems to accept it. Is it a leftover from earlier versions of this PR? Do you see where it is coming from?

JarettBakerDunn added 3 commits February 27, 2026 15:16

Make particle_distribution_score work correctly with multiple particl…

b93b94e

…e systems

make particle_distribution_statistics work with multiple particle man…

fbcffc4

…agers

apply diff

0612739

gassmoeller reviewed Mar 30, 2026

View reviewed changes

JarettBakerDunn added 2 commits March 30, 2026 12:55

address feedback

171f219

run make indent and apply diff (back to correct capitalization/spelling)

01e0dde

gassmoeller reviewed Mar 31, 2026

View reviewed changes

JarettBakerDunn added 2 commits March 31, 2026 12:34

apply feedback

a4e311e

make_indent

8061a38

gassmoeller added the ready to test label Apr 3, 2026

gassmoeller reviewed Apr 3, 2026

View reviewed changes

-                                                      output.str());
+          std::ostringstream output;
+          for (unsigned int particle_manager_index = 0; particle_manager_index < this->n_particle_managers(); ++particle_manager_index)
+            {
+              // Calculate the mean and standard deviation of cell scores across all processors
+              const std::pair<double, double> mean_and_standard_deviation =
+                Utilities::MPI::mean_and_standard_deviation(cell_scores[particle_manager_index].begin(),
+                                                            cell_scores[particle_manager_index].end(),
+                                                            this->get_mpi_communicator());
+              // get final values for min and max score from all processors
+              const double global_max_score = Utilities::MPI::max (local_max_scores[particle_manager_index], this->get_mpi_communicator());
+              const double global_min_score = Utilities::MPI::min (local_min_scores[particle_manager_index], this->get_mpi_communicator());
+              // write to statistics file for all particle managers
+              const std::string particle_manager_index_prefix = (particle_manager_index == 0) ? "" : "Particles " + std::to_string(particle_manager_index+1) + ": ";
+              statistics.add_value (particle_manager_index_prefix+"Minimal particle distribution score: ", global_min_score);
+              statistics.add_value (particle_manager_index_prefix+"Average particle distribution score: ", mean_and_standard_deviation.first);
+              statistics.add_value (particle_manager_index_prefix+"Maximal particle distribution score: ", global_max_score);
+              statistics.add_value (particle_manager_index_prefix+"Cell Score Standard Deviation: ", mean_and_standard_deviation.second);
+            // write screen output for the first particle manager
+            if (particle_manager_index == 0)
+              output << global_min_score << '/' << mean_and_standard_deviation.first << '/' << global_max_score << '/' << mean_and_standard_deviation.second;
+            }
+          return std::pair<std::string, std::string> ("Particle distribution score min/avg/max/stdev: ",
+                                                      output.str());

	unsigned int particles_in_cell = this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);
	const unsigned int particles_in_cell = this->get_particle_manager(particle_manager_index).get_particle_handler().n_particles_in_cell(cell);

	<< global_function_min_mean << "/" << global_function_max_mean << ", " << global_function_min_min << "/" << global_function_max_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): ",
	return std::pair<std::string, std::string> ("Particle Distribution Stats (stddev min/mean/max, mean min/max, absolute min/max):",

Conversation

JarettBakerDunn commented Mar 19, 2026

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JarettBakerDunn Mar 30, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

JarettBakerDunn commented Mar 31, 2026

Uh oh!

gassmoeller commented Apr 3, 2026

Uh oh!

gassmoeller left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

JarettBakerDunn Mar 30, 2026 •

edited

Loading