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
