Optimisation code de diffusion d'atomes

[produ13]

Bonsoir,

je suis étudiant et je dois faire un programme sur la diffusion. voici ce qu'on a codé pour l'instant (le programme marche bien). Pourriez vous s'il vous plait m'aider à l'améliorer en vitesse ou autre, je suis ouvert à toutes vos propositions. je vous remercie.

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
clear all;
n=input('entrer le nombre datomes ' );
t=input('entrer le nombre diterations ');
tic;
S1=zeros(n);
S1(round(n/2),:)=1;
L=ceil(4*rand(t,n));    
S2=S1;
A=zeros(1,n);
B=zeros(1,n);
x=zeros(1,n);
y=zeros(1,n);
a=(1:n);
for p=1:t
    S1=S2;
    [Y,X]=find(S1);
    S2=zeros(n);
    for m=1:n
        A(m)=X(m); 
        B(m)=Y(m); 
        if (L(p,m)==1 && Y(m)+1==n+1)
            B(m)=0;
        elseif (L(p,m)==2 && Y(m)-1==0)
            B(m)=n+1;
        end

        if (L(p,m)==3 && X(m)+1==n+1)
            A(m)=0;
        elseif (L(p,m)==4 && X(m)-1==0)
            A(m)=n+1;
        end

        if (L(p,m)==1 && S1(B(m)+1,X(m))==0)
            x(m)=X(m);
            y(m)=B(m)+1;
        elseif (L(p,m)==1 && S1(B(m)+1,X(m))==1)
            x(m)=X(m);
            y(m)=Y(m);
        end
       
        if (L(p,m)==2 && S1(B(m)-1,X(m))==0)
            x(m)=X(m);
            y(m)=B(m)-1;
        elseif (L(p,m)==2 && S1(B(m)-1,X(m))==1)
            x(m)=X(m);
            y(m)=Y(m);
        end
    
        if (L(p,m)==3 && S1(Y(m),A(m)+1)==0)
            x(m)=A(m)+1;
            y(m)=Y(m);
        elseif (L(p,m)==3 && S1(Y(m),A(m)+1)==1)
            x(m)=X(m);
            y(m)=Y(m);
        end

        if (L(p,m)==4 && S1(Y(m),A(m)-1)==0)
            x(m)=A(m)-1;
            y(m)=Y(m);
        elseif (L(p,m)==4 && S1(Y(m),A(m)-1)==1)
            x(m)=X(m);
            y(m)=Y(m);
        end
        S2(y(m),x(m))=S2(y(m),x(m))+1;
        if (S2(y(m),x(m))==2)
            S2(y(m),x(m))=1;
            x(m)=X(m);
            y(m)=Y(m);
            S2(y(m),x(m))=1;
        end
    end
end
figure(2);
imagesc(S2);
X2=zeros(n,1);
for u=1:n
X2(u,1)=sum(S2(u,:));
end
figure(3)
plot(a, X2, '*')
toc

Pourriez vous transfert le message dans le forum Environnements de développement> MATLAB svp

[ol9245]

J'ai essayé de regarder ton code pour voir si je pouvais te donner un coup de main. En fait, je n'ai trouvé qu'un seul conseil pour l'instant :

ton code fait 80 lignes.
avant de le faire lire à quelqu'un d'autre, écris 160 lignes de commentaire en plus pour expliquer ce que tu fais.

2 lignes de commentaire pour une ligne de code, c'est à peu près le minimum que tu peux trouver dans un code normal.

[produ13]

Merci pour votre conseil.
Dans ce message l'explication du programme ainsi que le code détaillé étape par étape.

Tout d'abord je suis étudiant en école d'ingénieur Le but de ce TD est de nous apprendre les rudiments de MATLAB tout en étudiant un phénomène physique particulier, ici la diffusion des atomes (migration des atomes sous l'effet de l'agitation thermique)

Voila le principe:
On crée une matrice de zéros qui sera l'espace de référence dans lequel les atomes diffusent. Pour modéliser les atomes, on ajoute des 1 dans la ligne centrale de la matrice. Tous les atomes se déplacent l'un après l'autre, une fois, soit en haut, à droite, bas ou gauche dans la matrice (déplacement aléatoire). Puis ça recommence...

Conditions:
Il ne peut pas y avoir deux atomes a une même position dans la matrice.
Les atomes ne peuvent pas sortir de la matrice.

Pour finir le programme édite un profil de l'espace avec la position finale des atomes et un graphe de répartition des atomes selon leur positions avec un lissage de courbe pour obtenir une belle gaussienne.

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
clear all;
n=input('entrer le nombre datomes ' );
t=input('entrer le nombre de déplacements ');
tic;
S1=zeros(n);
%génère l'espace de référence qui va accueillir les atomes
S1(round(n/2),:)=1;
%ligne d'atomes qui vont diffuser
L=ceil(4*rand(t,n));
%matrice de déplacements, les déplacement dont 1, 2, 3, 4. si le 1 alors l'atome se déplace droite, haut, , bas, , gauche
S2=S1;
A=zeros(1,n);
B=zeros(1,n);
%A et B sont les valeurs tampon de boucle
x=zeros(1,n);
%variable de positions verticale des atomes après leur déplacement
y=zeros(1,n);
%variable de positions horizontale des atomes après leur déplacement
a=(1:n);

for p=1:t
    %pour chaque déplacement
    S1=S2;
    [Y,X]=find(S1);
    %X et Y variables de positions des atomes avant leur déplacement

    %permet de localiser les atomes dans la matrice
    S2=zeros(n);

    for m=1:n
        %pour chaque atome
        A(m)=X(m);
        B(m)=Y(m);
        %A et B sont les valeurs tampon de la boucle

        if (L(p,m)==1 && Y(m)+1==n+1)
            B(m)=0;
            %si l'atome sort de la matrice par le bas, on change la valeur tampon de sa position pour lui permettre de rester dans la matrice après son déplacement. l'atome est donc placé au dessus de la matrice
        elseif (L(p,m)==2 && Y(m)-1==0)
            B(m)=n+1;
        end
        %si l'atome sort de la matrice sort par le haut,
        if (L(p,m)==3 && X(m)+1==n+1)
            A(m)=0;
            %si l'atome sort de la matrice sort par la droite,
        elseif (L(p,m)==4 && X(m)-1==0)
            A(m)=n+1;
        end
        %si l'atome sort de la matrice sort par la gauche,

        if (L(p,m)==1 && S1(B(m)+1,X(m))==0)
            x(m)=X(m);
            y(m)=B(m)+1;
            %déplacement de l'atome vers le bas si il n'y a pas d'atome en dessous
        elseif (L(p,m)==1 && S1(B(m)+1,X(m))==1)
            x(m)=X(m);
            y(m)=Y(m)
            %l'atome reste à sa position si il y a déplacement vers le bas et un atome en dessous
        end
        if (L(p,m)==2 && S1(B(m)-1,X(m))==0)
            x(m)=X(m);
            y(m)=B(m)-1;
        elseif (L(p,m)==2 && S1(B(m)-1,X(m))==1)
            x(m)=X(m);
            y(m)=Y(m);
        end
        if (L(p,m)==3 && S1(Y(m),A(m)+1)==0)
            x(m)=A(m)+1;
            y(m)=Y(m);
        elseif (L(p,m)==3 && S1(Y(m),A(m)+1)==1)
            x(m)=X(m);
            y(m)=Y(m);
        end
        if (L(p,m)==4 && S1(Y(m),A(m)-1)==0)
            x(m)=A(m)-1;
            y(m)=Y(m);
        elseif (L(p,m)==4 && S1(Y(m),A(m)-1)==1)
            x(m)=X(m);
            y(m)=Y(m);
        end

        S2(y(m),x(m))=S2(y(m),x(m))+1;
        if (S2(y(m),x(m))==2)
            S2(y(m),x(m))=1;
            x(m)=X(m);
            y(m)=Y(m);
        end

    end
end



figure(2);
imagesc(S2);
X2=zeros(n,1);

for u=1:n
    X2(u,1)=sum(S2(u,:));
end
%permet de sommer les atomes pour avoir la répartition dans les colonnes

figure(3)
plot(a, X2, '*')
toc

[ol9245]

Si tu veux optimiser ton code spécifiquement pour matlab, il faut virer toutes le boucles for...next. C'est un peu compliqué à imaginer quand on n'a pas l'habitude. Dans ton cas, l'idée est de séparer ton tableau d'atomes en 4 tableaux selon la direction où ils vont. Dans un des tableaux, tu peux faire tous les déplacement en une seule opération en faisant un shift sur les lignes ou les colonnes, selon la direction définie.

En plus, ton code ne simule pas un mouvement brownien car tu exclu la possibilité qu'un atome reste dans sa case. Soit ta case a une dimension physique très petite par rapport au déplacement moyen d'un atome dans un pas de temps, dans ce cas là tu dois envisager qu'un atome puisse parcourir plusieurs cases. Soit ta case est grande et tu dois accepter que ton atome puisse ne pas en sortir.

Bref, ce que tu as simulé, c'est plus proche du jeu des chaises musicales que de la diffusion. Perso, il m'arrive de bosser sur la diffusion (celle des grains de sable quand ils sont impactés par des gouttes de pluie) et je n'utiliserai jamais un code comme le tiens.

Autre point discutable : Le fait qu'un atome ne puisse pas atteindre les cases diagonales est également très discutable. Pour respeter le sens physique, tu devrais faire alternativement des mouvement aléatoires en X et en Y, et pas l'un ou l'autre dans le même pas de temps.

[ol9245]

Voici un code qui fait à peu près ce que tu as décrit sans jamais boucler sur les lignes et les colonnes de la matrice.

Ne jamais boucler sur les lignes et les colonnes de la matrice, C'est LA grande règle d'optimisation des codes matlab.

Si il y a des trucs que tu ne comprends pas dans mon petit code, explore l'aide de matlab ou repique les morceaux de code que tu comprends pas dans des bouts d'essai pour voir ce qu'ils font.

Si le code te plait, à toi de l'adapter à tes besoins (je le ferai pas pour toi ). Je t'ai écrit cet exemple surtout pour te mettre sur la voie d'un code matlab optimisé car, comme déja dit, c'est pas toujours très intuitif.

[ol9245]

et avec le code c'est mieux !

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
function brownien()
taille = 100 ; % longueur du côté de la matrice
M = rand(taille) ; % une carte de 100 x 100
n = 100 ; % nombre d'atomes
tmp = sort(M(:)) ; % trier les valeurs de M pour ne garder que les 100 plus grandes
seuil = tmp(numel(M)-n+1) ; % coupe à 100 atomes
M(M<seuil) = 0 ; % pas d'atome
M(M>=seuil) = 1 ; % un atome

direction = true ; % de haut en bas. On alterne avec false = de gauche à droite
for k = 1:1000 % pour chaque round
    if direction
        arrivee = deplacementh(M) ;
    else
        arrivee = deplacementv(M) ;
    end
    direction = ~direction ; % bascule la direction
    M = arrivee ; % repart pour le pas suivant
    imshow(M) ; % affiche la situation courante
    drawnow ; % force l'affichage maintenant pour voir l'animation
end

% fonction interne.
% cette fonction peut lire toutes les variables de la fonction mère. C'est pratique pour écrire vite du code bien cochon.
% déplacement des atomes de M horizontalement
    function A = deplacementv (M)
        
        D = distance() ; % tire les distances de déplacement entre -1 et 1
        
        % calcule le résultat des déplacements versle haut
        enhaut = M ;
        enhaut(D~=-1)=0 ; % détruit les atomes qui ne montent pas
        enhaut2 = [enhaut(2:end, :) ; zeros(1, taille)] ; % fait tout monter d'une ligne, complète d'une ligne vide en bas.
        
        % calcule le résultat des déplacements vers le bas
        enbas = M ;
        enbas(D~=1)=0 ; % détruit les atomes qui ne descendent pas
        enbas2 = [zeros(1, taille) ; enbas(1:end-1, :)] ; % fait tout descendre, complète d'une ligne vide en haut
        
        % calcul le résultat des atomes immobiles
        ici = M ;
        ici(D~=0) = 0 ; % détruit les atomes qui bougent
        
        % consolide tous les déplacements.
        % ici, les atomes multiples dans une même case sont purement et simplement fusionnés en un seul :-(
        A = min(1, enhaut2 + enbas2 + ici) ;
    end

    function A = deplacementh (M)
        A = deplacementv(M').' ; % transpose, déplace verticalement et transpose à nouveau
    end

    function D = distance()
        D = round(rand(taille)*3 -1.5) ; % de -1 case à + 1 case
    end

end

[produ13]

ok Merci beaucoup pour le code, je vais regarder ça.



En ce qui concerne mon modèle de déplacement aléatoire des atomes, j'obtiens des résultats qui sont très proches du modèle théorique de la 2ème loi de Fick:

C(x,t)=C0/√πDt*e^(-(x^2/Dt))

Puisque quand je fit mes résultats avec cette loi normale, voici ce que j'obtiens



(avec en ordonnées le nombre d'atomes par ligne dans le tableau et en abscisse le numéro de la ligne du tableau).

D'autre part le fait que les sauts dans les diagonales soit impossible me paraît justifiable dans le cas où on limite les déplacements aux cases contiguës. L'atome a une énergie qui peut le faire déplacer au maximum de la longueur qui le sépare d'une case adjacente.


Mais c'est vrai que il manque qq chose dans le programme...

[ol9245]

Envoyé par produ13

ok Merci beaucoup pour le code, je vais regarder ça.



En ce qui concerne mon modèle de déplacement aléatoire des atomes, j'obtiens des résultats qui sont très proches du modèle théorique de la 2ème loi de Fick:

C(x,t)=C0/√πDt*e^(-(x^2/Dt))

Mais c'est vrai que il manque qq chose dans le programme...

ça veut simplement dire que ton modèle théorique X et ton modèle numérique Y fournissent tous les deux la loi la plus répandue dans le monde (Gauss). Comme preuve que X=Y, on a déja vu mieux

pas besoin de chercher très loin une prévision un peu plus forte de X, que Y ne respecte pas.

X = le mouvement brownien
Y = ton code

D'une manière générale il est difficile de prouver que un code Y respecte un modèle théorique X. Le mieux qu'on puisse faire est de montrer qu'on a passé beaucoup de temps à essayer de montrer que X ~= Y et que on n'y est pas arrivé. Toi, tu trouves LA prédiction théorique que Y est capable de simuler correctement.... pas tout à fait la bonne approche.

[produ13]

Bonsoir,

effectivement, l'erreur réside dans le fait que j'applique une loi macroscopique (Loi de Fick) à un problème qui est microscopique.
En poussant un peu ma machine et après de longues minutes d'attentes, voila le profil que j'obtient avec 3000 déplacements pour chaque atome.



On est assez loin modèle théorique (la courbe) !

Pour améliorer ça il faut tenir compte de l'environnement de chaque atome c'est à dire si l'atome est entouré ou non d'un autre atome. En principe cela ramène le problème à l'échelle microscopique et cela nous rapproche plus de la réalité.

Nous avons aussi décidé d'ajouter une fonction pdepe (résolution d'une nouvelle équation différentielle). En gros ça devrait nous permettre de fixer un paramètre (D, coefficient de diffusion ou t, temps total de la diffusion) et d'utiliser notre programme pour simuler la diffusion d'atome dans des matériaux bien précis.

Désolé si il y a des redite par rapport au code que vous avez envoyé.

[ol9245]

Envoyé par produ13

Bonsoir,

Pour améliorer ça il faut tenir compte de l'environnement de chaque atome c'est à dire si l'atome est entouré ou non d'un autre atome. En principe cela ramène le problème à l'échelle microscopique et cela nous rapproche plus de la réalité.

Nous avons aussi décidé d'ajouter une fonction pdepe (résolution d'une nouvelle équation différentielle). En gros ça devrait nous permettre de fixer un paramètre (D, coefficient de diffusion ou t, temps total de la diffusion) et d'utiliser notre programme pour simuler la diffusion d'atome dans des matériaux bien précis.

Désolé si il y a des redite par rapport au code que vous avez envoyé.

Bonne chance dans vos travaux. (je suppose que tu n'es pas seul car tu dis "nous").
Faites attention à ne pas réinventer le mur mittoyen : il y a des étagères complètes de livres sur la diffusion. Je bosse avec un collègue du CEA qui a passé 15 ans rien que sur cette équation... Il y a donc fort à douter que vous puissiez créer du nouveau en quelques heures. => ces quelques heures vous seront plus profitables à feuilleter les bons livres d'une bibliothèque scienifique digne de ce nom (Eyrolles à Paris) où à piocher Google avec les mots clé qui vont bien. Quand vous aurez trouvé une solution déja publiée et proche de votre problème, codez-la et le tour est joué.

Pour répondre plus près de ta question :
* vision cartésienne : admettre plus que un atome par case et itérer sur chaque atome de chaque case.
* vision lagrangienne : pas de case : tu ne gères que des atomes et leur position.

Note : le mouvement brownien ne tiens pas nécéssairement compte des colisions. Le mouvement de la poussière en suspension dans l'air est brownien, et on peut le décrire en ignorant les colisions. Les chocs entre particules doivent être pris en compte uniquement si tu veux relier mouvement brownien et chaleur.