DISPATCH
patch_mod.f90
1 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2 !> Template module for patches, which adds pointers to memory and mesh, and
3 !> number of dimensions and variables. The number of dimensions is hardwired
4 !> to 3, since we can just choose the actual dimensions to get 2-D and 1-D.
5 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
6 MODULE patch_mod
7  USE iso_fortran_env, only: int8
8  USE io_mod
9  USE mpi_mod
10  USE omp_mod
11  USE omp_lock_mod
12  USE omp_timer_mod
13  USE timer_mod
14  USE mpi_mesg_mod
15  USE bits_mod
16  USE trace_mod
17  USE mpi_mod
18  USE kinds_mod
19  USE task_mod
20  USE link_mod
21  USE mesh_mod
22  USE boundaries_mod
23  USE index_mod
24  USE shared_mod
25  USE connect_mod
26  implicit none
27  private
28  integer, parameter:: ndim = 3
29  type, extends(task_t), public:: patch_t
30  real(kind=KindScalarVar), dimension(:,:,:,:,:,:), pointer:: mem => null()
31  real(kind=KindScalarVar), dimension(:,:,:,:,:), pointer:: vgas => null()
32  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: density => null()
33  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: pressure => null()
34  !dir$ attributes align: 64 :: mem, vgas, density, pressure
35  !---------------------------------------------------------------------------
36  ! These allocatables are intended to deposite external forces and heating in,
37  ! typically used by methods called from the extras_t data type procedures.
38  !---------------------------------------------------------------------------
39  real(kind=KindScalarVar), dimension(:,:,:,:), allocatable:: force_per_unit_mass
40  real(kind=KindScalarVar), dimension(:,:,:,:), allocatable:: force_per_unit_volume
41  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: heating_per_unit_mass
42  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: heating_per_unit_volume
43  integer(kind=int8) , dimension(:,:,:), allocatable:: filled
44  !dir$ attributes align: 64 :: force_per_unit_mass, force_per_unit_volume
45  !dir$ attributes align: 64 :: heating_per_unit_mass, heating_per_unit_volume
46  !dir$ attributes align: 64 :: filled
47  !--- anonymous pointer, which may be used to connect extra info to the solvers
48  type(connect_t):: connect
49  class(mesh_t), dimension(:), pointer:: mesh => null()
50  real:: courant ! courant number
51  real(8):: pdf_next=0.0 ! next PDF output time
52  real(8):: gamma=1.4 ! gas gamma
53  real(8):: dt_fixed ! for fixed time step Courant condition
54  real:: u_fixed, u_max ! for fixed time step Courant condition
55  integer:: ndim = ndim ! number of dimensions
56  integer:: nv=0 ! number of variables
57  integer:: nw=2 ! number of work indices
58  integer:: msplit=0 ! how to split a patch when refining
59  integer:: max_files ! prevent disk overrun when NaN happens
60  integer:: check_refine_next=0
61  integer:: time_derivs ! for restart changes
62  integer:: not_filled=0 ! guard zones not filled
63  integer:: refine_ratio=2 ! may differ btw patches
64  logical:: guard_zones ! for restart changes
65  logical:: mem_allocated=.false.
66  logical:: do_pebbles=.false.
67  logical:: staggered=.true.
68  logical:: no_mans_land=.false.
69  logical:: use_data_hub=.false.
70  logical:: manual_refine=.false.
71  logical:: get_umax_location=.false.
72  logical, dimension(:), pointer:: unsigned => null()
73  logical, dimension(:), pointer:: pervolume => null()
74  class(link_t), pointer:: link ! this could be useful
75  type(boundaries_t):: boundaries
76  character(len=32):: eos='ideal'
77  character(len=32):: opacity='none'
78  integer:: file_offset(4)
79  type(index_t):: idx
80  real:: safety_factor=1.0
81  !--- used in extras ---
82  logical, allocatable:: mask(:,:,:)
83  !--- used in refine ---
84  integer(int8), allocatable:: irefine(:,:,:)
85  type(lock_t), pointer:: mem_lock(:) => null()
86  !---------------------------------------------------------------------------
87  ! HACKs below: FIXME
88  !---------------------------------------------------------------------------
89  real(8):: theta=0d0 ! latitude
90  real(8):: phi=0d0 ! longitude at t=0.0
91  real:: csound ! isothermal sound speed, for tests
92  real:: grav ! gravitational consant (code units)
93  integer:: n(3)=0, ng(3)=0 ! defaults, to detect unchanged values
94  integer:: li(3), ui(3), gn(3) ! convenience -- avoid patch%mesg%li etc
95  integer:: ncell(3) ! number of cells per patch
96  integer:: ipos(3) ! integer position, 3D
97  real(8):: offset(3) ! this is used to compute indices, locations
98  real(8):: llc_cart(3)=0d0 ! the lower left corner of the patch in *Cartesian* coords.
99  real(8):: llc_nat(3)=0d0 ! the lower left corner of the patch in its *own* coords.
100  real(8):: centre_cart(3)=0d0 ! the geometric centre of the patch in *Cartesian* coords.
101  real(8):: centre_nat(3)=0d0 ! the geometric centre of the patch in its *own* coords.
102  integer :: mesh_type = 1 ! default: Cartesian mesh; FIXME (there must be a better way)
103  real(8):: grid_vel(3) = 0d0 ! velocity of the patch in its native coords. and
104  ! units (e.g. rad/s for the phi coord.); for
105  ! Cartesian coords., equal to `velocity`.
106  contains
107  procedure, nopass:: cast2patch
108  procedure:: init
109  procedure:: init_default
110  procedure:: init_bdries
111  procedure:: setup
112  procedure:: dealloc
113  procedure:: clone_mem_accounting
114  procedure:: rotate
115  procedure:: allocate_mesg ! allocate a mesg buffer
116  procedure:: pack ! pack patch into buffer
117  procedure, nopass:: unpack_id ! unpack patch id
118  procedure:: unpack ! unpack buffer into patch
119  procedure:: info ! print info to stdout
120  procedure, private:: contains1 ! patch contains a point
121  procedure, private:: contains2 ! patch contains a patch
122  generic:: contains => contains1, contains2
123  procedure:: overlaps ! patch overlaps test
124  procedure:: index ! compute the index for a single point
125  procedure:: nearest ! the nearest index that maps to patch point
126  procedure:: lowest ! the lowest index that maps to patch point
127  procedure:: highest ! the highest index that maps to patch point
128  procedure:: myposition
129  procedure:: myposition_staggered
130  procedure:: periodic_grid
131  procedure:: extrapolate_guards
132  procedure:: is_periodic
133  procedure, private :: make_periodic4
134  procedure, private :: make_periodic8
135  generic, public :: make_periodic => make_periodic4 , make_periodic8
136  procedure, private :: fsum4
137  procedure, private :: fsum8
138  generic, public :: fsum => fsum4, fsum8 ! interior sum
139  procedure, private :: faver4
140  procedure, private :: faver8
141  generic, public :: faver => faver4, faver8 ! interior average
142  procedure, private :: fminval4
143  procedure, private :: fminval8
144  procedure, private :: fminvali
145  generic, public :: fminval => fminval4, fminval8, fminvali
146  procedure, private :: fmaxval4
147  procedure, private :: fmaxval8
148  procedure, private :: fmaxvali
149  generic, public :: fmaxval => fmaxval4, fmaxval8, fmaxvali
150  generic, public :: fminvalvar => fminval4, fminval8 ! these two overloaded functions are
151  generic, public :: fmaxvalvar => fmaxval4, fmaxval8 ! used by `fminvali`/`fmaxvali` only.
152  procedure:: courant_condition
153  procedure:: check_density
154  procedure:: check_nan
155  procedure:: custom_refine
156  procedure:: patch_to_header
157  procedure:: header_to_patch
158  procedure, private:: stats_patch
159  procedure, private:: stats_scalar
160  procedure, private:: stats_vector
161  generic, public:: stats => stats_patch, stats_scalar, stats_vector
162  procedure:: mapvectors
163  procedure:: beyond_patch_edge
164  procedure:: myfloatindex_scalar
165  procedure:: log_interpolate
166  procedure:: time_interval
167  procedure:: distance => distance_curvilinear
168  procedure:: index_stagger
169  procedure:: interpolate
170  procedure:: interpolate_specific
171  procedure:: index_space ! full mapping to index space
172  procedure:: index_space_of
173  procedure:: index_space2
174  procedure:: index_only
175  procedure:: index_only2
176  procedure:: count_edges
177  procedure:: update_position
178  procedure:: get_overlap
179  procedure:: init_level
180  end type
181  !-----------------------------------------------------------------------------
182  ! Sequential header, for MPI and I/O
183  !-----------------------------------------------------------------------------
184  integer, public, save:: n_header = (14*4 + 4*8 + 5*8*3 + 1*8 + 2*8*8 + 1*4*8 + 2*4*3)/4
185  type, public:: header_t
186  sequence
187  integer:: version=1
188  integer:: n_header
189  integer:: n_data
190  integer:: id
191  integer:: status
192  integer:: rank
193  integer:: istep
194  integer:: iout
195  integer:: level
196  integer:: nt
197  integer:: nv ! number of variables
198  integer:: it
199  integer:: nq=0
200  integer:: npad=0
201  real(8):: time
202  real(8):: out_next
203  real(8):: print_next
204  real(8):: dtime
205  real(8):: position(3)
206  real(8):: velocity(3)
207  real(8):: size(3)
208  real(8):: box(3) ! the periodic wrapping size
209  real(8):: ds(3)
210  real(8):: send_time
211  real(8), dimension(8):: t
212  real(8), dimension(8):: dt
213  integer, dimension(8):: iit
214  integer:: n(3)=16 ! convenience
215  integer:: ng(3)=2 ! convenience
216  end type
217  !-----------------------------------------------------------------------------
218  ! Template task, with procedure to place it it MPI-space
219  !-----------------------------------------------------------------------------
220  type box_t
221  real:: size(3)=1d0
222  real:: num(3)
223  integer:: dim(3)
224  end type
225  integer, save:: print_every=0, iprint=0 ! print cadence
226  logical, save:: extrapolate_gz=.false. ! extrapolate guard zones?
227  logical, save:: do_check_nan=.false. ! check for NaN?
228  logical, save:: zero_order_extrap=.true. ! extrap to no-mans-land
229  public task2patch
230 CONTAINS
231 
232 !===============================================================================
233 !> Copy he most important fields to a sequenced header, for communication and I/O
234 !===============================================================================
235 SUBROUTINE patch_to_header (self, head)
236  class(patch_t):: self
237  type(header_t):: head
238  integer:: nt
239  !-----------------------------------------------------------------------------
240  head%n_header = n_header
241  head%n_data = product(self%gn)*self%nv
242  head%id = self%id
243  head%status = self%status
244  head%rank = self%rank
245  head%istep = self%istep
246  head%iout = self%iout
247  head%level = self%level
248  head%time = self%time
249  head%out_next = self%out_next
250  head%print_next = self%print_next
251  head%dtime = self%dtime
252  head%position = self%position
253  head%velocity = self%velocity
254  head%size = self%size
255  head%it = self%it
256  head%nt = self%nt
257  nt = self%nt
258  head%t(1:nt) = self%t
259  head%dt(1:nt) = self%dt
260  head%iit(1:nt) = self%iit
261  head%nq = self%nq
262  head%send_time = wallclock()
263  if (io%verbose > 1) &
264  write (io_unit%log,'(f12.6,3x,a,i7,i4,l3,2i4,f8.4)') &
265  wallclock(), 'patch_to_header: id, nq =', self%id, &
266  self%nq, self%is_set(bits%swap_request), &
267  mpi_mesg%sent_list%n, mpi_mesg%recv_list%n
268 END SUBROUTINE patch_to_header
269 
270 !===============================================================================
271 !> Copy he most important fields to a sequenced header, for communication and I/O
272 !===============================================================================
273 SUBROUTINE header_to_patch (self, head)
274  class(patch_t):: self
275  type(header_t):: head
276  integer:: nt, new
277  !-----------------------------------------------------------------------------
278  ! Swap the roles of boundary and virtual on incoming patches. Also clear the
279  ! ready bit, which is an indicator of a patch being in the ready queue, which
280  ! an incoming patch is not. By doing this here, we avoid having to do it in
281  ! a critical region, otherwise needed to protect an operation on the patch
282  !-----------------------------------------------------------------------------
283  if (iand(head%status, bits%boundary) /= 0) then
284  head%status = ior(head%status, bits%virtual)
285  head%status = iand(head%status, not(bits%boundary))
286  if (io%verbose>2) &
287  write (io_unit%mpi, *) head%id, 'swapping to virtual'
288  else if (iand(head%status, bits%virtual) /= 0) then
289  head%status = ior(head%status, bits%boundary)
290  head%status = iand(head%status, not(bits%virtual))
291  if (io%verbose>2) &
292  write (io_unit%mpi, *) head%id, 'swapping to boundary'
293  end if
294  head%status = iand(head%status, not(bits%ready+bits%busy))
295  if (head%id /= self%id) then
296  call io%abort('wrong message ID unpacked')
297  end if
298  !-----------------------------------------------------------------------------
299  ! Except for the first and last appearance, istep should increase be one
300  !-----------------------------------------------------------------------------
301  if ((self%istep > 0) .and. (iand(head%status, bits%remove) == 0)) then
302  if (head%istep > self%istep+1) then
303  write(io_unit%log,1) &
304  'WARNING: early arrival in time reversal id, rank =', &
305  self%id, mpi%rank, head%istep, self%istep, head%time, self%time
306  1 format(a,4i6,2f12.6)
307  flush (io_unit%log)
308  else if (head%istep < self%istep+1) then
309  write(io_unit%log,1) &
310  'WARNING: late arrival in time reversal id, rank =', &
311  self%id, mpi%rank, head%istep, self%istep, head%time, self%time
312  flush (io_unit%log)
313  end if
314  end if
315  self%id = head%id
316  self%status = head%status
317  self%rank = head%rank
318  self%istep = head%istep
319  self%iout = head%iout
320  self%level = head%level
321  self%time = head%time
322  self%out_next = head%out_next
323  self%print_next = head%print_next
324  self%dtime = head%dtime
325  self%position = head%position
326  self%velocity = head%velocity
327  self%size = head%size
328  !$omp atomic write
329  self%it = head%it
330  new = mod(head%it,head%nt) + 1
331  !$omp atomic write
332  self%new = new
333  self%nt = head%nt
334  nt = head%nt
335  if (nt > 0) then
336  self%t = head%t(1:nt)
337  self%dt = head%dt(1:nt)
338  self%iit = head%iit(1:nt)
339  !self%new = head%iit(nt)
340  else
341  write (io_unit%log,*) omp%thread,'header_to_patch WARNING: nt non-positive', nt
342  flush (io_unit%log)
343  end if
344  self%nq = head%nq
345  self%latency = wallclock() - head%send_time
346  !$omp atomic
347  timer%latency%max = max(timer%latency%max,self%latency)
348  !$omp atomic
349  timer%latency%aver = timer%latency%aver+self%latency
350  !$omp atomic
351  timer%latency%n = timer%latency%n+1
352  !.............................................................................
353 END SUBROUTINE header_to_patch
354 
355 !===============================================================================
356 !> Cast a generic task_t to patch_t
357 !===============================================================================
358 FUNCTION cast2patch (task) RESULT(patch)
359  class(task_t), pointer:: task
360  class(patch_t), pointer:: patch
361  !.............................................................................
362  select type (task)
363  class is (patch_t)
364  patch => task
365  class default
366  nullify(patch)
367  call io%abort ('patch_t%cast: failed to cast a task to patch_t')
368  end select
369 END FUNCTION cast2patch
370 
371 !===============================================================================
372 !> Public function (deprecated)
373 !===============================================================================
374 FUNCTION task2patch (task) RESULT (patch)
375  class(task_t), pointer:: task
376  class(patch_t), pointer:: patch
377  !.............................................................................
378  patch => cast2patch(task)
379 END FUNCTION task2patch
380 
381 !===============================================================================
382 !> Initialize the patch memory and patch dimensions. For this we need to be told the
383 !> patch dimensions (which include the number of variable and the number of time
384 !> slices), the patch position, and the patch size
385 !===============================================================================
386 SUBROUTINE init (self)
387  class(patch_t):: self
388  character(len=120):: id = &
389  '$Id: f01f80ce28846df45c29d0d9c93488206e197f17 $ tasks/patch_mod.f90'
390  !-----------------------------------------------------------------------------
391  call trace%begin ('patch_t%init')
392  call trace%print_id (id)
393  call self%setup
394  call trace%end()
395 END SUBROUTINE init
396 
397 !===============================================================================
398 !> Initialize default values for the patch memory and patch dimensions. The
399 !> order of setting values is that the calling code, typically experiment_mod,
400 !> may set preferred defaults that differ from the ones set here. The input
401 !> namelist then allows these to be changed. This is done (only) once, on first
402 !> call, after which the saved properties are reused on subsequent calls.
403 !>
404 !> However, in cases where not all patches should have identical properties it
405 !> is necessary be able to override the saved default. In the case of n and ng
406 !> -- the most likely properties to vary -- this can be done by changing the patch
407 !> values before the call, since only properties with default (zero) values
408 !> are set to the saved values. It is possible to change other parameters, by
409 !> instead of calling patch_t%init, first call patch_t%init_default, and then
410 !> call patch_t%setup (nt=..., nv=..., ...).
411 !>
412 !> Values that are expected to vary from patch to patch, such as size and
413 !> position, must be set before calling patch_t%init, or else must be set by
414 !> instead calling patch_t%init_defult + patch_t%setup.
415 !>
416 !> Any experiment that should start up with non-default values of parameters
417 !> given default values here (such as 'periodic'), should use the method:
418 !> call init_default + set values + call setup, since otherwise it would have
419 !> to rely on the values being set correctly in the input namelist.
420 !===============================================================================
421 SUBROUTINE init_default (self)
422  class(patch_t):: self
423  !.............................................................................
424  integer, save:: n(3)=16, ng(3)=2, nv=5, nt=5, nw, max_files=1000
425  logical, save:: no_mans_land=.false., use_data_hub=.false., limit_dtime=.false.
426  real, save:: u_fixed=-1.0, grace=0.0
427  real(8), save:: end_time=1.0, out_time=0.1, print_time=0.1, dt_fixed=-1.0
428  real(8):: position(3), size(3)
429  integer:: iostat, i, n_alloc=0, n_ext=0
430  namelist /patch_params/ n, nv, ng, nt, u_fixed, dt_fixed, &
431  grace, extrapolate_gz, no_mans_land, use_data_hub, do_check_nan, &
432  zero_order_extrap, limit_dtime
433  namelist /out_params/ out_time, end_time, print_time, print_every, &
434  max_files
435  logical, save:: first_time = .true.
436  !-----------------------------------------------------------------------------
437  call trace%begin ('patch_t%init_default')
438  call self%task_t%init ! sets unique id
439  !-----------------------------------------------------------------------------
440  ! Read and set patch parameters; the preferred defalts may have been set by
441  ! the calling code -- typically in experiment_mod.f90
442  !-----------------------------------------------------------------------------
443  !$omp critical (input_cr)
444  if (first_time) then
445  ! -- if values have been imposed from above, FI show them as default
446  if (any(self%n /= 0)) n = self%n ! patch dimensions
447  if (any(self%ng /= 0)) ng = self%ng ! patch guard zones
448  if ( self%nv /= 0 ) nv = self%nv ! number of variables
449  if ( self%nt /= 0 ) nt = self%nt ! number of time slices
450  first_time = .false.
451  rewind(io%input)
452  read (io%input, patch_params)
453  if (dt_fixed > 0.0) grace=max(grace,1e-6)
454  if (io%master) write (*, patch_params)
455  rewind(io%input)
456  read (io%input, out_params)
457  if (io%master) write (*, out_params)
458  if (mpi%rank > 0) &
459  self%track = .true.
460  iprint = print_every ! print 1st task
461  end if
462  !-----------------------------------------------------------------------------
463  ! To allow patches to have n, ng, nv, and nt values that differ, values that
464  ! differ from the data type zero defaults take precedence over saved values
465  !-----------------------------------------------------------------------------
466  if (any(self%n == 0)) self%n = n ! patch dimensions
467  if (any(self%ng == 0)) self%ng = ng ! patch guard zones
468  if ( self%nv == 0 ) self%nv = nv ! number of variables
469  if ( self%nt == 0 ) self%nt = nt ! number of time slices
470  !-----------------------------------------------------------------------------
471  ! Output related parameters
472  !-----------------------------------------------------------------------------
473  io%end_time = end_time
474  io%out_time = out_time
475  io%print_time = print_time
476  if (io%nv==0) &
477  io%nv = self%nv
478  !$omp end critical (input_cr)
479  self%restart = io%restart ! restart snapshot
480  self%grace = grace ! grace fraction of dtime
481  self%dt_fixed = dt_fixed
482  self%u_fixed = u_fixed
483  self%max_files = max_files ! prevent run-away NaN
484  self%use_data_hub = use_data_hub ! guard zone handling
485  self%no_mans_land = no_mans_land ! mesh centering
486  self%print_next = print_time
487  self%limit_dtime = limit_dtime
488  call trace%end()
489 END SUBROUTINE init_default
490 
491 SUBROUTINE init_bdries (self)
492  class(patch_t):: self
493  real(8):: position(3), limit(3)
494  !-----------------------------------------------------------------------------
495  ! If at edge, set boundary bits
496  !-----------------------------------------------------------------------------
497  call trace%begin ('patch_t%init_bdries')
498  position = self%position - self%origin - 0.5d0 * self%box
499  limit = 0.5d0 * (self%box - self%size - self%ds)
500  if (.not.self%periodic(1)) then
501  if (position(1) <= -limit(1)) then
502  call self%boundaries%set (bits%xl)
503  self%mesh(1)%lower_boundary = .true.
504  end if
505  if (position(1) >= +limit(1)) then
506  call self%boundaries%set (bits%xu)
507  self%mesh(1)%upper_boundary = .true.
508  end if
509  end if
510  if (.not.self%periodic(2)) then
511  if (position(2) <= -limit(2)) then
512  call self%boundaries%set (bits%yl)
513  self%mesh(2)%lower_boundary = .true.
514  end if
515  if (position(2) >= +limit(2)) then
516  call self%boundaries%set (bits%yu)
517  self%mesh(2)%upper_boundary = .true.
518  end if
519  end if
520  if (.not.self%periodic(3)) then
521  if (position(3) <= -limit(3)) then
522  call self%boundaries%set (bits%zl)
523  self%mesh(3)%lower_boundary = .true.
524  end if
525  if (position(3) >= +limit(3)) then
526  call self%boundaries%set (bits%zu)
527  self%mesh(3)%upper_boundary = .true.
528  end if
529  end if
530  call trace%end()
531 END SUBROUTINE init_bdries
532 
533 !===============================================================================
534 !> Setup a patch, possibly with non-standard parameters. First setup with
535 !> defaults, then possibly override, then do the actual setup
536 !===============================================================================
537 SUBROUTINE setup (self, position, size, n, ng, nt, nv, nw)
538  class(patch_t) :: self
539  real(8), optional :: position(3), size(3)
540  integer, optional :: n(3), ng(3)
541  integer, optional :: nt, nv, nw
542  !.............................................................................
543  type(header_t):: header
544  character(len=4) kind
545  integer:: i
546  !-----------------------------------------------------------------------------
547  call trace%begin ('patch_t%setup')
548  call self%init_default ! read default params
549  if (present(position)) self%position = position
550  if (present(size )) self%size = size
551  if (present(n )) self%n = n
552  if (present(ng )) self%ng = ng
553  if (present(nt )) self%nt = nt
554  if (present(nv )) self%nv = nv
555  if (present(nw )) self%nw = nw
556 !print '(a,3(2x,3i4),2(2x,1p,3e11.2))', &
557 !'patch_t%setup: n,ng,nt,nv,nw,size,pos=', &
558 !self%n, self%ng, self%nv, self%nt, self%nw, self%size, self%position
559  !-----------------------------------------------------------------------------
560  ! Pass the desired patch dimension via the shared_mod to refine_mod
561  !-----------------------------------------------------------------------------
562  shared%patch_n = self%n
563  !-----------------------------------------------------------------------------
564  self%llc_cart = self%position - 0.5 * self%size
565  self%llc_nat = self%llc_cart
566  !-----------------------------------------------------------------------------
567  ! Detect root grid
568  !-----------------------------------------------------------------------------
569  if (all(self%size==self%box)) then
570  !$omp critical (log_cr)
571  write (io_unit%log,*) mpi%rank, 'task', self%id, ' is root grid'
572  call self%set(bits%root_grid)
573  !$omp end critical (log_cr)
574  end if
575  !-----------------------------------------------------------------------------
576  ! Special handling of collapsed dimensions
577  !-----------------------------------------------------------------------------
578  where (self%n==1) ! for collapsed dimensions
579  self%ng = 0 ! no guard zones
580  self%box = 1 ! make sure it is non-zero
581  end where
582  self%llc_cart = self%position - 0.5_8*self%size
583  if (self%mesh_type == mesh_types%Cartesian) then
584  self%llc_nat = self%llc_cart
585  self%centre_nat = self%position
586  end if
587  !-----------------------------------------------------------------------------
588  ! Initialize and allocate the meshes, afterwards imposing the external values
589  ! box size and position
590  !-----------------------------------------------------------------------------
591  call omp%set_stacksize (self%nv)
592  call meshfactory(self%mesh_type, self%mesh, n=self%n, ng=self%ng, s=self%size, &
593  p=self%position, llc_native=self%llc_nat, nml=self%no_mans_land)
594  ! intialise staggering indices; it is the responsibility of the *solver* to
595  ! assign non-zero values to this variable (if any).
596  do i=1,3
597  allocate(self%mesh(i)%h(self%nv))
598  self%mesh(i)%h = 0.0
599  if (self%kind == 'zeus_mhd_patch' .and. &
600  self%staggered .and. self%n(i) /= 1) self%mesh(i)%gn = self%mesh(i)%gn + 1
601  if (self%kind == 'zeus_tw_mhd_patch' .and. &
602  self%staggered .and. self%n(i) /= 1) self%mesh(i)%gn = self%mesh(i)%gn + 1
603  end do
604  self%ds = self%mesh%d
605  self%ncell = self%mesh%nc
606  self%mesh%id = self%id
607  self%mesh%b = self%box
608  self%mesh%llc_cart = self%llc_cart
609  self%centre_nat = [self%mesh(1)%centre_nat,self%mesh(2)%centre_nat,self%mesh(3)%centre_nat]
610  self%centre_cart = currenttocartesian(self%mesh_type, &
611  self%centre_nat(1), self%centre_nat(2), self%centre_nat(3))
612  !-----------------------------------------------------------------------------
613  ! Hack for backwards compatibility
614  !-----------------------------------------------------------------------------
615  self%gn = self%mesh%gn
616  self%li = self%mesh%li
617  self%ui = self%mesh%ui
618  self%gsize = (self%size*self%gn/self%n)
619  self%offset = (self%mesh%lb+self%mesh%ub)/2.0 ! mid point index offset
620  if (self%id == io%id_debug) then
621  print 1,'patch_t%init: position =', self%position
622  print 1,'patch_t%init: m%p =', self%mesh%p
623  print 1,'patch_t%init: offset =', self%offset
624  print 1,'patch_t%init: m%o =', self%mesh%o
625  print 1,'patch_t%init: centre_nat =', self%centre_nat
626  print 1,'patch_t%init: centre_cart =', self%centre_cart
627  print 1,'patch_t%init: llc_nat =', self%llc_nat
628  print 1,'patch_t%init: m%llc_nat =', self%mesh%llc_nat
629  print 1,'patch_t%init: llc_cart =', self%llc_cart
630  print 1,'patch_t%init: m%llc_cart =', self%mesh%llc_cart
631  1 format(a,1p,3e15.6)
632  end if
633  !-----------------------------------------------------------------------------
634  ! Measure levels in terms of the refinement factor -- make sure this doesn't
635  ! come in conflict with levels set by e.g. rubiks_mod
636  !-----------------------------------------------------------------------------
637  call self%init_level
638  if (io%verbose>2) then
639  !$omp critical (log_cr)
640  write (io_unit%log,'(1x,a,i7,2(1p,3e12.4,2x))') 'patch_mod::init: id, size, pos =', &
641  self%id, self%size, self%position
642  if (io%verbose>3) &
643  write (io_unit%log,*) 'id, LEVEL: ', self%id, self%level, maxval(self%box), minval(self%ds)
644  !$omp end critical (log_cr)
645  end if
646  !-----------------------------------------------------------------------------
647  ! Allocate memory, if not already allocated
648  !-----------------------------------------------------------------------------
649  if (.not.self%mem_allocated) then
650  allocate (self%t(self%nt), self%dt(self%nt), self%iit(self%nt))
651 !print '(a,3(2x,3i4))', &
652 !'patch_t: gn, nv, nt, nw =', self%gn, self%mesh%gn, self%nv, self%nt, self%nw
653  allocate (self%mem(self%mesh(1)%gn, &
654  self%mesh(2)%gn, &
655  self%mesh(3)%gn, &
656  self%nv,self%nt,self%nw))
657  allocate (self%mem_lock(self%nt))
658  do i=1,self%nt
659  write (kind,'("mem",i1)') i
660  call self%mem_lock(i)%init (kind)
661  end do
662  self%mem_allocated = .true. ! mark
663  call io%bits_mem (storage_size(self%mem), product(shape(self%mem)), 'mem')
664  allocate (self%unsigned (self%nv))
665  allocate (self%pervolume(self%nv))
666  self%unsigned (:) = .false.
667  self%pervolume(:) = .false.
668  end if
669  self%dt = 0.0
670  self%t = -1.0 ! mark as "no guard zones"
671  self%time = -1.0 ! mark as "no guard zones"
672  self%iit = 1 ! initial memory slots
673  self%it = 1 ! short-hand (=1 initially)
674  self%new = min(2,self%nt) ! slot for updates
675  self%iit(self%nt) = self%new ! first output slot
676  self%new = self%iit(self%nt) ! short-hand
677  self%t(self%it) = 0.0 ! the only valid slot (no. 1)
678  self%wallclock = wallclock() ! starting time
679  n_header = storage_size(header)/32 ! 32 bits per word
680  if (n_header*4 /= storage_size(header)/8) then
681  write(io%output,*) n_header*4, storage_size(header)/8
682  call mpi%abort ('n_header is incorrect')
683  end if
684  if (all(self%size == self%box)) then
685  call self%set (bits%root_grid)
686  if (io%verbose > 1) write(io%output,*) 'set bits%root_grid id =', self%id
687  end if
688  call trace%end()
689 END SUBROUTINE setup
690 
691 !===============================================================================
692 !> Deallocate a patch and its allocated parts
693 !===============================================================================
694 SUBROUTINE dealloc (self)
695  class(patch_t):: self
696  integer :: i
697  !-----------------------------------------------------------------------------
698  call trace%begin ('patch_t%dealloc')
699  if (associated(self%mesh)) &
700  call meshrecycler (self%mesh)
701  call io%bits_mem (-storage_size(self%mem), product(shape(self%mem)),'-mem')
702  if (associated(self%mem)) deallocate (self%mem)
703  nullify (self%mem)
704  self%mem_allocated = .false.
705  if (allocated(self%force_per_unit_mass)) then
706  call io%bits_mem (-storage_size(self%force_per_unit_mass), &
707  product(shape(self%force_per_unit_mass)),'-fpm')
708  deallocate (self%force_per_unit_mass)
709  end if
710  if (allocated(self%force_per_unit_volume)) then
711  call io%bits_mem (-storage_size(self%force_per_unit_volume), &
712  product(shape(self%force_per_unit_volume)),'-fpv')
713  deallocate (self%force_per_unit_volume)
714  end if
715  if (allocated(self%heating_per_unit_mass)) then
716  call io%bits_mem (-storage_size(self%heating_per_unit_mass), &
717  product(shape(self%heating_per_unit_mass)),'-qpm')
718  deallocate (self%heating_per_unit_mass)
719  end if
720  if (allocated(self%heating_per_unit_volume)) then
721  call io%bits_mem (-storage_size(self%heating_per_unit_volume), &
722  product(shape(self%heating_per_unit_volume)),'-qpv')
723  deallocate (self%heating_per_unit_volume)
724  end if
725  if (associated(self%vgas)) then
726  call io%bits_mem (-storage_size(self%vgas), &
727  product(shape(self%vgas)),'-vgas')
728  deallocate (self%vgas)
729  end if
730  if (associated(self%density)) then
731  call io%bits_mem (-storage_size(self%density), &
732  product(shape(self%density)),'-density')
733  deallocate (self%density)
734  end if
735  if (associated(self%pressure)) then
736  call io%bits_mem (-storage_size(self%pressure), &
737  product(shape(self%pressure)),'-pressure')
738  deallocate (self%pressure)
739  end if
740  if (allocated(self%mask)) then
741  call io%bits_mem (-storage_size(self%mask), &
742  product(shape(self%mask)), '-mask')
743  deallocate (self%mask)
744  end if
745  if (allocated(self%irefine)) then
746  call io%bits_mem (-storage_size(self%irefine), &
747  product(shape(self%irefine)), '-irefine')
748  deallocate (self%irefine)
749  end if
750  call self%task_t%dealloc
751  call trace%end()
752 END SUBROUTINE dealloc
753 
754 SUBROUTINE clone_mem_accounting (self)
755  class(patch_t):: self
756  if (allocated(self%mask)) &
757  call io%bits_mem (storage_size(self%mask), &
758  product(shape(self%mask)), 'mask')
759  if (allocated(self%irefine)) &
760  call io%bits_mem (storage_size(self%irefine), &
761  product(shape(self%irefine)), 'irefine')
762  if (allocated(self%force_per_unit_mass)) &
763  call io%bits_mem (storage_size(self%force_per_unit_mass), &
764  product(shape(self%force_per_unit_mass)),'fpm')
765  if (allocated(self%force_per_unit_volume)) &
766  call io%bits_mem (storage_size(self%force_per_unit_volume), &
767  product(shape(self%force_per_unit_volume)),'fpv')
768  if (allocated(self%heating_per_unit_mass)) &
769  call io%bits_mem (storage_size(self%heating_per_unit_mass), &
770  product(shape(self%heating_per_unit_mass)),'qpm')
771  if (allocated(self%heating_per_unit_volume)) &
772  call io%bits_mem (storage_size(self%heating_per_unit_volume), &
773  product(shape(self%heating_per_unit_volume)),'qpv')
774 END SUBROUTINE clone_mem_accounting
775 
776 !===============================================================================
777 !> Rotate time slots. The initial conditions are in slot 1, and the first time
778 !> step puts new values in the 'new' slot 2, while saving the time step used in
779 !> dt(1). Then the current slot (it) becomes 2, and the new one becomes 3, etc.
780 !> This way, there is no need to copy memory btw time steps.
781 !===============================================================================
782 SUBROUTINE rotate (self)
783  class(patch_t):: self
784  integer:: iv, nv, new
785  real(8):: time
786  integer, save:: itimer=0
787  !-----------------------------------------------------------------------------
788  ! Periodic wrapping, time extrapolation
789  !-----------------------------------------------------------------------------
790  if (self%rotated) return
791  call trace%begin ('patch_t%rotate', itimer=itimer)
792  if (self%is_periodic()) then
793  call self%periodic_grid ! periodic root grid
794  else if (extrapolate_gz) then
795  call self%extrapolate_guards ! time-extrapolate guard cells
796  end if
797  !-----------------------------------------------------------------------------
798  ! One line summary
799  !-----------------------------------------------------------------------------
800  if (io%verbose>1) then
801  write (io_unit%log,'(a,i6,2i3,2e13.6,2x,a,1p,6e10.2)') 'update: id,j,k,l,n,dt,t', &
802  self%id, self%it, self%new, self%dtime, self%time, 'minmax:', &
803  minval(self%mem(:,:,:,1,self%it ,self%nw)), &
804  maxval(self%mem(:,:,:,1,self%it ,self%nw)), &
805  minval(self%mem(:,:,:,1,self%new,1)), &
806  maxval(self%mem(:,:,:,1,self%new,1))
807  end if
808  !-----------------------------------------------------------------------------
809  ! Update patch position for non-zero patch velocities;
810  ! This happens *before* `rotate` because `dtime` stores the time step just taken.
811  !-----------------------------------------------------------------------------
812  call self%update_position
813  !-----------------------------------------------------------------------------
814  ! Task time slot rotation
815  !-----------------------------------------------------------------------------
816  call self%task_t%rotate
817  call trace%end (itimer)
818 END SUBROUTINE rotate
819 
820 !===============================================================================
821 !> This function gives the physical position of an index triple (patch point)
822 !> in absolute Cartesian coords.
823 !===============================================================================
824 FUNCTION myposition (self, ii, iv, pp) RESULT(out)
825  class(patch_t) :: self
826  integer, intent(in) :: ii(3)
827  integer, optional :: iv
828  real, optional :: pp(3)
829  !.............................................................................
830  real(8) :: out(3), cellpos(3), h(3)
831  !.............................................................................
832  if (present(pp)) then
833  cellpos = self%position + (real(ii)-self%offset + pp)*self%ds
834  else
835  cellpos = self%position + (real(ii)-self%offset)*self%ds
836  end if
837  if (present(iv)) then
838  h = self%index_stagger(iv)
839  cellpos = cellpos + h*self%mesh%d
840  end if
841  if (self%id==io%id_debug .and. io%verbose>1) &
842  write(io%output,'(a,3i4,4(3x,3f10.2))') &
843  'mypos:',ii,self%position,self%offset,self%ds,cellpos
844  ! convert to Cartesian coords.
845  if (self%mesh_type == mesh_types%cylindrical) then
846  cellpos = cylindricaltocartesian(cellpos(1), cellpos(2), cellpos(3))
847  else if (self%mesh_type == mesh_types%spherical) then
848  cellpos = sphericaltocartesian(cellpos(1), cellpos(2), cellpos(3))
849  end if
850  out = cellpos ! absolute position
851 END FUNCTION myposition
852 
853 !===============================================================================
854 !> Scalar version of the previous function, but **including staggering factor**.
855 !> For the conversion to Cartesian coords., all three indices are required!
856 !===============================================================================
857 FUNCTION myposition_staggered (self, ii, idir, h) RESULT(out)
858  class(patch_t):: self
859  integer, intent(in):: ii(3), idir
860  real(8):: out, cellpos(3)
861  real(8), optional :: h
862  real(4) :: hlocal
863  !.............................................................................
864  if (present(h)) then
865  hlocal = h
866  else
867  hlocal = 0.0
868  end if
869 
870  cellpos = (real(ii-self%li,kind=8) + merge(0.5 + hlocal,0.0,self%n(idir)>1)) &
871  * self%ds + self%llc_nat ! in patch coords.
872  ! convert to Cartesian coords.
873  if (self%mesh_type == mesh_types%cylindrical) then
874  cellpos = cylindricaltocartesian(cellpos(1), cellpos(2), cellpos(3))
875  else if (self%mesh_type == mesh_types%spherical) then
876  cellpos = sphericaltocartesian(cellpos(1), cellpos(2), cellpos(3))
877  end if
878  out = cellpos(idir) ! absolute position
879 
880 END FUNCTION myposition_staggered
881 
882 !===============================================================================
883 FUNCTION myfloatindex_scalar (self, cpos, idir, h) RESULT(out)
884  class(patch_t) :: self
885  real(8), intent(in) :: cpos(3)
886  integer, intent(in) :: idir
887  real(8), intent(in) :: h
888  real(8) :: pp(3), ppCart(3), out, intermezzo
889  !.............................................................................
890  ! convert to native patch coords.
891  if (self%mesh_type == mesh_types%Cartesian) then
892  ! do nothing
893  pp = cpos
894  else if (self%mesh_type == mesh_types%spherical) then
895  ppcart = cpos - self%position ! Cart. coords to centre of patch's coord. system
896  pp = cartesiantospherical(ppcart(1), ppcart(2), ppcart(3))
897  else if (self%mesh_type == mesh_types%cylindrical) then
898  ppcart = cpos - self%position ! Cart. coords to centre of patch's coord. system
899  pp = cartesiantocylindrical(ppcart(1), ppcart(2), ppcart(3))
900  end if
901 
902  intermezzo = (pp(idir) - self%llc_nat(idir)) / self%ds(idir) + self%li(idir) &
903  - merge(0.5d0 + h,0.0d0,self%n(idir)>1)
904  out = intermezzo
905 
906 END FUNCTION myfloatindex_scalar
907 
908 !===============================================================================
909 !> Periodic mapping of a grid -- typically the root grid. This is typically
910 !> only used in tests, so does not need to be efficient.
911 !===============================================================================
912 SUBROUTINE periodic_grid (self, only)
913  class(patch_t):: self
914  integer, optional:: only
915  !.............................................................................
916  integer:: ix, iy, iz, i(3), j(3), lb(3), ub(3), li(3), ui(3), n(3)
917  logical:: periodic(3)
918  !-----------------------------------------------------------------------------
919  call trace%begin ('patch_t%periodic_grid')
920  periodic = self%periodic .and. self%size+self%ds > self%box
921  if (.not.any(periodic)) call mpi%abort ('patch_t%periodic_grid')
922  lb = self%mesh%lb
923  ub = self%mesh%ub
924  li = self%mesh%li
925  ui = self%mesh%ui
926  where (self%mesh%lower_boundary)
927  li = self%mesh%lb
928  end where
929  where (self%mesh%upper_boundary)
930  ui = self%mesh%ub
931  end where
932  n = self%ncell
933  do iz=lb(3),ub(3)
934  do iy=lb(2),ub(2)
935  do ix=lb(1),ub(1)
936  i = [ix,iy,iz]
937  if (all(i >= li .and. i <= ui)) cycle
938  where (n > 0)
939  j = modulo(i-li,n)+li
940  else where
941  j = 1
942  end where
943  !---------------------------------------------------------------------------
944  ! Prevent wrapping in non-periodic directions
945  !---------------------------------------------------------------------------
946  where (.not.periodic)
947  j = i
948  end where
949  !---------------------------------------------------------------------------
950  ! Make sure horizontal periodic wrapping continues over the whole vertical
951  ! range when physical boundary conditions are present. Note the use of
952  ! mesh%ui and mesh%li below; this is because the physical boundaries are
953  ! inside the normal patch boundaries.
954  !---------------------------------------------------------------------------
955  where (self%mesh%lower_boundary .and. i < self%mesh%li+4)
956  j = i
957  end where
958  where (self%mesh%upper_boundary .and. i > self%mesh%ui-4)
959  j = i
960  end where
961  if (present(only)) then
962  self%mem(ix,iy,iz,only,self%it ,1) = self%mem(j(1),j(2),j(3),only,self%it ,1)
963  self%mem(ix,iy,iz,only,self%new,1) = self%mem(j(1),j(2),j(3),only,self%new,1)
964  else
965  self%mem(ix,iy,iz,:,self%it ,1) = self%mem(j(1),j(2),j(3),:,self%it ,1)
966  self%mem(ix,iy,iz,:,self%new,1) = self%mem(j(1),j(2),j(3),:,self%new,1)
967  end if
968  end do
969  end do
970  end do
971  if (io%verbose > 1) write(io%output,*) &
972  'periodic_grid: id, it, new =', self%id, self%it, self%new
973  call trace%end()
974 END SUBROUTINE periodic_grid
975 
976 !===============================================================================
977 FUNCTION is_periodic(self) RESULT(out)
978  class(patch_t):: self
979  logical:: out
980  !.............................................................................
981  out = any(self%periodic .and. self%size+self%ds > self%box)
982 END FUNCTION is_periodic
983 
984 !=======================================================================
985 !> Optimized periodic boundary enforcer
986 !=======================================================================
987 SUBROUTINE make_periodic4 (patch, f)
988  class(patch_t):: patch
989  class(mesh_t), pointer:: m(:)
990  real(4):: f(:,:,:)
991  integer:: ix, iy, iz, n(3)
992  logical:: periodic(3)
993  !-----------------------------------------------------------------------------
994  periodic = patch%periodic .and. patch%size+patch%ds > patch%box
995  m => patch%mesh
996  n = patch%ncell
997  !-----------------------------------------------------------------------------
998  ! x periodicity; make periodic inside the yz patch limits
999  !-----------------------------------------------------------------------------
1000  if (n(1) > 1 .and. periodic(1)) then
1001  do iz=m(3)%li,m(3)%ui
1002  do ix=m(1)%lb,m(1)%lo
1003  do iy=m(2)%li,m(2)%ui
1004  f(ix,iy,iz) = f(ix+n(1),iy,iz)
1005  end do
1006  end do
1007  do ix=m(1)%uo,m(1)%ub
1008  do iy=m(2)%li,m(2)%ui
1009  f(ix,iy,iz) = f(ix-n(1),iy,iz)
1010  end do
1011  end do
1012  end do
1013  end if
1014  !-----------------------------------------------------------------------------
1015  ! y periodicity, full range in x now
1016  !-----------------------------------------------------------------------------
1017  if (n(2) > 1 .and. periodic(2)) then
1018  do iz=m(3)%li,m(3)%ui
1019  do iy=m(2)%lb,m(2)%lo
1020  do ix=m(1)%lb,m(1)%ub
1021  f(ix,iy,iz) = f(ix,iy+n(2),iz)
1022  end do
1023  end do
1024  do iy=m(2)%uo,m(2)%ub
1025  do ix=m(1)%lb,m(1)%ub
1026  f(ix,iy,iz) = f(ix,iy-n(2),iz)
1027  end do
1028  end do
1029  end do
1030  end if
1031  !-----------------------------------------------------------------------------
1032  ! z periodicity, full range in xy
1033  !-----------------------------------------------------------------------------
1034  if (n(3) > 1 .and. periodic(3)) then
1035  do iy=m(2)%lb,m(2)%ub
1036  do iz=m(3)%lb,m(3)%lo
1037  do ix=m(1)%lb,m(1)%ub
1038  f(ix,iy,iz) = f(ix,iy,iz+n(3))
1039  end do
1040  end do
1041  do iz=m(3)%uo,m(3)%ub
1042  do ix=m(1)%lb,m(1)%ub
1043  f(ix,iy,iz) = f(ix,iy,iz-n(3))
1044  end do
1045  end do
1046  end do
1047  end if
1048 END SUBROUTINE make_periodic4
1049 
1050 !=======================================================================
1051 !> Optimized periodic boundary enforcer
1052 !=======================================================================
1053 SUBROUTINE make_periodic8 (patch, f)
1054  class(patch_t):: patch
1055  class(mesh_t), pointer:: m(:)
1056  real(8):: f(:,:,:)
1057  integer:: ix, iy, iz, n(3)
1058  logical:: periodic(3)
1059  !-----------------------------------------------------------------------------
1060  periodic = patch%periodic .and. patch%size+patch%ds > patch%box
1061  m => patch%mesh
1062  n = patch%ncell
1063  !-----------------------------------------------------------------------------
1064  ! x periodicity; make periodic inside the yz patch limits
1065  !-----------------------------------------------------------------------------
1066  if (n(1) > 1 .and. periodic(1)) then
1067  do iz=m(3)%li,m(3)%ui
1068  do ix=m(1)%lb,m(1)%lo
1069  do iy=m(2)%li,m(2)%ui
1070  f(ix,iy,iz) = f(ix+n(1),iy,iz)
1071  end do
1072  end do
1073  do ix=m(1)%uo,m(1)%ub
1074  do iy=m(2)%li,m(2)%ui
1075  f(ix,iy,iz) = f(ix-n(1),iy,iz)
1076  end do
1077  end do
1078  end do
1079  end if
1080  !-----------------------------------------------------------------------------
1081  ! y periodicity, full range in x now
1082  !-----------------------------------------------------------------------------
1083  if (n(2) > 1 .and. periodic(2)) then
1084  do iz=m(3)%li,m(3)%ui
1085  do iy=m(2)%lb,m(2)%lo
1086  do ix=m(1)%lb,m(1)%ub
1087  f(ix,iy,iz) = f(ix,iy+n(2),iz)
1088  end do
1089  end do
1090  do iy=m(2)%uo,m(2)%ub
1091  do ix=m(1)%lb,m(1)%ub
1092  f(ix,iy,iz) = f(ix,iy-n(2),iz)
1093  end do
1094  end do
1095  end do
1096  end if
1097  !-----------------------------------------------------------------------------
1098  ! z periodicity, full range in xy
1099  !-----------------------------------------------------------------------------
1100  if (n(3) > 1 .and. periodic(3)) then
1101  do iy=m(2)%lb,m(2)%ub
1102  do iz=m(3)%lb,m(3)%lo
1103  do ix=m(1)%lb,m(1)%ub
1104  f(ix,iy,iz) = f(ix,iy,iz+n(3))
1105  end do
1106  end do
1107  do iz=m(3)%uo,m(3)%ub
1108  do ix=m(1)%lb,m(1)%ub
1109  f(ix,iy,iz) = f(ix,iy,iz-n(3))
1110  end do
1111  end do
1112  end do
1113  end if
1114 END SUBROUTINE make_periodic8
1115 
1116 !===============================================================================
1117 !> Time extrapolate guard cell values, to allow using also guard zone values to
1118 !> get boundary conditions for higher level patches. This is done before time
1119 !> slot rotation, so the slot that has just been updated is iit(nt).
1120 !===============================================================================
1121 SUBROUTINE extrapolate_guards (self)
1122  class(patch_t) :: self
1123  !.............................................................................
1124  integer :: it0, it1, it2, nt, iit(self%nt)
1125  integer, dimension(3) :: lb, lo, li, ub, uo, ui, i
1126  integer :: iv
1127  integer, save :: itimer=0
1128  real(8) :: t(self%nt)
1129  !-----------------------------------------------------------------------------
1130  call trace%begin ('patch_t%extrapolate_guards', 3, itimer=itimer)
1131  nt = self%nt
1132  call self%timeslots (iit, t)
1133  it0 = iit(nt) ! current time slot
1134  it1 = iit(nt-1) ! previous time slot
1135  it2 = iit(nt-2) ! the one before that
1136  !-----------------------------------------------------------------------------
1137  ! Before someone asks: Associate on the non-changing 6 doesn't please gfort
1138  !-----------------------------------------------------------------------------
1139  lb = self%mesh%lb
1140  li = self%mesh%li
1141  lo = self%mesh%lo
1142  ub = self%mesh%ub
1143  ui = self%mesh%ui
1144  uo = self%mesh%uo
1145  where (self%mesh%lower_boundary)
1146  li = self%mesh%lb
1147  end where
1148  where (self%mesh%upper_boundary)
1149  ui = self%mesh%ub
1150  end where
1151  if (io%id_debug>0.and.(self%id==io%id_debug.or.abs(self%id-io%id_debug)==27)) then
1152  write(io%output,'(i5,2x,a,6i7)') self%id, &
1153  'MK extrapolating guard zone values', it2, it1, it0, self%it, self%new
1154  end if
1155  do iv=1,self%nv
1156  if (self%unsigned(iv)) then
1157  !-------------------------------------------------------------------------
1158  ! XY-layers (Z varies the least)
1159  !-------------------------------------------------------------------------
1160  if (self%gn(3) > 1) then
1161  self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it0,1) = exp( &
1162  log(self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it1,1))*2 - &
1163  log(self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it2,1)) )
1164  self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it0,1) = exp( &
1165  log(self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it1,1))*2 - &
1166  log(self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it2,1)) )
1167  end if
1168  !-------------------------------------------------------------------------
1169  ! ZX-layers (Y varies the least)
1170  !-------------------------------------------------------------------------
1171  if (self%gn(2) > 1) then
1172  self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it0,1) = exp( &
1173  log(self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it1,1))*2 - &
1174  log(self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it2,1)) )
1175  self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it0,1) = exp( &
1176  log(self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it1,1))*2 - &
1177  log(self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it2,1)) )
1178  end if
1179  !-------------------------------------------------------------------------
1180  ! YZ-layers (X varies the least)
1181  !-------------------------------------------------------------------------
1182  if (self%gn(1) > 1) then
1183  self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = exp( &
1184  log(self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it1,1))*2 - &
1185  log(self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it2,1)) )
1186  self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = exp( &
1187  log(self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it1,1))*2 - &
1188  log(self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it2,1)) )
1189  end if
1190  else
1191  !-------------------------------------------------------------------------
1192  ! XY-layers (Z varies the least)
1193  !-------------------------------------------------------------------------
1194  if (self%gn(3) > 1) then
1195  self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it0,1) = &
1196  self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it1,1)*2 - &
1197  self%mem(lb(1):ub(1),lb(2):ub(2),lb(3):lo(3),iv,it2,1)
1198  self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it0,1) = &
1199  self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it1,1)*2 - &
1200  self%mem(lb(1):ub(1),lb(2):ub(2),uo(3):ub(3),iv,it2,1)
1201  end if
1202  !-------------------------------------------------------------------------
1203  ! ZX-layers (Y varies the least)
1204  !-------------------------------------------------------------------------
1205  if (self%gn(2) > 1) then
1206  self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it0,1) = &
1207  self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it1,1)*2 - &
1208  self%mem(lb(1):ub(1),lb(2):lo(2),li(3):ui(3),iv,it2,1)
1209  self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it0,1) = &
1210  self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it1,1)*2 - &
1211  self%mem(lb(1):ub(1),uo(2):ub(2),li(3):ui(3),iv,it2,1)
1212  end if
1213  !-------------------------------------------------------------------------
1214  ! YZ-layers (X varies the least)
1215  !-------------------------------------------------------------------------
1216  if (self%gn(1) > 1) then
1217  self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = &
1218  self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it1,1)*2 - &
1219  self%mem(lb(1):lo(1),li(2):ui(2),li(3):ui(3),iv,it2,1)
1220  self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it0,1) = &
1221  self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it1,1)*2 - &
1222  self%mem(uo(1):ub(1),li(2):ui(2),li(3):ui(3),iv,it2,1)
1223  end if
1224  end if
1225  end do
1226  call trace%end (itimer)
1227 END SUBROUTINE extrapolate_guards
1228 
1229 !===============================================================================
1230 !> Pack patch info into a buffer, for MPI communication
1231 !===============================================================================
1232 SUBROUTINE allocate_mesg (self)
1233  class(patch_t):: self
1234  !-----------------------------------------------------------------------------
1235  allocate (self%mesg)
1236  self%mesg%nbuf = product(self%gn) ! 4-bytes per element; `anonymous_copy` also uses 4-bytes.
1237  if (kind(self%mem)==8) &
1238  self%mesg%nbuf = self%mesg%nbuf * 2 ! 8-bytes per element, but `anonymous_copy` uses 4-bytes!
1239  self%mesg%nbuf = n_header + self%nv*self%mesg%nbuf
1240  allocate (self%mesg%buffer(self%mesg%nbuf))
1241  call io%bits_mem (storage_size(self%mesg%buffer), product(shape(self%mesg%buffer)))
1242 END SUBROUTINE allocate_mesg
1243 
1244 !===============================================================================
1245 !> Pack patch info into a buffer, for MPI communication
1246 !===============================================================================
1247 SUBROUTINE pack (self, mesg)
1248  class(patch_t):: self
1249  class(mesg_t), pointer:: mesg
1250  !.............................................................................
1251  type(header_t):: header
1252  integer:: n_buf, ibuf, iv, it
1253  real:: pt
1254  integer, save:: itimer=0
1255  !-----------------------------------------------------------------------------
1256  ! Compute size of buffer and allocate -- it will be freed by the writer
1257  !-----------------------------------------------------------------------------
1258  call trace%begin ('patch_t%pack', itimer=itimer)
1259  allocate (mesg)
1260  n_buf = product(self%gn) ! 4-bytes per element; `anonymous_copy` also uses 4-bytes.
1261  if (kind(self%mem)==8) &
1262  n_buf = n_buf * 2 ! 8-bytes per element, but `anonymous_copy` uses 4-bytes!
1263  mesg%nbuf = n_header + self%nv*n_buf
1264  if (self%is_set(bits%swap_request) .and. self%is_set(bits%boundary)) then
1265  mesg%nbuf = mesg%nbuf+ self%nv*n_buf
1266  end if
1267  allocate (mesg%buffer(mesg%nbuf))
1268  if (mpi_mesg%nbuf==0) then
1269  !$omp critical (nbuf_cr)
1270  if (mpi_mesg%nbuf==0) then
1271  mpi_mesg%nbuf = mesg%nbuf
1272  write (io_unit%log,*) &
1273  'patch_t%pack: setting mpi_mesg%nbuf =', mpi_mesg%nbuf
1274  end if
1275  !$omp end critical (nbuf_cr)
1276  end if
1277  call io%bits_mem (storage_size(mesg%buffer), product(shape(mesg%buffer)))
1278  allocate (mesg%reqs(self%n_nbors))
1279  mesg%id = self%id ! MPI message tag
1280  mesg%nreq = 0 ! reset counter
1281  !-----------------------------------------------------------------------------
1282  ! Copy relevant patch info to sequenced header, and copy that to the buffer
1283  !-----------------------------------------------------------------------------
1284  if (self%debug(2)) then
1285  write (io_unit%log,1) wallclock(), &
1286  'DBG patch_t%pack, id, nq, time, status:', &
1287  self%id, self%nq, self%time, &
1288  self%is_set (bits%internal), &
1289  self%is_set (bits%boundary), &
1290  self%is_set (bits%virtual), &
1291  self%is_set (bits%external), &
1292  self%is_set (bits%swap_request)
1293  1 format(f12.6,2x,a,2i9,f12.6,2x,5l1)
1294  flush (io_unit%log)
1295  end if
1296  call patch_to_header (self, header)
1297  ibuf = 1
1298  call anonymous_copy (n_header, header, mesg%buffer(ibuf))
1299  ibuf = ibuf + n_header
1300  !-----------------------------------------------------------------------------
1301  ! Copy the variables to the output buffer. If the bdry+swap bit are set,
1302  ! copy as a minimum also the previous time slot
1303  !-----------------------------------------------------------------------------
1304  do iv=1,self%nv
1305  call anonymous_copy (n_buf, self%mem(:,:,:,iv,self%it,1), mesg%buffer(ibuf))
1306  ibuf = ibuf + n_buf
1307  end do
1308  if (self%is_set(bits%swap_request) .and. self%is_set(bits%boundary)) then
1309  it = mod(self%it-1+self%nt-1,self%nt)+1
1310  if (io%verbose>1) then
1311  write (io_unit%log,*) self%id, 'extra time slot packed', it
1312  flush (io_unit%log)
1313  end if
1314  do iv=1,self%nv
1315  call anonymous_copy (n_buf, self%mem(:,:,:,iv,it,1), mesg%buffer(ibuf))
1316  ibuf = ibuf + n_buf
1317  end do
1318  end if
1319  if (io%verbose > 2) &
1320  write (io_unit%log,*) self%id, 'after pack', minval(self%mem(:,:,:,1,self%it,1))
1321  call trace%end (itimer)
1322 END SUBROUTINE pack
1323 
1324 !===============================================================================
1325 !> Unpack patch id from an MPI patch_t buffer
1326 !===============================================================================
1327 INTEGER FUNCTION unpack_id (mesg)
1328  class(mesg_t), pointer:: mesg
1329  !.............................................................................
1330  type(header_t):: header
1331  !-----------------------------------------------------------------------------
1332  call trace%begin ('patch_t%unpack_id')
1333  call anonymous_copy (n_header, mesg%buffer, header)
1334  unpack_id = header%id
1335  call trace%end ()
1336 END FUNCTION unpack_id
1337 
1338 !===============================================================================
1339 !> Unpack patch info from a buffer, for MPI communication
1340 !===============================================================================
1341 SUBROUTINE unpack (self, mesg)
1342  class(patch_t):: self
1343  class(mesg_t), pointer:: mesg
1344  !.............................................................................
1345  type(header_t):: header
1346  integer:: n_buf, ibuf, iv, it
1347  integer, dimension(:), pointer:: buffer
1348  integer, save:: itimer=0
1349  !-----------------------------------------------------------------------------
1350  ! Compute size of buffer
1351  !-----------------------------------------------------------------------------
1352  call trace%begin ('patch_t%unpack', itimer=itimer)
1353  call self%set (bits%busy)
1354  self%unpack_time = wallclock()
1355  buffer => mesg%buffer
1356  n_buf = product(self%gn) ! 4-bytes per element; `anonymous_copy` also uses 4-bytes.
1357  if (kind(self%mem)==8) &
1358  n_buf = n_buf * 2 ! 8-bytes per element, but `anonymous_copy` uses 4-bytes!
1359  !-----------------------------------------------------------------------------
1360  ! Copy relevant patch info from sequenced header
1361  !-----------------------------------------------------------------------------
1362  ibuf = 1
1363  call anonymous_copy (n_header, buffer(ibuf), header)
1364  ibuf = ibuf + n_header
1365  !-----------------------------------------------------------------------------
1366  ! Lock the %new slot, which will become the %it slot after header_to_patch
1367  !-----------------------------------------------------------------------------
1368  call self%header_to_patch (header)
1369  if (io%verbose > 1) &
1370  write (io_unit%log,'(f12.6,3x,a,i7,2i4,2x,3l1,2i4,f8.4,f11.6)') &
1371  wallclock(), 'header_to_patch: id, nq, BVS, nsnd, nrcv, latency, time =', &
1372  self%id, mesg%sender, self%nq, &
1373  self%is_set(bits%boundary), &
1374  self%is_set(bits%virtual), &
1375  self%is_set(bits%swap_request), &
1376  mpi_mesg%sent_list%n, mpi_mesg%recv_list%n, self%latency, self%time
1377  self%wc_last = wallclock()
1378  if (io%verbose > 1) write (io_unit%log,1) 'patch_mod::unpk id, nq, time after:', &
1379  self%id, self%nq, self%time, &
1380  self%is_set (bits%internal), &
1381  self%is_set (bits%boundary), &
1382  self%is_set (bits%virtual), &
1383  self%is_set (bits%external), &
1384  self%is_set (bits%swap_request), &
1385  self%wc_last
1386  1 format(a,2i9,f12.6,2x,5l1,f12.1)
1387  !-----------------------------------------------------------------------------
1388  ! Copy the variables from the output buffer. If this is swap of patch roles
1389  ! expect at least one more time slot.
1390  !-----------------------------------------------------------------------------
1391  call self%mem_lock(self%it)%set ('unpack')
1392  do iv=1,self%nv
1393  call anonymous_copy (n_buf, buffer(ibuf), self%mem(:,:,:,iv,self%it,1))
1394 !write(io_unit%log,'(a,i6,2i4,1p,e14.5)') &
1395 !'patch_t%unpack: id, iv, it, it, minval =', &
1396 !self%id, iv, self%it, minval(self%mem(:,:,:,iv,it,1))
1397  ibuf = ibuf + n_buf
1398  end do
1399  call self%mem_lock(self%it)%unset ('unpack')
1400  if (self%is_set(bits%swap_request) .and. self%is_set(bits%virtual)) then
1401  it = mod(self%it-1+self%nt-1,self%nt)+1
1402  call self%mem_lock(it)%set ('unpack')
1403  if (io%verbose>0) then
1404  write (io_unit%log,*) self%id, 'extra time slot unpacked', it
1405  flush (io_unit%log)
1406  end if
1407  do iv=1,self%nv
1408  call anonymous_copy (n_buf, buffer(ibuf), self%mem(:,:,:,iv,it,1))
1409  ibuf = ibuf + n_buf
1410  end do
1411  call self%mem_lock(it)%unset ('unpack')
1412  end if
1413  if (io%verbose > 2) &
1414  write (io_unit%log,*) self%id, 'after unpack', minval(self%mem(:,:,:,1,self%it,1))
1415  call self%clear (bits%busy)
1416  call trace%end (itimer)
1417 END SUBROUTINE unpack
1418 
1419 !===============================================================================
1420 !> Print info to stdout. If
1421 !> verbose==0: print only id=1 to stdout and io_unit%log
1422 !> verbose==1: print only id=1 to stdout but all to io_unit%log
1423 !> verbose==2: print iv=1 to stdout and io_unit%log for all id
1424 !> verbose==3: print iv=1,nv to stdout and io_unit%log for all id
1425 !===============================================================================
1426 SUBROUTINE info (self, nq, ntask, experiment_name)
1427  class(patch_t):: self
1428  integer, optional:: nq, ntask
1429  character(len=64), optional:: experiment_name
1430  !.............................................................................
1431  integer:: iv, nv, iprint_l
1432  integer:: it0, it1, it
1433  real:: fmin, fmax
1434  real, dimension(:,:,:), pointer:: df
1435  integer, parameter:: max_lines=50
1436  integer, save:: counts(6)=0, id_prv=0
1437  integer:: counts_l(6), i
1438  character(len=120):: fmt
1439  real(8):: time, print_next
1440  logical:: print_it
1441  type(lock_t), save:: lock
1442  integer, save:: itimer=0
1443  !-----------------------------------------------------------------------------
1444  call trace%begin ('patch_t%info', itimer=itimer)
1445  if (io_unit%do_validate) then
1446  it0 = merge(self%iit(self%nt-2), 1, self%nt>2)
1447  it1 = merge(self%iit(self%nt-1), 1, self%nt>1)
1448  if (self%nw==1) then
1449  allocate (df(size(self%mem,1),size(self%mem,2),size(self%mem,3)))
1450  end if
1451  do iv=1,self%nv
1452  if (self%nw==1) then
1453  df = self%mem(:,:,:,iv,it1,1) &
1454  - self%mem(:,:,:,iv,it0,1)
1455  fmin = self%fminval(df)
1456  fmax = self%fmaxval(df)
1457  else
1458  fmin = self%fminval(self%mem(:,:,:,iv,it0,self%nw))
1459  fmax = self%fmaxval(self%mem(:,:,:,iv,it0,self%nw))
1460  end if
1461  write(io%output,2) &
1462  'upd',self%id, self%level, iv, self%istep, self%time-self%dtime, &
1463  self%dtime, self%u_max, fmin, fmax, &
1464  self%fminval(self%mem(:,:,:,iv,it0,1)), &
1465  self%fmaxval(self%mem(:,:,:,iv,it0,1)), &
1466  real(timer%n_update), ntask, nq, self%is_set(bits%boundary), self%is_set(bits%virtual)
1467  2 format(a,i8,2i4,i3,1p,e14.5,e12.4,0p,f7.2,2x,1p,2e12.4,1x,2e12.4,2x,e10.2,2i5,l3,l1,o12)
1468  end do
1469  if (self%nw==1) then
1470  deallocate (df)
1471  end if
1472  call trace%end (itimer)
1473  return
1474  end if
1475  1 format(a,i8,2i4,i3,f12.6,1p,e12.4,0p,f7.2,2x,1p,2e10.2,2x,2e11.3,2x,e10.2,9i5,2i4,l3,l1,o12)
1476  if (self%time < 1d4) then
1477  fmt = '(a,i8,2i4,i3,f12.6,1p,e12.4,0p,f7.2,2x,1p,2e10.2,2x,2e11.3,2x,e10.2,i6,3i5,i6,4i5,2i4,l3,l1,o12)'
1478  else
1479  fmt = '(a,i8,2i4,i3,1p,2e12.5, 0p,f7.2,2x,1p,2e10.2,2x,2e11.3,2x,e10.2,i6,3i5,i6,4i5,2i4,l3,l1,o12)'
1480  end if
1481  !-----------------------------------------------------------------------------
1482  if ((io%verbose>=0).and.(self%is_clear(bits%no_density))) then
1483  !$omp atomic
1484  io%dmin = min(io%dmin, real(self%fminval(self%mem(:,:,:,self%idx%d,self%it,1)),kind=4))
1485  !$omp atomic
1486  io%dmax = max(io%dmax, real(self%fmaxval(self%mem(:,:,:,self%idx%d,self%it,1)),kind=4))
1487  !---------------------------------------------------------------------------
1488  ! After coming back to the saved ID, reset it, so we stop searching
1489  !---------------------------------------------------------------------------
1490  if (self%id == id_prv) then
1491  !$omp atomic write
1492  id_prv = 0
1493  end if
1494  !---------------------------------------------------------------------------
1495  ! Allow print every print_every update, if print_every > 0
1496  !---------------------------------------------------------------------------
1497  if (print_every > 0) then
1498  !$omp atomic read
1499  iprint_l = iprint
1500  if (iprint_l >= print_every .and. print_every > 0) then
1501  !$omp atomic write
1502  iprint = 1
1503  print_it = .true.
1504  else
1505  !$omp atomic update
1506  iprint = iprint+1
1507  print_it = .false.
1508  end if
1509  !---------------------------------------------------------------------------
1510  ! If print_every==0, then check if print_time controls cadence, using
1511  ! an if-lock-if sequence to ensure only one thread prints and updates
1512  !---------------------------------------------------------------------------
1513  else if (io%print_time > 0.0) then
1514  if (print_every == 0 .and. self%time >= io%print_next) then
1515  call lock%set ('print')
1516  if (self%time >= io%print_next) then
1517  print_it = .true.
1518  print_next = (nint(self%time/io%print_time)+1)*io%print_time
1519  io%print_next = print_next
1520  end if
1521  call lock%unset ('print')
1522  end if
1523  !---------------------------------------------------------------------------
1524  ! Otherwise, fall back on printing the task update with the smallest dt
1525  !---------------------------------------------------------------------------
1526  else if (self%id==io%id) then
1527  print_it = .true.
1528  end if
1529  if (print_it .or. io%verbose>0) then
1530  nv = 1
1531  if (io%verbose>2) nv = self%nv
1532  !-------------------------------------------------------------------------
1533  counts_l = &
1534  [mpi_mesg%n_check, mpi_mesg%n_ready, mpi_mesg%n_update, &
1535  mpi_mesg%n_send , mpi_mesg%n_recv , mpi_mesg%n_unpk ]
1536  do i=1,size(counts)
1537  !$omp atomic update
1538  counts(i) = counts(i) + counts_l(i)
1539  end do
1540  !$omp atomic update
1541  timer%n_lines = timer%n_lines-1
1542  if (timer%n_lines==0) then
1543  !$omp critical (info_cr2)
1544  !$omp atomic write
1545  timer%n_lines = max_lines
1546  !$omp end atomic
1547  write(io%output,'(a)') &
1548  '0....+....1....+....2....+....3....+....4....+....5....+....' &
1549  //'6....+....7....+....8....+....9....+....0....+....1 ntsk nq ' &
1550  //' chk rdy upd snt prb unp snq rcq unq BV'
1551  !$omp end critical (info_cr2)
1552  end if
1553  !-----------------------------------------------------------------------
1554  ! time-dtime is the task time when update started; it is guaranteed to
1555  ! be increasing monotomically, since the ready queue is time sorted.
1556  ! To allow using the column for plotting, we show this time when all
1557  ! times are printed.
1558  !-----------------------------------------------------------------------
1559  if (io%print_time == 0.0 .and. io%verbose > 0) then
1560  time = self%time - self%dtime
1561  else
1562  time = self%time
1563  end if
1564  do iv=1,nv
1565  write(io%output,fmt) &
1566  'upd',self%id, self%level, iv, omp_mythread, time, &
1567  self%dtime, log10(max(self%u_max,1e-30)), &
1568  io%dmin, io%dmax, &
1569  self%fminval(self%mem(:,:,:,iv,self%it ,1)), &
1570  self%fmaxval(self%mem(:,:,:,iv,self%it ,1)), &
1571  real(timer%n_update), ntask, nq, counts, &
1572  mpi_mesg%nq_send, mpi_mesg%nq_recv, mpi_mesg%unpk_list%n, &
1573  self%is_set(bits%boundary), self%is_set(bits%virtual)
1574  end do
1575  flush (io%output)
1576  !$omp atomic write
1577  io%dmin = huge(1.)
1578  !$omp atomic write
1579  io%dmax = tiny(1.)
1580  !$omp atomic write
1581  io%dtime = huge(1d0)
1582  !$omp atomic write
1583  id_prv = io%id
1584  counts(:) = 0
1585  !-----------------------------------------------------------------------
1586  ! Reset the message counters after every timestep, if relevant
1587  !-----------------------------------------------------------------------
1588  mpi_mesg%n_check=0; mpi_mesg%n_ready=0; mpi_mesg%n_update=0
1589  mpi_mesg%n_send=0 ; mpi_mesg%n_recv=0 ; mpi_mesg%n_unpk=0
1590  mpi_mesg%n_delay=0
1591  end if
1592  !---------------------------------------------------------------------------
1593  ! Keep track of the shortest times step between updates of the io%id local
1594  ! task. Use this task in patch%info, rather than the one with ip==1. As
1595  ! long as shorter dtime are being found, io%id will not match the next self%id.
1596  ! But if the same io%id survices one update cycle, it will match.
1597  !---------------------------------------------------------------------------
1598  if (self%dtime < io%dtime .and. id_prv /= 0) then
1599  !$omp atomic write
1600  io%id = self%id
1601  !$omp atomic write
1602  io%dtime = self%dtime
1603  end if
1604  end if
1605  call trace%end (itimer)
1606 END SUBROUTINE info
1607 
1608 !===============================================================================
1609 !> Return true if a point is inside the self interior
1610 !===============================================================================
1611 LOGICAL FUNCTION contains1 (self, point)
1612  class(patch_t):: self
1613  real(8):: point(3)
1614  contains1 = all(abs(point-self%position) < 0.5_8*self%size)
1615 END FUNCTION contains1
1616 
1617 !===============================================================================
1618 !> Return true if a patch position is inside the self interior
1619 !===============================================================================
1620 LOGICAL FUNCTION contains2 (self, patch)
1621  class(patch_t):: self, patch
1622  contains2 = contains1(self, patch%position)
1623 END FUNCTION contains2
1624 
1625 !===============================================================================
1626 !> Check if there is overlap, using the size plus an extra cell margin. The
1627 !> definition of overlap MUST be symmetric, since it is being used to define
1628 !> neighbor relations. This function overloads task_t%overlaps, in order to
1629 !> force use of the patch_t%distance function.
1630 !===============================================================================
1631 FUNCTION overlaps (self, task)
1632  class(patch_t):: self
1633  class(task_t):: task
1634  logical:: overlaps
1635  real(8):: distance, box, limit
1636  integer:: i
1637  !.............................................................................
1638  !call trace%begin ('patch_t%overlaps', 2)
1639  select type (task)
1640  class is (patch_t)
1641  overlaps = .true.
1642  do i=1,3
1643  distance = self%centre_nat(i) - task%centre_nat(i)
1644  box = self%box(i)
1645  if (self%periodic(i)) then
1646  distance = modulo(distance+0.5_8*box, box) - 0.5_8*box
1647  end if
1648  limit = 0.5_8*(self%size(i)+task%size(i)+self%ds(i)+task%ds(i))
1649  if (abs(distance) > limit) then
1650  overlaps = .false.
1651  exit
1652  end if
1653  end do
1654  class default
1655  overlaps = self%task_t%overlaps(task)
1656  end select
1657  !call trace%end ()
1658 END FUNCTION overlaps
1659 
1660 !===============================================================================
1661 !> Signed distance between the centers of two tasks in the periodic box
1662 !> This version is adjusted to account for curvilinear coords. (the position of
1663 !> a curvilinear patch is set to the coordinate *origin) and requires components
1664 !> defined in `patch_t`.
1665 !===============================================================================
1666 FUNCTION distance_curvilinear (self, task) RESULT (out)
1667  class(patch_t):: self
1668  class(task_t):: task
1669  real(8):: out(3)
1670  !.............................................................................
1671  call trace%begin ('patch_t%distance_curvilinear',2)
1672  !write (io_unit%log,*) 'distance'; flush(io_unit%log)
1673  select type (task)
1674  class is (patch_t)
1675  out = self%centre_nat-task%centre_nat
1676  !write (io_unit%log,*) out; flush(io_unit%log)
1677  class default
1678  out = self%position-task%position
1679  end select
1680  ! account for periodic wrapping
1681  !write (io_unit%log,*) self%box; flush(io_unit%log)
1682  where (self%periodic .and. self%box > 0d0)
1683  out = modulo(out+0.5_8*self%box, self%box) - 0.5_8*self%box
1684  end where
1685  if (io%verbose>1 .and. (task%id==io%id_debug .or. self%id==io%id_debug)) then
1686  write(io%output,1) &
1687  'distance: self', self%id, self%position, self%size, self%box
1688  write(io%output,1) &
1689  'distance: task', task%id, task%position, task%size, task%box, out
1690  1 format(a,i5,4(3x,3f10.2))
1691  end if
1692  call trace%end()
1693 END FUNCTION distance_curvilinear
1694 
1695 !===============================================================================
1696 !> The nearest smaller index position in the self patch coordinate system. The
1697 !> point with index self%mesh%lb is half a cell inside the patch guard size, so the
1698 !> floating point position is counted relative to that point, and is split into
1699 !> an integer part and a remainder. For patches smaller than the box, the index
1700 !> may become both smaller than self%mesh%lb and larger than self%mesh%ub. For the root
1701 !> grid patch, the index is in the range (li-1,ui) inside the interior.
1702 !===============================================================================
1703 FUNCTION index (self, other, oindex, f)
1704  class(patch_t):: self, other
1705  real, dimension(3):: position, f, pos, pos_lb
1706  integer, dimension(3):: index, oindex
1707  !.............................................................................
1708  pos = other%position + (oindex-other%offset)*other%ds ! position in other
1709  pos = pos-self%position ! relative
1710  pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box ! wrap
1711  f = self%offset + pos/self%ds ! floating integer
1712  index = f ! smaller integer
1713  f = f - index ! remainder
1714 END FUNCTION index
1715 
1716 !===============================================================================
1717 !> The nearest index position in the self patch of the point in another patch.
1718 !> The cost of this operation is 4 plus and 2 multiply per direction = 18 flops
1719 !===============================================================================
1720 FUNCTION nearest (self, other, index) RESULT (out)
1721  class(patch_t):: self, other
1722  integer, dimension(3):: out, index
1723  real, dimension(3):: position, pos_lb, f, pos
1724  !.............................................................................
1725  pos = other%position + (index-other%offset)*other%ds ! position in other
1726  pos = pos-self%position ! relative
1727  pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box ! wrap
1728  f = self%offset + pos/self%ds ! floating integer
1729  out = nint(f) ! nearest integer
1730 END FUNCTION nearest
1731 
1732 !===============================================================================
1733 !> The lowest index position in the self patch of the point in another patch
1734 !===============================================================================
1735 FUNCTION lowest (self, other, index) RESULT (out)
1736  class(patch_t):: self, other
1737  integer, dimension(3):: out, index
1738  real, dimension(3):: position, pos_lb, f, pos
1739  !.............................................................................
1740  pos = other%position + (index-other%offset-0.49)*other%ds ! position in other
1741  pos = pos-self%position ! relative
1742  pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box ! wrap
1743  f = self%offset + pos/self%ds ! floating integer
1744  out = nint(f) ! lowest integer
1745 END FUNCTION lowest
1746 !===============================================================================
1747 !> The highest index position in the self patch of the point in another patch
1748 !===============================================================================
1749 FUNCTION highest (self, other, index) RESULT (out)
1750  class(patch_t):: self, other
1751  integer, dimension(3):: out, index
1752  real, dimension(3):: position, pos_lb, f, pos
1753  !.............................................................................
1754  pos = other%position + (index-other%offset+0.49)*other%ds ! position in other
1755  pos = pos-self%position ! relative
1756  pos = modulo(pos+0.5*self%box,self%box) - 0.5*self%box ! wrap
1757  f = self%offset + pos/self%ds ! floating integer
1758  out = nint(f) ! highest integer
1759 END FUNCTION highest
1760 
1761 !===============================================================================
1762 !> Courant condition, one way or another
1763 !===============================================================================
1764 SUBROUTINE courant_condition (self, detailed_timer)
1765  class(patch_t):: self
1766  logical, optional:: detailed_timer
1767  real(8):: dtime, ds(3), courant
1768  integer:: i2
1769  integer, save:: itimer=0
1770  logical, parameter:: save_dbg=.false.
1771  !.............................................................................
1772  call trace%begin('patch_t%courant_condition', itimer=itimer)
1773  if (io%verbose>1) &
1774  write (io_unit%log,*) self%id, ' courant: u_max, ds', self%u_max, minval(self%ds)
1775  ds(1) = self%ds(1)
1776  ds(2) = self%ds(2)*self%mesh(1)%h2c(self%mesh(1)%li)
1777  ds(3) = self%ds(3)*self%mesh(1)%h31c(self%mesh(1)%li)*abs(self%mesh(2)%h32c(self%mesh(2)%li))
1778  if (trim(self%mesh(3)%type) == 'spherical_mesh') then
1779  do i2=self%mesh(2)%li,self%mesh(2)%ui
1780  ds(3) = min(ds(3),self%ds(3)*self%mesh(1)%h31c(i2)*abs(self%mesh(2)%h32c(i2)))
1781  end do
1782  end if
1783  if (self%dt_fixed > 0.0) then
1784  dtime = self%dt_fixed
1785  courant = dtime*self%u_max/minval(ds)
1786  else if (self%u_fixed > 0.0) then
1787  dtime = self%courant*minval(ds)/self%u_fixed
1788  courant = self%courant
1789  else
1790  dtime = self%courant*minval(ds)/self%u_max*self%safety_factor
1791  courant = self%courant
1792  endif
1793  !$omp atomic write
1794  self%dtime = dtime
1795  if (self%nt > 1 .and. self%istep > 5) then
1796  if (dtime < 0.5*self%dt(self%iit(self%nt-1)) .and. self%dt_fixed < 0.0) then
1797  self%safety_factor = 0.5
1798  end if
1799  if (self%safety_factor < 0.6 .and. save_dbg) then
1800  !$omp critical (fort88_cr)
1801  if (self%safety_factor == 0.5) then
1802  write (io_unit%dbg) self%id, self%gn, self%ipos
1803  write (io_unit%dbg) self%t(self%nt-3), self%mem(:,:,:,:,self%iit(self%nt-3),1)
1804  write (io_unit%dbg) self%id, self%gn, self%ipos
1805  write (io_unit%dbg) self%t(self%nt-2), self%mem(:,:,:,:,self%iit(self%nt-2),1)
1806  end if
1807  write (io_unit%dbg) self%id, self%gn, self%ipos
1808  write (io_unit%dbg) self%time, self%mem(:,:,:,:,self%it,1)
1809  flush (io_unit%dbg)
1810  if (io%verbose>2 .or. (self%istep>5)) &
1811  write(io%output,'(a,i6,i3,f12.6,1p,3g11.2,0p,1x,3f12.6,1x,1p,7g12.3)') &
1812  ' courant: ', self%id, self%it, &
1813  self%time, dtime, courant, self%u_max, self%position, self%dt(self%iit)
1814  !$omp end critical (fort88_cr)
1815  end if
1816  self%safety_factor = self%safety_factor**0.9
1817  self%get_umax_location = (self%safety_factor < 0.9)
1818  end if
1819  !-----------------------------------------------------------------------------
1820  ! Force the time step to give exactly sync_time if within reach
1821  !-----------------------------------------------------------------------------
1822  if (self%time < self%sync_time .and. self%time+self%dtime > self%sync_time) then
1823  self%dtime = self%sync_time - self%time
1824  self%syncing = .true.
1825  if (io%verbose>2) &
1826  write (io_unit%mpi,*) 'patch_t%courant_condition: syncing turned on for', &
1827  self%id, self%time
1828  end if
1829  !-----------------------------------------------------------------------------
1830  ! Prevent large (cheap) patches from taking time steps longer than the out_time
1831  !-----------------------------------------------------------------------------
1832  if (io%do_output .and. io%out_time > 0.0) &
1833  self%dtime = min(self%dtime, io%out_time)
1834  self%dt(self%it) = self%dtime
1835  if (io%verbose>3) &
1836  write(io%output,'(a,i7,f12.6,f8.4,1p,e10.2)') 'courant:', &
1837  self%id, self%dtime, self%courant, self%u_max
1838  call trace%end (itimer)
1839 END SUBROUTINE courant_condition
1840 
1841 !===============================================================================
1842 FUNCTION fsum4 (self, f) RESULT (out)
1843  class(patch_t):: self
1844  real:: out
1845  real, dimension(:,:,:), intent(in):: f
1846  integer:: n(3)
1847  !-----------------------------------------------------------------------------
1848  if (self%is_periodic()) then
1849  n = self%ncell-1
1850  out = sum(real(f(self%mesh(1)%li:self%mesh(1)%li+n(1), & self%mesh(2)%li:self%mesh(2)%li+n(2), & self%mesh(3)%li:self%mesh(3)%li+n(3)),kind=8))
1851  else
1852  out = sum(real(f(self%mesh(1)%li:self%mesh(1)%ui, & self%mesh(2)%li:self%mesh(2)%ui, & self%mesh(3)%li:self%mesh(3)%ui),kind=8))
1853  end if
1854 END FUNCTION fsum4
1855 !===============================================================================
1856 FUNCTION fsum8 (self, f) RESULT (out)
1857  class(patch_t):: self
1858  real(8):: out
1859  real(8), dimension(:,:,:), intent(in):: f
1860  integer:: n(3)
1861  !-----------------------------------------------------------------------------
1862  if (self%is_periodic()) then
1863  n = self%ncell-1
1864  out = sum(real(f(self%mesh(1)%li:self%mesh(1)%li+n(1), & self%mesh(2)%li:self%mesh(2)%li+n(2), & self%mesh(3)%li:self%mesh(3)%li+n(3)),kind=8))
1865  else
1866  out = sum(real(f(self%mesh(1)%li:self%mesh(1)%ui, & self%mesh(2)%li:self%mesh(2)%ui, & self%mesh(3)%li:self%mesh(3)%ui),kind=8))
1867  end if
1868 END FUNCTION fsum8
1869 
1870 !===============================================================================
1871 FUNCTION faver4 (self, f) RESULT (out)
1872  class(patch_t):: self
1873  real:: out
1874  real, dimension(:,:,:), intent(in):: f
1875  !-----------------------------------------------------------------------------
1876  out = self%fsum4(f)/product(self%ncell)
1877 END FUNCTION faver4
1878 !===============================================================================
1879 FUNCTION faver8 (self, f) RESULT (out)
1880  class(patch_t):: self
1881  real(8):: out
1882  real(8), dimension(:,:,:), intent(in):: f
1883  !-----------------------------------------------------------------------------
1884  out = self%fsum8(f)/product(self%ncell)
1885 END FUNCTION faver8
1886 
1887 !===============================================================================
1888 FUNCTION fminvali (self, iv) RESULT (fminval)
1889  class(patch_t):: self
1890  integer:: iv
1891  real:: fminval
1892  !-----------------------------------------------------------------------------
1893  fminval = self%fminvalvar (self%mem(:,:,:,iv,self%it,1))
1894 END FUNCTION fminvali
1895 !===============================================================================
1896 FUNCTION fmaxvali (self, iv) RESULT (fmaxval)
1897  class(patch_t):: self
1898  integer:: iv
1899  real:: fmaxval
1900  !-----------------------------------------------------------------------------
1901  fmaxval = self%fmaxvalvar (self%mem(:,:,:,iv,self%it,1))
1902 END FUNCTION fmaxvali
1903 !===============================================================================
1904 FUNCTION fmaxval4 (self, f, outer) RESULT (fmaxval)
1905  real:: fmaxval
1906  class(patch_t):: self
1907  real, dimension(:,:,:), intent(in):: f
1908  logical, optional:: outer
1909  logical:: outer_l
1910  !-----------------------------------------------------------------------------
1911  if (present(outer)) then
1912  outer_l = outer
1913  else
1914  outer_l = .false.
1915  end if
1916  if (self%boundaries%is_set (bits%spherical)) then
1917  fmaxval = maxval(f, mask=self%boundaries%mask(self%mesh))
1918  else if (outer_l) then
1919  fmaxval = maxval(f(self%mesh(1)%lo:self%mesh(1)%uo, &
1920  self%mesh(2)%lo:self%mesh(2)%uo, &
1921  self%mesh(3)%lo:self%mesh(3)%uo))
1922  else
1923  fmaxval = maxval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1924  self%mesh(2)%li:self%mesh(2)%ui, &
1925  self%mesh(3)%li:self%mesh(3)%ui))
1926  end if
1927 END FUNCTION fmaxval4
1928 !===============================================================================
1929 FUNCTION fminval4 (self, f) RESULT(fminval)
1930  real:: fminval
1931  class(patch_t):: self
1932  real, dimension(:,:,:), intent(in):: f
1933  !-----------------------------------------------------------------------------
1934  if (self%boundaries%is_set (bits%spherical)) then
1935  fminval = minval(f, mask=self%boundaries%mask(self%mesh))
1936  else
1937  fminval = minval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1938  self%mesh(2)%li:self%mesh(2)%ui, &
1939  self%mesh(3)%li:self%mesh(3)%ui))
1940  end if
1941 END FUNCTION fminval4
1942 !===============================================================================
1943 FUNCTION fmaxval8 (self, f, outer) RESULT (fmaxval)
1944  real(8):: fmaxval
1945  class(patch_t):: self
1946  real(8), dimension(:,:,:), intent(in):: f
1947  logical, optional:: outer
1948  logical:: outer_l
1949  !-----------------------------------------------------------------------------
1950  if (present(outer)) then
1951  outer_l = outer
1952  else
1953  outer_l = .false.
1954  end if
1955  if (self%boundaries%is_set (bits%spherical)) then
1956  fmaxval = maxval(f, mask=self%boundaries%mask(self%mesh))
1957  else if (outer_l) then
1958  fmaxval = maxval(f(self%mesh(1)%lo:self%mesh(1)%uo, &
1959  self%mesh(2)%lo:self%mesh(2)%uo, &
1960  self%mesh(3)%lo:self%mesh(3)%uo))
1961  else
1962  fmaxval = maxval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1963  self%mesh(2)%li:self%mesh(2)%ui, &
1964  self%mesh(3)%li:self%mesh(3)%ui))
1965  end if
1966 END FUNCTION fmaxval8
1967 !===============================================================================
1968 FUNCTION fminval8 (self, f) RESULT(fminval)
1969  real(8):: fminval
1970  class(patch_t):: self
1971  real(8), dimension(:,:,:), intent(in):: f
1972  !-----------------------------------------------------------------------------
1973  if (self%boundaries%is_set (bits%spherical)) then
1974  fminval = minval(f, mask=self%boundaries%mask(self%mesh))
1975  else
1976  fminval = minval(f(self%mesh(1)%li:self%mesh(1)%ui, &
1977  self%mesh(2)%li:self%mesh(2)%ui, &
1978  self%mesh(3)%li:self%mesh(3)%ui))
1979  end if
1980 END FUNCTION fminval8
1981 
1982 !===============================================================================
1983 !> Check for NaN
1984 !===============================================================================
1985 SUBROUTINE check_nan (self, panic, label, always)
1986  class(patch_t):: self
1987  logical:: panic
1988  character(len=*), optional:: label
1989  logical, optional:: always
1990  !.............................................................................
1991  integer:: ix, iy, iz, iv, ii(3)
1992  logical:: nan
1993  integer, save:: itimer=0
1994  !-----------------------------------------------------------------------------
1995  if (.not. do_check_nan) then
1996  if (io%verbose < 2 .and. .not. present(always)) return
1997  end if
1998  if (.not. associated(self%mem)) return
1999  call timer%begin ('patch_t%check_nan', itimer=itimer)
2000  if (io%verbose>1) then
2001  if (present(label)) then
2002  write(io%output,*) self%id, trim(label)//', min(d), max(d):', &
2003  minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2004  maxval(self%mem(:,:,:,self%idx%d,self%it,1))
2005  else
2006  write(io%output,*) self%id, 'min(d), max(d):', &
2007  minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2008  maxval(self%mem(:,:,:,self%idx%d,self%it,1))
2009  end if
2010  end if
2011  nan = .false.
2012  do iv=1,self%nv
2013  do iz=1,self%gn(3)
2014  do iy=1,self%gn(2)
2015  do ix=1,self%gn(1)
2016  if (isinf(self%mem(ix,iy,iz,iv,self%it,1))) then
2017  if (io%verbose>0) then
2018  if (present(label)) then
2019  write(io%output,1) trim(label)//' ERROR: NaN id,ix,iy,iz,iv:', &
2020  self%id,ix,iy,iz,iv
2021  1 format(a,i6,3i4,i5)
2022  else
2023  write(io%output,1) 'ERROR: NaN ix,iy,iz,iv:', ix,iy,iz,iv
2024  end if
2025  flush (io%output)
2026  end if
2027  ii = [ix,iy,iz]
2028  nan = .true.
2029  end if
2030  end do
2031  end do
2032  end do
2033  end do
2034  if (nan) then
2035  flush (io%output)
2036  panic = .true.
2037  if (present(label)) then
2038  write(io%output,*) self%id, trim(label)//' ERROR NaN, min(d), max(d):', &
2039  minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2040  maxval(self%mem(:,:,:,self%idx%d,self%it,1)), ii
2041  else
2042  write(io%output,*) self%id, ' ERROR NaN, min(d), max(d):', &
2043  minval(self%mem(:,:,:,self%idx%d,self%it,1)), &
2044  maxval(self%mem(:,:,:,self%idx%d,self%it,1)), ii
2045  end if
2046  write(io%output,*) self%unsigned
2047  flush (io%output)
2048  call io%abort ('found NaN -- check rank & thread logs!')
2049  end if
2050  call timer%end (itimer)
2051 END SUBROUTINE check_nan
2052 
2053 LOGICAL FUNCTION isinf(a)
2054  use ieee_arithmetic, only : ieee_is_nan
2055  real(kind=KindScalarVar):: a
2056  isinf = ieee_is_nan(a) .or. (abs(a) > huge(1.0_kindscalarvar))
2057 END FUNCTION isinf
2058 
2059 !===============================================================================
2060 !> Check if the density is negative somewhere, counting points in the interior
2061 !> and guard zones separately. Attempt to patch the guard zones, if applicable.
2062 !===============================================================================
2063 SUBROUTINE check_density (self, label)
2064  class(patch_t):: self
2065  character(len=*), optional:: label
2066  real(kind=KindScalarVar):: mind, posd
2067  real:: tiny = 1.0e-32
2068  integer:: ix, iy, iz, iv, ii(3), ntot, nins, l(3), u(3), it, jt, iw
2069  logical:: inside, nan, panic
2070  integer, save:: itimer=0
2071  !-----------------------------------------------------------------------------
2072  if ((io%verbose < 2 .or. io%print_time > 0.0 .or. print_every > 1) .and. &
2073  .not.do_check_nan) return
2074  if (self%is_set(bits%no_density+bits%frozen)) return
2075  if (.not. associated(self%mem)) return
2076  call timer%begin ('patch_t%check_density', itimer=itimer)
2077  l = self%mesh%lb
2078  u = self%mesh%ub
2079  mind = minval(self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%d,self%it,1))
2080  if (present(label) .and. io%verbose>1) write (io_unit%log,*) label, self%id, mind
2081  panic = .false.
2082  if (mind <= 0.0) then
2083  panic = .true.
2084  !---------------------------------------------------------------------------
2085  ! Negative densities appear to develop at a few points in the guard zones,
2086  ! so this tries to patch up, sampling densitites in the guard zone.
2087  !---------------------------------------------------------------------------
2088  posd = tiny
2089  ntot = 0
2090  nins = 0
2091  do iz=1,self%gn(3)
2092  do iy=1,self%gn(2)
2093  do ix=1,self%gn(1)
2094  ii = [ix,iy,iz]
2095  inside = all(ii>self%mesh%lo .and. ii<self%mesh%uo)
2096  posd = merge(posd, max(posd,self%mem(ix,iy,iz,self%idx%d,self%it,1)), inside)
2097  ntot = ntot + merge(1,0,self%mem(ix,iy,iz,self%idx%d,self%it,1)<=0.0)
2098  nins = nins + merge(1,0,self%mem(ix,iy,iz,self%idx%d,self%it,1)<=0.0 .and. inside)
2099  !-------------------------------------------------------------------------
2100  ! Pick the sampled positive value if a negative one is found in the
2101  ! guard zone
2102  !-------------------------------------------------------------------------
2103  !self%mem(ix,iy,iz,1,self%it,1) = merge (self%mem(ix,iy,iz,1,self%it,1), &
2104  ! posd, self%mem(ix,iy,iz,1,self%it,1) > 0.0 .or. inside)
2105  end do
2106  end do
2107  end do
2108  1 format(a,i8,3x,a,i8,a,i8,a)
2109  if (present(label)) then
2110  if (nins > 0) then
2111  write (io_unit%output,1) 'ID ',self%id, &
2112  'ERROR: density non-positive in',ntot,' points,',nins,' inside, '// &
2113  trim(label)
2114  else if (io%verbose>0) then
2115  write (io_unit%output,1) 'ID ',self%id, &
2116  'WARNING: density non-positive in',ntot,' points,',nins,' inside, '// &
2117  trim(label)
2118  end if
2119  else
2120  if (nins > 0) then
2121  write (io_unit%output,1) 'ID ',self%id, &
2122  'WARNING: density non-positive in',ntot,' points,',nins,' inside'
2123  else if (io%verbose>0) then
2124  write (io_unit%output,1) 'ID ',self%id, &
2125  'WARNING: density non-positive in',ntot,' points,',nins,' inside'
2126  end if
2127  end if
2128  if (nins > 0) &
2129  panic = .true.
2130  end if
2131  call self%check_nan (panic, label)
2132  !-----------------------------------------------------------------------------
2133  ! If error indicated by check_density or check_nan, dump the patch and abort
2134  !-----------------------------------------------------------------------------
2135  if (panic) then
2136  !$omp critical (panic_cr)
2137  write (io_unit%dump) self%id, self%gn, self%nv, self%nt, self%nw, self%istep
2138  write (io_unit%dump) self%time, self%dtime
2139  do iw=1,self%nw
2140  do it=1,self%nt
2141  !$omp atomic read
2142  jt = self%iit(it)
2143  !$omp end atomic
2144  write (io_unit%dump) self%mem(:,:,:,:,jt,iw)
2145  end do
2146  end do
2147  !$omp end critical (panic_cr)
2148  if (nins > 0 .or. io%verbose > 0) &
2149  print *, mpi%rank, self%id, &
2150  'WARNING: found non-positive density points -- check rank and thread logs'
2151  if (nins > 0) then
2152  write (stderr,*) 'check dump, rank and thread logs for rank', mpi%rank
2153  flush (stdout)
2154  flush (stderr)
2155  flush (io_unit%log)
2156  flush (io_unit%dump)
2157  call io%abort ('found non-positive internal density points')
2158  end if
2159  end if
2160  call timer%end (itimer)
2161 END SUBROUTINE check_density
2162 
2163 !===============================================================================
2164 !> Statistics for all vars in a patch
2165 !===============================================================================
2166 SUBROUTINE stats_patch (self, label)
2167  class(patch_t):: self
2168  character(len=*):: label
2169  integer:: iv
2170  real(8), dimension(4):: si, sb
2171  do iv=1,self%nv
2172  call scalar_stats (self%mem(:,:,:,iv,self%it,1), self%mesh%li, self%mesh%ui, si)
2173  call scalar_stats (self%mem(:,:,:,iv,self%it,1), self%mesh%lb, self%mesh%ub, sb)
2174  write (io_unit%log,'(i7,2x,a,i3,1p,2(3x,4g12.4))') &
2175  self%id, trim(label)//' stats: iv, min, aver, rms, max =', iv, si, sb
2176  end do
2177 END SUBROUTINE stats_patch
2178 
2179 !===============================================================================
2180 !> Statistics for a scalar field
2181 !===============================================================================
2182 SUBROUTINE stats_scalar (self, var, label)
2183  class(patch_t):: self
2184  real(kind=KindScalarVar), dimension(:,:,:):: var
2185  character(len=*):: label
2186  !.............................................................................
2187  real(8), dimension(4):: si, sb
2188  !-----------------------------------------------------------------------------
2189  call scalar_stats (var, self%mesh%li, self%mesh%ui, si)
2190  call scalar_stats (var, self%mesh%lb, self%mesh%ub, sb)
2191  write (io_unit%output,'(i7,2x,a,1p,2(3x,4g12.4))') &
2192  self%id, trim(label)//' stats: min, aver, rms, max =', si, sb
2193 END SUBROUTINE stats_scalar
2194 
2195 !===============================================================================
2196 !> Statistics for a vector field
2197 !===============================================================================
2198 SUBROUTINE stats_vector (self, var, label)
2199  class(patch_t):: self
2200  real(kind=KindScalarVar), dimension(:,:,:,:):: var
2201  character(len=*):: label
2202  !.............................................................................
2203  integer:: iv
2204  real(8), dimension(4):: si, sb
2205  !-----------------------------------------------------------------------------
2206  do iv=1,size(var,4)
2207  call scalar_stats (var(:,:,:,iv), self%mesh%li, self%mesh%ui, si)
2208  call scalar_stats (var(:,:,:,iv), self%mesh%lb, self%mesh%ub, sb)
2209  write (io_unit%output,'(i7,2x,a,i3,1p,2(3x,4g12.4))') &
2210  self%id, trim(label)//' stats: iv, min, aver, rms, max =', iv, si, sb
2211  end do
2212 END SUBROUTINE stats_vector
2213 
2214 !===============================================================================
2215 !===============================================================================
2216 SUBROUTINE scalar_stats (f, l, u, s)
2217  real(kind=KindScalarVar), dimension(:,:,:):: f
2218  integer:: l(3), u(3)
2219  real(8), dimension(4):: s
2220  integer:: ix, iy, iz
2221  s(2) = 0d0
2222  s(3) = 0d0
2223  s(1) = f(1,1,1)
2224  s(4) = f(1,1,1)
2225  do iz=l(3),u(3)
2226  do iy=l(2),u(2)
2227  do ix=l(1),u(1)
2228  s(2) = s(2) + f(ix,iy,iz)
2229  s(3) = s(3) + f(ix,iy,iz)**2
2230  s(1) = min(s(1),f(ix,iy,iz))
2231  s(4) = max(s(4),f(ix,iy,iz))
2232  end do
2233  end do
2234  end do
2235  s(2) = s(2)/product(u-l+1)
2236  s(3) = sqrt(max(0d0,s(3)/product(u-l+1)-s(2)**2))
2237 END SUBROUTINE scalar_stats
2238 
2239 !===============================================================================
2240 !> default version; only supports Cartesian coords.
2241 !> see the `solvers/zeus3d` version for curvilinear coord. suppport.
2242 !===============================================================================
2243 SUBROUTINE mapvectors (self, in, target, ii, jj, out)
2244  class(patch_t) :: self
2245  class(patch_t), pointer :: target
2246  real(kind=KindScalarVar), dimension(:), intent(in):: in
2247  real(kind=KindScalarVar), intent(out):: out(size(in))
2248  integer, intent(in) :: ii(3), jj(3)
2249  !.............................................................................
2250  call trace%begin('patch_t%MapVectors', 2)
2251 
2252  select case (self%mesh_type)
2253  case (1) ! Cartesian
2254  if (target%mesh_type == mesh_types%Cartesian) then
2255  ! no conversion necessary
2256  out = in
2257  else
2258  call mpi%abort('Not implemented error! patch_t%MapVectors')
2259  end if
2260  case default ! other; not implemented
2261  call mpi%abort('Not implemented error! patch_t%MapVectors')
2262  end select
2263 
2264  call trace%end()
2265 END SUBROUTINE mapvectors
2266 
2267 !===============================================================================
2268 FUNCTION beyond_patch_edge (self, idir, position) result (outside)
2269  class(patch_t) :: self
2270  integer, intent(in) :: idir
2271  real(8), intent(in) :: position(3)
2272  logical :: outside
2273  real(8) :: edge, curvpos(3)
2274  !.............................................................................
2275 
2276  if (self%n(idir) == 1) then
2277  outside = .false.
2278  return
2279  endif
2280 
2281  select case (self%mesh_type)
2282  case (1) ! Cartesian
2283  edge = 0.5 * self%gsize(idir)
2284  outside = abs(position(idir)-self%position(idir)) > edge
2285  case (2) ! spherical
2286  curvpos = position - self%centre_cart
2287  curvpos = cartesiantospherical(curvpos(1), curvpos(2), curvpos(3))
2288  select case (idir)
2289  case (1)
2290  edge = self%gsize(idir) + self%llc_nat(idir)
2291  outside = abs(curvpos(idir)) > edge .or. abs(curvpos(idir)) < self%llc_nat(idir)
2292  case (2)
2293  if (associated(self%mesh)) then
2294  edge = self%gsize(idir) + self%mesh(2)%xi
2295  else
2296  edge = self%gsize(idir)
2297  end if
2298  outside = abs(curvpos(idir)) > edge
2299  case (3)
2300  if (associated(self%mesh)) then
2301  edge = self%gsize(idir) + self%mesh(3)%zeta
2302  else
2303  edge = self%gsize(idir)
2304  end if
2305  outside = abs(curvpos(idir)) > edge
2306  end select
2307  outside = abs(curvpos(idir)) > edge
2308  case (3) ! cylindrical
2309  curvpos = position - self%centre_cart
2310  curvpos = cartesiantocylindrical(curvpos(1), curvpos(2), curvpos(3))
2311  select case (idir)
2312  case (1)
2313  edge = 0.5 * self%gsize(idir)
2314  outside = abs(curvpos(idir)) > edge
2315  case (2)
2316  edge = self%gsize(idir) + self%llc_nat(idir)
2317  outside = abs(curvpos(idir)) > edge .or. abs(curvpos(idir)) < self%llc_nat(idir)
2318  case (3)
2319  if (associated(self%mesh)) then
2320  edge = self%gsize(idir) + self%mesh(3)%zeta
2321  else
2322  edge = self%gsize(idir)
2323  end if
2324  outside = abs(curvpos(idir)) > edge
2325  end select
2326  end select
2327 END FUNCTION beyond_patch_edge
2328 
2329 !=======================================================================
2330 !> Interpolate -- only in time for now -- from the MHD time slots save in
2331 !> this module to the time requested from the solve procedure
2332 !=======================================================================
2333 SUBROUTINE log_interpolate (self, time, i, a)
2334  class(patch_t) :: self
2335  real(8) :: time
2336  integer :: i
2337  real(kind=KindScalarVar), dimension(:,:,:), pointer :: a, b
2338  !.....................................................................
2339  integer :: jt(2)
2340  real :: pt(2)
2341  integer, save :: itimer=0
2342  ! --------------------------------------------------------------------
2343  call trace%begin ('patch_t%log_interpolate', itimer=itimer)
2344  call self%time_interval (time, jt, pt)
2345  if (io%verbose > 1) then
2346  b => self%mem(:,:,:,i,jt(1),1)
2347  write(io%output,*) self%id, 'log_interp(1)', jt(1), minval(b), maxval(b)
2348  b => self%mem(:,:,:,i,jt(2),1)
2349  write(io%output,*) self%id, 'log_interp(2)', jt(2), minval(b), maxval(b)
2350  end if
2351 
2352  a = log(self%mem(:,:,:,i,jt(1),1))*pt(1) &
2353  + log(self%mem(:,:,:,i,jt(2),1))*pt(2)
2354  call trace%end (itimer)
2355 END SUBROUTINE log_interpolate
2356 
2357 !=======================================================================
2358 !> Search for a time interval to interpolate in, and return the time
2359 !> slot indices and weight factors
2360 !=======================================================================
2361 SUBROUTINE time_interval (self, time, jt, pt, all)
2362  class(patch_t) :: self
2363  real(8) :: time
2364  integer :: jt(2)
2365  real :: pt(2)
2366  logical, optional :: all
2367  !.....................................................................
2368  integer :: i, j, n, iit(self%nt-1)
2369  real(8) :: t(self%nt-1), dt
2370  integer, save :: itimer=0
2371  ! --------------------------------------------------------------------
2372  !call trace%begin ('patch_t%time_interval', 2, itimer=itimer)
2373  if (self%nt==1) then
2374  jt = 1
2375  pt(2) = 0.0
2376  else
2377  call self%timeslots (iit, t)
2378  j = 1
2379  n = self%nt-2
2380  !-------------------------------------------------------------------
2381  ! Seach and set the j index to the last point that has t(j)<time.
2382  ! Note that this does NOT guarantee that t(j+1) >= time. In cases
2383  ! where both t(nt-2) < time and t(n-1) < time, extrapolation in time
2384  ! will result. What should NOT happen is that a memory slot that
2385  ! has not yet received an update is chosen.
2386  !-------------------------------------------------------------------
2387  do i=1,n
2388  if (t(i) < time) j=i
2389  end do
2390  jt = iit(j:j+1)
2391  dt = t(j+1)-t(j)
2392  if (dt==0.0) then
2393  pt(2) = 0.0
2394  else
2395  pt(2) = (time-t(j))/dt
2396  end if
2397  end if
2398  pt(1) = 1.-pt(2)
2399  if ((io%verbose>3 .or. io_unit%do_validate) .and. pt(2)>0.0) then
2400  write(io%output,'(a,i6,i6,1p,e14.6,i3,0p,2f8.4,2x,1p,7e14.6)') &
2401  'time_interval: times =', self%id, self%istep, time, j, pt, t
2402  end if
2403  !call trace%end (itimer)
2404 END SUBROUTINE time_interval
2405 
2406 !===============================================================================
2407 !===============================================================================
2408 FUNCTION index_stagger (self, iv) RESULT (out)
2409  class(patch_t) :: self
2410  integer :: iv
2411  integer :: i
2412  real(8) :: out(3)
2413  !-----------------------------------------------------------------------------
2414  do i=1,3
2415  out(i) = self%mesh(i)%h(iv)
2416  end do
2417 END FUNCTION index_stagger
2418 
2419 !===============================================================================
2420 !> Return indices and fractions for 3D interpolation. Only called from
2421 !> patch_mod::interpolation, where a possible issue is to ensure the correct
2422 !> behavior when the source and target points are essentially the same, except
2423 !> for roundoff. We then want to make sure that the index becomes the same,
2424 !> with fraction vanishing, so there is no possible issue with extrapolations.
2425 !> Hence we should add a small fraction before taking the floor value, and then
2426 !> subtract the same.
2427 !===============================================================================
2428 SUBROUTINE index_space (self, pos, iv, i, p)
2429  class(patch_t):: self
2430  real(8) :: pos(3)
2431  integer :: iv
2432  integer :: i(3)
2433  real :: p(3)
2434  !.............................................................................
2435  real :: h(3)
2436  real, parameter :: eps=0.001
2437  integer, save :: itimer=0
2438  !-----------------------------------------------------------------------------
2439  !call trace%begin ('patch_t%index_space', itimer=itimer)
2440  h = self%index_stagger(iv)
2441  p = pos - self%position ! relative position
2442  p = modulo(p+0.5_8*self%mesh%b,self%mesh%b) & ! wrap
2443  -0.5_8*self%mesh%b
2444  p = p/self%mesh%d + self%mesh%o - h ! float index
2445  p = p + eps ! ensure same points ok
2446  i = p ! integer part
2447  i = max(self%mesh%li, &
2448  min(self%mesh%ui-1,i)) ! valid range
2449  p = p - eps - i ! remaining fraction
2450  !-----------------------------------------------------------------------------
2451  ! The slight allowance for extrapolation below is to detect it in tests
2452  !-----------------------------------------------------------------------------
2453  if (zero_order_extrap .and. .not.io_unit%do_validate) &
2454  p = min(p,max(p,-0.01),1.01) ! 0th order extrapolation
2455  !call trace%end (itimer)
2456 END SUBROUTINE index_space
2457 
2458 !===============================================================================
2459 !> Optimized case, where it is known that target and source (self) overlap, so
2460 !> no periodic mapping is needed, and iv is always present
2461 !===============================================================================
2462 SUBROUTINE index_space_of (self, target, ii, iv, jj, pp)
2463  class(patch_t) :: self
2464  class(patch_t), pointer :: target
2465  integer :: ii(3)
2466  integer :: iv
2467  real :: pp(3)
2468  integer :: jj(3)
2469  !.............................................................................
2470  class(mesh_t), pointer :: mt, ms
2471  integer :: i, j
2472  real :: p
2473  real, parameter :: eps=0.001
2474  integer, save :: itimer=0
2475  !-----------------------------------------------------------------------------
2476  do i=1,3
2477  mt => target%mesh(i) ! target mesh
2478  ms => self%mesh(i) ! source mesh
2479  p = mt%p - ms%p + (ii(i) - mt%o + mt%h(iv))*mt%d ! real position
2480  p = p/ms%d + ms%o - ms%h(iv) + eps ! float index
2481  j = p ! integer part
2482  j = max(ms%li, &
2483  min(ms%ui-1,j)) ! valid range
2484  p = p - eps - j ! remaining fraction
2485  if (zero_order_extrap .and. .not.io_unit%do_validate) &
2486  p = min(max(p,-0.01),1.01) ! 0th order extrapolation
2487  pp(i) = p
2488  jj(i) = j
2489  end do
2490 END SUBROUTINE index_space_of
2491 
2492 !===============================================================================
2493 !> Return only the index in self corresponding to a point in native space.
2494 !> "pos" is a relative position in target coordinate space, and since we are
2495 !> not informed about the target grid spacing it must be adjusted for staggering
2496 !> before the call.
2497 !===============================================================================
2498 FUNCTION index_only (self, pos, iv, roundup, nowrap) RESULT (ii)
2499  class(patch_t) :: self
2500  real(8), intent(in) :: pos(3)
2501  integer, optional :: iv
2502  logical, optional :: roundup, nowrap
2503  integer :: ii(3)
2504  !.............................................................................
2505  integer :: i
2506  real(8) :: p(3), h(3)
2507  class(mesh_t), pointer :: mesh
2508  !-----------------------------------------------------------------------------
2509  if (present(iv)) &
2510  h = self%index_stagger(iv)
2511  do i=1,3
2512  mesh => self%mesh(i)
2513  p(i) = pos(i) - self%position(i) ! relative position
2514  if (self%periodic(i)) then
2515  if (.not.present(nowrap)) &
2516  p(i) = modulo(p(i)+0.5_8*mesh%b, mesh%b)-0.5_8*mesh%b ! periodic wrap
2517  end if
2518  p(i) = p(i)/mesh%d + self%offset(i) ! float index
2519  if (present(iv)) then
2520  p(i) = p(i) + mesh%h(iv) ! staggering
2521  end if
2522  if (present(roundup)) then
2523  ii(i) = ceiling(p(i)) ! round up
2524  else
2525  ii(i) = p(i) + 0.0001 ! floor w margin
2526  end if
2527  end do
2528 END FUNCTION index_only
2529 
2530 !===============================================================================
2531 !> Return indices and fractions for 3D interpolation. Only called from
2532 !> patch_mod::interpolation, where a possible issue is to ensure the correct
2533 !> behavior when the source and target points are essentially the same, except
2534 !> for roundoff. We then want to make sure that the index becomes the same,
2535 !> with fraction vanishing, so there is no possible issue with extrapolations.
2536 !> Hence we should add a small fraction before taking the floor value, and then
2537 !> subtract the same.
2538 !===============================================================================
2539 SUBROUTINE index_space2 (self, nbor, pos, iv, ii, p)
2540  class(patch_t):: self
2541  class(patch_t), pointer :: nbor
2542  real(8) :: pos(3)
2543  integer :: iv
2544  integer :: ii(3)
2545  real :: p(3)
2546  !.....................................................................
2547  class(mesh_t), pointer :: m,mn
2548  integer :: i
2549  real(8) :: pp
2550  real(8), parameter:: eps=0.001_8
2551  integer, save :: itimer=0
2552  real, parameter:: sml = -1.0e-5
2553  !-----------------------------------------------------------------------------
2554  !call trace%begin ('patch_t%index_space', itimer=itimer)
2555  do i=1,3
2556  m => self%mesh(i)
2557  mn => nbor%mesh(i)
2558  pp = mn%p - m%p ! relative position
2559  pp = pp + pos(i) + m%h(iv)*(mn%d-m%d) ! position in the local frame
2560  if (self%periodic(i)) &
2561  pp = modulo(pp+0.5_8*m%b,m%b)-0.5_8*m%b ! wrap
2562  pp = pp/m%d + m%o ! float index
2563  !pp = round(pp,rnd)
2564  !pp = pp + eps ! ensure same points ok
2565  ii(i)= floor(pp+eps) ! integer part
2566  ii(i) = max(m%lb, &
2567  min(m%ub,ii(i))) ! valid range
2568  p(i) = pp - real(ii(i),kind=4)
2569  if ((p(i)<0.0).and.(p(i)>sml)) p(i) = 0.0
2570  !p(i) = round(p(i),rnd) ! remaining fraction
2571  end do
2572  !call trace%end (itimer)
2573 END SUBROUTINE index_space2
2574 
2575 !===============================================================================
2576 !> Return only the index in self corresponding to a point in native space
2577 !===============================================================================
2578 FUNCTION index_only2 (self, nbor, pos, iv, roundup, closest) RESULT (ii)
2579  class(patch_t) :: self
2580  class(patch_t), pointer :: nbor
2581  real(8) :: pos(3)
2582  integer :: iv
2583  integer :: ii(3)
2584  logical, optional :: roundup
2585  logical, optional :: closest
2586  !.............................................................................
2587  class(mesh_t), pointer :: m, mn
2588  integer :: i
2589  real(8), parameter:: eps=0.001_8
2590  real(8) :: p(3)
2591  !-----------------------------------------------------------------------------
2592  do i=1,3
2593  m => self%mesh(i)
2594  mn => nbor%mesh(i)
2595  p(i) = mn%p - m%p ! relative position
2596  p(i) = p(i) + pos(i) + m%h(iv)*(mn%d-m%d) ! position in the local frame
2597  if (self%periodic(i)) &
2598  p(i) = modulo(p(i)+0.5_8*m%b,m%b)-0.5_8*m%b ! wrap
2599  p(i) = p(i)/m%d+m%o ! float index
2600  !p(i) = round(p(i),m%rnd)
2601  if (present(roundup)) then
2602  ii(i) = ceiling(p(i)-eps) ! round up, with epsilon substracted
2603  else if (present(closest)) then
2604  ii(i) = nint(p(i))
2605  else
2606  ii(i) = floor(p(i)+eps) ! integer part with added epsilon
2607  end if
2608  ii(i) = max(m%lb,min(m%ub,ii(i))) ! valid range
2609  end do
2610 END FUNCTION index_only2
2611 
2612 !===============================================================================
2613 !> How many edges do `self` and its neighbor `source` share?
2614 !===============================================================================
2615 FUNCTION count_edges (self, source) RESULT (out)
2616  class(patch_t) :: self
2617  class(patch_t), pointer :: source
2618  integer:: out, i
2619  real(8):: dist(3)
2620  !-----------------------------------------------------------------------------
2621  dist = self%distance (source)
2622  out = 0
2623  do i=1,3
2624  if (abs(dist(i)) > self%ds(i)) out = out+1
2625  end do
2626 END FUNCTION count_edges
2627 
2628 !===============================================================================
2629 !> Interpolate in space and time, for a case where it is safe = the source point
2630 !> is inside the valid domain (inside a distance of half a cell from the source
2631 !> patch boundary for non-staggered points)
2632 !===============================================================================
2633 SUBROUTINE interpolate (self, target, ii, iv, jt, pt)
2634  class(patch_t) :: self
2635  class(patch_t), pointer :: target
2636  integer :: ii(3), iv, jt(2)
2637  real :: pt(2)
2638  !.............................................................................
2639  real(8) :: pos(3)
2640  integer :: j(3)
2641  real :: p(3), q(3)
2642  integer, save :: itimer=0
2643  !-----------------------------------------------------------------------------
2644  if (self%pervolume(iv)) then
2645  call self%interpolate_specific (target, ii, iv, jt, pt)
2646  else
2647  pos = target%myposition (ii, iv)
2648  call self%index_space (pos, iv, j, p)
2649  q = 1.0-p
2650  !if (io%verbose > 1) &
2651  ! write(io_unit%log,'(a,2i6,2x,3i4,2x,3f7.3)') &
2652  ! 'interpolate: target, source, j, p =', target%id, self%id, j, p
2653  !call trace%begin ('patch_t%interpolate', itimer=itimer)
2654  if (self%unsigned(iv)) then
2655  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = exp( &
2656  pt(1)*(q(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1)) + &
2657  p(1)*log(self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1))) + &
2658  p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1)) + &
2659  p(1)*log(self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1)))) + &
2660  p(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1)) + &
2661  p(1)*log(self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1))) + &
2662  p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1)) + &
2663  p(1)*log(self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1))))) + &
2664  pt(2)*(q(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1)) + &
2665  p(1)*log(self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1))) + &
2666  p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1)) + &
2667  p(1)*log(self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1)))) + &
2668  p(3)*(q(2)*(q(1)*log(self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1)) + &
2669  p(1)*log(self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1))) + &
2670  p(2)*(q(1)*log(self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1)) + &
2671  p(1)*log(self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))))
2672  else
2673  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
2674  pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3) ,iv,jt(1),1) + &
2675  ( p(1))*self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(1),1)) + &
2676  ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(1),1) + &
2677  ( p(1))*self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(1),1))) + &
2678  ( p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(1),1) + &
2679  ( p(1))*self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(1),1)) + &
2680  ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(1),1) + &
2681  ( p(1))*self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1)))) + &
2682  pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3) ,iv,jt(2),1) + &
2683  ( p(1))*self%mem(j(1)+1,j(2) ,j(3) ,iv,jt(2),1)) + &
2684  ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3) ,iv,jt(2),1) + &
2685  ( p(1))*self%mem(j(1)+1,j(2)+1,j(3) ,iv,jt(2),1))) + &
2686  ( p(3))*((1.-p(2))*((1.0-p(1))*self%mem(j(1) ,j(2) ,j(3)+1,iv,jt(2),1) + &
2687  ( p(1))*self%mem(j(1)+1,j(2) ,j(3)+1,iv,jt(2),1)) + &
2688  ( p(2))*((1.0-p(1))*self%mem(j(1) ,j(2)+1,j(3)+1,iv,jt(2),1) + &
2689  ( p(1))*self%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))
2690  end if
2691  end if
2692  if (target%id==io%id_debug .and. io%verbose>2) &
2693  write(io%output,1) target%id, ii, iv, self%id, j, jt, p, pt, &
2694  target%mem(ii(1),ii(2),ii(3),iv,target%it,1)
2695  1 format(i6,1x,3i3,i4,2x,"DBG patch_t%interpolate src:",i6, &
2696  " j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2," out:",1p,e15.6)
2697  !call trace%end (itimer)
2698 END SUBROUTINE interpolate
2699 
2700 !===============================================================================
2701 !> Interpolate in space and time, for a case where it is safe = the source point
2702 !> is inside the valid domain (inside a distance of half a cell from the source
2703 !> patch boundary for non-staggered points)
2704 !===============================================================================
2705 SUBROUTINE interpolate_specific (self, target, ii, iv, jt, pt)
2706  class(patch_t) :: self
2707  class(patch_t), pointer :: target
2708  integer :: ii(3), iv, jt(2)
2709  real :: pt(2)
2710  real :: tmp(2,2,2,2)
2711  !.............................................................................
2712  real(8) :: pos(3)
2713  integer :: j(3)
2714  real :: p(3)
2715  integer, save :: itimer=0
2716  !-----------------------------------------------------------------------------
2717  pos = target%myposition (ii, iv)
2718  call self%index_space (pos, iv, j, p)
2719  !call trace%begin ('patch_t%interpolate', itimer=itimer)
2720  !-----------------------------------------------------------------------------
2721  ! Pick up 16 values for tri-linear interpolation, divided by density
2722  !-----------------------------------------------------------------------------
2723  associate(id=>self%idx%d)
2724  tmp(:,:,:,:) = self%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,iv,jt(1:2),1) &
2725  / self%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,id,jt(1:2),1)
2726  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
2727  pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,1) + &
2728  ( p(1))*tmp(2,1,1,1)) + &
2729  ( p(2))*((1.0-p(1))*tmp(1,2,1,1) + &
2730  ( p(1))*tmp(2,2,1,1))) + &
2731  ( p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,1) + &
2732  ( p(1))*tmp(2,1,2,1)) + &
2733  ( p(2))*((1.0-p(1))*tmp(1,2,2,1) + &
2734  ( p(1))*tmp(2,2,2,1)))) + &
2735  pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,2) + &
2736  ( p(1))*tmp(2,1,1,2)) + &
2737  ( p(2))*((1.0-p(1))*tmp(1,2,1,2) + &
2738  ( p(1))*tmp(2,2,1,2))) + &
2739  ( p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,2) + &
2740  ( p(1))*tmp(2,1,2,2)) + &
2741  ( p(2))*((1.0-p(1))*tmp(1,2,2,2) + &
2742  ( p(1))*tmp(2,2,2,2))))
2743  !-----------------------------------------------------------------------------
2744  ! Multiply by target denssity
2745  !-----------------------------------------------------------------------------
2746  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
2747  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) * &
2748  target%mem(ii(1),ii(2),ii(3),id,target%it,1)
2749  end associate
2750  if (target%id==io%id_debug .and. io%verbose>2) &
2751  write(io%output,1) target%id, ii, iv, self%id, j, jt, p, pt, &
2752  target%mem(ii(1),ii(2),ii(3),iv,target%it,1)
2753  1 format(i6,1x,3i3,i4,2x,"DBG patch_t%interpolate_specific src:", &
2754  i6," j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2," out:",1p,e15.6)
2755  !call trace%end (itimer)
2756 END SUBROUTINE interpolate_specific
2757 
2758 SUBROUTINE update_position (self)
2759  class(patch_t):: self
2760  !-----------------------------------------------------------------------------
2761  if (all(self%velocity == 0.0)) return
2762  call trace%begin("patch_t%update_position")
2763  !
2764  self%position = self%position + self%velocity * self%dtime
2765  self%llc_cart = self%llc_cart + self%velocity * self%dtime
2766  self%centre_cart = self%centre_cart + self%velocity * self%dtime
2767  !
2768  self%mesh(:)%p = self%position(:)
2769  self%mesh(:)%llc_cart = self%llc_cart(:)
2770  !
2771  ! adjust the patch position for periodic wrapping
2772  where(self%position > 0.5*self%box)
2773  self%position = self%position - self%box
2774  self%mesh(:)%p = self%mesh(:)%p - self%box(:)
2775  self%llc_cart = self%llc_cart - self%box
2776  self%mesh(:)%llc_cart = self%mesh(:)%llc_cart - self%box
2777  self%centre_cart = self%centre_cart - self%box
2778  end where
2779  call trace%end()
2780 END SUBROUTINE update_position
2781 
2782 !===============================================================================
2783 SUBROUTINE custom_refine (self)
2784  class(patch_t):: self
2785  !-----------------------------------------------------------------------------
2786  call io%abort ("patch_t:: custom refine is triggered, but none provided. aborting.")
2787 END SUBROUTINE custom_refine
2788 
2789 !===============================================================================
2790 !> Check overlap and return overlap range in a patch, in terms self indices
2791 !===============================================================================
2792 FUNCTION get_overlap (self, patch, guards, l, u) RESULT(overlap)
2793  class(patch_t):: self
2794  class(patch_t), pointer:: patch
2795  integer:: l(3), u(3)
2796  logical:: guards, overlap
2797  !.............................................................................
2798  real(8):: dist(3), la(3), ua(3)
2799  !-----------------------------------------------------------------------------
2800  call trace%begin ('download_t%overlap_range')
2801  overlap = self%overlaps (patch)
2802  if (overlap) then
2803  !---------------------------------------------------------------------------
2804  ! First take care of mapping the distance between the patch and self
2805  ! centers, periodically, so this can be ignored in the next steps.
2806  !---------------------------------------------------------------------------
2807  dist = patch%position - self%position
2808  where (patch%periodic)
2809  dist = modulo(dist+0.5_8*patch%mesh%b,patch%mesh%b) - 0.5_8*patch%mesh%b
2810  end where
2811  !---------------------------------------------------------------------------
2812  ! Add the distances from the patch center to the start and end of guard
2813  ! zones, and convert to indices in self index space, making sure the range
2814  ! is legal
2815  !---------------------------------------------------------------------------
2816  if (guards) then
2817  la = dist - patch%size/2d0 - (patch%mesh%ng+1)*patch%mesh%d
2818  ua = dist + patch%size/2d0 + (patch%mesh%ng+1)*patch%mesh%d
2819  else
2820  la = dist - patch%size/2d0
2821  ua = dist + patch%size/2d0
2822  end if
2823  l = floor(la/self%mesh%d + self%mesh%o)
2824  u = ceiling(ua/self%mesh%d + self%mesh%o)
2825  l = max(l,self%mesh%lb)
2826  u = min(u,self%mesh%ub)
2827  !---------------------------------------------------------------------------
2828  ! Check that there really is overlap with the authoritative region
2829  !---------------------------------------------------------------------------
2830  if (any(l > self%mesh%ui .or. u < self%mesh%li .or. l > u)) then
2831  write (stderr,'(a,2(3x,a,i6,":",i2),3x,a,3f8.3,3x,a,2(3i4,2x))') &
2832  'patch_t%overlap_range false overlap, with ', &
2833  'self:', self%id, self%level, &
2834  'patch:', patch%id, patch%level, &
2835  'dist:', dist/self%ds, &
2836  'l, u:', l, u
2837  overlap = .false.
2838  end if
2839  !---------------------------------------------------------------------------
2840  ! Account for collapsed dimensions
2841  !---------------------------------------------------------------------------
2842  where (self%mesh%n == 1)
2843  l = min(l,self%mesh%ub)
2844  u = max(u,self%mesh%lb)
2845  end where
2846  end if
2847  call trace%end()
2848 END FUNCTION get_overlap
2849 
2850 !===============================================================================
2851 SUBROUTINE init_level (self)
2852  class(patch_t):: self
2853  !-----------------------------------------------------------------------------
2854  self%level = maxval(nint(log(self%box/self%ds) / &
2855  log(real(self%refine_ratio))), mask=self%n > 1)
2856 END SUBROUTINE init_level
2857 
2858 END MODULE patch_mod
2859 
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
Definition: timer_mod.f90:11
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
This module is a placeholder for shared information, solving the problem of making a piece of informa...
Definition: shared_mod.f90:7
Boundary conditions for centered variables, which have the physical boundary in-between mesh points l...
Module holding anonymous pointers back to extras features.
Definition: connect_mod.f90:4
Template module for mesh.
Definition: mesh_mod.f90:4
This index file has slot indices for all solver, all initially equal to zero It is the responsibility...
Definition: index_mod.f90:7
Definition: io_mod.f90:4
The lock module uses nested locks, to allow versatile use of locks, where a procedure may want to mak...
Template module for tasks.
Definition: task_mod.f90:4