Randomly disappearing particles in non-trivial geometries (hyper_ball) since the merge of #13202

[blaisb]

Dear all, I hope you are well.
The merge of #13202 has caused some significant issue on our DEM code base. I am having difficulties identifying what is the exact issue, but we have found that particles are starting to disappear randomly when simulations are carried out in an unstructured grid in which the topology is complicated (e.g. an hyperball). Everything seems to go as normal, but at random time-steps, particles will be erased from the simulation. As such, all of our unit tests that take place in such geometries have been "broken" and the simulations don't make much sense anymore. I presume something is going on with the sort_particles_into_subdomain_and_cells. So far I have managed to reproduce the bug on a single core with very very few particles (20ish).

The following two images higlight this:

This is at a time-step, particles are moving down

Next time-step (e.g a few 0.00001m further)

One particle has magically disappeared.

@gassmoeller or @kronbichler maybe one of you would have an idea of what's happening?
Maybe it is because we are not using the mapping correctly at one place and it means we had a lying bug in our code base, but I would just like to figure out how to fix this :).

Take care and be safe all :)
Bruno

[blaisb]

My hypothesis would be this:
https://github.com/dealii/dealii/pull/13202/files#diff-c3dbdb2659a6f4be5037f4b6539d758f47099ec4da8f3838a5fd70a21b0fa742R1404

I am not sure if searching all cell neighbours of the closest vertex is the best idea. I am pretty sure there are some edge cases that are missed by that.
Maybe a better fall back would to search all cell neighbours of closest vertex, then of the other vertices, then the particle would be lost?
The optimisation would still be there, but we would be a lot more robust.

[marcfehling]

@blaisb Can you provide a minimal test that reproduces the issue? Independent of the solution to the problem, it would be good to have it in the deal.II testsuite to prevent it from recurring in the future.

[blaisb]

@marcfehling ,
I will try to cook one up today or tomorrow.
You are right, the deal.II test suite for particles is a bit incomplete. It has happened a few time that some edge case was captured by an outside code (e.g. Lethe) instead of by the test suite. My goal for the Hackathon this year would be to write as many new particle tests as possible to make the library more robust.

[blaisb]

@marcfehling Here is a small application that is a very very very minimal step-68.
It generates a ball of particle in a ball and advects them in the y direction. The time step is relatively small.
Before #13202 you will have 48 particles during the entire simulation
After #13202 at step 159 you will start to loose particles. 4 particles will be lost during the simulation, even if they should not leave the domain.

I also added paraview output as this may help getting a visual on the problem.
Hope this helps. If you want, I will gladly turn this into a unit test once a PR fixing this problem is underway.

particle_disappearance.tar.gz

Example:
Step-158:

Step-159

[gassmoeller]

Hi @blaisb, first of all sorry that my change affected more than intended.

I agree that the behavior could be slightly different after the change, but I am not sure it is for the reason you suggested (but I could be convinced of course). The reason I do not think it is the simpler search option is that this option only triggers if particles move to a cell that is not a neighbor of the old cell. But in your model the particles move only very small amounts of space, so they should always end up in a cell that is a neighbor of the old cell. The only other case if this fallback is triggered is if they move out of the domain (for example if you use a forward euler integration scheme this could easily happen because of integration errors, but also for other integration schemes, just less often). Can you clarify if the disappearing particles are very close to the boundary of the domain? Because in this case what is different would not be the selection of which cells to search, but the actual check if a particle is in a cell. I saw that find_active_cell_around_point has a tolerance of 1e-10 (in reference cell space) in which it still considers the location to be inside the cell, while our new algorithm does not. Maybe this is the difference why some particles with the new algorithm are deleted as outside of the domain, while with the old algorithm they would not?

I agree that maybe a different search scheme could help (like checking not only all cells that use the closest vertex, but also all neighbor cells of these cells, or we could also check all cells whose bounding boxes overlap with the particle), but maybe the simpler test would be to add that tolerance to the check, i.e. change

dealii/source/particles/particle_handler.cc

Line 1412 in 040e788

if (GeometryInfo<dim>::is_inside_unit_cell(

into

              if (GeometryInfo<dim>::is_inside_unit_cell(
                    reference_locations[0], 1.e-10))

[blaisb]

@gassmoeller Please don't be sorry. It happens to the best of us and this is part of working with a community code right. These things are bound to happen :)! If

The particle that is disappearing is not in a cell which has a boundary, it is one layer above. They are also quite far from the outer wall of the domain. I think the issue is more related to the curvature of the cell, because this issue does not appear for meshes which do not have a curved manifold. I don't have any good hypothesis, but this is not a good behaviour.

I will gladly test what you gave me and report back the results. This small unit test reproduces the bug easily.

[blaisb]

@gassmoeller I tested setting line 1412 to :
if (GeometryInfo<dim>::is_inside_unit_cell( reference_locations[0], 1.e-10))

And the particle still disappear.

Looking at it very closely, the particle is not in the last layer of cell before disappearing..
Step-158 :

Step-159 :

I am trying to debug this, but I'm not exactly sure where to look.
Is there a way we could revert the patch temporarily and work on a PR together with additional tests? I don't think this is a normal behaviour. Particles within the domain should not disappear. That particle is very far away from outer walls of the domain.

Take care :)
Bruno