Salut,
(Désolé pour ceux qui on déjà lu mon intro plsuieurs fois)
J'ai une grille 2D (i,j) = [0..nx ; 0..ny] sur lequel est discrétisé l'équation :
u(x,y) - d^2Laplacien(u(x,y)) = f(x,y)
où u et f sont des vecteurs à trois composantes u_x,u_y, u_z, et f_x, f_y, f_z.
J'impose à u d'avoir toutes ces composantes périodiques selon la direction x, c'est à dire :
u_x(0,y) = u_x(nx,y)
u_y(0,y) = u_y(nx,y)
u_z(0,y) = u_z(nx,y)
Et sur les bords y=constante j'impose respectivement :
Condition de dirichlet aux composantes tangentielles à la frontière (u_x et u_z donc) et une condition de neumann à la condition normale à la frontière : u_y
c'est à dire :
u_x(x,0) = B1
u_x(x,ny) = B2
u_z(x,0) = 0
u_z(x,ny) = 0
et
du_y/dy (x,0) = du_y/dy(x,ny) = 0
Pour résoudre cette équation, j'utilise la méthode multigrille. Avec un smoother gauss-seidel dont l'ordering est un simple forward. Je fais 3 pre et post smoothing par cycle V et pour l'instant je n'ai pas de critère de convergence, j'itère sur 50 cycles pour voir ce qui se passe. Mon opérateur de restriction est le classique full weighting et mon je prolonge par interpolation bilinéaire. Je pars d'une grille (nx,ny) et descent jusqu'à là grille 2x2 où je fais 5 itérations de gauss seidel.
Afin de tester mon solver, je prend deux exemple de mon équation :
u_x - d^2Laplacien(u_x) = cos(a1*y)
u_y - d^2Laplacien(u_t) = cos(a2*x)
Je fixe B1=B2 = 0 (conditions de dirichlet homogène pour u_x)
La taille "physique" de mon domaine est xm=ym=1, le pas d'espace est donné par dx = xm/nx et dy = ym/ny.
Je fixe les parametres suivants :
a1 = 0.7*2pi/ym
a2 = 2pi/xm
50 cycles V
npressmooth = nPostSmooth=3
d^2 = 0.02
Voici mes questions :
Mon solver semble bien fonctionner pour nx=ny, quelque soit leur valeur (puissances de 2) j'ai testé de 16x16 à 1024x1024.
En revanche, j'ai parfois des problemes pour certains ratio de nx/ny. Toutes les combinaisons nx/ny=2 ou ny/nx=2 fonctionnent. J'ai systématiquement des ennuis lorsque ny/nx > 2 pour u_x et u_y, et parfois des ennuis lorsque nx/ny > 2 pour u_y.
Voici une liste de tests que j'ai fait tourner :
512² ok
512x256 ok
512x128 ok u_x mais pas u_y
512x64 ok u_x mais pas u_y
512x32 ok u_x mais pas u_y
512x16 ok (etrange !)
256x256 ok
256x128 ok
256x64 ok u_x mais pas u_y
256x32 ok u_x mais pas u_y
256x16 ok (étrange aussi j'aurais parié non pour u_y)
256x8 ok (encore étrange)
128^2 ok
128x64 ok
128x32 ok
128x16 ok
128x8 ok
64x64 ok
64x32 ok
64x16 ok
64x8 ok
32x32 ok
32x16 ok
32x8 ok
256x512 ok
128x512 aucun ne converge
64x512 aucun ne converge
32x512 aucun ne converge
32x512 aucun ne converge
16x512 aucun ne converge
128x256 ok
64x256 aucun ne converge
32x256 aucun ne converge
16x256 aucun ne converge
64x128 ok
32x128 aucun ne converge
16x128 aucun ne converge
32x64 ok
16x64 aucun ne converge
16x32 ok
16x16 ok
Voilà....
Je me demande si j'ai pas zappé quelque chose quant à la convergence de tout ça ? J'ai l'impression que mon code fonctionne mais que les paramètre ou ça diverge sont faux mathématiquement (i.e. si je connaissais les règles je demanderais pas ces parametres)
Je dis ça car, j'ai relancé les cas suivant avec xm=ym = 20 :
128x512 ; 64x512 ; 32x512 qui divergeaient pour xm=ym=1
et ils fonctionnent bien. Ce qui tend à me faire penser que je n'ai pas d'erreur dans la gestion de mes dimensions mais qu'il y a bien un critère de convergence sur la combinaison (dx,dy)
Quelqu'un peut-il m'éclairer ?
Merci
ps : je peux fournir mon code si besoin
Partager