MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_space_filling_curve.t
Go to the documentation of this file.
2  implicit none
3 
4 contains
5 
6  !> build Morton space filling curve for level 1 grid
7  subroutine level1_morton_order
8  ! use Morton curve to connect level 1 grid blocks
9  use mod_forest
11 
12  integer, allocatable :: gsq_sfc(:^d&),seq_sfc(:),seq_ig^D(:)
13  integer :: ig^D, ngsq^D, isq, total_number
14  !integer(kind=8), external :: mortonEncode
15  logical, allocatable :: in_domain(:)
16 
17  ! use the smallest square/cube to cover the full domain
18  ngsq^d=2**ceiling(log(real(ng^d(1)))/log(2.0));
19  {^nooned
20  {ngsq^d=max(ngsq^dd) \}
21  }
22  total_number={ngsq^d|*}
23  ! Morton number acquired by block numbers
24  allocate(gsq_sfc(ngsq^d))
25  ! Morton number in sequence
26  allocate(seq_sfc(total_number))
27  ! block numbers in sequence
28  {allocate(seq_ig^d(total_number))\}
29  allocate(in_domain(total_number))
30  in_domain=.true.
31  ! get Morton-order numbers in the square/cube
32  {do ig^db=1,ngsq^db\}
33  gsq_sfc(ig^d)=int(mortonencode(ig^d-1,ndim))+1
34  seq_sfc(gsq_sfc(ig^d))=gsq_sfc(ig^d)
35  {seq_ig^d(gsq_sfc(ig^dd))=ig^d \}
36  {end do\}
37  ! mark blocks that are out of the domain and change Morton number
38  do isq=1,total_number
39  if (seq_ig^d(isq)>ng^d(1)|.or.) then
40  seq_sfc(isq:total_number)=seq_sfc(isq:total_number)-1
41  in_domain(isq)=.false.
42  end if
43  end do
44  ! copy the modified Morton numbers to the blocks in the domain
45  allocate(iglevel1_sfc(ng^d(1)))
46  allocate(sfc_iglevel1(ndim,nglev1))
47  do isq=1,total_number
48  if(in_domain(isq)) then
49  iglevel1_sfc(seq_ig^d(isq))=seq_sfc(isq)
50  {sfc_iglevel1(^d,seq_sfc(isq))=seq_ig^d(isq) \}
51  end if
52  end do
53 
54  deallocate(gsq_sfc,seq_sfc,seq_ig^d,in_domain)
55 
56  end subroutine level1_morton_order
57 
58  integer(kind=8) function mortonencode(ig^D,ndim)
59  use iso_fortran_env, only : int64
60  implicit none
61  integer(kind=4), intent(in) :: ig^D,ndim
62  integer(kind=4) :: i
63  integer(kind=8) :: answer, lg^D
64 
65  ! Create a 64-bit version of ig^D
66  lg^d=ig^d;
67  answer = 0
68 
69  do i=0,64/ndim
70  {^ifoned answer=ig1}
71  {^iftwod
72  answer=ior(answer,ior(ishft(iand(lg1,ishft(1_int64,i)),i),&
73  ishft(iand(lg2,ishft(1_int64,i)),i+1)))
74  \}
75  {^ifthreed
76  answer=ior(answer,ior(ishft(iand(lg1,ishft(1_int64,i)),2*i),ior(ishft(&
77  iand(lg2,ishft(1_int64,i)),2*i+1),ishft(iand(lg3,ishft(1_int64,i)),2*i+2))))
78  \}
79  end do
80  mortonencode=answer
81  return
82  end function mortonencode
83 
84  !> Construct Morton-order as a global recursive lexicographic ordering.
85  subroutine amr_morton_order
88 
89  integer :: ig^D, Morton_no, isfc
90 
91  morton_no=0
92  nglev1={ng^d(1)*}
93  do isfc=1,nglev1
94  ig^d=sfc_iglevel1(^d,isfc)\
95  call get_morton_number(tree_root(ig^d))
96  end do
97 
98  if (morton_no/=nleafs) then
99  call mpistop("error in amr_Morton_order: Morton_no/=nleafs")
100  end if
101 
102  contains
103 
104  recursive subroutine get_morton_number(tree)
106  type(tree_node_ptr) :: tree
107 
108  integer :: ic^D
109 
110  if (tree%node%leaf) then
111  morton_no=morton_no+1
112  sfc(1,morton_no)=tree%node%igrid
113  sfc(2,morton_no)=tree%node%ipe
114  if (tree%node%active) then
115  sfc(3,morton_no)=1
116  else
117  sfc(3,morton_no)=0
118  end if
119  if(tree%node%ipe==mype) igrid_to_sfc(tree%node%igrid)=morton_no
120  else
121  {do ic^db=1,2\}
122  call get_morton_number(tree%node%child(ic^d))
123  {end do\}
124  end if
125 
126  end subroutine get_morton_number
127 
128  end subroutine amr_morton_order
129 
130  !> Set the Morton range for each processor
131  subroutine get_morton_range
134 
135  integer :: ipe, blocks_left, procs_left, num_blocks
136 
137  if (allocated(sfc_to_igrid)) deallocate(sfc_to_igrid)
138  {#IFDEF EVOLVINGBOUNDARY
139  if (allocated(sfc_phybound)) deallocate(sfc_phybound)
140  }
141 
142  blocks_left = nleafs
143 
144  do ipe = 0, npe-1
145  if (ipe == 0) then
146  morton_start(ipe) = 1
147  else
148  morton_start(ipe) = morton_stop(ipe-1) + 1
149  end if
150 
151  ! Compute how many blocks this cpu should take
152  procs_left = npe - ipe
153  num_blocks = ceiling(blocks_left / dble(procs_left))
154  morton_stop(ipe) = morton_start(ipe) + num_blocks - 1
155  blocks_left = blocks_left - num_blocks
156  end do
157 
159  {#IFDEF EVOLVINGBOUNDARY
160  allocate(sfc_phybound(nleafs))
161  sfc_phybound=0
162  }
163 
164  end subroutine get_morton_range
165 
166  subroutine get_morton_range_active
169 
170  ! Cut the sfc based on weighted decision.
171  ! Oliver Porth, 02.02.2012
172 
173  !!!Here you choose the weithts:!!
174  integer, parameter :: wa=3, wp=1
175  ! wp : Weight for passive block
176  ! wa : Weight for active block
177  ! wp=0 : balance load (active blocks) exactly and
178  ! don't care about memory imbalance
179  ! wp=wa : balance memory exactly and don't care about load
180  ! Maximum possible memory imbalance is X=wa/wp.
181 
182  ! If you run into memory issues, decrease this ratio.
183  ! Scaling should be better if you allow for higher ratio.
184  ! Note also that passive cells still do regridding and boundary swap.
185  ! I have best results with a ratio 2:1, but it is problem dependent.
186  ! It can't do magic though...
187  ! Best to make sure that the sfc is properly aligned with your problem.
188 
189  integer :: ipe, Morton_no
190  integer :: Mtot, Mstop, Mcurr
191  ! For debugging:
192  integer :: nactive(0:npe-1),npassive(0:npe-1)
193  !double precision, save :: ptasum=0
194 
195  if (allocated(sfc_to_igrid)) deallocate(sfc_to_igrid)
196  {#IFDEF EVOLVINGBOUNDARY
197  if (allocated(sfc_phybound)) deallocate(sfc_phybound)
198  }
199 
200  mtot = nleafs_active*wa+(nleafs-nleafs_active)*wp
201  ipe = 0
202  mcurr = 0
203 
204  nactive=0
205  npassive=0
206 
207  morton_start(0) = 1
208  do morton_no=1,nleafs
209  ! This is where we ideally would like to make the cuts:
210  mstop = (ipe+1)*int(mtot/npe)+min(ipe+1,mod(mtot,npe))
211  ! Build up mass:
212  mcurr = mcurr + (wa*sfc(3,morton_no)+wp*(1-sfc(3,morton_no)))
213 
214  if (sfc(3,morton_no)==1) then
215  nactive(ipe) = nactive(ipe) +1
216  else
217  npassive(ipe) = npassive(ipe) +1
218  end if
219 
220  if (mcurr >= mstop) then
221  morton_stop(ipe) = morton_no
222  ipe = ipe +1
223  if (ipe>=npe) exit
224  morton_start(ipe) = morton_no + 1
225  end if
226  end do
227 
228  xmemory=dble(maxval(npassive+nactive))/&
229  dble(minval(npassive+nactive))
230  xload=dble(maxval(nactive))/&
231  dble(minval(nactive))
232 
233  !ptasum = ptasum +dble(nleafs-nleafs_active)/dble(nleafs_active)
234 
235  !if (mype == 0) print*, 'nleafs_passive:',nleafs-nleafs_active, 'nleafs_active:',nleafs_active,'ratio:',dble(nleafs-nleafs_active)/dble(nleafs_active),'mean ratio:',ptasum/it
236 
237  if (morton_stop(mype)>=morton_start(mype)) then
239  {#IFDEF EVOLVINGBOUNDARY
240  allocate(sfc_phybound(nleafs))
241  sfc_phybound=0
242  }
243  end if
244 
245  end subroutine get_morton_range_active
246 
247 end module mod_space_filling_curve
integer, dimension(:^d &), allocatable, save iglevel1_sfc
iglevel1_sfc(ig^D) gives the Morton number for grid ig^D
Definition: mod_forest.t:50
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine get_morton_range
Set the Morton range for each processor.
integer nglev1
Definition: mod_forest.t:78
Module with basic grid data structures.
Definition: mod_forest.t:2
integer npe
The number of MPI tasks.
integer, dimension(:,:), allocatable, save sfc
Array to go from a Morton number to an igrid and processor index. Sfc(1:3, MN) contains [igrid...
Definition: mod_forest.t:43
integer(kind=8) function mortonencode(igD, ndim)
Pointer to a tree_node.
Definition: mod_forest.t:6
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
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
double precision xload
Stores the memory and load imbalance, used in printlog.
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(:^d &), allocatable, save tree_root
Pointers to the coarse grid.
Definition: mod_forest.t:29
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition: mod_forest.t:62
integer nleafs_active
Definition: mod_forest.t:78
integer, dimension(:), allocatable, parameter d
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition: mod_forest.t:65
integer mype
The rank of the current MPI task.
subroutine level1_morton_order
build Morton space filling curve for level 1 grid
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, parameter ndim
Number of spatial dimensions for grid variables.
recursive subroutine get_morton_number(tree)
subroutine amr_morton_order
Construct Morton-order as a global recursive lexicographic ordering.
integer, dimension(:), allocatable, save sfc_phybound
Space filling curve used for physical boundary blocks.
Definition: mod_forest.t:59
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, save nleafs
Number of leaf block.
Definition: mod_forest.t:76