module_boundary.f90
Go to the documentation of this file.
1 !*************************************************************************
2 ! This file is part of HERCULES.
3 !
4 ! HERCULES is free software; you can redistribute it and/or modify
5 ! it under the terms of the GNU General Public License as published by
6 ! the Free Software Foundation; either version 3 of the License, or
7 ! (at your option) any later version.
8 !
9 ! HERCULES is distributed in the hope that it will be useful,
10 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ! GNU General Public License for more details.
13 !
14 ! You should have received a copy of the GNU General Public License
15 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
16 !
17 ! Copyright 2015, Ping He, MEAS, NCSU
18 !*************************************************************************
19 
20 
21 
23 
25 
26  implicit none
27 
28 contains
29 
30 
31 
32  ! update boundary/ghost cell values for u,v,w in the z direction
33  subroutine update_boundary_uvw(u,v,w)
34 
35  implicit none
36 
37  real(mytype),dimension(:,:,:),intent(inout) :: u, v, w
38 
39  integer :: i,j
40 
41  !!!!!!!!!!!!!!!! U Eqn !!!!!!!!!!!!!!!!!
42  ! Top and bottom ghost cells for u
43  ! check if it is open channel or closed channel
44  if (ochannel .eq. 0) then
45 
46  do j=jstr3,jend3,1
47  do i=istr3,iend3,1
48 
49  ! if moving reference frame is used, we need to
50  ! give the non-slip walls with a moving speed
51  ! of -u_mrf in the x direction
52  u(i,j,kend3+1)=-u(i,j,kend3)-2.d0*u_mrf ! non-slip
53  u(i,j,kstr3-1)=-u(i,j,kstr3)-2.d0*u_mrf ! non-slip
54 
55  enddo
56  enddo
57 
58  elseif (ochannel .eq. 1) then
59 
60  do j=jstr3,jend3,1
61  do i=istr3,iend3,1
62 
63  ! if moving reference frame is used, we need to
64  ! give the non-slip walls with a moving speed
65  ! of -u_mrf in the x direction
66  u(i,j,kend3+1)= u(i,j,kend3) ! slip
67  u(i,j,kstr3-1)=-u(i,j,kstr3)-2.d0*u_mrf ! non-slip
68 
69  enddo
70  enddo
71 
72  endif
73 
74  !!!!!!!!!!!!!!! V Eqn !!!!!!!!!!!!!!!!!
75  ! Top and bottom ghost cells for v
76  ! check if it is open channel or closed channel
77  if (ochannel .eq. 0) then
78 
79  do j=jstr3,jend3,1
80  do i=istr3,iend3,1
81 
82  v(i,j,kend3+1)=-v(i,j,kend3) ! non-slip
83  v(i,j,kstr3-1)=-v(i,j,kstr3) ! non-slip
84 
85  enddo
86  enddo
87 
88  elseif (ochannel .eq. 1) then
89 
90  do j=jstr3,jend3,1
91  do i=istr3,iend3,1
92 
93  v(i,j,kend3+1)= v(i,j,kend3) ! slip
94  v(i,j,kstr3-1)=-v(i,j,kstr3) ! non-slip
95 
96  enddo
97  enddo
98 
99  endif
100 
101  !!!!!!!!!!!! W Equation !!!!!!!!!!!!!!
102  ! Top and bottom boundary for w
103  ! no matter it is open or closed channel, w=0 at the top and bottom
104  do j=jstr3,jend3,1
105  do i=istr3,iend3,1
106 
107  w(i,j,kend3)=0.d0
108  w(i,j,kend3+1)=-w(i,j,kend3-1)
109  w(i,j,kstr3-1)=0.d0
110 
111  enddo
112  enddo
113 
114  ! communication between processors
115  ! for horizontal ghost cells
116  ! we don't need this for spectral method
117  if (cds .ne. 3) then
118 
119  call comm_bound(u)
120  call comm_bound(v)
121  call comm_bound(w)
122 
123  endif
124 
125 
126  return
127 
128  end subroutine
129 
130 
131 
132  ! update boundary/ghost cell values for p
133  subroutine update_boundary_p ( p )
135  implicit none
136 
137  real(mytype),dimension(:,:,:),intent(inout) :: p
138 
139  integer :: i,j
140 
141  ! Top and bottom ghost cells for p, zero-gradient
142  do j=jstr3,jend3,1
143  do i=istr3,iend3,1
144 
145  p(i,j,kend3+1)=p(i,j,kend3)
146  p(i,j,kstr3-1)=p(i,j,kstr3)
147 
148  enddo
149  enddo
150 
151  ! communication between processors
152  ! for horizontal ghost cells
153  ! we don't need this for spectral method
154  if (cds .ne. 3) then
155 
156  call comm_bound(p)
157 
158  endif
159 
160  return
161 
162  end subroutine
163 
164 
165 
166  ! update boundary/ghost cell value for t
167  subroutine update_boundary_t ( t )
169  implicit none
170 
171  real(mytype),dimension(:,:,:),intent(inout) :: t
172 
173  integer :: i,j
174 
175  ! Top and bottom ghost cells for t, fixed temperature
176  do j=jstr3,jend3,1
177  do i=istr3,iend3,1
178 
179  t(i,j,kend3+1)=2.d0*ttop-t(i,j,kend3)
180  t(i,j,kstr3-1)=2.d0*tbot-t(i,j,kstr3)
181 
182  enddo
183  enddo
184 
185  ! communication between processors
186  ! for horizontal ghost cells
187  ! we don't need this for spectral method
188  if (cds .ne. 3) then
189 
190  call comm_bound(t)
191 
192  endif
193 
194  return
195 
196  end subroutine
197 
198 
199 
200  ! update boundary/ghost cell values for interpolated variables uhx, vhx, etc.
201  ! this is only needed for FD method
202  subroutine update_boundary_uvwhxyz(u,v,uhx,uhy,uhz,vhx,vhy,vhz,whx,why)
204  implicit none
205 
206  real(mytype),dimension(:,:,:), intent(in):: u, v
207  real(mytype),dimension(:,:,:),intent(inout):: uhx,uhy,uhz, &
208  vhx,vhy,vhz,whx,why
209 
210  integer :: i,j
211 
212  ! check if it is open channel or closed channel
213  if (ochannel .eq. 0) then
214 
215  ! Top boundary, non-slip
216  do j=jstr3,jend3,1
217  do i=istr3,iend3,1
218 
219  ! if moving reference frame is used, we need to
220  ! give the non-slip walls with a moving speed
221  ! of -u_mrf in the x direction
222  uhx(i,j,kend3+1)=-uhx(i,j,kend3)-2.d0*u_mrf
223  uhy(i,j,kend3+1)=-uhy(i,j,kend3)-2.d0*u_mrf
224  uhz(i,j,kend3) =-u_mrf
225  uhz(i,j,kend3+1)=-uhz(i,j,kend3-1)-2.d0*u_mrf
226 
227  vhx(i,j,kend3+1)=-vhx(i,j,kend3)
228  vhy(i,j,kend3+1)=-vhy(i,j,kend3)
229  vhz(i,j,kend3) = 0.d0
230  vhz(i,j,kend3+1)=-vhz(i,j,kend3-1)
231 
232  whx(i,j,kend3) = 0.d0
233  whx(i,j,kend3+1)=-whx(i,j,kend3-1)
234  why(i,j,kend3) = 0.d0
235  why(i,j,kend3+1)=-why(i,j,kend3-1)
236 
237  enddo
238  enddo
239 
240  elseif (ochannel .eq. 1) then
241 
242  ! open channel, slip for the top and non-slip for the bottom
243  do j=jstr3,jend3,1
244  do i=istr3,iend3,1
245 
246  uhx(i,j,kend3+1)= uhx(i,j,kend3)
247  uhy(i,j,kend3+1)= uhy(i,j,kend3)
248  uhz(i,j,kend3) = u(i,j,kend3)
249  uhz(i,j,kend3+1)= uhz(i,j,kend3)
250 
251  vhx(i,j,kend3+1)= vhx(i,j,kend3)
252  vhy(i,j,kend3+1)= vhy(i,j,kend3)
253  vhz(i,j,kend3) = v(i,j,kend3)
254  vhz(i,j,kend3+1)= vhz(i,j,kend3)
255 
256  whx(i,j,kend3) = 0.d0
257  whx(i,j,kend3+1)=-whx(i,j,kend3-1)
258  why(i,j,kend3) = 0.d0
259  why(i,j,kend3+1)=-why(i,j,kend3-1)
260 
261  enddo
262  enddo
263 
264  endif
265 
266  ! Bottom boundary, non-slip
267  do j=jstr3,jend3,1
268  do i=istr3,iend3,1
269 
270  ! if moving reference frame is used, we need to
271  ! give the non-slip walls with a moving speed
272  ! of -u_mrf in the x direction
273  uhx(i,j,kstr3-1)=-uhx(i,j,kstr3)-2.d0*u_mrf
274  uhy(i,j,kstr3-1)=-uhy(i,j,kstr3)-2.d0*u_mrf
275  uhz(i,j,kstr3-1)=-u_mrf
276 
277  vhx(i,j,kstr3-1)=-vhx(i,j,kstr3)
278  vhy(i,j,kstr3-1)=-vhy(i,j,kstr3)
279  vhz(i,j,kstr3-1)=0.d0
280 
281  whx(i,j,kstr3-1)=0.d0
282  why(i,j,kstr3-1)=0.d0
283 
284  enddo
285  enddo
286 
287  ! communication between processors
288  ! for horizontal ghost cells
289  ! we don't need this for spectral method
290  if (cds .ne. 3) then
291 
292  call comm_bound(uhx)
293  call comm_bound(uhy)
294  call comm_bound(uhz)
295  call comm_bound(vhx)
296  call comm_bound(vhy)
297  call comm_bound(vhz)
298  call comm_bound(whx)
299  call comm_bound(why)
300 
301  endif
302 
303  return
304 
305  end subroutine
306 
307 
308 
309  ! update boundary/ghost cell values for interpolated variables thx, thx, etc.
310  ! this is only needed for FD method
311  subroutine update_boundary_thxyz (thx,thy,thz)
313  implicit none
314 
315  real(mytype),dimension(:,:,:),intent(inout):: thx,thy,thz
316 
317  integer :: i,j
318 
319  ! Top ghost cells for t, fixed temperature
320  do j=jstr3,jend3,1
321  do i=istr3,iend3,1
322 
323  thx(i,j,kend3+1)=-thx(i,j,kend3)+2.d0*ttop
324  thy(i,j,kend3+1)=-thy(i,j,kend3)+2.d0*ttop
325  thz(i,j,kend3) = ttop
326  thz(i,j,kend3+1)=-thz(i,j,kend3-1)+2.d0*ttop
327 
328  enddo
329  enddo
330 
331  ! Bottom ghost cells for t, fixed temperature
332  do j=jstr3,jend3,1
333  do i=istr3,iend3,1
334 
335  thx(i,j,kstr3-1)=-thx(i,j,kstr3)+2.d0*tbot
336  thy(i,j,kstr3-1)=-thy(i,j,kstr3)+2.d0*tbot
337  thz(i,j,kstr3-1)=tbot
338 
339  enddo
340  enddo
341 
342  ! communication between processors
343  ! for horizontal ghost cells
344  ! we don't need this for spectral method
345  if (cds .ne. 3) then
346 
347  call comm_bound(thx)
348  call comm_bound(thy)
349  call comm_bound(thz)
350 
351  endif
352 
353  return
354 
355  end subroutine
356 
357 
358 
359  ! this function does the communication between processors for the
360  ! East and West boundary ghost cells.
361  subroutine comm_bound(var)
363  use mpi
364 
365  implicit none
366 
367  real(mytype),dimension(:,:,:),intent(inout) :: var
368  ! idsend: the processor id to send data to
369  ! isrecv: the processor id to receive data from
370  ! count_ns: how many cells to send/receive in the north and south directions
371  ! count_ew: how many cells to send/receive in the east and west directions
372  integer :: idsend,idrecv,count_ns,count_ew
373  integer,dimension(MPI_STATUS_SIZE):: status1
374 
375  ! calculate how many grids to exchange, note that for north and south,
376  ! we need to extend two ghost cells in the x direction
377  count_ew=(ny/p_col)*(nz+2)
378  count_ns=(nx/p_row+2*ghst)*(nz+2)
379 
380  ! DATA exchange
381  ! first step: send data to EAST and receive data from WEST
382  ! calculate send and receive id
383  ! here myid_rowindx, pij is calculated in initiate_indices
384  ! if it is the west boundary sub-domain, send data to
385  ! the east boundary sub-domain
386  ! otherwise, send data to the neighbor sub-domain to the west
387  if (myid_rowindx .eq. p_row) then
388  idsend=pij(1,myid_colindx)
389  else
390  idsend=myid+p_col
391  endif
392  if (myid_rowindx .eq. 1) then
393  idrecv=pij(p_row,myid_colindx)
394  else
395  idrecv=myid-p_col
396  endif
397  ! do communication for the 1st ghost cells
398  call mpi_sendrecv(var(iend3, jstr3:jend3,kstr3-1:kend3+1), &
399  count_ew,real_type,idsend,99, &
400  var(istr3-1,jstr3:jend3,kstr3-1:kend3+1), &
401  count_ew,real_type,idrecv,99, &
402  mpi_comm_world,status1,ierr)
403  ! do communication for the 2st ghost cells
404  call mpi_sendrecv(var(iend3-1,jstr3:jend3,kstr3-1:kend3+1), &
405  count_ew,real_type,idsend,98, &
406  var(istr3-2,jstr3:jend3,kstr3-1:kend3+1), &
407  count_ew,real_type,idrecv,98, &
408  mpi_comm_world,status1,ierr)
409 
410  ! second step: send data to WEST and receive data from EAST
411  ! calculate send and receive id
412  if (myid_rowindx .eq. 1) then
413  idsend=pij(p_row,myid_colindx)
414  else
415  idsend=myid-p_col
416  endif
417  if (myid_rowindx .eq. p_row) then
418  idrecv=pij(1,myid_colindx)
419  else
420  idrecv=myid+p_col
421  endif
422  ! do communication for the 1st ghost cells
423  call mpi_sendrecv(var(istr3, jstr3:jend3,kstr3-1:kend3+1), &
424  count_ew,real_type,idsend,97, &
425  var(iend3+1,jstr3:jend3,kstr3-1:kend3+1), &
426  count_ew,real_type,idrecv,97, &
427  mpi_comm_world,status1,ierr)
428  ! do communication for the 2st ghost cells
429  call mpi_sendrecv(var(istr3+1,jstr3:jend3,kstr3-1:kend3+1), &
430  count_ew,real_type,idsend,96, &
431  var(iend3+2,jstr3:jend3,kstr3-1:kend3+1), &
432  count_ew,real_type,idrecv,96, &
433  mpi_comm_world,status1,ierr)
434 
435  ! third step: send data to NORTH and receive data from SOUTH
436  ! calculate send and receive id
437  if (myid_colindx .eq. p_col) then
438  idsend=pij(myid_rowindx,1)
439  else
440  idsend=myid+1
441  endif
442  if (myid_colindx .eq. 1) then
443  idrecv=pij(myid_rowindx,p_col)
444  else
445  idrecv=myid-1
446  endif
447  ! do communication for the 1st ghost cells
448  call mpi_sendrecv(var(istr3-ghst:iend3+ghst,jend3, kstr3-1:kend3+1), &
449  count_ns,real_type,idsend,95, &
450  var(istr3-ghst:iend3+ghst,jstr3-1,kstr3-1:kend3+1), &
451  count_ns,real_type,idrecv,95, &
452  mpi_comm_world,status1,ierr)
453  ! do communication for the 2st ghost cells
454  call mpi_sendrecv(var(istr3-ghst:iend3+ghst,jend3-1,kstr3-1:kend3+1), &
455  count_ns,real_type,idsend,94, &
456  var(istr3-ghst:iend3+ghst,jstr3-2,kstr3-1:kend3+1), &
457  count_ns,real_type,idrecv,94, &
458  mpi_comm_world,status1,ierr)
459 
460  ! fourth step: send data to SOUTH and receive data from NORTH
461  ! calculate send and receive id
462  if (myid_colindx .eq. 1) then
463  idsend=pij(myid_rowindx,p_col)
464  else
465  idsend=myid-1
466  endif
467  if (myid_colindx .eq. p_col) then
468  idrecv=pij(myid_rowindx,1)
469  else
470  idrecv=myid+1
471  endif
472  ! do communication for the 1st ghost cells
473  call mpi_sendrecv(var(istr3-ghst:iend3+ghst,jstr3, kstr3-1:kend3+1), &
474  count_ns,real_type,idsend,93, &
475  var(istr3-ghst:iend3+ghst,jend3+1,kstr3-1:kend3+1), &
476  count_ns,real_type,idrecv,93, &
477  mpi_comm_world,status1,ierr)
478  ! do communication for the 2st ghost cells
479  call mpi_sendrecv(var(istr3-ghst:iend3+ghst,jstr3+1,kstr3-1:kend3+1), &
480  count_ns,real_type,idsend,92, &
481  var(istr3-ghst:iend3+ghst,jend3+2,kstr3-1:kend3+1), &
482  count_ns,real_type,idrecv,92, &
483  mpi_comm_world,status1,ierr)
484 
485  return
486 
487  end subroutine
488 
489 
490 
491 end module
subroutine update_boundary_p(p)
integer, save, protected myid
integer, save, protected ochannel
integer, save, protected ny
integer, dimension(:,:), allocatable, save, protected pij
integer, save, protected jend3
real(mytype), save, protected ttop
integer, save, protected nz
integer, save, protected p_row
subroutine update_boundary_uvwhxyz(u, v, uhx, uhy, uhz, vhx, vhy, vhz, whx, why)
integer, save, protected myid_rowindx
subroutine update_boundary_thxyz(thx, thy, thz)
integer, save, protected istr3
subroutine comm_bound(var)
integer, save, protected jstr3
real(mytype), save, protected u_mrf
integer, save, protected p_col
integer, parameter ghst
integer, save, protected cds
subroutine update_boundary_uvw(u, v, w)
integer, save, protected nx
subroutine update_boundary_t(t)
integer, save, protected iend3
real(mytype), save, protected tbot
integer, save, protected kstr3
integer, save, protected myid_colindx
integer, save, protected kend3