DISPATCH
amr_commons.f90
1 #define MAXLEVEL 30
2 
3 module amr_commons
4  use amr_parameters
5 
6  integer::mpi_comm_use ! communicator to nodes in use
7  integer::bitwise_test=0 ! =1 -> write, =2 -> compare
8  logical::idle_nodes=.false. ! test for bad nodes and idle them
9  logical::do_output=.true. ! Output just performed
10  logical::output_done=.false. ! Output just performed
11  logical::init=.false. ! Set up or run
12  logical::balance=.false. ! Load balance or run
13  logical::shrink=.false. ! Shrink mesh or run
14  logical::orphan=.false. ! Thread with myid > ncpu2
15  logical::extra_dump=.false. ! trigger extra dump
16  logical::extra_load_balance=.false. ! trigger extra load_balance
17  logical::tight_nbor=.false. ! require 3 neighbors in 3rd flag step
18  integer::nstep=0 ! Time step
19  integer::nstep_coarse=0 ! Coarse step
20  integer::nstep_coarse_old=0 ! Old coarse step
21  integer::nflag ! Refinements
22  integer::ncoarse ! nx.ny.nz
23  integer::ngrid_current ! Actual number of octs
24  integer(8)::n_cell ! count cells per main step
25  integer(kind=8)::allocate_size=0
26 
27  real(dp)::emag_tot=0.0d0 ! Total magnetic energy
28  real(dp)::ekin_tot=0.0d0 ! Total kinetic energy
29  real(dp)::eint_tot=0.0d0 ! Total internal energy
30  real(dp)::epot_tot=0.0d0 ! Total potential energy
31  real(dp)::epot_tot_old=0.0d0 ! Old potential energy
32  real(dp)::epot_tot_int=0.0d0 ! Time integrated potential
33  real(dp)::const=0.0d0 ! Energy conservation
34  real(dp)::aexp_old=1.0d0 ! Old expansion factor
35  real(dp)::rho_tot=0.0d0 ! Mean density in the box
36  real(dp)::t=0.0d0 ! Time variable
37  real(dp)::tset=0.0d0 ! Time variable
38  real(kind=8)::t8=0.0d0 ! Time variable used internally. 64 bit.
39 
40  logical::local_debug=.false.
41  real(kind=8)::dr_debug(3),r_debug(3)
42  logical::dbg_grid(nvector)
43  integer::dbg_ind(nvector)
44 
45  ! MPI variables
46  integer::ncpu,ncpu_dump,ndomain,myid,overload=1,ncore
47  integer, allocatable, dimension(:) :: isub_domain, cpu_domain ! Given domain, return local sub domain and cpu
48  integer, allocatable, dimension(:,:) :: map_domain ! Given local sub domain and cpu return domain
49 
50  ! Friedman model variables
51  integer::n_frw
52  real(dp),allocatable,dimension(:)::aexp_frw,hexp_frw,tau_frw,t_frw
53 
54  ! Initial conditions parameters from grafic
55  integer ::nlevelmax_part
56  real(dp) ::aexp_ini=10.
57  real(dp),dimension(1:MAXLEVEL)::dfact=1.0d0,astart
58  real(dp),dimension(1:MAXLEVEL)::vfact,dtlev
59  real(dp),dimension(1:MAXLEVEL)::xoff1,xoff2,xoff3,dxini
60  integer ,dimension(1:MAXLEVEL)::n1,n2,n3
61 
62  ! Level related arrays
63  real(dp) ::dtlevelmin ! Largest possible timestep at the coarsest level
64  real(dp),dimension(1:MAXLEVEL)::dtold,dtnew ! Time step at each level
65  real(dp),dimension(1:MAXLEVEL)::rho_max ! Maximum density at each level
66  integer ,dimension(1:MAXLEVEL)::nsubcycle=2 ! Subcycling at each level
67 
68  ! Pointers for each level linked list
69  integer,allocatable,dimension(:,:)::headl
70  integer,allocatable,dimension(:,:)::taill
71  integer,allocatable,dimension(:,:)::numbl
72  integer,allocatable,dimension(:,:)::numbtot
73 
74  ! Pointers for each level boundary linked list
75  integer,allocatable,dimension(:,:)::headb
76  integer,allocatable,dimension(:,:)::tailb
77  integer,allocatable,dimension(:,:)::numbb
78 
79  ! Pointers for free memory grid linked list
80  integer::headf,tailf,numbf,used_mem,used_mem_tot
81 
82  ! Tree arrays
83  real(dp),allocatable,dimension(:,:)::xg ! grids position
84  integer ,allocatable,dimension(:,:)::nbor ! neighboring father cells
85  integer ,allocatable,dimension(:) ::father ! father cell
86  integer ,allocatable,dimension(:) ::next ! next grid in list
87  integer ,allocatable,dimension(:) ::prev ! previous grid in list
88  integer ,allocatable,dimension(:) ::son ! sons grids
89  integer ,target, allocatable,dimension(:) ::flag1 ! flag for refine
90  integer ,target, allocatable,dimension(:) ::flag2 ! flag for expansion
91 
92  ! Global indexing
93  integer ,allocatable,dimension(:) ::cpu_map ! domain decomposition
94  integer ,allocatable,dimension(:) ::cpu_map2 ! new domain decomposition
95 
96  ! Hilbert key
97  real(qdp),allocatable,dimension(:)::hilbert_key
98  real(qdp),allocatable,dimension(:)::bound_key,bound_key2
99  real(qdp) ::order_all_min,order_all_max
100 
101  ! Recursive bisection
102  real(dp),allocatable,dimension(:) ::bisec_wall ! bisection wall positions
103  integer ,allocatable,dimension(:,:) ::bisec_next ! next 2 child cells in bisection
104  integer::bisec_root ! root of bisection tree
105 
106  integer,allocatable,dimension(:) ::bisec_indx ! map from leaf cell id to cpu id
107  real(dp),allocatable,dimension(:,:) ::bisec_cpubox_min ! cpu domains boxes
108  real(dp),allocatable,dimension(:,:) ::bisec_cpubox_max
109  real(dp),allocatable,dimension(:,:) ::bisec_cpubox_min2 ! cpu domains boxes for new decomp
110  real(dp),allocatable,dimension(:,:) ::bisec_cpubox_max2
111 
112  integer,allocatable,dimension(:) ::bisec_cpu_load ! CPU loads (for stats)
113  integer,allocatable,dimension(:,:) ::bisec_hist ! histograms for load computation
114  integer,allocatable,dimension(:) ::bisec_hist_bounds ! histogram splitting boundaries
115  integer,allocatable,dimension(:) ::new_hist_bounds
116  integer,allocatable,dimension(:) ::bisec_ind_cell ! histo swap id -> cell id map (big)
117  integer,allocatable,dimension(:) ::cell_level ! store the level of the cells (big)
118 
119  real(dp)::bisec_res ! resolution parameters
120  integer(kind=8)::bisec_nres
121 
122  ! Communication structure
124  integer ::ngrid=0
125  integer ::npart=0
126  integer ,dimension(:) ,pointer::igrid => null()
127  integer ,dimension(:,:),pointer::f
128  real(kind=8),dimension(:,:),pointer::u
129  integer ,dimension(:,:),pointer::fp
130  real(kind=8),dimension(:,:),pointer::up
131 #ifdef ATON
132  real(kind=8),dimension(:,:),pointer::u_radiation
133 #endif
134  end type communicator
135 
136  ! Active grid, emission and reception communicators
137  type(communicator),allocatable,dimension(:) ::active
138  type(communicator),allocatable,dimension(:,:)::boundary
139  type(communicator),allocatable,dimension(:,:)::emission
140  type(communicator),allocatable,dimension(:,:)::reception
141 
142  ! Types for physical boundary conditions
143  CHARACTER(LEN=20)::type_hydro ='hydro'
144  CHARACTER(LEN=20)::type_accel ='accel'
145  CHARACTER(LEN=20)::type_flag ='flag'
146 
147  ! Units specified by the user in the UNITS_PARAMS namelist for non-cosmo runs.
148  ! These values shouldn't be used directly. Instead call units() in amr/units.f90.
149  real(dp)::units_density=1.0 ! [g/cm^3]
150  real(dp)::units_time=1.0 ! [seconds]
151  real(dp)::units_length=1.0 ! [cm]
152  real(dp)::units_velocity=0.0 ! [cm/s] (if set, replaces units_time)
153 
154  ! Counters for refinement
155  integer n_refine_d, n_refine_p, n_refine_u, n_refine_c, n_refine_s, n_refine_o, &
156  n_refine_m, n_refine_l, n_refine_q, n_refine_b
157  integer n_gfloor_u, n_gfloor_b
158 
159  ! Type for supporting restart from arbitrary type of files
161  ! Switches detailing the input data file format
162 #ifndef NPRE
163  logical :: double_precision_run = .false.
164 #else
165 #if NPRE==4
166  logical :: double_precision_run = .false.
167 #else
168  logical :: double_precision_run = .true.
169 #endif
170 #endif
171 #ifdef QUADHILBERT
172  logical :: quad_precision_run = .true.
173 #else
174  logical :: quad_precision_run = .false.
175 #endif
176  logical :: big_endian, double_precision, quad_precision
177  end type
178  type(io_restart_t) :: io_restart
179 !
180  interface load_float
181  module procedure load_float_scalar, load_float_1d
182  end interface
183 contains
184 !==============================================================================!
185 subroutine print_id (id)
186  implicit none
187  character(len=*) id
188  !real(kind=8) wallclock
189 !-----------------------------------------------------------------------
190  !$omp critical
191  if (myid==1 .and. id .ne. ' ') then
192  !print '(1x,a,f12.2)', trim(id), wallclock()
193  id = ' '
194  end if
195  !$omp end critical
196 end subroutine
197 !==============================================================================!
198 subroutine load_float_scalar(ilun,x,n)
199  implicit none
200  integer, intent(in) :: ilun, n
201  real(dp) :: x
202  real(kind=8) :: x_dp
203  real(kind=4) :: x_sp
204  if (io_restart%double_precision) then
205  read(ilun) x_dp
206  x = x_dp
207  else
208  read(ilun) x_sp
209  x = x_sp
210  end if
211 end subroutine load_float_scalar
212 !==============================================================================!
213 subroutine load_float_1d(ilun,x,n)
214  implicit none
215  integer, intent(in) :: ilun, n
216  real(dp), dimension(:) :: x
217  real(kind=8), dimension(:), allocatable :: x_dp
218  real(kind=4), dimension(:), allocatable :: x_sp
219  if (io_restart%double_precision .eqv. io_restart%double_precision_run) then ! no need for scratch variable
220  read(ilun) x(1:n)
221  else
222  if (io_restart%double_precision) then
223  allocate(x_dp(n))
224  read(ilun) x_dp
225  x(1:n) = x_dp
226  deallocate(x_dp)
227  else
228  allocate(x_sp(n))
229  read(ilun) x_sp
230  x(1:n) = x_sp
231  deallocate(x_sp)
232  end if
233  end if
234 end subroutine load_float_1d
235 !==============================================================================!
236 end module amr_commons
237 !==============================================================================!
238 module reduce_comm_m
239  USE amr_commons
240  implicit none
241  ! Generic type for encapsulating information for a level-communicator
243  logical :: inuse=.false. ! is it allocated ?
244  integer :: comm ! communicator handle
245  integer :: nodes ! #nodes in comm
246  integer :: rank ! rank in communicator
247  integer, dimension(:), allocatable :: g2l ! mapping from global ranks to communicator
248  integer, dimension(:), allocatable :: l2g ! mapping from communicator ranks to global
249  integer, dimension(:), allocatable :: nemission ! nr of grids to emit in global coordinates
250  end type
251 
252  ! communicator for exchanging ghostzones between nodes on a given level
253  type(level_communicator_t), dimension(:), allocatable :: full_level
254 
255  ! communicator for exchanging ghostzones between nodes that have active grids on a given level
256  type(level_communicator_t), dimension(:), allocatable :: virtual_level
257 end module reduce_comm_m
258 !==============================================================================!
259 module openmp_support
260  use amr_commons
261  implicit none
262  integer :: omp_nthreads, omp_mythread,mpi_omp_threads
263  logical :: omp_master
264  !$omp threadprivate(omp_mythread,omp_master)
265  character(len=80), save:: id='amr_commons.f90 $Id: 0082ecb9f612def63f2fb80a39ec9d2d3b71bcf6 $'
266 contains
267 ! Function that resets the stacksize of OpenMP threads to avoid segfaults
268 #if defined (_OPENMP)
269 subroutine set_openmp_stacksize
270 #ifdef __INTEL_COMPILER
271  use hydro_parameters, only : nvar
272  use omp_lib, only : kmp_set_stacksize_s, kmp_get_stacksize_s, kmp_size_t_kind
273 #endif
274  implicit none
275 #ifdef __INTEL_COMPILER
276  integer(kind=kmp_size_t_kind) :: omp_stacksize, new_omp_stacksize
277 
278  omp_stacksize = kmp_get_stacksize_s()
279  new_omp_stacksize = (4*1024**2 * nvar) / 8 ! Adjust according to number of passive scalars, with baseline being 4M for MHD
280  if (dp==kind(1.0e0_8)) new_omp_stacksize = new_omp_stacksize * 2 ! Double for double precision runs
281 
282  if (new_omp_stacksize > omp_stacksize) then
283  call kmp_set_stacksize_s(new_omp_stacksize) ! Set stacksize to needed value
284  if (myid==1) write(*,*) 'WARNING! OpenMP stacksize has been reset to ', new_omp_stacksize/1024**2, 'MB to avoid segfault'
285  endif
286 #endif
287 end subroutine set_openmp_stacksize
288 #endif
289 !
290 subroutine init_openmp
291 #if defined (_OPENMP)
292  use omp_lib, only : omp_get_thread_num, omp_get_num_threads
293 #endif
294  implicit none
295 #if defined (_OPENMP)
296  ! Make sure stack size is large enough -- if compiler supports setting it on-the-fly
297  call set_openmp_stacksize
298  !
299 !$omp parallel
300 !$omp master
301  omp_nthreads = omp_get_num_threads()
302 !$omp end master
303  omp_mythread = omp_get_thread_num()
304 #else
305  omp_mythread = 0
306  omp_nthreads = 1
307 #endif
308  omp_master = omp_mythread .eq. 0
309 !$omp end parallel
310  mpi_omp_threads = (ncpu-1)/omp_nthreads+1 ! Ratio of mpi to openmp threads. For hybrid communication.
311  ncore = ncpu*omp_nthreads
312 !
313  call print_id(id)
314  if (myid==1) print'(1x,a,4i8)','init_openmp: ncpu, ncore, omp_nthreads, mpi_omp_threads =', &
315  ncpu, ncore, omp_nthreads, mpi_omp_threads
316 end subroutine init_openmp
317 end module openmp_support
318 !==============================================================================!