DISPATCH
refine_mod.f90
1 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
2 !> This module handles checking max change between neighboring points. Each
3 !> instance of it needs an index to the variable, whether to take the log, and
4 !> a value for the max allowed change. It might be applied to log density, log
5 !> pressure, log opacity, velocity amplitude, magnetic field amplitude, etc.
6 !>
7 !> task_list_t%update
8 !> check_current
9 !> refine_needed
10 !> selective_refine
11 !> if ... refine
12 !> make_child_patch
13 !> tasklist%lock%set
14 !> init_nbors ...
15 !> check_nbors ...
16 !> tasklist%lock%unset
17 !> if ... derefine
18 !> remove_and_reset
19 !> tasklist%lock%set
20 !> init_nbors ...
21 !> check_nbors ...
22 !> tasklist%lock%unset
23 !> download_link
24 !> update
25 !> ...
26 !>
27 !> Support criterion
28 !> -----------------
29 !> A patch should be able to load guard zone values from nbors differing at
30 !> most by one AMR-level, to avoid too abrupt changes of resolution at patch
31 !> boundaries. The criterion to test must be formulated as a criterion on
32 !> the patch that should potentially be created, and not as a criterion on
33 !> the that needs the support (to avoid having to find a sub-region of the
34 !> parent of the parent -- and thus have to maintain parent relationships)
35 !>
36 !> On the other hand, the criterion needs to be checked at the point in time
37 !> when a new, higher level patch is to be created, to check if this requires
38 !> that an nbor patch should also be refined. The typical situation to look
39 !> for is when creating a patch at level L+1 requires of the patch whose sub-
40 !> region is being investigated requires a new patch at level L -- a so far
41 !> non-existing nbor of the patch at level L under investigation. The logic
42 !> is, in brief:
43 !>
44 !> for candidate in sub_regions_of (patch):
45 !> if needed(candidate):
46 !> attach_preliminary_nbor_list_to (candidate)
47 !> for nbor in candidate.nbors:
48 !> if nbor.patch.level==patch_candidate+2:
49 !> for nbor_candidate in sub_region_of (nbor.patch):
50 !> if overlap_between (candidate,nbor_candidate):
51 !> create(nbor_candidate)
52 !> candidate.add_nbor (nbor_candidate)
53 !> update_nbor_lists
54 !>
55 !> Note that this procedure does not assume a specific (e.g. side-by-side)
56 !> arrangement of patches, and hence allows for moving patches. Thus, it also
57 !> does not rely on the existence of a parent-child relation between patches.
58 !>
59 !> OpenMP synchronization is implemented with two (nested) OMP locks; one
60 !> (check_lock) responsible for synchronizing calls to list_t%check_nbors,
61 !> and one (download_lock) responsible for synchronizing filling of guard
62 !> zones (as well as restrict and prolong operations).
63 !>
64 !> When this module decides to refine or derefine a task, it acquires both
65 !> locks, and hence blocks check_nbors and download_link() calls while the
66 !> operation is on-going. Since no time advance is involved in the refine
67 !> / derefine, the blocking cannot cause a hang.
68 !>
69 !> Since refine/derefine is not done every timestep, the impact on running
70 !> tasks from the locking is not severe, and since two different locks are
71 !> used, one avoids the significant impact that would be the result of
72 !> locking out check_nbors and download_link from occuring at the same time.
73 !>
74 !> The impact on download_link could be entirely removed by always using
75 !> the sorted nbor list, and by collecting patches to be removed on a GC
76 !> list, removing them a few AMR cycles later.
77 !>
78 !> The absolute minimum amount of synchronization would be to just make sure
79 !> that an nbor list 1) can never change exactly when it is accessed, 2) to
80 !> immediately make a copy of it, 3) to make sure no to remove any task on the
81 !> nbor list while it is being used, and 4) to make sure that rotating the
82 !> memory slots will not invalidate using the nbor list, as is.
83 !>
84 !> An nbor list relevant for a particular task update is accessed in the
85 !> following order:
86 !>
87 !> 1) During a link_t%check_nbors() call, made after an nbor of the task in
88 !> question has been updated, link_t%check_ready() is called, and runs
89 !> through the tasks nbor list, comparing task%time values, finding that
90 !> all nbors of the task are advanced enough to provide boundary zone data,
91 !> or other data that the task needs. The task is then added to the ready
92 !> queue. The comparisons of times requires only that task_t%time is updated
93 !> atomically. These updates happen in task_t%rotate, which also manipulates
94 !> the task_t%iit(:) array, which holds the indices to the memory slots.
95 !>
96 !> 2) After a thread picks that task off the head of the ready queue, the first
97 !> action (in task_list_t%update) is to check if the task should be refined
98 !> or derefined. The test itself does not use the nbor lists, but if either
99 !> the current thread, or a thread working on an nbor decides to refine or
100 !> derefine, then the nbor list of the task may become modified.
101 !>
102 !> a) If it is refined, the data upon which the decision to use it as a
103 !> source of guard zone data still exists (in the parent task), and the
104 !> refined data in the child task(s) is initially just a prolongation of the
105 !> parent data, so does not yet provide any improvement. Therefore, there
106 !> is no problem, and the insertion of the new nbor list should be delayed
107 !> until it is entirely complete, at which point the insertion can be done
108 !> using a very brief locking of the link, while changeing task%link%nbor to
109 !> point to the new first nbor.
110 !>
111 !> b) If it is derefined, one just has to leave the nbor list and the task
112 !> there until after they have been used, in guard zone download, or in
113 !> check_ready().
114 !>
115 !> 3) The next step is download of data, which relies on data from all source
116 !> tasks, one of which could possible be in the process of being derefined.
117 !> To minimize the acess time needed, the download procedure should copy the
118 !> nbor list, a step that is anyway needed to sort the list into level order.
119 !> Then, to guarantee existence of the derefined patch, it has to either be
120 !> prevented from going away by locking, or else it must be put in a garbage
121 !> collector, which only removes patches when they are no longer needed.
122 !>
123 !> The latter option is the easier one, since it avoid messy synchronization.
124 !> A sync would require that derefinement is prevented in the interval btw
125 !> the task has been put in the ready queue, and until all nbors that might
126 !> possibly need it has done their downloads; a tall order.
127 !>
128 !> As part of putting the task in the ready queue, a copy of the nbor list
129 !> should be made -- this could be exactly the sorted copy that will be needed
130 !> in the download step in any case.
131 !>
132 !> Keeping track of "no longer needed" can be done with an atomically updated
133 !> counter in the task, which starts out with a value of one, and is incremented
134 !> when a task is put on a sorted, temporary nbor list, pending use in download.
135 !> When used as a source in download, the counter is decremented, endingg up on
136 !> one again, unless the task has been derefined away in the meantime, in which
137 !> case the counter ends on zero, and the task is then deallocated / deleted.
138 !>
139 !> Updating an nbor list with list_t%init_nbors() should be done in two steps:
140 !> first the head of the list is allocated, without linking it to link%nbor,
141 !> and the rest of the tasks are appended (or prepended). Then, in an lock
142 !> protected assignment, link%nbor is pointed to the head of the new nbor list.
143 !> Before that, a link is saved to the head of the existing nbor list, and after
144 !> the atomic reassignment, the old list is deleted.
145 !>
146 !> A derefined task is immediately removed from the task list, and hence can no
147 !> longer be added to re-initialized nbor lists. It thus ceases to be needed
148 !> when the last task that still has it on its nbor list copy deallocates its
149 !> nbor list, after using it in dowload_link().
150 !>
151 !> Summary: A sorted copy of the nbor list is created in connection with the
152 !> queue_by_time() call in check_ready(). At the same time a counter task_t%
153 !> n_needed is incremented for all tasks in the temporary nbor list. When the
154 !> task is used as a source in download_link() the counter is decremented
155 !> (atomically), and if the task is derefined the counter is decremented there
156 !> as well, leading to a zero cound when the task is no longer on an nbor list
157 !> copy. Then, and only then, the task may be deleted.
158 !>
159 !> FIXME: Care should be extended also to updates over MPI.
160 !> FIXME: Use of the deprecated shared_mod should be removed
161 !>
162 !> MPI specific handling: Boundary patches are always sent to vnbors, and if
163 !> they are new on the recieving rank, new nbor lists are generated, so no new
164 !> actions should be needed there. If a boundary task is removed, the vnbors
165 !> are send a copy with the "remove" bit set, and they should react correspond-
166 !> ingly, by doing essentially what remove_and_reset does.
167 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
169  USE iso_fortran_env, only: int8
170  USE io_mod
171  USE trace_mod
172  USE omp_mod
173  USE omp_timer_mod
174  USE omp_lock_mod
175  USE mpi_mod
176  USE bits_mod
177  USE kinds_mod
178  USE link_mod
179  USE list_mod
180  USE task_mod
181  USE patch_mod
182  USE solver_mod
183  USE extras_mod
184  USE shared_mod
185  USE download_mod
186  USE timer_mod
187  USE data_io_mod
188  implicit none
189  private
190  !
191  type:: refine_t
192  logical:: on ! allow external access
193  integer:: levelmin=0 ! Min level of refinment
194  integer:: levelmax=0 ! Max level of refinment
195  integer:: ratio=2 ! default
196  real:: safe_ratio=2.1 ! safe derefine ratio
197  contains
198  procedure:: init
199  procedure:: refine_needed
200  procedure:: check_current
201  procedure:: check_support
202  procedure:: need_to_support
203  procedure:: check_cubical
204  procedure:: check_jeans
205  procedure:: shock_and_cd_detector
206  procedure:: refine_region
207  procedure:: gradient_detector
208  procedure:: sanity_check
209  end type
210  !
211  integer:: check_interval=10 ! Number of steps betwen checks
212  real:: min_dx=0.0, max_dx=1e32 ! Min and max grid size
213  real:: refine_from_t=0.0 ! Time limits for occurence of refinement
214  real:: refine_until_t=1e30 ! Time limits for occurence of refinement
215  real:: max_jeans=0.0 ! Jeans number below which a patch is flagged for refinement
216  real:: min_jeans=0.0 ! Jeans number above which a patch is flagged for derefinement
217  real:: max_shock=0.0, max_contact=0.0 ! Tolerance values for shock/discontinuity detector refinement
218  real:: min_shock=0.0, min_contact=0.0 ! Tolerance values for shock/discontinuity detector derefinement
219  real:: min_compress=0.0, max_compress=0.0 ! Tolerance values for compression (negative divergence)
220  real:: min_vorticity=0.0, max_vorticity=0.0 ! Tolerance values for vorticity
221  integer, parameter:: ngradvar=10 ! Number of gradient variables
222  real, dimension(3):: region_llc=0.0 ! Sometimes you want to refine a region (manually)
223  real, dimension(3):: region_urc=0.0
224  real, dimension(ngradvar):: max_grad=0.0 ! Tolerance values for gradient detector refinement
225  real, dimension(ngradvar):: min_grad=0.0 ! Tolerance values for gradient detector derefinement
226  character(len=32), dimension(ngradvar):: grad_var=""
227  logical:: on=.false. ! Flag to turn the entire process on or off
228  logical:: force_cubical=.false. ! Force cubical patches
229  logical:: detailed_timer=.false. ! Lump all AMR together when false
230  integer:: n_locks=1 ! number of locks to use
231  integer:: verbose=0 ! Noise level
232  type(refine_t), public:: refine
233 CONTAINS
234 
235 !===============================================================================
236 !> Initialize refinement parameters (called from task_list_t%init_levels)
237 !===============================================================================
238 SUBROUTINE init (self)
239  class(refine_t):: self
240  integer:: iostat
241  integer, save:: levelmin=0, levelmax=0
242  integer, save:: ratio=2 ! How does a patch split when it is refined
243  logical, save:: first_time = .true.
244  namelist /refine_params/ verbose, on, levelmin, levelmax, check_interval, &
245  ratio, refine_from_t, refine_until_t, min_dx, max_dx, max_jeans, min_jeans, &
246  max_shock, min_shock, max_contact, min_contact, max_grad, min_grad, &
247  max_vorticity, min_vorticity, max_compress, min_compress, &
248  grad_var, region_llc, region_urc, force_cubical, n_locks, detailed_timer
249  character(len=120):: id = &
250  '$Id: aeec0874d565557a5ada56517633008870dafc5f $ tasks/refine_mod.f90'
251  !-----------------------------------------------------------------------------
252  call trace%begin ('refine_t%init')
253  call trace%print_id (id)
254  !$omp critical (input_cr)
255  if (first_time) then
256  first_time=.false.
257  levelmin = shared%levelmin
258  levelmax = shared%levelmax
259  rewind(io%input)
260  read (io%input, refine_params, iostat=iostat)
261  write (io%output, refine_params)
262  if (on) then
263  omp_lock%links = .true.
264  else
265  omp_lock%links = .false.
266  end if
267  end if
268  !-----------------------------------------------------------------------------
269  ! If AMR is active, we must use check_filled, since nbor levels are in
270  ! decreasing order, to save computing time
271  !-----------------------------------------------------------------------------
272  !$omp end critical (input_cr)
273  self%levelmin = levelmin
274  self%levelmax = levelmax
275  shared%levelmin = levelmin
276  shared%levelmax = levelmax
277  self%ratio = ratio
278  self%on = on
279  !-----------------------------------------------------------------------------
280  ! Initialize the lock in the garbage collector (shared by link_mod)
281  !-----------------------------------------------------------------------------
282  call garbage%init
283  call trace%end()
284 END SUBROUTINE init
285 
286 !===============================================================================
287 !> Check if refinement is desired and needed for the current task
288 !===============================================================================
289 FUNCTION refine_needed (self, tasklist, link) RESULT (refine)
290  class(refine_t):: self
291  class(list_t):: tasklist
292  class(link_t), pointer:: link
293  integer:: refine
294  !.............................................................................
295  class(solver_t), pointer:: patch
296  class(extras_t), pointer:: extras
297  real:: dxmin, dxmax
298  integer, save:: itimer=0
299  !----------------------------------------------------------------------------
300  refine = 0
301  if (.not. on) return
302  call trace%begin('refine_t%needed', itimer=itimer, detailed_timer=detailed_timer)
303  !----------------------------------------------------------------------------
304  ! Refinement criteria, starting with assumed derefinement
305  !------------------------------------------------------------------------
306  patch => cast2solver(link%task)
307  refine = -1
308  !
309  ! Extras refinement
310  extras => patch
311  refine = max(refine, extras%check_refine (extras))
312  if (verbose > 3) write (io_unit%output,*) &
313  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'extras'
314  !
315  ! Ensure consistent, cubical patch size (e.g. after RAMSES startup)
316  if (force_cubical) refine = max(refine, self%check_cubical(patch))
317  if (verbose > 3) write (io_unit%output,*) &
318  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'cubical'
319  !
320  ! Jeans criterion
321  if (max_jeans > 0.0) refine = max(refine, self%check_Jeans(patch))
322  if (verbose > 3) write (io_unit%output,*) &
323  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'Jeans'
324  !
325  ! shock and/or contact discontinuity detector
326  if (max_shock > 0.0 .or. max_contact > 0.0) &
327  refine = max(refine, self%shock_and_cd_detector(patch))
328  if (verbose > 3) write (io_unit%output,*) &
329  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'shock'
330  !
331  ! compression detector
332  if (max_compress > 0.0) &
333  refine = max(refine, check_compress(self, patch))
334  if (verbose > 3) write (io_unit%output,*) &
335  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'compress'
336  !
337  ! vorticity detector
338  if (max_vorticity > 0.0) &
339  refine = max(refine, check_vorticity(self, patch))
340  if (verbose > 3) write (io_unit%output,*) &
341  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'vorticity'
342  !
343  ! gradient detector
344  if (any(max_grad > 0.0)) refine = max(refine, self%gradient_detector(patch))
345  if (verbose > 3) write (io_unit%output,*) &
346  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'gradient'
347  !
348  ! refine on (a) user-specified region(s)
349  if (any(region_llc /= 0.0) .and. any(region_urc /= 0.0)) &
350  refine = max(refine, self%refine_region(patch))
351  if (verbose > 3) write (io_unit%output,*) &
352  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'user'
353  !
354  ! Check if support constraints require refinement
355  refine = max(refine, need_to_support(self, tasklist, link))
356  if (verbose > 3) write (io_unit%output,*) &
357  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'need to support'
358  !
359  ! Check if support constraints allow derefinement
360  refine = max(refine, check_derefine_support(self, tasklist, link))
361  if (verbose > 3) write (io_unit%output,*) &
362  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'derefine support'
363  !
364  call trace%end (itimer, detailed_timer=detailed_timer)
365 END FUNCTION refine_needed
366 
367 !===============================================================================
368 !> Check if refinement is desired and needed for the current task; if so, push
369 !> new tasks onto the (ready) queue with the same time as the parent, i.e., to
370 !> the head of the queue. If the task is a leaf task, and the level of refinement
371 !> does not appear to be needed, request derefinement.
372 !===============================================================================
373 SUBROUTINE check_current (self, tasklist, link, was_refined, was_derefined)
374  class(refine_t):: self
375  class(list_t):: tasklist
376  class(link_t), pointer:: link, nbor
377  class(solver_t), pointer:: patch
378  logical, intent(out):: was_refined, was_derefined
379  real:: dxmin, dxmax
380  integer:: refine, id, n_added, l(3), u(3)
381  integer, save:: itimer=0
382  !.............................................................................
383  was_refined = .false.
384  was_derefined = .false.
385  !-----------------------------------------------------------------------------
386  ! Punt on refinement if the task is very new, or has the init_nbors bit set
387  !-----------------------------------------------------------------------------
388  if (.not.on .or. &
389  link%task%istep < 3 .or. &
390  link%task%is_set (bits%init_nbors)) &
391  return
392  !---------------------------------------------------------------------------
393  call trace%begin('refine_t%check_current', itimer=itimer)
394  patch => cast2solver(link%task)
395  id = patch%id
396  if (verbose > 3) &
397  write (io_unit%output,'(i6,2x,a,f12.6,2i4)') &
398  link%task%id, 'refine_t%check_current, at t =', &
399  link%task%time, link%task%level, self%levelmax
400  !----------------------------------------------------------------------------
401  dxmin = minval(patch%ds)
402  dxmax = maxval(patch%ds)
403  if (verbose > 3) then
404  write (io_unit%output,'(a,i6,2i4,l4,1p,4e11.2)') &
405  'refine_t%refine_needed: id, istep, refine_next =', &
406  patch%id, patch%istep, patch%check_refine_next, &
407  patch%is_set(bits%virtual), min_dx, dxmin, max_dx, dxmax
408  end if
409  !----------------------------------------------------------------------------
410  ! Filter out patches that do not match the necessary criteria
411  !----------------------------------------------------------------------------
412  if (patch%is_clear(bits%virtual) .and. &
413  dxmin > min_dx .and. &
414  dxmax <= max_dx .and. &
415  patch%time >= refine_from_t .and. &
416  patch%time < refine_until_t .and. &
417  patch%istep >= patch%check_refine_next) then
418  if (verbose > 3) &
419  tasklist%verbose = verbose
420  !---------------------------------------------------------------------------
421  ! Before testing any criteria, check if patch%irefine is set
422  !---------------------------------------------------------------------------
423  if (.not.allocated(patch%irefine)) &
424  allocate (patch%irefine(patch%gn(1),patch%gn(2),patch%gn(3)))
425  patch%irefine = -1
426  l = patch%mesh%li
427  u = patch%mesh%ui
428  refine = maxval(patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)))
429  if (verbose > 1 .and. &
430  any(patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) > 0)) then
431  write (stdout,'(3x,a,i6,":",i2,2x,a)') &
432  'refine' , patch%id, patch%level, &
433  'requested for support'
434  end if
435  !---------------------------------------------------------------------------
436  ! Check all refinement criteria. If patch%irefine is already marked for
437  ! refinement, then it has been done to request support of new level+2 patches
438  !---------------------------------------------------------------------------
439  refine = max(refine,self%refine_needed (tasklist, link))
440  !---------------------------------------------------------------------------
441  if (verbose > 3) write (io_unit%output,*) &
442  patch%id, refine, minval(patch%irefine), maxval(patch%irefine), 'mask'
443  !---------------------------------------------------------------------------
444  ! Apply refinement criteria on potential child patch parts, but only if
445  ! either refinement has been indicated, or if derfinement is allowed
446  !---------------------------------------------------------------------------
447  if (refine /= 0) then
448  call selective_refine (self, tasklist, patch, refine, was_refined, was_derefined)
449  if (was_derefined) then
450  call trace%end (itimer)
451  return
452  end if
453  end if
454  !---------------------------------------------------------------------------
455  ! Update next refinement step counter and deallocate local irefine
456  !---------------------------------------------------------------------------
457  patch%check_refine_next = patch%istep + check_interval
458  deallocate (patch%irefine)
459  end if
460  call trace%end (itimer)
461 END SUBROUTINE check_current
462 
463 !===============================================================================
464 !> Loop over the ratio x ratio x ratio positions and determine whether to refine
465 !> or derefine.
466 !===============================================================================
467 SUBROUTINE selective_refine (self, tasklist, patch, refine, was_refined, was_derefined)
468  class(refine_t):: self
469  class(list_t):: tasklist
470  class(solver_t), pointer:: patch
471  integer:: refine
472  logical:: was_refined, was_derefined
473  !.............................................................................
474  class(solver_t), pointer:: child
475  class(link_t), pointer:: nbor, link, nbors, tail
476  class(patch_t), pointer:: nbpatch
477  integer:: i, j, k
478  integer, dimension(3):: l, u, n, ii
479  real(8), dimension(3):: size, pos
480  logical:: not_already_refined
481  integer:: refine_it, n_added
482  integer(int8), allocatable:: refined(:,:,:)
483  integer:: itimer=0
484  !-----------------------------------------------------------------------------
485  call trace%begin ('refine_t%selective_refine', itimer=itimer, detailed_timer=detailed_timer)
486  was_refined = .false.
487  was_derefined = .false.
488  link => patch%link
489  call link%lock%set ('selective_refine')
490  nbors => link%nbor
491  tail => tasklist%tail
492  nbor => patch%link%nbor
493  !-----------------------------------------------------------------------------
494  ! As preparation, mark all points that already are refined
495  !-----------------------------------------------------------------------------
496  allocate (refined(patch%gn(1),patch%gn(2),patch%gn(3)))
497  refined = 0
498  do while (associated(nbor))
499  nbpatch => nbpatch%cast2patch(nbor%task)
500  if (patch%contains(nbpatch%position) .and. &
501  patch%level == nbpatch%level-1) then
502  l = patch%index_only (nbpatch%position + nbpatch%mesh%lf)
503  u = patch%index_only (nbpatch%position + nbpatch%mesh%uf)
504  !-------------------------------------------------------------------------
505  ! A patch with existing child patches should not be derefined. To prevent
506  ! this, refined regions inside are set to at least 0 in irefine, and to
507  ! prevent creating doubles, the refined array is set to 1
508  !-------------------------------------------------------------------------
509  refined(l(1):u(1),l(2):u(2),l(3):u(3)) = 1
510  patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) = &
511  max(0,patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)))
512  !-------------------------------------------------------------------------
513  if (verbose > 2) then
514  write (io_unit%output,'(a,i6,a,i6,2(3x,3i4))') &
515  'derefine: patch', patch%id, ' contains', nbpatch%id, l, u
516  end if
517  !-------------------------------------------------------------------------
518  end if
519  nbor => nbor%next
520  end do
521  call link%lock%unset ('selective_refine')
522  !-----------------------------------------------------------------------------
523  ! First, consider if any part should be refined
524  !-----------------------------------------------------------------------------
525  if (patch%level < self%levelmax .and. refine > 0) then
526  if (any(patch%irefine > 0)) then
527  do k=1,self%ratio
528  do j=1,self%ratio
529  do i=1,self%ratio
530  !-------------------------------------------------------------------
531  ! Size and position
532  !-------------------------------------------------------------------
533  size = patch%size/self%ratio
534  pos = patch%position + ([i,j,k]-(self%ratio+1)*0.5)*size
535  !-------------------------------------------------------------------
536  ! Lower and upper index in parent patch
537  !-------------------------------------------------------------------
538  n = patch%mesh%n/self%ratio
539  l = patch%mesh%li + ([i,j,k]-1)*n
540  u = l + n - 1
541  !-------------------------------------------------------------------
542  ! Should this part be refined or not?
543  !-------------------------------------------------------------------
544  refine_it = maxval(patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)))
545  !-------------------------------------------------------------------
546  ! Is this part refined before, or should it be now?
547  !-------------------------------------------------------------------
548  ii = (l+u)/2
549  not_already_refined = refined(ii(1),ii(2),ii(3)) == 0
550  if (not_already_refined .and. refine_it > 0) then
551  !io%do_trace = .true.
552  !-----------------------------------------------------------------
553  ! Create a refined patch here
554  !-----------------------------------------------------------------
555  if (verbose > 1) &
556  write (io_unit%output,'(a,i6,1p2g16.5,2x,3(i4,":",i2))') 'refine', &
557  patch%id, patch%fmaxval(patch%mem(:,:,:,patch%idx%d,patch%it,1)), &
558  maxval(patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),patch%idx%d,patch%it,1)), &
559  l(1),u(1),l(2),u(2),l(3),u(3)
560  call make_child_patch (self, tasklist, patch, pos, size, child)
561  !-----------------------------------------------------------------
562  if (verbose > 0) then
563  associate(unit => merge(stdout, io_unit%output, verbose>1))
564  write (unit,'(2(3x,a,i6,":",i2),2x,a,2l1)') &
565  'refine' , patch%id, patch%level, &
566  'to' , child%id, child%level, &
567  'BV:', child%is_set(bits%boundary), child%is_set(bits%virtual)
568  flush (unit)
569  end associate
570  end if
571  !-----------------------------------------------------------------
572  was_refined= .true.
573  end if
574  end do
575  end do
576  end do
577  end if
578  end if
579  !-----------------------------------------------------------------------------
580  ! Then, consider if the whole patch should be derefined. This requires that
581  ! it has no higher level patch inside (would trigger patch%refine > 0 there),
582  ! and secondly that irefine(:,:,:) indicates that derefine would be OK
583  !-----------------------------------------------------------------------------
584  if (refine == -1 .and. & ! derefine possible
585  patch%level > self%levelmin .and. & ! relevant?
586  patch%check_refine_next > 0 .and. & ! not immediately!
587  patch%istep >= 4 .and. & ! must take 4 steps
588  .not.was_refined) then ! not after a refine!
589  l = patch%mesh%li
590  u = patch%mesh%ui
591  if (all(patch%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) == -1)) then
592  !-----------------------------------------------------------------------
593  if (verbose > 0) then
594  associate(unit => merge(stdout, io_unit%output, verbose > 1))
595  write (io_unit%output,'(1x,a,i6,":",i2,2x,"B:",l1,3x,a,i6,3x,a,f12.6)') &
596  'derefine', patch%id, patch%level, patch%is_set (bits%boundary), &
597  'istep:', patch%istep, 'time:', patch%time
598  flush (unit)
599  end associate
600  end if
601  !-------------------------------------------------------------------------
602  ! If task is to be removed, call the comprehensive remove procedure, which
603  ! includes calling set_init_nbors(), queueing all nbors for immediate
604  ! update of their nbor lists.
605  !-------------------------------------------------------------------------
606  call data_io%update_counters (patch, -1)
607  call patch%log ('remove')
608  call tasklist%remove_and_reset (link)
609  was_derefined = .true.
610  end if
611  end if
612  deallocate (refined)
613  call trace%end (itimer, detailed_timer=detailed_timer)
614 END SUBROUTINE selective_refine
615 
616 !===============================================================================
617 !> Check if new%task is a patch, and if it overlaps with other level L-2 patches,
618 !> with points not covered by L-1 patches, indicating that a new L-1 patch should
619 !> be created from the L-2 patch.
620 !>
621 !> To implement this, we need to look at the regions that overlap with an L-2
622 !> patch, and then run through the nbor list L-1 patches, to see if they cover
623 !> the area. If not, "make an impression" on the L-2 patch, by setting the
624 !> corresponding cells in its %xrefine map.
625 !>
626 !> NOTE: This mechanism is currently NOT used; it is left for now, to be removed
627 !> when need_to_support() has been validated as being sufficient -- it looks
628 !> at the same issue from the point of view of the L-2 patch, which is better,
629 !> since it eliminates the need for synchronizing access to %xrefine.
630 !===============================================================================
631 SUBROUTINE check_support (self, tlist, new, n_added)
632  class(refine_t):: self
633  class(list_t):: tlist
634  class(link_t), pointer:: new
635  integer:: n_added
636  !.............................................................................
637  class(link_t), pointer:: old, nbor
638  class(patch_t), pointer:: old_patch, new_patch, nbor_patch
639  type(patch_t):: patch
640  integer:: l(3), u(3)
641  logical:: overlap, guards
642  integer, allocatable, dimension(:,:,:):: trefine
643  integer:: itimer=0
644  !-----------------------------------------------------------------------------
645  ! Return if there are less than two levels down to levelmin
646  !-----------------------------------------------------------------------------
647  n_added = 0
648  if (new%task%level <= self%levelmin+1) return
649  !-----------------------------------------------------------------------------
650  call trace%begin ('refine_t%check_support', itimer=itimer, detailed_timer=detailed_timer)
651 !write (io_unit%log,*) 'unexpected call to check_support for', new%task%id, &
652 !new%task%is_set(bits%boundary), new%task%is_set(bits%virtual)
653  !-----------------------------------------------------------------------------
654  ! Select patch-based tasks, and assume for now it is OK for updates
655  !-----------------------------------------------------------------------------
656  new_patch => patch%cast2patch (new%task)
657  if (associated(new_patch)) then
658  if (verbose > 2) &
659  write (stdout,*) 'check_support: id =', &
660  new_patch%id, new_patch%is_set (bits%support)
661  !---------------------------------------------------------------------------
662  ! Search task list for patch-based overlapping L-2 patches
663  !---------------------------------------------------------------------------
664  call new_patch%clear (bits%support)
665  call tlist%lock%set ('check_support')
666  old => tlist%head
667  do while (associated(old))
668  if (old%task%level <= new%task%level-2) then
669  old_patch => patch%cast2patch (old%task)
670  guards = .true.
671  if (old_patch%get_overlap (new_patch, guards, l, u)) then
672  if (verbose > 2) &
673  write (stdout,'(a,i6,3(2x,i2,":",i2))') &
674  ' level-2 overlap: id =', old_patch%id, &
675  l(1),u(1), l(2),u(2), l(3),u(3)
676  !---------------------------------------------------------------------
677  ! Allocate a temporary integer mask, to assemble evidence in
678  !---------------------------------------------------------------------
679  allocate (trefine(old_patch%gn(1),old_patch%gn(2),old_patch%gn(3)))
680  trefine = 0
681  trefine(l(1):u(1), l(2):u(2), l(3):u(3)) = 1
682  !---------------------------------------------------------------------
683  ! Find all L-1 nbors of the new patch, and cancel the refinement
684  !---------------------------------------------------------------------
685  nbor => new%nbor
686  do while (associated(nbor))
687  if (nbor%task%level == new%task%level-1) then
688  nbor_patch => patch%cast2patch (nbor%task)
689  guards = .false.
690  if (old_patch%get_overlap (nbor_patch, guards, l, u)) then
691  if (verbose > 2) &
692  write (stdout,'(a,i6,3(2x,i2,":",i2))') &
693  ' level-1 overlap: id =', nbor_patch%id, &
694  l(1),u(1), l(2),u(2), l(3),u(3)
695  trefine(l(1):u(1), l(2):u(2), l(3):u(3)) = 0
696  end if
697  end if
698  nbor => nbor%next
699  end do
700  !---------------------------------------------------------------------
701  ! If any refinement marking remains, emit messages, set bits%support
702  !---------------------------------------------------------------------
703  l = old_patch%li
704  u = old_patch%ui
705  if (any(trefine(l(1):u(1), l(2):u(2), l(3):u(3)) == 1)) then
706  call new_patch%set (bits%support)
707  if (verbose > 1) then
708  write (stdout,'(2(a,i6,":",i2,2x))') &
709  'id =', old_patch%id, old_patch%level, 'support request from', &
710  new_patch%id, new_patch%level
711  write (io_unit%log,'(2(a,i6,":",i2,2x))') &
712  'id =', old_patch%id, old_patch%level, 'support request from', &
713  new_patch%id, new_patch%level
714  end if
715  else if (verbose > 3) then
716  write (stdout,'(2(a,i6,":",i2,2x))') &
717  'id =', old_patch%id, old_patch%level, 'support not requested from', &
718  new_patch%id, new_patch%level
719  end if
720  deallocate (trefine)
721  end if
722  end if
723  old => old%next
724  end do
725  call tlist%lock%unset ('check_support')
726  end if
727  !-----------------------------------------------------------------------------
728  call trace%end (itimer, detailed_timer=detailed_timer)
729 END SUBROUTINE check_support
730 
731 !===============================================================================
732 !> Decide if there are level L+2 tasks we need to support be refining
733 !===============================================================================
734 FUNCTION need_to_support (self, tlist, link) RESULT (refine)
735  class(refine_t):: self
736  class(list_t):: tlist
737  class(link_t), pointer:: link
738  integer:: refine
739  !-----------------------------------------------------------------------------
740  class(link_t), pointer:: nbor1, nbor2, nbor3
741  class(patch_t), pointer:: task, nbtask
742  type(patch_t):: patch
743  logical:: guards
744  integer:: l(3), u(3)
745  !-----------------------------------------------------------------------------
746  call trace%begin ('refine_t%need_to_support')
747  task => patch%cast2patch(link%task)
748  !-----------------------------------------------------------------------------
749  ! Find L+2 and higher tasks, and check for overlap of their guard zones
750  !-----------------------------------------------------------------------------
751  call tlist%lock%set ('need_to_support')
752  nbor1 => tlist%head
753  do while (associated(nbor1))
754  if (nbor1%task%level >= link%task%level+2) then
755  nbtask => patch%cast2patch (nbor1%task)
756  guards = .true.
757  if (task%get_overlap (nbtask, guards, l, u)) then
758  task%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) = 1
759  end if
760  end if
761  nbor1 => nbor1%next
762  end do
763  !-----------------------------------------------------------------------------
764  ! Find L+1 nbor-nbors, check for overlap of their interiors
765  !-----------------------------------------------------------------------------
766  call link%lock%set ('need_to_support')
767  nbor1 => link%nbor
768  do while (associated(nbor1))
769  if (nbor1%task%level == link%task%level+1) then
770  nbtask => patch%cast2patch (nbor1%task)
771  guards = .false.
772  if (task%get_overlap (nbtask, guards, l, u)) then
773  task%irefine(l(1):u(1),l(2):u(2),l(3):u(3)) = 0
774  end if
775  end if
776  nbor1 => nbor1%next
777  end do
778  call link%lock%unset ('need_to_support')
779  call tlist%lock%unset ('need_to_support')
780  l = task%mesh%li
781  u = task%mesh%ui
782  refine = maxval(task%irefine(l(1):u(1),l(2):u(2),l(3):u(3)))
783  call trace%end ()
784 END FUNCTION need_to_support
785 
786 !===============================================================================
787 !> To check if a patch at level L may be removed w/o patches at level L+1 loosing
788 !> support it is enough to check if there are nbors at level L+1, because if
789 !> there are they would have, at the place where it had level L support, instead
790 !> only a level L-1 patch.
791 !===============================================================================
792 FUNCTION check_derefine_support (self, tlist, link, check_only) RESULT (refine)
793  USE task_mod
794  class(refine_t):: self
795  class(list_t):: tlist
796  class(link_t), pointer:: link
797  logical, optional:: check_only
798  integer:: refine
799  !.............................................................................
800  class(link_t), pointer:: nbor, nbors
801  integer:: n_nbor1, n_nbor2, n_added
802  integer, save:: itimer=0
803  !-----------------------------------------------------------------------------
804  refine = -1
805  if (link%task%level == self%levelmin) &
806  return
807  call trace%begin('refine_t%check_derefine_support', itimer=itimer, detailed_timer=detailed_timer)
808  !-----------------------------------------------------------------------------
809  call tlist%lock%set ('check_derefine_support')
810  nbor => link%nbor
811  do while (associated(nbor))
812  if (nbor%task%level == link%task%level+1) then
813  refine = 0
814  exit
815  end if
816  nbor => nbor%next
817  end do
818  if (verbose > 2) &
819  write (stdout,*) link%task%id, 'check_derefine_support =', refine
820  !-----------------------------------------------------------------------------
821  ! Release lock
822  !-----------------------------------------------------------------------------
823  call tlist%lock%unset ('check_derefine_support')
824  !-----------------------------------------------------------------------------
825  call trace%end (itimer, detailed_timer=detailed_timer)
826 END FUNCTION check_derefine_support
827 
828 !===============================================================================
829 !> Create a new patch using an existing "parent" patch.
830 !> This procedure is used frequently during refinement.
831 !===============================================================================
832 SUBROUTINE make_child_patch (self, tlist, parent, pos, child_size, &
833  child)
834  class(refine_t):: self
835  class(list_t) :: tlist
836  class(solver_t), pointer :: parent, child
837  real(8), dimension(3), intent(in) :: pos, child_size
838  !.............................................................................
839  class(link_t), pointer :: child_link, child_nbor
840  integer:: j, count, l(3), u(3), n_added
841  integer, save:: itimer=0
842  !-----------------------------------------------------------------------------
843  call trace%begin('refine_t%make_child_patch', itimer=itimer, detailed_timer=detailed_timer)
844  if (parent%rank /= mpi%rank) then
845  print *, 'make_child_patch: mpi%rank, parent%rank =', mpi%rank, parent%rank
846  call io%abort ('make_child_patch: wrong MPI rank')
847  end if
848  !-----------------------------------------------------------------------------
849  ! Allocate a new task, 1/nsplit the size of original one, and adjust its
850  ! position in each direction with "pos", in units of the new size.
851  !-----------------------------------------------------------------------------
852  ! Using sourced allocation below effectively clones the parent. This is
853  ! necessary, to pick up properties that may have been set at the experiment
854  ! level, which would otherwise be left undefined, or in inconsistent state
855  !-----------------------------------------------------------------------------
856  call parent%lock%set ('parent')
857  allocate (child, source=parent)
858  call parent%lock%unset ('parent')
859  call child%clone_mem_accounting
860  !-----------------------------------------------------------------------------
861  child%mem_allocated = .false. ! make sure to get ..
862  nullify (child%mem) ! .. new mem
863  nullify (child%mesh) ! .. new mesh
864  nullify (child%lock) ! .. new lock
865  nullify (child%mem_lock) ! .. new mem_lock
866  child%id = 0 ! ensure a new ID
867  child%size = child_size ! new size
868  child%position = pos ! patch position
869  child%llc_cart = child%position - & ! mid point position
870  merge(0.5,0.0,child%n>1) * child%size ! lower left corner
871  child%llc_nat = child%llc_cart ! FIXME!: this form only applies for Cartesian coords.
872  child%n_needed = 1 ! reset to initial
873  child%istep = 0 ! for diagnostics
874  !-----------------------------------------------------------------------------
875  ! Force consistent resolution in new patches, overriding the resolution
876  ! inherited from the parent patch.
877  !-----------------------------------------------------------------------------
878  child%n = shared%patch_n
879  !-----------------------------------------------------------------------------
880  ! Clear flags -- explictly to show which ones are relevant
881  !-----------------------------------------------------------------------------
882  !child%status = 0
883  call child%clear(bits%root_grid)
884  call child%clear(bits%ready)
885  call child%clear(bits%busy)
886  call child%clear(bits%frozen)
887  child%boundaries%bits = 0
888  !-----------------------------------------------------------------------------
889  ! Append links to the new patch in the patch list, sorted by quality
890  ! Neighbour lists are constructed at the end of `split_patch`.
891  !-----------------------------------------------------------------------------
892  allocate (child_link)
893  call child_link%init ! initialise the lock
894  child_link%task => child ! associate task with link
895  child_link%parent => parent%link ! used in remove_patch
896  child%link => child_link ! associate link with task
897  child%parent => parent ! task relation
898  !-----------------------------------------------------------------------------
899  ! Initialise the new patch
900  !-----------------------------------------------------------------------------
901  child%level = parent%level + 1 ! must be set before init!
902  call child%set (bits%frozen)
903  call child%init() ! patch, for a new id
904  call child%clear (bits%frozen)
905  child%level = parent%level + 1 ! override patch_t%init
906  child%time = parent%time ! same time
907  child%it = 1 ! time index
908  child%new = 2 ! new time slot
909  child%iit = 1 ! single time slice
910  child%iit(child%nt) = child%new ! new slot
911  child%dt(:) = parent%dt(parent%it) / self%ratio ! reset the time-step
912  child%t(:) = child%time ! sinle time
913  child%parentid = parent%id ! for debugging
914  child%check_refine_next = 0 ! allow recursive refine
915  child%iout = child%time/io%out_time+1 ! previous output index
916  child%out_next = child%iout*io%out_time ! next output time
917  child%dnload_time = -1.0 ! for duplicate call test
918  !-----------------------------------------------------------------------------
919  ! Increment I/O write counter, to avoid repeated writing of metadata
920  !-----------------------------------------------------------------------------
921  call data_io%update_counters (child, +1)
922  call child%log ('created')
923  !-----------------------------------------------------------------------------
924  if (verbose > 0) &
925  write (io_unit%log,'(f12.5,2x,a,2i6,1p,g14.5)') &
926  wallclock(), 'refine_t%make_child_patch, parent, child, time =', &
927  parent%id, child%id, child%time
928  !-----------------------------------------------------------------------------
929  ! Check if the new patch, which might have been created to support another
930  ! patch, possibly should be immediately frozen. If that is the case, it still
931  ! needs to be added to the task list, in case it is needed for support, but by
932  ! setting bits%frozen it is prevented from being added to the ready queue.
933  !-----------------------------------------------------------------------------
934  if (child%time > io%end_time) then
935  call child%set (bits%frozen)
936  end if
937  !-----------------------------------------------------------------------------
938  ! Set the bits%support here, which causes it to be set on the rank nbors
939  ! if this is a boundary task. There, a check for support is then triggered,
940  ! where the nbor lists also in the virtual nbors are refreshed. On the owner
941  ! rank the bits%support will remain until next time the task is updated, when
942  ! code in dispatcher0_t%update() sees the bit, call init_nbors(), and clears
943  ! the bit.
944  !-----------------------------------------------------------------------------
945  call child%set (bits%support)
946  !-----------------------------------------------------------------------------
947  ! Make parent a neighbor of the new child patch. This will be overwritten
948  ! in `init_nbor_nbors` but is necessary here to permit interpolation of
949  ! values from parent to child.
950  !-----------------------------------------------------------------------------
951  allocate (child_nbor)
952  child_nbor%link => parent%link
953  child_nbor%task => parent
954  nullify (child_link%nbor)
955  child_link%nbors_by_level => child_nbor
956  !------------------------------------------------------------------------
957  ! Download values from parent and neighbors, except at time = 0, to allow
958  ! recursive refinement there. If time = 0, patch should have been filled
959  ! with values during `child%init()` call.
960  !------------------------------------------------------------------------
961  !$omp atomic update
962  parent%n_needed = parent%n_needed + 2
963  if (parent%time > 0.0 .and. child%restart <= 0) then
964  call download%download_link (child_link, all_cells=.true.)
965  end if
966  !-----------------------------------------------------------------------------
967  ! Add the child task to the task list, including giving it an nbor list etc,
968  ! and add it also to the ready queue
969  !-----------------------------------------------------------------------------
970  call tlist%add_new_link (child_link)
971  call tlist%queue_by_time (child_link)
972  !-----------------------------------------------------------------------------
973  if (verbose > 1) then
974  l = parent%index_only (child%position + child%mesh%lf)
975  u = parent%index_only (child%position + child%mesh%uf)
976  l = max(l,parent%mesh%li)
977  u = min(u,parent%mesh%ui)
978  print 1,'make_child_patch: parent, level, iout, time, size(3), position(3) =', &
979  parent%id, parent%level, parent%time, parent%size, parent%position , &
980  maxval(parent%mem(l(1):u(1),l(2):u(2),l(3):u(3),parent%idx%d,parent%it,1)), &
981  l, u, parent%contains(child)
982  1 format(a,i6,i3,f9.6,2(1x,3f8.4),1pe12.2,2(1x,3i3),l2)
983  l(:) = parent%mesh(:)%li
984  u(:) = parent%mesh(:)%ui
985  print 1,' child, level, iout, time, size(3), position(3) =', &
986  child%id, child%level, child%time, child%size, child%position , &
987  maxval(child%mem(l(1):u(1),l(2):u(2),l(3):u(3),child%idx%d,child%it,1))
988  flush (io%output)
989  end if
990  !-----------------------------------------------------------------------------
991  call trace%end (itimer, detailed_timer=detailed_timer)
992 END SUBROUTINE make_child_patch
993 
994 !===============================================================================
995 !> Check that the patch dimensions are the desired cubical ones
996 !===============================================================================
997 FUNCTION check_cubical (self, patch) RESULT (refine_it)
998  class(refine_t):: self
999  class(solver_t), pointer:: patch
1000  integer:: refine_it
1001  !.............................................................................
1002  if (all(patch%n==shared%patch_n)) then
1003  refine_it = -1
1004  else
1005  refine_it = 1
1006  patch%msplit = 1
1007  end if
1008 END FUNCTION check_cubical
1009 
1010 !===============================================================================
1011 !> Refine / derefine on ratio of Jeans length to cell size (max_jeans/min_jeans)
1012 !===============================================================================
1013 FUNCTION check_jeans (self, patch) RESULT(refine_it)
1014  USE math_mod
1015  USE scaling_mod
1016  class(refine_t):: self
1017  class(solver_t), pointer:: patch
1018  integer:: refine_it
1019  !.............................................................................
1020  integer:: imin(3), l(3), u(3)
1021  real:: Jeans_num, dmin, c
1022  real, dimension(:,:,:), allocatable:: pg
1023  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d, e
1024  real, allocatable:: test(:,:,:)
1025  logical:: masked
1026  !-----------------------------------------------------------------------------
1027  call trace%begin ('refine_t%check_jeans')
1028  call self%sanity_check (min_jeans, max_jeans, 'check_Jeans')
1029  l = 1
1030  u = patch%gn
1031  d => patch%mem(l(1):u(1),l(2):u(2),l(3):u(3),patch%idx%d,patch%it,1)
1032  allocate(pg(patch%gn(1),patch%gn(2),patch%gn(3)))
1033  select type (patch)
1034  class is (solver_t)
1035  pg = patch%gas_pressure() + tiny(1.0)
1036  class default
1037  call mpi%abort("Unrecognised patch type. Abort!")
1038  end select
1039  dmin = patch%fminval(d)
1040  if (dmin <= 0.0) then
1041  imin = minloc(d)
1042  write (io%output,*) &
1043  'check_Jeans ERROR: id, dmin, imin =', patch%id, dmin, imin
1044  jeans_num = min_jeans
1045  else
1046  allocate (test(patch%gn(1),patch%gn(2),patch%gn(3)))
1047  c = minval(patch%mesh%d)*sqrt(scaling%grav/(math%pi*patch%gamma))
1048  !print *, patch%id, scaling%grav, 'grav'
1049  test = sqrt(pg)/(c*d)
1050  !print *, patch%id, patch%fminval(d), patch%fmaxval(d), 'min/max d'
1051  !print *, patch%id, patch%fminval(test), patch%fmaxval(test), 'min/max L_J/ds'
1052  !---------------------------------------------------------------------------
1053  ! The Jeans number in the unrefined part
1054  !---------------------------------------------------------------------------
1055  jeans_num = minval(test, mask=patch%irefine <= 0)
1056  !---------------------------------------------------------------------------
1057  ! Compute a per-cell criterion, which sets patch%irefine to +1 where it
1058  ! needs refinement, and sets it to 0 where the current resolution is needed.
1059  ! By implication, where it does not set patch%irefine, and where it remains
1060  ! at the intial -1 value, derefinement is allowed
1061  !---------------------------------------------------------------------------
1062  where (test < min_jeans)
1063  patch%irefine = +1
1064  !---------------------------------------------------------------------------
1065  ! Jeans number not small enough, or overlying refinement
1066  !---------------------------------------------------------------------------
1067  else where (test < max_jeans)
1068  patch%irefine = max(patch%irefine, 0)
1069  end where
1070  deallocate (test)
1071  refine_it = 0
1072  if (jeans_num > max_jeans) refine_it = -1
1073  if (jeans_num < min_jeans) refine_it = +1
1074  !---------------------------------------------------------------------------
1075  if (verbose > 2 .or. &
1076  (verbose > 1 .and. patch%level < self%levelmax .and. refine_it > 0) .or. &
1077  (verbose > 1 .and. patch%level > self%levelmin .and. refine_it < 0)) then
1078  masked = any(patch%irefine >= 0)
1079  write (io_unit%output,'(a,i6,1p,3e9.2,i4,e9.2,l3)') &
1080  'Jeans criterion: id, max, num, min, refine, maxval(d)', &
1081  patch%id, max_jeans, jeans_num, min_jeans, refine_it, maxval(d), &
1082  masked
1083  end if
1084  !---------------------------------------------------------------------------
1085  end if
1086  deallocate (pg)
1087  call trace%end()
1088 END FUNCTION check_jeans
1089 
1090 !===============================================================================
1091 !> Refine / derefine on vorticity (velocity difference per cell)
1092 !===============================================================================
1093 FUNCTION check_compress (self, patch) RESULT(refine_it)
1095  USE scaling_mod
1096  class(refine_t):: self
1097  class(solver_t), pointer:: patch
1098  integer:: refine_it
1099  !.............................................................................
1100  integer:: imin(3), l(3), u(3)
1101  real:: w_num, wmin
1102  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: test
1103  logical:: masked
1104  integer, save:: itimer=0
1105  !-----------------------------------------------------------------------------
1106  call trace%begin ('refine_t%check_compress', itimer=itimer, detailed_timer=detailed_timer)
1107  l = 1
1108  u = patch%gn
1109  allocate (test(patch%gn(1),patch%gn(2),patch%gn(3)))
1110  select type (patch)
1111  class is (solver_t)
1112  call patch%compression_magnitude (test)
1113  class default
1114  call mpi%abort("Unrecognised patch type. Abort!")
1115  end select
1116  test = minval(patch%mesh%d)*test
1117  w_num = maxval(test, mask=patch%irefine <= 0)
1118  !---------------------------------------------------------------------------
1119  ! Evaluate a per-cell criterion, which sets patch%irefine to +1 where it
1120  ! needs refinement, and sets it to 0 where the current resolution is needed.
1121  ! By implication, where it does not set patch%irefine, and where it remains
1122  ! at the intial -1 value, derefinement is allowed
1123  !---------------------------------------------------------------------------
1124  where (test > max_compress)
1125  patch%irefine = +1
1126  !---------------------------------------------------------------------------
1127  ! compression number not small enough, or overlying refinement
1128  !---------------------------------------------------------------------------
1129  else where (test > min_compress)
1130  patch%irefine = max(patch%irefine, 0)
1131  end where
1132  refine_it = 0
1133  if (w_num > max_compress) refine_it = +1
1134  if (w_num < min_compress) refine_it = -1
1135  !---------------------------------------------------------------------------
1136  if (verbose > 2 .or. &
1137  (verbose > 1 .and. patch%level < self%levelmax .and. refine_it > 0) .or. &
1138  (verbose > 1 .and. patch%level > self%levelmin .and. refine_it < 0)) then
1139  masked = any(patch%irefine >= 0)
1140  write (io_unit%output,'(a,i6,1p,3e9.2,i4,e9.2,l3)') &
1141  'compression criterion: id, max, num, min, refine, maxval(test)', &
1142  patch%id, max_compress, w_num, min_compress, refine_it, maxval(test), &
1143  masked
1144  end if
1145  deallocate (test)
1146  call trace%end (itimer, detailed_timer=detailed_timer)
1147 END FUNCTION check_compress
1148 
1149 !===============================================================================
1150 !> Refine / derefine on vorticity (velocity difference per cell)
1151 !===============================================================================
1152 FUNCTION check_vorticity (self, patch) RESULT(refine_it)
1154  USE scaling_mod
1155  class(refine_t):: self
1156  class(solver_t), pointer:: patch
1157  integer:: refine_it
1158  !.............................................................................
1159  integer:: imin(3), l(3), u(3)
1160  real:: w_num, wmin
1161  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: test
1162  logical:: masked
1163  integer, save:: itimer=0
1164  !-----------------------------------------------------------------------------
1165  call trace%begin ('refine_t%check_vorticity', itimer=itimer, detailed_timer=detailed_timer)
1166  l = 1
1167  u = patch%gn
1168  allocate (test(patch%gn(1),patch%gn(2),patch%gn(3)))
1169  select type (patch)
1170  class is (solver_t)
1171  call patch%vorticity_magnitude (test)
1172  class default
1173  call mpi%abort("Unrecognised patch type. Abort!")
1174  end select
1175  test = minval(patch%mesh%d)*test
1176  w_num = maxval(test, mask=patch%irefine <= 0)
1177  !---------------------------------------------------------------------------
1178  ! Evaluate a per-cell criterion, which sets patch%irefine to +1 where it
1179  ! needs refinement, and sets it to 0 where the current resolution is needed.
1180  ! By implication, where it does not set patch%irefine, and where it remains
1181  ! at the intial -1 value, derefinement is allowed
1182  !---------------------------------------------------------------------------
1183  where (test > max_vorticity)
1184  patch%irefine = +1
1185  !---------------------------------------------------------------------------
1186  ! vorticity number not small enough, or overlying refinement
1187  !---------------------------------------------------------------------------
1188  else where (test > min_vorticity)
1189  patch%irefine = max(patch%irefine, 0)
1190  end where
1191  refine_it = 0
1192  if (w_num > max_vorticity) refine_it = +1
1193  if (w_num < min_vorticity) refine_it = -1
1194  !---------------------------------------------------------------------------
1195  if (verbose > 2 .or. &
1196  (verbose > 1 .and. patch%level < self%levelmax .and. refine_it > 0) .or. &
1197  (verbose > 1 .and. patch%level > self%levelmin .and. refine_it < 0)) then
1198  masked = any(patch%irefine >= 0)
1199  write (io_unit%output,'(a,i6,1p,3e9.2,i4,e9.2,l3)') &
1200  'vorticity criterion: id, max, num, min, refine, maxval(test)', &
1201  patch%id, max_vorticity, w_num, min_vorticity, refine_it, maxval(test), &
1202  masked
1203  end if
1204  deallocate (test)
1205  call trace%end (itimer, detailed_timer=detailed_timer)
1206 END FUNCTION check_vorticity
1207 
1208 !===============================================================================
1209 !> Refine on shocks, with threshold set by tolerance parameter `tol_shock`.
1210 !===============================================================================
1211 FUNCTION shock_and_cd_detector (self, patch) RESULT (refine)
1213  class(refine_t):: self
1214  class(solver_t), pointer:: patch
1215  integer:: refine
1216  !.............................................................................
1217  class(mesh_t), pointer:: mx, my, mz
1218  integer:: ix, iy, iz, ione=0, jone=0, kone=0
1219  real(kind=KindScalarVar):: dp1, dp2, dp3, dv1, dv2, dv3, dd1, dd2, dd3
1220  real, parameter:: small = 1.0e-10
1221  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: pg
1222  real(kind=KindScalarVar), dimension(:,:,:,:), allocatable:: v
1223  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d, e
1224  logical:: derefine(3)
1225  !-----------------------------------------------------------------------------
1226  call trace%begin('refine_t%shocks_and_cd_detector')
1227  mx => patch%mesh(1)
1228  my => patch%mesh(2)
1229  mz => patch%mesh(3)
1230  if (mx%n > 1) ione = 1
1231  if (my%n > 1) jone = 1
1232  if (mz%n > 1) kone = 1
1233  allocate(pg(patch%mesh(1)%gn,patch%mesh(2)%gn,patch%mesh(3)%gn), &
1234  v(patch%mesh(1)%gn,patch%mesh(2)%gn,patch%mesh(3)%gn,3))
1235  ! we need the gas pressure, density and velocity
1236  select type (patch)
1237  class is (solver_t)
1238  d => patch%mem(:,:,:,patch%idx%d,patch%it,1)
1239  e => patch%mem(:,:,:,patch%idx%e,patch%it,1)
1240  pg = patch%gas_pressure()
1241  v = patch%gas_velocity()
1242  class default
1243  call mpi%abort("Unrecognised patch type. Abort!")
1244  end select
1245  !
1246  refine = 0
1247  do iz=mz%li,mz%ui
1248  do iy=my%li,my%ui
1249  do ix=mx%li,mx%ui
1250  dp1 = abs(pg(ix+ione,iy,iz) - pg(ix-ione,iy,iz)) &
1251  / (min(pg(ix+ione,iy,iz), pg(ix-ione,iy,iz)) + small )
1252  dp2 = abs(pg(ix,iy+jone,iz) - pg(ix,iy-jone,iz)) &
1253  / (min(pg(ix,iy+jone,iz), pg(ix,iy-jone,iz)) + small )
1254  dp3 = abs(pg(ix,iy,iz+kone) - pg(ix,iy,iz-kone)) &
1255  / (min(pg(ix,iy,iz+kone), pg(ix,iy,iz-kone)) + small )
1256  !
1257  if (max_shock > 0.0) then
1258  dv1 = v(ix ,iy,iz,1) - v(ix+ione,iy,iz,1)
1259  dv2 = v(ix,iy ,iz,2) - v(ix,iy+jone,iz,2)
1260  dv3 = v(ix,iy,iz ,3) - v(ix,iy,iz+kone,3)
1261  !
1262  if ( dp1 > max_shock .and. dv1 > small ) refine = 1
1263  if ( dp2 > max_shock .and. dv2 > small ) refine = 1
1264  if ( dp3 > max_shock .and. dv3 > small ) refine = 1
1265  !
1266  if ( dp1 < min_shock .and. dv1 > small ) derefine(1) = .true.
1267  if ( dp2 < min_shock .and. dv2 > small ) derefine(2) = .true.
1268  if ( dp3 < min_shock .and. dv3 > small ) derefine(3) = .true.
1269  if (all(derefine)) refine = -1
1270  end if
1271  !
1272  if (max_contact > 0.0) then
1273  dd1 = abs(d(ix+ione,iy,iz) - d(ix-ione,iy,iz)) &
1274  / (min(d(ix+ione,iy,iz), d(ix-ione,iy,iz)) + small )
1275  dd2 = abs(d(ix,iy+jone,iz) - d(ix,iy-jone,iz)) &
1276  / (min(d(ix,iy+jone,iz), d(ix,iy-jone,iz)) + small )
1277  dd3 = abs(d(ix,iy,iz+kone) - d(ix,iy,iz-kone)) &
1278  / (min(d(ix,iy,iz+kone), d(ix,iy,iz-kone)) + small )
1279  !
1280  if ( dd1 > max_contact .and. dp1 < 0.01*max_contact ) refine = 1
1281  if ( dd2 > max_contact .and. dp2 < 0.01*max_contact ) refine = 1
1282  if ( dd3 > max_contact .and. dp3 < 0.01*max_contact ) refine = 1
1283  !
1284  if ( dd1 < min_contact .and. dp1 < 0.01*max_contact ) derefine(1) = .true.
1285  if ( dd2 < min_contact .and. dp2 < 0.01*max_contact ) derefine(2) = .true.
1286  if ( dd3 < min_contact .and. dp3 < 0.01*max_contact ) derefine(3) = .true.
1287  if (all(derefine)) refine = -1
1288  end if
1289  end do
1290  end do
1291  end do
1292  !-----------------------------------------------------------------------------
1293  ! Cancel actions that would go outside the interval [levelmin...levelmax]
1294  !-----------------------------------------------------------------------------
1295  if (patch%level == self%levelmax) refine = min(refine,0)
1296  if (patch%level == self%levelmin) refine = max(refine,0)
1297  if (verbose > 1) then
1298  if (refine > 0) then
1299  write (io%output,*) 'shock_and_cd_detector: refine', patch%id, patch%level, self%levelmax
1300  else if (refine < 0) then
1301  write (io%output,*) 'shock_and_cd_detector: derefine', patch%id, patch%level, self%levelmin
1302  end if
1303  end if
1304  !
1305  nullify (mx, my, mz, d)
1306  deallocate (pg, v)
1307  !
1308  call trace%end()
1309 END FUNCTION shock_and_cd_detector
1310 
1311 !===============================================================================
1312 !> Refine on shocks, with threshold set by tolerance parameter `max_shock`.
1313 !===============================================================================
1314 FUNCTION gradient_detector (self, patch) RESULT (refine)
1316  USE solver_mod
1317  class(refine_t):: self
1318  class(solver_t), pointer:: patch
1319  integer:: refine
1320  !.............................................................................
1321  class(mesh_t), pointer:: mx, my, mz
1322  integer:: i, j, k, ione=0, jone=0, kone=0, im1, jm1, km1, ip1, jp1, kp1, n
1323  real, parameter:: small=1e-10
1324  real(kind=KindScalarVar):: dqmax
1325  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: q, dq1, dq2, dq3
1326  real, parameter:: rsmall = 1.0e-10
1327  logical:: is_unsigned
1328  character(len=1):: unsigned(3) = ['d', 'p', 'e']
1329  !-----------------------------------------------------------------------------
1330  call trace%begin('refine_t%gradient_detector')
1331  allocate (dq1(patch%gn(1),patch%gn(2),patch%gn(3)), &
1332  dq2(patch%gn(1),patch%gn(2),patch%gn(3)), &
1333  dq3(patch%gn(1),patch%gn(2),patch%gn(3)), &
1334  q(patch%gn(1),patch%gn(2),patch%gn(3)))
1335  mx => patch%mesh(1)
1336  my => patch%mesh(2)
1337  mz => patch%mesh(3)
1338  if (mx%n > 1) ione = 1
1339  if (my%n > 1) jone = 1
1340  if (mz%n > 1) kone = 1
1341  refine = 0
1342  do n=1,size(grad_var)
1343  if (trim(grad_var(n)) == '') exit
1344  if (trim(grad_var(n)) == 'd') q = patch%mem(:,:,:,patch%idx%d,patch%it,1)
1345  if (trim(grad_var(n)) == 'e') q = patch%mem(:,:,:,patch%idx%e,patch%it,1)
1346  if (trim(grad_var(n)) == 'p1') q = patch%mem(:,:,:,patch%idx%px,patch%it,1)
1347  if (trim(grad_var(n)) == 'p2') q = patch%mem(:,:,:,patch%idx%py,patch%it,1)
1348  if (trim(grad_var(n)) == 'p3') q = patch%mem(:,:,:,patch%idx%pz,patch%it,1)
1349  if (trim(grad_var(n)) == 'b1') q = patch%mem(:,:,:,patch%idx%bx,patch%it,1)
1350  if (trim(grad_var(n)) == 'b2') q = patch%mem(:,:,:,patch%idx%by,patch%it,1)
1351  if (trim(grad_var(n)) == 'b3') q = patch%mem(:,:,:,patch%idx%bz,patch%it,1)
1352  select type(patch)
1353  class is (solver_t)
1354  if (trim(grad_var(n)) == 'p') q = patch%gas_pressure()
1355  if (trim(grad_var(n)) == 'v1') q = patch%gas_velocity(1)
1356  if (trim(grad_var(n)) == 'v2') q = patch%gas_velocity(2)
1357  if (trim(grad_var(n)) == 'v3') q = patch%gas_velocity(3)
1358  class default
1359  call mpi%abort("Unrecognised patch type. Abort!")
1360  end select
1361  is_unsigned = .false.
1362  do i=1,size(unsigned)
1363  if (grad_var(n) == unsigned(i)) then
1364  is_unsigned = .true.
1365  exit
1366  end if
1367  end do
1368  ! calculate the gradient in each direction
1369  do k=mz%lo,mz%uo
1370  km1 = k - kone
1371  kp1 = k + kone
1372  do j=my%lo,my%uo
1373  jm1 = j - jone
1374  jp1 = j + jone
1375  if (is_unsigned) then
1376  do i=mx%lo,mx%uo
1377  im1 = i - ione
1378  ip1 = i + ione
1379  dq1(i,j,k) = abs( q(ip1,j,k) - q(im1,j,k) ) &
1380  / ( max( q(ip1,j,k), q(im1,j,k) ) &
1381  + small )
1382  dq2(i,j,k) = abs( q(i,jp1,k) - q(i,jm1,k) ) &
1383  / ( max( q(i,jp1,k), q(i,jm1,k) ) &
1384  + small )
1385  dq3(i,j,k) = abs( q(i,j,kp1) - q(i,j,km1) ) &
1386  / ( max( q(i,j,kp1), q(i,j,km1) ) &
1387  + small )
1388  end do
1389  else
1390  do i=mx%lo,mx%uo
1391  im1 = i - ione
1392  ip1 = i + ione
1393  dq1(i,j,k) = abs( q(ip1,j,k) - q(im1,j,k) )
1394  dq2(i,j,k) = abs( q(i,jp1,k) - q(i,jm1,k) )
1395  dq3(i,j,k) = abs( q(i,j,kp1) - q(i,j,km1) )
1396  end do
1397  end if
1398  end do
1399  end do
1400  ! choose max. gradient; if the max. is greater than the threshold, flag
1401  ! for refinement
1402  do k=mz%li+kone,mz%ui-kone
1403  km1 = k - kone
1404  kp1 = k + kone
1405  do j=my%li+jone,my%ui-jone
1406  jm1 = j - jone
1407  jp1 = j + jone
1408  do i=mx%li+ione,mx%ui-ione
1409  im1 = i - ione
1410  ip1 = i + ione
1411  dqmax = max( dq1(i,j,k), dq1(ip1,j,k), dq1(im1,j,k) &
1412  , dq2(i,j,k), dq2(i,jp1,k), dq2(i,jm1,k) &
1413  , dq3(i,j,k), dq3(i,j,kp1), dq3(i,j,km1) )
1414  if (patch%level < self%levelmax .and. dqmax > max_grad(n)) refine = 1
1415  if (patch%level > self%levelmin .and. dqmax < min_grad(n)) refine = -1
1416  end do
1417  end do
1418  end do
1419  if (verbose > 1) then
1420  if (refine > 0) then
1421  write (io%output,*) 'gradient_detector: refine ', grad_var(n), patch%id
1422  else if (refine < 0) then
1423  write (io%output,*) 'gradient_detector: derefine ', grad_var(n), patch%id
1424  end if
1425  end if
1426  end do
1427  deallocate (dq1, dq2, dq3, q)
1428  call trace%end()
1429 END FUNCTION gradient_detector
1430 
1431 !===============================================================================
1432 !> Given a box defined by its lower and upper corners, flag a patch for refinement
1433 !> if it lies inside this box, up to a maximum level.
1434 !> This particular criteria is primarily used for testing and debugging.
1435 !===============================================================================
1436 FUNCTION refine_region (self, patch) RESULT (refine)
1438  class(refine_t) :: self
1439  class(solver_t), pointer :: patch
1440  integer:: refine
1441  class(mesh_t), pointer:: m1, m2, m3, m(:)
1442  logical :: intersect
1443  logical :: refine_it
1444  integer :: ix, iy, iz
1445  real :: x, y, z
1446  real, dimension(3) :: p_llc, p_urc, b_llc, b_urc
1447  !.............................................................................
1448  call trace%begin('refine_t%refine_on_box')
1449  refine_it = .false.
1450  m1 => patch%mesh(1)
1451  m2 => patch%mesh(2)
1452  m3 => patch%mesh(3)
1453  m => patch%mesh(:)
1454  ! check if the patch (as a whole) intersects the box-shaped region of interest
1455  do ix=1,3
1456  associate(mr => patch%mesh(ix)%r)
1457  p_llc(ix) = m(ix)%p + mr(m(ix)%li) - 0.5 * m(ix)%d
1458  p_urc(ix) = m(ix)%p + mr(m(ix)%ui) + 0.5 * m(ix)%d
1459  end associate
1460  end do
1461  b_llc = region_llc
1462  b_urc = region_urc
1463  intersect = all((p_urc - b_llc) > 0) .and. all((b_urc - p_llc) > 0)
1464  if (.not. intersect) then
1465  call trace%end ()
1466  return
1467  endif
1468  do iz=m3%li,m3%ui
1469  z = m3%p + m3%r(iz)
1470  do iy=m2%li,m2%ui
1471  y = m2%p + m2%r(iy)
1472  do ix=m1%li,m1%ui
1473  x = m1%p + m1%r(ix)
1474  refine_it = refine_it .or. (x >= b_llc(1) .and. x <= b_urc(1) .and. &
1475  y >= b_llc(2) .and. y <= b_urc(2) .and. &
1476  z >= b_llc(3) .and. z <= b_urc(3) )
1477  end do
1478  end do
1479  end do
1480  refine = merge(1,0,refine_it)
1481  call trace%end()
1482 END FUNCTION refine_region
1483 
1484 !===============================================================================
1485 !> Check that the derefinement criterion is sane
1486 !===============================================================================
1487 SUBROUTINE sanity_check (self, min, max, label, reverse)
1488  class(refine_t):: self
1489  real:: min, max
1490  character(len=*):: label
1491  logical, optional:: reverse
1492  !.............................................................................
1493  if (max == 0.0) then
1494  if (present(reverse)) then
1495  min = max/self%safe_ratio
1496  else
1497  max = min*self%safe_ratio
1498  end if
1499  else if (max < min*2.0) then
1500  if (present(reverse)) then
1501  min = max/2.0
1502  else
1503  max = min*2.0
1504  end if
1505  if (io%master) then
1506  write (stdout,*) 'WARNING: '//trim(label)//' max value set to', max
1507  end if
1508  end if
1509 END SUBROUTINE sanity_check
1510 
1511 !===============================================================================
1512 !> Generic task to patch function
1513 !===============================================================================
1514 FUNCTION cast2solver (task) RESULT (solver)
1515  class(task_t), pointer:: task
1516  class(solver_t), pointer:: solver
1517  select type (task)
1518  class is (solver_t)
1519  solver => task
1520  class default
1521  nullify (solver)
1522  end select
1523 END FUNCTION cast2solver
1524 
1525 END MODULE refine_mod
Each thread uses a private timer data type, with arrays for start time and total time for each regist...
Definition: timer_mod.f90:11
download_link: takes care of downloads to linktask same: called for patches on the same level differ:...
Template version of a module intended for adding extra features, at a level between the basic task an...
Definition: extras_mod.f90:54
Support tic/toc timing, as in MATLAB, and accurate wallclock() function. The timing is generally much...
Module with list handling for generic class task_t objects.
Definition: list_mod.f90:4
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Interface from gpatch_mod to a choice of binary data I/O methods, controlled by the iomethod text str...
Definition: data_io_mod.f90:17
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
This module handles checking max change between neighboring points. Each instance of it needs an inde...
Definition: refine_mod.f90:168
This module is a placeholder for shared information, solving the problem of making a piece of informa...
Definition: shared_mod.f90:7
Fundamental constants in CGS and SI units.
Definition: math_mod.f90:4
Template module for mesh.
Definition: mesh_mod.f90:4
Definition: io_mod.f90:4
The lock module uses nested locks, to allow versatile use of locks, where a procedure may want to mak...
This module contains all experiment specific information necessary to solve the heat diffusion proble...
Definition: solver_mod.f90:5
Template module for tasks.
Definition: task_mod.f90:4