Hello,
>
> i am still trying to understand the parallelized version of the heat
> equation 2D solving that we saw at school. In order to explain my problem, i
> need to list the main code :
>
> 9 program heat
> 10
> !**************************************************************************
> 11 !
> 12 ! This program solves the heat equation on the unit square
> [0,1]x[0,1]
> 13 ! | du/dt - Delta(u) = 0
> 14 ! | u/gamma = cste
> 15 ! by implementing a explicit scheme.
> 16 ! The discretization is done using a 5 point finite difference scheme
> 17 ! and the domain is decomposed into sub-domains.
> 18 ! The PDE is discretized using a 5 point finite difference scheme
> 19 ! over a (x_dim+2)*(x_dim+2) grid including the end points
> 20 ! correspond to the boundary points that are stored.
> 21 !
> 22 ! The data on the whole domain are stored in
> 23 ! the following way :
> 24 !
> 25 ! y
> 26 ! ------------------------------------
> 27 ! d | |
> 28 ! i | |
> 29 ! r | |
> 30 ! e | |
> 31 ! c | |
> 32 ! t | |
> 33 ! i | x20 |
> 34 ! o /\ | |
> 35 ! n | | x10 |
> 36 ! | | |
> 37 ! | | x00 x01 x02 ... |
> 38 ! | ------------------------------------
> 39 ! -------> x direction x(*,j)
> 40 !
> 41 ! The boundary conditions are stored in the following submatrices
> 42 !
> 43 !
> 44 ! x(1:x_dim, 0) ---> left temperature
> 45 ! x(1:x_dim, x_dim+1) ---> right temperature
> 46 ! x(0, 1:x_dim) ---> top temperature
> 47 ! x(x_dim+1, 1:x_dim) ---> bottom temperature
> 48 !
> 49
> !**************************************************************************
> 50 implicit none
> 51 include 'mpif.h'
> 52 ! size of the discretization
> 53 integer :: x_dim, nb_iter
> 54 double precision, allocatable :: x(:,:),b(:,:),x0(:,:)
> 55 double precision :: dt, h, epsilon
> 56 double precision :: resLoc, result, t, tstart, tend
> 57 !
> 58 integer :: i,j
> 59 integer :: step, maxStep
> 60 integer :: size_x, size_y, me, x_domains,y_domains
> 61 integer :: iconf(5), size_x_glo
> 62 double precision conf(2)
> 63 !
> 64 ! MPI variables
> 65 integer :: nproc, infompi, comm, comm2d, lda, ndims
> 66 INTEGER, DIMENSION(2) :: dims
> 67 LOGICAL, DIMENSION(2) :: periods
> 68 LOGICAL, PARAMETER :: reorganisation = .false.
> 69 integer :: row_type
> 70 integer, parameter :: nbvi=4
> 71 integer, parameter :: S=1, E=2, N=3, W=4
> 72 integer, dimension(4) :: neighBor
> 73
> 74 !
> 75 intrinsic abs
> 76 !
> 77 !
> 78 call MPI_INIT(infompi)
> 79 comm = MPI_COMM_WORLD
> 80 call MPI_COMM_SIZE(comm,nproc,infompi)
> 81 call MPI_COMM_RANK(comm,me,infompi)
> 82 !
> 83 !
> 84 if (me.eq.0) then
> 85 call readparam(iconf, conf)
> 86 endif
> 87 call MPI_BCAST(iconf,5,MPI_INTEGER,0,comm,infompi)
> 88 call MPI_BCAST(conf,2,MPI_DOUBLE_PRECISION,0,comm,infompi)
> 89 !
> 90 size_x = iconf(1)
> 91 size_y = iconf(1)
> 92 x_domains = iconf(3)
> 93 y_domains = iconf(4)
> 94 maxStep = iconf(5)
> 95 dt = conf(1)
> 96 epsilon = conf(2)
> 97 !
> 98 size_x_glo = x_domains*size_x+2
> 99 h = 1.0d0/dble(size_x_glo)
> 100 dt = 0.25*h*h
> 101 !
> 102 !
> 103 lda = size_y+2
> 104 allocate(x(0:size_y+1,0:size_x+1))
> 105 allocate(x0(0:size_y+1,0:size_x+1))
> 106 allocate(b(0:size_y+1,0:size_x+1))
> 107 !
> 108 ! Create 2D cartesian grid
> 109 periods(:) = .false.
> 110
> 111 ndims = 2
> 112 dims(1)=x_domains
> 113 dims(2)=y_domains
> 114 CALL MPI_CART_CREATE(MPI_COMM_WORLD, ndims, dims, periods, &
> 115 reorganisation,comm2d,infompi)
> 116 !
> 117 ! Identify neighbors
> 118 !
> 119 NeighBor(:) = MPI_PROC_NULL
> 120 ! Left/West and right/Est neigbors
> 121 CALL MPI_CART_SHIFT(comm2d,0,1,NeighBor(W),NeighBor(E),infompi)
> 122
> 123 print *,'mpi_proc_null=', MPI_PROC_NULL
> 124 print *,'rang=', me
> 125 print *, 'ici premier mpi_cart_shift : neighbor(w)=',NeighBor(W)
> 126 print *, 'ici premier mpi_cart_shift : neighbor(e)=',NeighBor(E)
> 127
> 128 ! Bottom/South and Upper/North neigbors
> 129 CALL MPI_CART_SHIFT(comm2d,1,1,NeighBor(S),NeighBor(N),infompi)
> 130
> 131
> 132 print *, 'ici deuxieme mpi_cart_shift : neighbor(s)=',NeighBor(S)
> 133 print *, 'ici deuxieme mpi_cart_shift : neighbor(n)=',NeighBor(N)
> 134
> 135
> 136
> 137 !
> 138 ! Create row data type to coimmunicate with South and North neighbors
> 139 !
> 140 CALL MPI_TYPE_VECTOR(size_x, 1, size_y+2, MPI_DOUBLE_PRECISION,
> row_type,infompi)
> 141 CALL MPI_TYPE_COMMIT(row_type, infompi)
> 142 !
> 143 ! initialization
> 144 !
> 145 call initvalues(x0, b, size_x+1, size_x )
> 146 !
> 147 ! Update the boundaries
> 148 !
> 149 call updateBound(x0,size_x,size_x, NeighBor, comm2d, row_type)
> 150
> 151 step = 0
> 152 t = 0.0
> 153 !
> 154 tstart = MPI_Wtime()
> 155 ! REPEAT
> 156 10 continue
> 157 !
> 158 step = step + 1
> 159 t = t + dt
> 160 ! perform one step of the explicit scheme
> 161 call Explicit(x0,x,b, size_x+1, size_x, size_x, dt, h, resLoc)
> 162 ! update the partial solution along the interface
> 163 call updateBound(x0,size_x,size_x, NeighBor, comm2d, row_type)
> 164 ! Check the distance between two iterates
> 165 call MPI_ALLREDUCE(resLoc,result,1, MPI_DOUBLE_PRECISION,
> MPI_SUM,comm,infompi)
> 166 result= sqrt(result)
> 167 !
> 168 if (me.eq.0) write(*,1002) t,result
> 169 !
> 170 if ((result.gt.epsilon).and.(step.lt.maxStep)) goto 10
> 171 !
> 172 ! UNTIL "Convergence"
> 173 !
> 174 tend = MPI_Wtime()
> 175 if (me.eq.0) then
> 176 write(*,*)
> 177 write(*,*) ' Convergence after ', step,' steps '
> 178 write(*,*) ' Problem size ',
> size_x*x_domains*size_y*y_domains
> 179 write(*,*) ' Wall Clock ', tend-tstart
> 180
> 181 !
> 182 ! Print the solution at each point of the grid
> 183 !
> 184 write(*,*)
> 185 write(*,*) ' Computed solution '
> 186 write(*,*)
> 187 do 30, j=size_x+1,0,-1
> 188 write(*,1000)(x0(j,i),i=0,size_x+1)
> 189 30 continue
> 190 endif
> 191 !
> 192 call MPI_FINALIZE(infompi)
> 193 !
> 194 deallocate(x)
> 195 deallocate(x0)
> 196 deallocate(b)
> 197 !
> 198 ! Formats available to display the computed values on the grid
> 199 !
> 200 1000 format(100(1x, f7.3))
> 201 1001 format(100(1x, e7.3))
> 202 1002 format(' At time ',E8.2,' Norm ', E8.2)
> 203
> 204 !
> 205 stop
> 206 end
> 207 !
>
> ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
> The routine "Explicit" at line 161 allows to compute the domain 2D ( the
> values of Temperature are stocked in vector "x" : x(1:size_x,1:size_y)) with
> the following scheme :
>
> ! The stencil of the explicit operator for the heat equation
> ! on a regular rectangular grid using a five point finite difference
> ! scheme in space is
> !
> ! | + 1
> |
> ! |
> |
> !dt/(h*h) * | +1 -4 + h*h/dt + 1 |
> ! |
> |
> ! | + 1
> |
> !
> diag = - 4.0 + h*h/dt
> weight = dt/(h*h)
> !
> ! perform an explicit update on the points within the domain
> do 20, j=1,size_x
> do 30, i=1,size_y
> x(i,j) = weight *( x0(i-1,j) + x0(i+1,j) &
> + x0(i,j-1) + x0(i,j+1) +x0(i,j)*diag)
> 30 continue
> 20 continue
>
>
> ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
> The routine "updateBound" updates the bound values like this :
>
> 1 !*******************************************************************
> 2 SUBROUTINE updateBound(x, size_x, size_y, NeighBor, comm2d, row_type)
> 3
> 4 !*****************************************************************
> 5 include 'mpif.h'
> 6 !
> 7 INTEGER size_x, size_y
> 8 !Iterate
> 9 double precision, dimension(0:size_y+1,0:size_x+1) :: x
> 10 !Type row_type
> 11 INTEGER :: row_type
> 12 integer, parameter :: nbvi=4
> 13 integer, parameter :: S=1, E=2, N=3, W=4
> 14 integer, dimension(4) :: neighBor
> 15 !
> 16 INTEGER :: infompi, comm2d
> 17 INTEGER :: flag
> 18 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status
> 19
> 20 !********* North/South communication
> ************************************
> 21 flag = 1
> 22 !Send my boundary to North and receive from South
> 23 CALL MPI_SENDRECV(x(size_y, 1), 1, row_type ,neighBor(N),flag, &
> 24 x(0,1), 1, row_type,neighbor(S),flag, &
> 25 comm2d, status, infompi)
> 26
> 27 !Send my boundary to South and receive from North
> 28 CALL MPI_SENDRECV(x(1,1), 1, row_type ,neighbor(S),flag, &
> 29 x(size_y+1,1), 1, row_type ,neighbor(N),flag, &
> 30 comm2d, status, infompi)
> 31
> 32 !********* Est/West communication ************************************
> 33 flag = 2
> 34 !Send my boundary to West and receive from Est
> 35 CALL MPI_SENDRECV(x(1,1), size_y, MPI_DOUBLE_PRECISION, neighbor(W),
> flag, &
> 36 x(1, size_x+1), size_y,
> MPI_DOUBLE_PRECISION,neighbor(E), flag, &
> 37 comm2d, status, infompi)
> 38
> 39 !Send my boundary to Est and receive from West
> 40 CALL MPI_SENDRECV( x(1,size_x), size_y, MPI_DOUBLE_PRECISION,
> neighbor(E), flag, &
> 41 x(1,0),size_y, MPI_DOUBLE_PRECISION, neighbor(W),
> flag, &
> 42 comm2d, status, infompi)
> 43
> 44 END SUBROUTINE updateBound
>
>
> --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
> I am confused between the shift of the values near to the bounds done by
> the "updateBound" routine and the main loop (at line 161 in main code)
> which calls the routine "Explicit".
> For a given process (say number 1) ( i use 4 here for execution), i send to
> the east process (3) the penultimate column left column, to the north
> process (0) the penultimate row top ,and to the others (mpi_proc_null=-2)
> the penultimate right column and the bottom row. But how the 4 processes
> are synchronous ? I don't understand too why all the processes go through
> the solving piece of code calling the "Explicit" routine.
>
> Here'i the sequential version of this simulation :
>
> 49 program heat
> 50 implicit none
> 51 ! size of the discretization
> 52 integer size_x, maxStep
> 53 integer lda
> 54 double precision, allocatable:: x(:,:),b(:,:), x0(:,:)
> 55 double precision dt, h, result, epsilon
> 56
> 57 ! index loop
> 58 integer i,j, step
> 59 double precision t
> 60 !
> 61 intrinsic abs
> 62 !
> 63 !
> 64 print *,' Size of the square '
> 65 read(*,*) size_x
> 66 h = 1.0/(size_x+1)
> 67 dt = 0.25*h*h
> 68 maxStep = 100
> 69 print *, 'Max. number of steps '
> 70 read(*,*) maxStep
> 71 epsilon = 1.d-8
> 72 !
> 73 allocate(x(0:size_x+1,0:size_x+1))
> 74 allocate(x0(0:size_x+1,0:size_x+1))
> 75 allocate(b(0:size_x+1,0:size_x+1))
> 76 !
> 77 !
> 78 ! initialization
> 79 !
> 80 call initvalues(x0, b, size_x+1, size_x )
> 81 !
> 82 step = 0
> 83 t = 0.0
> 84 ! REPEAT
> 85 10 continue
> 86 !
> 87 step = step + 1
> 88 t = t + dt
> 89 !
> 90 call Explicit(x0, x, b, size_x+1, size_x, size_x, dt, &
> 91 h, result)
> 92 result = sqrt(result)
> 93 if ((result.gt.epsilon).and.(step.lt.maxStep)) goto 10
> 94 !
> 95 write(*,*)
> 96 write(*,*) ' Convergence after ', step,' steps '
> 97 write(*,*) ' Problem size ', size_x*size_x
> 98 !
>
>
> ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
> Sorry for this long post but i would like to clarify my problem as much as
> it was possible.
>
> Any explanation would help me.
>
> Thanks in advance
>
>
