MPI-AMRVAC  3.0
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  if(.not. allocated(iglevel1_sfc)) allocate(iglevel1_sfc(ng^d(1)))
46  if(.not. allocated(sfc_iglevel1)) 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
86  use mod_forest
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)
105 
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
132  use mod_forest
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 
167  use mod_forest
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
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
recursive subroutine get_morton_number(tree)
Module with basic grid data structures.
Definition: mod_forest.t:2
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, dimension(:), allocatable, save sfc_phybound
Space filling curve used for physical boundary blocks.
Definition: mod_forest.t:59
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition: mod_forest.t:62
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
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, 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
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision xload
Stores the memory and load imbalance, used in printlog.
integer, parameter ndim
Number of spatial dimensions for grid variables.
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 npe
The number of MPI tasks.
subroutine get_morton_range
Set the Morton range for each processor.
integer(kind=8) function mortonencode(igD, ndim)
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