Forum |  HardWare.fr | News | Articles | PC | Prix | S'identifier | S'inscrire | Aide | Shop Recherche
1508 connectés 

  FORUM HardWare.fr
  Programmation
  C++

  interaction à grand nombre de corps et quadtree (barnes hut)

 


 Mot :   Pseudo :  
 
Bas de page
Auteur Sujet :

interaction à grand nombre de corps et quadtree (barnes hut)

n°2149297
Profil sup​primé
Posté le 13-07-2012 à 22:47:41  answer
 

Bonjour à tous,
 
Mon projet
 
Je travaille actuellement sur un programme devant simuler l'évolution d'un système de n corps en interaction (gravitationnelle dans notre cas...) avec n très grand. (càd, minimum 1000 voire plus si possible sur mon PC :D)
Du coup la méthode qui consiste à calculer les n(n-1) interactions n'est pas envisageable (bien trop lente).
D'autre part je me fiche de la précision (enfin il faut que ça reste un peu réaliste quand même) le but étant d'avoir une idée assez grossière du comportement de grands ensemble de particules, voire de ce qui se produit lors de collisions entre deux amas.
 
La méthode retenue
 
Du coup, j'ai opté pour la méthode de Barnes-Hut (http://en.wikipedia.org/wiki/Barnes%E2%80%93Hut_simulation) qui est tout à fait adaptée à ce genre de calculs.
Le principe est de créer un arbre contenant chaque particule. Dans mon cas, dans un premier temps, je travaille en 2D, et j'ai donc à manipuler des quad tree.
On a donc un quad tree, avec un arbre principal (root), et un certain nombre de nœuds tels que chaque particule soit comprise dans un nœud et que chaque nœud contienne une seule particule (au maximum).
Les noeuds sont caractérisés par :
 - leur dimension
 - leur éventuelle particule
 - leurs éventuels noeuds "fils"
 
Après avoir "construit" l'arbre, le centre de masse de chaque nœud est calculé.  
Enfin, les forces exercées sur chaque particules sont calculées, mais cette fois si au lieu de considérer les interactions entre chaque particule, on procède ainsi pour chaque particule :
 On parcourt la liste des noeuds. Dès qu'un noeud de taille s est suffisamment éloigné de la particule (d'une distance d) on calcule l'interaction entre le centre de masse du nœud et la particule considérée. (On vérifie que le noeud est suffisamment éloigné en s'assurant que le rapport s/d est inférieur à un certain "thresold" défini arbitrairement.)
 
Vous pouvez avoir accès à des infos plus détaillées sur cet algo ici : http://www.cs.berkeley.edu/~demmel [...] ure26.html
 
Ma question
 
J'ai déjà codé quelque chose, mais ça ne m'a pas l'air de fonctionner tout à fait correctement pour l'instant :D  
Mais ma vraie question est : quelles précautions dois je prendre pour coder cela en C++ ? Comment le coder de façon optimale en terme de performances ?
Quel algorithme vous parait le plus efficace, et notamment pour la création de l'arbre ? Pour cette étape, qui me parait être assez couteuse, je me suis inspiré d'une méthode trouvée dans ce code : http://pastebin.com/4vL1AKSX (javascript)
 
Pour calculer les centre de masse, j'utilise cette fonction (récursive) :
 

Code :
  1. // compute_cm(root); browses the whole tree
  2. void compute_cm(QuadTree *q)
  3. {
  4.    if(!q) return;
  5.    double total_mass = 0;
  6.    vec cm = vec(0, 0, 0);
  7.    if(!q->leaf) loop(i, 4) if(q->children[i])
  8.    {
  9.        QuadTree *c = q->children[i];
  10.        // the subtree contains a particle
  11.        if(c->leaf && c->b)
  12.        {
  13.            total_mass += c->b->mass;
  14.            cm += c->b->pos * c->b->mass;
  15.        }
  16.        // the substree contains children
  17.        else if(!c->leaf)
  18.        {
  19.            compute_cm(c);
  20.            total_mass += c->mass;
  21.            cm += c->cm * c->mass;
  22.        }
  23.    }
  24.    else if(q->b)
  25.    {
  26.        total_mass = q->b->mass;
  27.        cm = q->b->pos;
  28.    }
  29.    q->cm = cm;
  30.    if(total_mass != 0) { q->cm.div(total_mass); }
  31.    q->mass = total_mass;
  32. }


 
Et enfin pour calculer les accélérations :
 

Code :
  1. // recursive function to calculate the force on each particle
  2. void apply_force(Body *b, QuadTree *q)
  3. {
  4.    if(q->b)
  5.    {
  6.        if(q->b != b) // can occur ?!
  7.        {
  8.            vec dir = (q->b->pos-b->pos);
  9.            dir.normalize();
  10.            if(!dir.iszero()) b->a += dir*G*q->b->mass/(b->pos-q->b->pos).squaredlen();
  11.        }
  12.    }
  13.    else
  14.    {
  15.        vec dir = (q->cm-b->pos);
  16.        double dist_square = dir.squaredlen();
  17.        dir.normalize();
  18.        if(q->size*q->size/dist_square < thresold) { b->a += (q->cm-b->pos).normalize()*G*q->mass/dist_square; }
  19.        else loop(i, 4) if(q->children[i] != NULL) { apply_force(b, q->children[i]); }
  20.    }
  21. }
  22.  
  23. void apply_force(double period)
  24. {
  25.    loopvector(i, bodies)
  26.    {
  27.        bodies[i]->a = vec(0, 0, 0);
  28.        apply_force(bodies[i], root);
  29.        bodies[i]->v += bodies[i]->a * period;
  30.    }
  31.  
  32.    loopvector(i, bodies)
  33.    {
  34.        bodies[i]->pos += bodies[i]->v * period;
  35.    }
  36. }


 
 
Avez vous des idées pour optimiser ce code ? Passer par des fonctions récursives me parait inévitable, mais il me semble aussi que c'est assez couteux en terme de performances et qu'il y a surement des moyens d'améliorer tout ça.
 
Merci d'avance !
 
Au revoir.

mood
Publicité
Posté le 13-07-2012 à 22:47:41  profilanswer
 

n°2149408
Profil sup​primé
Posté le 15-07-2012 à 20:15:42  answer
 

Rebonjour,
 
Pas d'idées ? Si j'ai pas été très clair, s'il manque des informations, n'hésitez pas à me le signaler !
 
Merci

n°2149824
bjone
Insert booze to continue
Posté le 19-07-2012 à 00:57:18  profilanswer
 

A première vue, à part mettre de l'OpenMP dans apply_force() (le quadtree est en "lecture seule" ) comme ça, je vois pas trop.
Je me pose la question si l'actualisation/reconstruction du quadtree est facilement parallélisable efficacement....


Message édité par bjone le 19-07-2012 à 01:03:08
n°2150143
in_your_ph​ion
Posté le 20-07-2012 à 19:51:48  profilanswer
 


 
tu devrais peut être regarder la complexité du quad tree (n log n) et des autres algos que tu pourrais utiliser, c'est ce qui fera la plus grande différence.
 
Une fois choisis, pour le reste c'est de l'optimisation C++, donc ça dépend de ta manière de coder (éviter les copies inutiles, etc) et autrement tu peux utiliser un profiler pour voir où ton algo pêche.

n°2150423
bjone
Insert booze to continue
Posté le 23-07-2012 à 15:39:19  profilanswer
 

Et voir comment améliorer la localité au niveau allocation pour utiliser au mieux le cache.  
Allouer les childrens par pavé de quatre, stocker que le pointeur sur le premier.
Des trucs comme ça.

n°2163230
Profil sup​primé
Posté le 08-11-2012 à 19:10:03  answer
 

Re bonjour et merci beaucoup pour vos réponses !
 
Désolé de ne pas avoir répondu plus tôt. Faute de temps jai du laissé ce projet de côté un petit moment.
 
j'ai changé l'algo de zéro et je suis passé à la 3D (octrees).
J'en suis au début mais les perfs ont l'air "correctes" pour l'instant même si ce n'est pas fini.
 
J'ai cependant un problème : je dois indexer toutes les entités placées dans chaque "octant". Ca c'est bon (je place les pointeurs vers ces entités dans un vector stl par octant).
 
Quand je calcule la force subie par chaque entité, je dois bien sur l'influence que cette entité à sur elle même. Je dois donc traiter de façon particulière les octants contenant cette entité.
Il me suffit en théorie de recherche si le pointeur vers l'entité est indexé dans chaque octant considéré. Mais je pense que si cette recherche est mal codée, l'impact sur les perfs sera très négatif. Donc, j'aurais aimé savoir quelle méthode de recherche était la plus appropriée ? D'après cette discussion : http://stackoverflow.com/questions [...] -stdvector
il y a un risque que "find" soit trop lent si la taille des tableaux est trop grande. Pour l'octant principal le nombre d'élements est ~ 1000 (de 500 à 10 000 en gros) mais ça se réduit vite pour les éléments fils.
 
Puis je utiliser find ? existe-t-il des outils plus performants ? Dois je changer de méthode de stockage (map plutôt que vector ?)
 
Ou dois je avoir carrément une autre approche ? Par exemple, stocker pour chaque entité l'octant le plus petit qui le contient, puis remonter la chaine...
 
Encore merci et désolé d'avoir mis temps de temps à répondre.
 
PS : en ce qui concerne vos remarques, j'ai essayé de les prendre en compte. Pour la reconstruction de l'octree, je me demande si on peut pas le mettre à jour sans reconstruction complète sur des cours intervalles de temps, puis le reconstruire from scratch de temps en temps histoire d'être bien synchro. Ca peut être efficace si le calcul des centres de masse est assez performant.
 
A+

n°2163243
Profil sup​primé
Posté le 08-11-2012 à 21:00:41  answer
 

Re. Après réflexion, mieux vaut indexer, pour chaque entité, les octants qui la contiennent. Si je fais ça la recherche elle même ne diminue pas les perfs de façon significative. Par contre, je sèche pour la suite. Quand je change le code pour que l'influence d'un corps sur lui même ne soit pas calculé, less perfs tombent carrément. Et pas à cause de la recherche donc. Voici le code de cette partie :
 

Code :
  1. void apply_gravity(Entity *e, OctreeNode *o, double step, double thresold)
  2. {
  3.    double ratio = (o->sup-o->inf).squaredlen()/(e->p-(o->barycenter/o->mass)).squaredlen();
  4.  
  5.    bool in_octreenode = in_vector(o, e->nodes); // est ce que l'octant contient l'entité courante ?
  6.    if(o->leaf) // si l'octant n'a pas de noeud fils
  7.    {
  8.        if(in_octreenode) // si l'entité courante est à l'intérieur, on calcule l'influence des corps séparés sauf elle m^me
  9.        {
  10.            for(int i = 0; i < o->ents.size(); ++i) if(o->ents[i] != e)
  11.            {
  12.                e->v += (((o->ents[i]->p)-e->p).normalize() * step * CONST_GRAVITY / (e->p-(o->ents[i]->p)).squaredlen())*o->ents[i]->mass;
  13.            }
  14.        }
  15.        else // sinon, on se sert du barycentre (+rapide)
  16.        {
  17.            e->v += (((o->barycenter/o->mass)-e->p).normalize() * step * CONST_GRAVITY / (e->p-(o->barycenter/o->mass)).squaredlen())*o->mass;
  18.        }
  19.    }
  20.    else if(ratio < thresold && !in_octreenode) // si le corps est suffisamment éloigné de l'octant considéré, et qu'il n'y appartient pas, calcul avec barycentre...
  21.    {
  22.        e->v += (((o->barycenter/o->mass)-e->p).normalize() * step * CONST_GRAVITY / (e->p-(o->barycenter/o->mass)).squaredlen())*o->mass;
  23.    }
  24.    else // sinon, on continue avec les octants fils...
  25.    {
  26.        for(int i = 0; i < 8; ++i)
  27.        {
  28.            if(o->children[i]->entities_count > 0) apply_gravity(e, o->children[i], step, thresold);
  29.        }
  30.    }
  31. }
  32.  
  33. void move_entities(double step, double thresold, OctreeNode *o, std::vector<Entity *> &entities)
  34. {
  35.    int size = entities.size();
  36.    for(int i = 0; i < size; ++i)
  37.    {
  38.        if(o->entities_count > 0) apply_gravity(entities[i], o, step, thresold);
  39.    }
  40.  
  41.    for(int i = 0; i < size; ++i)
  42.    {
  43.        entities[i]->p += entities[i]->v * step;
  44.    }
  45. }

n°2163345
Profil sup​primé
Posté le 09-11-2012 à 23:43:59  answer
 

Bon j'ai pris une vidéo histoire de montrer le fonctionnement. :
http://www.youtube.com/watch?v=MmT [...] e=youtu.be
 
Mais j'arrive pas à savoir si les perfs obtenues sont normales ou non. Je vois des simul avec 10k voire 100k corps sur le net alors que là j'en ai que 1k... Peut être que je dois coupler mon algo avec d'autres algos. Je cherche pas la précision mais j'aimerais pas que ça soit aberrant non plus !
 
Merci

n°2164423
Profil sup​primé
Posté le 17-11-2012 à 17:14:35  answer
 

Je m'autorise un honteux up :o

n°2164432
bjone
Insert booze to continue
Posté le 17-11-2012 à 19:06:24  profilanswer
 

Honnêtement sans avoir tester le truc par moi même, je peux pas trop te dire.
 
Dans ta vidéo tu indiques "4 corps par octant", ça me parait peu.
 
Tu peux aussi regard les implémentation OpenCl & Cuda, peut-être que tu identifiera des erreurs d'approche dans ton code.

mood
Publicité
Posté le 17-11-2012 à 19:06:24  profilanswer
 

n°2164444
Joel F
Real men use unique_ptr
Posté le 17-11-2012 à 20:47:27  profilanswer
 

ca sert a rien d'avoir des structures de données hyepr compelxes si tu passe t a vie a peter ton cache et à faire des opérations non vectorizables.  
 
Avec du SSE et de 'lopenMP tu tape facilement du 50K particules, un peu de tunign fin monte a 100K. Du GPGPU pour ca c'ets de l'overkill.


---------------
[NumScale] [ MKM | EBay ]
n°2164469
Profil sup​primé
Posté le 18-11-2012 à 11:46:00  answer
 

Bonjour !
Merci beaucoup pour vos réponses.
J'ai donc activé l'SSE et j'essaie de mettre l'openMP en place (mais j'ai quelques manips à faire apparemment avec la version express de visual studio)
 
Je vous tiens au courant
A+

n°2164480
Profil sup​primé
Posté le 18-11-2012 à 12:32:37  answer
 

Bon bah j'ai mis un #pragma omp parallel for devant la boucle ligne 36 du code posté ici. (seule boucle sur laquelle c'est "rentable" je crois). Le gain est significatif mais elle prend quand même  un temps important à s'exécuter et du coup ça reste très lent. Pour moi cest vraiment le code ci dessus qui pose problème. La construction de l'octree est assez performante.
J'essaierai des bidouiller quelques param (genre max de corps par octant) mais bon.

 

Merci !


Message édité par Profil supprimé le 18-11-2012 à 12:33:00
n°2164507
Joel F
Real men use unique_ptr
Posté le 18-11-2012 à 16:09:50  profilanswer
 

SSE na se s'active pas, ca se code. Les autovectoriseurs ne font pas de la magie.
et si ton cache est toujours pourri par tes acces memoire, openMP ne changera rien.

 

y a plus de cours d'architecture des ordinateurs au primaire :o

Message cité 1 fois
Message édité par Joel F le 18-11-2012 à 16:10:53

---------------
[NumScale] [ MKM | EBay ]
n°2164510
Profil sup​primé
Posté le 18-11-2012 à 16:42:31  answer
 

Joel F a écrit :

SSE na se s'active pas, ca se code.

oui d'accord mais laissez moi le temps de découvrir vos pistes :D

Joel F a écrit :

Les autovectoriseurs ne font pas de la magie.
et si ton cache est toujours pourri par tes acces memoire, openMP ne changera rien.


y a un net gain de perfs quand même (x2) sur la boucle dont j'ai parlé.
Maintenant pour la construction de l'octree je vois pas trop comment la rendre parallélisable.
Justement, la question c'est aussi comment adapter l'algo pour que je puisse utiliser ces solutions.
 

Joel F a écrit :

y a plus de cours d'architecture des ordinateurs au primaire :o

:sleep:  
 

n°2164538
Profil sup​primé
Posté le 18-11-2012 à 21:55:52  answer
 

bon j'essaie de modifier ma classe vec pour utiliser SSE mais j'ai encore quelques prob parce que c'est un peu chiant de gérer des triplets de double (alors qu'avec des triplets de float c'est beaucoup plus simple) enfin je me comprends.
Je vous tiens au jus

 

Merci encore


Message édité par Profil supprimé le 18-11-2012 à 21:56:07
n°2164661
Profil sup​primé
Posté le 19-11-2012 à 21:54:58  answer
 

Jbricole comm eje peux, par exemple :

 
Code :
  1. vec vec::sse_mul(const vec &a, const double &b) const
  2. {
  3.        vec ret;
  4.        // obliger de faire en 2 fois parce que _mm256_store_pd plante systématiquement.
  5.        __m128d B = _mm_set1_pd(b);
  6.        __m128d R1 = _mm_mul_pd(_mm_load_pd(a.v), B);
  7.        __m128d R2 = _mm_mul_pd(_mm_load_pd(a.v+2), B);
  8.        _mm_store_pd(ret.v, R1);
  9.        _mm_store_pd(ret.v+2, R2);
  10.        return ret;
  11. }
 

avec

Code :
  1. vec  operator* (double f) const { return sse_mul(*this, f); }
 

Ca a l'air de marcher correctement quand je l'appelle de façon normale au début du code pour tester mais ça plante lamentablement par la suite.
Je vois pas trop d'où ça peut venir.


Message édité par Profil supprimé le 19-11-2012 à 21:55:30
n°2164664
bjone
Insert booze to continue
Posté le 19-11-2012 à 23:03:05  profilanswer
 

Problème d'alignement ;)
C'est le truc chiant avec le SSE, au lieu de ralentir ça sblouche.


Message édité par bjone le 19-11-2012 à 23:03:18
n°2164665
bjone
Insert booze to continue
Posté le 19-11-2012 à 23:04:52  profilanswer
 

Mais de toutes façons, ce sont des structures ou la manière de les manipuler qui posent problème.

n°2164813
Joel F
Real men use unique_ptr
Posté le 21-11-2012 à 07:08:57  profilanswer
 

faire du Structure of Array plutot que du Array of Structure deja
ensuite _m256 c'ets du AVX, donc alignement sur 32 octets de la memoire obligatoire.


---------------
[NumScale] [ MKM | EBay ]
n°2164958
Profil sup​primé
Posté le 21-11-2012 à 17:44:08  answer
 

'soir,
 
Merci pour vos réponses, je corrigerai ça dès que je peux ! faut que je me documente un peu.
 
:jap:

n°2165259
Profil sup​primé
Posté le 23-11-2012 à 12:55:16  answer
 

Bon, bah rien que ça je m'en sors pas... Pourtant je manipule des objets simples :
 

Code :
  1. class vec
  2. {
  3. public:
  4.    union
  5.    {
  6.        struct { double x, y, z, w; };
  7.        double v[4];
  8.    };
  9. // [...]


 
Mais impossible de m'en sortir avec ça, même en jouant avec __declspec(align()). Auriez vous des idées et surtout des liens à consulter traitant de ce problème ? J'ai fouillé sur le net mais rien trouvé de bien éclairant.
 
@JoelF : qu'entends tu par des "stucture of array" dans mon cas ? Tu veux dire me débarrasser carrément de cette classe vec ?  
 
Merci.

mood
Publicité
Posté le   profilanswer
 


Aller à :
Ajouter une réponse
  FORUM HardWare.fr
  Programmation
  C++

  interaction à grand nombre de corps et quadtree (barnes hut)

 

Sujets relatifs
Intégrer un tchat à un site web avec intéraction sur iphoneComment compter le nombre de lignes dans un tableau croisé dynamique ?
interaction entre une application web mobile et un serveur LAMPnombre d'occurrences dans un XML avec PHP
compter le nombre de phrases dans un paragrapheTableau - revenir au début si déplacement trop grand
compter le nombre de champs vides dans 1 enregistrement SQLGrand moche!
Excel macro, grille de nombre, changer couleur fond selon choixCompter le nombre de résultats d'une requête
Plus de sujets relatifs à : interaction à grand nombre de corps et quadtree (barnes hut)



Copyright © 1997-2016 Hardware.fr SARL (Signaler un contenu illicite) / Groupe LDLC / Shop HFR