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

  FORUM HardWare.fr
  Programmation
  Python

  problème avec de manager.canvas.draw de matplotlib et gobject.idle_add

 



 Mot :   Pseudo :  
 
Bas de page
Auteur Sujet :

problème avec de manager.canvas.draw de matplotlib et gobject.idle_add

n°2150614
Fauster
Posté le 24-07-2012 à 11:48:01  profilanswer
 

Bonjour à tous, j'ai besoin de votre aide. Je suis en stage dans le domaine de la mécanique des fluides (je ne suis pas informaticien de formation) et j'ai un problème sur mon code de résolution de l'équation de la chaleur à 2 dimensions avec terme advectif.
 
Le problème est le suisant:
 
J'ai un fluide (de l'eau) qui entre dans un canal à 20°C et qui est en contacte avec une paroi chauffante. Je résoud l'équation de la chaleur par différence finis en coordonnées cartésiennes. Mes conditions aux limites sont:
 
Une rampe de température de la paroi de 20°C à 100°C à en une seconde puis maintient à 100°C le reste du temps (bord gauche du maillage).
Une température d'entrée de 20°C
Une conditions de flux nul pour le bord droit du maillage.
En sortie je fais une extrapolation de la température gràce aux valeurs de températures précédentes
 
Voila mon code, je me suis inspiré d'un code que j'ai trouvé sur le net résolvant un problème de diffision de particule
 

Code :
  1. import scipy as sp
  2. import matplotlib
  3. matplotlib.use('GTKAgg')
  4. import gobject
  5. from pylab import *
  6. #variables du probleme
  7. dx=0.0001
  8. dy=0.001
  9. phof=1000.0
  10. cpf=4185.0
  11. lambdaf=0.6
  12. alphaf = lambdaf/(phof*cpf)
  13. timesteps=4000
  14. nx = int(0.01/dx)
  15. ny = int(0.2/dy)
  16. dx2=dx**2
  17. dy2=dy**2
  18. dtm=dx2*dy2/(2*alphaf*(dx2+dy2))
  19. vo = 1.0
  20. dtm_diff = min(dx2,dy2)*0.5/alphaf
  21. dtm_adv = dy/abs(vo)
  22. dt = 0.9 * min(dtm_diff, dtm_adv)
  23. ui = 20.0*ones([nx,ny])
  24. u = 20.0*ones([nx,ny])
  25. def updatefig(*args):
  26.    global u, ui, m
  27.    im.set_array(ui)
  28.    manager.canvas.draw()
  29.    print m
  30.    if (m*dt<1.0):
  31. #  conditions de bords
  32.       ui[0,:] = 20.0+80.0*m*dt
  33.       ui[:,0] = 20.0
  34.       u[0,:] = ui[0,:]
  35.       u[:,0] = ui[:,0]
  36.    
  37.       u[1:-1, 1:-1] = ui[1:-1, 1:-1] + alphaf*dt*( (ui[2:, 1:-1] - 2*ui[1:-1, 1:-1] + ui[:-2, 1:-1])/dx2 + (ui[1:-1, 2:] - 2*ui[1:-1, 1:-1] + ui[1:-1, :-2])/dy2 )-dt*vo*(ui[1:-1, 1:-1] - ui[1:-1, :-2])/dy
  38.  
  39.    else:
  40. #  conditions de bords
  41.       ui[0,:] = 100.00
  42.       ui[:,0] = 20.00
  43.       u[0,:] = ui[0,:]
  44.       u[:,0] = ui[:,0]
  45.       u[1:-1, 1:-1] = ui[1:-1, 1:-1] + alphaf*dt*( (ui[2:, 1:-1] - 2*ui[1:-1, 1:-1] + ui[:-2, 1:-1])/dx2 + (ui[1:-1, 2:] - 2*ui[1:-1, 1:-1] + ui[1:-1, :-2])/dy2 )-dt*vo*(ui[1:-1, 1:-1] - ui[1:-1, :-2])/dy        
  46.              
  47.    print u
  48.    ui = sp.copy(u)
  49.    m+=1
  50.    if m >= timesteps:
  51.       return False
  52.    return True
  53. fig = plt.figure()
  54. img = fig.add_subplot(1,1,1)
  55. ticks_at = [0,1]
  56. im = img.imshow( ui/100, interpolation='nearest', origin='lower')
  57. cbar = fig.colorbar(im,ticks=ticks_at)
  58. manager = get_current_fig_manager()
  59. m=1
  60. gobject.idle_add(updatefig)
  61. show()


 
Je n'arrive pas à visualiser correctement le champ de température. L 'écran est totalement noir. J'ai fait des prints pour regarder l'allure de mon tableau ui et les résultats sont physiques (le fluide se chauffe bien au contacte de la paroi) pourtant je n'arrive pas à les visualiser. Il faudrait que les zones chaudes se colorent en rouge au cours du calcul. Sauriez vous d'où vient le problème ?
 
Merci d'avance,
 
Fauster
 
edit du 25 Juillet : 8h13  prise en compte des remarques de Tentacle sur la stabilité et le schéma pour la dérivée d'ordre 1 (terme d'advection)
                                    Modification sur le calcul de u[i,j]. u[i,j] remplace u[1:-1, 1:-1] etc..
                                    Bonne prise en compte des conditions aux limites.
 
J'ai réussi à faire apparaître une image du champs de température (en bidouillant un peu)  mais je ne parvient toujours pas à itérer cette image au cours du calcul:
 
http://nsa29.casimages.com/img/2012/07/26/120726084000222156.png
 
PS: Le résultat est bien physique: apparition d'une couche limite de température pendant le transitoire. (champs de température normalisé, le fluide rentre à gauche et va vers la droite)
 
http://nsa29.casimages.com/img/2012/07/25/120725091918954233.png


Message édité par Fauster le 26-07-2012 à 13:58:08
mood
Publicité
Posté le 24-07-2012 à 11:48:01  profilanswer
 

n°2150629
Fauster
Posté le 24-07-2012 à 12:57:13  profilanswer
 

Hum, j'ai peut être des problèmes de convergences...j'éclaircis cela

n°2150638
theshockwa​ve
I work at a firm named Koslow
Posté le 24-07-2012 à 13:27:42  profilanswer
 

Si tu pouvais surtout encadrer ton code de la balise appropriée (balise "code" ), ce serait bien :)


---------------
last.fm
n°2150653
Fauster
Posté le 24-07-2012 à 14:01:22  profilanswer
 

Quoiqu'il en soit, lorsque le calcul converge, ça ne marche toujours pas. Je suis pas sur que ce soit ça la balise "code"


Message édité par Fauster le 24-07-2012 à 14:35:17
n°2150783
Tentacle
Posté le 25-07-2012 à 02:27:00  profilanswer
 

Bonjour,
 

Fauster a écrit :

Hum, j'ai peut être des problèmes de convergences...j'éclaircis cela


Si j'ai bien compris ton problème, tu as effectivement des problèmes de convergence mais ça ne résoudra peut-être pas ton problème d'affichage (quoique, en cas d'instabilité, les valeurs explosent ce qui donne un affichage saturé) :
 

  • dtm représente apparemment ton pas de temps maximal assurant la stabilité de ton schéma. Pour la diffusion, j'aurais plutôt mis min(dx2,dy2)*0.5/alphaf mais ton expression en est un minorant donc ok ; Le problème est que tu ne prends pas en compte la condition de stabilité de ton schéma d'advection qui devrait être (après correction du schéma) dt <= dy/abs(vo). Vu les paramètres de ton problème, tu devrais avoir dt <= 0.001 ce qui n'est pas le cas. Tu peux faire un truc du genre :
Code :
  1. dtm_diff = min(dx2,dy2)*0.5/alphaf
  2. dtm_adv = dy/abs(vo)
  3. dt = 0.9 * min(dtm_diff, dtm_adv)


le coefficient 0.9 étant là pour contrôler le pas de temps par rapport à sa limite maximale (sachant que plus il est bas, plus ton schéma sera artificiellement diffusif).
 

  • concernant ton schéma d'advection, il correspond à un schéma d'Euler explicite décentré mais il est décentré du mauvais côté : pour une vitesse positive, il faut décentrer à gauche (en amont par rapport au vecteur vitesse) et ici tu décentres à droite (le schéma est alors inconditionnellement instable). Il faudrait plutôt mettre :
Code :
  1. -dt*vo*(ui[1:-1, 1:-1] - ui[1:-1, :-2])/dy


à la place du dernier terme de ton schéma.
 
En espérant que ça fera avancer le smilblick ;)
 
PS : tu affiches l'ensemble du vecteur ui or sur les bords du domaine, il vaudra toujours 0 au moment de l'affichage car il est copié à partir de u qui lui-même n'est calculé que sur les indices [1:-1, 1:-1]. Soit tu n'affiches que cette sous-matrice, soit il faut compléter ces valeurs avec les contraintes de ton problème.
 
PS2 : Je ne m'y connais pas beaucoup en python mais si j'ai bien compris, la méthode idle_add de gobject appelle updatefig tant que celle-ci renvoie True. Tu devrais peut-être faire en sorte de pouvoir contrôler ton résultat numérique toutes les n itérations ?

n°2150789
Fauster
Posté le 25-07-2012 à 08:20:24  profilanswer
 

Merci à toi Tentacle. J'ai également lu un document qui expliquait pourquoi le terme advectif engendrait une instabilité vu son évaluation.  
 
J'ai un peu modifié le calcul de u|i,j] pour que cela soit plus clair, cela reste équivalent à l'ancienne version, bien que le calcul soit moins rapide.
 
Je cherche toujours à résoudre le problème d'affichage.  
 
J'ai tenté de normaliser mon profil de température-->ui/umax = ui/100 mais ecla ne fonctionne pas

Message cité 1 fois
Message édité par Fauster le 25-07-2012 à 08:30:28
n°2150968
Tentacle
Posté le 25-07-2012 à 15:22:24  profilanswer
 

Fauster a écrit :

Merci à toi Tentacle. J'ai également lu un document qui expliquait pourquoi le terme advectif engendrait une instabilité vu son évaluation.  
 
J'ai un peu modifié le calcul de u|i,j] pour que cela soit plus clair, cela reste équivalent à l'ancienne version, bien que le calcul soit moins rapide.


Le calcul vectoriel est probablement plus rapide que faire des boucles sous python. D'ailleurs pourquoi l'avoir modifié ?  
Sur le coup, tu mélanges des parties vectorielles (les conditions aux bords) avec le schéma numérique calculé en chaque i,j.
De plus, tu t'es trompé en convertissant le terme advectif :

Code :
  1. -dt*vo*(ui[i,j]- ui[i,j-1])/dy


Autre point : il y a un problème d'indentation dans ton code (ou du moins pour l'affichage sur le forum) ce qui le rend ambigu (surtout pour du python).
Comme on est dans les considérations de codage, tu pourrais ne mettre que la définition des conditions aux bords dans le test m*dt<0.1, le schéma numérique étant indépendant du temps.
 

Citation :


Je cherche toujours à résoudre le problème d'affichage.  
 
J'ai tenté de normaliser mon profil de température-->ui/umax = ui/100 mais ecla ne fonctionne pas


Tu devrais peut-être mettre plus en avant, dans ton post, la partie affichage (problème avec l'utilisation de manager.canvas.draw de matplotlib et gobject.idle_add) et reléguer au second plan la partie numérique. Il te reste quoi comme problème d'affichage ? C'est toujours noir ou c'est juste que tu ne vois pas l'évolution du calcul ?
Le titre n'est d'ailleurs pas du tout clair.
 
À part ça, l'image que tu obtiens ressemble bien à la solution stationnaire que tu devrais obtenir. En attendant de résoudre tes soucis d'affichage, tu devrais d'ailleurs diminuer la vitesse vo pour obtenir une solution stationnaire moins proche de l'état initial et ainsi mieux visualiser l'évolution de ta solution.

n°2151079
Fauster
Posté le 26-07-2012 à 08:16:36  profilanswer
 

Citation :

De plus, tu t'es trompé en convertissant le terme advectif :
 
-dt*vo*(ui[i,j]- ui[i,j-1])/dy


 
C'est corrigé, merci.
 

Citation :

Sur le coup, tu mélanges des parties vectorielles (les conditions aux bords) avec le schéma numérique calculé en chaque i,j.


 
Oui en effet...je modifie ça
 

Citation :

C'est toujours noir ou c'est juste que tu ne vois pas l'évolution du calcul ?


 
Toujours noir

Message cité 1 fois
Message édité par Fauster le 26-07-2012 à 11:55:02
n°2151130
Fauster
Posté le 26-07-2012 à 13:42:32  profilanswer
 

Le code suivant fonctionne:
 

Code :
  1. import scipy as sp
  2. import matplotlib
  3. matplotlib.use('GTKAgg')
  4. import gobject
  5. from pylab import *
  6. #variables du probleme
  7. dx=0.0001
  8. dy=0.001
  9. phof=1000.0
  10. cpf=4185.0
  11. lambdaf=0.6
  12. alphaf = lambdaf/(phof*cpf)
  13. timesteps=4000
  14. nx = int(0.01/dx)
  15. ny = int(0.2/dy)
  16. dx2=dx**2
  17. dy2=dy**2
  18. dtm=dx2*dy2/(2*alphaf*(dx2+dy2))
  19. vo = 1.0
  20. dtm_diff = min(dx2,dy2)*0.5/alphaf
  21. dtm_adv = dy/abs(vo)
  22. dt = 0.9 * min(dtm_diff, dtm_adv)
  23. ui = 20.0*ones([nx,ny])
  24. u = 20.0*ones([nx,ny])
  25. ui[0,:] = 100
  26. ui[:,0] = 20.0
  27. u[0,:] = ui[0,:]
  28. u[:,0] = ui[:,0]
  29. def updatefig(*args):
  30.    global u, ui, m
  31.    im.set_array(ui)
  32.    manager.canvas.draw()
  33.    print m
  34.    
  35.    u[1:-1, 1:-1] = ui[1:-1, 1:-1] + alphaf*dt*( (ui[2:, 1:-1] - 2*ui[1:-1, 1:-1] + ui[:-2, 1:-1])/dx2 + (ui[1:-1, 2:] - 2*ui[1:-1, 1:-1] + ui[1:-1, :-2])/dy2 )-dt*vo*(ui[1:-1, 1:-1] - ui[1:-1, :-2])/dy
  36.  
  37.    print u
  38.    ui = sp.copy(u)
  39.    m+=1
  40.    if m >= timesteps:
  41.       return False
  42.    return True
  43. fig = plt.figure()
  44. img = fig.add_subplot(1,1,1)
  45. ticks_at = [0,100]
  46. im = img.imshow( ui, interpolation='nearest', origin='lower')
  47. cbar = fig.colorbar(im,ticks=ticks_at)
  48. manager = get_current_fig_manager()
  49. m=1
  50. gobject.idle_add(updatefig)
  51. show()


 
Par contre ça ne fait pas ce que je veux. Je n'ai pas des conditions aux limites qui changent avec le pas de temps. Manifestement il y a un problème avec ça. Je rappel que je veux pour ui[0,:] une rampe de température de 0 à 1 seconde (20°C à 100) puis maintient en température à 100 °C le reste deu temps


Message édité par Fauster le 26-07-2012 à 14:04:48
n°2151140
Tentacle
Posté le 26-07-2012 à 14:32:42  profilanswer
 

Fauster a écrit :

Citation :

C'est toujours noir ou c'est juste que tu ne vois pas l'évolution du calcul ?


 
Toujours noir


 
En regardant la doc de imshow, la palette de couleur est automatiquement adaptée au min/max des données fournies mais uniquement au moment de l'appel de cette méthode. Dans ton code, quand imshow est appelée, ui est constant et vaut 20.0 ce qui donne une palette de couleurs réduite à sa plus petite valeur (noir).
 
Pour éviter ça, tu peux spécifier les bornes de tes données :

Code :
  1. im = img.imshow( ui, norm=matplotlib.colors.Normalize(vmin=20, vmax=100), cmap=cm.hot, interpolation='nearest', origin='lower')


 
PS: faut attendre un peu que la température monte pour voir quelquechose.


Message édité par Tentacle le 26-07-2012 à 14:35:44
mood
Publicité
Posté le 26-07-2012 à 14:32:42  profilanswer
 

n°2151145
Fauster
Posté le 26-07-2012 à 14:43:15  profilanswer
 

WOUOU ca marche. Merci encore à toi ;)

n°2151152
Tentacle
Posté le 26-07-2012 à 14:59:16  profilanswer
 

Fauster a écrit :

WOUOU ca marche. Merci encore à toi ;)


Au passage, je te conseille *fortement* ( pour l'avoir testé ) de repasser à la version vectorielle de ton code, tu devrais voir la différence de performance :p (tu peux aussi n'afficher le résultat que toutes les n itérations)
 
 :hello:

n°2151171
Fauster
Posté le 26-07-2012 à 15:58:08  profilanswer
 

Oui c'est franchement bcp bcp plus rapide


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

  problème avec de manager.canvas.draw de matplotlib et gobject.idle_add

 

Sujets relatifs
Variable Tableau qui ne passe pas en condition ( ! )[RESOLU] php script affichage fichier
[MS SQL Server] Problème de volumétrie : changement de type de colonneProblème formulaire vb/access 2010
Affichage title avec balise .attr[BATCH] probleme dans la redaction de mon batch
Problème lien animation flash / HTMLProblème d'affichage - taille des bordures
Plus de sujets relatifs à : problème avec de manager.canvas.draw de matplotlib et gobject.idle_add


Copyright © 1997-2018 Hardware.fr SARL (Signaler un contenu illicite / Données personnelles) / Groupe LDLC / Shop HFR