re-generate particles in initial adaptive refinement

[jdannberg]

In large application models with adaptive refinement, one often wants to initialize the simulation with a similar number of particles in each cell (i.e. based on the adaptive refinement). This strategy is useful in models where the number of particles is controlled via a Minimum/Maximum number of particles per cell, because it avoids having to create/delete additional particles later on (which would be added at random locations).

So far, ASPECT only called the particle generator once after the global refinements. This PR changes this to delete and re-generate particles in each initial adaptive refinement step. The result is show in the image below for the reference cell particle generator:

Left: this PR. Right: main.

• I have followed the instructions for indenting my code.
• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[gassmoeller]

I like the approach and I think it is universally more accurate, and the additional cost (generating particles multiple times) is likely small compared to the total model cost, because generation still only happens a few times in timestep 0. Can you:

• List for future reference how much time the generation takes for your test between main and this PR
• Update the test cases that are affected by this change, and give them a look if they still behave as expected?

The code itself looks ready.

[gassmoeller]

It looks like the test particle_interpolator_empty_cells is crashing with this change. Can you take a look? It is likely just an issue that now a cell and all of its neighbors is empty of particles. An adjustment of the test should do the trick (so that individual cells might be empty of particles, but never a cell and all of its neighbors).

[jdannberg]

I started to look at the tests that changed to see if there was any case where one might want to keep the original method. I did not get very far, since already for the first test that failed (grain_size_crossed_transition) I got some confusing results:

Top is new, bottom is old behavior.
In this test, the particle generator is a probability density function using a sine with the maximum in the center. The new version seemed to show some odd behavior at cell boundaries, so I increased the number of particles to see what was happening:

and with the mesh:

Is this expected behavior (i.e., within once cell, there is just one value for the probability density function)? Would one want to keep the old behavior in this case, i.e., should I make it an input parameter which for which behavior to use?

[tjhei]

I got some confusing results

It makes sense to me that the new behavior correctly correlates the density with the final cell size, while the old version does not.

[gassmoeller]

It makes sense to me that the new behavior correctly correlates the density with the final cell size, while the old version does not.

I think the confusing part is more that the new form shows more clearly that the particle density within one cell is always constant. This is because the particle density function is integrated (in some form) over the cell volume and that will be the chance to generate a particle somewhere within that cell. Because the particle positions are now generated on the final adaptive grid cell boundaries now feature more prominently in the particle distribution. I think this is as expected, it just wasnt as visible before.

In a discussion with @jdannberg we couldnt really come up with a situation where one would use the old implementation, but we couldnt completely rule it out either. In order to avoid more confusing input parameters, my suggestion would be to switch the algorithm to the new form as suggested in this PR. If we encounter use cases in the future that would benefit from the old implementation we can introduce an input parameter to switch between the options then.

[jdannberg]

Thanks for looking into this!

I went through all the tests that changed (which all had initial adaptive refinement) to make sure the new test results make sense, and I think they do. Most of the changes are small and because the particles are at different locations now.

There is a large change in the number of particles in the particles_with_AMR_on test, but that is because it uses the reference cell generator and there are way more cells after adaptive refinement, so we create way more particles now. The particle_load_balancing_repartition_multiple_systems now has exactly one particle per cell, but @gassmoeller mentioned that this is not a problem.

I also updated the particle_interpolator_empty_cells test.

Now, for the timing:

For the test case I added (averaged over 4 model runs):
New: Particles: Generate | 2 | 0.0001275s | 0%
Old: Particles: Generate | 1 | 6.785e-05s | 0%

(We now have two calls, so we take twice as long, but it’s still zero percent of the runtime.)

So I thought I’d run another test where we have about as many particles as we would have in an actual model (15 million) and got the following (same test, just more cells, and 100x100 particles in each cell, averaged over 3 runs):

New: | Particles: Generate | 2 | 1.49s | 3% |
Old: | Particles: Generate | 1 | 0.615s | 1.6% |

Generating is still cheap compared to everything else. For example, initialize properties takes about 10% of the time, and that’s something we already do in every initial refinement step.

[gassmoeller]

Can you address my one comment and then merge yourself?

Thanks for the improvement 👍

[gassmoeller]

Can you modify this to make clear this only refers to initial adaptive and global refinement cycles here or in the second sentence? At the moment it sounds like we are regenerating particles during every adaptive refinement.