Storage for particles

[chrisrichardson]

On each process, we need some storage for the particle data. Each particle should have the following:

• a cell (index) to which is it attached
• a position
• some values which it carries

All particles should have the same set of attributes (with different values).

Possible implementations:

struct Particle {
int cell;
Eigen::Vector3d x;
std::vector<double> values;
} Particle;

// Storage of particles
std::vector<Particle> particles;

// Offsets for each value in values
std::vector<int> value_layout;

(looks attractive, but requires new std::vector for each particle)

or

std::vector<int> cell_indices;
std::vector<double> particle_values;
std::vector<int> value_layout;

(should be more efficient).

How should the values be set up? Previously, the first value was always the position. Should there be a method to set up the values (e.g. with string names for debugging?) at construction time.

[chrisrichardson]

Looking at this more carefully, it is a bit difficult to decide the best way. In many cases, we want the particles collected by cell (e.g. for projection/interpolation), so we could have a method e.g. compute_cell_particles() to get that, or maybe it is better to just use std::vector<std::vector<int>> cell_to_particles for indexing them? Previously in original LEoPart, it was std::vector<std::vector<particle>> which probably resulted in some needless copying...

[chrisrichardson]

@nate-sime @jmmaljaars - any thoughts?

[jmmaljaars]

Seems that decisions here are pretty crucial, since we want a good balance between flexibility and efficiency. Just a few first thoughts in view of this:

• most of the time, a particle doesn't really have to be aware of its hosting cell. Rather, a process and a cell should be aware which particles it's hosting. So this pleads in favor of a std::vector<std::vector<int>> cell_to_particles. (Apart from the needless copying when a particle migrates from one cell to the other on the same process, I was happy with the choice std::vector<std::vector<particle>> in the original LEoPart)
• so at first sight it seems that std::vector<double> particle_values; is the way to go for storing particle properties.
• however, we want to conveniently access the different particle properties - in particular the position x - for the particle advection and updating of properties. Can we do something with a Particle struct here which just provides a view to the underlying particle_values vector?

Just an additional question:

• with std::vector<int> value_layout; (maybe name it std::vector<int> particle_layout;?): how are we going to enable tensor-valued properties? Can we say without loss of generality that in case the value of int is n_d**2 (with n_d>1 the geometric dimension) we are dealing with a tensor-valued particle property?

[chrisrichardson]

Yes, I think std::vector<std::vector<int>> can be OK - it is not really efficient, because the vector for each cell requires allocation. Imagine if we have 10^6 cells on a process... But I suppose, once it is done, it is done.

As for storing the actual values, I think we can have one large vector std::vector<double> particle_values. If we need to access into it, we can use Eigen::Map to pull out e.g. the position, or a tensor value etc. The only thing we might want to consider is whether it should be possible to use complex numbers. That could make it more complicated (as x is obviously not complex, but the other "values" could be).

The value_layout is also quite a small thing. It would be easy to e.g. make it a "class" and add more properties, like tensor shape for each item. I was wondering whether to attach a FunctionSpace for each value? But maybe it is better to keep all dolfin stuff a bit more outside.