MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_functions_forest.t
Go to the documentation of this file.
2 
3  implicit none
4  private
5 
6  public :: write_forest
7  public :: read_forest
8  public :: coarsen_tree_leaf
9  public :: refine_tree_leaf
10  public :: change_ipe_tree_leaf
11  public :: init_forest_root
12 
13 
14 contains
15 
16 
17  !> build root forest
18  subroutine init_forest_root
19  use mod_forest
22  use mod_amr_solution_node, only: getnode
24 
25  integer :: ig^d, level, igrid, ipe
26  integer :: iside, i^d, morton_no, isfc
27 
28 
29  level=1
30  morton_no=0
31  ipe=0
32  nparents=0
33  nleafs={ng^d(1)*}
35  nleafs_level(1)={ng^d(1)*}
37  call get_morton_range
38  ! Generate Morton-order space-filling curve to connect level 1 blocks
40  do isfc=1,nglev1
41  ig^d=sfc_iglevel1(^d,isfc)\
42  morton_no=morton_no+1
43  if (morton_no>morton_stop(ipe)) ipe=ipe+1
44  igrid=getnode(ipe)
45  if (ipe==mype) then
46  sfc_to_igrid(morton_no)=igrid
47  igrid_to_sfc(igrid)=morton_no
48  end if
49  call init_tree_leaf(tree_root(ig^d),ig^d,level,igrid,ipe,.true.)
50  end do
51 
52  ! update root neighbor
53  {do ig^db=1,ng^db(1)\}
54  {do iside=1,2
55  i^dd=kr(^dd,^d)*(2*iside-3);
56  call find_root_neighbor(tree_root(ig^dd)%node%neighbor(iside,^d), &
57  tree_root(ig^dd),i^dd)
58  end do\}
59  {end do\}
60 
61  ! This call is here to ensure the sfc array is initialized
62  call amr_morton_order()
63 
64  end subroutine init_forest_root
65 
66  subroutine init_tree_leaf(tree,ig^D,level,igrid,ipe,active)
67  use mod_forest
68  implicit none
69 
70  type(tree_node_ptr) :: tree
71  integer, intent(in) :: ig^D, level, igrid, ipe
72  logical, intent(in) :: active
73  integer :: ic^D
74 
75  allocate(tree%node)
76 
77  tree%node%ig^d=ig^d;
78  tree%node%level=level
79  tree%node%igrid=igrid
80  tree%node%ipe=ipe
81 
82  tree%node%leaf=.true.
83  tree%node%active=active
84 
85  nullify(tree%node%parent%node)
86  {do ic^db=1,2\}
87  nullify(tree%node%child(ic^d)%node)
88  {end do\}
89 
90  call add_to_linked_list(level,tree)
91 
92  ! initialize neighbor pointers
93  nullify({tree%node%neighbor(1,^d)%node},{tree%node%neighbor(2,^d)%node})
94 
95  igrid_to_node(igrid,ipe)%node => tree%node
96 
97  end subroutine init_tree_leaf
98 
99  subroutine coarsen_tree_leaf(igrid,ipe,child_igrid,child_ipe,active)
100  use mod_forest
101  implicit none
102 
103  integer, intent(in) :: igrid, ipe
104  integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
105  logical, intent(out) :: active
106 
107 
108  integer :: level, ic^d, child_level, iside, iotherside, vote
109  type(tree_node_ptr) :: tree, child, child_neighbor
110 
111  tree%node => igrid_to_node(child_igrid(1^d&),child_ipe(1^d&))%node%parent%node
112  level=tree%node%level
113 
114  call add_to_linked_list(level,tree)
115 
116  child_level=level+1
117  vote=0
118 
119  {do ic^db=1,2\}
120  child%node => tree%node%child(ic^d)%node
121 
122  ! vote for active:
123  if(child%node%active) vote=vote+1
124 
125  call delete_from_linked_list(child_level,child)
126 
127  ! update neighbor pointers
128  {iside=ic^d
129  child_neighbor%node => child%node%neighbor(iside,^d)%node
130  if (associated(child_neighbor%node)) then
131  if (child%node%ig^d==child_neighbor%node%ig^d) then ! pole
132  nullify(child_neighbor%node%neighbor(iside,^d)%node)
133  else
134  iotherside=3-iside
135  nullify(child_neighbor%node%neighbor(iotherside,^d)%node)
136  end if
137  end if\}
138 
139  nullify(tree%node%child(ic^d)%node)
140  deallocate(igrid_to_node(child_igrid(ic^d),child_ipe(ic^d))%node)
141  {end do\}
142 
143  tree%node%leaf=.true.
144  tree%node%igrid=igrid
145  tree%node%ipe=ipe
146  igrid_to_node(igrid,ipe)%node => tree%node
147 
148  ! Count the vote and set active/passive state:
149 
150  if (vote /= 2**^nd) then
151  !if (vote == 0) then
152  tree%node%active = .false.
153  nleafs_active = nleafs_active - vote
154  else
155  tree%node%active = .true.
156  nleafs_active = nleafs_active - vote + 1
157  end if
158  active = tree%node%active
159 
160  nleafs=nleafs-2**^nd+1
161  nparents=nparents-1
162  nleafs_level(child_level)=nleafs_level(child_level)-2**^nd
163  nleafs_level(level)=nleafs_level(level)+1
164 
165  end subroutine coarsen_tree_leaf
166 
167  subroutine refine_tree_leaf(child_igrid,child_ipe,igrid,ipe,active)
168  use mod_forest
170  use mod_amr_neighbors, only: find_neighbor
171 
172  integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
173  integer, intent(in) :: igrid, ipe
174  logical, intent(out):: active
175 
176  integer :: ig^d, level, i^d, ic^d, child_ig^d, child_level, iside
177  integer :: my_neighbor_type
178  logical, dimension(ndim) :: pole
179  type(tree_node_ptr) :: tree, child, my_neighbor
180 
181  tree%node => igrid_to_node(igrid,ipe)%node
182  ig^d=tree%node%ig^d;
183  level=tree%node%level
184  active=tree%node%active
185 
186  tree%node%ipe=-1
187  tree%node%igrid=0
188  tree%node%leaf=.false.
189  tree%node%active=.true.
190 
191  call delete_from_linked_list(level,tree)
192 
193  child_level=level+1
194 
195  {do ic^db=1,2\}
196  child_ig^d=2*(ig^d-1)+ic^d;
197  call init_tree_leaf(child,child_ig^d,child_level, &
198  child_igrid(ic^d),child_ipe(ic^d),active)
199 
200  igrid_to_node(child_igrid(ic^d),child_ipe(ic^d))%node => child%node
201 
202  tree%node%child(ic^d)%node => child%node
203  child%node%parent%node => tree%node
204  {end do\}
205 
206  ! update neighbor pointers
207  {do ic^db=1,2\}
208  child%node => tree%node%child(ic^d)%node
209  {iside=ic^d
210  i^dd=kr(^dd,^d)*(2*iside-3);
211  call find_neighbor(my_neighbor,my_neighbor_type,child,i^dd,pole)
212  select case (my_neighbor_type)
213  case (neighbor_sibling, neighbor_fine)
214  child%node%neighbor(iside,^d)%node => my_neighbor%node
215  if (pole(^d)) then
216  my_neighbor%node%neighbor(iside,^d)%node => child%node
217  else
218  my_neighbor%node%neighbor(3-iside,^d)%node => child%node
219  end if
220  case default
221  nullify(child%node%neighbor(iside,^d)%node)
222  end select
223  child%node%neighbor(3-ic^d,^d)%node=>tree%node%child(3-ic^d^d%ic^dd)%node\}
224  {end do\}
225 
226  nleafs=nleafs+2**^nd-1
227  nparents=nparents+1
228  nleafs_level(child_level)=nleafs_level(child_level)+2**^nd
229  nleafs_level(level)=nleafs_level(level)-1
230 
231  if (active) nleafs_active = nleafs_active + 2**^nd-1
232 
233  end subroutine refine_tree_leaf
234 
235  subroutine change_ipe_tree_leaf(recv_igrid,recv_ipe,send_igrid,send_ipe)
236  use mod_forest
237  implicit none
238 
239  integer, intent(in) :: recv_igrid, recv_ipe, send_igrid, send_ipe
240 
241  type(tree_node_ptr) :: tree
242 
243  tree%node => igrid_to_node(send_igrid,send_ipe)%node
244 
245  tree%node%igrid=recv_igrid
246  tree%node%ipe=recv_ipe
247 
248  nullify(igrid_to_node(send_igrid,send_ipe)%node)
249  igrid_to_node(recv_igrid,recv_ipe)%node => tree%node
250 
251  end subroutine change_ipe_tree_leaf
252 
253  subroutine add_to_linked_list(level,tree)
254  use mod_forest
255  implicit none
256 
257  integer, intent(in) :: level
258  type(tree_node_ptr) :: tree
259 
260  nullify(tree%node%next%node)
261  if (associated(level_head(level)%node)) then
262  tree%node%prev%node => level_tail(level)%node
263  level_tail(level)%node%next%node => tree%node
264  level_tail(level)%node => tree%node
265  else
266  level_head(level)%node => tree%node
267  level_tail(level)%node => tree%node
268  nullify(tree%node%prev%node)
269  end if
270 
271  end subroutine add_to_linked_list
272 
273  subroutine delete_from_linked_list(level,tree)
274  use mod_forest
275  implicit none
276 
277  integer, intent(in) :: level
278  type(tree_node_ptr) :: tree
279 
280  type(tree_node_ptr) :: next, prev
281 
282  prev%node => tree%node%prev%node
283  next%node => tree%node%next%node
284  if (associated(next%node).and.associated(prev%node)) then
285  prev%node%next%node => next%node
286  next%node%prev%node => prev%node
287  else if (associated(prev%node)) then
288  level_tail(level)%node => prev%node
289  nullify(prev%node%next%node)
290  else if (associated(next%node)) then
291  level_head(level)%node => next%node
292  nullify(next%node%prev%node)
293  else
294  nullify(level_head(level)%node)
295  nullify(level_tail(level)%node)
296  end if
297 
298  end subroutine delete_from_linked_list
299 
300  subroutine write_forest(file_handle)
301  use mod_forest
303 
304  integer, intent(in) :: file_handle
305 
306  integer, dimension(MPI_STATUS_SIZE) :: status
307  integer :: ig^d,isfc
308 
309  do isfc=1,nglev1
310  ig^d=sfc_iglevel1(^d,isfc)\
311  call write_node(tree_root(ig^d))
312  end do
313 
314  contains
315 
316  recursive subroutine write_node(tree)
317  implicit none
318 
319  type(tree_node_ptr) :: tree
320 
321  integer :: ic^d
322 
323  call mpi_file_write(file_handle,tree%node%leaf,1,mpi_logical,status,ierrmpi)
324 
325  if (.not.tree%node%leaf) then
326  {do ic^db=1,2\}
327  call write_node(tree%node%child(ic^d))
328  {end do\}
329  end if
330 
331  end subroutine write_node
332 
333  end subroutine write_forest
334 
335  subroutine read_forest(file_handle)
336  use mod_forest
339  use mod_amr_solution_node, only: getnode
342 
343  integer, intent(in) :: file_handle
344 
345  integer, dimension(MPI_STATUS_SIZE) :: status
346  !integer :: ig^D, level, size_logical, Morton_no, igrid, ipe
347  integer :: ig^d, level, morton_no, igrid, ipe, isfc
348 
349  morton_no=0
350  ipe=0
351  level=1
352  nleafs_level(1:nlevelshi) = 0
353  nparents = 0
354 
355  call get_morton_range
357  do isfc=1,nglev1
358  ig^d=sfc_iglevel1(^d,isfc)\
359  allocate(tree_root(ig^d)%node)
360  nullify(tree_root(ig^d)%node%parent%node)
361  call read_node(tree_root(ig^d),ig^d,level)
362  end do
363 
364  call get_level_range
365 
366  ! Rebuild tree connectivity
367  call getigrids
368  call build_connectivity
369 
370  ! This call is here to ensure the sfc array is initialized
371  call amr_morton_order()
372 
373  contains
374 
375  recursive subroutine read_node(tree,ig^D,level)
376  implicit none
377 
378  type(tree_node_ptr) :: tree
379  integer, intent(in) :: ig^d, level
380 
381  logical :: leaf
382  integer :: ic^d, child_ig^d, child_level
383 
384  if (mype==0) then
385  call mpi_file_read(file_handle,leaf,1,mpi_logical, &
386  status,ierrmpi)
387  end if
388  if (npe>1) call mpi_bcast(leaf,1,mpi_logical,0,icomm,ierrmpi)
389 
390  tree%node%leaf=leaf
391  tree%node%ig^d=ig^d;
392  tree%node%level=level
393  tree%node%active=.true. .and. leaf
394 
395  {do ic^db=1,2\}
396  nullify(tree%node%child(ic^d)%node)
397  {end do\}
398  nullify({tree%node%neighbor(1,^d)%node},{tree%node%neighbor(2,^d)%node})
399  nullify(tree%node%next%node,tree%node%prev%node)
400 
401  call asign_tree_neighbor(tree)
402 
403  if (leaf) then
404  call add_to_linked_list(level,tree)
405  nleafs_level(level) = nleafs_level(level) + 1
406 
407  morton_no=morton_no+1
408  if (morton_no>morton_stop(ipe)) ipe=ipe+1
409  igrid=getnode(ipe)
410  tree%node%igrid=igrid
411  tree%node%ipe=ipe
412  igrid_to_node(igrid,ipe)%node => tree%node
413  if (ipe==mype) sfc_to_igrid(morton_no)=igrid
414  else
415  nparents = nparents + 1
416  tree%node%igrid=0
417  tree%node%ipe=-1
418  child_level=level+1
419  {do ic^db=1,2\}
420  child_ig^d=2*(ig^d-1)+ic^d;
421  allocate(tree%node%child(ic^d)%node)
422  tree%node%child(ic^d)%node%parent%node => tree%node
423  call read_node(tree%node%child(ic^d),child_ig^d,child_level)
424  {end do\}
425  end if
426 
427  end subroutine read_node
428 
429  end subroutine read_forest
430 
431 
432 end module mod_functions_forest
recursive subroutine write_node(tree)
recursive subroutine read_node(tree, igD, level)
subroutine, public find_root_neighbor(tree_neighbor, tree, iD)
find neighors of level-one root blocks
subroutine, public asign_tree_neighbor(tree)
asign tree node neighor
subroutine, public find_neighbor(my_neighbor, my_neighbor_type, tree, iD, pole)
find neighors of all blocks
integer function, public getnode(ipe)
Get first available igrid on processor ipe.
Module with basic grid data structures.
Definition: mod_forest.t:2
integer, dimension(:), allocatable, save nleafs_level
How many leaves are present per refinement level.
Definition: mod_forest.t:81
integer, dimension(:), allocatable, save sfc_to_igrid
Go from a Morton number to an igrid index (for a single processor)
Definition: mod_forest.t:53
integer nleafs_active
Definition: mod_forest.t:78
integer, save nleafs
Number of leaf block.
Definition: mod_forest.t:76
integer, dimension(:), allocatable, save igrid_to_sfc
Go from a grid index to Morton number (for a single processor)
Definition: mod_forest.t:56
integer nglev1
Definition: mod_forest.t:78
integer, dimension(:,:), allocatable, save sfc_iglevel1
Space filling curve for level 1 grid. sfc_iglevel1(^D, MN) gives ig^D (the spatial index of the grid)
Definition: mod_forest.t:47
type(tree_node_ptr), dimension(:), allocatable, save level_tail
The tail pointer of the linked list per refinement level.
Definition: mod_forest.t:38
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition: mod_forest.t:65
type(tree_node_ptr), dimension(:^d &), allocatable, save tree_root
Pointers to the coarse grid.
Definition: mod_forest.t:29
integer, save nparents
Number of parent blocks.
Definition: mod_forest.t:73
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
type(tree_node_ptr), dimension(:), allocatable, save level_head
The head pointer of the linked list per refinement level.
Definition: mod_forest.t:35
subroutine, public change_ipe_tree_leaf(recv_igrid, recv_ipe, send_igrid, send_ipe)
subroutine delete_from_linked_list(level, tree)
subroutine, public write_forest(file_handle)
subroutine, public read_forest(file_handle)
subroutine, public refine_tree_leaf(child_igrid, child_ipe, igrid, ipe, active)
subroutine add_to_linked_list(level, tree)
subroutine, public coarsen_tree_leaf(igrid, ipe, child_igrid, child_ipe, active)
subroutine init_tree_leaf(tree, igD, level, igrid, ipe, active)
subroutine, public init_forest_root
build root forest
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter nlevelshi
The maximum number of levels in the grid refinement.
integer icomm
The MPI communicator.
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
subroutine get_morton_range
Set the Morton range for each processor.
subroutine amr_morton_order
Construct Morton-order as a global recursive lexicographic ordering.
subroutine level1_morton_order
build Morton space filling curve for level 1 grid
Pointer to a tree_node.
Definition: mod_forest.t:6