FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
fstr_contact_comm.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
7 
8  use hecmw
9 
10  implicit none
11 
12  private
13  public :: fstrst_contact_comm
14  public :: fstr_contact_comm_init
22 
23  type fstrst_contact_comm
24  private
25  integer(kind=kint) :: n_neighbor_pe
26  integer(kind=kint), pointer :: neighbor_pe(:)
27  integer(kind=kint) :: MPI_COMM
28  integer(kind=kint), pointer :: ext_index(:)
29  integer(kind=kint), pointer :: ext_item(:)
30  integer(kind=kint), pointer :: int_index(:)
31  integer(kind=kint), pointer :: int_item(:)
32  end type fstrst_contact_comm
33 
34  integer(kind=kint), parameter :: op_overwrite = 46810
35 
36  integer(kind=kint), parameter :: DEBUG = 0
37 
38 contains
39 
40  subroutine fstr_contact_comm_init(conComm, hecMESH, ndof, n_contact_dof, contact_dofs)
41  implicit none
42  type (fstrst_contact_comm), intent(inout) :: concomm
43  type (hecmwst_local_mesh), intent(in) :: hecmesh
44  integer(kind=kint), intent(in) :: ndof, n_contact_dof
45  integer(kind=kint), intent(in) :: contact_dofs(:)
46  integer(kind=kint), pointer :: ext_index(:), ext_item(:), int_index(:), int_item(:)
47  integer(kind=kint), allocatable :: n_ext_per_dom(:), n_int_per_dom(:), ext_item_remote(:)
48  integer(kind=kint), allocatable :: statuses(:,:), requests(:)
49  integer(kind=kint) :: nn_int, np, ilag, icontact, irow, irank, idom, tag, idof, idx, irow_remote
50  integer(kind=kint) :: n_send, is, ie, len
51  if (hecmesh%n_neighbor_pe == 0) then
52  concomm%n_neighbor_pe = 0
53  return
54  endif
55  nn_int = hecmesh%nn_internal
56  np = hecmesh%n_node
57  ! count external contact dofs
58  allocate(n_ext_per_dom(hecmesh%n_neighbor_pe))
59  n_ext_per_dom(:) = 0
60  do ilag = 1, n_contact_dof
61  icontact = contact_dofs(ilag)
62  if (icontact <= nn_int*ndof) cycle ! skip internal contact dof
63  irow = (icontact+ndof-1) / ndof
64  irank = hecmesh%node_ID(2*irow)
65  call rank_to_idom(hecmesh, irank, idom)
66  n_ext_per_dom(idom) = n_ext_per_dom(idom) + 1
67  enddo
68  ! send external / recv internal contact dofs
69  allocate(statuses(hecmw_status_size, hecmesh%n_neighbor_pe))
70  allocate(requests(hecmesh%n_neighbor_pe))
71  do idom = 1, hecmesh%n_neighbor_pe
72  irank = hecmesh%neighbor_pe(idom)
73  tag = 5001
74  call hecmw_isend_int(n_ext_per_dom(idom), 1, irank, tag, &
75  hecmesh%MPI_COMM, requests(idom))
76  enddo
77  allocate(n_int_per_dom(hecmesh%n_neighbor_pe))
78  do idom = 1, hecmesh%n_neighbor_pe
79  irank = hecmesh%neighbor_pe(idom)
80  tag = 5001
81  call hecmw_recv_int(n_int_per_dom(idom), 1, irank, tag, &
82  hecmesh%MPI_COMM, statuses(:,1))
83  enddo
84  call hecmw_waitall(hecmesh%n_neighbor_pe, requests, statuses)
85  ! make index
86  allocate(ext_index(0:hecmesh%n_neighbor_pe))
87  allocate(int_index(0:hecmesh%n_neighbor_pe))
88  ext_index(0) = 0
89  int_index(0) = 0
90  do idom = 1, hecmesh%n_neighbor_pe
91  ext_index(idom) = ext_index(idom-1) + n_ext_per_dom(idom)
92  int_index(idom) = int_index(idom-1) + n_int_per_dom(idom)
93  enddo
94  ! make ext_item
95  allocate(ext_item(ext_index(hecmesh%n_neighbor_pe)))
96  allocate(ext_item_remote(ext_index(hecmesh%n_neighbor_pe)))
97  n_ext_per_dom(:) = 0
98  do ilag = 1, n_contact_dof
99  icontact = contact_dofs(ilag)
100  if (icontact <= nn_int*ndof) cycle ! skip internal contact dof
101  irow = (icontact+ndof-1) / ndof
102  idof = icontact - ndof*(irow-1)
103  irank = hecmesh%node_ID(2*irow)
104  call rank_to_idom(hecmesh, irank, idom)
105  n_ext_per_dom(idom) = n_ext_per_dom(idom) + 1
106  idx = ext_index(idom-1)+n_ext_per_dom(idom)
107  ext_item(idx) = icontact
108  irow_remote = hecmesh%node_ID(2*irow-1)
109  ext_item_remote(idx) = ndof*(irow_remote-1)+idof
110  enddo
111  deallocate(n_ext_per_dom)
112  deallocate(n_int_per_dom)
113  ! send ext_item_remote and recv int_item
114  n_send = 0
115  do idom = 1, hecmesh%n_neighbor_pe
116  irank = hecmesh%neighbor_pe(idom)
117  is = ext_index(idom-1)+1
118  ie = ext_index(idom)
119  len = ie-is+1
120  if (len == 0) cycle
121  n_send = n_send + 1
122  tag = 5002
123  call hecmw_isend_int(ext_item_remote(is:ie), len, irank, tag, &
124  hecmesh%MPI_COMM, requests(n_send))
125  enddo
126  allocate(int_item(int_index(hecmesh%n_neighbor_pe)))
127  do idom = 1, hecmesh%n_neighbor_pe
128  irank = hecmesh%neighbor_pe(idom)
129  is = int_index(idom-1)+1
130  ie = int_index(idom)
131  len = ie-is+1
132  if (len == 0) cycle
133  tag = 5002
134  call hecmw_recv_int(int_item(is:ie), len, irank, tag, &
135  hecmesh%MPI_COMM, statuses(:,1))
136  enddo
137  call hecmw_waitall(n_send, requests, statuses)
138  deallocate(statuses, requests)
139  if (debug >= 2) then
140  write(0,*) ' DEBUG2: ext_index',ext_index(:)
141  write(0,*) ' DEBUG2: ext_item',ext_item(:)
142  write(0,*) ' DEBUG2: ext_item_remote',ext_item_remote(:)
143  write(0,*) ' DEBUG2: int_index',int_index(:)
144  write(0,*) ' DEBUG2: int_item',int_item(:)
145  endif
146  deallocate(ext_item_remote)
147  !
148  concomm%n_neighbor_pe = hecmesh%n_neighbor_pe
149  allocate(concomm%neighbor_pe(concomm%n_neighbor_pe))
150  concomm%neighbor_pe(:) = hecmesh%neighbor_pe(:)
151  concomm%MPI_COMM = hecmesh%MPI_COMM
152  concomm%ext_index => ext_index
153  concomm%ext_item => ext_item
154  concomm%int_index => int_index
155  concomm%int_item => int_item
156  end subroutine fstr_contact_comm_init
157 
158  subroutine fstr_contact_comm_finalize(conComm)
159  implicit none
160  type (fstrst_contact_comm), intent(inout) :: concomm
161  if (concomm%n_neighbor_pe == 0) return
162  if (associated(concomm%neighbor_pe)) deallocate(concomm%neighbor_pe)
163  if (associated(concomm%ext_index)) deallocate(concomm%ext_index)
164  if (associated(concomm%ext_item)) deallocate(concomm%ext_item)
165  if (associated(concomm%int_index)) deallocate(concomm%int_index)
166  if (associated(concomm%int_item)) deallocate(concomm%int_item)
167  concomm%n_neighbor_pe = 0
168  concomm%MPI_COMM = 0
169  end subroutine fstr_contact_comm_finalize
170 
171  subroutine fstr_contact_comm_reduce_r(conComm, vec, op)
172  implicit none
173  type (fstrst_contact_comm), intent(in) :: concomm
174  real(kind=kreal), intent(inout) :: vec(:)
175  integer(kind=kint), intent(in) :: op
176  if (concomm%n_neighbor_pe == 0) return
177  call send_recv_contact_info_r(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
178  concomm%ext_index, concomm%ext_item, concomm%int_index, concomm%int_item, vec, op)
179  end subroutine fstr_contact_comm_reduce_r
180 
181  subroutine fstr_contact_comm_bcast_r(conComm, vec)
182  implicit none
183  type (fstrst_contact_comm), intent(in) :: concomm
184  real(kind=kreal), intent(inout) :: vec(:)
185  integer(kind=kint) :: op
186  if (concomm%n_neighbor_pe == 0) return
187  op = op_overwrite
188  call send_recv_contact_info_r(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
189  concomm%int_index, concomm%int_item, concomm%ext_index, concomm%ext_item, vec, op)
190  end subroutine fstr_contact_comm_bcast_r
191 
192  subroutine fstr_contact_comm_reduce_i(conComm, vec, op)
193  implicit none
194  type (fstrst_contact_comm), intent(in) :: concomm
195  integer(kind=kint), intent(inout) :: vec(:)
196  integer(kind=kint), intent(in) :: op
197  if (concomm%n_neighbor_pe == 0) return
198  call send_recv_contact_info_i(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
199  concomm%ext_index, concomm%ext_item, concomm%int_index, concomm%int_item, vec, op)
200  end subroutine fstr_contact_comm_reduce_i
201 
202  subroutine fstr_contact_comm_bcast_i(conComm, vec)
203  implicit none
204  type (fstrst_contact_comm), intent(in) :: concomm
205  integer(kind=kint), intent(inout) :: vec(:)
206  integer(kind=kint) :: op
207  if (concomm%n_neighbor_pe == 0) return
208  op = op_overwrite
209  call send_recv_contact_info_i(concomm%n_neighbor_pe, concomm%neighbor_pe, concomm%MPI_COMM, &
210  concomm%int_index, concomm%int_item, concomm%ext_index, concomm%ext_item, vec, op)
211  end subroutine fstr_contact_comm_bcast_i
212 
213  subroutine fstr_contact_comm_allreduce_r(conComm, vec, op)
214  implicit none
215  type (fstrst_contact_comm), intent(in) :: concomm
216  real(kind=kreal), intent(inout) :: vec(:)
217  integer(kind=kint), intent(in) :: op
218  call fstr_contact_comm_reduce_r(concomm, vec, op)
219  call fstr_contact_comm_bcast_r(concomm, vec)
220  end subroutine fstr_contact_comm_allreduce_r
221 
222  subroutine fstr_contact_comm_allreduce_i(conComm, vec, op)
223  implicit none
224  type (fstrst_contact_comm), intent(in) :: concomm
225  integer(kind=kint), intent(inout) :: vec(:)
226  integer(kind=kint), intent(in) :: op
227  call fstr_contact_comm_reduce_i(concomm, vec, op)
228  call fstr_contact_comm_bcast_i(concomm, vec)
229  end subroutine fstr_contact_comm_allreduce_i
230 
231  !
232  ! private subroutines
233  !
234 
235  subroutine rank_to_idom(hecMESH, rank, idom)
236  implicit none
237  type (hecmwst_local_mesh), intent(in) :: hecmesh
238  integer(kind=kint), intent(in) :: rank
239  integer(kind=kint), intent(out) :: idom
240  integer(kind=kint) :: i
241  do i = 1, hecmesh%n_neighbor_pe
242  if (hecmesh%neighbor_pe(i) == rank) then
243  idom = i
244  return
245  endif
246  enddo
247  stop 'ERROR: exp_rank not found in neighbor_pe'
248  end subroutine rank_to_idom
249 
250  subroutine send_recv_contact_info_r(n_neighbor_pe, neighbor_pe, MPI_COMM, &
251  send_index, send_item, recv_index, recv_item, vec, op)
252  implicit none
253  integer(kind=kint), intent(in) :: n_neighbor_pe
254  integer(kind=kint), intent(in) :: neighbor_pe(:)
255  integer(kind=kint), intent(in) :: mpi_comm
256  integer(kind=kint), pointer, intent(in) :: send_index(:), send_item(:), recv_index(:), recv_item(:)
257  real(kind=kreal), intent(inout) :: vec(:)
258  integer(kind=kint), intent(in) :: op
259  real(kind=kreal), allocatable :: send_buf(:), recv_buf(:)
260  integer(kind=kint) :: i, n_send, idom, irank, is, ie, len, tag
261  integer(kind=kint), allocatable :: requests(:), statuses(:,:)
262  if (n_neighbor_pe == 0) return
263  allocate(requests(n_neighbor_pe))
264  allocate(statuses(hecmw_status_size, n_neighbor_pe))
265  allocate(send_buf(send_index(n_neighbor_pe)))
266  allocate(recv_buf(recv_index(n_neighbor_pe)))
267  do i = 1, send_index(n_neighbor_pe)
268  send_buf(i) = vec(send_item(i))
269  enddo
270  n_send = 0
271  do idom = 1, n_neighbor_pe
272  irank = neighbor_pe(idom)
273  is = send_index(idom-1)+1
274  ie = send_index(idom)
275  len = ie-is+1
276  if (len == 0) cycle
277  n_send = n_send + 1
278  tag = 5011
279  call hecmw_isend_r(send_buf(is:ie), len, irank, tag, &
280  mpi_comm, requests(n_send))
281  enddo
282  do idom = 1, n_neighbor_pe
283  irank = neighbor_pe(idom)
284  is = recv_index(idom-1)+1
285  ie = recv_index(idom)
286  len = ie-is+1
287  if (len == 0) cycle
288  tag = 5011
289  call hecmw_recv_r(recv_buf(is:ie), len, irank, tag, &
290  mpi_comm, statuses(:,1))
291  enddo
292  call hecmw_waitall(n_send, requests, statuses)
293  if (op == hecmw_sum) then
294  do i = 1, recv_index(n_neighbor_pe)
295  vec(recv_item(i)) = vec(recv_item(i)) + recv_buf(i)
296  enddo
297  elseif (op == hecmw_prod) then
298  do i = 1, recv_index(n_neighbor_pe)
299  vec(recv_item(i)) = vec(recv_item(i)) * recv_buf(i)
300  enddo
301  elseif (op == hecmw_max) then
302  do i = 1, recv_index(n_neighbor_pe)
303  vec(recv_item(i)) = max(vec(recv_item(i)), recv_buf(i))
304  enddo
305  elseif (op == hecmw_min) then
306  do i = 1, recv_index(n_neighbor_pe)
307  vec(recv_item(i)) = min(vec(recv_item(i)), recv_buf(i))
308  enddo
309  else ! overwrite
310  do i = 1, recv_index(n_neighbor_pe)
311  vec(recv_item(i)) = recv_buf(i)
312  enddo
313  endif
314  deallocate(requests)
315  deallocate(statuses)
316  if (debug >= 2) then
317  write(0,*) ' DEBUG2: send_buf',send_buf(:)
318  write(0,*) ' DEBUG2: recv_buf',recv_buf(:)
319  endif
320  deallocate(send_buf)
321  deallocate(recv_buf)
322  end subroutine send_recv_contact_info_r
323 
324  subroutine send_recv_contact_info_i(n_neighbor_pe, neighbor_pe, MPI_COMM, &
325  send_index, send_item, recv_index, recv_item, vec, op)
326  implicit none
327  integer(kind=kint), intent(in) :: n_neighbor_pe
328  integer(kind=kint), intent(in) :: neighbor_pe(:)
329  integer(kind=kint), intent(in) :: mpi_comm
330  integer(kind=kint), pointer, intent(in) :: send_index(:), send_item(:), recv_index(:), recv_item(:)
331  integer(kind=kint), intent(inout) :: vec(:)
332  integer(kind=kint), intent(in) :: op
333  integer(kind=kint), allocatable :: send_buf(:), recv_buf(:)
334  integer(kind=kint) :: i, n_send, idom, irank, is, ie, len, tag
335  integer(kind=kint), allocatable :: requests(:), statuses(:,:)
336  if (n_neighbor_pe == 0) return
337  allocate(requests(n_neighbor_pe))
338  allocate(statuses(hecmw_status_size, n_neighbor_pe))
339  allocate(send_buf(send_index(n_neighbor_pe)))
340  allocate(recv_buf(recv_index(n_neighbor_pe)))
341  do i = 1, send_index(n_neighbor_pe)
342  send_buf(i) = vec(send_item(i))
343  enddo
344  n_send = 0
345  do idom = 1, n_neighbor_pe
346  irank = neighbor_pe(idom)
347  is = send_index(idom-1)+1
348  ie = send_index(idom)
349  len = ie-is+1
350  if (len == 0) cycle
351  n_send = n_send + 1
352  tag = 5011
353  call hecmw_isend_int(send_buf(is:ie), len, irank, tag, &
354  mpi_comm, requests(n_send))
355  enddo
356  do idom = 1, n_neighbor_pe
357  irank = neighbor_pe(idom)
358  is = recv_index(idom-1)+1
359  ie = recv_index(idom)
360  len = ie-is+1
361  if (len == 0) cycle
362  tag = 5011
363  call hecmw_recv_int(recv_buf(is:ie), len, irank, tag, &
364  mpi_comm, statuses(:,1))
365  enddo
366  call hecmw_waitall(n_send, requests, statuses)
367  if (op == hecmw_sum) then
368  do i = 1, recv_index(n_neighbor_pe)
369  vec(recv_item(i)) = vec(recv_item(i)) + recv_buf(i)
370  enddo
371  elseif (op == hecmw_prod) then
372  do i = 1, recv_index(n_neighbor_pe)
373  vec(recv_item(i)) = vec(recv_item(i)) * recv_buf(i)
374  enddo
375  elseif (op == hecmw_max) then
376  do i = 1, recv_index(n_neighbor_pe)
377  vec(recv_item(i)) = max(vec(recv_item(i)), recv_buf(i))
378  enddo
379  elseif (op == hecmw_min) then
380  do i = 1, recv_index(n_neighbor_pe)
381  vec(recv_item(i)) = min(vec(recv_item(i)), recv_buf(i))
382  enddo
383  else ! overwrite
384  do i = 1, recv_index(n_neighbor_pe)
385  vec(recv_item(i)) = recv_buf(i)
386  enddo
387  endif
388  deallocate(requests)
389  deallocate(statuses)
390  if (debug >= 2) then
391  write(0,*) ' DEBUG2: send_buf',send_buf(:)
392  write(0,*) ' DEBUG2: recv_buf',recv_buf(:)
393  endif
394  deallocate(send_buf)
395  deallocate(recv_buf)
396  end subroutine send_recv_contact_info_i
397 
398 end module m_fstr_contact_comm
Definition: hecmw.f90:6
subroutine, public fstr_contact_comm_init(conComm, hecMESH, ndof, n_contact_dof, contact_dofs)
subroutine, public fstr_contact_comm_bcast_i(conComm, vec)
subroutine, public fstr_contact_comm_finalize(conComm)
subroutine, public fstr_contact_comm_allreduce_i(conComm, vec, op)
subroutine, public fstr_contact_comm_reduce_r(conComm, vec, op)
subroutine, public fstr_contact_comm_reduce_i(conComm, vec, op)
subroutine, public fstr_contact_comm_allreduce_r(conComm, vec, op)
subroutine, public fstr_contact_comm_bcast_r(conComm, vec)