DISPATCH
remesh_mod.f90
1 !===============================================================================
2 !> This module handles remeshing, eg from lower resolution to higher resolution,
3 !> with special attention to maximizing speed, by identifying cases where cells
4 !> are aligned (static meshes). FIXME: initially assumed to always be the case.
5 !===============================================================================
6 MODULE remesh_mod
7  USE io_mod
8  USE trace_mod
9  USE omp_timer_mod
10  USE patch_mod
11  implicit none
12  private
13  type, public:: remesh_t
14  contains
15  procedure:: init
16  procedure:: prolong
17  procedure:: scatter_0th
18  procedure:: scatter_1st
19  procedure, nopass:: weights
20  end type
21  integer:: verbose=0
22  integer:: order=0
23  type(remesh_t), public:: remesh
24 CONTAINS
25 
26 !===============================================================================
27 !> Initialize
28 !===============================================================================
29 SUBROUTINE init (self)
30  class(remesh_t):: self
31  !.............................................................................
32  integer:: iostat
33  logical, save:: first_time=.true.
34  namelist /remesh_params/ verbose, order
35  !-----------------------------------------------------------------------------
36  call trace%begin ('remesh_t%init')
37  if (first_time) then
38  !$omp critical (input_cr)
39  if (first_time) then
40  rewind(io_unit%input)
41  read (io_unit%input, remesh_params, iostat=iostat)
42  write (io_unit%output, remesh_params)
43  first_time = .false.
44  end if
45  !$omp end critical (input_cr)
46  end if
47  call trace%end ()
48 END SUBROUTINE init
49 
50 !===============================================================================
51 !> Interface btw prolong and scatter. In the basic case with no_mans_land=T,
52 !> the lower left corner in source space is the index point to the left of the
53 !> target index point; i.e., the source%index_only point corresponding to
54 !> target%myposition (target%mesh%li).
55 !===============================================================================
56 SUBROUTINE prolong (self, target, source, jt, pt, all_cells)
57  class(remesh_t):: self
58  class(patch_t), pointer:: target, source
59  integer:: jt(2)
60  real:: pt(2)
61  logical:: all_cells
62  !.............................................................................
63  real(8):: pnt(3)
64  integer:: ll(3), nn(3)
65  !-----------------------------------------------------------------------------
66  call trace%begin ('remesh_t%prolong')
67  call self%init
68  if (all_cells) then
69  pnt = target%myposition (target%mesh%li)
70  if (order==0) then
71  ll = source%index_only (pnt, roundup=.true.)
72  nn = source%mesh%n/2
73  call self%scatter_0th (source, target, ll, nn)
74  else
75  ll = source%index_only (pnt)
76  nn = source%mesh%n/2+1
77  call self%scatter_1st (source, target, ll, nn)
78  end if
79  if (verbose > 0) then
80  write (io_unit%output,'(a,2i6,3f10.6,2(2x,3i4))') &
81  'remesh_t%prolong: source, target, pnt, ll, nn =', &
82  source%id, target%id, pnt, ll, nn
83  end if
84  else
85  call io%abort('remesh_t%prolong: cannot handle all_cells=.false.')
86  end if
87  call trace%end ()
88 END SUBROUTINE prolong
89 
90 !===============================================================================
91 !> Scatter = remesh from coarse to fine mesh. The index corner ll gives the
92 !> lower left corner in index space. If this is, for example, mesh%li, the
93 !> region covered will start at li1 for the coarse mesh, and at li for the
94 !> fine mesh.
95 !===============================================================================
96 SUBROUTINE scatter_0th (self, coarse, fine, ll, nn)
97  class(remesh_t):: self
98  class(patch_t), pointer:: coarse, fine
99  integer:: ll(3), nn(3)
100  optional:: nn
101  !.............................................................................
102  integer:: i1, i2, i3, j1, j2, j3, iv, lc(3), uc(3), lf(3), uf(3)
103  real, dimension(:,:,:), pointer:: c, f
104  real:: s(2)
105  integer, save:: itimer=0
106  !-----------------------------------------------------------------------------
107  call trace%begin ('remesh_t%scatter_0th', itimer=itimer)
108  lc = ll
109  if (present(nn)) then
110  uc = ll + nn - 1
111  else
112  uc = ll + coarse%mesh%n/2 - 1
113  end if
114  lf = fine%mesh%li
115  uf = fine%mesh%ui
116  !-----------------------------------------------------------------------------
117  ! Loop over variables, duplicating the coarse value in all three directions
118  !-----------------------------------------------------------------------------
119  do iv=1,fine%nv
120  c => coarse%mem(:,:,:,iv,coarse%it,1)
121  f => fine%mem(:,:,:,iv, fine%it,1)
122  j3 = lf(3)
123  do i3=lc(3),uc(3)
124  j2 = lf(2)
125  do i2=lc(2),uc(2)
126  !-----------------------------------------------------------------------
127  !-----------------------------------------------------------------------
128  j1 = lf(1)
129  do i1=lc(1),uc(1)
130  f(j1 ,j2 ,j3 ) = c(i1,i2,i3)
131  f(j1+1,j2 ,j3 ) = c(i1,i2,i3)
132  f(j1 ,j2+1,j3 ) = c(i1,i2,i3)
133  f(j1+1,j2+1,j3 ) = c(i1,i2,i3)
134  f(j1 ,j2 ,j3+1) = c(i1,i2,i3)
135  f(j1+1,j2 ,j3+1) = c(i1,i2,i3)
136  f(j1 ,j2+1,j3+1) = c(i1,i2,i3)
137  f(j1+1,j2+1,j3+1) = c(i1,i2,i3)
138  j1 = j1+2
139  end do
140  j2 = j2+2
141  end do
142  j3 = j3+2
143  end do
144  if (verbose > 1) then
145  s(1) = sum(c(lc(1):uc(1),lc(2):uc(2),lc(3):uc(3)))*product(coarse%ds)
146  s(2) = sum(f(lf(1):uf(1),lf(2):uf(2),lf(3):uf(3)))*product( fine%ds)
147  write (io_unit%output,'(a,i4,1p,3e12.3,4(2x,3i4))') &
148  'remesh_t%scatter_0th: iv, sum_c, sum_f =', &
149  iv, s, fine%fmaxval(f), lc, uc, lf, uf
150  end if
151  end do
152  call trace%end (itimer)
153 END SUBROUTINE scatter_0th
154 
155 !===============================================================================
156 !> Scatter = remesh from coarse to fine mesh. The index corner ll gives the
157 !> lower left corner in index space. If this is, for example, mesh%li, the
158 !> region covered will start at li-1 for the coarse mesh, and at li-2 for the
159 !> fine mesh. The li-2 point will lack contributions, so the the fine mesh will
160 !> only be valid from lo=li-1 to uo=ui+1.
161 !===============================================================================
162 SUBROUTINE scatter_1st (self, coarse, fine, ll, nn)
163  class(remesh_t):: self
164  class(patch_t), pointer:: coarse, fine
165  integer:: ll(3), nn(3)
166  optional:: nn
167  !.............................................................................
168  integer:: i1, i2, i3, j1, j2, j3, iv, j, lc(3), uc(3), lf(3), ic(3)
169  real, dimension(:,:,:), pointer:: c, f
170  real, allocatable:: out(:,:), w(:,:)
171  integer, save:: itimer=0
172  !-----------------------------------------------------------------------------
173  call trace%begin ('remesh_t%scatter_1st', itimer=itimer)
174  ic = coarse%mesh%li
175  lc = ll
176  if (present(nn)) then
177  uc = ll + nn
178  else
179  uc = ll + coarse%mesh%n/2 - 1
180  end if
181  lf = fine%mesh%li
182  allocate (out(fine%gn(1),8), w(8,fine%gn(1)))
183  !-----------------------------------------------------------------------------
184  ! Loop over variables, computing remeshing weights for each type of centering
185  !-----------------------------------------------------------------------------
186  do iv=1,fine%nv
187  c => coarse%mem(:,:,:,iv,coarse%it,1)
188  f => fine%mem(:,:,:,iv, fine%it,1)
189  w = weights(iv)
190  f = 0.0
191  do i3=lc(3)-1,uc(3)+1
192  do i2=lc(2)-1,uc(2)+1
193  !-------------------------------------------------------------------------
194  ! Loop over the 8 fine points for each coarse point, and over coarse points
195  ! in the x-direction (this loop vectorizes)
196  !-------------------------------------------------------------------------
197  do j=1,8
198  do i1=lc(1)-1,uc(1)+1
199  out(i1,j) = w(1,j)*c(i1,i2 ,i3 ) + w(2,j)*c(i1+1,i2 ,i3 ) &
200  + w(3,j)*c(i1,i2+1,i3 ) + w(4,j)*c(i1+1,i2+1,i3 ) &
201  + w(5,j)*c(i1,i2 ,i3+1) + w(6,j)*c(i1+1,i2 ,i3+1) &
202  + w(7,j)*c(i1,i2+1,i3+1) + w(8,j)*c(i1+1,i2+1,i3+1)
203  end do
204  end do
205  !-------------------------------------------------------------------------
206  ! Add up the contributions to the 8 fine points for each coarse point.
207  ! This loop does not vectorize, but has a significant load per iteration
208  !-------------------------------------------------------------------------
209  i1 = lc(1)-1
210  j1 = lf(1) + (i1-ic(1))*2
211  j2 = lf(2) + (i2-ic(2))*2
212  j3 = lf(3) + (i3-ic(3))*2
213  do i1=lc(1)-1,uc(1)+2
214  !if (i1 >= lc(1) .and. i2 >= lc(2) .and. i3 >= lc(3)) &
215  f(j1 ,j2 ,j3 ) = f(j1 ,j2 ,j3 ) + out(i1,1)
216  !if (i1 <= uc(2) .and. i2 >= lc(2) .and. i3 >= lc(3)) &
217  f(j1+1,j2 ,j3 ) = f(j1+1,j2 ,j3 ) + out(i1,2)
218  !if (i1 >= lc(1) .and. i2 <= uc(2) .and. i3 >= lc(3)) &
219  f(j1 ,j2+1,j3 ) = f(j1 ,j2+1,j3 ) + out(i1,3)
220  !if (i1 <= uc(1) .and. i2 <= uc(2) .and. i3 >= lc(3)) &
221  f(j1+1,j2+1,j3 ) = f(j1+1,j2+1,j3 ) + out(i1,4)
222  !if (i1 >= lc(1) .and. i2 >= lc(2) .and. i3 <= uc(3)) &
223  f(j1 ,j2 ,j3+1) = f(j1 ,j2 ,j3+1) + out(i1,5)
224  !if (i1 <= uc(1) .and. i2 >= lc(2) .and. i3 <= uc(3)) &
225  f(j1+1,j2 ,j3+1) = f(j1+1,j2 ,j3+1) + out(i1,6)
226  !if (i1 >= lc(1) .and. i2 <= uc(2) .and. i3 <= uc(3)) &
227  f(j1 ,j2+1,j3+1) = f(j1 ,j2+1,j3+1) + out(i1,7)
228  !if (i1 <= uc(1) .and. i2 <= uc(2) .and. i3 <= uc(3)) &
229  f(j1+1,j2+1,j3+1) = f(j1+1,j2+1,j3+1) + out(i1,8)
230  j1 = j1+2
231  end do
232  end do
233  end do
234  end do
235  deallocate (out)
236  call trace%end (itimer)
237 END SUBROUTINE scatter_1st
238 
239 !===============================================================================
240 !> Weight factors for RAMSES type arrangement of coarse / fine (no_mans_land=t).
241 !> The weight w(i,j) is the contribution of coarse point i to the value at fine
242 !> point j. Both i and j are numbered in increasing order of first x, then y,
243 !> then z.
244 !>
245 !> FIXME: There should be a solver-specific weight function defined for each
246 !> solver -- this one is only intended for testing.
247 !===============================================================================
248 FUNCTION weights (iv) RESULT (w)
249  integer:: iv
250  real:: w(8,8), w1, w2, w3
251  integer:: i1, i2, i3, i, j1, j2, j3, j
252  !-----------------------------------------------------------------------------
253  w = 0.0
254  i = 1
255  do i3=0,1
256  do i2=0,1
257  do i1=0,1
258  j = 1
259  do j3=0,1
260  w3 = merge(0.75,0.25,i3==j3)
261  do j2=0,1
262  w2 = merge(0.75,0.25,i2==j2)
263  !-------------------------------------------------------------------
264  ! i1 is the coarse point location in x: 0 for 0.0, 1 for 1.0
265  ! j1 is the fine point location in x: 0 for .25, 1 for .75
266  !-------------------------------------------------------------------
267  do j1=0,1
268  w1 = merge(0.75,0.25,i1==j1)
269  w(i,j) = w1*w2*w3
270  j = j+1
271  end do
272  end do
273  end do
274  i = i+1
275  end do
276  end do
277  end do
278 END FUNCTION weights
279 
280 END MODULE remesh_mod
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
This module handles remeshing, eg from lower resolution to higher resolution, with special attention ...
Definition: remesh_mod.f90:6
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Definition: io_mod.f90:4