Système de Particules

Version imprimable

Voir 40 message(s) de cette discussion en une page

Bonjour,

Je suis actuellement entrain d'implémenter un système de particules dans un environement en 2D. Pour le moment le programme peut traiter environ 2 000 particules avant de commencer a ramer. Je commence a être a cours d'idée sur les optimisations que je peut faire pour augmenter le nombre de particules. J'ai donc plusieurs questions :

- Le conteneur choisis pour les particules vous semble t-il correct ?
- La répartition des particules dans une grille est-elle correct ?
- Est que une ligne de code vous parait bizarre ? :D

En vous remerciant d'avance pour vos aides :)

Le code:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
void ParticleSystem::computeParticles(Map *pMap)
{
    int widthCell = 32;
    int heightCell = 32;
 
    int nbCellWidth = pMap->WIDTH * 32 / widthCell + 1;
    int nbCellHeight = pMap->HEIGHT * 32 / heightCell + 1;
 
     std::vector<int> grid[nbCellHeight][nbCellWidth];
 
    for(int i = 0; i<mParticles.size(); i++){
 
        int x = mParticles[i]->mPositionY/heightCell;
        int y = mParticles[i]->mPositionX/widthCell;
 
        grid[x][y].push_back(i);
 
    }
 
    // trouver les particules voisines
    for(int i = 0; i<nbCellHeight; i++)
        for(int j = 0; j<nbCellWidth; j++)
        {
            while(!grid[i][j].empty()){
                int particule = grid[i][j].back(); // remove particle par la fin pour O(1)
                grid[i][j].pop_back();
 
                for(int x = -1; x<=1; x++){
                    for(int y = -1; y<=1; y++){
                        int a = i + x;
                        int b = j + y;
 
                        if(a < 0 || a >=nbCellHeight || b < 0 || b >= nbCellWidth)
                            continue;
 
                        for(int z = 0; z<grid[a][b].size(); z++)
                        {
                            int voisin = grid[a][b][z];
 
                            float dx = mParticles[voisin]->mPositionX - mParticles[particule]->mPositionX;
                            float dy = mParticles[voisin]->mPositionY - mParticles[particule]->mPositionY;
                            float distance2 = dx * dx + dy * dy;
                            if(distance2 < mParticleWidth * mParticleWidth){
                                float distance = sqrt(distance2);
                                float weight = (1.0f-distance/mParticleWidth);
                                ParticleContact contact(particule, voisin, distance, weight);
                                mParticleContacts.push_back(contact);
 
                            }
                        }
                    }
                }
                mParticles[particule]->mDensity=0.0f;
                mParticles[particule]->mPressure=0.0f;
                mParticles[particule]->mVelocityY -= GRAVITY;
            }
        }
 
    // compute density
    for(int i = 0; i<mParticleContacts.size(); i++){
 
        int p1 = mParticleContacts[i].mParticule1;
        int p2 = mParticleContacts[i].mParticule2;
 
        mParticles[p1]->mDensity += mParticleContacts[i].mWeight;
        mParticles[p2]->mDensity += mParticleContacts[i].mWeight;
 
    }
 
    while(!mParticleContacts.empty()){
 
        ParticleContact contact = mParticleContacts.back();
        mParticleContacts.pop_back();
 
        int p1 = contact.mParticule1;
        int p2 = contact.mParticule2;
 
        //compute pressure
        mParticles[p1]->mPressure = std::max(0.0f, n * (mParticles[p1]->mDensity - w0));
        mParticles[p2]->mPressure = std::max(0.0f, n * (mParticles[p2]->mDensity - w0));
 
        //compute force
        float nX = (mParticles[p2]->mPositionX - mParticles[p1]->mPositionX)/contact.mDistance;
        float nY = (mParticles[p2]->mPositionY - mParticles[p1]->mPositionY)/contact.mDistance;
 
        mParticles[p1]->mVelocityX += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * nX;
        mParticles[p1]->mVelocityY += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * nY;
 
        mParticles[p2]->mVelocityX += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * -nX;
        mParticles[p2]->mVelocityY += TIME_STEP * a * (mParticles[p1]->mPressure + mParticles[p2]->mPressure)
                                        * contact.mWeight * -nY;
 
        mParticles[p1]->mVelocityX += TIME_STEP * VISCOSITY * (mParticles[p2]->mVelocityX - mParticles[p1]->mVelocityX);
        mParticles[p1]->mVelocityY += TIME_STEP * VISCOSITY * (mParticles[p2]->mVelocityY - mParticles[p1]->mVelocityY);
 
        mParticles[p2]->mVelocityX += TIME_STEP * VISCOSITY * (mParticles[p1]->mVelocityX - mParticles[p2]->mVelocityX);
        mParticles[p2]->mVelocityY += TIME_STEP * VISCOSITY * (mParticles[p1]->mVelocityY - mParticles[p2]->mVelocityY);
 
    }
 
    // update position
    for (int i = 0; i < mParticles.size(); ++i) {
        mParticles[i]->updatePosition(pMap ,TIME_STEP);
    }
}

Les conteneurs dans ParticleSystem.h

Code:

1
2
std::vector<std::shared_ptr<Particle>> mParticles;
        std::vector<ParticleContact> mParticleContacts;

et la classe ParticleContact

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
class ParticleContact
{
    public:
        ParticleContact() {}
        ParticleContact(int p1, int p2, float distance, float weight) {
            mParticule1 = p1;
            mParticule2 = p2;
            mDistance = distance;
            mWeight = weight;
 
        }
        virtual ~ParticleContact() {}
 
        int mParticule1;
        int mParticule2;
        float mDistance;
        float mWeight;
};

17/08/2017, 10h52
ternel

Première réaction: Pourquoi un vecteur de shared_ptr? Ca te tue les performances car les données ne sont plus contigües en mémoire.
17/08/2017, 11h25
oodini

Je n'ai pas le temps de rentrer dans le détail de ton code, notamment l'algorithmique, mais utiliser des push_back, pour les perfs, ce n'est pas génial.
Pour ce genre de trucs, il vaut mieux coder dans le style C. Cela te permettra en outre de plus facilement faire un portage en CUDA/OpenCL, qui sont des technos bien adaptés à ce genre de problèmes.
17/08/2017, 11h42
dragonjoker59
Salut !

Code:

std::vector<int> grid[nbCellHeight][nbCellWidth];

Ce conteneur n'est pas cache friendly, n'est pas valide (on ne peut pas déclarer un conteneur statique avec des tailles dynamiques), et comme je suppose que ta map ne doit pas changer de taille souvent, tu devrais le mettre en membre de la classe, pour l'initialiser une fois pour toutes, plutôt que de faire les allocations/désallocations qu'il implique, à chaque calcul de tes particules.
Je l'écrirais comme suit (pour l'instant) :
Code:

1 2 3 Map const & mMap; using Particles = std::vector< int >; std::vector< Particles > mGrid;
et dans le constructeur de ParticleSystem :
Code:

1 2 3 4 5 ParticleSystem::ParticleSystem( Map const & map ) : mGrid( map.WIDTH * map.HEIGHT, Particles() ) , mMap( map ) { }
Il faudra donc, pour te faciliter la vie, écrire les accesseurs te permettant d'accéder à l'élément aux coordonnées X et Y :
Code:

1 2 3 4 Particles & getParticles( size_t x, size_t y ) { return mGrid[x * mMap.HEIGHT + y ); }
Ensuite, si tu connais le nombre de particules max que peut contenir Particles, tu peux le définir comme suit, afin de gagner encore au niveau du parcours du tableau :
Code:

1 2 size_t constexpr MaxParticles = 12; //par exemple using Particles = std::array< int, MaxParticles >;
Ca te permettrait de n'avoir qu'une seule allocation en création du ParticleSystem, le reste n'étant que du parcours.
De plus, cette linéarisation des données te permettrait de supprimer la double boucle au niveau de la recherche des voisins (les double boucles c'est mauvais pour les perfs) :
Code:

1 2 3 4 for ( auto & particles : mGrid ) { // Recherche voisins... }
Bon, mais tout ça, ça montre une chose : tu as mal choisi tes structures de données.
Parce que ta recherche de voisins consiste en 6 BOUCLES IMBRIQUEES et si 2 boucles imbriquées c'est mal, je te laisse imaginer avec 6.
Ton soucis de perfs vient exactement de là, car au pire des cas, tu te retrouves avec une complexité de O(N^6) :aie:.
17/08/2017, 15h53
guppyworld

Citation:

Première réaction: Pourquoi un vecteur de shared_ptr? Ca te tue les performances car les données ne sont plus contigües en mémoire.

Oui c'est en faisant des test j'ai oublié de remettre comme avant. Donc pour le moment ce sont des unique_ptr.

dragonjoker59:

Donc en faite si j'ai bien compris tu veut que pour la grille contenant l'indice des particules, je passe d'un tableau 2D à un tableau 1D ?

Et qu'il est préférable d'allouer directement un tableau pouvant contenir 10 000 particules par exemple plutôt que de faire un tableau dynamique ?

Citation:

Parce que ta recherche de voisins consiste en 6 BOUCLES IMBRIQUEES et si 2 boucles imbriquées c'est mal, je te laisse imaginer avec 6.
Ton soucis de perfs vient exactement de là, car au pire des cas, tu te retrouves avec une complexité de O(N^6) .

La j'ai un doute sur le O(N^6) car les 2 premiers for ne dépendent pas du nombre de particules mais plutôt d'un nombre constant. Pareil pour les 2 for avec x et y qui sont également constant.
17/08/2017, 15h59
Bousk

Des unique_ptr fragmentent toujours ta mémoire.
17/08/2017, 16h08
ternel
En fait, la technique est même plutot recommandée, parce que tu limites les collisions aux collisions d'un même sous-espace, chacun contenant en général beaucoup moins de particules.
Dans un espace quadrillés en B blocs, pour N particules, tu ne calcules que B * (N/B)² = N²/B au lieu de N². À condition d'avoir une répartition uniforme dans la grille.

C'est en partant sur un tel raisonnement que tu arriveras sur le quad-tree (et son cousin en 3D, l'oct-tree), qui a pour vocation d'assurer une bonne distribution dans la grille.

Ta grille peut tout à fait être linéarisée.

D'un point de vue algorithmique, la boucle pourrait être pour toute particule p, pour toute autre particule q de la meme case, calculer l'effet de la collision de q sur la p.

Comme "pour toute particule p" est mathématiquement équivalent à "pour toute case c, pour toute particule p de c", c'est ce que tu as fait.

Par contre, tes calculs m'ont l'air très lourds, et surtout, ton assignation de case.
Je ne pense pas qu'il soit nécessaire de faire des allocations (et réallocation) à chaque pas de temps.

Avec un peu de "mise en cache", ta particule pourrait ressembler à
Code:

1 2 3 4 5 6 struct particule { coord_type coord; //une structure {x, y} vector_type vitesse; // une structure {dx, dy} vector_type acceleration; masse_type masse; };
La raison pour laquelle je différencie systématiquement coord_type de vector_type, c'est que les opérateurs n'ont pas la même signification. Notamment coord_type + coord_type n'a pas de sens.

Et à chaque pas de temps, ton code devrait globalement être:
```
pour chaque pas de temps
    pour chaque particule,
        calculer la nouvelle accélération
            remise à 0
            calcul des collisions
        appliquer l'accélération,
        puis la vitesse.
```
éventuellement, tu peux stocker deux vector/array de particules, au début identiques. A chaque pas de temps, tu utilises l'un comme le temps actuel, et l'autre comme le temps suivant. Et au pas suivant, tu réutilises l'ancien temps actuel comme nouveau temps futur.
Il suffit de prendre les deux listes en références sur ta fonction effectuant un pas de temps.
17/08/2017, 16h28
ternel

Par contre, quel symptome qualifies-tu de "ramer"?

Combien as-tu de mémoire vive d'installée sur ton pc?

Normalement, avec un vector sans pointeurs, tu devrais avoir une consommation mémoire de plus ou moins une cacahuète (les quelques constantes de ton programme) + une cacahuète (la pile pour atteindre ta fonction de pas de temps) + sizeof <Particule> * N octets + N² * sizeof<Collision>.

Et cette dernière partie peut couter cher. N² sizeof collision, ca peut rapidement monter si tes particules sont dans un espace petit pour leur vitesse.
Il n'est pas nécessaire de stocker chaque collision, elle peuvent s'additionner au fur et à mesure.

Citation:

Par contre, tes calculs m'ont l'air très lourds, et surtout, ton assignation de case.

Tu veux parler des lignes de code suivante ?
Code:

1 2 3 4 5 6 7 8 for(int i = 0; i<mParticles.size(); i++){ int x = mParticles[i]->mPositionY/heightCell; int y = mParticles[i]->mPositionX/widthCell; grid[x][y].push_back(i); }
Sinon pour le moment je ne compte pas utiliser le QuadTree, je reste sur une carte divisé en case de 32 pixels sur 32.

Par contre je sais pas si vous avez compris mais en faite la grille contient les indices des particules pour chaque cellule. Je ne vois donc pas trop comment je peut faire autrement que la structure

Code:

std::vector<int> mGrid[NB_CELL_HEIGHT][NB_CELL_WIDTH];

. Ou au mieux je peut faire

Code:

std::vector<int> mGrid[NB_CELL_HEIGHT * NB_CELL_WIDTH];

.

17/08/2017, 18h09
ternel

Tu recalcules à chaque itération l'intégralité de ces listes.
Je te rappelles que vector est sucesptible de réallouer son tableau interne à chaque allocation.
Donc tu passes ton temps à faire des recopies de tableaux (certes d'entier) à chaque itération.

Alors que tu pourrais qu'avec deux grilles, tu n'aurais pratiquement plus ce probleme.
Utilise deux "grid", prises en argument, et tu devrais déjà récupérer pas mal.

De plus, ce que je te suggère, c'est de tout de suite appliquer tes collisions, pour ne pas avoir à les stocker.
17/08/2017, 18h42
guppyworld

Citation:

Tu recalcules à chaque itération l'intégralité de ces listes.

Je suis obligé car une particule peut changer de cellule entre 2 itérations.

Citation:

Je te rappelles que vector est sucesptible de réallouer son tableau interne à chaque allocation.

Je ne vois pas quoi d'autre d'utiliser qu'un vecteur car je ne sais pas combien de particules vont etre dans la cellule.

Concernant l'utilisation des pointeurs intélligents, je les oubli et utilise uniquement le new et delete pour les particules ?
18/08/2017, 10h41
Bousk

Non. new/delete, unique_ptr, shared_ptr, tout ça c'est du pointeur et ça fragmente la mémoire, du coup tu perds quasi tout l'intérêt du vector et la contiguité mémoire. Bonjour les cache-miss et les perfs plombées.
- pourquoi stocker la collision ?
- pourquoi trier les particules par case ?
-> pourquoi les particules ne sont pas dans leur collection globale avec une simple position et affichées tel quel ?

L'idée générale de la grille est assez bonne, j'ai expliqué pourquoi plus haut.
Ce qui ne l'est pas, c'est de la recréer à chaque tour.

C'est la raison pour laquelle je te conseille vivement deux grilles, les réallocations seront beaucoup moins nombreuses.

Tu prends la première particule.
tu réinitialises son accélération à la valeur infligée par l'environnement extérieur (la gravité par exemple)
tu cherches si chaque autre particule (pas trop lointaine) la percute
le cas échéant, tu ajoutes a son accélération l'accélération infligée par cette collision
Lorsque tu as vu toutes ces autres particules, tu peux appliquer cette accélération à la particule, pour obtenir sa nouvelle vitesse, puis sa nouvelle position.
Ce faisant, tu mets à jour cette particule dans la liste "futur" (ainsi que sa position dans la grille si elle change).

Et tu passes à la suivante.

Le fait d'avoir deux listes (et deux grilles) permet de ne pas mélanger le futur et le présent.

Quant à ton code appelant cette fonction, il devrait ressembler à :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 etat simulation(int n) { etat base = initialiser(); etat autre = present; //juste histoire d'initialiser tous les conteneurs, les valeurs en elles etat *e1 = present; etat *e2 = autre; for (int i = 0; i < n; ++i) { simuler_un_temps(*e1, *e2); std::swap(e1, e2); } return (n%2==0) ? *e1 : *e2; }
Ca manque de noms correctes, mais ca résume l'utilisation.
Il te faut bien sur une classe état utile, qui contiennent à la fois la description en grille et la liste des particules.

De toute façon, tu as besoin d'utilitaire pour clarifier ton code: classe de point, de vecteur, de particule, de grille, de liste (std::vector pouvant suffire), d'état, de force, éventuellement distinguer le vecteur vitesse, le vecteur accélération, un type de temps, etc.

Citation:

pourquoi stocker la collision ?

J'ai fait ca pour éviter de devoir refaire les calculs de distance et poids entre les particules.

Citation:

-> pourquoi les particules ne sont pas dans leur collection globale avec une simple position et affichées tel quel ?

Dans ParticuleSystem.h je l'ai déja fait non ?

Code:

std::vector<std::unique_ptr<Particle>> mParticles;

Je suis obligé d'alloué dynamiquement les particules si je ne veut pas qu'elles soit détruites à la fin de la portée il me semble.
Code:

1 2 3 4 5 void ParticleSystem::createParticle(int type, float positionX, float positionY, float velocityX, float velocityY, int lifeTime, float width) { std::unique_ptr<Particle> particle(new Particle(type, positionX, positionY, velocityX, velocityY, lifeTime, width)); mParticles.push_back(std::move(particle)); }
Je vais essayer de faire avec 2 états pour voir.

18/08/2017, 16h56
dragonjoker59

Si ton intention n'est pas d'utiliser le polymorphisme sur tes particules, alors les pointeurs sont contre productifs au possible.
Ma devise : "Ne jamais utiliser les pointeurs, sauf si j'ai une bonne raison"

Un vecteur de unique_ptr ou un vecteur d'objets, c'est identique du point de vue de la portée.
Comme je le dis dans ma signature, "un pointeur de moins, c'est une montagne d'erreurs en moins".

Que craignais-tu précisément?
Sans pointeur, tu crées une particule temporaire, et tu la passe en argument à push_back. Cette fonction crée dans le vecteur une copie de la particule fournie.
Ta variable temporaire sortant de portée, elle est détruite, mais pas la copie dans le vecteur.

Tu te doutes bien que si ce n'était pas le cas, vector ne servirait pratiquement jamais.

Par contre, une fois le pointeur supprimé, ta fonction peut être légèrement simplifiée avec emplace_back.
Code:

1 2 3 4 void ParticleSystem::createParticle(int type, float positionX, float positionY, float velocityX, float velocityY, int lifeTime, float width) { mParticles.emplace_back(type, positionX, positionY, velocityX, velocityY, lifeTime, width); }
Remarque comme je n'ai pas besoin de créer de variable temporaire.
En fait, emplace_back appelle directement le constructeur "normal" de Particule, plutot que le constructeur de copie.

18/08/2017, 18h00
guppyworld

OK merci pour les précisions.

En effet il n'y a pas besoin de pointeurs intelligents et cela permet de passer de 2000 particules max a 2600 :)

Par contre le principe des deux grilles j'ai pas trop compris l'utilité et comment le mettre en place puisqu'il faut à chaque iteration déplacer les particules dans des cellules differentes. Donc si on fait comme tu le suggère une copie de la grille entre 2 itérations il va falloir faire des removes au milieu du vector si il faut deplacer la particules vers une autre cellule =/
18/08/2017, 20h40
ternel

Si l'objet vector en lui-même n'est pas détruit, un vidage puis remplissage ne coute pas de réallocation, puisque la mémoire est déjà là.

Oui, tu vas la reremplir, mais ca coutera beaucoup moins de temps que de tout recommencer.
Autre question, pourquoi ta grille est-elle mesurée en case de 32 pixels? le pixel est une unitée d'affichage, qui n'a pas lieu d'être l'unité de position.
18/08/2017, 20h51
guppyworld

Citation:

Autre question, pourquoi ta grille est-elle mesurée en case de 32 pixels?

J'ai pris cette valeur car cela permet d'avoir peu de particules par cellule.

Après j'aurais peut être dû prendre des cellules de 20x20 car le diamètre des particules est de 20.
19/08/2017, 11h26
ternel

En fait, la "bonne" valeur, c'est en fonction de la taille de l'espace disponible, et du nombre de particules.
Théoriquement, tu peux calculer une distance moyenne de parcours sans collision, et adapter ta taille de cellule en fonction.
En pratique, ca devrait être parmi les paramètres de ton modèle.

Cela dit, ma question portait plus sur "pixels" que sur 32.

Voir 40 message(s) de cette discussion en une page