DISPATCH
data_hub_mod.f90
1 !===============================================================================
2 !> A communications hub between patches. Three methods are used: download_same,
3 !> download_higher and download_lower
4 !===============================================================================
6  USE io_mod
7  USE omp_mod
8  USE bits_mod
9  USE kinds_mod
10  USE mesh_mod
11  USE trace_mod
12  USE task_mod
13  USE patch_mod
14  USE link_mod
15  implicit none
16  private
17  ! --------------------------------------------------------------------
18  type, public:: data_hub_t
19  class(patch_t), pointer:: target => null() ! pointer to a target patch
20  class(patch_t), pointer:: source => null() ! pointer to a source patch
21  class(data_hub_t), pointer :: head => null()
22  class(data_hub_t), pointer :: next => null()
23  class(data_hub_t), pointer :: prev => null()
24  integer, dimension(:,:), pointer :: t_l => null(), t_u => null()
25  integer, dimension(:,:), pointer :: s_l => null(), s_u => null()
26  logical :: first_time = .true. ! if first time, hib has to be initialized.
27  logical :: rt_target = .false.
28  logical :: all_cells = .false.
29  procedure(dnload_from_same), pointer :: dnload => null()
30  integer :: int_order = -1
31  integer :: level
32  integer :: id
33  contains
34  procedure:: add_hub
35  procedure:: check_in_nbors
36  procedure:: check_in_target
37  procedure:: dnload_from_higher
38  procedure:: dnload_from_lower
39  procedure:: dnload_from_same
40  procedure:: init
41  procedure:: init_higher
42  procedure:: init_lower
43  procedure:: init_same
44  procedure:: interpolate
45  procedure:: interpolate_unsigned
46  procedure:: interpolate_one
47  procedure:: interpolate_unsigned_one
48  procedure:: remove_hub
49  procedure:: select_dnloader
50  procedure:: update
51  end type
52 CONTAINS
53 
54 !=======================================================================
55 !> Initialise the data communication hub.
56 !=======================================================================
57 SUBROUTINE init (self, link, all_cells)
58  class(data_hub_t):: self
59  class(link_t), pointer :: link
60  logical, optional :: all_cells
61  ! ....................................................................
62  class(patch_t), pointer :: patch
63  ! --------------------------------------------------------------------
64  call trace%begin ('data_hub_t%init')
65  if (self%first_time) then
66  patch => task2patch(link%task) ! convert to patch
67  self%target => patch ! save a pointer to target patch
68  call self%select_dnloader(link, all_cells) ! map the dh-cell shell
69  self%first_time = .false. ! should not be repeated again.
70  end if
71  call trace%end()
72 END SUBROUTINE init
73 
74 !=======================================================================
75 !> Select the download method for each nbor.
76 !=======================================================================
77 SUBROUTINE select_dnloader (self, link, all_cells)
78  class(data_hub_t) :: self
79  class(link_t), pointer :: link
80  logical, optional :: all_cells
81  ! ....................................................................
82  class(link_t), pointer :: nbor
83  class(patch_t), pointer :: source, patch
84  class(data_hub_t), pointer :: hub
85  ! --------------------------------------------------------------------
86  patch => self%target
87  nbor => link%nbor
88  do while (associated(nbor))
89  source => task2patch(nbor%task)
90  if (nbor%download) then
91  if (trim(patch%kind)/='rt_solver') then
92  if ((source%is_clear(bits%no_download))) then
93  if (abs(source%level-patch%level) < 2) then
94  call self%add_hub (patch, source, hub, all_cells)
95  if (source%level == patch%level) then
96  hub%dnload => dnload_from_same ! target and source are of the same level
97  call hub%init_same
98  else if (source%level > patch%level) then
99  hub%dnload => dnload_from_higher ! source is of higher level than the target
100  call hub%init_higher
101  else
102  hub%dnload => dnload_from_lower ! source is of lower level than the target
103  call hub%init_lower
104  end if
105  end if
106  end if
107  else if (trim(source%kind)=='rt_solver') then
108  if (abs(source%level-patch%level) < 2) then
109  call self%add_hub (patch, source, hub, all_cells)
110  hub%rt_target = .true.
111  if (source%level == patch%level) then
112  hub%dnload => dnload_from_same ! target and source are of the same level
113  call hub%init_same
114  else if (source%level > patch%level) then
115  hub%dnload => dnload_from_higher ! source is of higher level than the target
116  call hub%init_higher
117  else
118  hub%dnload => dnload_from_lower ! source is of lower level than the target
119  call hub%init_lower
120  end if
121  end if
122  end if
123  end if
124  nbor => nbor%next ! move to next nbor
125  end do
126 END SUBROUTINE select_dnloader
127 
128 !=======================================================================
129 !> Add a new data hub, responsible for the source in question. Data hubs
130 !> are arranged in the source level of refinement from highest to lowest,
131 !> as for levels, lower than the target, additional cells are needed for
132 !> interpolation.
133 !=======================================================================
134 SUBROUTINE add_hub (self, target, source, hub, all_cells)
135  class(data_hub_t) :: self
136  class(patch_t), pointer :: target, source
137  class(data_hub_t), pointer :: hub
138  logical, optional:: all_cells
139  ! ....................................................................
140  class(data_hub_t), pointer :: next, prev
141  logical :: found
142  ! --------------------------------------------------------------------
143  allocate (hub)
144  hub%target => target
145  hub%source => source
146  hub%level = source%level
147  hub%ID = source%id
148  if (present(all_cells)) hub%all_cells = all_cells
149  if (.not.associated(self%head)) then
150  self%head => hub
151  else
152  nullify (prev)
153  found = .false.
154  next => self%head
155  do while (associated(next).and.(.not.found))
156  if (.not.associated(next%source,hub%source)) then
157  if (next%source%level > hub%source%level) then
158  found = .true.
159  hub%next => next
160  if (associated(prev)) then
161  prev%next => hub
162  else
163  self%head => hub
164  end if
165  end if
166  end if
167  prev => next
168  next => next%next
169  end do
170  if (.not.found) then
171  if (associated(prev)) then
172  prev%next => hub
173  else
174  call io%abort("data_hub_t::add_hub:: location not found?")
175  end if
176  end if
177  end if
178 END SUBROUTINE add_hub
179 
180 !=======================================================================
181 !> Remove a particular data hub. If clear is set to be .true., all hubs
182 !> in the list are removed.
183 !=======================================================================
184 SUBROUTINE remove_hub (self, id, clear)
185  class(data_hub_t) :: self
186  integer, optional :: id
187  logical, optional :: clear
188  ! ....................................................................
189  class(data_hub_t), pointer :: this, next
190  logical:: cleared
191  ! --------------------------------------------------------------------
192  if (present(id)) then
193  cleared = .false.
194  if (associated(self%head)) then
195  this => self%head
196  do while(associated(this).and.(.not.cleared))
197  if (this%id == id) then
198  if (associated(this%prev)) then
199  this%prev%next => this%next
200  this%next%prev => this%prev
201  deallocate(this%s_l, this%s_u)
202  deallocate(this%t_l, this%t_u)
203  deallocate(this)
204  cleared = .true.
205  else
206  if (associated(this%next)) then
207  self%head => this%next
208  else
209  self%first_time = .true.
210  end if
211  deallocate(this%s_l, this%s_u)
212  deallocate(this%t_l, this%t_u)
213  deallocate(this)
214  cleared = .true.
215  end if
216  end if
217  this => this%next
218  end do
219  end if
220  end if
221  if (present(clear)) then
222  this => self%head
223  do while(associated(this))
224  next => this%next
225  deallocate(this%s_l, this%s_u)
226  deallocate(this%t_l, this%t_u)
227  deallocate(this)
228  this => next
229  end do
230  end if
231 END SUBROUTINE remove_hub
232 
233 !=======================================================================
234 !> Initialize the method when both the target and the source are of the
235 !> same level of refinement. Approach is identical to the dnload_same in
236 !> download_mod.
237 !> If "no no-mans land" mode is used, m%n has to be changed into (m%n-1)
238 !=======================================================================
239 SUBROUTINE init_same (self)
240  class(data_hub_t) :: self
241  ! ....................................................................
242  class(mesh_t), pointer ::m
243  real(8) :: dist(3)
244  integer :: i
245  ! --------------------------------------------------------------------
246  dist = self%source%distance(self%target)
247  allocate (self%t_l(3,1))
248  allocate (self%t_u(3,1))
249  allocate (self%s_l(3,1))
250  allocate (self%s_u(3,1))
251  do i = 1,3
252  m => self%target%mesh(i)
253  if (dist(i) > self%source%ds(i)) then
254  self%t_l(i,1) = m%uo
255  self%t_u(i,1) = m%ub
256  self%s_l(i,1) = self%t_l(i,1) - m%n
257  self%s_u(i,1) = self%t_u(i,1) - m%n
258  else if (dist(i) < -self%source%ds(i)) then
259  self%t_l(i,1) = m%lb
260  self%t_u(i,1) = m%lo
261  self%s_l(i,1) = self%t_l(i,1) + m%n
262  self%s_u(i,1) = self%t_u(i,1) + m%n
263  else
264  self%t_l(i,1) = m%li
265  self%t_u(i,1) = m%ui
266  self%s_l(i,1) = self%t_l(i,1)
267  self%s_u(i,1) = self%t_u(i,1)
268  end if
269  end do
270 END SUBROUTINE init_same
271 
272 !=======================================================================
273 !> Download guard cell values from the source, which has the same level
274 !> of refinement, as the target patch. Download range is pre-determined
275 !> by the initialization subroutine above.
276 !=======================================================================
277 SUBROUTINE dnload_from_same (self, only)
278  class(data_hub_t) :: self
279  integer, optional :: only
280  ! ....................................................................
281  integer :: jt(2), iv, sl(3), su(3), tl(3), tu(3), iv1, iv2
282  real :: pt(2)
283  integer, save:: itimer=0
284  ! --------------------------------------------------------------------
285  call trace%begin ('data_hub_t%dnload_same',itimer=itimer)
286  call self%source%time_interval (self%target%time, jt, pt)
287  sl = self%s_l(:,1); su = self%s_u(:,1)
288  tl = self%t_l(:,1); tu = self%t_u(:,1)
289  if (self%rt_target) then
290  self%target%mem(tl(1):tu(1),tl(2):tu(2),tl(3):tu(3),:,self%target%it,1)= &
291  pt(1) *self%source%mem(sl(1):su(1),sl(2):su(2),sl(3):su(3),:,jt(1),1) + &
292  pt(2) *self%source%mem(sl(1):su(1),sl(2):su(2),sl(3):su(3),:,jt(2),1)
293  else
294  if (present(only)) then
295  iv1 = only
296  iv2 = only
297  else
298  iv1 = 1
299  iv2 = self%source%nv
300  end if
301  do iv=iv1, iv2
302  if (self%source%unsigned(iv)) then
303  self%target%mem(tl(1):tu(1),tl(2):tu(2),tl(3):tu(3),iv,self%target%it,1) = exp &
304  (pt(1) *log(self%source%mem(sl(1):su(1),sl(2):su(2),sl(3):su(3),iv,jt(1),1)) + &
305  pt(2) *log(self%source%mem(sl(1):su(1),sl(2):su(2),sl(3):su(3),iv,jt(2),1)))
306  else
307  self%target%mem(tl(1):tu(1),tl(2):tu(2),tl(3):tu(3),iv,self%target%it,1)= &
308  pt(1) *self%source%mem(sl(1):su(1),sl(2):su(2),sl(3):su(3),iv,jt(1),1) + &
309  pt(2) *self%source%mem(sl(1):su(1),sl(2):su(2),sl(3):su(3),iv,jt(2),1)
310  end if
311  end do
312  end if
313  call trace%end(itimer)
314 END SUBROUTINE dnload_from_same
315 
316 !=======================================================================
317 !> Initialize the method when the source is of a higher level of refinement
318 !> than the target patch.
319 !=======================================================================
320 SUBROUTINE init_higher (self)
321  class(data_hub_t) :: self
322  ! ....................................................................
323  class(mesh_t), pointer:: m, m2
324  class(patch_t), pointer:: source, target
325  real(8) :: pos(3)
326  integer :: iv, i, j(3), nv
327  real(8) :: eps, pl(3), pu(3), dp
328  ! --------------------------------------------------------------------
329  source => self%source
330  target => self%target
331  if (self%rt_target) then
332  allocate (self%t_l(3,1))
333  allocate (self%t_u(3,1))
334  nv = 1
335  else
336  allocate (self%t_l(3,target%nv))
337  allocate (self%t_u(3,target%nv))
338  nv = target%nv
339  end if
340  eps = 0.001_8
341  do iv = 1, nv
342  ! --------------------------------------------------------------------
343  ! First and last points inside the source domain
344  ! --------------------------------------------------------------------
345  do i = 1,3
346  m => source%mesh(i)
347  pl(i) = m%r(m%li) - eps
348  pu(i) = m%r(m%ui) + eps
349  end do
350  !---------------------------------------------------------------------
351  ! Map to the first and last point inside the source domain
352  !---------------------------------------------------------------------
353  self%t_l(:,iv) = target%index_only2 (source, pl, iv, roundup=.true.)
354  self%t_u(:,iv) = target%index_only2 (source, pu, iv)
355  ! --------------------------------------------------------------------
356  ! Make sure the target range is legal
357  ! --------------------------------------------------------------------
358  if (.not.self%all_cells) then
359  do i=1,3
360  m => target%mesh(i)
361  m2=> source%mesh(i)
362  self%t_l(i,iv) = max(m%lb,min(m%uo,self%t_l(i,iv)))
363  self%t_u(i,iv) = max(m%lo,min(m%ub,self%t_u(i,iv)))
364  dp = round(m%p-m2%p,12)+m%h(iv)*(m%d-m2%d)
365  if (self%t_u(i,iv)>m%lo) then
366  if (m%r(m%li)+dp > m2%r(m2%ui)) &
367  self%t_u(i,iv) = min(self%t_u(i,iv),m%lo)
368  end if
369  if (self%t_l(i,iv)<m%uo) then
370  if (m%r(m%ui)+dp < m2%r(m2%li)) &
371  self%t_l(i,iv) = max(self%t_l(i,iv),m%uo)
372  end if
373  if (self%t_u(i,iv)>m%ui) then
374  if (m%r(m%uo)+dp > m2%uf) then
375  self%t_u(i,iv) = min(self%t_u(i,iv),m%ui)
376  else
377  self%t_u(i,iv) = max(self%t_l(i,iv),m%ub)
378  end if
379  end if
380  if (self%t_l(i,iv)<m%li) then
381  if (m%r(m%lo)+dp < m2%lf) then
382  self%t_l(i,iv) = max(self%t_l(i,iv),m%li)
383  else
384  self%t_l(i,iv) = min(self%t_l(i,iv),m%lb)
385  end if
386  end if
387  end do
388  end if
389  end do
390 END SUBROUTINE init_higher
391 
392 !=======================================================================
393 !> Download guard cell values from the source, which has a larger level
394 !> of refinement, than the target patch. A temporary scratch cube is
395 !> created and values the right entries in the source memory slots are
396 !> copied. Then an interpolation is done. Can be easily extended to
397 !> accommodate the van Leer interpolation scheme.
398 !=======================================================================
399 SUBROUTINE dnload_from_higher (self, only)
400  class(data_hub_t) :: self
401  integer, optional :: only
402  ! ....................................................................
403  integer :: jt(2), iv, tl(3), tu(3), i, j, k, iv1, iv2
404  integer :: d, l(3), u(3), ii
405  real(kind=KindScalarVar), dimension(:,:,:,:,:), pointer :: mem
406  real(kind=KindScalarVar), dimension(:,:,:,:), pointer :: cube
407  real(kind=KindScalarVar) :: tcell
408  integer, save:: itimer=0
409  real :: pt(2), p(3)
410  class(patch_t), pointer :: source, target
411  class(mesh_t), pointer :: m
412  real(8) :: pos(3)
413  ! --------------------------------------------------------------------
414  call trace%begin ('data_hub_t%dload_higher',itimer=itimer)
415  source => self%source
416  target => self%target
417  allocate(cube(2,2,2,2))
418  call source%time_interval (target%time, jt, pt)
419  if (self%rt_target) then
420  tl = self%t_l(:,1); tu = self%t_u(:,1)
421  do k = tl(3), tu(3)
422  do j = tl(2), tu(2)
423  do i = tl(1), tu(1)
424  l=[i,j,k]
425  do ii = 1,3
426  m => target%mesh(ii)
427  pos(ii) = m%r(l(ii))
428  end do
429  !
430  if ((self%all_cells) .or. &
431  any((l-target%mesh%li) < 0 .or. (l-target%mesh%ui) > 0)) then
432  call source%index_space2(target,pos, 1, l, p)
433  u=l+1
434  mem => source%mem(l(1):u(1),l(2):u(2),l(3):u(3),:,:,1)
435  do iv = 1, target%nv
436  cube(:,:,:,1) = source%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1)
437  cube(:,:,:,2) = source%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1)
438  call self%interpolate(cube, pt, p,tcell)
439  target%mem(i,j,k,iv,target%it,1) = tcell
440  end do
441  end if
442  end do
443  end do
444  end do
445  else
446  if (present(only)) then
447  iv1 = only
448  iv2 = only
449  else
450  iv1 = 1
451  iv2 = target%nv
452  end if
453  d = self%source%idx%d
454  do iv = iv1, iv2
455  tl = self%t_l(:,iv); tu = self%t_u(:,iv)
456  do k = tl(3), tu(3)
457  do j = tl(2), tu(2)
458  do i = tl(1), tu(1)
459  l=[i,j,k]
460  do ii = 1,3
461  m => target%mesh(ii)
462  pos(ii) = m%r(l(ii))
463  end do
464  if ((self%all_cells) .or. &
465  any((l-target%mesh%li) < 0 .or. (l-target%mesh%ui) > 0)) then
466  call source%index_space2(target,pos, iv, l, p)
467  u=l+1
468  mem => source%mem(l(1):u(1),l(2):u(2),l(3):u(3),:,:,1)
469  if (source%pervolume(iv)) then
470  cube(:,:,:,1) = source%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1) / &
471  source%mem(l(1):u(1),l(2):u(2),l(3):u(3), d,jt(1),1)
472  cube(:,:,:,2) = source%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1) / &
473  source%mem(l(1):u(1),l(2):u(2),l(3):u(3), d,jt(2),1)
474  else
475  cube(:,:,:,1) = source%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1)
476  cube(:,:,:,2) = source%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1)
477  end if
478  if (source%unsigned(iv)) then
479  call self%interpolate_unsigned(cube, pt, p,tcell)
480  target%mem(i,j,k,iv,target%it,1) = tcell
481  else if (source%pervolume(iv)) then
482  call self%interpolate(cube, pt, p,tcell)
483  target%mem(i,j,k,iv,target%it,1) = tcell * target%mem(i,j,k,d,target%it,1)
484  else
485  call self%interpolate(cube, pt, p,tcell)
486  target%mem(i,j,k,iv,target%it,1) = tcell
487  end if
488  end if
489  end do
490  end do
491  end do
492  end do
493  end if
494  deallocate(cube)
495  call trace%end(itimer)
496 END SUBROUTINE dnload_from_higher
497 
498 !=======================================================================
499 !> Initialize the method when the source is of a lower level of refinement
500 !> than the target patch.
501 !=======================================================================
502 SUBROUTINE init_lower (self)
503  class(data_hub_t) :: self
504  ! ....................................................................
505  class(mesh_t), pointer:: m, m2
506  class(patch_t), pointer:: source, target
507  integer :: iv, i, j(3), nv
508  real(8) :: eps, pl(3), pu(3), dp
509  ! --------------------------------------------------------------------
510  source => self%source
511  target => self%target
512  eps = 0.001_8
513  if (self%rt_target) then
514  allocate (self%t_l(3,1))
515  allocate (self%t_u(3,1))
516  allocate (self%s_l(3,1))
517  allocate (self%s_u(3,1))
518  nv = 1
519  else
520  allocate (self%t_l(3,target%nv))
521  allocate (self%t_u(3,target%nv))
522  allocate (self%s_l(3,target%nv))
523  allocate (self%s_u(3,target%nv))
524  nv = target%nv
525  end if
526  ! --------------------------------------------------------------------
527  ! Corners of the source domain
528  ! --------------------------------------------------------------------
529  do iv = 1, nv
530  do i = 1,3
531  m => source%mesh(i)
532  pl(i) = m%lf - eps
533  pu(i) = m%uf + eps
534  end do
535  !---------------------------------------------------------------------
536  ! Map to the range inside the target domain
537  !---------------------------------------------------------------------
538  self%t_l(:,iv) = target%index_only2 (source, pl, iv, roundup = .true.)
539  self%t_u(:,iv) = target%index_only2 (source, pu, iv)
540  ! --------------------------------------------------------------------
541  ! Enforce proper limits
542  ! --------------------------------------------------------------------
543  if (.not.self%all_cells) then
544  do i=1,3
545  m => target%mesh(i)
546  m2=> source%mesh(i)
547  self%t_l(i,iv) = max(m%lb,min(m%uo,self%t_l(i,iv)))
548  self%t_u(i,iv) = max(m%lo,min(m%ub,self%t_u(i,iv)))
549  dp = round(m%p-m2%p,12)+m%h(iv)*(m%d-m2%d)
550  if (self%t_u(i,iv)>m%lo) then
551  if (m%r(m%li)+dp > m2%r(m2%ui)) &
552  self%t_u(i,iv) = min(self%t_u(i,iv),m%lo)
553  end if
554  if (self%t_l(i,iv)<m%uo) then
555  if (m%r(m%ui)+dp < m2%r(m2%li)) &
556  self%t_l(i,iv) = max(self%t_l(i,iv),m%uo)
557  end if
558  if (self%t_u(i,iv)>m%ui) then
559  if (m%r(m%uo)+dp > m2%uf) then
560  self%t_u(i,iv) = min(self%t_u(i,iv),m%ui)
561  else
562  self%t_u(i,iv) = max(self%t_l(i,iv),m%ub)
563  end if
564  end if
565  if (self%t_l(i,iv)<m%li) then
566  if (m%r(m%lo)+dp < m2%lf) then
567  self%t_l(i,iv) = max(self%t_l(i,iv),m%li)
568  else
569  self%t_l(i,iv) = min(self%t_l(i,iv),m%lb)
570  end if
571  end if
572  end do
573  end if
574  ! --------------------------------------------------------------------
575  ! Map the acquired indices back to the source
576  ! --------------------------------------------------------------------
577  do i = 1,3
578  m => target%mesh(i)
579  pl(i) = m%r(self%t_l(i,iv))
580  pu(i) = m%r(self%t_u(i,iv))
581  end do
582  self%s_l(:,iv) = source%index_only2 (target,pl, iv)
583  self%s_u(:,iv) = source%index_only2 (target,pu, iv, roundup=.true.)
584  if (.not.self%all_cells) then
585  do i = 1,3
586  m => source%mesh(i)
587  self%s_l(i,iv) = max(m%lo,min(m%ui,self%s_l(i,iv)))
588  self%s_u(i,iv) = min(m%uo,max(m%li,self%s_u(i,iv)))
589  end do
590  end if
591  end do
592 END SUBROUTINE init_lower
593 
594 !=======================================================================
595 !> Download guard cell values from the source, which has a lower level
596 !> of refinement, than the target patch. A temporary scratch array is
597 !> created and are taken from the source patch. If the needed values lay
598 !> outside of the source%li/%ui, then a search in the target is done. If
599 !> the right location for the data is not in the target, a search in the
600 !> neighbors is done. All values are initially interpolated in time to the
601 !> target%time and then interpolated in space to get the required quantities.
602 !> Can be extended to accommodate the van Leer interpolation scheme.
603 !=======================================================================
604 SUBROUTINE dnload_from_lower (self, only)
605  class(data_hub_t) :: self
606  integer, optional :: only
607  ! ....................................................................
608  integer :: jt(2),jt2(2), iv, sl(3), su(3), tl(3), tu(3), i, j, k
609  integer :: d, cl(3), ii, l(3), u(3), id, n(3), iv1, iv2
610  class(mesh_t), pointer :: m
611  class(patch_t), pointer:: source, target
612  real(kind=KindScalarVar), dimension(:,:,:,:,:), pointer :: mem
613  real(kind=KindScalarVar), dimension(:,:,:), pointer :: scr
614  real(kind=KindScalarVar), dimension(:,:,:,:), pointer :: rt_scr, rt_cube
615  real(kind=KindScalarVar), dimension(:,:,:), pointer :: cube
616  real(kind=KindScalarVar), pointer :: scr3
617  real(kind=KindScalarVar), dimension(:), pointer :: rt_scr3
618  real(kind=KindScalarVar) :: tcell
619  real :: pt(2), pt2(2), p(3)
620  real(8) :: pos(3), time, pos2(3)
621  logical :: outside, found
622  integer, save:: itimer=0
623  ! --------------------------------------------------------------------
624  call trace%begin ('data_hub_t%dnload_lower',itimer=itimer)
625  source => self%source
626  target => self%target
627  d = source%idx%d
628  mem => self%source%mem(:,:,:,:,:,1)
629  if (present(only)) then
630  iv1 = only
631  iv2 = only
632  else
633  iv1 = 1
634  iv2 = source%nv
635  end if
636  ! --------------------------------------------------------------------
637  ! interpolate in time and fill in a scratch cube. When downloading from
638  ! a larger patch, not all needed cells are usually available inside the
639  ! source patch, and have to be retrieved either from the target or from
640  ! other nbors.
641  ! --------------------------------------------------------------------
642  call source%time_interval (target%time, jt2, pt2)
643  time = target%time
644  id = -1
645  if (self%rt_target) then
646  sl = self%s_l(:,1); su = self%s_u(:,1)
647  tl = self%t_l(:,1); tu = self%t_u(:,1)
648  n = su - sl +2
649  allocate (rt_scr(n(1), n(2), n(3),target%nv))
650  do k = 1, n(3)
651  do j = 1, n(2)
652  do i = 1, n(1)
653  cl = [sl(1)+i-1,sl(2)+j-1,sl(3)+k-1]
654  outside = .false.
655  do ii = 1,3
656  if ((cl(ii)<source%li(ii)).or.(cl(ii)>source%ui(ii))) outside = .true.
657  end do
658  if (outside) then
659  found = .false.
660  do ii = 1,3
661  m => source%mesh(ii)
662  pos(ii) = m%r(cl(ii))
663  pos2(ii) = m%r(cl(ii))+m%p+m%h(1)*m%d
664  end do
665  rt_scr3 => rt_scr(i,j,k,:)
666  call self%check_in_target(pos, 1, time, jt, pt, id, found, rt_scr=rt_scr3) ! many cells are needed from the target for good interpolation
667  if (.not.found) &
668  call self%check_in_nbors (pos, 1, time, jt, pt, id, found, rt_scr=rt_scr3) ! worst case.. at least necessary nbors will be same size or smaller.
669  if (.not.found) then
670  print *, iv, cl, pos
671  print *, pos2
672  call io%abort("data_hub::dnload_lower:: cell not found") ! very unusual, but a necessary check and a stop.
673  end if
674  else
675  rt_scr(i,j,k,:) = pt2(1) * mem(cl(1),cl(2),cl(3),:, jt2(1)) + &
676  pt2(2) * mem(cl(1),cl(2),cl(3),:, jt2(2))
677  end if
678  end do
679  end do
680  end do
681  ! -------------------------------------------------------------------
682  ! Now that the scratch cube is ready, we can spatially interpolate
683  ! back to target
684  ! -------------------------------------------------------------------
685  do k=tl(3), tu(3)
686  do j=tl(2), tu(2)
687  do i=tl(1), tu(1)
688  cl = [i,j,k]
689  if ((self%all_cells) .or. &
690  any((cl-target%mesh%li) < 0 .or. (cl-target%mesh%ui) > 0)) then
691  do ii = 1,3
692  m => target%mesh(ii)
693  pos(ii) = m%r(cl(ii))
694  end do
695  call source%index_space2(target,pos, 1, cl, p)
696  l = cl - sl + 1
697  u = l +1
698  do iv = 1, target%nv
699  cube => rt_scr(l(1):u(1),l(2):u(2),l(3):u(3), iv)
700  call self%interpolate_one(cube, p,tcell)
701  target%mem(i,j,k,iv,target%it,1) = tcell
702  end do
703  end if
704  end do
705  end do
706  end do
707  deallocate(rt_scr)
708  else
709  do iv = iv1, iv2
710  sl = self%s_l(:,iv); su = self%s_u(:,iv)
711  tl = self%t_l(:,iv); tu = self%t_u(:,iv)
712  n = su - sl +2
713  allocate (scr(n(1), n(2), n(3)))
714  do k = 1, n(3)
715  do j = 1, n(2)
716  do i = 1, n(1)
717  cl = [sl(1)+i-1,sl(2)+j-1,sl(3)+k-1]
718  outside = .false.
719  do ii = 1,3
720  if ((cl(ii)<source%li(ii)).or.(cl(ii)>source%ui(ii))) outside = .true.
721  end do
722  if (outside) then
723  found = .false.
724  do ii = 1,3
725  m => source%mesh(ii)
726  pos(ii) = m%r(cl(ii))
727  pos2(ii) = m%r(cl(ii))+m%p+m%h(iv)*m%d
728  end do
729  scr3 => scr(i,j,k)
730  call self%check_in_target(pos, iv, time, jt, pt, id, found, scr=scr3) ! many cells are needed from the target for good interpolation
731  if (.not.found) &
732  call self%check_in_nbors (pos, iv, time, jt, pt, id, found, scr=scr3) ! worst case.. at least necessary nbors will be same size or smaller.
733  if (.not.found) then
734  print *, iv, cl, pos
735  print *, pos2
736  call io%abort("data_hub::dnload_lower:: cell not found") ! very unusual, but a necessary check and a stop.
737  end if
738  else
739  if (source%unsigned(iv)) then
740  scr(i,j,k) = exp(pt2(1) * log(mem(cl(1),cl(2),cl(3),iv, jt2(1))) + &
741  pt2(2) * log(mem(cl(1),cl(2),cl(3),iv, jt2(2))))
742  else if (source%pervolume(iv)) then
743  scr(i,j,k) = pt2(1) * (mem(cl(1),cl(2),cl(3),iv, jt2(1)) / &
744  mem(cl(1),cl(2),cl(3), d, jt2(1))) + &
745  pt2(2) * (mem(cl(1),cl(2),cl(3),iv, jt2(2)) / &
746  mem(cl(1),cl(2),cl(3), d, jt2(2)))
747  else
748  scr(i,j,k) = pt2(1) * mem(cl(1),cl(2),cl(3),iv, jt2(1)) + &
749  pt2(2) * mem(cl(1),cl(2),cl(3),iv, jt2(2))
750  end if
751  end if
752  end do
753  end do
754  end do
755  ! -------------------------------------------------------------------
756  ! Now that the scratch cube is ready, we can spatially interpolate
757  ! back to target
758  ! -------------------------------------------------------------------
759  do k=tl(3), tu(3)
760  do j=tl(2), tu(2)
761  do i=tl(1), tu(1)
762  cl = [i,j,k]
763  if ((self%all_cells) .or. &
764  any((cl-target%mesh%li) < 0 .or. (cl-target%mesh%ui) > 0)) then
765  do ii = 1,3
766  m => target%mesh(ii)
767  pos(ii) = m%r(cl(ii))
768  end do
769  call source%index_space2(target,pos, iv, cl, p)
770  l = cl - sl + 1
771  u = l +1
772  cube => scr(l(1):u(1),l(2):u(2),l(3):u(3))
773  if (source%unsigned(iv)) then
774  call self%interpolate_unsigned_one(cube, p,tcell)
775  target%mem(i,j,k,iv,target%it,1) = tcell
776  else if (self%source%pervolume(iv)) then
777  call self%interpolate_one(cube, p,tcell)
778  target%mem(i,j,k,iv,target%it,1) = tcell * target%mem(i,j,k,d,target%it,1)
779  else
780  call self%interpolate_one(cube, p,tcell)
781  target%mem(i,j,k,iv,target%it,1) = tcell
782  end if
783  end if
784  end do
785  end do
786  end do
787  deallocate(scr)
788  end do
789  end if
790  call trace%end(itimer)
791 END SUBROUTINE dnload_from_lower
792 
793 !=======================================================================
794 !> Some support data is needed, when download_lower is done. Usually it
795 !> resides in the target patch. This data is spatially interpolated to
796 !> the source resolution and added to the scratch array.
797 !=======================================================================
798 SUBROUTINE check_in_target (self, pos, iv, time, jt, pt, id, found, scr, rt_scr)
799  class(data_hub_t) :: self
800  real(8) :: pos(3), time
801  real(kind=KindScalarVar), pointer, optional :: scr
802  real(kind=KindScalarVar), dimension(:),pointer, optional :: rt_scr
803  integer :: iv, jt(2), id
804  real :: pt(2)
805  logical :: found
806  ! ....................................................................
807  integer :: sl(3), su(3), tl(3), tu(3), i, j, k, aa, bb, cc
808  integer :: d, cl(3), ii, l(3), u(3), v
809  real :: p(3)
810  real(8) :: eps = 0.001_8
811  class(mesh_t), pointer :: m
812  class(patch_t), pointer:: source, target
813  real(kind=KindScalarVar), dimension(:,:,:,:), pointer :: cube
814  logical :: inside
815  ! --------------------------------------------------------------------
816  source => self%source
817  target => self%target
818  d = source%idx%d
819  inside = .true.
820  call target%index_space2(source,pos,iv, cl, p)
821  inside = all(cl>=target%mesh%li .and. cl<=target%mesh%ui)
822  if (inside) then ! do the interpolation from smaller
823  ! --------------------------------------------------------------------
824  ! do spatial interpolation.
825  ! --------------------------------------------------------------------
826  if (id /=target%id) then
827  call target%time_interval (time, jt, pt)
828  id = target%id
829  end if
830  allocate(cube(2,2,2,2))
831  l = cl
832  u = cl+1
833  if (self%rt_target) then
834  do v=1,target%nv
835  cube(:,:,:,1) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),v,jt(1),1)
836  cube(:,:,:,2) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),v,jt(2),1)
837  call self%interpolate(cube, pt, p,rt_scr(v))
838  end do
839  else
840  if (source%unsigned(iv)) then
841  cube(:,:,:,1) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1)
842  cube(:,:,:,2) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1)
843  call self%interpolate_unsigned(cube, pt, p,scr)
844  else if (source%pervolume(iv)) then
845  cube(:,:,:,1) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1) / &
846  target%mem(l(1):u(1),l(2):u(2),l(3):u(3), d,jt(1),1)
847  cube(:,:,:,2) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1) / &
848  target%mem(l(1):u(1),l(2):u(2),l(3):u(3), d,jt(2),1)
849  call self%interpolate(cube, pt, p,scr)
850  else
851  cube(:,:,:,1) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1)
852  cube(:,:,:,2) = target%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1)
853  call self%interpolate(cube, pt, p,scr)
854  end if
855  end if
856  deallocate(cube)
857  found = .true.
858  end if
859 END SUBROUTINE check_in_target
860 
861 !=======================================================================
862 !> In some cases, e.g. corners, ends of edges, some support data is located
863 !> not in the target patch, but rather in adjacent nbors. They have the
864 !> same level of refinement as the target or the source. This data is
865 !> spatially interpolated to the source resolution, if needed as well as
866 !> interpolated in time to the needed target time. Then it is added
867 !> to the scratch array.
868 !=======================================================================
869 SUBROUTINE check_in_nbors (self, pos, iv, time, jt, pt, id, found, scr, rt_scr)
870  class(data_hub_t) :: self
871  real(8) :: pos(3), time
872  real(kind=KindScalarVar), pointer, optional :: scr
873  real(kind=KindScalarVar), dimension(:),pointer, optional :: rt_scr
874  integer :: jt(2), id
875  real :: pt(2)
876  logical :: found
877  ! ....................................................................
878  integer :: iv, sl(3), su(3), tl(3), tu(3), i, j, k
879  integer :: d, cl(3), ii, l(3), u(3), v
880  real :: p(3)
881  real(8) :: eps = 0.001_8
882  class(mesh_t), pointer :: m
883  class(link_t), pointer :: nbor
884  class(patch_t), pointer:: source, target, nbpatch
885  real(kind=KindScalarVar), dimension(:,:,:,:), pointer :: cube
886  logical :: inside
887  ! --------------------------------------------------------------------
888  source => self%source
889  target => self%target
890  d = source%idx%d
891  nbor => source%link%nbor
892  do while((associated(nbor)).and..not.found)
893  if ((nbor%task%id /= source%id).and.(nbor%task%id /= target%id)) then ! source is already checked, target should not be there
894  nbpatch => task2patch(nbor%task)
895  if (abs(nbpatch%level-source%level)<2) then
896  if (nbpatch%level == source%level) then ! same level
897  cl = nbpatch%index_only2 (source,pos, iv, closest=.true.)
898  inside = all(cl>=nbpatch%mesh%li .and. cl<=nbpatch%mesh%ui)
899  if (inside) then ! same size, just take the value by temporal interpolation
900  if (id /=nbpatch%id) then
901  call nbpatch%time_interval (time, jt, pt)
902  id = nbpatch%id
903  end if
904  if(self%rt_target) then
905  rt_scr = &
906  pt(1) * nbpatch%mem(cl(1),cl(2),cl(3),:,jt(1),1) + &
907  pt(2) * nbpatch%mem(cl(1),cl(2),cl(3),:,jt(2),1)
908  found = .true.
909  else
910  if (source%unsigned(iv)) then
911  scr = exp( &
912  pt(1) * log(nbpatch%mem(cl(1),cl(2),cl(3),iv,jt(1),1)) + &
913  pt(2) * log(nbpatch%mem(cl(1),cl(2),cl(3),iv,jt(2),1)))
914  else if (source%pervolume(iv)) then
915  scr = &
916  pt(1) * (nbpatch%mem(cl(1),cl(2),cl(3),iv,jt(1),1) / &
917  nbpatch%mem(cl(1),cl(2),cl(3), d,jt(1),1))+ &
918  pt(2) * (nbpatch%mem(cl(1),cl(2),cl(3),iv,jt(2),1) / &
919  nbpatch%mem(cl(1),cl(2),cl(3), d,jt(2),1))
920  else
921  scr = &
922  pt(1) * nbpatch%mem(cl(1),cl(2),cl(3),iv,jt(1),1) + &
923  pt(2) * nbpatch%mem(cl(1),cl(2),cl(3),iv,jt(2),1)
924  end if
925  found = .true.
926  end if
927  end if
928  else
929  call nbpatch%index_space2(source,pos,iv, cl, p)
930  inside = all(cl>=nbpatch%mesh%li .and. cl<=nbpatch%mesh%ui)
931  if (inside) then ! different size, spatial and temporal interpolation
932  if (id /=nbpatch%id) then
933  call nbpatch%time_interval (time, jt, pt)
934  id = nbpatch%id
935  end if
936  allocate(cube(2,2,2,2))
937  l = cl
938  u = cl+1
939  if (self%rt_target) then
940  do v = 1, target%nv
941  cube(:,:,:,1) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),v,jt(1),1)
942  cube(:,:,:,2) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),v,jt(2),1)
943  call self%interpolate(cube, pt, p,rt_scr(v))
944  end do
945  else
946  if (source%unsigned(iv)) then
947  cube(:,:,:,1) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1)
948  cube(:,:,:,2) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1)
949  call self%interpolate_unsigned(cube, pt, p,scr)
950  else if (source%pervolume(iv)) then
951  cube(:,:,:,1) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1) / &
952  nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3), d,jt(1),1)
953  cube(:,:,:,2) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1) / &
954  nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3), d,jt(2),1)
955  call self%interpolate(cube, pt, p,scr)
956  else
957  cube(:,:,:,1) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(1),1)
958  cube(:,:,:,2) = nbpatch%mem(l(1):u(1),l(2):u(2),l(3):u(3),iv,jt(2),1)
959  call self%interpolate(cube, pt, p,scr)
960  end if
961  end if
962  deallocate(cube)
963  found = .true.
964  end if
965  end if
966  end if
967  end if
968  nbor => nbor%next
969  end do
970 END SUBROUTINE check_in_nbors
971 
972 !=======================================================================
973 !> Update the guard zone values
974 !=======================================================================
975 SUBROUTINE update (self, link, all_cells, only)
976  class(data_hub_t):: self
977  class(link_t), pointer :: link
978  logical, optional :: all_cells
979  integer, optional :: only
980  ! ....................................................................
981  class(data_hub_t), pointer :: hub
982  integer, save:: itimer=0
983  ! --------------------------------------------------------------------
984  !call trace%begin ('data_hub_t%update',itimer=itimer)
985  if (self%first_time) call self%init(link, all_cells) ! map data shell if it is the first time
986  !
987  hub => self%head
988  do while(associated(hub))
989  call hub%dnload(only=only)
990  hub => hub%next
991  end do
992  !call trace%end(itimer)
993 END SUBROUTINE update
994 
995 SUBROUTINE interpolate (self, src, pt, p, res)
996  class(data_hub_t) :: self
997  real(kind=KindScalarVar), dimension(:,:,:,:):: src
998  real :: pt(2), p(3)
999  real(kind=KindScalarVar) :: res
1000  !---------------------------------------------------------------------
1001  if (any(p < 0.0).or.any(p>1.0)) then
1002  print *, p
1003  call io%abort("p negative")
1004  end if
1005  res = &
1006  pt(1)*((1.-p(3))*((1.-p(2))*((1.-p(1))*src(1,1,1,1) + &
1007  ( p(1))*src(2,1,1,1)) + &
1008  ( p(2))*((1.-p(1))*src(1,2,1,1) + &
1009  ( p(1))*src(2,2,1,1))) + &
1010  ( p(3))*((1.-p(2))*((1.-p(1))*src(1,1,2,1) + &
1011  ( p(1))*src(2,1,2,1)) + &
1012  ( p(2))*((1.-p(1))*src(1,2,2,1) + &
1013  ( p(1))*src(2,2,2,1)))) + &
1014  pt(2)*((1.-p(3))*((1.-p(2))*((1.-p(1))*src(1,1,1,2) + &
1015  ( p(1))*src(2,1,1,2)) + &
1016  ( p(2))*((1.-p(1))*src(1,2,1,2) + &
1017  ( p(1))*src(2,2,1,2))) + &
1018  ( p(3))*((1.-p(2))*((1.-p(1))*src(1,1,2,2) + &
1019  ( p(1))*src(2,1,2,2)) + &
1020  ( p(2))*((1.-p(1))*src(1,2,2,2) + &
1021  ( p(1))*src(2,2,2,2))))
1022 END SUBROUTINE interpolate
1023 
1024 SUBROUTINE interpolate_unsigned (self, src, pt, p, res)
1025  class(data_hub_t) :: self
1026  real(kind=KindScalarVar), dimension(:,:,:,:):: src
1027  real :: pt(2), p(3)
1028  real(kind=KindScalarVar) :: res
1029  !---------------------------------------------------------------------
1030  if (any(p < 0.0).or.any(p>1.0)) then
1031  print *, p
1032  call io%abort("p negative")
1033  end if
1034  res = exp( &
1035  pt(1)*((1.-p(3))*((1.-p(2))*((1.-p(1))*log(src(1,1,1,1)) + &
1036  ( p(1))*log(src(2,1,1,1))) + &
1037  ( p(2))*((1.-p(1))*log(src(1,2,1,1)) + &
1038  ( p(1))*log(src(2,2,1,1)))) + &
1039  ( p(3))*((1.-p(2))*((1.-p(1))*log(src(1,1,2,1)) + &
1040  ( p(1))*log(src(2,1,2,1))) + &
1041  ( p(2))*((1.-p(1))*log(src(1,2,2,1)) + &
1042  ( p(1))*log(src(2,2,2,1))))) + &
1043  pt(2)*((1.-p(3))*((1.-p(2))*((1.-p(1))*log(src(1,1,1,2)) + &
1044  ( p(1))*log(src(2,1,1,2))) + &
1045  ( p(2))*((1.-p(1))*log(src(1,2,1,2)) + &
1046  ( p(1))*log(src(2,2,1,2)))) + &
1047  ( p(3))*((1.-p(2))*((1.-p(1))*log(src(1,1,2,2)) + &
1048  ( p(1))*log(src(2,1,2,2))) + &
1049  ( p(2))*((1.-p(1))*log(src(1,2,2,2)) + &
1050  ( p(1))*log(src(2,2,2,2))))))
1051 END SUBROUTINE interpolate_unsigned
1052 
1053 SUBROUTINE interpolate_one (self, src, p, res)
1054  class(data_hub_t) :: self
1055  real(kind=KindScalarVar), dimension(:,:,:):: src
1056  real :: p(3)
1057  real(kind=KindScalarVar) :: res
1058  !---------------------------------------------------------------------
1059  if (any(p < 0.0).or.any(p>1.0)) then
1060  print *, p
1061  call io%abort("p negative")
1062  end if
1063  res = &
1064  ((1.-p(3))*((1.-p(2))*((1.-p(1))*src(1,1,1) + &
1065  ( p(1))*src(2,1,1)) + &
1066  ( p(2))*((1.-p(1))*src(1,2,1) + &
1067  ( p(1))*src(2,2,1))) + &
1068  ( p(3))*((1.-p(2))*((1.-p(1))*src(1,1,2) + &
1069  ( p(1))*src(2,1,2)) + &
1070  ( p(2))*((1.-p(1))*src(1,2,2) + &
1071  ( p(1))*src(2,2,2))))
1072 END SUBROUTINE interpolate_one
1073 
1074 SUBROUTINE interpolate_unsigned_one (self, src, p, res)
1075  class(data_hub_t) :: self
1076  real(kind=KindScalarVar), dimension(:,:,:):: src
1077  real :: p(3)
1078  real(kind=KindScalarVar) :: res
1079  !---------------------------------------------------------------------
1080  if (any(p < 0.0).or.any(p>1.0)) then
1081  print *, p
1082  call io%abort("p negative")
1083  end if
1084  res = exp &
1085  ((1.-p(3))*((1.-p(2))*((1.-p(1))*log(src(1,1,1)) + &
1086  ( p(1))*log(src(2,1,1))) + &
1087  ( p(2))*((1.-p(1))*log(src(1,2,1)) + &
1088  ( p(1))*log(src(2,2,1)))) + &
1089  ( p(3))*((1.-p(2))*((1.-p(1))*log(src(1,1,2)) + &
1090  ( p(1))*log(src(2,1,2))) + &
1091  ( p(2))*((1.-p(1))*log(src(1,2,2)) + &
1092  ( p(1))*log(src(2,2,2)))))
1093 END SUBROUTINE interpolate_unsigned_one
1094 
1095 END MODULE data_hub_mod
A communications hub between patches. Three methods are used: download_same, download_higher and down...
Definition: data_hub_mod.f90:5
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4
Template module for tasks.
Definition: task_mod.f90:4