Clearing a particle properties without getting properties for each particle

[blais...@gmail.com]

Dear all,

I hope you are well :).

We are currently facing a weird bottle neck. There are some operation which require us to set a particle property to zero and I think we are doing something which is not optimal at all. Right now, we loop over the particles, then we get the properties array view and set the property to zero.

I was wondering if there was a more optimal way to do this? I have profiled it with callgrind and it seems that getting the properties array view from the particle can be really expensive.

Since what we want to achieve is just to set one of these properties to zero, I was wondering if there was not a possible optimisation? For example, could we just get the array view once and assuming its continuity, set all of the same property to zero?

Sorry if the question seems dumb, I am really less familiar with this type of data structure.

Thank you so much!
Bruno

[Wolfgang Bangerth]

> I was wondering if there was a more optimal way to do this? I have profiled it
> with callgrind and it seems that getting the properties array view from the
> particle can be really expensive.
> Since what we want to achieve is just to set one of these properties to zero,
> I was wondering if there was not a possible optimisation? For example, could
> we just get the array view once and assuming its continuity, set all of the
> same property to zero?

That seems strange to me. The get_properties() function is in essence a
one-liner with a few indirections. Can you show the code in which you set
properties to zero?

Best
W.

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

[blais...@gmail.com]

Maybe it is callgrind that is playing tricks on me.

Here is the code :

for (auto particle = particle_handler.begin();
particle != particle_handler.end();
++particle)
{
// Getting properties of particle as local variable
auto particle_properties = particle->get_properties();

// Reinitializing forces and momentums of particles in the system
particle_properties[DEM::PropertiesIndex::force_x] = 0;
particle_properties[DEM::PropertiesIndex::force_y] = 0;

particle_properties[DEM::PropertiesIndex::M_x] = 0;
particle_properties[DEM::PropertiesIndex::M_y] = 0;

if (dim == 3)
{
particle_properties[DEM::PropertiesIndex::force_z] = 0;
particle_properties[DEM::PropertiesIndex::M_z] = 0;
}
}  

With dim a template parameter of the class.

The index for the particle properties are part of an enum:

enum PropertiesIndex : int
{
type = 0,
dp = 1,
rho = 2,
v_x = 3,
v_y = 4,
v_z = 5,
acc_x = 6,
acc_y = 7,
acc_z = 8,
force_x = 9,
force_y = 10,
force_z = 11,
omega_x = 12,
omega_y = 13,
omega_z = 14,
mass = 15,
mom_inertia = 16,
M_x = 17,
M_y = 18,
M_z = 19,
displacement_x = 20,
displacement_y = 21,
displacement_z = 22,
n_properties = 23,
};  

I was not sure on what would be the best data structure to identify the data in the particle properties, hence this type of enum (which I know is maybe not ideal but...)

Do you have any suggestions? Could it just be the cost of iterating through the map that callgrinds wrongly affects to the zeroing of the variables?

[Wolfgang Bangerth]

The code looks ok. I have no idea why that shows up so highly in your
profiles, but would check this with a TimerOutput section like we use in other
tutorial programs. This also tells you how often the code is actually called.
I suspect that valgrind is giving you inaccurate information.

[blais...@gmail.com]

The issue is that this code is being called a lot (in our case every iteration, so 100k to 10M times) and each iteration is very fast (1k iteration is a second or so). Consequently, we've had some severe issue when timing in parallel.

If nothing is wrong with the code I guess what I am seeing is just the cost of iterating through the particle data structure and valgrind is confusing me :). I'll try to see if I can combine this loop with another one.

Thank you very much for the time and the help!

[Wolfgang Bangerth]

On 10/19/20 4:58 AM, blais...@gmail.com wrote:
> The issue is that this code is being called a lot (in our case every
> iteration, so 100k to 10M times) and each iteration is very fast (1k iteration
> is a second or so). Consequently, we've had some severe issue when timing in
> parallel.
> If nothing is wrong with the code I guess what I am seeing is just the cost of
> iterating through the particle data structure and valgrind is confusing me :).
> I'll try to see if I can combine this loop with another one.

That's one option. Option 2 is to write an interface to the PropertyPool class
that allows you to set to zero certain components of all particle properties.

[blais...@gmail.com]

I am not sur if I fully understand the memory structure of the Property Pool (or at least I am still confused by it).

If I understand correctly, the properties of all the particles are stored in a single Handle array and the get_properties of a single particle returns an array view pointing to this large array.

However, the n_properties_per_slot() returns the number of properties associated with a single particle. Is there a way to know from within the property pool the size of the handle array?

If the properties are stored in a continuous block, one could easily write a function that takes into argument a vector of property index, the value to set and the number of particles for example. This function would then be called from the particle handler which in turn would have an interface to this function with the vector of property index and the value to set.

Would that be the right approach? If so, I can work on a PR.

[Wolfgang Bangerth]

On 10/19/20 9:03 AM, blais...@gmail.com wrote:
> I am not sur if I fully understand the memory structure of the Property Pool
> (or at least I am still confused by it).
> If I understand correctly, the properties of all the particles are stored in a
> single Handle array and the get_properties of a single particle returns an
> array view pointing to this large array.
> However, the n_properties_per_slot() returns the number of properties
> associated with a single particle. Is there a way to know from within the
> property pool the size of the handle array?

Maybe
https://github.com/dealii/dealii/pull/11057
clarifies the intent of the PropertyPool? The idea is that the property pool
is responsible for all of the memory allocated for particle properties.

I've come to realize that what you want to do wasn't possible in the current
design (because we don't keep track of which memory we have allocated), so wrote
https://github.com/dealii/dealii/pull/11059

Once that is merged, I think what you want to do is maybe something along the
lines of the following addition to the PropertyPool class:

/**
* Return the first of all of the memory pool handles that are managed
* by this class. Iterating from begin() to end() then returns the handles
* associated with all particles in the ParticleHandler class that owns
* this property pool; using the get_properties() function, these handles
* can be converted to ArrayView objects that then allow access to the
* actual properties.
*
* This function can be used to access and modify the properties of all
* particles, though possibly in a different order than if one accessed
* them using ParticleHandler::begin() and ParticleHandler::end().
*
* @note PropertyPools also store the properties of ghost particles,
* and consequently iterating over all handles also reaches the
* properties of ghost particles.
*/
std::set<Handle>::iterator begin ()
{ return currently_open_handler.begin() }

+ same for the corresponding end() function.

> If the properties are stored in a continuous block, one could easily write a
> function that takes into argument a vector of property index, the value to set
> and the number of particles for example. This function would then be called
> from the particle handler which in turn would have an interface to this
> function with the vector of property index and the value to set.
> Would that be the right approach? If so, I can work on a PR.

PropertyPool doesn't advertise how it stores properties. Currently, it's not
in one large block, though that's something we might want to change at some
future point. The interface above avoids exposing how exactly they are stored.

[shahab.g...@gmail.com]

Hi,

Can you please write an example of accessing memory pool handles, iterating them and converting the handles to particle properties?

Thanks in advance.

Best regards,

Shahab

[Wolfgang Bangerth]

On 12/3/20 1:19 PM, shahab.g...@gmail.com wrote:
> Can you please write an example of accessing memory pool handles, iterating
> them and converting the handles to particle properties?

The way we think about it is that the MemoryPool stores the data in some way
that should not be of any concern to you. You access these properties through
the particles they are attached to. As a user you don't have contact with the
PropertyPool at all.

At least that was the original idea. Can you explain what you're trying to do?

[shahab.g...@gmail.com]

Thanks for the reply.

I want to iterate through all the particles in each iteration of my solver and update some of the properties. The problem is iterating through particles using particle_handler() is rather computationally expensive. I am looking for a way to update the properties without iterating through the particle_handler.

Best regards,

Shahab

[Wolfgang Bangerth]

On 12/3/20 4:13 PM, shahab.g...@gmail.com wrote:
> I want to iterate through all the particles in each iteration of my solver and
> update some of the properties. The problem is iterating through particles
> using particle_handler() is rather computationally expensive. I am looking for
> a way to update the properties without iterating through the particle_handler.

I see.

So you'd need to iterate through the the blocks of properties stored in the
PropertyPool? We could write such an interface, but let me ask a couple of
questions first:
* Do you need to map back from block of memory to the owning particle?
* Will you somehow need to treat ghost particles differently than the locally
owned ones? How will you make sure that locally owned particles on one
processor stay in sync with the corresponding ghost particles on the other
processor(s)?

[shahab.g...@gmail.com]

Not here, but it other functions it would be very useful to be able to map back to the owning particle.

In fact, I have to treat local and ghost particles differently. In terms of updating properties, I only need to update the properties of local particles. When the maximum displacement of a particles reaches half of the cell length in the triangulation, I call sort_particles_into_subdomains_and_cells() to update the list of local and ghost particles. I can explain the algorithm in more details if you want.

Thanks again.

Shahab

[blais...@gmail.com]

Hello Shahab,

I would wait before trying to do some changes. Rene has started looking more deeply into the particle structure (see : https://github.com/dealii/dealii/pull/11314) and this will, most likely, fix your problem.

[Wolfgang Bangerth]

On 12/4/20 6:18 AM, shahab.g...@gmail.com wrote:
> Not here, but it other functions it would be very useful to be able to map
> back to the owning particle.

I suspect that in that case, it wouldn't be cheaper to loop over particles and
map back to particles, but that you could as well loop over particles to begin
with.

> In fact, I have to treat local and ghost particles differently. In terms of
> updating properties, I only need to update the properties of local particles.

This is where things get difficult. The PropertyPool doesn't know which
particle properties correspond to locally owned or ghost cells -- and I'd
argue that it shouldn't know that either. PropertyPool is really just a
slightly more intelligent way to allocate memory on the heap. What the
*meaning* of what it stores is not something it knows about, and I don't think
it should.

If you needed to do something different between locally owned or ghost
particles, you'd have to store information in duplicate in some way, and that
presents the potential for things to go out of sync.