Cell search strategies: reducing compute_point_locations search horizon?

[Benjamin Bercovici]

Dear all, 

I am working on a particle-in-cell code that uses the function compute_point_locations(P,whole_mesh) to find the cell to which each particle belongs. 

• P is a shortcut for the particle passed as an argument, and means either the particle or its location in what follows.
  
• whole_mesh is an instance of the class dealii::DofHandler<2> that stores the dofs of the entire mesh.

This solution works but is a bit slow (I have about 60,000 quads in my geometry). I am assuming that compute_point_locations(P,whole_mesh) does not have any a-priori knowledge of the location of P, and therefore loops over the entire mesh until it finds a cell that contains P. 

I was therefore trying to devise an alternative, for which I would need to have access to the neighbors of each cell.

1. At t=0, for each particle P
1. compute_point_locations(P,whole_mesh) is called and searches the entire mesh for P.
2. An iterator to the cell containing P at t=0 in then stored (as an attribute of the Particle class). Let's call this iterator it_cell_init
3. I would create a new DoFHandler object now restricted to the dofs of the cells surrounding *it_cell_init. A little bit of programming magic would be needed here, but let's assume that I can come up with a dealii::DofHandler<2> neighboring_mesh that does the job.
4. The particle then undergoes a displacement.
2. for each t>0, for each particle P
1. compute_point_locations(P, neighboring_mesh) could then be called, and would be expected to perform much faster since they would be only a limited number of inspected cells.
2. neighboring_mesh is updated to reflect the neighboring of the cell where P is in.
3. The particle is moved
4. ...

This is a little naive because 

1.  I do not have much details about compute_point_locations, and there might be a better way to solve my problem...
2.  I don't know if it is possible to extract the neighbors of a given cell and construct neighboring_mesh while keeping track of the correspondence between the original cell numbering and the resulting numbering of the restricted mesh.
   
3. Since my particles are moving around, I must make sure that 
1. The neighboring area in which they are searched for in is large enough, otherwise P won't be in the neighboring_mesh
2. It's not too large either

Anyways, let me know if this could be a valid way to help compute_point_locations to perform faster.

Thank you, 
Benjamin Bercovici

http://bbercovici.github.io/

[Matthias Maier]

Am 23. Sep 2014, 20:39 schrieb Benjamin Bercovici <benjamin....@gmail.com>:

> Dear all,
>
> 2. I don't know if it is possible to extract the neighbors of a given

> cell and construct neighboring_mesh while keeping track of the
> correspondence between the original cell numbering and the resulting
> numbering of the restricted mesh.

> 3. Since my particles are moving around, I must make sure that
> 1. The neighboring area in which they are searched for in is large

> enough, otherwise P won't be in the neighboring_mesh

> 2. It's not too large either

If the number of particles is relatively small it is a good idea to
store the previous location in form of a cell iterator and construct a
patch with something like the following as search candidates.

Have a look at cell->neighbor(i) [1]. You can do something like

57 std::set<CELL_IT> patch_around_cell;
58
59 patch_around_cell.insert(cell);
60
61 for (unsigned int i = 0; i < depth; ++i) {
63 for (unsigned int j = 0; j < dealii::GeometryInfo<dim>::faces_per_cell;
64 ++j) {
65
66 std::set<CELL_IT> tmp;
67 for (auto it = patch_around_cell.begin(); it != patch_around_cell.end(); ++it) {
68 if (!(*it)->face(j)->at_boundary())
69 tmp.insert((*it)->neighbor(j));
70 }
71 patch_around_cell.insert(tmp.begin(), tmp.end());
72 }
73 }

To get a patch around a cell (assuming a uniformly refined mesh and
without complex geometry - otherwise the iteration has to be a bit more
sophisticated).

You can always fall back to something more expensive if this fails.

Furthermore, it is likely that you need a small time step anyway, so the
particle should be quite close to the old cell (i.e. a patch with
depth==1 may already be sufficient).

Best,
Matthias

[1] http://www.dealii.org/8.1.0/doxygen/deal.II/classDoFCellAccessor.html#ab3fead7b17d7e990f40bbeeecbbdf3e4

[Benjamin Bercovici]

Matthias, 

this is exactly what I want to do.

The non-uniformity of the mesh might be a problem. If a more expensive scheme is needed, I don't think that it can be any worse than looping over 60,000 cells for nothing!

[Benjamin Bercovici]

And before I forget, 

Thank you!

I'll post the results here when i'll be done.

[Benjamin Bercovici]

Dear Matthias, 

I have adapted the code snippet that you sent, and I am now able to create a patch out of the neighbors of a given cell.

The only missing brick is the function which, when given a point (a set of two coordinates), searches the patch for the cell containing the point and returns the iterator to this cell (under the assumption of a neighborhood that it big enough to contain the point)

Actually, this would just boil down to this: a function associated to a cell which takes the point as an argument, and returns true/false whether the point belongs/does not belong to the cell.

I would either need such a function, or find a way to construct a DoFHandler object from the patch, and provide it to compute_point_locations().

I don't know if such a magic function exists ( a naive research for cell->contains did not return any result). 

What would you suggest?

Thank you, 

Benjamin Bercovici

Le mardi 23 septembre 2014 14:11:53 UTC-5, Matthias Maier a écrit :

[Luca Heltai]

Benjamin,

this problem is very relevant also for some of the things we have developed in the past.

The current implementation of find_active_cell_around_point is thought for very large coarse grids. If this is your case, then I think what has been discussed already is a good solution.

An alternative is to use a quadtree-octree splitting of space, which has a log(n) cost, if your coarse mesh is very small (ideally, one square/cube as a starting grid).

When you look for a point, start from the mesh at level 0, transform the point to unit cell, and check if the point is inside the unit cell (this is what is done internally in one part of the find_active_cell_around_point) (skip this if your coarse mesh has only one cell...).

This will not give you the *active* cell, but the coarsest cell at level 0 that contains the point. Now cycle over the 2^dim children of that cell, and repeat.

When a point has moved out, do the opposite: go to the father of this cell, see if the point is still there. If yes, cycle over children (maybe excluding the one you already know), and you are done. Otherwise go to the father of that cell. *Never* query the function find_active_cell_around_point as a wild guess, since the way this function works is

1. Find closest vertex (cost = N vertices)
2. Loop over all cells, to see who this vertex belongs to (average cost = N cells/2)
3. Loop over all cells that contain the given vertex, and
3.a tranform point to unit cell
3.b check if point is inside unit cell
3.c iterate

An even better solution would be:

1. Find closest *center* in coarse mesh (cost = N vertices **coarse mesh**)
2. Find closest *center* in children, until you get to the active child (cost = N levels * 2^dim)
3. 3.a, 3.b

If your coarse mesh is not coarse at all, create an external cartesian octree grid (starting from an hypercube containing your full domain), and use this as a "support mesh" to make your search algorithms efficient...

Let me know if you need help with this...

Luca.

--
Luca Heltai <luca....@gmail.com>
http://people.sissa.it/~heltai/
Scuola Internazionale Superiore di Studi Avanzati
Phone: +39 040 3787 449, Office: 622
--
There are no answers, only cross references

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Benjamin Bercovici]

Luca, 

thank you for your response. 
I had already started working on another solution when I saw your mail.

Here is what I have come up with:

• compute_point_location is still used at the initial time step to locate the cell where each particle is at t=0.
  
• A patch of cells surrounding this initial_cell is built and stored as an attribute by each particle.
• then, particles are moved to different cells. In order to find the cell which now contains each particle, I have come up with the following function Particle::find_current_cell

void Particle::find_current_cell(unsigned int depth){

    bool cell_found=false;

    for (auto it = patch_around_cell.begin(); it != patch_around_cell.end(); ++it){

      if ((*it)->point_inside(pos)){

        current_cell=*it;

        old_cell=current_cell;

        cell_found=true;

        build_patch(1);

        break;

      }

    }

    if(!cell_found){

      ++depth;

      build_patch(depth);

      this->find_current_cell(depth);

    }

  }

Note that find_current_cell() can be called recursively, if the function fails to find the particle within the immediate neighbors (the smallest value of depth is 1). If so, a bigger patch of cells is created (as depth is increased by one) and the search continues through this new patch, until the particle is found.

• when the particle is found, the cell it belongs to (is it what you call the "active" cell? I have to say I do not really see the difference) becomes the center of a new patch of depth 1, and so forth.

Searching for "active" cells this way, instead of relying on compute_point_locations at each iteration, seems much more effective: 

• At startup, compute_point_locations called on each particle (2000 particles, mesh comprised of ~ 60,000 quadrangles) takes 50s to run. 
• distributing the charge of each particle (i.e, finding the cell they belong to, and performing some weighing) takes about 5s to run.

In the previous version of my code, compute_point_locations was called at each timestep.

Further improvements could be made though (for instance, the search starts from scratch when a bigger patch is created which is very time consuming since a majority of the cells are known to be empty...). Another thing to keep in mind is that the depth of the patches can vary as the computation progresses so they might swell a little. But I don't think it can be any worse than before.

What is your opinion on this approach?

[Benjamin Bercovici]

Dear all, 

I have run several testcases to measure the performance improvement of patch searching over global searching.

50+50 computational cycles (each cycle equals 4 consecutive steps: distributing the space charge - > solving Poisson's equation -> computing E -> moving the particles) with ~60,000 cells in the mesh take (for a population of 2000 particles)

• 3374 s to run (wall clock time) with global searching at each timestep
• 1653 s to run (wall clock time) with global searching at the initial timestep, and patch searching at each subsequent timestep.

Attached to this message is the average patch depth at each timestep. Namely, if there were only two particles and if the depth of patches allowing the particles to be tracked was (1,1), (2,4) and (2,3), the average depth of  timestep # 1,2,3 would be 1, 3 and 2.5 respectively.

There appears to be no "explosion" of the average patch depth. Ions move a little farther away over one timestep than electrons do (the timestep used in the integrator scheme is not the same for ions and electrons), but the patch depth seems bounded...

However, this is not as satisfying as it looks. I ran another test case with 4000 particles for 400+400 timesteps, and the computation got stuck somewhere before the end of the first 400 iterations which deal with to electron motion only (ions are kept stationary during these 400 iterations, so it is impossible that they affect the computation speed at this point).

I am suspecting a brutal acceleration of the electrons motion to happen somewhere around 300 iterations. Which is physically acceptable but could nevertheless cause the patching method to fail. 

Yet, if the patch depth is indeed going to "infinity" (according to the implementation of find_current_cell hereabove), I would rather expect the patch to grow bigger until it includes the cell which contains the particle.

I'll keep thinking about it.

Benjamin 

[Bruno Turcksin]

Benjamin,

I am suspecting a brutal acceleration of the electrons motion to happen somewhere around 300 iterations. Which is physically acceptable but could nevertheless cause the patching method to fail. 

Yet, if the patch depth is indeed going to "infinity" (according to the implementation of find_current_cell hereabove), I would rather expect the patch to grow bigger until it includes the cell which contains the particle.

what could happen it's that you may spend a lot of time building bigger and bigger patch. If you need to build 10 patches for each electron, it could be slow. Maybe if the final patch for the previous particle had a depth of 10, you can start with a patch of depth 9 for the new particle.

Best,

Bruno

[Benjamin Bercovici]

Dear all,

I fixed a number of issues and the patching method is now operational.
It works significantly faster than the global search I had implemented before, and I think that I can still improve the code.

More specifically, each patch around a particle only contains the cell where the particle was found at the last iteration.
If the particle is no longer inside this cell, the patch grows until it eventually contains the particle.

The issue here is that when the patch grows in size, it is comprised that are mostly empty ( since they already were in the old patch)

Therefore, I would need a way to compare the old patch and the new patch in order to derive their difference ( in the same meaning as in std::set_difference), so that their difference new_cell can be searched for the particle.

However, a naive implementation of std::set_difference (old_patch.begin (), old_patch.end (), new_patch.begin (), new_patch.end (), std::inserter (new_cell)) did not yield the expected result, as new_cells remained empty.

A patch is defined as a std::set <dealii <DoFHandler<2>::cell_iterator > ( see line 410 of the linked page)

Is there an alternative way to extract the new cells from the knowledge of the old_patch and the new patch? If I had a way to compare two cell_iterators, I could come up with something doing the job. But I would prefer to use std::set_difference instead...

Thank you,
Benjamin Bercovici

Ps: you can find the implementation of Particle::build_patch () at the following location https://github.com/bbercovici/uiuc-cpp/blob/master/pic7.0/mycode/source/Particle.cc

[Bruno Turcksin]

Benjamin,

On 09/29/2014 04:32 PM, Benjamin Bercovici wrote:
> The issue here is that when the patch grows in size, it is comprised that are mostly empty ( since they already were in the old patch)
>
> Therefore, I would need a way to compare the old patch and the new patch in order to derive their difference ( in the same meaning as in std::set_difference), so that their difference new_cell can be searched for the particle.
>

you could create the new patch with only the new cells using deal.II
user_flag
http://dealii.org/8.1.0/doxygen/deal.II/DEALGlossary.html#GlossUserFlags
Once a cell is added to a patch, you set the flag. When you build the
next patch, you accept only the cells that have the flag set.

Best,

Bruno

[Benjamin Bercovici]

Bruno,
Thank you so much for pointing this out!

It works great right now.

10,000 particles over 100 timesteps take less than half an hour to run.

Attached are a few plots showing the patch depth, full patch size and sparse patch size averaged over the particle population at each timestep. The simulation itself was not physicaly accurate, but it was not the point here.

Thanks to all of you!

[Timo Heister]

I haven't read through the whole thread (sorry), but my two cents:

If you need to repeatedly find the cells a point belongs to, it is
worth constructing a structure to accelerate the process. Bounding box
trees come to mind (google for it). I think an AABB tree is the right
thing to do. With this you can identify a collection of a few cells
very very rapidly that you can then test using the expensive but exact
method.

This is something that might be useful for other users, so it would be
awesome if it made its way into deal.II.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
Timo Heister
http://www.math.clemson.edu/~heister/

[Wolfgang Bangerth]

Benjamin,
you've already mostly solved your problem, and I know I'm late to the
discussion, but I still have a suggestion:

> More specifically, each patch around a particle only contains the cell
> where the particle was found at the last iteration. If the particle is no
> longer inside this cell, the patch grows until it eventually contains the
> particle.

I assume once you find the cell you are looking for, then you just throw the
whole patch away and again only store that one cell? That is how I would have
implemented it as well.

The suggestion I wanted to make is that right now you're first building the
patch by including all neighbors of cells in the patch into the set (the
operation automatically throws away duplicates), and only then walking over
the cells to find which one the particle is in. That's inefficient. You should
of course terminate the construction of the patch once you've found the particle.

But even then, I think you can improve this. Instead of just walking over all
neighbors of all cells that are already in the patch, consider which
*direction* the particle has gone. If you know the old and the new position of
the particle, you can at least guess which neighbor it is in. This guess is
going to be correct most of the time, and if it isn't you can still search all
neighbors. What I want to say is that if you are smart about it, I'm sure you
can reduce the number of tests whether a particle is in a cell by a very
significant factor.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/