Discussion :

[Hyunkel07]

Bonjour à tous,

Il y quelque jours j´ai voulu commencer à apprendre à coder en C++ en réécrivant avec la bonne syntax un petit programme MATLAB tiré d´un livre de modélisation.

Il s´agit d´un simple calcul de champs de vitesse en 2D, d´un manteau (W = 1000km x H = 1500km) subissant une convection.
Le champs de vitesse est décrit par: Vx = -vx0 * sin(2*pi(x/W))*cos(pi(y/H)) (composante horizontale) et Vy = vy0 * cos(2*pi(x/W))*sin(pi(y/H)) (composante verticale). x et y étant les coordonées horizontales et verticales du modèle.
Le modèle est discritisé de la manière suivante: grille en 2D (36 x 36) regulièrement distribuer.

J´ai réécris le programme en C++ assez facilement (il marche et me donne les mêmes résultats que le programme MATLAB) mais n´étant qu un débutant je me suis demandé si il était vraiment bien écrit et optimiser (Je ne pense pas ...) . Auriez-vous quelques conseils/ idées ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
#include <iostream>
#include <cmath>
#include <fstream>
 
using namespace std;
 
int main()
{
    const long xsize = 1000000,ysize = 1500000; //Horizontal and Vertical size (m)
    const int xnum = 36,ynum = 36; // Horizontal and Vertical resolution
    const double vx0 = 1.0e-9 * xsize / 2 / ysize,vy0 = 1.0e-9,pi = 3.14159; //Scale for Horizontal and Vertical Velocity (m/s)
    const double xstp = xsize/(xnum-1),ystp = ysize/(ynum-1); // Horizontal and Vertical grid step
    double vx[xnum][ynum],vy[xnum][ynum],dvxdx[ynum][xnum],dvydy[ynum][xnum],divv[ynum][xnum]; // Declaration of the Horizontal and Vertical velocity component and Partial DerivativesdVx/dx and dVy/dy
    double x[xnum]={0.0}, y[ynum]={0.0}; // Declaration of x and y variables (Coordinates m)
    int n=0; //counter for the first "for" loop
 
    //Generate x and y coordinates of the nodal point (Grid 36x36)
    for ( n = 0; n < xnum; n++)
     {
       x[n] = n * xstp;
       y[n] = n * ystp;
     }
 
    /*for ( m = 0; m < ynum; m++)
     {
       y[m] = y[m-1] + ystp;
       cout << "y " << m << " equal to : " << y[m] <<"\n" << endl;
     }*/
 
    // Create and test files to stock results
    ofstream velox;
    velox.open("./vx.xy", ios::out); //on ouvrre le fichier en ecriture
    if (velox.bad()) //permet de tester si le fichier s'est ouvert sans probleme
    return 1;
    ofstream veloy;
    veloy.open("./vy.xy", ios::out); //on ouvrre le fichier en ecriture
    if (veloy.bad()) //permet de tester si le fichier s'est ouvert sans probleme
    return 1;
    ofstream velodvdy;
    velodvdy.open("./dvdy.xy", ios::out); //on ouvrre le fichier en ecriture
    if (velodvdy.bad()) //permet de tester si le fichier s'est ouvert sans probleme
    return 1;
    ofstream velodvdx;
    velodvdx.open("./dvdx.xy", ios::out); //on ouvrre le fichier en ecriture
    if (velodvdx.bad()) //permet de tester si le fichier s'est ouvert sans probleme
    return 1;
    ofstream velodivv;
    velodivv.open("./divv.xy", ios::out); //on ouvrre le fichier en ecriture
    if (velodivv.bad()) //permet de tester si le fichier s'est ouvert sans probleme
    return 1;
 
 
    int i(0),j(0); //counters for the two followings loops
 
    for (i = 0; i < ynum; i++)
       {
       for(j = 0; j < xnum; j++)
        {
 
        //Computing Velocity Fields Vx and Vy of the convecting mantle (Horizontal and Vertical component)
 
        vx[i][j] = - vx0 * sin(2*pi*(x[j]/xsize)) * cos(pi*(y[i]/ysize));
 
        vy[i][j]= vy0 * cos(2*pi*(x[j]/xsize)) * sin(pi*(y[i]/ysize));
 
        // Computing partial derivatives
 
        dvxdx[i][j]=-vx0*pi/xsize*2*cos(pi*x[j]/xsize*2)*cos(pi*y[i]/ysize);
 
        dvydy[i][j]=vy0*pi/ysize*cos(pi*y[i]/ysize)*cos(pi*x[j]/xsize*2);
 
        //Computing div(v)=dVx/dx+Dvy/dy
 
        divv[i][j]=dvxdx[i][j]+dvydy[i][j];
 
        //Write results in files
        velox << vx[i][j] << "\t";
        veloy << vy[i][j] << "\t";
        velodvdx << dvxdx[i][j] << "\t";
        velodvdy << dvydy[i][j] << "\t";
        velodivv << divv[i][j] << "\t";
 
        }
        velox << vx[i][j] << endl;
        veloy << vy[i][j] << endl;
        velodvdx << dvxdx[i][j] << endl;
        velodvdy << dvydy[i][j] << endl;
        velodivv << divv[i][j] << endl;
       }
        // Close files
        velox.close(); //on ferme le fichier pour liberer la mémoire
        veloy.close(); //on ferme le fichier pour liberer la mémoire
        velodvdx.close(); //on ferme le fichier pour liberer la mémoire
        velodvdy.close(); //on ferme le fichier pour liberer la mémoire
        velodivv.close(); //on ferme le fichier pour liberer la mémoire
    return 0;
}

Merci

Rémi

[Iradrille]

Hello,

C'est correct, juste 2 détails :
De manière générale on réduit la portée des variables au maximum et donc on limite la portée des compteurs de boucles dans les boucles

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
for(int i=0; i<42; ++i) {
   // ...
}
// Au lieu de 
int i=0;
for(; i<42; ++i) {
   // ...
}

Pour l'ouverture de tes fichiers : un std::ofstream est automatiquement ouvert en écriture, pas besoin de la préciser, tu peux ouvrir aussi les fichiers lors de la création de l'objet, et std::ofstream::operator bool() est préférable pour tester la validité du flux.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
std::ofstream velox("./vx.xy");
if(!velox) {
   return 1;
}
// Au lieu de
ofstream velox;
velox.open("./vx.xy", ios::out); //on ouvrre le fichier en ecriture
if (velox.bad()) { //permet de tester si le fichier s'est ouvert sans probleme
   return 1;
}

Le fichier est aussi automatiquement fermé quand l'objet std::ofstream est détruit, tes velox.close(); sont donc inutiles à la fin de ton main (mais c'est pas faux).

Niveau optimisation, il y à bien sur quelques améliorations possibles, parmi les plus évidentes / simples à mettre en places :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
vx[i][j] = - vx0 * sin(2*pi*(x[j]/xsize)) * cos(pi*(y[i]/ysize)); // ligne 62
vy[i][j]= vy0 * cos(2*pi*(x[j]/xsize)) * sin(pi*(y[i]/ysize)); // ligne 64
 
// sin(2*pi*(x[j]/xsize)) et cos(pi*(y[i]/ysize)) peuvent être calculés 1 fois avant ta double boucle
// pas besoin de les calculer plusieurs fois chacun à chaque itération

Çà te permet de virer ces sin / cos de la boucle principale, ce qui est probablement ce qui prend le plus de temps dans les calculs.

Il est probablement rentable de faire les traitements en 2 étapes : calculs puis écriture des résultats sur le disque pour ne pas pourrir le temps de calcul avec des I/O.

Le code à l'air facilement parallélisable (vectorisation et/ou multithreading) car il n'y à aucune dépendances entre les itérations, mais vu le peu de calculs face à la quantité d'I/O; faut vraiment avoir besoin de ce probable faible gain de performance pour se prendre la tête avec ça.