+ Répondre à la discussion
Affichage des résultats 1 à 3 sur 3
  1. #1
    Candidat au titre de Membre du Club
    Inscrit en
    avril 2009
    Messages
    46
    Détails du profil
    Informations forums :
    Inscription : avril 2009
    Messages : 46
    Points : 10
    Points
    10

    Par défaut Bug avec un code de prolongement vers le haut

    Bonjour à tous

    Je me permets de solliciter votre aide parce que j'ai un souci avec un programme en Fortran que j'ai essayé de faire en me basant sur un autre programme existant.

    Je souhaite un code qui me calcule un prolongement vers le haut d'un champ potentiel (En utilisant des transformées de Fourier). Le code que j'ai écrit est le suivant:

    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
          program upward
    
    c Déclarations
    
          character name*70
          real rft(512,512), flo(512,512), fla(512,512)
    
        1 format (a)
        2 format (f11.6,x,f10.6,x,f10.2)
    
          pi = 3.1415926535
          tpi = 2. * pi
          rad = pi / 180.
    
          dimension grid(nx*ny),store(2*nx*ny),nn(2)      
          real kx, ky, k
          complex cgrid, cmplx
          data pi/3.14159265/
          index(i,j,ncol)=(j-1)*ncol+i
          nn(1) = ny
          nn(2) = nx
          ndim = 2
          dkx = 2.*pi/(nx*dx)
          dky = 2.*pi/(ny*dy)
    
    c Lecture des valeurs
    
          print*,'Nom de la grille d''anomalies magnetiques?'
          read (5,1) name
          print*,'Nombre de colonnes de la grille, en 2**n?'
          read (5,*) ny
          ny = 2 ** ny
          print*,'Nombre de lignes de la grille, en 2**n?'
          read (5,*) nx
          nx = 2 ** nx
          print*,'Pas d''echantillonnage en x?'
          read (5,*) dx
          print*,'Pas d''échantillonnage en y?'
          read (5,*) dy
          print*,'Hauteur du prolongement vers le haut?'
          read (5,*) dz
          
          store = 2. * nx * ny
    
    c Programme
    
          eu = tpi / float(nx)
          ev = tpi / float(ny)
          nx2 = nx / 2 + 1
          ny2 = ny / 2 + 1
    
          open (70, file = name)
          do 3 i = 1, nx
          do 3 j = 1, ny
          read (70,*) flo(nx-i+1,j), fla(nx-i+1,j), rft(nx-i+1,j)
        3 continue
          close (70)
    
          do 4 i = 1, nx
          do 4 j = 1, ny
    
        4 continue
    
          do 10 j=1,nx
          do 10 i=1,ny
          ij=index(i,j,ny)
          store(2*ij-1)=grid(ij)
       10 store(2*ij)=0.
          call fourn(store,nn,ndim,-1)
          do 20 j=1,nx
          do 20 i=1,ny
          ij=index(i,j,ny)
          call kvalue(i,j,nx,ny,dkx,dky,kx,ky)
          k=sqrt(kx**2+ky**2)
          cont=exp(-k*dz)
          cgrid=cmplx(store(2*ij-1),store(2*ij))*cont
          store(2*ij-1)=real(cgrid)
       20 store(2*ij)=aimag(cgrid)
          endif
        5 continue
    
          call fourn(store,nn,ndim,+1)
          do 30 j=1,nx
          do 30 i=1,ny
          ij=index(i,j,ny)
       30 grid(ij)=store(2*ij-1)/(nx*ny)
    
          print*,'Nom de la grille prolongée vers le haut?'
          read (5,1) name
          open (70, file = name)
          do 6 i = 1, nx
          do 6 j = 1, ny
          rft(i,j) = rft(i,j) / acof * 2.
        6 continue
          do 7 i = 1, nx
          do 7 j = 1, ny
          write (70,2) flo(nx-i+1,j), fla(nx-i+1,j), rft(nx-i+1,j)
        7 continue
          close (70)
    
          stop
          end
    Bien entendu, comme tout programme qui se respecte, il bugge

    Je précise que je suis complètement nul en programmation. J'ai essayé de faire ce que je pouvais mais il continue de m'afficher toute une gamme de messages d'erreurs quand je le compile et je ne comprends pas forcément ce que je dois faire pour que ça fonctionne.

    Je vous mets les subroutines avec lesquelles il fonctionne:

    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
    subroutine fourn(data,nn,ndim,isign)
    c
    c  Replaces DATA by its NDIM-dimensional discrete Fourier transform, 
    c  if ISIGN is input as 1.  NN is an integer array of length NDIM, 
    c  containing the lengths of each dimension (number of complex values), 
    c  which must all be powers of 2.  DATA is a real array of length twice 
    c  the product of these lengths, in which the data are stored as in a 
    c  multidimensional complex Fortran array.  If ISIGN is input as -1, 
    c  DATA is replaced by its inverse transform times the product of the
    c  lengths of all dimensions.  From Press, W.H., Flannery, B.P., 
    c  Teukolsky, S.A., and Vetterling, W.T., 1986, Numerical Recipes, 
    c  Cambridge Univ. Press, p. 451-453.
    c
          real*8 wr,wi,wpr,wpi,wtemp,theta
          dimension nn(ndim),data(*)
          ntot=1
          do 11 iidim=1,ndim
          ntot=ntot*nn(iidim)
       11    continue
          nprev=1
          do 18 iidim=1,ndim
            n=nn(iidim)
            nrem=ntot/(n*nprev)
            ip1=2*nprev
            ip2=ip1*n
            ip3=ip2*nrem
            i2rev=1
            do 14 i2=1,ip2,ip1
              if(i2.lt.i2rev)then
                do 13 i1=i2,i2+ip1-2,2
                  do 12 i3=i1,ip3,ip2
                    i3rev=i2rev+i3-i2
                    tempr=data(i3)
                    tempi=data(i3+1)
                    data(i3)=data(i3rev)
                    data(i3+1)=data(i3rev+1)
                    data(i3rev)=tempr
                    data(i3rev+1)=tempi
       12         continue
       13       continue
              endif
              ibit=ip2/2
       1      if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then
                i2rev=i2rev-ibit
                ibit=ibit/2
              go to 1
              endif
              i2rev=i2rev+ibit
       14   continue
            ifp1=ip1
       2    if(ifp1.lt.ip2)then
              ifp2=2*ifp1
              theta=isign*6.28318530717959d0/(ifp2/ip1)
              wpr=-2.d0*dsin(0.5d0*theta)**2
              wpi=dsin(theta)
              wr=1.d0
              wi=0.d0
              do 17 i3=1,ifp1,ip1
                do 16 i1=i3,i3+ip1-2,2
                  do 15 i2=i1,ip3,ifp2
                    k1=i2
                    k2=k1+ifp1
                    tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1)
                    tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2)
                    data(k2)=data(k1)-tempr
                    data(k2+1)=data(k1+1)-tempi
                    data(k1)=data(k1)+tempr
                    data(k1+1)=data(k1+1)+tempi
       15         continue
       16       continue
                wtemp=wr
                wr=wr*wpr-wi*wpi+wr
                wi=wi*wpr+wtemp*wpi+wi
       17     continue
              ifp1=ifp2
            go to 2
            endif
            nprev=n*nprev
       18 continue
          return
          end
    Et

    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
    subroutine kvalue(i,j,nx,ny,dkx,dky,kx,ky)
    c  Subroutine KVALUE finds the wavenumber coordinates of one 
    c  element of a rectangular grid from subroutine FOURN.  
    c
    c  Input parameters:
    c    i  - index in the ky direction.
    c    j  - index in the kx direction.
    c    nx - dimension of grid in ky direction (a power of two).
    c    ny - dimension of grid in kx direction (a power of two).
    c    dkx - sample interval in the kx direction.
    c    dky - sample interval in the ky direction.
    c
    c  Output parameters:
    c    kx - the wavenumber coordinate in the kx direction.
    c    ky - the wavenumber coordinate in the ky direction.
    c
          real kx,ky
          nyqx=nx/2+1
          nyqy=ny/2+1
          if(j.le.nyqx)then
             kx=(j-1)*dkx
             else
                kx=(j-nx-1)*dkx
                end if
          if(i.le.nyqy)then
            ky=(i-1)*dky
            else
                ky=(i-ny-1)*dky
                end if
          return
          end
    Voilà voilà.

    Je vous remercie par avance pour toute réponse qui pourrait me sortir de ma galère actuelle

  2. #2
    Modérateur

    Profil pro
    Inscrit en
    août 2006
    Messages
    844
    Détails du profil
    Informations personnelles :
    Localisation : Canada

    Informations forums :
    Inscription : août 2006
    Messages : 844
    Points : 1 178
    Points
    1 178

    Par défaut

    Je n'ai pas le temps de garder à fond, mais je vois rapidement certaines erreurs :
    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
          program upward
    
    c Déclarations
    
          character name*70
          real rft(512,512), flo(512,512), fla(512,512)
    
        1 format (a)
        2 format (f11.6,x,f10.6,x,f10.2)
    
          pi = 3.1415926535
          tpi = 2. * pi
          rad = pi / 180.
    
          dimension grid(nx*ny),store(2*nx*ny),nn(2)      
          real kx, ky, k
          complex cgrid, cmplx
          data pi/3.14159265/
    ...
    Les lignes "1 format" à "rad = " doivent être déplacées juste après la ligne "data" car ce sont des lignes exécutables et en Fortran, les lignes exécutables doivent suivre les lignes de déclaration. De plus, pi est assigné 2 fois (pi = et data pi...).

    Ensuite, tu déclares des tableaux (grid et store) dont tu ne connais pas la taille. Tu dois les délacer comme suit :
    Code :
    real, allocatable :: grid(:),store(:)
    ou (au choix) :
    Code :
    real, allocatable, dimension(:) :: grid,store
    Quand nx et ny seront connus (donc après les "read"), tu devras allouer les tableaux :
    Code :
    allocate (grid(nx*ny),store(2*nx*ny))
    Finalement, comme ton programme est petit et que tu utilises 2 routines, je te conseille d'ajouter une section "contains" à la fin de ton programme, juste avant le "end" et y insérer les 2 routines (attention: les end à la fin des subroutines doivent être remplacés par des end subroutine) . Le compilateurs fera alors un meilleur diagnostic des erreurs :
    Code :
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
          close (70)
    
          stop
    
          contains
    
          subroutine fourn(data,nn,ndim,isign)
    c
    c  Replaces DATA by its NDIM-dimensional discrete Fourier transform, 
    ...
       18 continue
          return
          end subroutine
    
          subroutine kvalue(i,j,nx,ny,dkx,dky,kx,ky)
    c  Subroutine KVALUE finds the wavenumber coordinates of one 
    ...
                end if
          return
          end subroutine
    
          end

  3. #3
    Membre éprouvé
    Inscrit en
    mars 2007
    Messages
    373
    Détails du profil
    Informations forums :
    Inscription : mars 2007
    Messages : 373
    Points : 492
    Points
    492

    Par défaut

    Bonjour,
    Citation Envoyé par Flo Flo Voir le message

    Bien entendu, comme tout programme qui se respecte, il bugge

    Je précise que je suis complètement nul en programmation. J'ai essayé de faire ce que je pouvais mais il continue de m'afficher toute une gamme de messages d'erreurs quand je le compile et je ne comprends pas forcément ce que je dois faire pour que ça fonctionne.
    Je ne peux que te recommander chaudement d'utiliser les options de compilation du compilateur qui permettent justement d'identifier les problèmes et erreurs "de base" (variable utilisée alors que non initialisée, tableau non alloué, dépassement d'indice dans un tableau, etc.).
    Par exemple avec le compilateur gfortran: "-g -O2 -Wall -fbounds-check -ffpe-trap=invalid,zero,overflow".
    Consulter leur doc pour les autres compilateurs.

    Bonne continuation.

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
  •