Publicité
+ Répondre à la discussion
Affichage des résultats 1 à 3 sur 3
  1. #1
    Membre du Club
    Inscrit en
    mars 2003
    Messages
    191
    Détails du profil
    Informations forums :
    Inscription : mars 2003
    Messages : 191
    Points : 61
    Points
    61

    Par défaut Boucles imbriquées, intégrale double, tableaux numpy

    Salut !


    J'ai le code suivant, dans lequel je n'arrive pas à remplacer la double boucle imbriquée par une opération sur des tableaux numpy. La double boucle est beaucoup trop lente (je dois appeler la fonction des centaines de fois et les tableaux font 1000x800 :

    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
     
        def GetFlux(self, time):
     
                        bx = self.GetField("bx", time) * self.wpewce
                        by = self.GetField("by", time) * self.wpewce
                        bz = self.GetField("bz", time) * self.wpewce
     
     
     
                        flux  = np.zeros((self.ncells[0]+1,self.ncells[1]+1),"float32", order='FORTRAN')
                        flux2  = np.zeros((self.ncells[0]+1,self.ncells[1]+1),"float32", order='FORTRAN')
     
                        dx = self.dl[0]
                        dz = self.dl[1]
     
                        nx = self.ncells[0]
                        nz = self.ncells[1]
     
                        j = 0
     
                       for i in np.arange(1, nx):
                               flux2[i,0] = flux2[i-1,0] + bz[i-1,0]*dx
     
                        flux[1:,0] = flux[0,0] + np.cumsum(bz[:-1,0]*dx)
     
     
     
                       for j in np.arange(1,nz):
                               flux2[0,j] = flux2[0,j-1] - bx[0,j-1]*dz
     
                        flux[0,1:] = flux[0,0] - np.cumsum(bx[0,:-1]*dz)
     
     
                       for i in np.arange(1,nx):
                               for j in np.arange(1,nz):
                                       flux2[i,j] = 0.5*(flux2[i-1,j] + bz[i-1,j]*dx) + 0.5*(flux2[i,j-1] - bx[i,j-1]*dz)
     
     
                        # flux = ???
     
     
                        return flux2

    (voir la me chose ici : http://pastebin.com/0ZEr6hKL)


    Le but est de calculer l'intégrale double donnant F a partir de F=dB_z/dx et F=-dB_x/dz.


    'flux2' est calculé à partir de simples boucles, comme on pourrait le faire en C, et contient le bon résultat.

    'flux' devrait contenir la meme chose, mais fait a partir d'opérations sur tableaux avec un fancy indexing. Seulement je ne parviens pas a faire la double boucle.


    Attention, c'est une intégrale, donc dans les boucles;, c'est bien flux2[i,j] qui dépend de flux2[i-1,j] et flux2[i,j-1], i.e. le tableau dépend de lui même à litération d'avant... c'est ce qui me pose problème pour la boucle 2D.

    Dans les boucles simples plus haut, j'ai trouvé cumsum() pour bien faire le travail, et les cases de tableau concernées sont égales à celles obtenues dans les boucles simples.

    Je me dis qu'il doit exister une méthode basée sur cumsum pour la boucle imbriquée ?


    Merci d'avance bcp pour votre aide
    --
    Heimdall

  2. #2
    Membre éprouvé
    Homme Profil pro Thomas Pegot
    Ingénieur développement logiciels
    Inscrit en
    janvier 2012
    Messages
    274
    Détails du profil
    Informations personnelles :
    Nom : Homme Thomas Pegot
    Localisation : France, Rhône (Rhône Alpes)

    Informations professionnelles :
    Activité : Ingénieur développement logiciels
    Secteur : High Tech - Éditeur de logiciels

    Informations forums :
    Inscription : janvier 2012
    Messages : 274
    Points : 424
    Points
    424

    Par défaut

    Bonjour,

    Avez-vous pensez à un produit de convolution ( convolve) avec dans ce cas un masque 3x3 pout flux2 3x1 et 1x3 pour bz et bx.

  3. #3
    Membre expérimenté
    Homme Profil pro arnaud
    Ingénieur développement logiciels
    Inscrit en
    janvier 2013
    Messages
    280
    Détails du profil
    Informations personnelles :
    Nom : Homme arnaud
    Localisation : France

    Informations professionnelles :
    Activité : Ingénieur développement logiciels
    Secteur : Conseil

    Informations forums :
    Inscription : janvier 2013
    Messages : 280
    Points : 502
    Points
    502

    Par défaut équation de Maxwell-Ampère

    Bonjour.

    Je n'ai pas compris comment Heimdall a obtenu l'équation sur flux2. Alors j'ai déterminé le flux en intégrant l'équation de Maxwell-Ampère.
    En discrétisant, on montre que le flux est proportionnel à :
    Code :
    dz * np.cumsum( bz, axis=0) - dx * np.cumsum( bx, axis=1)
    Cette expression est bien évidemment différente de celle donnée pour calculer flux2.

    Pour les comparer, je fais l'hypothèse que le champ magnétique est suivant X (Bz=0).
    Je modélise le champ par une gaussienne 2D, de révolution autour de l'axe Y.
    La gaussienne est placée au centre du domaine. Soient larg l'écart type, et maxi la valeur maximale des abscisses. Je note r le rapport larg sur maxi.

    Pour r=3, la gaussienne est "étendue" par rapport au domaine, le champ tend à être uniforme (Fig. ChampPresqueUniforme).
    Lorsque le champ est rigoureusement uniforme, les résultats sont les mêmes : Bx a une valeur constante dans les plans perpendiculaires à l'axe X.
    Un léger décalage apparait à r=3 avec la formule de Heimdall vers les Z>0. Avec ma formule, les courbes d'iso-valeurs sont de révolution autour de l'axe X.

    Pour r=1/3, la gaussienne est peu "étendue" sur le domaine, avec un comportement asymptotique de champ nul. (Fig. ChampPresqueNul).
    Avec ma méthode les courbes d'iso-valeurs restent de révolution autour de X, tandis qu'avec celle de Heimdall, elles sont clairement orientées suivant un axe faisant un angle de 45° avec X.

    Qui a raison ?
    Quelqu'un connaitrait-il un cas théorique permettant de trancher ?
    Images attachées Images attachées

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •