!************************************************************************! !Resolution de l'equation de la chaleur sur le domaine [0,1]x[0,1]x[0,1] ! discretisation spatiale : differences finies 2nd ordre ! discretisation temporelle : Crank-Nicholson, 2nd ordre ! avec comme solveur le Gradient Conjugue. ! d(u)/dt - \lambda * Delta u = f(x,y,z) ! u sur les bords vaut 0 ! La solution exacte est u = x*y*z*(x-1)*(y-1)*(z-1) ! ! Dans cette version, on se donne ! le nombre de points interieurs total suivant x (ntx), ! le nombre de points inferieurs total suivant y (nty) et ! le nombre de points inferieurs total suivant z (ntz). ! ntx =? nty =? ntz =? ! hx represente le pas suivant x, hx = 1 / (ntx+1), ! hy represente le pas suivant y, hy = 1 / (nty+1), ! hz represente le pas suivant z, hz = 1 / (ntz+1). ! ! Pour chaque processeur : ! 1) decomposer le domaine ! 2) connaitre ses 6 voisins ! 3) echanger les points aux interfaces ! 4) calculer ! ! Auteur : Isabelle DUPAYS (CNRS/IDRIS - France) ! ! Cree le : Fri nov 15 10:32:45 1996 ! !!! !!! ! ! Modifie le : Mon Sept 17 9:32:41 2001 ! ! Nature : passage du code en 3D ! conversion Poisson -> Chaleur ! conversion jacobi -> gradient conjugue ! ! Auteur : Guy Moebs (CRIHAN - France) ! ! !***********************************************************************! ! PROGRAM Chaleur ! USE MPI ! USE numerics ! IMPLICIT NONE ! INCLUDE 'comthd.h' ! INTERFACE ! INCLUDE 'i_communication.h' ! INCLUDE 'i_cree_type.h' ! INCLUDE 'i_domaine.h' ! INCLUDE 'i_fctfx.h' ! INCLUDE 'i_scdmb.h' ! INCLUDE 'i_gradconj.h' ! INCLUDE 'i_nrmerr.h' ! INCLUDE 'i_voisinage.h' ! END INTERFACE ! !---Declaration des variables ! !Solution u REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: u ! !Solution exacte REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: u_exact ! !Second membre REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: b ! !Forcage exterieur REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: f1 REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: f2 ! !Vecteurs de travail REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: p REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: q REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: r ! ! Tableaux des repartitions INTEGER, DIMENSION(:), ALLOCATABLE :: npts INTEGER, DIMENSION(:), ALLOCATABLE :: ntproc ! !Nombre iterations en temps INTEGER :: itg ! !Nombre d'iterations maximum gradient conjugue INTEGER :: it_max ! !Nombre d'iterations en temps INTEGER :: nbiter ! !Nombre de threads INTEGER :: ir, iq INTEGER :: iproc ! !$ EXTERNAL OMP_GET_NUM_THREADS !$ INTEGER :: OMP_GET_NUM_THREADS ! !$ EXTERNAL OMP_GET_THREAD_NUM !$ INTEGER :: OMP_GET_THREAD_NUM ! !Scalaires gradient conjugue REAL(rp) :: rnorm_k INTEGER :: itgc INTEGER :: nbitgc ! !Scalaire Crank Nicholson REAL(rp) :: rnudt2 ! !Normes diverses REAL(rp) :: rnrm2 REAL(rp) :: rnrmoo REAL(rp) :: rerr2 REAL(rp) :: rerroo ! !Temps courant REAL(rp) :: tps ! !Frequence des sorties INTEGER :: ifreq ! !Consommation memoire REAL(rp) :: gmem ! !Convergence gradient conjugue REAL(rp) :: prec ! !Consommation memoire REAL(rp), DIMENSION(:), ALLOCATABLE :: igmem ! !Tableau contenant les voisins du domaine courant INTEGER, PARAMETER :: NB_VOISINS = 6 ! INTEGER, DIMENSION(NB_VOISINS) :: voisin INTEGER, PARAMETER :: N=1, E=2, S=3 INTEGER, PARAMETER :: W=4, AV=5, AR=6 ! !Compteur INTEGER :: j ! !Mesure du temps REAL(rp) :: t1, t2 ! !Volume memoire REAL(rp) :: ilmem ! !Test de convergence LOGICAL :: convergence ! !Constantes MPI ! INTEGER :: imyrank INTEGER :: comm3d INTEGER :: ip_xy, ip_xz, ip_yz INTEGER, DIMENSION(3) :: dims LOGICAL, DIMENSION(3) :: periods ! !***********************************************************************! ! !Initialisation de MPI CALL MPI_INIT( ierr ) ! !Savoir quel processeur je suis CALL MPI_COMM_RANK( MPI_COMM_WORLD, myrank, ierr ) ! !Connaitre le nombre total de processeurs CALL MPI_COMM_SIZE( MPI_COMM_WORLD, nbproc, ierr ) ! !Lecture des parametres machep = MAX( EPSILON(machep), 1E-15_rp ) ! IF ( myrank == 0 ) THEN OPEN(19, FILE='C3D.dat', STATUS='OLD', FORM='FORMATTED') READ(19, 101) nbthd READ(19, 101) ntx READ(19, 101) nty READ(19, 101) ntz READ(19, 101) nbiter READ(19, 102) dt READ(19, 102) rlambda READ(19, 102) prec READ(19, 101) it_max READ(19, 101) ifreq CLOSE(19) ! 101 FORMAT ( I5 ) 102 FORMAT ( E14.4 ) ! WRITE (6,201) 'Equation de la chaleur 3D' WRITE (6,202) 'machep = ', machep WRITE (6,201) 'Nombre de threads OpenMP : ', nbthd WRITE (6,201) 'Nombre de points en X : ', ntx WRITE (6,201) 'Nombre de points en Y : ', nty WRITE (6,201) 'Nombre de points en Z : ', ntz ! WRITE (6,201) 'Nombre d iterations en temps : ', nbiter WRITE (6,202) 'Pas de temps : ', dt WRITE (6,202) 'Viscosite cinematique : ', rlambda ! WRITE (6,202) 'Precision gradient : ', prec WRITE (6,201) 'Nbre max d iterations : ', it_max WRITE (6,201) 'Frequence des sorties : ', ifreq ! 201 FORMAT ( A, I5 ) 202 FORMAT ( A, E12.4 ) ! WRITE (6,1010) 'Nombre de processus MPI : ', nbproc ! IF ( rlambda < machep ) THEN WRITE (6,*) 'scalaire lambda positif strictement' CALL MPI_ABORT ( MPI_COMM_WORLD, 1, ierr ) STOP END IF ! END IF ! !Type de donnees IF ( rp == dp ) THEN ITYPE_REAL = MPI_DOUBLE_PRECISION ELSE ITYPE_REAL = MPI_REAL END IF ! !Envoi de la valeur de ntx a tous les autres processeurs CALL MPI_BCAST( nbthd, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( ntx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( nty, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( ntz, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( nbiter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( it_max, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( ifreq, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( dt, 1, ITYPE_REAL, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( rlambda, 1, ITYPE_REAL, 0, MPI_COMM_WORLD, ierr ) CALL MPI_BCAST( prec, 1, ITYPE_REAL, 0, MPI_COMM_WORLD, ierr ) ! !Connaitre ! le nombre de processeurs selon x ! le nombre de processeurs selon y ! le nombre de processeurs selon z en fonction du nombre total de processeurs dims(:) = 0 IF ( ntx == 1 ) dims(1) = 1 IF ( nty == 1 ) dims(2) = 1 IF ( ntz == 1 ) dims(3) = 1 CALL MPI_DIMS_CREATE( nbproc, 3, dims, ierr ) ! !Creation de la grille de processeurs 3D sans periodicite periods(:) = .FALSE. ! CALL MPI_CART_CREATE( MPI_COMM_WORLD, 3, dims, periods, .TRUE., & comm3d, ierr) ! !Recherche de ses 6 voisins pour chaque processeur CALL voisinage( comm3d, voisin ) ! !Determiner les indices de chaque sous-domaine CALL domaine( comm3d ) ! ! Consommation memoire en ko ilmem = 0.0625_rp * REAL(ex-sx+3, rp) * REAL(ey-sy+3, rp) * REAL(ez-sz+3, rp) ALLOCATE( igmem (0:nbproc-1), STAT=ierr) igmem (0:nbproc-1) = 0 ! CALL MPI_GATHER( ilmem, 1, ITYPE_REAL, & igmem, 1, ITYPE_REAL, & 0, comm3d, ierr ) ! CALL MPI_COMM_RANK( comm3d, imyrank, ierr) IF ( imyrank == 0 ) THEN WRITE (6,300) WRITE (6,301) WRITE (6,300) DO j = 0, nbproc - 1 WRITE (6,302) j, igmem(j) END DO WRITE (6,300) WRITE (6,303) SUM( igmem(0:nbproc-1) ) WRITE (6,300) WRITE (6,*) 300 FORMAT ( '+-----------+------------+') 301 FORMAT ( '| Processus | Mem. (ko) |') 302 FORMAT ( '| ', I6, ' |', F12.1, '|') 303 FORMAT ( '| TOTAL |', F12.1, '|') END IF ! !Allocation dynamique des tableaux ALLOCATE( u (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & b (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & f1 (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & f2 (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & p (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & q (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & r (sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & u_exact(sx-1:ex+1, sy-1:ey+1, sz-1:ez+1), & npts (0:nbthd-1), ntproc (0:nbthd), STAT = ierr ) IF ( ierr /= 0 ) THEN WRITE (6,1010) 'rang = ', myrank WRITE (6,1010) 'ierr = ', ierr WRITE (6,1010) 'Probleme allocation' CALL MPI_ABORT ( MPI_COMM_WORLD, 1, ierr ) 1010 FORMAT ( A, I6 ) STOP END IF ! !Repartition des iterations selon z iq = (ez - sz + 1) / nbthd ! DO iproc = 0, nbthd - 1 npts(iproc) = iq END DO ! ir = (ex - sx + 1) - nbthd * iq IF ( ir /= 0 ) THEN DO iproc = 0, ir - 1 npts(iproc) = npts(iproc) + 1 END DO END IF ! ntproc (0) = 1 DO iproc = 0, nbthd - 1 ntproc(iproc+1) = ntproc(iproc) + npts(iproc) END DO ! write (6,*) 'npts : ', npts(0:nbthd-1) write (6,*) 'ntproc : ', ntproc(0:nbthd) ! rnudt2 = 0.5_rp * rlambda * dt ! itgc = 0 nbitgc = 0 ! rnorm_k = zero tps = zero ! !Creation du type type_ligne pour echanger les points a l'ouest et a l'est CALL cree_type( ip_xy, ip_xz, ip_yz ) ! !************************************************************************! ! !Debut de la region parallele ! CALL OMP_SET_NUM_THREADS ( nbthd ) ! !$OMP PARALLEL DEFAULT( NONE ) & !$OMP SHARED( u, b, u_exact, f1, f2, p, q, r ) & !$OMP SHARED( tps, dt, rnorm_k, itgc, nbitgc ) & !$OMP SHARED( rnrm2, rnrmoo, rerr2, rerroo ) & !$OMP SHARED( nbproc, ntproc, nbthd, nbiter, prec, machep ) & !$OMP SHARED( ifreq, it_max, rnudt2, rlambda ) & !$OMP SHARED( ip_xy, ip_xz, ip_yz, voisin, comm3d ) & !$OMP SHARED( t2, t1, myrank ) & !$OMP PRIVATE( itg ) ! !$OMP MASTER nbthd = 1 !$ nbthd = OMP_GET_NUM_THREADS () WRITE (6,1040) 'Processus MPI : ', myrank, ' Nombre de threads OpenMP : ', nbthd !$OMP END MASTER 1040 FORMAT ( 2(A, I5), /) ! mythd = 0 !$ mythd = OMP_GET_THREAD_NUM () ! kmin = ntproc(mythd) kmax = ntproc(mythd+1) - 1 ! !Initialisation des matrices CALL fctfx( u, f1, f2, p, u_exact, tps ) ! !Erreurs initiales CALL nrmerr( rnrm2, rnrmoo, rerr2, rerroo, u, u_exact, tps, comm3d ) ! IF ( myrank == 0 ) THEN !$OMP MASTER WRITE (6, 401) 'TPS', 'uL2', 'uLoo', 'eL2', 'eLoo', 'itGC', 'res' WRITE (6, 402) tps, rnrm2, rnrmoo, rerr2, rerroo, itgc, SQRT(rnorm_k) !$OMP END MASTER 401 FORMAT ( 4x, A3, 7x, A3, 8x, A4, 9x, A3, 8x, A4, 7x, A4, 5x, A3 ) 402 FORMAT ( F8.4, 4E12.4, 4x, I4, E12.4 ) END IF ! !Mesure du temps en seconde dans la boucle en temps !$OMP MASTER t1 = MPI_WTIME() !$OMP END MASTER ! !================================= ! Debut de la boucle en temps !================================= ! DO itg = 1, nbiter ! !Echange des points aux interfaces pour u !$OMP MASTER CALL communication( u, ip_xy, ip_xz, ip_yz, voisin, comm3d ) ! nbitgc = nbitgc + itgc itgc = 0 rnorm_k = zero ! !$OMP END MASTER !$OMP BARRIER ! !Calcul second membre Crank Nicholson CALL scdmb( b, f1, f2, u, tps, rnudt2 ) ! !Avancement en temps !$OMP MASTER tps = tps + dt !$OMP END MASTER ! !Resolution gradient conjugue CALL gradconj( u, b, p, q, r, it_max, prec, itgc, rnorm_k, & rnudt2, ip_xy, ip_xz, ip_yz, voisin, comm3d ) ! !Calcul des normes IF ( MOD(itg,ifreq) == 0 ) THEN ! CALL nrmerr( rnrm2, rnrmoo, rerr2, rerroo, u, u_exact, tps, comm3d ) ! IF ( myrank == 0 ) THEN !$OMP MASTER WRITE (6, 402) tps, rnrm2, rnrmoo, rerr2, rerroo, itgc, SQRT(rnorm_k) !$OMP END MASTER ! END IF ! END IF ! !================================= ! Fin de la boucle en temps !================================= ! END DO ! !Mesure du temps a la sortie de la boucle !$OMP MASTER t2 = MPI_WTIME() - t1 !$OMP END MASTER ! !************************************************************************! ! !Fin de la region parallele !$OMP END PARALLEL ! nbitgc = nbitgc + itgc ! !Affichage du temps de convergence par le processeur 0 IF ( myrank == 0 ) THEN t1 = REAL(ntx,rp) * REAL(nty,rp) * REAL(ntz,rp) * ( & 23.0E-06_rp * REAL(nbitgc,rp) & ! gradconj + 30.0E-06_rp * REAL(nbiter) & ! scdmb + 12.0E-06_rp * REAL(nbiter/ifreq + 1, rp) ) / t2 ! nrmerr WRITE (6,202) 'Temps ecoule : ', t2 WRITE (6,203) 'Performances cumulees : ', t1, ' MFlops' WRITE (6,203) 'Performances par processus-thread : ', t1 / REAL(nbproc*nbthd,rp), ' MFlops' ! 203 FORMAT ( A, F8.2, A ) ! END IF ! !Desallocation des tableaux DEALLOCATE( u, u_exact, b, f1, f2, p, q, r, igmem, npts, ntproc ) ! !Liberation du type type_ligne et du communicateur comm2d CALL MPI_TYPE_FREE( ip_xy, ierr ) CALL MPI_TYPE_FREE( ip_xz, ierr ) CALL MPI_TYPE_FREE( ip_yz, ierr ) CALL MPI_COMM_FREE( comm3d, ierr ) ! !Desactivation de MPI CALL MPI_FINALIZE(ierr) ! END PROGRAM Chaleur