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 )
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  
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) :
 

// compute_cm(root); browses the whole tree
void compute_cm(QuadTree *q)
{
    if(!q) return;
    double total_mass = 0;
    vec cm = vec(0, 0, 0);
    if(!q->leaf) loop(i, 4) if(q->children[i])
    {
        QuadTree *c = q->children[i];
        // the subtree contains a particle
        if(c->leaf && c->b)
        {
            total_mass += c->b->mass;
            cm += c->b->pos * c->b->mass;
        }
        // the substree contains children
        else if(!c->leaf)
        {
            compute_cm(c);
            total_mass += c->mass;
            cm += c->cm * c->mass;
        }
    }
    else if(q->b)
    {
        total_mass = q->b->mass;
        cm = q->b->pos;
    }
    q->cm = cm;
    if(total_mass != 0) { q->cm.div(total_mass); }
    q->mass = total_mass;
}

 
Et enfin pour calculer les accélérations :
 

// recursive function to calculate the force on each particle
void apply_force(Body *b, QuadTree *q)
{
    if(q->b)
    {
        if(q->b != b) // can occur ?!
        {
            vec dir = (q->b->pos-b->pos);
            dir.normalize();
            if(!dir.iszero()) b->a += dir*G*q->b->mass/(b->pos-q->b->pos).squaredlen();
        }
    }
    else
    {
        vec dir = (q->cm-b->pos);
        double dist_square = dir.squaredlen();
        dir.normalize();
        if(q->size*q->size/dist_square < thresold) { b->a += (q->cm-b->pos).normalize()*G*q->mass/dist_square; }
        else loop(i, 4) if(q->children[i] != NULL) { apply_force(b, q->children[i]); }
    }
}
 
void apply_force(double period)
{
    loopvector(i, bodies)
    {
        bodies[i]->a = vec(0, 0, 0);
        apply_force(bodies[i], root);
        bodies[i]->v += bodies[i]->a * period;
    }
 
    loopvector(i, bodies)
    {
        bodies[i]->pos += bodies[i]->v * period;
    }
}

 
 
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.

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

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....

 
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.

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.

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+

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 :
 

void apply_gravity(Entity *e, OctreeNode *o, double step, double thresold)
{
    double ratio = (o->sup-o->inf).squaredlen()/(e->p-(o->barycenter/o->mass)).squaredlen();
 
    bool in_octreenode = in_vector(o, e->nodes); // est ce que l'octant contient l'entité courante ?
    if(o->leaf) // si l'octant n'a pas de noeud fils
    {
        if(in_octreenode) // si l'entité courante est à l'intérieur, on calcule l'influence des corps séparés sauf elle m^me
        {
            for(int i = 0; i < o->ents.size(); ++i) if(o->ents[i] != e)
            {
                e->v += (((o->ents[i]->p)-e->p).normalize() * step * CONST_GRAVITY / (e->p-(o->ents[i]->p)).squaredlen())*o->ents[i]->mass;
            }
        }
        else // sinon, on se sert du barycentre (+rapide)
        {
            e->v += (((o->barycenter/o->mass)-e->p).normalize() * step * CONST_GRAVITY / (e->p-(o->barycenter/o->mass)).squaredlen())*o->mass;
        }
    }
    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...
    {
        e->v += (((o->barycenter/o->mass)-e->p).normalize() * step * CONST_GRAVITY / (e->p-(o->barycenter/o->mass)).squaredlen())*o->mass;
    }
    else // sinon, on continue avec les octants fils...
    {
        for(int i = 0; i < 8; ++i)
        {
            if(o->children[i]->entities_count > 0) apply_gravity(e, o->children[i], step, thresold);
        }
    }
}
 
void move_entities(double step, double thresold, OctreeNode *o, std::vector<Entity *> &entities)
{
    int size = entities.size();
    for(int i = 0; i < size; ++i)
    {
        if(o->entities_count > 0) apply_gravity(entities[i], o, step, thresold);
    }
 
    for(int i = 0; i < size; ++i)
    {
        entities[i]->p += entities[i]->v * step;
    }
}

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

Je m'autorise un honteux up

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.

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.

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+

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 !

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

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

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.
 

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

Jbricole comm eje peux, par exemple :

vec vec::sse_mul(const vec &a, const double &b) const
{
        vec ret;
        // obliger de faire en 2 fois parce que _mm256_store_pd plante systématiquement.
        __m128d B = _mm_set1_pd(b);
        __m128d R1 = _mm_mul_pd(_mm_load_pd(a.v), B);
        __m128d R2 = _mm_mul_pd(_mm_load_pd(a.v+2), B);
        _mm_store_pd(ret.v, R1);
        _mm_store_pd(ret.v+2, R2);
        return ret;
}

avec

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.

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

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

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.

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

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

class vec
{
public:
    union
    {
        struct { double x, y, z, w; };
        double v[4];
    };
// [...]

 
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.