DISPATCH
guard_zones_mod.f90
1 !===============================================================================
2 !> Optimized restrict operations, intended for guard zones (not conservative)
3 !>
4 !> 1 2 3 4 5 6 7 8 9 0 11 17 18 19 20 21 22
5 !> +---+---+-|-o---o---o---o---o---o---0---0-| --0---o---o---o-|-+---+---+
6 !> -o-o-o-o-o-o|+-+-+ | +-+-+|o-o-o-o-o-o
7 !> 4 5 6 7 8 9 0 1 2 1 2 3 4 5 6 7 8 9
8 !> o-o-o|o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o|o-o-o
9 !> 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2
10 !>
11 !> When boundaries are lined up for 16x16x16 patches, restrict operations should
12 !> involve interpolations over index 14-15, 16-17, and 18-19 on the high source
13 !> side, and interpolations over index 4-5, 6-7, and 8-9 on the low source side,
14 !> while of coure 1-3 and 20-22 on the low and high target side. When the
15 !> source overlaps with the target it covers target index 4-11 left and 12-19
16 !> right, with its own range 4-19.
17 !===============================================================================
19  USE io_mod
20  USE io_unit_mod
21  USE trace_mod
22  USE patch_mod
23  USE mesh_mod
24  USE omp_timer_mod
25  implicit none
26  private
27  type, public:: guard_zones_t
28  procedure(restrict_simple), pointer:: method => null()
29  contains
30  procedure:: init
31  procedure:: restrict_simple
32  procedure:: different_opt
33  procedure:: restrict_vect
34  procedure, nopass:: info
35  end type
36  real(8):: update_t(3)=0d0
37  integer:: verbose=0
38  logical:: detailed_timer=.false.
39  type(guard_zones_t), public:: guard_zones
40 CONTAINS
41 
42 !===============================================================================
43 !> Choose method and verbosity
44 !===============================================================================
45 SUBROUTINE init (self)
46  class(guard_zones_t):: self
47  character(len=6):: method='opt'
48  namelist /guard_zone_params/ method, detailed_timer, verbose
49  integer:: iostat
50  logical, save:: first_time=.true.
51  !-----------------------------------------------------------------------------
52  call trace%begin ('guard_zones_t%init')
53  if (first_time) then
54  !$omp critical (init_cr)
55  if (first_time) then
56  rewind(io%input)
57  read (io%input, guard_zone_params, iostat=iostat)
58  write (io%output, guard_zone_params)
59  first_time = .false.
60  end if
61  !$omp end critical (init_cr)
62  end if
63  select case (trim(method))
64  case ('simple')
65  self%method => restrict_simple
66  case ('opt')
67  self%method => different_opt
68  case ('vect')
69  self%method => restrict_vect
70  case default
71  call io%abort ('guard zone method unknown')
72  end select
73  call trace%end ()
74 END SUBROUTINE init
75 
76 !===============================================================================
77 !> Compute interpolated guard zone values by setting up interface to the
78 !> slow but correct patch_t%interpolate() point-by-point interpolation.
79 !===============================================================================
80 SUBROUTINE restrict_simple (self, target, source, jt, pt, same)
81  class(guard_zones_t):: self
82  class(patch_t), pointer:: target, source
83  integer:: jt(2)
84  real(4):: pt(2)
85  logical, optional:: same
86  !.............................................................................
87  real, dimension(:,:,:), pointer:: tmem, smem, td, sd
88  class(mesh_t), pointer:: mesh
89  integer:: i1, i2, i3, iv, jv, tl(3), tu(3), sl(3), su(3)
90  real(8):: tlc(3), tuc(3), slc(3), suc(3), dist(3)
91  real(4):: pl(3), pu(3)
92  real(8):: used
93  integer, save:: itimer=0
94  !-----------------------------------------------------------------------------
95  ! Interpolate in time and space, point-by-point
96  !-----------------------------------------------------------------------------
97  call trace%begin ('guard_zones_t%restrict_simple', itimer=itimer)
98  if (detailed_timer) &
99  used = wallclock()
100  call source%time_interval (target%time, jt, pt)
101  do iv=1,target%nv
102  jv = 0
103  !---------------------------------------------------------------------------
104  ! Target lower and upper corner in real space
105  !---------------------------------------------------------------------------
106  tlc = target%myposition (target%mesh%lb, jv)
107  tuc = target%myposition (target%mesh%ub, jv)
108  !---------------------------------------------------------------------------
109  ! Map to source index space and limit to source interior
110  !---------------------------------------------------------------------------
111  sl = source%index_only (tlc, jv, nowrap=.true.)
112  su = source%index_only (tuc, jv, nowrap=.true.)
113  sl = min(max(sl,source%mesh%li),source%mesh%ui-1)
114  su = min(max(su,source%mesh%li),source%mesh%ui-1)
115  !---------------------------------------------------------------------------
116  ! Map source points to target legal index space
117  !---------------------------------------------------------------------------
118  slc = source%myposition(sl, jv)
119  suc = source%myposition(su, jv)
120  tl = target%index_only (slc, jv, nowrap=.true., roundup=.true.)
121  tu = target%index_only (suc, jv, nowrap=.true., roundup=.true.)
122  tl = min(max(tl,target%mesh%lb),target%mesh%ub)
123  tu = min(max(tu,target%mesh%lb),target%mesh%ub)
124  if (verbose > 0) &
125  write(io_unit%log,'(a,i3,2(2x,3i4),2(2x,3f7.3))') &
126  'guard_zones_t%restrict_simple1: iv, tl, tu =', iv, tl, tu
127  !call dnload_range (target, source, iv, tl, tu)
128  !if (verbose > 0) &
129  ! write(io_unit%log,'(a,i3,2(2x,3i4),2(2x,3f7.3))') &
130  ! 'guard_zones_t%restrict_simple2: iv, tl, tu =', iv, tl, tu
131  do i3=tl(3),tu(3); do i2=tl(2),tu(2); do i1=tl(1),tu(1)
132  call interpolate (source, target, [i1,i2,i3], iv, jt, pt)
133  end do; end do; end do
134  end do
135  if (detailed_timer) then
136  used = wallclock()-used
137  write (stdout,'(a,6p,f6.3," mus/cell")') &
138  "guard_zone_t%restrict_simple:", used/product(tu-tl+1)
139  end if
140  call trace%end (itimer)
141 END SUBROUTINE restrict_simple
142 
143 !===============================================================================
144 !> Compute interpolated guard zone values by setting up interface to the
145 !> slow but correct patch_t%interpolate() point-by-point interpolation.
146 !===============================================================================
147 SUBROUTINE different_opt (self, target, source, jt, pt, same)
148  class(guard_zones_t):: self
149  class(patch_t), pointer:: target, source
150  integer:: jt(2)
151  real(4):: pt(2)
152  logical, optional:: same
153  !.............................................................................
154  real, parameter:: eps=1e-4
155  real, dimension(:,:,:), pointer:: tmem, smem, td, sd
156  class(mesh_t), pointer:: mesh
157  integer:: i, i1, i2, i3, iv, jv, tl(3), tu(3), sl(3), su(3), cells
158  real:: pi(16)
159  real(8):: tlc, tuc, slc, suc, dist, box, tc
160  class(mesh_t), pointer:: tm, sm
161  real(8):: start
162  integer, save:: itimer=0
163  !-----------------------------------------------------------------------------
164  ! Interpolate in time and space, point-by-point
165  !-----------------------------------------------------------------------------
166  call trace%begin ('guard_zones_t%different_opt', itimer=itimer)
167  if (detailed_timer) then
168  start = wallclock()
169  end if
170  call source%time_interval (target%time, jt, pt)
171  do iv=1,target%nv
172  jv = 0
173  do i=1,3
174  !-------------------------------------------------------------------------
175  ! Source distance from target center
176  !-------------------------------------------------------------------------
177  tm => target%mesh(i)
178  sm => source%mesh(i)
179  box = tm%b
180  dist = sm%p - tm%p
181  if (target%periodic(i)) &
182  dist = modulo(dist+0.5d0*box,box) - 0.5d0*box
183  !-------------------------------------------------------------------------
184  ! Target corner distances from center
185  !-------------------------------------------------------------------------
186  tlc = tm%r(tm%lb)
187  tuc = tm%r(tm%ub)
188  !-------------------------------------------------------------------------
189  ! Source lower and upper li and ui points, relative to target center
190  !-------------------------------------------------------------------------
191  slc = sm%r(sm%li) + dist
192  suc = sm%r(sm%ui) + dist
193  !-------------------------------------------------------------------------
194  ! Lower and upper index distance from target center
195  !-------------------------------------------------------------------------
196  tlc = max(tlc,slc)
197  tuc = min(tuc,suc)
198  !-------------------------------------------------------------------------
199  ! Target indices to lower and upper guard zone limits
200  !-------------------------------------------------------------------------
201  tc = 1d0/tm%d
202  tl(i) = ceiling(tlc*tc + tm%o - eps)
203  tu(i) = ceiling(tuc*tc + tm%o - eps)
204  if (verbose > 0) &
205  write(io_unit%log,'(a,i3,2x,5f8.3,2x,2i4)') &
206  'guard_zones_t%different_opt: i,dist,tlc,tuc,slc,suc,tl,tu =', &
207  i, dist, tlc, tuc, slc, suc, tl(i), tu(i)
208  end do
209  !call dnload_range (target, source, iv, tl, tu)
210  !if (verbose > 0) &
211  ! write(io_unit%log,'(a,i3,2(2x,3i4),2(2x,3f7.3))') &
212  ! 'guard_zones_t%restrict_simple2: iv, tl, tu =', iv, tl, tu
213  if (present(same)) then
214  call interpolate (source, target, [i1,i2,i3], iv, jt, pt, pi=pi, set=.true.)
215  do i3=tl(3),tu(3); do i2=tl(2),tu(2); do i1=tl(1),tu(1)
216  call interpolate (source, target, [i1,i2,i3], iv, jt, pt, pi=pi)
217  end do; end do; end do
218  else
219  do i3=tl(3),tu(3); do i2=tl(2),tu(2); do i1=tl(1),tu(1)
220  call interpolate (source, target, [i1,i2,i3], iv, jt, pt)
221  end do; end do; end do
222  end if
223  end do
224  if (detailed_timer) then
225  cells = product(tu-tl+1)
226  write (stdout,'(a,6p,2(2x,i6,i3),i6," cells,",f9.3," mus/cell")') &
227  "guard_zones_t%different_opt: target, level, source, level,", &
228  target%id, target%level, source%id, source%level, cells, &
229  (wallclock()-start)/max(1,cells)
230  end if
231  call trace%end (itimer)
232 END SUBROUTINE different_opt
233 
234 !===============================================================================
235 !> Compute interpolated guard zone values
236 !===============================================================================
237 SUBROUTINE restrict_vect (self, target, source, jt, pt, same)
238  class(guard_zones_t):: self
239  class(patch_t), pointer:: target, source
240  integer:: jt(2)
241  real(4):: pt(2)
242  logical, optional:: same
243  !.............................................................................
244  real, dimension(:,:,:), pointer:: tmem, smem, td, sd
245  class(mesh_t), pointer:: mesh
246  integer:: i, i1, i2, i3, j, k, iv, tl(3), tu(3), sl(3), su(3), j0(3)
247  real(8):: tlc(3), tuc(3), slc(3), suc(3)
248  real(4):: p(3), p1, p2, p3, w(8)
249  real(8):: used
250  integer, save:: itimer=0
251  !-----------------------------------------------------------------------------
252  ! Map the source corners into target index space
253  !-----------------------------------------------------------------------------
254  call trace%begin ('guard_zones_t%restrict_vect', itimer=itimer)
255  if (detailed_timer) &
256  used = wallclock()
257  !-----------------------------------------------------------------------------
258  ! Interpolate in time and space
259  !-----------------------------------------------------------------------------
260  call source%time_interval (target%time, jt, pt)
261  if (verbose > 0) &
262  write(io_unit%log,'(a,2(2x,i6,i3,2x,3f7.3),2x,3f7.3)') &
263  'target, level, lc, source, level, lc, relative =', &
264  target%id, target%level, target%llc_cart, &
265  source%id, source%level, source%llc_cart, &
266  source%distance(target)/(0.5d0*(target%size+source%size))
267  do iv=1,target%nv
268  !---------------------------------------------------------------------------
269  ! Target lower and upper corner in real space; these are the positions of
270  ! the first and last target points.
271  !---------------------------------------------------------------------------
272  tlc = target%myposition (target%mesh%lb, iv)
273  tuc = target%myposition (target%mesh%ub, iv)
274  !---------------------------------------------------------------------------
275  ! Map to source index space and limit to source interior. These are the last
276  ! source points to the left of a target point, with the next source point to
277  ! the right of the target point, and still inside the source interior.
278  !---------------------------------------------------------------------------
279  sl = source%index_only (tlc, 1, nowrap=.true.)
280  su = source%index_only (tuc, 1, nowrap=.true.)
281  sl = min(max(sl,source%mesh%li),source%mesh%ui-1)
282  su = min(max(su,source%mesh%li),source%mesh%ui-1)
283  !---------------------------------------------------------------------------
284  ! Map those source points back to the target index space, rounding up to
285  ! find the target points between the source points.
286  !---------------------------------------------------------------------------
287  slc = source%myposition(sl, iv)
288  suc = source%myposition(su, iv)
289  !---------------------------------------------------------------------------
290  ! Compensate for periodic wrapping of corner points
291  !---------------------------------------------------------------------------
292  !where (target%position < source%position .and. slc > source%position)
293  !slc = slc - target%mesh%b
294  !end where
295  !where (target%position > source%position .and. suc < source%position)
296  !suc = suc + target%mesh%b
297  !end where
298  tl = target%index_only (slc, 1, roundup=.true., nowrap=.true.)
299  tu = target%index_only (suc, 1, roundup=.true., nowrap=.true.)
300  !---------------------------------------------------------------------------
301  ! The brackets below should not really be needed, since the points should
302  ! all be inside the target legal range
303  !---------------------------------------------------------------------------
304  tl = min(max(tl,target%mesh%lb),target%mesh%ub)
305  tu = min(max(tu,target%mesh%lb),target%mesh%ub)
306  if (verbose > 0) then
307  write(io_unit%log,1) &
308  'guard_zones_t%restrict_vect: iv, tlc, sl, slc, tl =', &
309  iv, tlc, sl, slc, tl
310  1 format(a,i3,2(2(2x,3f8.4,2x,3i4)))
311  write(io_unit%log,1) &
312  'guard_zones_t%restrict_vect: iv, tuc, su, suc, tu =', &
313  iv, tuc, su, suc, tu
314  end if
315  !---------------------------------------------------------------------------
316  ! Compute indices and interpolation weights for the first target point, and
317  ! save the weights for the 8 nearest neighbor points in source space
318  !---------------------------------------------------------------------------
319  call source%index_space (target%myposition(tl), 1, j0, p)
320  j = 0
321  do i3=0,1
322  p3 = merge(p(3), 1.0-p(3), i3==1)
323  do i2=0,1
324  p2 = merge(p(2), 1.0-p(2), i2==1)
325  do i1=0,1
326  p1 = merge(p(1), 1.0-p(1), i1==1)
327  j = j+1
328  w(j) = p1*p2*p3
329  end do
330  end do
331  end do
332  !---------------------------------------------------------------------------
333  if (verbose > 0) &
334  write(io_unit%log,'(a,i3,4(2x,3i4))') &
335  'guard_zones_t%restrict_vect: iv, tl, tu, j0 =', &
336  iv, tl, tu, j0
337  !---------------------------------------------------------------------------
338  ! Map target and source memory to contiguous 3-D arrays
339  !---------------------------------------------------------------------------
340  tmem => target%mem(:,:,:,iv,target%it,1)
341  td => target%mem(:,:,:,target%idx%d,target%it,1)
342  do k=1,2
343  smem => source%mem(:,:,:,iv,jt(k),1)
344  sd => source%mem(:,:,:,source%idx%d,jt(k),1)
345  call update (target, tmem, td, tl, tu, iv, source, smem, sd, j0, k, &
346  pt(k), w, verbose)
347  end do
348  end do
349  if (detailed_timer) then
350  used = wallclock()-used
351  write (stdout,'(a,6p,f6.3," mus/cell")') &
352  "guard_zone_t%restrict_vect:", used/product(tu-tl+1)
353  end if
354  call trace%end (itimer)
355 END SUBROUTINE restrict_vect
356 
357 !===============================================================================
358 !> Update interpolated guard zone values by looping over target indices,
359 !> incrementing source indices by two (FIXME: this assumes factor two AMR)
360 !> for each step in target indices.
361 !>
362 !> For a 16x16x16 test case, the expected source index ranges are:
363 !> case source target j0 l u
364 !> -----------------------------------
365 !> left: 4-9 1-3 4 1 3
366 !> right: 14-19 19-22 14 19 22
367 !> full: 4-19 4-11 4 4 11
368 !> full: 4-19 12-19 4 12 19
369 !===============================================================================
370 SUBROUTINE update (target, tmem, td, l, u, iv, source, smem, sd, j0, k, &
371  pt, w, verbose)
372  class(patch_t), pointer:: target, source
373  real, dimension(:,:,:), contiguous, intent(out):: tmem
374  real, dimension(:,:,:), contiguous, intent(in):: smem
375  real, dimension(:,:,:), contiguous, intent(inout):: td, sd
376  integer:: l(3), u(3), iv, j0(3), k, verbose
377  real:: pt, w(8)
378  !.............................................................................
379  real(8):: wc(4)
380  real, allocatable:: f(:,:), tot(:)
381  integer:: i, j, n, i1 , i2, i3, j1, j2, j3, jj(3)
382  !-----------------------------------------------------------------------------
383  call trace%begin ('guard_zones_t%update')
384  if (detailed_timer) then
385  wc(1) = wallclock()
386  end if
387  n = product(u-l+1)
388  allocate (f(n,8), tot(n))
389  !-----------------------------------------------------------------------------
390  ! Gather
391  !-----------------------------------------------------------------------------
392  i = 0
393  j3 = j0(3)
394  do i3=l(3),u(3)
395  j2 = j0(2)
396  do i2=l(2),u(2)
397  j1 = j0(1)
398  if (target%pervolume(iv)) then
399  do i1=l(1),u(1)
400  i = i+1
401 !define VALIDATE
402 #ifdef VALIDATE
403  jj = [j1,j2,j3]
404  if (any(jj < source%li .or. jj > source%ui)) then
405  write(stderr ,1) iv, jj, j0, l, u
406  write(io_unit%log,1) iv, jj, j0, l, u
407  1 format("guard_zones_t%update: ERROR iv,jj,j0,l,u =",i3,4(2x,3i4))
408  end if
409 #endif
410  f(i,1) = smem(j1 ,j2 ,j3 )/sd(j1 ,j2 ,j3 )
411  f(i,2) = smem(j1+1,j2 ,j3 )/sd(j1+1,j2 ,j3 )
412  f(i,3) = smem(j1 ,j2+1,j3 )/sd(j1 ,j2+1,j3 )
413  f(i,4) = smem(j1+1,j2+1,j3 )/sd(j1+1,j2+1,j3 )
414  f(i,5) = smem(j1 ,j2 ,j3+1)/sd(j1 ,j2 ,j3+1)
415  f(i,6) = smem(j1+1,j2 ,j3+1)/sd(j1+1,j2 ,j3+1)
416  f(i,7) = smem(j1 ,j2+1,j3+1)/sd(j1 ,j2+1,j3+1)
417  f(i,8) = smem(j1+1,j2+1,j3+1)/sd(j1+1,j2+1,j3+1)
418  j1 = j1+2
419  end do
420  else
421  do i1=l(1),u(1)
422  i = i+1
423 #ifdef VALIDATE
424  jj = [j1,j2,j3]
425  if (any(jj < source%li .or. jj > source%ui)) then
426  write(stderr ,1) iv, jj, j0, l, u
427  write(io_unit%log,1) iv, jj, j0, l, u
428  end if
429 #endif
430  f(i,1) = smem(j1 ,j2 ,j3 )
431  f(i,2) = smem(j1+1,j2 ,j3 )
432  f(i,3) = smem(j1 ,j2+1,j3 )
433  f(i,4) = smem(j1+1,j2+1,j3 )
434  f(i,5) = smem(j1 ,j2 ,j3+1)
435  f(i,6) = smem(j1+1,j2 ,j3+1)
436  f(i,7) = smem(j1 ,j2+1,j3+1)
437  f(i,8) = smem(j1+1,j2+1,j3+1)
438  j1 = j1+2
439  end do
440  end if
441  j2 = j2+2
442  end do
443  j3 = j3+2
444  end do
445  if (detailed_timer) then
446  wc(2) = wallclock()
447  end if
448  !-----------------------------------------------------------------------------
449  ! Compute weighted results, vectorized
450  !-----------------------------------------------------------------------------
451  tot = 0.0
452  if (target%unsigned(iv)) then
453  do j=1,8
454  do i=1,n
455  tot(i) = tot(i) + pt*w(j)*log(f(i,j))
456  end do
457  end do
458  tot = exp(tot)
459  else
460  do j=1,8
461  do i=1,n
462  tot(i) = tot(i) + pt*w(j)*f(i,j)
463  end do
464  end do
465  end if
466  if (detailed_timer) then
467  wc(3) = wallclock()
468  end if
469  !-----------------------------------------------------------------------------
470  ! Scatter, vectorized
471  !-----------------------------------------------------------------------------
472  i = 0
473  do i3=l(3),u(3)
474  do i2=l(2),u(2)
475  if (k==1) then
476  do i1=l(1),u(1)
477  i = i+1
478  tmem(i1,i2,i3) = tot(i)
479  end do
480  else
481  do i1=l(1),u(1)
482  i = i+1
483  tmem(i1,i2,i3) = tmem(i1,i2,i3) + tot(i)
484  end do
485  end if
486  end do
487  end do
488  !-----------------------------------------------------------------------------
489  ! For per-volume quantities, multiply with mass density
490  !-----------------------------------------------------------------------------
491  if (target%pervolume(iv) .and. k==2) then
492  do i3=l(3),u(3)
493  do i2=l(2),u(2)
494  do i1=l(1),u(1)
495  tmem(i1,i2,i3) = tmem(i1,i2,i3)*td(i1,i2,i3)
496  end do
497  end do
498  end do
499  end if
500  deallocate (f, tot)
501  if (detailed_timer) then
502  wc(4) = wallclock()
503  do i=1,3
504  !$omp atomic update
505  update_t(i) = update_t(i) + (wc(i+1)-wc(i))
506  end do
507  end if
508  call trace%end()
509 END SUBROUTINE update
510 
511 !===============================================================================
512 !> Interpolate in space and time, for a case where it is safe = the source point
513 !> is inside the valid domain (inside a distance of half a cell from the source
514 !> patch boundary for non-staggered points)
515 !===============================================================================
516 SUBROUTINE interpolate (source, target, ii, iv, jt, pt, pi, set)
517  class(patch_t) :: source
518  class(patch_t), pointer :: target
519  integer :: ii(3), iv, jt(2)
520  real :: pt(2)
521  real, optional :: pi(:)
522  logical, optional :: set
523  !.............................................................................
524  real(8) :: pos(3)
525  integer :: j(3)
526  real :: p(3)
527  integer, save :: itimer=0
528  !-----------------------------------------------------------------------------
529  if (present(set)) then
530  pi(01) = pt(1)*(1.-p(3))*(1.-p(2))*(1.0-p(1))
531  pi(02) = pt(1)*(1.-p(3))*(1.-p(2))*( p(1))
532  pi(03) = pt(1)*(1.-p(3))*( p(2))*(1.0-p(1))
533  pi(04) = pt(1)*(1.-p(3))*( p(2))*( p(1))
534  pi(05) = pt(1)*( p(3))*(1.-p(2))*(1.0-p(1))
535  pi(06) = pt(1)*( p(3))*(1.-p(2))*( p(1))
536  pi(07) = pt(1)*( p(3))*( p(2))*(1.0-p(1))
537  pi(08) = pt(1)*( p(3))*( p(2))*( p(1))
538  pi(09) = pt(2)*(1.-p(3))*(1.-p(2))*(1.0-p(1))
539  pi(10) = pt(2)*(1.-p(3))*(1.-p(2))*( p(1))
540  pi(11) = pt(2)*(1.-p(3))*( p(2))*(1.0-p(1))
541  pi(12) = pt(2)*(1.-p(3))*( p(2))*( p(1))
542  pi(13) = pt(2)*( p(3))*(1.-p(2))*(1.0-p(1))
543  pi(14) = pt(2)*( p(3))*(1.-p(2))*( p(1))
544  pi(15) = pt(2)*( p(3))*( p(2))*(1.0-p(1))
545  pi(16) = pt(2)*( p(3))*( p(2))*( p(1))
546  else if (source%pervolume(iv)) then
547  call interpolate_specific (source, target, ii, iv, jt, pt)
548  else
549  pos = target%myposition (ii, iv)
550  call source%index_space (pos, iv, j, p)
551  !if (io%verbose > 1) &
552  ! write(io_unit%log,'(a,2i6,2x,3i4,2x,3f7.3)') &
553  ! 'interpolate: target, source, j, p =', target%id, source%id, j, p
554  !call trace%begin ('patch_t%interpolate', itimer=itimer)
555  if (source%unsigned(iv)) then
556  if (present(pi)) then
557  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = exp( &
558  pi(01)*log(source%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1)) + &
559  pi(02)*log(source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1)) + &
560  pi(03)*log(source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1)) + &
561  pi(04)*log(source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1)) + &
562  pi(05)*log(source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1)) + &
563  pi(06)*log(source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1)) + &
564  pi(07)*log(source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1)) + &
565  pi(08)*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1)) + &
566  pi(09)*log(source%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1)) + &
567  pi(10)*log(source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1)) + &
568  pi(11)*log(source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1)) + &
569  pi(12)*log(source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1)) + &
570  pi(13)*log(source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1)) + &
571  pi(14)*log(source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1)) + &
572  pi(15)*log(source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1)) + &
573  pi(16)*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1)))
574  else
575  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = exp( &
576  pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1)) + &
577  ( p(1))*log(source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1))) + &
578  ( p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1)) + &
579  ( p(1))*log(source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1)))) + &
580  ( p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1)) + &
581  ( p(1))*log(source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1))) + &
582  ( p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1)) + &
583  ( p(1))*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1))))) + &
584  pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1)) + &
585  ( p(1))*log(source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1))) + &
586  ( p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1)) + &
587  ( p(1))*log(source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1)))) + &
588  ( p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1)) + &
589  ( p(1))*log(source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1))) + &
590  ( p(2))*((1.0-p(1))*log(source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1)) + &
591  ( p(1))*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))))
592  end if
593  else
594  if (present(pi)) then
595  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
596  pi(01)*source%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1) + &
597  pi(02)*source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1) + &
598  pi(03)*source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1) + &
599  pi(04)*source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1) + &
600  pi(05)*source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1) + &
601  pi(06)*source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1) + &
602  pi(07)*source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1) + &
603  pi(08)*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1) + &
604  pi(09)*source%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1) + &
605  pi(10)*source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1) + &
606  pi(11)*source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1) + &
607  pi(12)*source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1) + &
608  pi(13)*source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1) + &
609  pi(14)*source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1) + &
610  pi(15)*source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1) + &
611  pi(16)*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1)
612  else
613  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
614  pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1) + &
615  ( p(1))*source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1)) + &
616  ( p(2))*((1.0-p(1))*source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1) + &
617  ( p(1))*source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1))) + &
618  ( p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1) + &
619  ( p(1))*source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1)) + &
620  ( p(2))*((1.0-p(1))*source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1) + &
621  ( p(1))*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1)))) + &
622  pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1) + &
623  ( p(1))*source%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1)) + &
624  ( p(2))*((1.0-p(1))*source%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1) + &
625  ( p(1))*source%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1))) + &
626  ( p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1) + &
627  ( p(1))*source%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1)) + &
628  ( p(2))*((1.0-p(1))*source%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1) + &
629  ( p(1))*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))
630  end if
631  end if
632  end if
633  if (target%id==io%id_debug .and. io%verbose>2) &
634  write(io%output,1) target%id, ii, iv, source%id, j, jt, p, pt, &
635  target%mem(ii(1),ii(2),ii(3),iv,target%it,1)
636  1 format(i6,1x,3i3,i4,2x,"DBG patch_t%interpolate src:",i6, &
637  " j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2," out:",1p,e15.6)
638  !call trace%end (itimer)
639 END SUBROUTINE interpolate
640 
641 !===============================================================================
642 !> Interpolate in space and time, for a case where it is safe = the source point
643 !> is inside the valid domain (inside a distance of half a cell from the source
644 !> patch boundary for non-staggered points)
645 !===============================================================================
646 SUBROUTINE interpolate_specific (source, target, ii, iv, jt, pt)
647  class(patch_t) :: source
648  class(patch_t), pointer :: target
649  integer :: ii(3), iv, jt(2)
650  real :: pt(2)
651  real :: tmp(2,2,2,2)
652  !.............................................................................
653  real(8) :: pos(3)
654  integer :: j(3)
655  real :: p(3)
656  integer, save :: itimer=0
657  !-----------------------------------------------------------------------------
658  pos = target%myposition (ii, iv)
659  call source%index_space (pos, iv, j, p)
660  !call trace%begin ('patch_t%interpolate', itimer=itimer)
661  !-----------------------------------------------------------------------------
662  ! Pick up 16 values for tri-linear interpolation, divided by density
663  !-----------------------------------------------------------------------------
664  associate(id=>source%idx%d)
665  tmp(:,:,:,:) = source%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,iv,jt(1:2),1) &
666  / source%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,id,jt(1:2),1)
667  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
668  pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,1) + &
669  ( p(1))*tmp(2,1,1,1)) + &
670  ( p(2))*((1.0-p(1))*tmp(1,2,1,1) + &
671  ( p(1))*tmp(2,2,1,1))) + &
672  ( p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,1) + &
673  ( p(1))*tmp(2,1,2,1)) + &
674  ( p(2))*((1.0-p(1))*tmp(1,2,2,1) + &
675  ( p(1))*tmp(2,2,2,1)))) + &
676  pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,2) + &
677  ( p(1))*tmp(2,1,1,2)) + &
678  ( p(2))*((1.0-p(1))*tmp(1,2,1,2) + &
679  ( p(1))*tmp(2,2,1,2))) + &
680  ( p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,2) + &
681  ( p(1))*tmp(2,1,2,2)) + &
682  ( p(2))*((1.0-p(1))*tmp(1,2,2,2) + &
683  ( p(1))*tmp(2,2,2,2))))
684  !-----------------------------------------------------------------------------
685  ! Multiply by target denssity
686  !-----------------------------------------------------------------------------
687  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
688  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) * &
689  target%mem(ii(1),ii(2),ii(3),id,target%it,1)
690  end associate
691  if (target%id==io%id_debug .and. io%verbose>2) &
692  write(io%output,1) target%id, ii, iv, source%id, j, jt, p, pt, &
693  target%mem(ii(1),ii(2),ii(3),iv,target%it,1)
694  1 format(i6,1x,3i3,i4,2x,"DBG patch_t%interpolate_specific src:", &
695  i6," j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2," out:",1p,e15.6)
696  !call trace%end (itimer)
697 END SUBROUTINE interpolate_specific
698 
699 SUBROUTINE info
700  if (verbose > 0) then
701  write (stdout,*) 'guard_zones_t%info: update-t =', update_t
702  end if
703 END SUBROUTINE info
704 
705 END MODULE guard_zones_mod
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
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
Optimized restrict operations, intended for guard zones (not conservative)